<figure><img src="image_1.png" /></figure>

HAL
open science

# Wave propagation in random 2D turbulence: a multi-scale approach

Valentin Resseguier, Erwan Hascoët, Bertrand Chapron

► To cite this version:

Valentin Resseguier, Erwan Hascoët, Bertrand Chapron. Wave propagation in random 2D turbulence: a multi-scale approach. 2024. hal-04549029

HAL Id: hal-04549029

https://hal.science/hal-04549029

Preprint submitted on 16 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Wave propagation in random 2D turbulence: a multi-scale approach

Valentin Resseguier, Erwan Hascoët, Bertrand Chapron

► To cite this version:

Valentin Resseguier, Erwan Hascoët, Bertrand Chapron. Wave propagation in random 2D turbulence: a multi-scale approach. Journal of Fluid Mechanics, 2024, 997, pp.A51. 10.1017/jfm.2024.769. hal-04549029v3

HAL Id: hal-04549029

https://hal.science/hal-04549029v3

Submitted on 22 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Banner appropriate to article type will appear here in typeset article

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

# Wave propagation in random 2D turbulence: a multi-scale approach

Valentin Resseguier¹,²†, Erwan Hascoët³ and Bertrand Chapron⁴

¹ACTA, OPAALE, INRAE

²Lab, SCALIAN DS

³Oceandatalab

⁴LOPS, Ifremer

(Received xx; revised xx; accepted xx)

To study two-dimensional dispersive waves propagating over turbulent flows, a new and less restrictive fast waves approximation is proposed using a multiscale setting. In this ansatz, correlation lengths of the random small scale turbulence components can be considered negligible in the wave packet propagating frame. Still, 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 robustly be calibrated 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.

19

Key words:

20

MSC Codes 76B15,76M35, 82B31,82C31, 37H10, 37L55

21

## 1. Introduction

22

23

24

25

26

27

28

29

30

31

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 et al. 2019; Kafiabad et al. 2019; Bôas & Young 2020; Garnier et al. 2020). Some of these approaches can be traced back to wave-wave interactions models, e.g. McComas & Bretherton (1977) (see also Kafiabad et al. 2019, and reference therein). In most cases, the central assumption is either delta-correlated turbulence (Voronovich 1976; Klyatskin 2005; Klyatskin & Koshel 2015) and/or fast waves in comparison to fluid flow

† Email address for correspondence: valentin.resseguier@inrae.fr
Banner appropriate to article type will appear here in typeset article

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

# Wave propagation in random 2D turbulence: a multi-scale approach

Valentin Resseguier¹,²†, Erwan Hascoët³ and Bertrand Chapron⁴

¹INRAE, UR OPAALE, 17 avenue de Cucillé, F-35044 Rennes, France

²LAB, SCALIAN DS, 2 Rue Antoine Becquerel, 35700 Rennes, France

³Oceandatalab, 870 Rte de Deolen, 29280 Locmaria-Plouzané, France

⁴LOPS, Ifremer, 1625 Rte de Sainte-Anne, 29280 Plouzané

(Received xx; revised xx; accepted xx)

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 scale and small scale 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. Still, 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 robustly be calibrated 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.

## Key words:

21

MSC Codes 76B15,76M35, 82B31,82C31, 37H10, 37L55

22

## 1. Introduction

23

24

25

26

27

28

29

30

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 et al. 2019; Kafiabad et al. 2019; Bôas & Young 2020; Garnier et al. 2020). Some of these approaches can be traced back to wave-wave interactions models, e.g. McComas & Bretherton (1977) (see also Kafiabad et al. 2019, and reference therein). In most

† Email address for correspondence: valentin.resseguier@inrae.fr
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 et al. 2023; Wang et al. 2023). Medium variations may be slow and delta-correlations are hardly justifiable in a fixed frame. Though, 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 implies conservation along ray of intrinsic frequency, wavenumber, and group velocity magnitude in two dimension (Boury et al. 2023). Subsequently, most of wave dynamics closures 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 solely been derived for a zero large-scale current. More generally, fast wave models mostly rely either on zero or constant current components at larger scales. West (1978), for instance, discussed acoustic waves in a 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 et al. 2023), including dynamics under Location Uncertainty (LU) (Mémin 2014; Resseguier et al. 2017a) and Stochastic Advection by Lie Transport (SALT) (Holm 2015). Under this framework, resulting conservative stochastic closures build on time delta-correlation for the small-scale turbulence component (Cotter et al. 2017). 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 et al. 2023). Considering a single-wavevector current, solutions for a monochromatic shallow water wave were developed 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 principle of ray tracing method, we present the multiscale framework for fast waves dynamics, its physical ground and a calibration method for the closure. Simplified stochastic equations are then derived for ray dynamics and 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 2D 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 train, it follows a form $h(x, t)e^{i\phi(x, t)}$ + c.c., for most properties. If a packet is to be followed, the phase, $\phi(x, t)$, must smoothly vary along the propagation, i.e. $\phi(x, t)$ is differentiable. The relative frequency is then $\omega = -\partial_t \phi(x, t)$, and the wave number vector $k = \nabla\phi(x, t)$, with wavenumber $k = \|k\|$ and direction given by the normalized wave-vector, $\tilde{k} = k/k = (\cos\theta_k, \sin\theta_k)$. To first order, such a train of waves is dispersive and the intrinsic frequency reads

