## Abstract

The propagation of ocean swells from generating regions to remote coastlines is affected by submesoscale turbulence in the surface flow field. The presence of submesoscale velocity variations results in random scattering of wave rays. While the interactions with these flow fields are weak, cumulative effects over oceanic scales are significant and result in observable changes in the wave field. Using geometrical optics and statistical mechanics we derive a framework to express these scattering effects on the mean wave statistics directly in terms of the variance spectrum of the submesoscale current field. The theoretical results are presented in Lagrangian and Eulerian forms, where the latter takes the form of a radiative transport equation augmented with a diffusive term in directional space. The theoretical results are verified through Monte Carlo simulations with a geometrical optics model. We show that including submesoscale scattering on ocean wave evolution can explain observed delays in swell arrivals, accelerated wave height decay, and much larger directional spreading of the wave field than predicted by geometrical spreading alone.

## 1. Introduction

Pioneering work in the mid-twentieth century (Munk and Traylor 1947; Barber and Ursell 1948; Munk et al. 1963; Snodgrass et al. 1966) explored the foundations of swell dynamics over large propagation distances in the open ocean, and explained observed narrowing (in frequency and directional space) of the wave spectrum by frequency dispersion and geometric spreading. The evolution of ocean waves over long distances is important for many physical processes, including wave-driven flows (e.g., Stokes drift, Langmuir circulation; Xu and Bowen 1994; McWilliams et al. 2004), air–sea interaction (transfer of momentum and upper ocean mixing; Sullivan and McWilliams 2010; Cavaleri et al. 2012), and cross-shelf mass transport (e.g., Lentz et al. 2008) and affects nearly every aspect of wave-driven coastal dynamics (e.g., longshore currents, rip currents, sediment transport; Longuet-Higgins 1970; Battjes 1974; MacMahan et al. 2006).

The cumulative effects of geometric spreading and dispersion successfully explain some characteristics of oceanic swell evolution. However, considerable differences between geometrical optics and observed evolution remain unaccounted for and are poorly understood. Specifically, directional spreading values (a measure of spectral width in directional space) are considerably larger than predicted by geometrical optics (e.g., Ewans 2002), travel times of swell fields are often incorrectly predicted (with bias toward late arrivals; Jiang et al. 2016a), and observed decay of swell fields remains unexplained (e.g., Ardhuin et al. 2009; Young et al. 2013; Stopa et al. 2016; Jiang et al. 2016b).

As an example we can consider wave arrivals along the California coast originating from a storm in the Southern Ocean. In this example, the storm is schematized as a directionally isotropic source with a 1500-km width traveling eastward at 30 km h^{−1} along a line of constant latitude (Fig. 1A). Assuming that the waves generated in the storm propagate along great circles we find that the width of the wave field approaching the California coast ranges from 1° to 10° (see Fig. 1). However, the historical record (spanning 10 years of data) at Coastal Data Information Program (CDIP) site 191 (moored in deep water at 1143-m depth near San Diego, CA) shows directional spread values ranging mostly between 15° and 25° and events narrower than 5° spreading are never observed. The considerable discrepancy between these observations and geometric-dispersive theory are not due to instrumental bias, as the instrument-type used (Datawell WaveRider) has been shown to yield unbiased estimates for directional spreading (O’Reilly et al. 1996). Further, similar directional widths are observed at other CDIP sites (not shown), and have been reported for southwestern swell at the west coast of New Zealand (Ewans 2002). Although multimodality of the wave field (multiple sources) can explain some of the bias toward larger spreads, the virtual absence of sub-10° events suggests that something is missing.

As ocean waves travel over long distances in open ocean, the refraction associated with mesoscale and submesoscale eddies will drive deviations from great-circle paths, extending pathlength, and delaying arrivals (Gallet and Young 2014; Ardhuin et al. 2017). Although mesoscale flow fields are generally more energetic than submesoscale (1–100 km) eddies^{1} (Callies et al. 2015), observations and numerical simulations indicate that submesoscale eddies are more energetic than previously presumed (Callies et al. 2013, 2015). Moreover, the shorter scales result in larger shear and vorticity values, so that the effective inhomogeneity presented to the wave field by submesoscale currents is similar to the mesoscale variability. During a transoceanic crossing, swell waves encounter hundreds of eddies of varying scale, which affects the wave characteristics in the far field.

The propagation of waves through random media has received considerable interest in the context of extreme wave statistics. Refraction creates constructive and destructive interference, that combined with nonlinear effects, may drive rapid (on the scale of individual waves) variability in the statistics and increase the probability of extreme wave events (e.g., Janssen and Herbers 2009; Fabrikant and Raevsky 1994; Metzger et al. 2014; Heller et al. 2008). These are typically near-field effects and thus important locally, but are averaged out in the far field (Smit and Janssen 2013; Heller et al. 2008). Refraction of waves over submesoscale flows will drive variations in phase-averaged properties, such as wave height, on the scale of the submesoscale flow (Ardhuin et al. 2017). On this scale, crossing waves originate from statistically independent regions, and will appear as a superposition of multidirectional uncorrelated waves, without coherent interference effects (Smit and Janssen 2013). These submesoscale effects however do affect the mean properties on the larger scales, such as directional spreading and pathlengths (and thus arrival times). Technically, if the submesoscale flow field is known, the resulting wave effects (submesoscale and larger) can be determined with operational third-generation wave models. Unfortunately ocean-wide submesoscale flow fields are not generally available, and resolving the global wave field at these scales is not feasible due to computational efforts needed. If we assume that the individual submesoscale eddies are statistically uncorrelated, and that submesoscale refraction does not affect energy conservation (e.g., no enhanced dissipation), describing mesoscale swell fields with mean properties is reasonable. However, neglecting the interaction with the small-scale flow altogether ignores the directional broadening of the wave field.

In the present work we investigate the large-scale effects of submesoscale turbulence on the directional properties of swell and how we might incorporate these effects in statistical (Eulerian) models without explicitly resolving the underlying flow fields. Our goal is to approximate the long-term mean effects of scattering due to submesoscale flows on the evolution of mesoscale properties of the wave field. We propose a theoretical framework based on geometrical optics to predict the average directional properties of a large ensemble of wave rays that propagate through a medium with turbulent velocity fluctuations. Specifically, we demonstrate that under the assumption of isotropic, homogeneous turbulence the long-range random scattering effects can be accounted for by including a directional diffusion term in the action balance equation. We will start our exploration using geometrical optics (section 2) to consider the average statistical properties of an ensemble of wave rays (section 3). Subsequently, we use this Lagrangian framework to derive additional diffusive contributions to the action balance (section 4). We validate predictions from stochastic theory using Monte Carlo simulations (section 5), discuss implications (section 6), and summarize our main results (section 7).

## 2. Deterministic framework

We consider propagation of surface waves in deep water, in the absence of dissipation or source inputs (conservative), disconnected from local winds (swell) and characterized by low steepness. The medium is represented by a surface current field where the ratio *ε* between a characteristic surface velocity and characteristic wave celerity is assumed small (*ε* ≪ 1). The influence of vertical shear on the wave motion is assumed sufficiently small so that we can represent the surface velocity field by a representative velocity field **u** = [*u*_{x}, *u*_{y}]^{T}, with *u*_{x}, *u*_{y} representing orthogonal horizontal velocities as a function of spatial coordinates *x*, *y* and time *t*. Further spatial changes in the flow field occur on temporal and spatial scales larger than a typical wave period and wavelength (see appendix A for scaling details). Since waves exchange energy with the flow wave energy density is not conserved. However, under the stated assumptions the wave action is (approximately) an adiabatic invariant (Bretherton and Garrett 1968; Andrews and Mcintyre 1978; Ardhuin et al. 2008), and our starting point will be conservation of wave action for ocean swells,