$$ \omega - k \cdot v = \omega_0 = \begin{cases} \text{cst.} \frac{1}{\alpha} k^{\alpha}, & \alpha \neq 0 \\ \text{cst.} \log(k), & \alpha = 0 \end{cases} \qquad (2.1) $$
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 et al. 2023; Wang et al. 2023). Medium variations may be slow and delta-correlations are hardly justifiable in a fixed frame. Though, 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 implies conservation along ray of intrinsic frequency, wavenumber, and group velocity magnitude in two dimension (Boury et al. 2023). Subsequently, most of 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 solely been derived for a zero large-scale current. More generally, fast wave models mostly rely either on zero or constant current components at larger scales. West (1978), for instance, discussed acoustic waves in a 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 et al. 2023), including dynamics under Location Uncertainty (LU) (Mémin 2014; Resseguier et al. 2017a) and Stochastic Advection by Lie Transport (SALT) (Holm 2015). Under this framework, the small-scale velocity component is delta-correlated in time (Cotter et al. 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 et al. 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 2D 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 train, it follows a form $h(\mathbf{x}, t)e^{i\phi(\mathbf{x},t)} + \text{c.c.}$, for most properties. Typically, $h$ would be the local wave height with 70 meters units. If a packet is to be followed, the phase, $\phi(\mathbf{x}, t)$, must smoothly vary along the 71 propagation, i.e. $\phi(\mathbf{x}, t)$ is differentiable. The relative frequency is then $\omega = -\partial_t\phi(\mathbf{x}, t)$, and 72 the wave number vector $\mathbf{k} = \nabla\phi(\mathbf{x}, t)$, with wavenumber $k = \|\mathbf{k}\|$ and direction given by 73 the normalized wave-vector, $\tilde{\mathbf{k}} = \mathbf{k}/k = \left(\frac{\cos\theta_k}{\sin\theta_k}\right)$. To first order, such a train of waves is 74
and propagates with its group velocity $v_g = \nabla_k \omega$, constantly modified by the local velocity of the currents $v$, 73 74

$$ \frac{dx_r}{dt} = v_g = v_g^0 + v, \qquad (2.2) $$

where $x_r$ is the centroid of a wave group, $v_g^0 = \frac{\partial \omega_0(k)}{\partial k} \tilde{k}$ is the group velocity without currents, i.e. solely depending on the wave vector. For $\alpha = 1$, the medium is non-dispersive (e.g. acoustic waves). $\alpha = 1/2$ corresponds to gravity waves over deep ocean ($\omega_0 = \sqrt{gk}$). The dominant wave-vector $k$ within the group evolves according to 76 77 78 79

$$ \frac{dk}{dt} = -\nabla v^T k. \qquad (2.3) $$

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

To complete the wave field description, the wave action, $A(x, t) = h^2(x, t)/\omega_0(k(x, t))$, is considered to be an adiabatic invariant in 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 wave-vectors of the action spectrum, $N$, also related to the wave energy spectrum, $E$: 89 90 91 92 93 94

$$ A(x, t) = \int dk \, N(x, k, t) = \int dk \, \frac{E(x, k, t)}{\omega_0(k, t)}. \qquad (2.4) $$

Action and energy spectrum quantify action and energy by unit of surface (unit of $x$) and by unit of wave-vector surface (unit of $k$). Consider the $(x, k)$ variable change between different times $t_i$ and $t_f$ integrating the characteristic eikonal equations (2.2)-(2.3) 96 97 98

$$ \begin{pmatrix} x_r(t_i) \\ k(t_i) \end{pmatrix} \mapsto \begin{pmatrix} x_r(t_f) \\ k(t_f) \end{pmatrix} \qquad (2.5) $$

By the Liouville theorem and since $\nabla_x \cdot \frac{dx_r}{dt} + \nabla_k \cdot \frac{dk}{dt} = 0$, the state-space of the "packet-by-packet" approach (the $(x, k)$ space) does not contract nor dilates along time. Therefore, if wave dissipation is neglected, the wave action spectrum $N$ is conserved Lavrenov (2013), i.e. 100 101 102 103

$$ N(x_r(t_i), k(t_i), t_i) = N(x_r(t_f), k(t_f), t_f). \qquad (2.6) $$

This result is extremely useful because it only involves 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 105 106 107

$$ E(x_r(t_f), k(t_f), t_f) = \frac{\omega_0(k(t_f))}{\omega_0(k(t_i))} E(x_r(t_i), k(t_i), t_i). \qquad (2.7) $$

starting with an initial incoming wave spectrum $E(x_r(t_i), k(t_i), t_i)$ for every wave-vectors $k(t_i)$, starting from a small set of spatial points $x_r(t_i)$. 109 110
76 dispersive and the intrinsic frequency reads

$$
77 \qquad \omega - \mathbf{k} \cdot \mathbf{v} = \omega_0 = \begin{cases} \operatorname{cst.} \frac{1}{\alpha} k^{\alpha}, & \alpha \neq 0 \\ \operatorname{cst.} \log(k), & \alpha = 0 \end{cases} \qquad (2.1)
$$

78 and propagates with its group velocity $\mathbf{v}_g = \nabla_k \omega$, constantly modified by the local velocity
79 of the currents $\mathbf{v}$,

$$
80 \qquad \frac{d\mathbf{x}_r}{dt} = \mathbf{v}_g = \mathbf{v}_g^0 + \mathbf{v}, \qquad (2.2)
$$

81 where $\mathbf{x}_r$ is the centroid of a wave group, $\mathbf{v}_g^0 = \frac{\partial \omega_0(k)}{\partial k} \tilde{\mathbf{k}}$ is the group velocity without
82 currents, i.e. solely depending on the wave vector. For $\alpha = 1$, the medium is non-dispersive
83 (e.g. acoustic waves). $\alpha = 1/2$ corresponds to gravity waves over deep ocean ($\omega_0 = \sqrt{g\tilde{k}}$).
84 The dominant wave-vector $\mathbf{k}$ within the group evolves according to

$$
85 \qquad \frac{d\mathbf{k}}{dt} = -\nabla \mathbf{v}^T \mathbf{k}. \qquad (2.3)
$$

86 Equations (2.2)-(2.3) are the Hamilton’s eikonal equations. Along the propagating ray,
87 velocity gradients induce linear variations. Decelerating currents will, for instance, shorten
88 waves, and reduce the group velocity. Traveling over fields of random velocities **v**, the wave-
89 vector **k** will also become randomly distributed. Scattering of ocean surface wave packets
90 by random currents can generally be assumed to be weak, with ||**v**|| of order 0.5 m.s⁻¹,
91 much smaller than v_g⁰ = ||v_g⁰|| of order 10 m.s⁻¹. Yet, cumulative effects of these random
92 surface currents can lead to strong convergence or divergence between initially nearby ray
93 trajectories.
94

94 To complete the wave field description, $E(\mathbf{x}, t) = \frac{1}{2}\rho g h^2(\mathbf{x}, t)$ and $A(\mathbf{x}, t) = E(\mathbf{x}, t)/\omega_0(k(\mathbf{x}, t))$ denote energy and action by unit of surface. $E$ is expressed in $J/m^2$ and $A$ in $J.s/m^2$. To avoid spurious notations, we set the multiplicative constant $\frac{1}{2}\rho g$ to unity. The wave action is considered to be an adiabatic invariant in 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 wave-vectors of the action spectrum, $N$, also related to the wave energy spectrum, $E$:

$$
102 \qquad A(\mathbf{x}, t) = \int d\mathbf{k} \, N(\mathbf{x}, \mathbf{k}, t) = \int d\mathbf{k} \, \frac{E(\mathbf{x}, \mathbf{k}, t)}{\omega_0(\mathbf{k}, t)}. \qquad (2.4)
$$

103 Action and energy spectrum quantify action and energy by unit of surface (unit of x) and by
104 unit of wave-vector surface (unit of k). Consider the (x, k) variable change between different
times tᵢ and tⱼ integrating the characteristic eikonal equations (2.2)-(2.3)
105

$$
106 \qquad \begin{pmatrix} \mathbf{x}_r(t_i) \\ \mathbf{k}(t_i) \end{pmatrix} \mapsto \begin{pmatrix} \mathbf{x}_r(t_f) \\ \mathbf{k}(t_f) \end{pmatrix}. \qquad (2.5)
$$

107 According to the Liouville theorem for Hamiltonian mechanics (Landau & Lifshits 1960,
108 §46), the state-space of the "packet-by-packet" approach (the $(\mathbf{x}, \mathbf{k})$ space) does not contract
109 nor dilates along time. Readers not familiar with Hamiltonian dynamics may see the
110 divergence free of the 4-dimensional flow (2.5) – i.e. $\nabla_{\mathbf{x}} \cdot \frac{d\mathbf{x}_r}{dt} + \nabla_{\mathbf{k}} \cdot \frac{d\mathbf{k}}{dt} = 0$ – as the
111 divergence free of incompressible flow velocities, leading naturally to volume-preserving
dynamics. Therefore, if wave dissipation is neglected, the wave action spectrum $N$ is
112 conserved (Lavrenov (2013)), i.e.
113

$$
114 \qquad N(\mathbf{x}_r(t_i), \mathbf{k}(t_i), t_i) = N(\mathbf{x}_r(t_f), \mathbf{k}(t_f), t_f). \qquad (2.6)
$$
### 3. A new fast wave assumption

Commonly, the Eulerian current $v$ is decomposed into a low-frequency large-scale component

$\bar{v}$ and a transient small-scale unresolved component $v'$:

$$v = \bar{v} + v'. \quad (3.1)$$

From now, we shall consider divergence-free two-dimensional 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 $v'$, its Eulerian spatio-temporal covariance is considered, assuming

statistical homogeneity and stationarity for the Eulerian velocity $v'_E(t, \mathbf{x}) = v'(t, \mathbf{x})$

$$C_{ij}^{v'_E}(\delta t, \delta \mathbf{x}) = \mathbb{E} \left( v'_i(t, \mathbf{x}) v'_j(t + \delta t, \mathbf{x} + \delta \mathbf{x}) \right) = \mathbb{E} \left( v'_i(t, \mathbf{x}_r(t)) v'_j(t + \delta t, \mathbf{x}_r(t) + \delta \mathbf{x}) \right) \quad (3.2)$$

where $\mathbf{x}_r$ is solution of (2.2) with an arbitrary initial position $\mathbf{x}_r^0$. Then, we define, $v'_R(t) =$

$v'(t, \mathbf{x}_r(t))$, the Lagrangian velocity along the ray $\mathbf{x}_r(t)$. The temporal covariance of the

small-scale component $v'$ – in the wave group frame – is the covariance of that Lagrangian velocity:

$$C_{ij}^{v'_R}(\delta t) = \mathbb{E} \left( v'_i(t, \mathbf{x}_r(t)) v'_j(t + \delta t, \mathbf{x}_r(t + \delta t)) \right) = C_{ij}^{v'_E}(\delta t, \mathbf{x}_r(t + \delta t) - \mathbf{x}_r(t)), \quad (3.3)$$

Assume a typical isotropic form for this covariance:

$$C^{v'_E}(\delta t, \delta \mathbf{x}) = C \left( \frac{|\delta t|}{\tau_{v'}} + \frac{\|\delta \mathbf{x}\|}{l_{v'}} \right), \quad (3.4)$$

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

$$C^{v'_R}(\delta t) = C \left( \frac{|\delta t|}{\tau_{v'}} + \frac{\|\mathbf{x}_r(t' + t) - \mathbf{x}_r(t')\|}{l_{v'}} \right) = C \left( \left( \frac{1}{\tau_{v'}} + \frac{\|\mathbf{v}_g\|}{l_{v'}} \right) |\delta t| + O(\delta t^2) \right). \quad (3.5)$$

Therefore, $(\frac{1}{\tau_{v'}} + \frac{\|\mathbf{v}_g\|}{l_{v'}})^{-1}$ is the correlation time of $v'(t, \mathbf{x}_r(t))$. The same derivation is valid

for $\nabla(v')^\mathrm{T}(t, \mathbf{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$. The ratio $\epsilon$, between this along-ray correlation time and

the characteristic time of the wave group properties evolution, will then control the time

decorrelation assumption of $v'$:

$$\epsilon = \frac{l_{v'}}{v_g^0} \|\nabla v^\mathrm{T}\| \sim \frac{l_{v'}}{l_v} \frac{\|v\|}{v_g^0}. \quad (3.6)$$

Note the Eulerian small-scale velocity $v'$ is not necessarily time uncorrelated, as assumed

in Voronovich (1976); Klyatskin & Koshel (2015). Yet, for small enough $\epsilon$, the Lagrangian

small-scale velocity along the ray can be considered time uncorrelated. From the expression

of $\epsilon$, such a condition depends upon:

* $v_g^0$, the fast wave group velocity

* $\|v\|$, often slow but not always negligible compared to the intrinsic wave group, $v_g^0$.

* $l_{v'}/l_v$, related to the separation between large scales $\bar{v}$ and small scales $v'$, e.g. the spatial filtering cutoff of the large-scale velocity $\bar{v}$, but also related to its kinetic energy (KE)

distribution over spatial scales, typically the spectrum slope.
This along-ray partially 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 $\frac{\|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 $\frac{l_{v'}}{l_v} \ll 1$. Similarly, this last validity criterion can be obtained replacing in (3.2)-(3.6) $\mathbf{x}_r$ by the fluid particle Lagrangian path $\mathbf{x}$ (solution of $\frac{dx}{dt} = \mathbf{v}$) and thus $\mathbf{v}_g^0$ by $\mathbf{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 shall 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 $\mathbf{v}'$.

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 $\frac{dx(t)}{dt} = \mathbf{v}_L(t) = \mathbf{v}(t, \mathbf{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)

$$ \frac{1}{2} a^L = \int_0^\infty d\delta t \, C^{v'_L}(\delta t) = \int_0^\infty d\delta t \, C^{v'_E}(\delta t, \mathbf{x}(t+\delta t) - \mathbf{x}(t)) \approx \frac{1}{2} \tau_{v'} C(0). \quad (3.7) $$

Along a propagating wave group, a ray absolute diffusivity occurs and slightly differs from the usual absolute diffusivity to become

$$ \frac{1}{2} a^R = \int_0^\infty d\delta t \, C^{v'_R}(\delta t) \approx \frac{1}{2} \left( \frac{1}{\tau_{v'}} + \frac{\|v_g\|}{l_{v'}} \right)^{-1} C(0) \approx \frac{1}{2} \frac{l_{v'}}{v_g^0} C(0). \quad (3.8) $$

The absolute diffusivity sets the amplitude of the small scale velocity $\mathbf{v}'$. To calibrate its spatial correlations, we may focus on its Fourier transform, $\hat{v}'(\mathbf{k}, t)$, denoting by $\mathbf{k} = \mathbf{\kappa} \begin{pmatrix} \cos \theta_{\mathbf{k}} \\ \sin \theta_{\mathbf{k}} \end{pmatrix}$, the surface current wave-vector. By analogy with the current kinetic energy spectra $E_{\mathbf{k}} = \frac{1}{2} \oint_0^{2\pi} d\theta_{\mathbf{k}} \kappa \frac{\|\hat{v}(\mathbf{k}, t)\|^2}{(2\pi)^2}$, Resseguier et al. (2017b, 2020b) decompose the absolute diffusivity scale by scale. Referring it to Absolute Diffusivity Spectral Density (ADSD), they defined it by

$$ \frac{1}{2} A^R(\mathbf{k}) = \frac{1}{2} \frac{1/\kappa}{v_g^0(\mathbf{k})} E_k(\mathbf{k}), \quad (3.9) $$

where $k$ denotes the wave wavenumber and $\kappa$ the current wavenumber. To calibrate an equivalent noise, we model $\mathbf{v}'$ by $\sigma dB_t/dt$, where $dB_t/dt$ is a spatio-temporal white noise and $\sigma$ denotes a spatial filtering operator which encodes spatial correlations through $A_{v'}^R$ and the horizontal incompressibility condition $(\nabla \cdot \mathbf{\sigma}) = 0$. Specifically, we assume an ADSD power law, $A^R(\mathbf{k}) \approx A_0^R \kappa^{-\mu}$, and we choose a divergence-free isotropic spatial filter $\nabla^\perp \tilde{\psi}_\sigma$ such that

$$ \mathbf{v}' = \sigma dB_t/dt = \nabla^\perp \tilde{\psi}_\sigma \star dB_t/dt, \quad (3.10) $$
148 assumption of $v'$:

$$ \epsilon = \frac{l_{v'}}{v_g^0} \|\nabla v^T\| \sim \frac{l_{v'}}{l_v} \frac{\|v\|}{v_g^0}. \qquad (3.6) $$

150 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 the Eulerian small-scale velocity $v'$ is not necessarily time uncorrelated, as assumed in Voronovich (1991); 152 Klyatskin & Koshel (2015). Yet, for small enough $\epsilon$, the Lagrangian small-scale velocity 153 along the ray can be considered time uncorrelated. From the expression of $\epsilon$, such a condition 154 depends upon: 155

• $v_g^0$, the fast wave group velocity

• $\|v\|$, often slow but not always negligible compared to the intrinsic wave group, $v_g^0$.

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

161 This along-ray partial time-correlation assumption is less restrictive than the usual fast 162 wave approximation (White & Fornberg 1998; Dysthe 2001; Bal & Chou 2002; Borcea 163 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 $\|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; 165 Resseguier et al. 2020a) – say $l_{v'}/l_v \ll 1$. Similarly, this last validity criterion can be obtained 166 replacing in (3.2)-(3.6) $x_r$ by the fluid particle Lagrangian path $x$ (solution of $\frac{dx}{dt} = v$) 167 and thus $v_g^0$ by $v$. These asymptotic models often rely on averaging or homogenization 168 techniques (Papanicolaou & Kohler 1974; White & Fornberg 1998) to derive Markovian 169 dynamics involving various types of diffusivity. 170

### 3.2. Ray absolute diffusivity and turbulence statistics: calibration

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

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

$$ \frac{1}{2} a^L = \int_0^{+\infty} d\delta t \, C^{v'L} (\delta t) = \int_0^{+\infty} d\delta t \, C^{v'E} (\delta t, x(t+\delta t) - x(t)) \approx \frac{1}{2} \tau_{v'} C(0). \quad (3.7) $$

183 This diffusivity well describes effects of fast-varying eddies but is not appropriate in our case. 184 Indeed, along a propagating wave group, $\frac{dx_r(t)}{dt} = v_g^0(t) + v_R(t)$, a ray absolute diffusivity 185 occurs and slightly differs from the usual absolute diffusivity to become 186

$$ \frac{1}{2} a^R = \int_0^{+\infty} d\delta t \, C^{v'R} (\delta t) \approx \frac{1}{2} \left( \frac{1}{\tau_{v'}} + \frac{\|v_g\|}{l_{v'}} \right)^{-1} C(0) \approx \frac{1}{2} \frac{l_{v'}}{v_g^0} C(0). \quad (3.8) $$

187 The absolute diffusivity sets the amplitude of the small scale velocity $v'$. Indeed, since the 188
<figure><img src="image_2.png" /><figcaption>Figure 1: KE spectrum ($m^2 \cdot s^{-2}/(\text{rad} \cdot m^{-1})$) (left) and ADSD ($m^2 \cdot s^{-1}/(\text{rad} \cdot m^{-1})$) (right) of the resolved high-resolution velocity, $A^R$, in brown, low-resolution velocity, $A_{\tilde{v}}^R$, in blue, and modeled stochastic velocity, $A_{v'}^R(\kappa) = A_0^R \kappa^{-\mu} - A_{\tilde{v}}^R(\kappa)$, in red. For the ADSD power law, $A^R(\kappa) \approx A_0^R \kappa^{-\mu}$, we impose the theoretical KE spectrum slope $-5/3$ (black solid line), coherently with homogeneous SQG dynamics (see "Numerical results" section). The residual ADSD (red line) is set to extrapolate that power law at small scales.</figcaption></figure>

and $\frac{\partial \tilde{\kappa}^3}{(2\pi)^2} |\tilde{\psi}_\sigma(\kappa)|^2 = \frac{1}{2} \oint_0^{2\pi} d\theta_\kappa \kappa \frac{\|\overline{\sigma dB_t(\kappa)}\|^2}{(2\pi)^2 dt} = A_{v'}^R(\kappa)$. The power law assumption enables automatic closure calibration: $A_{v'}^R(\kappa) = A_0^R \kappa^{-\mu} - A_{\tilde{v}}^R(\kappa)$, from instantaneous large-scale current statistics $A_{\tilde{v}}^R$ only (Resseguier et al. 2020b) as illustrated in figure 1.

## 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 $v$ by $\tilde{v} + \sigma \circ dB_t/dt$. Then, the stochastic transport of phase, $d\tilde{v} = \omega_0(\|\nabla\phi\|)$, i.e. – up to that velocity replacement – the Stratonovich dispersion relation is exactly (2.1). The method of characteristics also applies. Note, one can switch from Stratonovich to Itō notations, where $v'$ corresponds to $\sigma dB_t/dt$. The characteristics equations (2.2)-(2.3) also remain unchanged for homogeneous and isotropic $v'$:

$$ \left\{ \begin{aligned} dx_r &= (v_g^0 + \tilde{v})dt + \sigma dB_t, \\ dk &= -\nabla(\tilde{v}dt + \sigma dB_t)^T k. \end{aligned} \right. \qquad (4.1) $$

### 4.1. Single-ray stochastic differential equations

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

Notably, dynamics of wave-vectors (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 et al. (1999) and references therein. Without loss of generality, the large-scale velocity can be parameterized as

$$ \bar{v} = \bar{v} \begin{pmatrix} \cos \bar{\theta} \\ \sin \bar{\theta} \end{pmatrix} \quad \text{and} \quad \nabla \bar{v}^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}. \qquad (4.2) $$

The influence of the large-scale currents on the ray dynamics (4.1), expressed in the local crest-oriented frame ($\tilde{k}, \tilde{k}^\perp$), is straightforward (Lapeyre et al. 1999). The small-scale
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, $\hat{v}'(\kappa, t)$, denoting by $\kappa = \kappa \begin{pmatrix} \cos \theta_{\kappa} \\ \sin \theta_{\kappa} \end{pmatrix}$, the surface current

wave-vector. By analogy with the current kinetic energy spectra $E_{\kappa} = \frac{1}{2} \oint_{0}^{2\pi} d\theta_{\kappa} \kappa \frac{\|\hat{v}(\kappa, t)\|^2}{(2\pi)^2}$,

Resseguier et al. (2017b, 2020b) decompose the absolute diffusivity scale by scale:

$$a^R = \int_{0}^{+\infty} A_{v'}^R(\kappa) d\kappa. \qquad (3.9)$$

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) = \frac{1/\kappa}{v_g^0(\kappa)}$ and then

$$\frac{1}{2} A^R(\kappa) = \frac{1}{2} \tau^R(\kappa) E_{\kappa}(\kappa) = \frac{1/\kappa}{v_g^0(\kappa)} E_{\kappa}(\kappa), \qquad (3.10)$$

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

To calibrate an equivalent noise, we model $v'$ by $\sigma dB_t/dt$, where $dB_t/dt$ is a spatio-

temporal white noise and $\sigma$ denotes a spatial filtering operator which encodes spatial correlations through its ADSD, $A_{v'}^R$, and the horizontal incompressibility condition $(\nabla \cdot \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 $\hat{B}$ with a filter $\tilde{\psi}_{\sigma}$

(Resseguier et al. 2017b), that is $\tilde{\psi}_{\sigma} \star \hat{B}$ where $\star$ denotes a spatial convolution. The velocity

field is hence:

$$v' = \sigma dB_t/dt = \nabla^{\perp} \tilde{\psi}_{\sigma} \star dB_t/dt, \qquad (3.11)$$

with $\nabla^{\perp}$ the two-dimensional curl. That formula is easily written and implementable

in Fourier space (see equation (A 2)). To define the streamfunction filter, we note that

$$\frac{\pi \kappa^3}{(2\pi)^2} |\tilde{\psi}_{\sigma}(\kappa)|^2 = \frac{1}{2} \oint_{0}^{2\pi} d\theta_{\kappa} \kappa \frac{\|\sigma dB_t(\kappa)\|^2}{(2\pi)^2 dt} = A_{v'}^R(\kappa), \text{ i.e. the filter can be fully defined by the}$$

small-scale ADSD, $A_{v'}^R$. To close our model, we assume an ADSD power law:

$$A^R(\kappa) \approx A_0^R \kappa^{-\mu}. \qquad (3.12)$$

It enables automatic closure calibration: $A_{v'}^R(\kappa) = A_0^R \kappa^{-\mu} - A_{\bar{v}}^R(\kappa)$, from instantaneous

large-scale current statistics $A_{\bar{v}}^R$ only (Resseguier et al. 2020b) as illustrated in figure 1.

## 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 $v$ by $\bar{v} + \sigma \circ dB_t/dt$. Then, the stochastic transport of phase,

$\frac{d}{dt} \phi = \omega_0 (\|\nabla \phi\|)$, i.e. – up to that velocity replacement – the Stratonovich dispersion

relation is exactly (2.1). The method of characteristics also applies. Note, one can switch

from Stratonovich to Itō notations, where $v'$ corresponds to $\sigma dB_t/dt$. The characteristics
currents force the ray dynamics through a stochastic noise. For a single ray $(x_r, k) =$ $(x_r, y_r, k \cos \theta_k, k \sin \theta_k)$, this noise can be rigorously described by four independent one-dimensional white noises only (see Appendix), $\dot{B}_t^{(1)}, \dot{B}_t^{(2)}, \dot{B}_t^{(3)}$, and $\dot{B}_t^{(4)}$, and:

$$ \frac{d}{dt} x_r = v_g^0 \cos \theta_k + \bar{v} \cos \bar{\theta} + \sqrt{a_0} \dot{B}_t^{(1)}, \quad (4.3) $$

$$ \frac{d}{dt} y_r = v_g^0 \sin \theta_k + \bar{v} \sin \bar{\theta} + \sqrt{a_0} \dot{B}_t^{(2)}, \quad (4.4) $$

$$ \frac{d}{dt} \log k = -\bar{\sigma} \sin(\zeta) + \gamma_0 + \sqrt{\gamma_0} \dot{B}_t^{(3)}, \quad (4.5) $$

$$ \frac{d}{dt} \theta_k = \frac{1}{2}(\bar{\omega} - \bar{\sigma} \cos(\zeta)) + \sqrt{3\gamma_0} \dot{B}_t^{(4)}, \quad (4.6) $$

where $\zeta = 2(\theta_k + \bar{\phi})$ and

$$ a_0 = \frac{1}{2\bar{d}} \mathbb{E} \| \sigma dB_t \|^2 = \int_0^{+\infty} A_{\bar{v}}^R(k) dk, \quad (4.7) $$

$$ \gamma_0 = \frac{1}{8\bar{d}} \mathbb{E} \| \nabla_x (\sigma dB_t)^T \|^2 = \frac{1}{4} \int_0^{+\infty} k^2 A_{\bar{v}}^R(k) dk \quad (4.8) $$

Diffusivity constants depend through (3.9) on both the correlation length and the spectrum slope of the small-scale velocity. Contrasting 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.10). This isotropy assumption and its implication are discussed in Appendix. Note that neither large-scale nor small-scale component is assumed to be steady, even though the small-scale time variations are mainly driven by the wave speed and not by their Eulerian correlation time. The current unsteadiness can also lead to wavenumber variations (Dong et al. 2020; Cox et al. 2023; Boury et al. 2023). Given a known wavevector angle, it leads to a wavenumber evolution

$$ k(t) = k(0) \exp \left( - \int_{0}^{t} \bar{\sigma} \sin(2(\theta_k + \bar{\phi})) \right) \exp \left( \gamma_0 t + \sqrt{\gamma_0} B_t^{(3)} \right), \quad (4.9) $$

and hence 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 $v'$ tends to shorten the wavelength (Voronovich 1991; Boury et al. 2023) leading to this 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{v}$.

For completeness, the action distribution over space and wave vector can be derived. Some approaches consider finite-size wave trains either through additional equations (Jonsson 1990; White & Fornberg 1998) or re-meshing (Hell et al. in preparation). 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(x, k|x_r^0, k_r^0, t)$, the distribution of the ray $(x, k)$ at time $t$ given initial conditions $(x_r^0, k^0)$. Then, by analogy with tracers in incompressible turbulence (Piterbarg & Ostrovskii 1997, equation (1.31)) we can evaluate the wave action spectrum mean – or any point-wise statistics – as follows

$$ \mathbb{E}N(x, k, t) = \iint dx_r^0 dk^0 N^0(x_r^0, k_r^0) p(x, k|x_r^0, k^0, t) \quad (4.10) $$

where $N^0$ is the initial wave action spectrum. Integrating this expression over wavevectors,
<figure><img src="image_2.png" /><figcaption>Figure 1: KE spectrum ($m^2 \cdot s^{-2}/(\text{rad} \cdot m^{-1})$) (left) and ADSD ($m^2 \cdot s^{-1}/(\text{rad} \cdot m^{-1})$) (right) of the resolved high-resolution velocity, $A^R$, in red, low-resolution velocity, $A_v^R$, in blue, and modeled stochastic velocity, $A_{v'}^R(\kappa) = A_0^R \kappa^{-\mu} - A_v^R(\kappa)$, in green. For the ADSD power law, $A^R(\kappa) \approx A_0^R \kappa^{-\mu}$, we impose the theoretical KE spectrum slope $-5/3$ (black solid line), coherently with homogeneous SQG dynamics (see "Numerical results" section). The residual ADSD (green line) is set to extrapolate that power law at small scales.</figcaption></figure>

equations (2.2)-(2.3) also remain unchanged for homogeneous and isotropic $v'$:

$$ \left\{ \begin{aligned} \mathrm{d}\mathbf{x}_r &= (\mathbf{v}_g^0 + \bar{\mathbf{v}})\mathrm{d}t + \sigma \mathrm{d}B_t, \\ \mathrm{d}\mathbf{k} &= -\nabla(\bar{\mathbf{v}}\mathrm{d}t + \sigma \mathrm{d}B_t)^\mathrm{T}\mathbf{k}. \end{aligned} \right. \qquad (4.1) $$

## 4.1. Single-ray stochastic differential equations

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

Notably, dynamics of wave-vectors (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 et al. (1999) and references therein. Without loss of generality, the large-scale velocity can be parameterized as

$$ \bar{\mathbf{v}} = \bar{\mathbf{v}} \begin{pmatrix} \cos \bar{\theta} \\ \sin \bar{\theta} \end{pmatrix} \quad \text{and} \quad \nabla \bar{\mathbf{v}}^{\mathrm{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}. \qquad (4.2) $$

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

$$ \frac{\mathrm{d}}{\mathrm{d}t}x_r = v_g^0 \cos \theta_k + \bar{\mathbf{v}} \cos \bar{\theta} + \sqrt{a_0} \dot{\mathbf{B}}_t^{(1)}, \qquad (4.3) $$

$$ \frac{\mathrm{d}}{\mathrm{d}t}y_r = v_g^0 \sin \theta_k + \bar{\mathbf{v}} \sin \bar{\theta} + \sqrt{a_0} \dot{\mathbf{B}}_t^{(2)}, \qquad (4.4) $$

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

$$ \frac{\mathrm{d}}{\mathrm{d}t} \theta_k = \frac{1}{2}(\bar{\omega} - \bar{\sigma} \cos(\zeta)) + \sqrt{3\gamma_0 \dot{\mathbf{B}}_t^{(4)}}, \qquad (4.6) $$
we note that the wave action mean solely depends of group positions distribution

$$ \mathbb{E}A(\mathbf{x}, t) = \iint d\mathbf{x}_r^0 d\mathbf{k}^0 N^0(\mathbf{x}_r^0, \mathbf{k}_r^0) p(\mathbf{x}|\mathbf{x}_r^0, \mathbf{k}^0, t) \quad (4.11) $$

Multi-point action statistics – e.g. focusing $\mathbb{E}\|\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 wave vector.

## 4.2. Eulerian dynamics and action diffusion

Wave action spectrum is transported along a 4-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

$$ \begin{aligned} & \partial_t N + (\mathbf{v}_g^0 + \bar{\mathbf{v}} + \sigma \frac{dB_t}{dt}) \cdot \nabla_x N + \left(-\nabla_x (\bar{\mathbf{v}} + \sigma \frac{dB_t}{dt})^\mathrm{T} \mathbf{k}\right) \cdot \nabla_k N \\ & = \begin{bmatrix} \nabla_x \\ \nabla_k \end{bmatrix} \cdot \left( \mathbf{D} \begin{bmatrix} \nabla_x \\ \nabla_k \end{bmatrix} N \right) = \frac{1}{2} a_0 \Delta_x N + \frac{1}{2} \gamma_0 \frac{1}{k} \partial_k (k^3 \partial_k N) + \frac{3}{2} \gamma_0 \partial_{\theta_k}^2 N. \end{aligned} \quad (4.12) $$

The RHS is reminiscent to Eq. (3.16) in Bôas & Young (2020) and Eq. (36) in Smit & Janssen (2019), and more generally to rapid wave models. Nevertheless, equation (4.12) is not averaged and explicitly involves large-scale currents and noise terms (terms with factor $\frac{dB_t}{dt}$). Differences with Smit & Janssen (2019); Bôas & Young (2020) for the diffusivity estimates and the detailed computation of the 4 × 4 diffusion matrix D can be found in Appendix. Itō notations of (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{3}{2}\gamma_0$ act to diffuse wave action in space, wavenumber and wave-vector 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.10) 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:

$$ (\partial_t + \mathbf{v} \cdot \nabla) \Theta = 0 \quad \text{with} \quad \mathbf{v} = -\nabla^\perp (-\Delta)^{-\xi} \Theta. \quad (5.1) $$

Two extreme cases: the Surface Quasi-Geostrophic dynamics ($\xi = \frac{1}{2}$) (Pierrehumbert 1994; Lapeyre 2017), abbreviated SQG, and the two-dimensional Euler dynamics ($\xi = 1$), abbreviated 2D Euler. SQG is a dynamics with an extreme locality (KE spectrum slope -5/3) whereas 2D Euler has an extreme non-locality (KE spectrum slope -3). The objective is to test how the proposed closures apply for both dynamics – and possibly beyond validity of rapid wave models – to also be 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 medium. First, surface fast waves travel in a homogeneous and isotropic SQG turbulence, then waves propagate in a spatially localized 2D Euler turbulence, mimicking an oceanic jet. For both SQG and 2D Euler dynamics, a reference simulation is
<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic view of vectors and angles involved in single-ray dynamics. $\bar{S}_-$ and $\bar{S}_+$ are respectively compression and dilatation axes associated with the large-scale velocity gradient $\nabla \bar{v}^T$.</figcaption></figure>

where $\zeta = 2(\theta_k + \bar{\phi})$ and

$$a_0 = \frac{1}{2dt} \mathbb{E} \| \sigma dB_t \|^2 = \int_0^{+\infty} A_{v'}^R(\kappa) d\kappa, \qquad (4.7)$$

$$\gamma_0 = \frac{1}{8dt} \mathbb{E} \| \nabla_x (\sigma dB_t)^T \|^2 = \frac{1}{4} \int_0^{+\infty} k^2 A_{v'}^R(\kappa) d\kappa. \qquad (4.8)$$

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 large-scale nor 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 mainly driven by the large wave speed and not by the Eulerian velocity unsteadiness. The current unsteadiness can also lead to wavenumber variations (Dong et al. 2020; Cox et al. 2023; Boury et al. 2023). Given a known wavevector angle, it leads to a wavenumber evolution

$$k(t) = k(0) \exp \left( - \int_{0}^{t} \bar{\sigma} \sin(2(\theta_k + \bar{\phi})) dt' \right) \exp \left( \gamma_0 t + \sqrt{\gamma_0} B_t^{(3)} \right), \qquad (4.9)$$

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 $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{v}$.

For completeness, the action distribution over space and wave vector can be derived. Some approaches consider finite-size wave trains either through additional equations (Jonsson 1990; White & Fornberg 1998) or re-meshing (Hell et al. in preparation). 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(x, k|x_r^0, k_r^0, t)$, the distribution of the ray $(x, k)$ at time $t$ given initial conditions
<figure><img src="image_3.png" /><figcaption>Figure 2: Current velocity norm of the SQG homogeneous turbulence (left) and of the 2D Euler jet current at high-resolution (512 × 512).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Swell interacting with a high-resolution (512 × 512) deterministic SQG current (left), a low-resolution (32 × 32) deterministic SQG current (middle) and a low-resolution (32 × 32) deterministic SQG current plus (one realization of) the time-uncorrelated stochastic model (right) – colored by the corresponding wave amplitude (right-hand side colorbar) – computed by forward advection and superimposed on the current vorticity $\omega = \nabla^{\perp} \cdot \mathbf{v}$. The red cross indicate where the bidirectional wave spectra of figure 4 are computed.</figcaption></figure>

obtained at a resolution 512 × 512 for a 1000-km squared domain, through a pseudo-spectral code (Resseguier et al. 2017b, 2020b). Once initialized, the current velocity **v** is about 0.1 m.s⁻¹ for the homogeneous turbulence and 1 m.s⁻¹ for the jet (see figure 2).

## 5.2. Rays scattering in homogeneous SQG turbulence

A wave system enters the bottom boundary, propagating to the top. The carrier incident wave has an intrinsic wave group velocity of 10 m/s, i.e. wavelength $\lambda = 250$ m. Its envelope is Gaussian with an isotropic spatial extension of $30\lambda$. The left panels of figures 3 and 4 illustrates the resulting branched regime, spreading the wavevectors (figure 4) of the incoming waves. From bottom to top, spectral diffusion occurs (figure 4), in the direction orthogonal (here $k_x$) to the propagation (here $k_y$), in line with equation (4.6). This accelerates – along the propagation – the wave position spread (figure 3). This acceleration is explained by the ray equation (4.3) dominated by the intrinsic wave group velocity.

To mimic a badly resolved $\bar{v}$ field, **v** is smoothed at a resolution 32 × 32. Using this coarse-scale current, middle panels of figures 3 and 4, the branched regime is strongly depleting, i.e. the spectral small-scale turbulence diffusion is missing. A stochastic current $\mathbf{v}'$ is then
(x_r^0, k^0). Then, by analogy with tracers in incompressible turbulence (Piterbarg & Ostrovskii 1997, equation (1.31), see also Appendix D) we can evaluate the wave action spectrum mean – or any point-wise statistics – as follows

$$ \mathbb{E}N(\mathbf{x}, \mathbf{k}, t) = \iint \mathrm{d}\mathbf{x}_r^0 \mathrm{d}\mathbf{k}^0 N^0(\mathbf{x}_r^0, \mathbf{k}_r^0) p(\mathbf{x}, \mathbf{k} | \mathbf{x}_r^0, \mathbf{k}^0, t), \quad (4.10) $$

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

$$ \mathbb{E}A(\mathbf{x}, t) = \iint \mathrm{d}\mathbf{x}_r^0 \mathrm{d}\mathbf{k}^0 N^0(\mathbf{x}_r^0, \mathbf{k}_r^0) p(\mathbf{x} | \mathbf{x}_r^0, \mathbf{k}^0, t). \quad (4.11) $$

The wave action mean solely depends of group positions distribution. Multi-point action statistics – e.g. focusing $\mathbb{E}\|\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 wave vector.

## 4.2. Eulerian dynamics and action diffusion

Wave action spectrum is transported along a 4-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

$$ \begin{aligned} \partial_t N + (\nu_g^0 + \bar{\nu} + \sigma \frac{dB_t}{dt}) \cdot \nabla_x N + \left( -\nabla_x (\bar{\nu} + \sigma \frac{dB_t}{dt})^\top \mathbf{k} \right) \cdot \nabla_k N \\ = \begin{bmatrix} \nabla_x \\ \nabla_k \end{bmatrix} \cdot \begin{pmatrix} \mathbf{D} \\ \begin{bmatrix} \nabla_x \\ \nabla_k \end{bmatrix} \end{pmatrix} N = \frac{1}{2} a_0 \Delta_x N + \frac{1}{2} \gamma_0 \frac{1}{k} \partial_k (k^3 \partial_k N) + \frac{3}{2} \gamma_0 \partial_{\theta_k}^2 N. \end{aligned} \quad (4.12) $$

The RHS is reminiscent to Eq. (3.16) in Bôas & Young (2020) and Eq. (36) in Smit & Janssen (2019), and more generally to rapid wave models. Nevertheless, equation (4.12) is not averaged and explicitly involves large-scale currents and noise terms (terms with factor $\frac{dB_t}{dt}$). Differences with Smit & Janssen (2019); Bôas & Young (2020) for the diffusivity estimates and the detailed computation of the 4 × 4 diffusion matrix **D** can be found in Appendix A. Itō notations of (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{3}{2}\gamma_0$ act to diffuse wave action in space, wavenumber and wave-vector 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:

$$ (\partial_t + \nu \cdot \nabla) \Theta = 0 \quad \text{with} \quad \nu = -\nabla^{\perp}(-\Delta)^{-\xi} \Theta, \quad (5.1) $$

where $\Theta$ stands for the buoyancy, $\nabla^{\perp}$ the curl and $\Delta$ the Laplacian. Two extreme cases: the Surface Quasi-Geostrophic dynamics ($\xi = \frac{1}{2}$) (Held et al. 1995; Lapeyre 2017), abbreviated
<figure><img src="image_5.png" /><figcaption>Figure 4: Bidirectional wave spectra, computed by backward advection, at 8 locations along a vertical axis (the mean wave propagation direction) resulting from a swell interacting with a high-resolution (512 × 512) deterministic SQG current (left), a low-resolution (32 × 32) deterministic SQG current (middle) and a low-resolution (32 × 32) deterministic SQG current plus (one realization of) the stochastic model (3.10) (right). The spatial locations where the spectra are calculated are specified on top of each image. We point out these locations on figure 3 by the red crosses.</figcaption></figure>

300 added for ray tracing (4.1). This stochastic component is divergence-free and has a self-301 similar distribution of energy across spatial scales (3.10) (see figure 1). The resulting spatial 302 and spectral spreads are now comparable to simulations with high-resolution currents. For 303 this setting, the stochastic closure provides satisfying results for a sufficiently well-resolved 304 large-scale current. The key decorrelation ratio $\epsilon = \frac{l_{v'} \|v\|}{l_v v_g^0}$ indeed depends on the resolution 305 through $l_{v'}$. The large-scale current $\bar{v}$ is resolved on a 32×32 grid, i.e. with a resolution
<figure><img src="image_4.png" /><figcaption>Figure 3: Current velocity norm of the SQG homogeneous turbulence (left) and of the 2D Euler jet current at high-resolution (512 × 512).</figcaption></figure>

SQG, and the two-dimensional Euler dynamics ($\xi = 1$), abbreviated 2D Euler. SQG is a
dynamics with an extreme locality (KE spectrum slope $-5/3$) whereas 2D Euler has an
extreme non-locality (KE 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 medium. 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 2D Euler turbulence, mimicking an oceanic jet. For
both SQG and 2D Euler dynamics, a reference simulation is obtained at a resolution 512×512
for a 1000-km squared domain, with the help of a pseudo-spectral code (Resseguier et al.
2017b, 2020b). Once initialized, the current velocity **v** is about $0.1 \ \text{m.s}^{-1}$ for the homogeneous
turbulence and $1 \ \text{m.s}^{-1}$ for the jet (see figure 3).

5.2. *Rays scattering in homogeneous SQG turbulence*

A wave system enters the bottom boundary, propagating to the top. The carrier incident wave
has an intrinsic wave group velocity of 10 m/s, i.e. wavelength $\lambda = 250$ m. Its envelope
is Gaussian with an isotropic spatial extension of $30\lambda$. The left panels of figures 4 and
5 illustrates 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 equation
(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.

To mimic a badly resolved $\bar{v}$ field, $v$ is smoothed at a resolution 32 × 32. Using this
coarse-scale current, middle panels of figures 4 and 5, 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 $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
306 $l_{v'} = \frac{\|\nabla v^T\|}{\|\nabla v^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
307 $\|v\| \approx 0.12$ m.s$^{-1}$, so $\frac{\|v\|}{v_g^0} \approx 1.2 \times 10^{-2}$, which is sufficiently small to make the proposed
308 model applicable.

309 From the ADSD estimate (equation (3.9) illustrated by figure 1) and (4.7)-(4.8), evaluations
310 of the diffusivity coefficients $a_0$ and $\gamma_0$ are straightforward. Previously discussed Smit &
311 Janssen (2019), the spatial diffusivity is extremely weak: $a_0 = 6.4 \times 10^{-1}$ m$^2$.s$^{-1}$ (spatial
312 variations in ray equations (4.3)-(4.4) of about $\sqrt{a_0}t = 230$ m during 1 day). In contrast, the
313 spectral angle diffusivity is large: $3\gamma_0 = 3.0 \times 10^{-8}$ rad$^2$.s$^{-1}$. Along our 1-day simulation,
314 neglecting large-scale velocity influence, (4.6) leads to a Brownian wave vector angle
315 variations $\delta\theta_k = \theta_k - \theta_k(0) = \sqrt{3\gamma_0}B_t^{(4)}$ with a standard deviation $\sigma_{\delta\theta_k} = \sqrt{3\gamma_0}t =
316 5.2 \times 10^{-2}$ rad $\approx 3.0^\circ$, eventually increasing the wave group spectral maximal extension from
317 $\pm 2\sigma_{k_x} = \pm 2\frac{2\pi}{30\lambda} = \pm 1.7 \times 10^{-3}$ rad.m$^{-1}$ to $\pm 2\sigma_{k_x} \approx \pm 2\sqrt{\left(\frac{2\pi}{30\lambda}\right)^2 + (k\sigma_{\delta\theta_k})^2} = \pm 3.1 \times 10^{-3}$
318 rad.m$^{-1}$, confirmed by figure 4. This figure also illustrates the wave action diffusion induced
319 by diffusivity $\gamma_0$, well predicted by the Eulerian wave action model (4.12). In the branched
320 regime, the increased angle variability leads, by advection, to a spatial spread. The simplified
321 ray equation (4.3) gives $\delta x \approx \int_0^t v_g^0 \cos \theta_k \approx v_g^0 \int_0^t \delta \theta_k \approx v_g^0 \sqrt{3\gamma_0} \int_0^t B_{t'}^{(4)}$ with maximal
322 extension $\pm 2\sigma_x \approx \pm 2v_g^0 \sqrt{\gamma_0}t^3 \approx \pm 52$ km, in agreement with figure 3. Finally, we estimate a
323 first-order delay along the propagation $\delta t = t - (y-y(0))/v_g^0 \approx \int_0^t (1-\sin \theta_k) \approx \frac{1}{2} \int_0^t \delta \theta_k^2 \approx
324 \frac{3}{2}\gamma_0 \int_0^t (B_{t'}^{(4)})^2$, with mean value $\mathbb{E}\delta t = \frac{3}{4}\gamma_0 t^2$.

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

325 Tests are now performed for rays traveling in fast and strongly heterogeneous 2D Euler flows.
326 Classical fast wave models – assuming flows of weak amplitude and often uniform statistics –
327 are expected to fail here. Jets exhibit strong current gradients (e.g. Kudryavtsev et al. 2017)
328 creating strong ray focusing and possibly rogue events. Passing through localized spatial
329 structures, caustics can appear but solely from unrealistically collimated wave trains (White
330 & Fornberg 1998; Heller et al. 2008; Wang et al. 2023). Occurrences strongly reduce for
331 finite directional spread (Slunyaev & Shrira 2023). Here, wave groups are trapped in a jet, but
332 nonlinear wave interactions are neglected. The high-resolution numerical simulations (see
333 figure 5) reveal that even linear wave trains are well trapped in adversarial currents. Freund
334 & Fleischman (2002) observed a similar behavior for acoustic waves in a 3D turbulent jet.
335 Note that during our simulation, rays cross the domain several times (because of the doubly
336 periodic boundary conditions, see appendix for technical details). Top (resp. bottom) of the
337 jet, the vorticity and thus – at first order – rays curvatures (Dysthe 2001) are negative (resp.
338 positve). Therefore, rays oscillate around the jet. A toy model can explain this behavior.
339 Following the multiscale stochastic approach (4.3)-(4.6), wave scattering is also taken into
340 account.
341

342 For very-coarsed-grained (4 × 4) current $\bar{v}$, oscillation remains, but most of the scattering
343 vanishes, as illustrated by figure 6. Moreover, the curvature of the jet creates artificial wave
344 focusing at $t = 8$ and 10 days. Introducing a time-uncorrelated model (3.10) corrects the
345 resolution issue on figure 7. Figure 8 plots the current ADSD. The current is strong ($\|v\| \approx 1.4$
346 m.s$^{-1}$), and the usual fast wave approximation cannot be applied ($\frac{\|v\|}{v_g^0} \approx 1.2 \times 10^{-1}$). However,
347 the proposed modified fast wave model is valid, even at the very coarse 4 × 4 resolution.
348 Indeed, 2D Euler spectra are steeper than for SQG dynamics, and the length scale ratio is
<figure><img src="image_5.png" /><figcaption>Figure 4: Swell (wavelength $\lambda$ = 250 m) interacting with a high-resolution (512 $\times$ 512) deterministic SQG current (left), a low-resolution (32 $\times$ 32) deterministic SQG current (middle) and a low-resolution (32 $\times$ 32) deterministic SQG current plus (one realization of) the time-uncorrelated stochastic model (right) – colored by the corresponding wave amplitude, $h(t) = \sqrt{\omega_0(k(t))N(t=0)}$ (right-hand side colorbar) – computed by forward advection and superimposed on the current vorticity $\omega = \nabla^{\perp} \cdot \mathbf{v}$. The red cross indicate where the bidirectional wave spectra of figure 5 are computed.</figcaption></figure>

currents. For this setting, the stochastic closure provides satisfying results for a sufficiently well-resolved large-scale current. The key decorrelation ratio $\epsilon = \frac{l_{v'}}{l_v} \frac{\|\mathbf{v}\|}{v_g^0}$ indeed depends on the resolution through $l_{v'}$. The large-scale current $\bar{\mathbf{v}}$ is resolved on a 32×32 grid, i.e. with a resolution $l_{v'} = \frac{\|\nabla \mathbf{v}'\|}{\|\nabla \mathbf{v}'\|} l_v = 0.33 l_v$. As such $\epsilon = 4.1 \times 10^{-3}$, computed with $v_g^0 \approx 10 \text{ m.s}^{-1}$ and $\|\mathbf{v}\| \approx 0.12 \text{ m.s}^{-1}$, so $\frac{\|\mathbf{v}\|}{v_g^0} \approx 1.2 \times 10^{-2}$, which is sufficiently small to make the proposed model applicable.

From the ADSD estimate (equation (3.10) illustrated by figure 1) and (4.7)-(4.8), evaluations of the diffusivity coefficients $a_0$ and $\gamma_0$ are straightforward. Previously discussed Smit & Janssen (2019), the spatial diffusivity is extremely weak: $a_0 = 6.4 \times 10^{-1} \text{ m}^2.\text{s}^{-1}$ (spatial variations in ray equations (4.3)-(4.4) of about $\sqrt{a_0 t} = 230 \text{ m}$ during 1 day). In contrast, the spectral angle diffusivity is large: $3\gamma_0 = 3.0 \times 10^{-8} \text{ rad}^2.\text{s}^{-1}$. Along our 1-day simulation, neglecting large-scale velocity influence, (4.6) leads to a Brownian wave vector angle variations $\delta\theta_k = \theta_k - \theta_k(0) = \sqrt{3\gamma_0}B_t^{(4)}$ with a standard deviation $\sigma_{\delta\theta_k} = \sqrt{3\gamma_0 t} = 5.2 \times 10^{-2} \text{ rad} \approx 3.0^\circ$, eventually increasing the wave group spectral maximal extension from $\pm 2\sigma_{k_x} = \pm 2\frac{2\pi}{30\lambda} = \pm 1.7 \times 10^{-3} \text{ rad.m}^{-1}$ to $\pm 2\sigma_{k_x} \approx \pm 2\sqrt{\left(\frac{2\pi}{30\lambda}\right)^2 + \left(k\sigma_{\delta\theta_k}\right)^2} = \pm 3.1 \times 10^{-3} \text{ 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 dt' \approx v_g^0 \int_0^t \delta\theta_k dt' \approx v_g^0 \sqrt{3\gamma_0} \int_0^t B_{t'}^{(4)} dt'$ with maximal extension $\pm 2\sigma_x \approx \pm 2v_g^0 \sqrt{\gamma_0 t^3} \approx \pm 52 \text{ km}$, in agreement with figure 4. Finally, we estimate a first-order delay along the propagation $\delta t = t - (y-y(0))/v_g^0 \approx \int_0^t (1-\sin \theta_k) dt' \approx \frac{1}{2} \int_0^t \delta\theta_k^2 dt' \approx \frac{3}{2}\gamma_0 \int_0^t (B_{t'}^{(4)})^2 dt'$, with mean value $\mathbb{E}\delta t = \frac{3}{4}\gamma_0 t^2$.
<figure><img src="image_6.png" /><figcaption>Figure 5: Rays facing a high-resolution (512 × 512) deterministic 2D Euler jet current – colored by the corresponding wave amplitude (right-hand side colorbar) – computed by forward advection and superimposed on the current vorticity $\omega = \nabla^{\perp} \cdot \mathbf{v}$.</figcaption></figure>

349 already significant at this resolution, $\frac{l_v'}{l_v} = 0.14$, and the derived time-decorrelation ratio is

small: $\epsilon = \frac{l_v'}{l_v} \frac{\|\mathbf{v}\|}{\bar{v}_g^0} = 1.6 \times 10^{-2}$.

351 Furthermore, by approximating the under-resolved current $\bar{v}$, an analytic stochastic solution

352 can be obtained for a ray traveling against the current. The large-scale pattern of the jet takes
<figure><img src="image_6.png" /><figcaption>Figure 5: Bidirectional wave spectra, computed by backward advection, at 8 locations along a vertical axis (the mean wave propagation direction) resulting from a swell interacting with a high-resolution (512 × 512) deterministic SQG current (panel (a)), a low-resolution (32 × 32) deterministic SQG current (panel (b)) and a low-resolution (32 × 32) deterministic SQG current plus (one realization of) the stochastic model (3.11) (panel (c)). The spatial locations where the spectra are calculated are highlighted on figure 4 by the red crosses.</figcaption></figure>

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

Tests are now performed for rays traveling in fast and strongly heterogeneous 2D 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 et al. 2008; Wang et al. 2023). Occurrences strongly reduce for

365

366

367

368

369

370

371
<figure><img src="image_7.png" /><figcaption>Figure 6: Rays facing a low-resolution (4 × 4) deterministic 2D Euler jet current – colored by the corresponding wave amplitude (right-hand side colorbar) – computed by forward advection and superimposed on the current vorticity $\omega = \nabla^{\perp} \cdot \mathbf{v}$.</figcaption></figure>

353 a quadratic form

$$ \overline{u} \approx \overline{U}_0 - \frac{1}{2} \overline{\beta} \left( y - \frac{L_y}{2} \right)^2 \quad \text{and} \quad \overline{v} \approx 0, \quad \text{with} \quad \overline{U}_0, \overline{\beta} < 0 \qquad (5.2) $$

355 Note, the toy model (5.2) simply considers a straight jet, neglecting its curvature. For weak 356 subgrid currents and a ray, $(x_r, y'_r + \frac{L_y}{2}, k, \theta_k)$, propagating mainly to the right, $\theta_k$ is small
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 behavior for acoustic waves in a 3D 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). 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 behavior. Following the multiscale stochastic approach (4.3)-(4.6), wave scattering is also taken into account.

For very-coarsed-grained (4 × 4) current $\bar{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 on figure 8. Figure 9 plots the current ADSD. The current is strong ($\|\bar{v}\| \approx 1.4$ m.s$^{-1}$), and the usual fast wave approximation cannot be applied ($\frac{\|\bar{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 × 4 resolution. Indeed, 2D Euler spectra are steeper than for SQG dynamics, and the length scale ratio is already significant at this resolution, $\frac{l_{v'}}{l_v} = 0.14$, and the derived time-correlation ratio is small: $\epsilon = \frac{l_{v'}}{l_v} \frac{\|\bar{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 traveling against the current. The large-scale pattern of the jet takes a quadratic form

$$ \bar{u} \approx \bar{U}_0 - \frac{1}{2}\bar{\beta} \left(y - \frac{L_y}{2}\right)^2 \quad \text{and} \quad \bar{v} \approx 0, \quad \text{with} \quad \bar{U}_0, \bar{\beta} < 0. \qquad (5.2) $$

Note, the toy model (5.2) simply considers a straight jet, neglecting its curvature. For weak subgrid currents and a ray, $(x_r, y'_r + \frac{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$.

$$ \frac{d}{dt} y'_r \approx v_g^0 \sin(\theta_k) = v_g^0 \theta_k + O(\theta_k^2). \qquad (5.3) $$

For frozen turbulence, the wave-number and hence $v_g^0$ will not significantly vary. 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)$ term from now on. Moreover, $\tilde{k}^{\perp} \cdot \nabla \tilde{v}^{\top} \tilde{k} \approx -\partial_y \bar{u}$ and the dynamics of wave vector angle (4.6) simplifies to a stochastic oscillator equation:

$$ \frac{d^2}{dt^2} y'_r = v_g^0 \frac{d}{dt} \theta_k = -\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)}. \qquad (5.4) $$

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. Solution of this linear equation is known (e.g. Resseguier et al. 2017a, Eq.(51)-(55)):

$$ y_r(t) = \underbrace{\frac{L_y}{2} + y'_r(0) \cos(\bar{\omega}_r t) + \frac{v_g^0}{\bar{\omega}_r} \theta_k(0) \sin(\bar{\omega}_r t)}_{=\mathbb{E}(y_r(t))} + \underbrace{Y_{\gamma_0} \sqrt{\bar{\omega}_r} \int_0^t \sin(\bar{\omega}_r (t-r)) dB_r^{(4)}}_{=y_r''(t)}, \qquad (5.5) $$

with $Y_{\gamma_0} = v_g^0 \sqrt{3\gamma_0/\bar{\omega}_r^3}$. The wavevector angle solution is similar. The solution ensemble
<figure><img src="image_8.png" /><figcaption>Figure 7: Rays facing a low-resolution (4 × 4) deterministic 2D Euler jet current plus (one realization of) the time-uncorrelated stochastic model – colored by the corresponding wave amplitude (right-hand side colorbar) – computed by forward advection and superimposed on the current vorticity $\omega = \nabla^{\perp} \cdot \mathbf{v}$.</figcaption></figure>

357 and the simplified ray equation (4.4) determines the group position with respect to the jet $y'_r$

$$ \frac{d}{dt} y'_r \approx v_g^0 \sin(\theta_k) = v_g^0 \theta_k + O(\theta_k^2). \qquad (5.3) $$

359 For frozen turbulence, the wave-number and hence $v_g^0$ will not significantly vary. The other
<figure><img src="image_7.png" /><figcaption>Figure 6: Rays facing a high-resolution (512 × 512) deterministic 2D Euler jet current – colored by the corresponding wave amplitude, $h(t) = \sqrt{\omega_0(k(t))}N(t=0)$ (right-hand side colorbar) – computed by forward advection and superimposed on the current vorticity $\omega = \nabla^{\perp} \cdot \mathbf{v}$ (top colorbar).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: ADSD ($m^2 \cdot s^{-1} / (\text{rad} \cdot m^{-1})$) of the resolved high-resolution jet velocity in brown, low-resolution jet velocity in blue, and modeled stochastic velocity, in red. The theoretical spectrum slope -3 (black solid line) is imposed, consistent with homogeneous 2D Euler dynamics. The residual ADSD (red line) is set to extrapolate that power law at small scales.</figcaption></figure>

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{k}^{\perp} \cdot \nabla \tilde{v}^{\perp} \tilde{k} \approx -\partial_y \bar{u}$ and the dynamics of wave vector angle (4.6) simplifies to a stochastic oscillator equation:

$$ \frac{d^2}{dt^2} y_r' = v_g^0 \frac{d}{dt} \theta_k = -\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)}. \quad (5.4) $$

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. Solution of this linear equation is known (e.g. Resseguier et al. 2017a, Eq.(51)-(55)):

$$ y_r(t) = \underbrace{\frac{L_y}{2} + y_r'(0) \cos(\bar{\omega}_r t) + \frac{v_g^0}{\bar{\omega}_r} \theta_k(0) \sin(\bar{\omega}_r t)}_{=\mathbb{E}(y_r(t))} + Y_{\gamma_0} \sqrt{\bar{\omega}_r} \int_0^T \sin(\bar{\omega}_r (t-r)) dB_r^{(4)}, \quad (5.5) $$

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 6. From the coarse-scale vorticity shear plotted in figure 9 in the vicinity of the jet, we can estimate $\bar{\beta} = -2.7 \times 10^{-11} \text{m}^{-1} \cdot \text{s}^{-1}$. It yields an oscillation frequency $\bar{\omega}_r = 1.3 \times 10^{-5} \text{rad} \cdot \text{s}^{-1}$ i.e. a period of $2\pi / \bar{\omega}_r = 5.7$ days, in agreement with the ray tracing simulations. Note that the high-resolution vorticity shear, left panel of figure 9, 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 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) are entirely defined by their mean and covariance functions. Specifically,

$$ 384 \mathbb{E} (y_r''(t) y_r''(t+\tau)) = \frac{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))). \quad (5.6) $$
<figure><img src="image_8.png" /></figure>

Figure 7: Rays facing a low-resolution (4 × 4) deterministic 2D Euler jet current – colored by the corresponding wave amplitude, $h(t) = \sqrt{\omega_0(k(t))}N(t=0)$ (right-hand side colorbar) – computed by forward advection and superimposed on the current vorticity $\omega = \nabla^{\perp} \cdot \mathbf{v}$ (top colorbar).
<figure><img src="image_10.png" /><figcaption>Figure 9: Vorticity shear $\partial_y^2 u$ of the deterministic 2D Euler jet current at high-resolution (512 × 512) (left), at low-resolution (4 × 4) (middle), and the corresponding swell system period $2\pi/\bar{\omega}_r$. Far from the jet (±200 km away), the vorticity shear becomes zero or even positive, so period larger than 10 days are cropped.</figcaption></figure>

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\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{\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 branched regime with a much faster vertical location spreading: $\pm 2\sigma_y = \pm 2\sqrt{(2\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 $t^3$ scaling of the branched regime.

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) \approx 2 \int_0^t \bar{\omega} \theta_k = 2\bar{\beta} \int_0^t y'_r \theta_k$ 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)$ oscillates at frequency $2\bar{\omega}_r$ and the oscillations modulate the wave amplitude: $h = \sqrt{E} = \sqrt{\omega_0 N} = \text{cst.} k^{1/4}$. The modulations are associated with wave-current energy exchanges Boury et al. (2023), visible in the colored rays, figures 5, 6 and 7, when the groups enter and exit the jet.

Finally, the conditional ray distribution, $p(\mathbf{x}, \mathbf{k} | \mathbf{x}_r^0, \mathbf{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

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

with $\mathcal{N}(\bullet|\tilde{\sigma}_y^2(t))$, a Gaussian function with variance $\tilde{\sigma}_y^2(t) = \sigma_y^2(t) + \left(\frac{v_g^0}{\bar{\omega}_r} \sin(\bar{\omega}_r t) \sigma_{\delta\theta_k}^0\right)^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 an 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
<figure><img src="image_9.png" /></figure>

Figure 8: Rays facing a low-resolution (4 × 4) deterministic 2D Euler jet current plus (one realization of) the time-uncorrelated stochastic model – colored by the corresponding wave amplitude, $h(t) = \sqrt{\omega_0(k(t))N(t=0)}$ (right-hand side colorbar) – computed by forward advection and superimposed on the low-resolution current vorticity $\overline{\omega} = \nabla^\perp \cdot \overline{v}$ (top colorbar).
413 decorrelate faster, leading to a broader validity range of fast wave approximations. The
414 proposed framework encodes both large-scale refraction and random scattering effects on
415 wave statistical properties. The mean wave-action statistics are directly linked to resolved
416 strain-rate and vorticity, but also to unresolved KE spectral properties. Both Eulerian and
417 Lagrangian views are presented. A convenient calibration method is also proposed for the
418 subgrid parametrization.

419 Anticipated, random horizontal currents delay wave arrival and augment the initial
420 radiative transport equation with a directional diffusive term. These phenomena are illustrated
421 with numerical simulations, analytical solutions, and quantitative proxies describing weak
422 homogeneous turbulence. Using these proxies, measured delays in ray arrivals, estimated
423 wave energy spectral characteristics and decays, and/or varying directional spread shall then
424 be used to more quantitatively interpreted the upper the turbulent underlying flow properties.

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

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

439 **Funding.** This work is supported by the R&T CNES R-S19/OT-0003-084, the ERC project 856408-STUOD,
440 the European Space Agency World Ocean Current project (ESA Contract No. 4000130730/20/I-NB), and
441 SCALIAN DS.

442 Declaration of interests. The authors report no conflict of interest.

443 **Data availability statement.** SCALIAN DS owns a portion of the developed code intellectual property.
444 For commercial reasons, that code will remains private.

445 Author ORCIDs. V. Resseguier, https://orcid.org/0000-0002-9301-9493; E. Hascoët https://orcid.org/0009-0002-0760-6534; B. Chapron, https://orcid.org/0000-0001-6088-8775

447 Author contributions. VR developed the theory. VR and EH wrote the code and performed numerical
448 experiments. VR and BC wrote the paper.

449 Appendix A. Stochastic forcing covariance

450 In this appendix, we will compute the conditionnal covariance of the stochastic forcing of
451 our eikonal characteristic equations (4.1), that is:

$$
2\mathbf{D} \triangleq \frac{1}{dt} \mathbb{E}_t \left\{ \begin{pmatrix} \sigma \mathrm{d}B_t \\ \mathrm{d}\eta_t \end{pmatrix} \begin{pmatrix} \sigma \mathrm{d}B_t \\ \mathrm{d}\eta_t \end{pmatrix}^\top \right\} = \begin{bmatrix} \mathbf{a} & \Sigma_{\eta, \sigma} \\ \Sigma_{\eta, \sigma}^\top & \Sigma_{\eta} \end{bmatrix} \quad (A1)
$$

453 where, $d\eta_t = -\nabla(\sigma dB_t)^\top k$, denotes the wave-vector stochastic forcing and, $\Sigma_\eta dt$, its
454 covariance, and $\mathbb{E}_t\{\cdot\}$ stands for the conditional expectation evaluated
455 with given characteristics $(x_r(t), k(t))$ at the current time $t$. Note that in the appendix we
456 use $\text{Itô}$ notations only.
<figure><img src="image_10.png" /><figcaption>Figure 9: ADSD ($m^2.s^{-1}/(rad.m^{-1})$) of the resolved high-resolution jet velocity in red, low-resolution jet velocity in blue, and modeled stochastic velocity, in green. The theoretical spectrum slope -3 (black solid line) is imposed, consistent with homogeneous 2D Euler dynamics. The residual ADSD (green line) is set to extrapolate that power law at small scales.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Vorticity shear $\partial_z^2 u$ of the deterministic 2D Euler jet current at high-resolution (512 x 512) (left), at low-resolution (4 x 4) (middle), and the corresponding swell system period $2\pi/\bar{\omega}_r$. Far from the jet (±200 km away), the vorticity shear becomes zero or even positive, so period larger than 10 days are cropped.</figcaption></figure>

mean, $Ey_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} \text{m}^{-1} \text{s}^{-1}$. It yields an oscillation frequency $\bar{\omega}_r = 1.3 \times 10^{-5} \text{rad.s}^{-1}$ i.e. a period of $2\pi/\bar{\omega}_r = 5.7$ days, in agreement with the ray tracing simulations. Note that the high-resolution vorticity shear, left panel of figure 10, 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 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) are entirely defined by their mean and covariance functions. Specifically,

$$424 \mathbb{E} (y_r''(t)y_r''(t+\tau)) = \frac{1}{4} Y_{\gamma_0}^2 \left( \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)) \right). \quad (5.6)$$

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
The subgrid velocity, $v' = \sigma dB_t/dt$ is constructed in Fourier space with a divergence-free isotropic spatial filter $\nabla^\perp \psi_\sigma$ (see (3.10)).

$$ \hat{v}'(\boldsymbol{\kappa}, t) = \int d\mathbf{x} \, v'(\mathbf{x}, t) e^{-i\boldsymbol{\kappa} \cdot \mathbf{x}} = \overline{\sigma dB_t/dt}(\boldsymbol{\kappa}, t) = i\boldsymbol{\kappa}^{\perp} \hat{\psi}_{\sigma}(\boldsymbol{\kappa}) \overline{dB_t}(\boldsymbol{\kappa})/dt, \quad (A2) $$

The computation of the variance tensor $\mathbf{a}$ is classical and straightforward from the definition of the inverse Fourier transform and the identity $\mathbb{E}\{\overline{dB_t}(\boldsymbol{\kappa}_1)\overline{dB_t}(\boldsymbol{\kappa}_2)\} = (2\pi)^2\delta(\boldsymbol{\kappa}_1 - \boldsymbol{\kappa}_2)dt$.

We simply need to split the integral of the stochastic forcing spectrum over the current wavevector $\boldsymbol{\kappa} = \kappa(\cos \theta_\kappa, \sin \theta_\kappa)$:

$$ \mathbf{a} = \frac{1}{(2\pi)^4 dt} \iint d\boldsymbol{\kappa}_1 d\boldsymbol{\kappa}_1 \mathbb{E}_t \left\{ \overline{\sigma dB_t(\boldsymbol{\kappa}_1, \mathbf{k})} \overline{\sigma dB_t^\top(\boldsymbol{\kappa}_2, \mathbf{k})} \right\} e^{i(\boldsymbol{\kappa}_1 - \boldsymbol{\kappa}_2) \cdot \mathbf{x}}, \quad (A3) $$

$$ = \frac{1}{(2\pi)^2} \int d\boldsymbol{\kappa} \, \kappa^2 |\hat{\psi}_{\sigma}(\boldsymbol{\kappa})|^2 \begin{pmatrix} -\sin \theta_{\boldsymbol{\kappa}} \\ \cos \theta_{\boldsymbol{\kappa}} \end{pmatrix} \begin{pmatrix} -\sin \theta_{\boldsymbol{\kappa}} \\ \cos \theta_{\boldsymbol{\kappa}} \end{pmatrix}^{\top}, \quad (A4) $$

$$ = \frac{1}{(2\pi)^2} \int_0^\infty \oint_0^{2\pi} d\boldsymbol{\kappa} d\theta_{\boldsymbol{\kappa}} \kappa^3 |\hat{\psi}_{\sigma}(\boldsymbol{\kappa})|^2 \begin{pmatrix} \sin^2 \theta_{\boldsymbol{\kappa}} & -\sin \theta_{\boldsymbol{\kappa}} \cos \theta_{\boldsymbol{\kappa}} \\ -\sin \theta_{\boldsymbol{\kappa}} \cos \theta_{\boldsymbol{\kappa}} & \cos^2 \theta_{\boldsymbol{\kappa}} \end{pmatrix}, \quad (A5) $$

$$ = \frac{2}{2\pi} a_0 \oint_0^{2\pi} d\theta_{\boldsymbol{\kappa}} \begin{pmatrix} \sin^2 \theta_{\boldsymbol{\kappa}} & -\sin \theta_{\boldsymbol{\kappa}} \cos \theta_{\boldsymbol{\kappa}} \\ -\sin \theta_{\boldsymbol{\kappa}} \cos \theta_{\boldsymbol{\kappa}} & \cos^2 \theta_{\boldsymbol{\kappa}} \end{pmatrix}, \quad (A6) $$

$$ = a_0 \mathbb{I}_d. \quad (A7) $$

Now, the Fourier transform of the wave-vector stochastic forcing is

$$ d\hat{\eta}_t = -\nabla(\sigma dB_t)^\top \mathbf{k} = -i\boldsymbol{\kappa}(i\boldsymbol{\kappa}^{\perp}\hat{\psi}_{\sigma}\overline{dB_t}) \cdot \mathbf{k} = -\boldsymbol{\kappa}(\mathbf{k}^{\perp} \cdot \boldsymbol{\kappa})\hat{\psi}_{\sigma}\overline{dB_t}. \quad (A8) $$

Then, applying the crest-oriented rotation matrix, $M_k = [\tilde{\mathbf{k}} \ \tilde{\mathbf{k}}^{\perp}]$, leads to

$$ d\hat{Z}_t = M_k^\top d\hat{\eta}_t = - \begin{pmatrix} \tilde{\mathbf{k}} \cdot \boldsymbol{\kappa} \\ \tilde{\mathbf{k}}^{\perp} \cdot \boldsymbol{\kappa} \end{pmatrix} (\mathbf{k}^{\perp} \cdot \boldsymbol{\kappa}) \hat{\psi}_{\sigma} \overline{dB_t} = - \begin{pmatrix} \cos \delta\theta \sin \delta\theta \\ \sin^2 \delta\theta \end{pmatrix} \kappa^2 k \hat{\psi}_{\sigma} \overline{dB_t}, \quad (A9) $$

with $\delta\theta = \theta_\kappa - \theta_k$. From there, we can evaluate the conditional covariance matrix $\Sigma_Z = \frac{1}{dt} \mathbb{E}_t\{dZ_t dZ_t^\top\}$ of $dZ_t$ as before:

$$ \Sigma_Z = \frac{1}{(2\pi)^4 dt} \iint d\boldsymbol{\kappa}_1 d\boldsymbol{\kappa}_1 \mathbb{E}_t \left\{ d\hat{Z}_t(\boldsymbol{\kappa}_1, \mathbf{k}) d\overline{\hat{Z}_t^\top}(\boldsymbol{\kappa}_2, \mathbf{k}) \right\} e^{i(\boldsymbol{\kappa}_1 - \boldsymbol{\kappa}_2) \cdot \mathbf{x}}, \quad (A10) $$

$$ = \frac{1}{(2\pi)^2} \int_0^\infty \oint_0^{2\pi} d\boldsymbol{\kappa} d\delta\theta \, \kappa^5 k^2 |\hat{\psi}_{\sigma}(\boldsymbol{\kappa})|^2 \begin{pmatrix} \cos^2 \delta\theta \sin^2 \delta\theta & \cos \delta\theta \sin^3 \delta\theta \\ \cos \delta\theta \sin^3 \delta\theta & \sin^4 \delta\theta \end{pmatrix}, \quad (A11) $$

$$ = \gamma_0 k^2 \begin{bmatrix} 1 & 0 \\ 0 & 3 \end{bmatrix}. \quad (A12) $$

Finally, we come back to the canonical frame to get $\Sigma_\eta$

$$ \Sigma_\eta = \mathbb{E}_t \{d\eta_t d\eta_t^\top\} = M_k \Sigma_Z M_k^\top = \gamma_0 k^2 \left[ \tilde{\mathbf{k}} \tilde{\mathbf{k}}^\top + 3 \tilde{\mathbf{k}}^{\perp} (\tilde{\mathbf{k}}^{\perp})^\top \right] \quad (A13) $$

For noises cross-correlations, by isotropy, it is also straightforward to show that

$$ \Sigma_{\eta, \sigma} = 0. \quad (A14) $$

The stochastic forcings of $\mathbf{x}_r$ and $\mathbf{k}$ are hence (conditionally) independent from one another.
t = 2π/ωr, the group has oscillated once around the jet and the maximal position extension reaches ±2σy = ±2√π Yγ₀ = ±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: ±2σy = ±2√(2π)³/³ Yγ₀ = ±217 km. For large time, our multiscale approach predicts a scaling in t, much slower than the usual scattering t³ scaling.

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) dt' \approx 2 \int_0^t \bar{\omega} \theta_k dt' = 2\bar{\beta} \int_0^t y'_r \theta_k dt'$ and (4.9) together with the analytic solutions for y'_r and θ_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 dt')$ oscillates at frequency $2\bar{\omega}_r$ and the oscillations modulate the wave amplitude: $h = \sqrt{E} = \sqrt{\omega_0 N} = \text{cst.} k^{1/4}$. The modulations are associated with wave-current energy exchanges Boury et al. (2023), visible in the colored rays, figures 6, 7 and 8, when the groups enter and exit the jet.

Finally, the conditional ray distribution, $p(x, k|x^0, 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

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

$$ \text{with } \mathcal{N}(\bullet|\tilde{\sigma}_y^2(t)), \text{ a Gaussian function with variance } \tilde{\sigma}_y^2(t) = \sigma_y^2(t) + \left(\frac{v_g^0}{\bar{\omega}_r} \sin(\bar{\omega}_r t) \sigma_{\delta\theta_k}^0\right)^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 an 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 KE 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 shall then be used to more quantitatively interpreted the upper the turbulent 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 fast wave literature, isotropic diffusion and, hence, wavenumber diffusion may (e.g. Voronovich 1991) or may not (e.g. Bôas & Young 2020) comes into play (see Appendix
## 483 Appendix B. Single ray dynamics

The Itō noise ($\sigma dB_t / d\eta_t$) is white in time and conditionally Gaussian. Its conditional single-point distribution is fully determined by its zero mean and its local covariance matrix (given by equations (A1), (A7), (A13) and (A14)). In particular, we can replace this noise by another zero-mean Gaussian vector with the same covariance without changing the single-ray dynamics – typically replacing $\sigma dB_t$ by $\sqrt{a_0} \left( dB_t^{(1)} / dB_t^{(2)} \right)$ and $dZ_t$ by $-\sqrt{\gamma_0}k \left( dB_t^{(3)} / \sqrt{3}dB_t^{(4)} \right)$. It yields the simplified ray equations (4.3)-(4.4).

Then note that from Itō lemma (Oksendal 1998) $d\tilde{k} = d(\cos\theta_k / \sin\theta_k) = \tilde{k}^\perp d\theta_k - \frac{1}{2}\tilde{k}d<\theta_k, \theta_k>_t$ where $<\cdot, \cdot>_t$ denotes the quadratic covariation. Thus, $d\tilde{k} = dk\tilde{k} + kd\tilde{k} + d<\cdot, \cdot>_t = (dk - \frac{1}{2}kd<\theta_k, \theta_k>_t)\tilde{k} + (kd\theta_k + d<\cdot, \cdot>_t)\tilde{k}^\perp$ (B 1). Projecting this equation and $d\tilde{k} = -\nabla\tilde{v}^\top\tilde{k}dt + d\eta_t$ on $\tilde{k}$ and $\tilde{k}^\perp$, we have

$$494 \qquad \left\{ \begin{aligned} dk &= -\tilde{k} \cdot \nabla \tilde{v}^\top \tilde{k} dt + (dZ_t)_1 + \frac{1}{2} k d<\theta_k, \theta_k>_t \\ k d\theta_k &= -\tilde{k}^\perp \cdot \nabla \tilde{v}^\top \tilde{k} dt + (dZ_t)_2 - d<\cdot, \cdot>_t \end{aligned} \right. , \qquad (B 2)$$

$$495 \qquad \left\{ \begin{aligned} dk &= -\tilde{k} \cdot \nabla \tilde{v}^\top \tilde{k} dt + (dZ_t)_1 + \frac{1}{2} k^{-1} d<Z_2, Z_2>_t \\ d\theta_k &= -\tilde{k}^\perp \cdot \nabla \tilde{v}^\top \tilde{k} dt + k^{-1} (dZ_t)_2 + \frac{1}{2} k^{-2} d<Z_1, Z_2>_t \end{aligned} \right. . \qquad (B 3)$$

The treatment of the large-scale terms $\tilde{k} \cdot \nabla \tilde{v}^\top \tilde{k}$ and $\tilde{k}^\perp \cdot \nabla \tilde{v}^\top \tilde{k}$ is classical. Interested readers can refer to Lapeyre et al. (1999) for details. From Ito lemma again, $d\log k = dk/k - \frac{1}{2}d<\cdot, \cdot>_t/k^2$ leading to the simplified wave-vector dynamics (4.5)-(4.6).

## 499 Appendix C. Subgrid flow anisotropy and comparison with other works

Throughout this paper, we have considered an isotropic model for the stochastic subgrid velocity (3.10). The isotropic diffusivity matrix $a = a_0 \mathbb{I}_d$ is a good illustration of this. In contrast, many authors (e.g. White & Fornberg 1998; Bôas & Young 2020; Smit & Janssen 2019) assume isotropic and homogeneous turbulence and obtain anisotropic stochastic subgrid models for $||v||_g \to 0$. In these approaches, the integral over $\delta\theta$ in diffusivity matrix computations (A4) and (A12) involve singular integrations over the direction $v_g^0 = v_g^0 \tilde{k}$. It makes appear a Dirac delta function, $2\pi\delta(\kappa \cdot v_g^0) = \frac{2\pi}{\kappa v_g^0} (\delta(\theta_\kappa - \theta_k - \frac{\pi}{2}) + \delta(\theta_\kappa - \theta_k + \frac{\pi}{2}))$ (see Appendix in Bôas & Young 2020). This precision imposes a statistical anisotropy for $\sigma dB_t$ (oriented along $k$) and $d\eta_t$ (oriented along $k^\perp$), eventually leading to a covariance $\Sigma_Z = \gamma_0 k^2 \begin{bmatrix} 0 & 0 \\ 0 & 16 \end{bmatrix}$ (Eq. (3.17) in Bôas & Young (2020) and Eq. (24) in Smit & Janssen (2019)), no noise $dZ_1$, and no Brownian motion $B_t^{(3)}$. Moreover, because of scaling assumption, Bôas & Young (2020) neglect the spatial diffusivity matrix, $a$, while Smit & Janssen (2019) find $a = 4a_0 (\mathbb{I}_d + \frac{5}{4}\tilde{k}\tilde{k}^\top)$ (Eq. (22)-(23)). In this anisotropic framework, the Stratonovich wavevector equation (2.3), $dk = -\nabla(\bar{v}dt + \sigma \circ dB_t)^\top k$, would involve an additional drift term in Itō notations.

Further developing this anisotropic stochastic closure is an interesting avenue. Neverthe-
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 analyze wave dynamics with subsequent fast simulations of ensembles. Beside comprehension and analysis, our stochastic simulation tools aim at eventually facilitate future ensemble-based data assimilation algorithms (Smit et al. 2021).

**Funding.** This work is supported by the R&T CNES R-S19/OT-0003-084, the ERC project 856408-STUOD, the European Space Agency World Ocean Current project (ESA Contract No. 4000130730/20/I-NB), and SCALIAN DS.

**Declaration of interests.** The authors report no conflict of interest.

**Data availability statement.** SCALIAN DS owns a portion of the developed code intellectual property. For commercial reasons, that code will remains private.

**Author ORCIDs.** V. Resseguier, https://orcid.org/0000-0002-9301-9493; E. Hascoët https://orcid.org/0009-0002-0760-6534; B. Chapron, https://orcid.org/0000-0001-6088-8775

**Author contributions.** VR developed the theory. VR and EH wrote the code and performed numerical experiments. VR and BC wrote the paper.

## 488 Appendix A. Stochastic forcing covariance

In this appendix, we will compute the conditional covariance of the stochastic forcing of our eikonal characteristic equations (4.1), that is:

$$2\mathbf{D} \triangleq \frac{1}{dt} \mathbb{E}_t \left\{ \left( \frac{\sigma dB_t}{d\eta_t} \right) \left( \frac{\sigma dB_t}{d\eta_t} \right)^{\top} \right\} = \begin{bmatrix} \mathbf{a} & \Sigma_{\eta, \sigma} \\ \Sigma_{\eta, \sigma}^{\top} & \Sigma_{\eta} \end{bmatrix} \quad (A 1)$$

where, $d\eta_t = -\nabla(\sigma dB_t)^{\top}k$, denotes the wave-vector stochastic forcing and, $\Sigma_{\eta}dt$, its covariance, and $\mathbb{E}_t\{\bullet\} \triangleq \mathbb{E}\{\bullet|x_r(t), k(t)\}$ stands for the conditional expectation evaluated with given characteristics $(x_r(t), k(t))$ at the current time $t$. Note that in the appendix we use Itō notations only.

The subgrid velocity, $v' = \sigma dB_t/dt$ is constructed in Fourier space with a divergence-free isotropic spatial filter $\nabla^{\perp}\tilde{\psi}_{\sigma}$ (see (3.11)).

$$\tilde{v}'(\kappa, t) = \int dx \, v'(x, t) e^{-i\kappa \cdot x} = \overline{\sigma dB_t/dt}(\kappa, t) = i\kappa^{\perp} \tilde{\psi}_{\sigma}(\kappa) \overline{dB_t}(\kappa), \quad (A 2)$$

where $\kappa^{\perp}$ is the vector directly orthogonal to $\kappa$. The computation of the variance tensor $\mathbf{a}$ is classical and straightforward from the definition of the inverse Fourier transform and the identity $\mathbb{E}\{\overline{dB_t}(\kappa_1)\overline{dB_t}^*(\kappa_2)\} = (2\pi)^2\delta(\kappa_1 - \kappa_2)dt$, where * denotes complex conjugate. We simply need to split the integral of the stochastic forcing spectrum over the current wavevector
516 less, in the present study, we adopt the isotropic model for $\sigma d B_t$, which is much more

517 convenient for multi-ray numerical simulations.

## 518 Appendix D. Jet simulation

519 Again, currents are simulated at a resolution 512 × 512 on a 1000-km-width squared domain

520 through the same code. A backward velocity $v_{Bk}$ forces a leftward jet structure.

$$521 \qquad \partial_t \omega + \mathbf{v} \cdot \nabla \omega = S_{\omega} \quad \text{with} \quad \mathbf{v} = \nabla^{\perp} \Delta^{-1}(\omega + \omega_{Bk}). \qquad (\text{D 1})$$

522 $S_\omega$ encompasses the linear drag and the hyperviscosity, and the background vorticity, $\omega_{Bk}$,

523 is a smooth step function with a wavy interface at $y = Y_{Bk}(x)$:

$$524 \qquad \omega_{Bk}(x, y) = \Omega_{Bk} \left( \frac{1}{2} - \operatorname{erf} \left( \frac{y - Y_{Bk}(x)}{L_y^\omega} \right) \right) \quad \text{with } Y_{Bk}(x) = L_y \left( \frac{1}{2} + \frac{1}{30} \cos \left( \frac{2\pi}{L_x} x \right) \right) \qquad (\text{D 2})$$

525 To better highlight the interplay between ray oscillations and scattering, we consider
526 very-collimated swells, with a spatial extension of $100\lambda = 25$ km.

527 Besides, the curvature of the simulated jet can force an additional faster oscillation around

528 the jet for small enough wavevector angle. Indeed, a wave group traveling exactly rightward

529 would cross an alternation of positive and negative vorticity regions with a period $L_x/(v_g^0 - \overline{U}_0) \approx 1$ day $< 2\pi/\overline{\omega}_r$. Here, we set an initial wavevector angle large enough to prevent the

530 $\overline{U}_0 \approx 1$ day $< 2\pi/\overline{\omega}_r$. Here, we set an initial wavevector angle large enough to prevent the

531 additional harmonics.

## REFERENCES

532 BAL, GUILLAUME & CHOU, TOM 2002 Capillary-gravity wave transport over spatially random drift. *Wave*

533 *Motion* **35** (2), 107–124.

534 BAUER, WERNER, CHANDRAMOULI, PRANAV, CHAPRON, BERTRAND, LI, LONG & MÉMIN, ETIENNE 2020

535 Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic

536 approach. *Journal of Physical Oceanography* **50** (4), 983–1003.

537 BÔAS, ANA B VILLAS & YOUNG, WILLIAM R 2020 Directional diffusion of surface gravity wave action by

538 ocean macroturbulence. *Journal of Fluid Mechanics* **890**.

539 BORCEA, LILIANA, GARNIER, JOSSELIN & SOLNA, KNUT 2019 Wave propagation and imaging in moving

540 random media. *Multiscale Modeling & Simulation* **17** (1), 31–67.

541 BOURY, SAMUEL, BÜHLER, OLIVER & SHATAH, JALAL 2023 Fast-slow wave transitions induced by a random

542 mean flow. *Physical Review E* **108** (5), 055101.

543 BÜHLER, OLIVER 2009 *Waves and mean flows*. Cambridge University Press.

544 COTTER, COLIN J, GOTTWALD, GEORG A & HOLM, DARRYL D 2017 Stochastic partial differential fluid

545 equations as a diffusive limit of deterministic lagrangian multi-time dynamics. *Proceedings of the*

546 *Royal Society A: Mathematical, Physical and Engineering Sciences* **473** (2205), 20170388.

547 COX, MICHAEL R, KAFIABAD, HOSSEIN A & VANNESTE, JACQUES 2023 Inertia-gravity-wave diffusion by

548 geostrophic turbulence: the impact of flow time dependence. *Journal of Fluid Mechanics* **958**, A21.

549 CRISAN, DAN & HOLM, DARRYL D 2018 Wave breaking for the stochastic camassa–holm equation. *Physica*

550 *D: Nonlinear Phenomena* **376**, 138–143.

551 DINVAY, EVGUENI & MÉMIN, ETIENNE 2022 Hamiltonian formulation of the stochastic surface wave problem.

552 *Proceedings of the Royal Society A* **478** (2265), 20220050.

553 DONG, WENJING, BÜHLER, OLIVER & SMITH, K SHAFER 2020 Frequency diffusion of waves by unsteady

554 flows. *Journal of Fluid Mechanics* **905**, R3.

555 DYSTHE, KRISTIAN B 2001 Refraction of gravity waves by weak current gradients. *Journal of Fluid Mechanics*

556 **442**, 157–159.

557 FREUND, JB & FLEISCHMAN, TG 2002 Ray traces through unsteady jet turbulence. *International Journal of*

558 *Aeroacoustics* **1** (1), 83–96.

559 GARNIER, JOSSELIN, GAY, ETIENNE & SAVIN, ERIC 2020 Multiscale analysis of spectral broadening of acoustic

560 waves by a turbulent shear layer. *Multiscale Modeling & Simulation* **18** (2), 798–823.
$503 \quad \kappa = \kappa(\cos \theta_{\kappa}, \sin \theta_{\kappa}):$

$504 \quad \mathbf{a} = \frac{1}{(2\pi)^4 dt} \iint d\boldsymbol{\kappa}_1 d\boldsymbol{\kappa}_2 \mathbb{E}_t \left\{ (\overline{\sigma dB_t})(\boldsymbol{\kappa}_1, \mathbf{k}) (\overline{\sigma dB_t}^T)^*(\boldsymbol{\kappa}_2, \mathbf{k}) \right\} e^{i(\boldsymbol{\kappa}_1 - \boldsymbol{\kappa}_2) \cdot \mathbf{x}}, \quad (A3)

$505 \quad = \frac{1}{(2\pi)^2} \int d\boldsymbol{\kappa} \, \kappa^2 |\tilde{\psi}_{\sigma}(\boldsymbol{\kappa})|^2 \begin{pmatrix} -\sin \theta_{\boldsymbol{\kappa}} \\ \cos \theta_{\boldsymbol{\kappa}} \end{pmatrix} \begin{pmatrix} -\sin \theta_{\boldsymbol{\kappa}} \\ \cos \theta_{\boldsymbol{\kappa}} \end{pmatrix}^T, \quad (A4)$

$506 \quad = \frac{1}{(2\pi)^2} \int_0^{+\infty} \oint_0^{2\pi} d\boldsymbol{\kappa} d\theta_{\boldsymbol{\kappa}} \kappa^3 |\tilde{\psi}_{\sigma}(\boldsymbol{\kappa})|^2 \begin{pmatrix} \sin^2 \theta_{\boldsymbol{\kappa}} & -\sin \theta_{\boldsymbol{\kappa}} \cos \theta_{\boldsymbol{\kappa}} \\ -\sin \theta_{\boldsymbol{\kappa}} \cos \theta_{\boldsymbol{\kappa}} & \cos^2 \theta_{\boldsymbol{\kappa}} \end{pmatrix}, \quad (A5)$

$507 \quad = \frac{2}{2\pi} a_0 \oint_0^{2\pi} d\theta_{\boldsymbol{\kappa}} \begin{pmatrix} \sin^2 \theta_{\boldsymbol{\kappa}} & -\sin \theta_{\boldsymbol{\kappa}} \cos \theta_{\boldsymbol{\kappa}} \\ -\sin \theta_{\boldsymbol{\kappa}} \cos \theta_{\boldsymbol{\kappa}} & \cos^2 \theta_{\boldsymbol{\kappa}} \end{pmatrix}, \quad (A6)$

$508 \quad = a_0 \mathbb{I}_d, \quad (A7)$

where $a_0$ is defined by (4.7).

Now, the Fourier transform of the wave-vector stochastic forcing is

$511 \quad d\hat{\eta}_t = -\nabla(\sigma dB_t)^T \mathbf{k} = -i\boldsymbol{\kappa}(i\boldsymbol{\kappa}^{\perp}\tilde{\psi}_{\sigma}\overline{dB_t}) \cdot \mathbf{k} = \boldsymbol{\kappa}(\boldsymbol{\kappa}^{\perp} \cdot \mathbf{k})\tilde{\psi}_{\sigma}\overline{dB_t} = -\boldsymbol{\kappa}(\mathbf{k}^{\perp} \cdot \boldsymbol{\kappa})\tilde{\psi}_{\sigma}\overline{dB_t}. \quad (A8)$

Then, applying the crest-oriented rotation matrix, $M_k = [\tilde{k} \ \tilde{k}^{\perp}]$, leads to

$513 \quad d\hat{Z}_t = M_k^T d\hat{\eta}_t = - \begin{pmatrix} \tilde{k} \cdot \boldsymbol{\kappa} \\ \tilde{k}^{\perp} \cdot \boldsymbol{\kappa} \end{pmatrix} (\mathbf{k}^{\perp} \cdot \boldsymbol{\kappa}) \tilde{\psi}_{\sigma} \overline{dB_t} = - \begin{pmatrix} \cos \delta\theta \sin \delta\theta \\ \sin^2 \delta\theta \end{pmatrix} \kappa^2 \mathbf{k} \tilde{\psi}_{\sigma} \overline{dB_t}, \quad (A9)$

with $\delta\theta = \theta_{\kappa} - \theta_k$. From there, we can evaluate the conditional covariance matrix $\Sigma_Z = \frac{1}{dt} \mathbb{E}_t \{dZ_t dZ_t^T\}$ of $dZ_t$ as before:

$516 \quad \Sigma_Z = \frac{1}{(2\pi)^4 dt} \iint d\boldsymbol{\kappa}_1 d\boldsymbol{\kappa}_2 \mathbb{E}_t \left\{ (d\hat{Z}_t)(\boldsymbol{\kappa}_1, \mathbf{k})(d\hat{Z}_t^T)^*(\boldsymbol{\kappa}_2, \mathbf{k}) \right\} e^{i(\boldsymbol{\kappa}_1 - \boldsymbol{\kappa}_2) \cdot \mathbf{x}}, \quad (A10)$

$517 \quad = \frac{1}{(2\pi)^2} \int_0^{+\infty} \oint_0^{2\pi} d\boldsymbol{\kappa} d\delta\theta \, \kappa^5 k^2 |\tilde{\psi}_{\sigma}(\boldsymbol{\kappa})|^2 \begin{pmatrix} \cos^2 \delta\theta \sin^2 \delta\theta & \cos \delta\theta \sin^3 \delta\theta \\ \cos \delta\theta \sin^3 \delta\theta & \sin^4 \delta\theta \end{pmatrix}, \quad (A11)$

$518 \quad = \gamma_0 k^2 \begin{bmatrix} 1 & 0 \\ 0 & 3 \end{bmatrix}. \quad (A12)$

Finally, we come back to the canonical frame to get $\Sigma_{\eta}$

$520 \quad \Sigma_{\eta} = \mathbb{E}_t \{d\eta_t d\eta_t^T\} = M_k \Sigma_Z M_k^T = \gamma_0 k^2 \left[ \tilde{k} \tilde{k}^T + 3 \tilde{k}^{\perp} (\tilde{k}^{\perp})^T \right]. \quad (A13)$

For noises cross-correlations, by isotropy, it is also straightforward to show that

$522 \quad \Sigma_{\eta, \sigma} = 0. \quad (A14)$

The stochastic forcings of $x_r$ and $k$ are hence (conditionally) independent from one another.

## Appendix B. Single ray dynamics

The Itō noise $(\sigma dB_t) / d\eta_t$ is white in time and conditionally Gaussian. Its conditional single-point distribution is fully determined by its zero mean and its local covariance matrix (given by equations (A1), (A7), (A13) and (A14)). In particular, we can replace this noise by another zero-mean Gaussian vector with the same covariance without changing the single-
561 HELL, MOMME C, FOX-KEMPER, BAYLOR & CHAPRON, BERTRAND in preparation An efficient wave model for surface wave growth and propagation in coupled climate models. *Journal of Advances in Modeling Earth Systems*.

562
563
564 HELLER, EJ, KAPLAN, L & DAHLEN, A 2008 Refraction of a gaussian seaway. *Journal of Geophysical Research: Oceans* **113** (C9).

566 HOLM, D. 2015 Variational principles for stochastic fluid dynamics. *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences* **471** (2176).

568
569
570 HOLM, DARRYL D 2021 Stochastic variational formulations of fluid wave-current interaction. *Journal of nonlinear science* **31** (1), 4.

571 HOLM, DARRYL D, HU, RUIAO & STREET, OLIVER D 2023 On the interactions between mean flows and inertial gravity waves. *arXiv preprint arXiv:2302.04838*.

572
573 HOLM, DARRYL D & LUESINK, ERWIN 2021 Stochastic wave-current interaction in thermal shallow water dynamics. *Journal of Nonlinear Science* **31**, 1–56.

574 JONSSON, IVAR G 1990 Wave-current interactions. *The sea* **9**, 65–120.

575 KAFIABAD, HOSSEIN A, SAVVA, MILES AC & VANNESTE, JACQUES 2019 Diffusion of inertia-gravity waves by geostrophic turbulence. *Journal of Fluid Mechanics* **869**, R7.

577 KLYATSKIN, V. 2005 *Stochastic equations through the eye of the physicist: Basic concepts, exact results and asymptotic approximations*. Elsevier.

579 KLYATSKIN, VI & KOSHEL, KV 2015 Anomalous sea surface structures as an object of statistical topography. *Physical Review E* **91** (6), 063003.

581 KUDRYAVTSEV, VLADIMIR, YUROVSKAYA, MARIA, CHAPRON, BERTRAND, COLLARD, FABRICE & DONLON, CRAIG 2017 Sun glitter imagery of ocean surface waves. part 1: Directional spectrum retrieval and validation. *Journal of Geophysical Research: Oceans* **122** (2), 1369–1383.

584 KUNITA, H. 1997 *Stochastic flows and stochastic differential equations*,, vol. 24. Cambridge university press.

585 LAPEYRE, GUILLAUME 2017 Surface quasi-geostrophy. *Fluids* **2** (1), 7.

586 LAPEYRE, G., KLEIN, P. & HUA, B. 1999 Does the tracer gradient vector align with the strain eigenvectors in 2D turbulence? *Physics of Fluids* **11** (12), 3729–3737.

588 LAVRENOV, IGOR 2013 *Wind-waves in oceans: dynamics and numerical simulations*. Springer Science & Business Media.

590 MCCOMAS, C HENRY & BRETHERTON, FRANCIS P 1977 Resonant interaction of oceanic internal waves. *Journal of Geophysical Research* **82** (9), 1397–1412.

592 MÉMIN, E. 2014 Fluid flow dynamics under location uncertainty. *Geophysical & Astrophysical Fluid Dynamics* **108** (2), 119–146.

594 MÉMIN, ETIENNE, LI, LONG, LAHAYE, NOÉ, TISSOT, GILLES & CHAPRON, BERTRAND 2022 Linear wave solutions of a stochastic shallow water model. In *Stochastic Transport in Upper Ocean Dynamics Annual Workshop*, pp. 223–245. Springer Nature Switzerland Cham.

597 MIKULEVICIUS, R. & ROZOVSKII, B. 2004 Stochastic Navier–Stokes equations for turbulent flows. *SIAM Journal on Mathematical Analysis* **35** (5), 1250–1310.

599 OKSENDAL, B. 1998 *Stochastic differential equations*. Spinger-Verlag.

600 PAPANICOLAU, G. & KOHLER, W. 1974 Asymptotic theory of mixing stochastic ordinary differential equations. *Communications on Pure and Applied Mathematics* **27** (5), 641–668.

602 PIERREHUMBERT, RT 1994 Tracer microstructure in the large-eddy dominated regime. *Chaos, Solitons & Fractals* **4** (6), 1091–1110.

604 PITERBARG, LEONID & OSTROVSKII, A 1997 *Advection and diffusion in random media: implications for sea surface temperature anomalies*. Kluwer Academic.

606 PLOUGONVEN, RIWAL & ZHANG, FUQING 2014 Internal gravity waves from atmospheric jets and fronts. *Reviews of Geophysics* **52** (1), 33–76.

608 RESSEGUIER, VALENTIN, LI, LONG, JOUAN, GABRIEL, DÉRIAN, PIERRE, MÉMIN, ETIENNE & BERTRAND, CHAPRON 2020a New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics. *Archives of Computational Methods in Engineering* pp. 1–82.

611 RESSEGUIER, VALENTIN, MÉMIN, ETIENNE & CHAPRON, BERTRAND 2017a Geophysical flows under location uncertainty, part I random transport and general models. *Geophysical & Astrophysical Fluid Dynamics* **111** (3), 149–176.

614 RESSEGUIER, VALENTIN, MÉMIN, ETIENNE & CHAPRON, BERTRAND 2017b Geophysical flows under location uncertainty, part II quasi-geostrophy and efficient ensemble spreading. *Geophysical & Astrophysical Fluid Dynamics* **111** (3), 177–208.

617 RESSEGUIER, VALENTIN, PAN, WEI & FOX-KEMPER, BAYLOR 2020b Data-driven versus self-similar
529 ray dynamics – typically replacing $\sigma dB_t$ by $\sqrt{a_0} \left( \frac{dB_t^{(1)}}{dB_t^{(2)}} \right)$ and $dZ_t$ by $-\sqrt{\gamma_0}k \left( \frac{dB_t^{(3)}}{\sqrt{3}dB_t^{(4)}} \right)$. It
530 yields the simplified ray equations (4.3)-(4.4).

531 Then note that from Itō lemma (Oksendal 1998) $d\tilde{k} = d(\cos\theta_k) = \tilde{k}^{\perp}d\theta_k - \frac{1}{2}\tilde{k}dt <$
532 $\theta_k, \theta_k)_t$ where $< \bullet, \bullet>_t$ denotes the quadratic covariation. Thus,
533 $dk = dk\tilde{k} + kd\tilde{k} + d\langle k, \tilde{k} \rangle = (dk - \frac{1}{2}kd\langle \theta_k, \theta_k \rangle_t)\tilde{k} + (kd\theta_k + d\langle k, \theta_k \rangle_t)\tilde{k}^{\perp}$ (B 1)
534 Projecting this equation and $dk = -\nabla\tilde{v}^T kdt + d\eta_t$ on $\tilde{k}$ and $\tilde{k}^{\perp}$, we have

$$
535 \begin{cases} dk = -\tilde{k} \cdot \nabla \tilde{v}^T k dt + (dZ_t)_1 + \frac{1}{2} k d \langle \theta_k, \theta_k \rangle_t \\ k d \theta_k = -\tilde{k}^{\perp} \cdot \nabla \tilde{v}^T k dt + (dZ_t)_2 - d \langle k, \theta_k \rangle_t \end{cases} , \qquad (B 2)
$$

$$
536 \begin{cases} dk = -\tilde{k} \cdot \nabla \tilde{v}^T k dt + (dZ_t)_1 + \frac{1}{2} k^{-1} d \langle Z_2, Z_2 \rangle_t \\ d\theta_k = -\tilde{k}^{\perp} \cdot \nabla \tilde{v}^T \tilde{k} dt + k^{-1} (dZ_t)_2 + \frac{1}{2} k^{-2} d \langle Z_1, Z_2 \rangle_t \end{cases} . \quad (B 3)
$$

537 The treatment of the large-scale terms $\tilde{k} \cdot \nabla \tilde{v}^T \tilde{k}$ and $\tilde{k}^{\perp} \cdot \nabla \tilde{v}^T \tilde{k}$ is classical. Interested
538 readers can refer to Lapeyre et al. (1999) for details. From Ito lemma again, $d \log k =
539 dk/k - \frac{1}{2} d \langle k, k \rangle_t / k^2$ leading to the simplified wave-vector dynamics (4.5)-(4.6).

540 Appendix C. Subgrid flow anisotropy and comparison with other works

541 Throughout this paper, we have considered an isotropic model for the stochastic subgrid
542 velocity (3.11). The isotropic diffusivity matrix **a** = $a_0 \mathbb{I}_d$ is a good illustration of this. In
543 contrast, many authors (e.g. White & Fornberg 1998; Bôas & Young 2020; Smit & Janssen
544 2019) assume isotropic and homogeneous turbulence and obtain anisotropic stochastic
545 subgrid models for $\frac{\|\mathbf{v}\|}{v_g^0} \to 0$. In these approaches, the integral over $\delta\theta$ in diffusivity matrix
546 computations (A 4) and (A 12) involve singular integrations over the direction $\mathbf{v}_g^0 = v_g^0 \tilde{\mathbf{k}}$. It
547 makes appear a Dirac delta function, $2\pi\delta(\mathbf{\kappa} \cdot \mathbf{v}_g^0) = \frac{2\pi}{\kappa v_g^0} (\delta(\theta_\kappa - \theta_k - \frac{\pi}{2}) + \delta(\theta_\kappa - \theta_k + \frac{\pi}{2}))$
548 (see Appendix in Bôas & Young 2020). This precision imposes a statistical anisotropy for
549 $\sigma dB_t$ (oriented along **k**) and $d\eta_t$ (oriented along $\mathbf{k}^{\perp}$), eventually leading to a covariance
550 $\Sigma_Z = \gamma_0 k^2 \begin{bmatrix} 0 & 0 \\ 0 & 16 \end{bmatrix}$ (Eq. (3.17) in Bôas & Young (2020) and Eq. (24) in Smit & Janssen
551 (2019)), no noise $dZ_1$, and no Brownian motion $B_t^{(3)}$. Moreover, because of scaling
552 assumption, Bôas & Young (2020) neglect the spatial diffusivity matrix, **a**, while Smit
553 & Janssen (2019) find **a** = $4a_0 (\mathbb{I}_d + \frac{5}{4} \tilde{\mathbf{k}} \tilde{\mathbf{k}}^T)$ (Eq. (22)-(23)). In this anisotropic framework,
554 the Stratonovich wavevector equation (2.3), $dk = -\nabla(\bar{v}dt + \sigma \circ dB_t)^T k$, would involve an
555 additional drift term in Itô notations.

556 Further developing this anisotropic stochastic closure is an interesting avenue. A multiscale
557 anisotropic stochastic closure would involve wavenumber variations but no wavenumber
558 diffusion. Nevertheless, in the present study, we adopt the isotropic model for $\sigma dB_t$, which
559 is much more convenient for multi-ray numerical simulations.

560 Appendix D. Action spectra and ray distribution

561 Here we highlight the link between mean action spectral density and the ray distribution.
562 We denote by $N^0$ the initial wave action spectrum. We first use the definition of the Dirac
618 parameterizations for stochastic advection by lie transport and location uncertainty. *Nonlinear*

619 *Processes in Geophysics* **27** (2), 209–234.

620 SLUNYAEV, AV & SHRIRA, VI 2023 Extreme dynamics of wave groups on jet currents. *Physics of Fluids*

621 **35** (12).

622 SMIT, PB, HOUGHTON, IA, JORDANOVA, K, PORTWOOD, T, SHAPIRO, E, CLARK, D, SOSA, M & JANSSEN, TT

623 2021 Assimilation of significant wave height from distributed ocean wave sensors. *Ocean Modelling*

624 **159**, 101738.

625 SMIT, PIETER B & JANSSEN, TIM T 2019 Swell propagation through submesoscale turbulence. *Journal of*

626 *Physical Oceanography* **49** (10), 2615–2630.

627 VORONOVICH, AG 1976 Propagation of internal and surface gravity waves in the approximation of geometrical

628 optics. *Izv. Atmos. Ocean. Phys* **12**, 850–857.

629 VORONOVICH, A. 1991 The effect of shortening of waves on random currents. In *Proceedings of nonlinear*

630 water waves. Bristol.

631 WANG, HAN, BÔAS, ANA B VILLAS, YOUNG, WILLIAM R & VANNESTE, JACQUES 2023 Scattering of swell by

632 currents. *arXiv preprint arXiv:2305.12163*.

633 WEST, BRUCE J 1978 Ray paths in a fluctuating environment. *Physical Review A* **18** (4), 1646.

634 WHITE, BENJAMIN S 1999 Wave action on currents with vorticity. *Journal of Fluid Mechanics* **386**, 329–344.

635 WHITE, BENJAMIN S & FORNBERG, BENGT 1998 On the chance of freak waves at sea. *Journal of fluid*

636 mechanics **355**, 113–138.

637 ZHEN, YICUN, RESSEGUIER, VALENTIN & CHAPRON, BERTRAND 2023 Physically constrained covariance

638 inflation from location uncertainty. *EGUsphere* **2023**, 1.
measure then we perform a variable change corresponding to the characteristic (2.5) from

$t = t_i$ to $t = t_f$:

$$565 \qquad \mathbb{E}N(\mathbf{x}, \mathbf{k}, t) = \mathbb{E} \iint d\mathbf{x}_r d\mathbf{k}_r \ N(\mathbf{x}_r, \mathbf{k}_r, t)\delta(\mathbf{x}_r - \mathbf{x})\delta(\mathbf{k}_r - \mathbf{k}), \qquad (\text{D 1})$$

$$566 \qquad = \mathbb{E} \iint d\mathbf{x}_r^0 d\mathbf{k}_r^0 \ N(\mathbf{x}_r^0, \mathbf{k}_r^0, 0)\delta(\mathbf{x}_r(\mathbf{x}_r^0, \mathbf{k}_r^0, t) - \mathbf{x})\delta(\mathbf{k}_r(\mathbf{x}_r^0, \mathbf{k}_r^0, t) - \mathbf{k}), \quad (\text{D 2})$$

$$567 \qquad = \iint d\mathbf{x}_r^0 d\mathbf{k}_r^0 \ N^0(\mathbf{x}_r^0, \mathbf{k}_r^0) p(\mathbf{x}, \mathbf{k} | \mathbf{x}_r^0, \mathbf{k}_r^0, t), \qquad (\text{D 3})$$

where the standard relation between the Dirac measure and the probability distribution

function has been used.

## 570 Appendix E. Jet simulation

Again, currents are simulated at a resolution 512 × 512 on a 1000-km-width squared domain

$[0, L_x] \times [0, L_y]$ through the same code. A backward velocity $v_{Bk}$ forces a leftward jet

structure.

$$574 \qquad \partial_t \omega + \mathbf{v} \cdot \nabla \omega = S_{\omega} \quad \text{with} \quad \mathbf{v} = \nabla^{\perp} \Delta^{-1}(\omega + \omega_{Bk}). \qquad (\text{E 1})$$

$S_{\omega}$ encompasses the linear drag and the hyperviscosity with coefficient $1/\tau_F = 3.22 \times 10^{-8} s^{-1}$

and $v_{HV}/dx^8 = 3.33 \times 10^{-9} s^{-1}$ respectively. The background vorticity, $\omega_{Bk}$, is a smooth step

function with a wavy interface at $y = Y_{Bk}(x)$:

$$578 \qquad \omega_{Bk}(x, y) = \Omega_{Bk} \left( \frac{1}{2} - \text{erf} \left( \frac{y - Y_{Bk}(x)}{L_y^\omega} \right) \right) \quad \text{with} \quad Y_{Bk}(x) = L_y \left( \frac{1}{2} + \frac{1}{30} \cos \left( \frac{2\pi}{L_x} x \right) \right). \quad (\text{E 2})$$

To better highlight the interplay between ray oscillations and scattering, we consider very-collimated swells, with a spatial extension of $100\lambda = 25$ km.

Besides, the curvature of the simulated jet can force an additional faster oscillation around

the jet for small enough wavevector angle. Indeed, a wave group traveling exactly rightward

would cross an alternation of positive and negative vorticity regions with a period $L_x/(v_g^0 - \bar{U}_0) \approx 1$ day $< 2\pi/\bar{\omega}_r$. Here, we set an initial wavevector angle large enough to prevent the

additional harmonics.

## 585 REFERENCES

586 BAL, GUILLAUME & CHOU, TOM 2002 Capillary-gravity wave transport over spatially random drift. *Wave*

587 *Motion* **35** (2), 107–124.

588 BAUER, WERNER, CHANDRAMOULI, PRANAV, CHAPRON, BERTRAND, LI, LONG & MÉMIN, ETIENNE 2020

589 Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic

590 approach. *Journal of Physical Oceanography* **50** (4), 983–1003.

591 BÔAS, ANA B VILLAS & YOUNG, WILLIAM R 2020 Directional diffusion of surface gravity wave action by

592 ocean macroturbulence. *Journal of Fluid Mechanics* **890**.

593 BORCEA, LILIANA, GARNIER, JOSSELIN & SOLNA, KNUT 2019 Wave propagation and imaging in moving

594 random media. *Multiscale Modeling & Simulation* **17** (1), 31–67.

595 BOURY, SAMUEL, BÜHLER, OLIVER & SHATAH, JALAL 2023 Fast-slow wave transitions induced by a random

596 mean flow. *Physical Review E* **108** (5), 055101.

597 BÜHLER, OLIVER 2009 *Waves and mean flows*. Cambridge University Press.

598 COTTER, COLIN J, GOTTWALD, GEORG A & HOLM, DARRYL D 2017 Stochastic partial differential fluid

599 equations as a diffusive limit of deterministic lagrangian multi-time dynamics. *Proceedings of the*

600 Royal Society A: Mathematical, Physical and Engineering Sciences **473** (2205), 20170388.
601 COX, MICHAEL R, KAFIABAD, HOSSEIN A & VANNESTE, JACQUES 2023 Inertia-gravity-wave diffusion by geostrophic turbulence: the impact of flow time dependence. *Journal of Fluid Mechanics* **958**, A21.

602 CRISAN, DAN & HOLM, DARRYL D 2018 Wave breaking for the stochastic camassa–holm equation. *Physica D: Nonlinear Phenomena* **376**, 138–143.

605 DINVAY, EVGUENI & MÉMIN, ETIENNE 2022 Hamiltonian formulation of the stochastic surface wave problem. *Proceedings of the Royal Society A* **478** (2265), 20220050.

607 DONG, WENJING, BÜHLER, OLIVER & SMITH, K SHAFER 2020 Frequency diffusion of waves by unsteady flows. *Journal of Fluid Mechanics* **905**, R3.

609 DYSTHE, KRISTIAN B 2001 Refraction of gravity waves by weak current gradients. *Journal of Fluid Mechanics* **442**, 157–159.

611 FREUND, JB & FLEISCHMAN, TG 2002 Ray traces through unsteady jet turbulence. *International Journal of Aeroacoustics* **1** (1), 83–96.

613 GARNIER, JOSSELIN, GAY, ETIENNE & SAVIN, ERIC 2020 Multiscale analysis of spectral broadening of acoustic waves by a turbulent shear layer. *Multiscale Modeling & Simulation* **18** (2), 798–823.

615 HELD, I., PIERREHUMBERT, R., GARNER, S. & SWANSON, K. 1995 Surface quasi-geostrophic dynamics. *Journal of Fluid Mechanics* **282**, 1–20.

617 HELL, MOMME C, FOX-KEMPER, BAYLOR & CHAPRON, BERTRAND in preparation An efficient wave model for surface wave growth and propagation in coupled climate models. *Journal of Advances in Modeling Earth Systems*.

620 HELLER, EJ, KAPLAN, L & DAHLEN, A 2008 Refraction of a gaussian seaway. *Journal of Geophysical Research: Oceans* **113** (C9).

622 HOLM, D. 2015 Variational principles for stochastic fluid dynamics. *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences* **471** (2176).

624 HOLM, DARRYL D 2021 Stochastic variational formulations of fluid wave–current interaction. *Journal of nonlinear science* **31** (1), 4.

626 HOLM, DARRYL D, HU, RUIAO & STREET, OLIVER D 2023 On the interactions between mean flows and inertial gravity waves. *arXiv preprint arXiv:2302.04838*.

628 HOLM, DARRYL D & LUESINK, ERWIN 2021 Stochastic wave–current interaction in thermal shallow water dynamics. *Journal of Nonlinear Science* **31**, 1–56.

630 JONSSON, IVAR G 1990 Wave-current interactions. *The sea* **9**, 65–120.

631 KAFIABAD, HOSSEIN A, SAVVA, MILES AC & VANNESTE, JACQUES 2019 Diffusion of inertia-gravity waves by geostrophic turbulence. *Journal of Fluid Mechanics* **869**, R7.

633 KLYATSKIN, V. 2005 *Stochastic equations through the eye of the physicist: Basic concepts, exact results and asymptotic approximations*. Elsevier.

635 KLYATSKIN, VI & KOSHEL, KV 2015 Anomalous sea surface structures as an object of statistical topography. *Physical Review E* **91** (6), 063003.

637 KUDRYAVTSEV, VLADIMIR, YUROVSKAYA, MARIA, CHAPRON, BERTRAND, COLLARD, FABRICE & DONLON, CRAIG 2017 Sun glitter imagery of ocean surface waves. part 1: Directional spectrum retrieval and validation. *Journal of Geophysical Research: Oceans* **122** (2), 1369–1383.

640 KUNITA, H. 1997 *Stochastic flows and stochastic differential equations*, vol. 24. Cambridge university press.

641 LANDAU, LEV DAVIDOVICH & LIFSHITS, EVGENIİ MIKHAİLOVICH 1960 *Mechanics*, vol. 1. CUP Archive.

642 LAPEYRE, GUILLAUME 2017 Surface quasi-geostrophy. *Fluids* **2** (1), 7.

643 LAPEYRE, G., KLEIN, P. & HUA, B. 1999 Does the tracer gradient vector align with the strain eigenvectors in 2D turbulence? *Physics of Fluids* **11** (12), 3729–3737.

645 LAVRENOV, IGOR 2013 *Wind-waves in oceans: dynamics and numerical simulations*. Springer Science & Business Media.

647 MCCOMAS, C HENRY & BRETHERTON, FRANCIS P 1977 Resonant interaction of oceanic internal waves. *Journal of Geophysical Research* **82** (9), 1397–1412.

649 MÉMIN, E. 2014 Fluid flow dynamics under location uncertainty. *Geophysical & Astrophysical Fluid Dynamics* **108** (2), 119–146.

651 MÉMIN, ETIENNE, LI, LONG, LAHAYE, NOÉ, TISSOT, GILLES & CHAPRON, BERTRAND 2022 Linear wave solutions of a stochastic shallow water model. In *Stochastic Transport in Upper Ocean Dynamics Annual Workshop*, pp. 223–245. Springer Nature Switzerland Cham.

654 MIKULEVICIUS, R. & ROZOVSKII, B. 2004 Stochastic Navier–Stokes equations for turbulent flows. *SIAM Journal on Mathematical Analysis* **35** (5), 1250–1310.

656 OKSENDAL, B. 1998 *Stochastic differential equations*. Spinger-Verlag.
657 PAPANICOLAOU, G. & KOHLER, W. 1974 Asymptotic theory of mixing stochastic ordinary differential equations. *Communications on Pure and Applied Mathematics* **27** (5), 641–668.

658 PITERBARG, LEONID & OSTROVSKII, A 1997 *Advection and diffusion in random media: implications for sea surface temperature anomalies*. Kluwer Academic.

661 PLOUGONVEN, RIWAL & ZHANG, FUQING 2014 Internal gravity waves from atmospheric jets and fronts. *Reviews of Geophysics* **52** (1), 33–76.

663 RESSEGUIER, VALENTIN, LI, LONG, JOUAN, GABRIEL, DÉRIAN, PIERRE, MÉMIN, ETIENNE & BERTRAND, CHAPRON 2020a New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics. *Archives of Computational Methods in Engineering* pp. 1–82.

666 RESSEGUIER, VALENTIN, MÉMIN, ETIENNE & CHAPRON, BERTRAND 2017a Geophysical flows under location uncertainty, part I random transport and general models. *Geophysical & Astrophysical Fluid Dynamics* **111** (3), 149–176.

669 RESSEGUIER, VALENTIN, MÉMIN, ETIENNE & CHAPRON, BERTRAND 2017b Geophysical flows under location uncertainty, part II quasi-geostrophy and efficient ensemble spreading. *Geophysical & Astrophysical Fluid Dynamics* **111** (3), 177–208.

672 RESSEGUIER, VALENTIN, PAN, WEI & FOX-KEMPER, BAYLOR 2020b Data-driven versus self-similar parameterizations for stochastic advection by lie transport and location uncertainty. *Nonlinear Processes in Geophysics* **27** (2), 209–234.

675 SLUNYAEV, AV & SHRIRA, VI 2023 Extreme dynamics of wave groups on jet currents. *Physics of Fluids* **35** (12).

677 SMIT, PB, HOUGHTON, IA, JORDANOVA, K, PORTWOOD, T, SHAPIRO, E, CLARK, D, SOSA, M & JANSSEN, TT 2021 Assimilation of significant wave height from distributed ocean wave sensors. *Ocean Modelling* **159**, 1017–1038.

680 SMIT, PIETER B & JANSSEN, TIM T 2019 Swell propagation through submesoscale turbulence. *Journal of Physical Oceanography* **49** (10), 2615–2630.

682 VORONOVICH, A. 1991 The effect of shortening of waves on random currents. In *Proceedings of nonlinear water waves*. Bristol.

684 WANG, HAN, BÔAS, ANA B VILLAS, YOUNG, WILLIAM R & VANNESTE, JACQUES 2023 Scattering of swell by currents. *arXiv preprint arXiv:2305.12163*.

686 WEST, BRUCE J 1978 Ray paths in a fluctuating environment. *Physical Review A* **18** (4), 1646.

687 WHITE, BENJAMIN S 1999 Wave action on currents with vorticity. *Journal of Fluid Mechanics* **386**, 329–344.

688 WHITE, BENJAMIN S & FORNBERG, BENGT 1998 On the chance of freak waves at sea. *Journal of fluid mechanics* **355**, 113–138.

690 ZHEN, YICUN, RESSEGUIER, VALENTIN & CHAPRON, BERTRAND 2023 Physically constrained covariance inflation from location uncertainty. *EGUsphere* **2023**, 1.