Here *E*_{k}(*t*) = *E*[**k**(*t*), **x**(*t*), *t*] is the energy density of a wave packet with wavenumber **k**(*t*) = [*k*_{x}(*t*), *k*_{y}(*t*)]^{T} at location **x**(*t*) = [*x*(*t*), *y*(*t*)]^{T} and time *t*. Further, *σ*(*t*) is the intrinsic frequency related to the wavenumber through the deep-water dispersion relation |**k**| = *σ*^{2}/*g*, with *g* the gravitational acceleration.

We express the wavenumber in terms of the wave direction *θ* and intrinsic frequency *σ* {using **k** = *k*(*σ*)[cos *θ*, sin *θ*]^{T}}, and denote the location of a wave packet in geographical and spectral space as **X** = [*x*, *y*, *θ*, *σ*]^{T}. In this case **X** evolves along the rays of geometrical optics according to

with

Here *u*_{n}[**x**(*t*), *θ*(*t*), *t*], *u*_{m}[**x**(*t*), *θ*(*t*), *t*] are velocities directed parallel and perpendicular to the local wave direction, and ∂_{n} and ∂_{m} are the derivatives along and perpendicular to the wave direction. Specifically,

with $A$ the rotation matrix for rotation angle *θ* defined as

Note that *θ*(*t*) is a Lagrangian (dependent) variable, whereas *θ* in Eqs. (4) and (5) is an independent variable, representing some angle of rotation. For convenience (and brevity) we consider wave motion in a Euclidean geometry since Earth’s curvature does not affect our results to the order considered. Equations (1)–(3) constitute a geometric optics approximation (e.g., Mei et al. 2005) to account for the effect of horizontal variations in depth uniform flows on the surface wave motion. We do not include effects due to, for example, velocity shear, higher-order wave–current or wave–wave interactions, and dissipation or generation. These may contribute to long-term swell evolution (e.g., swell dissipation; Ardhuin et al. 2009), but will be excluded here because our focus is on the leading-order scattering due to spatial variations in the surface flow.

Given the surface velocity field **u**, Eq. (2) can be integrated directly. However, global deterministic knowledge of the turbulent velocity field at submesoscales is unavailable, and the computational effort involved in resolving *O*(10) km scales for cross-basin propagation is formidable. We therefore split the near-surface horizontal flow field **u** in large-scale (mesoscale and larger scale) flows **U** and smaller-scale (submesoscale) fluctuations **u**′, written as

In operational models, typically only the large-scale flow field that varies on time and space scales *T*_{U} and *L*_{U} (mesoscale and larger) is resolved. To capture the effects of the smaller-scale dynamics on the wave propagation, we assume that **u′** can be represented as a zero-mean ergodic stochastic process that decorrelates on spatial and temporal scales *L*_{u′} and *T*_{u′} (with *L*_{u′} ≪ *L*_{U} and *T*_{u′} ≪ *T*_{U}). We further assume that the underlying random process is quasi-stationary and quasi-homogeneous, so that the variance of *u*′, changes on the *T*_{U} and *L*_{U} mesoscales. Further, we assume that waves travel fast relative to submesoscale time and space scales, expressed as

so that the submesoscale velocity field appears frozen from the waves perspective. This means that changes in the surface flow field as observed by an observer moving with the wave field are principally due to the propagation of the waves through the (frozen) flow field.

The nonlinear dependence of ray velocities on direction and intrinsic frequency complicates statistical analysis of the influence of **u**′ on the coupled ray equations. Because in the absence of the flow field the problem has rotational symmetry, we first consider *n* and *m* which measure the distances traveled in the crest normal and along-crest directions, respectively. Specifically,

with $A$(*t*) = $A$[*θ*(*t*)] the local rotation matrix.

We define the ray starting location as **x**(*t*_{0}) and *m*(*t*_{0}) = *n*(*t*_{0}) = 0. For a given ray starting at **X**_{0}, we can then express the local position **x**[**M**(*t*; **X**_{0}); **X**_{0}] as a function of **M** = [*m*, *n*, *θ*, *σ*]^{T} and the starting label **X**_{0}. The time evolution of **M**(*t*) then follows

with

and where *U*_{m}, *U*_{n}, and $um\u2032,un\u2032$ are the velocities directed along the local *m*, *n* directions as a function of position **x**, *θ*, and *t*. To leading order the dependence of **c** on directions disappears, but the nonlinear dependence of **c** on *σ* remains due to the group velocity *C*_{g} = *g*(2*σ*)^{−1}.

To further simplify we note that the absolute frequency *ω* is locally related to the intrinsic frequency and flow field through the dispersion relation

We decompose absolute and relative frequencies as $\sigma =\sigma \xaf+\sigma \u2032$ and $\omega =\omega \xaf+\omega \u2032$, where here (and in what follows) the bar indicates a mean contribution averaged over submesoscale (but varying in the mesoscale) and primed variables are fluctuating components about that mean. Over the time scales of *O*(10) days relevant for swell propagation we find that the fluctuating component to the absolute frequency may be neglected (see appendix A for details) and to *O*(*ε*) we find that $\omega \u2248\omega \xaf$ so that *ω* is directly tied to the large scale flow and mean intrinsic frequency as

In addition, absolute frequency fluctuations *σ*′ are directly tied to velocity fluctuations ( appendix A)

Since $\sigma \u2032/\sigma \xaf=O\u2061(\epsilon )$ we can now approximate *C*_{g} as a first order Taylor series, so that the along-ray velocity $cn=Cg+Un+un\u2032$ becomes

with $C\xafg=g\u2061(2\sigma \xaf)\u22121$. From this we see that the group velocity is affected by the presence of fluctuating motions. With these approximations in place, **c** can be separated as

where $c\xaf$ is the local mean ray velocity that depends on the mean properties of the medium,

and $c\xaf\u2032$ represents fluctuations in **c** due to the presence of small-scale surface flow fluctuations,

As a consequence, for small fluctuations, the ray propagation speed varies linearly with **u**′ and gradients of **u**′.

## 3. Ensemble statistics

To develop a statistical framework, we consider the properties of a large ensemble of wave packets that propagate through given mean flow superposed by random flow fluctuations. We express energy *E* and ray position **M** in terms of an ensemble mean and fluctuating component

where ⟨…⟩ denotes the ensemble average. Initially, the location and action density *N* of each ensemble member is known and set (without loss of generality) to *m* = *n* = *θ* = 0, with $\sigma \xaf=\omega $ and *σ*′ defined by the local $un\u2032$ as above. As a consequence, the initial *σ*′ is randomly distributed. Further, because *c*′ and *σ*′ are linearly dependent on $un\u2032,um\u2032$ (to the order considered) we find that $\u2329c\u232a=c\xaf$ and $\u2329\sigma \u232a=\sigma \xaf$. The objective in what follows is to predict the mean energy, the position, and the scatter around the mean (variance) for a large ensemble of rays (see, e.g., Fig. 2).

### a. Ensemble average ray path and energy density

To find the mean path ⟨**M**⟩(*t*) traced by the ensemble we take the ensemble average of the ray equations

We further assume that **M**′ does not induce *O*(*L*_{U}) changes in position, expressed as

so that in Eq. (17) we may approximate $\u2329c\xaf\u2061(M)\u232a=c\xaf\u2061(\u2329M\u232a)$. Further, invoking quasi-stationarity and quasi-homogeneity, and noting that velocity perturbations **u**′ (or gradients thereof) across the ensemble are uncorrelated with **M**′, we have $\u2329c\u2032\u2061(M)\u232a\u2248\u2329c\u2032\u2061(\u2329M\u232a)\u232a=0$. With these assumptions in place, and upon integration in time of Eq. (17), the average ray path traced out by the ensemble mean ⟨**M**⟩ can be written as

where $\Gamma =\omega /\sigma \xaf\u22121$ and |Γ| ≪ 1. Along-ray mean displacement ⟨*n*⟩ is dominated by displacement due to the wave motion proper. Propagation effects driven by the mean flow only become significant over long integration times. In what follows, given **U**, the mean properties in Eq. (18) are treated as known (e.g., through numerical integration).

The random fluctuations in the surface velocity field drive random fluctuations in the intrinsic frequency *σ* around its mean. Due to the conservation of wave action, *N* = *E*/*σ* such random variations in *σ* result in random variations in the energy density. The ensemble-averaged energy density ⟨*E*⟩ can be written as

which shows that along a wave ray, ⟨*E*⟩_{k} is unaffected by small-scale fluctuations. In contrast, the frequency-direction energy density $(E\sigma ,\theta =cg\u22121kEk)$ is slightly enhanced in the presence of such zero-mean fluctuations. Specifically, making use of the binomial theorem we find

However, since ⟨*σ*′⟩ = 0, and $\sigma \u2032/\sigma \xaf=O\u2061(\epsilon )$ (see appendix A) corrections only enter at *O*(*ε*^{2}) and this bias can be neglected.

### b. Variance and diffusivity estimates

To characterize the spreading around the mean location in each of the coordinates we consider the ensemble (co)variance defined as

where *α* and *β* denote *m*, *n*, *σ*, or *θ*. Due to the explicit dependence of *σ*′ on $un\u2032$, *V*_{σ′σ′} is constant in time,

As a consequence, since $u\u2032/C\xafg=O\u2061(\epsilon )$, the relative variance due to velocity fluctuations is *O*(*ε*^{2}), so that conservation of absolute frequency ensures that the intrinsic frequency remains constrained within limited bands around the mean.

For the remaining variables no such constrains apply. In general, individual realizations will follow random walks with growing variance in time. Specifically, since *u*′ decorrelates over spatial scales *L*_{u′}, and the waves propagate at a characteristic velocity *c*_{0}, we assume that velocity fluctuations decorrelate on time scale *τ*_{c} = *L*_{u′}/*c*_{0}, such that the correlation Γ_{αβ}(*t*, *τ*) ≈ 0 for *τ* > *τ*_{c}, with

The evolution of variances *V*_{αβ} follows from random walk theory. Specifically, since $\alpha \u2032=\u222b0tc\alpha \u2032\u2061(t1)dt1$, we have

where $t\xaf=t1$, *τ* = *t*_{1} − *t*_{2}, and where we used that Γ_{αβ} ≈ 0 for *τ* > *τ*_{c}. We assume that *t* ≫ *t*_{c}, so that the integration domain of the last two terms in (22) is $\u221d\tau c2$, and the first term on the right side integrates over 2*τ*_{c}*t*. If we assume that Γ_{αβ} integrated over *τ*_{c} is nonzero, the last two terms on the right of Eq. (22) can be neglected.^{2} Further, making use of $\Gamma \u2061(t\xaf,\tau )\u2248\Gamma \u2061(t\xaf,\u2212\tau )$ the rate of change of the variance is given by

where *D*_{αβ} is the diffusivity which may change on slow scales. For stationary statistics (*D*_{αβ} constant) the variance grows linearly in time (*V*_{αβ} = 2*D*_{αβ}*t*) akin to Brownian motion, and in general we find

where $D\u02dc\alpha \beta $ is the effective mean diffusion. Since $um\u2032un\u2032\xaf=0$ and $\u2202nun\u2032$ is in quadrature with $un\u2032$, nondiagonal contributions to the (mean) diffusion and covariance tensors vanish (*D*_{αβ} = *V*_{αβ} = 0 if *α* ≠ *β*). The diagonal entries *D*_{αα} can be nonzero, but are expressed in terms of the correlation in time of the fluctuations $c\alpha \u2032$ as observed in a coordinate system following a wave group. The leading order contribution to ray propagation occurs along *n* (see appendix A), so that

with $\xi n=g\tau /\u2061(2\sigma \xaf)$, *ξ*_{m} = 0, and

Invoking quasi-stationarity and quasi-homogeinity Γ_{α,α}(*t*, *τ*) can then be expressed in terms of the spatial correlation *γ*_{α,α} as

with

As a consequence, the diffusion coefficient can be expressed in terms of integration over the spatial correlation function *γ* as

with $\xi c=g\tau c/\u2061(2\sigma \xaf)$. Alternatively, for a quasi-homogeneous field (statistics do not change within spatial correlation scale) this may be approximated in terms of the local velocity wavenumber variance spectrum as

with *e*_{αα}(*λ*_{m}, *λ*_{n}) the variance density spectrum of velocity fluctuations $c\alpha \u2032$ as a function of the along-crest *λ*_{m} and along-ray *λ*_{n} wavenumbers. Given the definition of **c**′, and under the assumption that velocity fluctuations are isotropic such that $eum\u2032=eun\u2032=e/2$, with *e*(*λ*_{m}, *λ*_{n}) the total variance density, we find that

Expressed in terms of wavenumber magnitude *λ* and flow direction Θ (unrelated to the wave direction), and with the assumption that there is no preferential decorrelation direction, we have *e*(*λ*, Θ) = (2*π*)^{−1}*e*(*λ*), and since *e*(*λ*_{m}, *λ*_{n}) = *e*(*λ*, Θ)/*λ* we find after substitution that

### c. Geographical mean and variance in fixed coordinates

Predicting variance and mean values in a fixed coordinate frame **X** is complicated due to the nonlinear dependence of the leading order geographical ray velocities on wave orientation *θ*, as seen in

Estimates of spatial mean location and variance in fixed coordinates differ from estimates in wave-oriented coordinates. However, statistical estimates for *θ* and *σ* are unchanged. To estimate spatial mean and variance we neglect mean flows (*U*_{x} = *U*_{y} = 0), and assume $\theta \xaf=0$. Further, deviations from the mean are assumed small, and approximating sin*θ* and cos*θ* by truncated second-order series approximations then yields,

The angle *θ* depends on the integrated history over velocity gradients, but is not correlated to instantaneous velocities. As a consequence, ensemble-averaged velocities are

so that ensemble average velocity estimates depend on the directional variance *V*_{θθ}. This reduces the mean propagation speed of the ensemble. Physically, this speed reduction is due to path extension by the random fluctuation around a straight path. As a consequence, the mean ensemble location is affected by small scale fluctuations since

and ⟨*y*⟩ = 0. Similarly, spatial variances are influenced by directionality since (see appendix B for details)

Although initially spatial variance growth is dominated by *D*_{mm}, *D*_{nn}, as time progresses and directional variance of the field increases, the growth is dominated by additional spatial spreading due to enhanced directionality of the field. Further, for time scales such that the small-angle approximation remains valid (*t* < 1/*D*_{θθ}), we have that *V*_{yy} > *V*_{xx}, or lateral spreading dominates during initial evolution. Essentially, the process is comparable to shear-induced dispersion, since a small perturbation in the angle between two particles will result in large lateral differences in position due to difference in the lateral propagation velocity. Dispersion in the *y* direction dominates initially because we start along *θ* = 0 for which ∂_{θ}*c*_{x} ≈ 0.

## 4. The diffusive action balance equation

The Lagrangian wave-packet viewpoint is convenient to estimate diffusive rates, but is not easily amendable to calculating the evolution of swell fields. Our goal here is to describe the evolution of the mean Eulerian action density $N\xaf\u2061(x,y,\sigma ,\theta ,t)$ in the presence of random scattering. Assuming the action density is known at a time *t*, where *NJ* (with *J* = *C*_{g}/*k* the Jacobian) is conserved along rays, we consider the action integrated over a small volume *δV* = *δxδyδσδθ* at time *t* + *τ*. Here *τ* is a time interval long compared to the decorrelation time scale, but short compared to other time scales. Assuming that the statistics are quasi-homogeneous and quasi-stationary, the mean wave action in a small volume may be approximated as

During the time interval *τ* rays undergo displacement Δ(*τ*, **X**) such that **X** − Δ denotes the point of origin at time *t*. Consequently (similar to Einstein 1905), the expected value of the total action density *t* + *τ* is

where the displacement Δ during *τ* is given by

Since the displacements and time interval *τ* are assumed small, we approximate both sides of Eq. (34) by Taylor approximation, which results in

Substitution of the definition of Δ, dropping terms proportional to *τ*^{2} and higher, and using that $NJ\xaf=NJ$ then results in

Here *D*_{αβ} are the Lagrangian diffusivities. To leading order, the LHS of Eq. (35) represents an energy conservation equation (e.g., Tolman and Booij 1998). Using *D*_{σσ} ≈ 0, ∂_{θ}*J* = 0, and neglecting spatial gradients of *J* due to slow-scale variations in the diffusive terms we obtain

which is a conventional action balance augmented with diffusive terms on the right hand side. Further, the spatial diffusion coefficients are given by

with *D*_{mm} defined as in Eq. (27), whereas the directional diffusion coefficient *D*_{θθ} is given by Eq. (29). Note that the dependency of spatial diffusion on directional variance, which explicitly is accounted for in *D*_{xx} in the previous section, does not explicitly enter the spatial diffusion coefficients (due to assumed small *τ*). Instead, under the present Eulerian approximation, the effects of directional variance on spatial diffusion enter as a shear-dispersion effect.

## 5. Monte Carlo simulations

To verify the theoretical Lagrangian predictions, and to estimate the significance of small-scale flow effects on the wave statistics, we conduct Monte Carlo simulations for a large ensemble of wave rays propagating through a flow field with submesoscale variability. Specifically, following Callies et al. (2013, 2015), we parameterize the variance in the background flow *e*(*λ*) as

where *n* is the power of spectral decay and *e*_{0} is the kinetic energy density^{3} at reference wavenumber *λ*_{0}. Depending on season and location, the exponent *n* typically varies between 2 and 3 for oceanic surface flows, and using 10 km as a reference wavelength, *e*_{0} ranges from 1 to 10 m^{3} s^{−2}. We use *n* = 2.5 and set *e*_{0} = 1 or 10 m^{3} s^{−2}. Further, *λ*_{low} and *λ*_{up} denote the upper and lower wavenumbers considered. We assume that scales above 50 km are part of the mean flow (which we will neglect here) and use 1 km to define the upper wavenumber limit, to avoid excessively small integration steps in the ray tracing procedure.

To construct realizations of the flow field we assume that the flow statistics are isotropic. We use a discrete 2D spatial Fourier representation for $ux\u2032$ and $uy\u2032$ on a domain of 500 km × 500 km. For each realization the amplitude for each wavenumber component is set according to the variance contained in the binned interval around that wavenumber, and the phase is drawn randomly from a uniform distribution.

For a single realization of the flow field a set of 40 rays is initialized along *x* = 0 spaced at 12.5-km intervals along *y* and initialized with a given absolute radian frequency *ω* = 2*πf* and direction *θ* = 0. For each ray *j*, the evolution of **X**_{j} is subsequently calculated by numeric integration of the ray using Eqs. (2) with a fourth-order Runge–Kutta solver. For efficiency, use is made of the periodic nature of Fourier series approximation, and rays that leave the domain continue at the opposite boundary. We consider a total of 50 independent realizations with 40 rays each and ensemble-averaged statistics are based on 2000 rays.

### Results

The theoretical values for the standard deviation $S\alpha =V\alpha \alpha $ are in excellent agreement with estimates from Monte Carlo simulations, for a range of frequencies and different values of *e*_{0}. Here we specifically consider results for waves with periods of 10 and 20 s and *e*_{0} = 1 or 10 m^{3} s^{−2} (Fig. 3). Spatial spreading in *m* and *n* coordinates (*S*_{m} and *S*_{n}) typically grows to *O*(1) km over 5 days (Fig. 3, left and center panels). Although distances increase for large *e*_{0} and *f* (since *D*_{mm} ∝ *fe*_{0} the standard deviations *S*_{m} and *S*_{n} are proportional to $fe0t$), in general *S*_{m}, *S*_{n} are small relative to the mean distance traveled (a 10-s wave travels ~3000 km in this time) and the typical spatial size of a storm (~100 km). These findings imply that the direct spatial dispersion due to propagation of the waves by random flows is negligible.

The directional spreading *S*_{θ} is also $\u221dfe0t$, which implies that directional variance increases with increased flow variance and increased wave frequency. However *S*_{θ} varies from ~ 1° to 10° after 1 day (Fig. 3, right panel), so that even the lower limit is significant compared with observed directional spreading of ocean swell fields of *O*(10°).

In Earth-fixed coordinates, the increased pathlengths induce a spatial lag $C\xafgt\u2212x\xaf$ in the *x* coordinate along the mean wave direction. This spatial lag is frequency independent (results with different periods coincide), and scales as ∝ *e*_{0}*t*^{2}. After 5 days, the spatial lag grows to 10–100 km (Fig. 4, right panel). Although this is relatively small compared to the total travel distance this will have significant influence on swell arrival times (we return to this in the discussion).

Spreading along the mean wave direction *S*_{x} for small *t* is determined by the fluctuating flow (dependent on frequency). However, as directionality grows, shear-dispersion dominates variance in Eq. (33) and *S*_{x} becomes frequency independent and ∝ *e*_{0}*t*^{2} (Fig. 4, left panel). As a consequence *S*_{x} values are typically an order of magnitude larger, and comparable to the spatial lag in magnitude [in fact, $Sx\u22482\u2061(C\xafgt\u2212x\xaf)/3$ for large *t*]. However, spatial spreading on the order of days is typically dominated by lateral spreading. The *S*_{y} is $\u221dte0ft$ and consequently increases with increased variance of the flow or increased frequency, and after 1 day *S*_{y} ~ 20–100 km (Fig. 4, center panel).

## 6. Discussion

The statistical framework developed here describes the mean effects on a wave packet propagating over long distances through a field of submesoscale eddies. The predictions made by the model are in excellent agreement with direct Monte Carlo simulations. Since oceanic submesoscale turbulence is present everywhere, these effects are important for operational modeling of wave propagation over long distances. Direct empirical verification of the theoretical framework presented here can be obtained by observations of the directional evolution of ocean swell across the ocean along its route from generation area to remote coastline. Such observations are becoming more and more available from drifting wave and weather buoys. Once the theoretical predictions are empirically validated, the form of the directional scattering term in the radiative transfer equation can be immediately incorporated in operational models to account for scattering of waves on surface flow turbulence.

### a. Delayed swell arrival

Deviations from the great-circle path due to refraction in a random medium increases travel time for oceanic swells. Although the spatial lag is small compared to the total travel distance, the resulting delay in arrival times (compared with arrival times of swells in a quiescent ocean) is significant and easily observable. Specifically, from Eq. (32), we find that the expected delay Δ*t* = *t*_{0} − *t* is

where *t*_{0} is the travel time for a great-circle route (assuming that Δ*t*/*t*_{0} ≪ 1) and $x\xaf$ is the mean distance from the source to the observation point. Hence, expected delays rapidly increase as a function of increased source distance and frequency. The delay for swells with *f* > 0.05 that traveled ~3000 km can easily be in the order of hours (Fig. 5). In the context of operational forecasting applications this is a significant difference.

The order of magnitude of time delay and bias to late arrivals (travel times never decrease relative to the shortest path) is consistent with Jiang et al. (2016a) who reported mean arrival delays of several hours for swell at Pacific sites (compared with model results). While other factors may contribute, and the exact influence depends on the strength of oceanic surface flows, the present mechanism is a plausible explanation for delayed arrivals.

### b. Balance between geometric and diffusive spreading

Arguably, the most significant result of the present work is that within an Eulerian framework the effect of random scattering of waves can be represented by diffusive terms that can be directly related to the variance spectrum of the flow turbulence. Since spatial diffusivities are small, the leading contribution comes from the diffusion term in directional space, so that—augmented with source terms *S* representing generation, dissipation, and nonlinear interaction—the action balance becomes

where *D*_{θθ} is determined from Eq. (29). The inclusion of a diffusive term implies that the evolution of the directional shape of the spectrum will be influenced by a balance between geometric spreading, which decreases directional width, and the diffusive term, which increases directional width. Initially, near the source region the directional distribution is relatively broad and geometric spreading dominates the evolution. As the directional aperture of the wave field increasingly narrows, diffusion becomes more important and eventually arrests the geometrical narrowing, so that a stable directional shape develops.

To illustrate this, we consider an initial Gaussian wave field (period 15 s, cosine squared distribution, mean direction due north) with a 500-km-diameter footprint. We solve Eq. (38) numerically (finite differences) for a single frequency, in the absence of mean ambient flow [fourth and fifth term of LHS of Eq. (38) do not contribute] and sources and sinks (*S* = 0). The initial development of the wave field is similar to a reference solution without diffusion (Fig. 6). However, when including diffusion, after 2 days the narrowing of the distribution due to geometric spreading is arrested and a balance between diffusive scattering and geometric spreading is established. The resulting directional aperture is approximately 15° (Fig. 6, center panel), whereas the aperture of the reference solution (no diffusion) continues to narrow. Further, compared to the reference solution the wave height is significantly reduced after 10 days (by 50%). This simple modeling example shows that the scattering of surface waves on currents results in directional apertures in agreement with what is commonly observed (see Fig. 1). Further, these results are very similar to findings by Gallet and Young (2014), who used Monte Carlo ray tracing to estimate the probability density function of the far-field directional deflection angle. Moreover, the observed swell decay is weak but significant and may explain or contribute to widely observed, but unexplained, oceanic swell decay rates (see, e.g., Ardhuin et al. 2009; Young et al. 2013; Stopa et al. 2016).

## 7. Conclusions

We derived a stochastic model to account for the mean effects of near-surface submesoscale turbulent flows on the evolution of ocean surface wave fields. Specifically, we argued that due to the ubiquitous presence of submesoscale surface flow fields in the ocean, a wave packet’s trajectory deviates from a straight (great circle) propagation path, and—analogous to Brownian motion—this deviation takes the form of a random walk in geographical, intrinsic-frequency, and directional spaces. Based on Lagrangian and random walk theory we determined the ensemble mean and variance of the wavepacket position and frequency and direction for a wave ray bundle propagating through a random medium. We derived diffusion coefficients in each of the coordinate spaces, expressed in terms of the variance spectrum of the turbulent flows, and verified our theoretical results with Monte Carlo simulations. The scattering of surface waves on (sub)mesoscale turbulence can be accounted for in stochastic wave models by including a directional diffusion term in the action balance equation. Accounting for these subgrid refraction effects in operational forecast models will improve our ability to estimate swell arrivals, wave height decay (due to directional dispersion of swell fields) and the aperture of the directional distribution (due to a balance between geometric narrowing and diffusive spreading). Modeling errors in swell arrivals, wave heights, and—in particular—the directional aperture, have been widely observed, and some of these errors could likely be explained by accounting for submesoscale refraction effects. The changes in the wave statistics predicted by this theory are observable with standard oceanographic instruments, and direct empirical validation of these theoretical results is possible.

## Acknowledgments

This study was funded by the Office of Naval Research (N00014-16-1-2856). We thank David Clark for comments on a draft version of this work. The quality of this manuscript improved significantly from constructive feedback by two anonymous reviewers.

### APPENDIX A

#### Order *O*(*ε*) Approximations for *σ*′ and Group Velocity

The starting point are the previously introduced ray equations together with the along ray time variation of the absolute frequency

where the absolute frequency *ω* is related to the intrinsic frequency through the dispersion relation

The wave field is characterized by characteristic period *T*_{w}, wavelength *L*_{w}, (absolute) frequency *σ*_{0} = 2*π*/*T*_{w}, wavenumber *k*_{0} = *σ*_{0}/*c*_{0}, celerity *c*_{0} = *g*/*σ*_{0}, and we introduce the scaled variables

where $\u2026\u02dc$ denote scaled variables. Further, the flow and flow gradients are scaled as

with *u*_{0} the characteristic flow scale; *L*_{U} the spatial scale and *T*_{U} the time scale of the large scale flow; and *L*_{u′} and *T*_{u′} the spatial scale of the small scale flow. We define the scale parameters

which are all assumed to be small (≪ 1) so that the current speed is assumed small compared to the wave propagation and the characteristic scales of the current are large compared with the wave scale. Further, we introduce the scaling assumption

with *O*(*α*_{T}) = *O*(*α*_{L}) ≪ 1. This formalizes the assumption that small and large scales vary by at least an order of magnitude. With this scaling in place the ray equation and dispersion relation reduce to

with

and

At the leading order we thus have $\omega \u02dc=\sigma \u02dc$ and to *O*(*ε*) we have

We now decompose absolute and intrinsic frequencies into mean and fluctuating components as $\omega \u02dc=\omega \u02dc\xaf+\omega \u02dc\u2032$ and $\sigma \u02dc=\sigma \u02dc\xaf+\sigma \u02dc\u2032$. We further introduce the scaling assumption $\omega \u02dc\u2032=O\u2061(\epsilon \mu u\u2032)$ so that $\omega \u02dc\u2032$ can be neglected to *O*(*ε*) (we motivate this assumption below).

The average $\sigma \xaf$ and fluctuating part of the intrinsic frequency *σ*′ now follow immediately from taking the average of Eq. (A12), and subtracting the result from Eq. (A12). Noting that to the leading order $\sigma \u02dc\xaf=\omega \u02dc$, and since accurate up to *O*(*ε*) we have $\omega \u02dc=\omega \u02dc\xaf$, we then obtain

To the order considered changes in absolute frequency are thus associated with changes in mean intrinsic frequency, whereas $\sigma \u02dc\u2032$ balances with the Doppler shift due to local velocity fluctuations.

From Eq. (A13) we find that $\u2061(\sigma \u2032\u02dc)\u22121\sigma \u02dc\xaf=O\u2061(\epsilon )$ and we can thus approximate the propagation velocity $\u2061(2\sigma \u02dc)\u22121$ with a Taylor series approximation, which at the first order reduces to

As a consequence, in the leading order approximation particle displacement only occurs in the normal direction, with $cn=\u2061(2\sigma \u02dc\xaf)\u22121+O\u2061(\epsilon )$.

##### a. Order of $\omega \u02dc\u2032$

To obtain the evolution equation for $\omega \u02dc\u2032$ we take the mean of Eq. (A9), and subsequently subtract it from Eq. (A9). To the leading order this results in

As a consequence $\omega \u02dc\u2032=O\u2061(\epsilon \mu u\u2032)$ on short time scales, and since we only consider up to *O*(*ε*) in Eq. (A12) may be neglected on *O*(1) time scales. Further, due to fluctuating nature of $u\u2032\u02dc$, growth is not monotonic; instead $\sigma \u02dc\u2032$ oscillates and a mean drift only accumulates slowly. To estimate the time scale when $\sigma \u02dc\u2032$ becomes *O*(*ε*), we assume (analogous with the other variables in the main text) random walk for the perturbed variable so that

with $S\omega \u02dc\u2032$ the standard deviation and $D\omega \u02dc\u2032\omega \u02dc\u2032$ a coefficient characterizing the diffusive strength. If the initial scale $S\omega \u2032$ is order *O*(*εμ*) this implies that $S\omega \u2032$ becomes *O*(*ε*) after $t\u02dc\u2248O\u2061(\mu \u22122)$. Setting *T*_{w} = *O*(10) s and *T*_{u′} = *O*(10^{4}) we find that *μ*^{−2} represents *O*(10^{6}) wave periods, or *O*(100) days. As a consequence, we may neglect the $\omega \u02dc\u2032$ contribution on scales of *O*(10) days that are of interest for swell propagation.

### APPENDIX B

#### Estimation of *V*_{xx} and *V*_{yy}

To estimate *V*_{xx} and *V*_{yy} in Eq. (33) we have to approximate

where we made use of symmetry of the integrands about *t*_{1} = *t*_{2} so that now *t*_{2} ≤ *t*_{1}. To approximate *V*_{yy} (which we will make use of in approximating *V*_{xx} below) we note that $\u2329cy\u2032\u2061(t1)cy\u2032\u2061(t2)\u232a=C\xafg2\u2329\theta \u2061(t1)\theta \u2061(t2)\u232a+\u2061(9/4)\u2329ux\u2032\u2061(t1)ux\u2032\u2061(t2)\u232a$ [neglecting in Eq. (30) assumed small contributions originating from products between *u*′ and *θ*′], and assume quasi-stationary statistics so that

where we introduced the coordinate transformation *τ* = *t*_{3} − *t*_{4} and $t\xaf=t4$ and integrated over $t\xaf$ assuming stationary statistics. Assuming that *t* is much larger than the coherence length scale, we may approximate the upper bounds of the integrands over *τ* with the coherence length scale *τ*_{c} (thus neglecting contributions to *V*_{yy} when |*t*_{1} − *t*_{2}| ≤ *τ*_{c}, *t*_{1} ≤ *τ*_{c} and *t*_{2} ≤ *τ*_{c}). And since we then have *τ* ≪ *t*_{1} and *τ* ≪ *t*_{2} we find

Substitution into *V*_{yy} of this approximation we then find after integration that

To approximate *V*_{xx}, we substitute $cx\u2032=cx\u2212c\xafx$ so that

If the statistics are (and remain) Gaussian, we can neglect fourth-order cumulant contributions to the fourth-order moment,

Consequently, with Eq. (B1) we have $\u2329cx\u2032\u2061(t1)cx\u2032\u2061(t2)\u232a=2D\theta \theta 2t22$, so that after substitution and integration of the results we find

## REFERENCES

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

More restrictive definitions of submesoscale based on flow dynamics exist in the literature. Since only flow kinematics will enter our discussion, we use a scale-based (literal) definition.

^{3}

The factor *π*^{−1} results from a change from inverse wavelength density as in Callies et al. (2015) to wavenumber density, and transforming from kinetic energy to variance.