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.

Fig. 1.

(a) Directional width at a buoy site near San Diego (CDIP site 191) due to waves originating from a hypothetical storm south of New Zealand traveling due east. The storm is parameterized as an isotropic radiator traveling eastward at 30 km h−1 with the locally generated sea represented by a JONSWAP spectrum (peak at 14 s and 8-m significant wave height). Directional widths at CDIP site 191 are estimated through great-circle ray tracing from the moving storm. (b) Histogram for the directional width of swell arriving at the same site. The site is in deep water (1143 m), and we consider 10 years’ worth of data (October 2007–January 2018, in total 39 832 events). To select relevant swell fields only and exclude more local fields, which would bias the spreading high, only spectral peaks below 0.1 Hz with mean directions that originate from the Southern or Indian Ocean are considered. Directional properties are calculated according to Kuik et al. (1988) 

Fig. 1.

(a) Directional width at a buoy site near San Diego (CDIP site 191) due to waves originating from a hypothetical storm south of New Zealand traveling due east. The storm is parameterized as an isotropic radiator traveling eastward at 30 km h−1 with the locally generated sea represented by a JONSWAP spectrum (peak at 14 s and 8-m significant wave height). Directional widths at CDIP site 191 are estimated through great-circle ray tracing from the moving storm. (b) Histogram for the directional width of swell arriving at the same site. The site is in deep water (1143 m), and we consider 10 years’ worth of data (October 2007–January 2018, in total 39 832 events). To select relevant swell fields only and exclude more local fields, which would bias the spreading high, only spectral peaks below 0.1 Hz with mean directions that originate from the Southern or Indian Ocean are considered. Directional properties are calculated according to Kuik et al. (1988) 

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) eddies1 (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 = [ux, uy]T, with ux, uy 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,

 
ddt(Ekσ)=0.
(1)

Here Ek(t) = E[k(t), x(t), t] is the energy density of a wave packet with wavenumber k(t) = [kx(t), ky(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

 
dXdt=c(X,t),
(2)

with

 
c(X,t)=[cxcycθcσ]=g2σ[cosθ(t)sinθ(t)00]+[uxuymun12σnun].
(3)

Here un[x(t), θ(t), t], um[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,

 
[un(x,θ,t)um(x,θ,t)]=A[ux(x,t)uy(x,t)],[nm]=A[xy],
(4)

with A the rotation matrix for rotation angle θ defined as

 
A(θ)=[cosθsinθsinθcosθ].
(5)

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

 
u=U+u.

In operational models, typically only the large-scale flow field that varies on time and space scales TU and LU (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 Lu and Tu (with LuLU and TuTU). We further assume that the underlying random process is quasi-stationary and quasi-homogeneous, so that the variance of u′, changes on the TU and LU mesoscales. Further, we assume that waves travel fast relative to submesoscale time and space scales, expressed as

 
LuCgTu,
(6)

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,

 
[n(t)m(t)]=[n(t0)m(t0)]+t0tA(t)[cx(X(t),t)cy(X(t),t)]dt,
(7)

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

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

 
dMdt=c[M(t),t],
(8)

with

 
c=[cmcncθcσ]=[0g(2σ)100]+[UmUnmUn12σnUn]+[umunmun12σnun],
(9)

and where Um, Un, and um,un 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 Cg = 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

 
ω=σ+σ2g(Un+un).
(10)

We decompose absolute and relative frequencies as σ=σ¯+σ and ω=ω¯+ω, 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 ωω¯ so that ω is directly tied to the large scale flow and mean intrinsic frequency as

 
σ¯=ω(σ¯)2gUn.
(11)

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

 
σ(σ¯)2gun.
(12)

Since σ/σ¯=O(ε) we can now approximate Cg as a first order Taylor series, so that the along-ray velocity cn=Cg+Un+un becomes

 
cng2σ¯gσ2(σ¯)2+Un+un=C¯g+Un+32un,

with C¯g=g(2σ¯)1. 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

 
c=c¯+c,
(13)

where c¯ is the local mean ray velocity that depends on the mean properties of the medium,

 
c¯=[Um,C¯g+Un,mUn,12σ¯nUn]T,
(14)

and c¯ represents fluctuations in c due to the presence of small-scale surface flow fluctuations,

 
(c)=[um,3un2,mun,12σ¯nun]T.
(15)

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

 
Ek=Ek+Ek,M=M+M,
(16)

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 σ¯=ω and σ′ defined by the local un as above. As a consequence, the initial σ′ is randomly distributed. Further, because c′ and σ′ are linearly dependent on un,um (to the order considered) we find that c=c¯ and σ=σ¯. 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).

Fig. 2.

(left) Ray paths through a turbulent medium for a wavefield with initially (along x = 0) parallel (eastward directed) waves. The background flow, calculated as described below in section 5, is characterized with the vorticity (positive red, negative blue). Due to the irregular refraction individual rays deviate from a straight line, and the deviation takes the form of a random walk. (right) For example, the ensemble mean wave direction remains approximately 0, but the wave direction for an individual ray (e.g., black line) develops erratically.

Fig. 2.

(left) Ray paths through a turbulent medium for a wavefield with initially (along x = 0) parallel (eastward directed) waves. The background flow, calculated as described below in section 5, is characterized with the vorticity (positive red, negative blue). Due to the irregular refraction individual rays deviate from a straight line, and the deviation takes the form of a random walk. (right) For example, the ensemble mean wave direction remains approximately 0, but the wave direction for an individual ray (e.g., black line) develops erratically.

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

 
dMdt=c¯(M)+c(M).
(17)

We further assume that M′ does not induce O(LU) changes in position, expressed as

 
|x(M,t)x(M,t)|LU1,

so that in Eq. (17) we may approximate c¯(M)=c¯(M). 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 c(M)c(M)=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

 
m=0tUmdt,n=g2ω[t0tΓdt]+0tUndt,θ=0tUnmdt,
(18)

where Γ=ω/σ¯1 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

 
Ek=Nk(σ¯+σ)=Nkσ¯,

which shows that along a wave ray, ⟨Ek is unaffected by small-scale fluctuations. In contrast, the frequency-direction energy density (Eσ,θ=cg1kEk) is slightly enhanced in the presence of such zero-mean fluctuations. Specifically, making use of the binomial theorem we find

 
Eσ,θ=2σ¯4g2Nk(σ¯+σ)4σ¯4=2σ¯4g2Nk[1+4σσ¯+6(σσ¯)2+4(σσ¯)3+(σσ¯)4]

However, since ⟨σ′⟩ = 0, and σ/σ¯=O(ε) (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

 
Vαβ=(αα)(ββ),
(19)

where α and β denote m, n, σ, or θ. Due to the explicit dependence of σ′ on un, Vσσ is constant in time,

 
Vσσσ¯2=(un)24C¯g2.
(20)

As a consequence, since u/C¯g=O(ε), 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 Lu, and the waves propagate at a characteristic velocity c0, we assume that velocity fluctuations decorrelate on time scale τc = Lu/c0, such that the correlation Γαβ(t, τ) ≈ 0 for τ > τc, with

 
Γαβ(t,τ)=cα[M(t),t]cβ[M(t+τ),t+τ].
(21)

The evolution of variances Vαβ follows from random walk theory. Specifically, since α=0tcα(t1)dt1, we have

 
Vαβ=0t0tcα(t1)cβ(t2)dt1dt2=0ttt¯t¯Γαβ(t¯,τ)dτdt¯=0tτcτcΓαβ(t¯,τ)dτdt¯0τct¯τcΓαβ(t¯,τ)dτdt¯tτctt¯t0Γαβ(t¯,τ)dτdt¯,
(22)

where t¯=t1, τ = t1t2, and where we used that Γαβ ≈ 0 for τ > τc. We assume that ttc, so that the integration domain of the last two terms in (22) is τc2, and the first term on the right side integrates over 2τct. 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 Γ(t¯,τ)Γ(t¯,τ) the rate of change of the variance is given by

 
dVαβdt=2Dαβ(t),Dαβ(t)=0τcΓαβ(t,τ)dτ,t>τc,
(23)

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

 
Vαβ=2D˜αβt,D˜αβ=0tDαβ(t)dtt,
(24)

where D˜αβ is the effective mean diffusion. Since umun¯=0 and nun is in quadrature with un, 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α 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

 
M(t+τ)M(t)+ΔM[ξm(t),ξn(t)],
(25)

with ξn=gτ/(2σ¯), ξm = 0, and

 
ΔM(ξm,ξn)=[ξm,ξn,0,0]T.

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

 
Γαα(t,τ)γαα(0,gτ2σ¯,M,t)

with

 
γαα(ξm,ξn,M,t)=cα(M¯,t)cα[M(t)+ΔM(ξm,ξn),t].

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

 
Dαα=2σ¯g0ξcγαα(0,ξ,M¯,t)dξ

with ξc=gτc/(2σ¯). 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

 
Dαα=2πσ¯geαα(λm,0,M,t)dλm,
(26)

with eαα(λm, λn) the variance density spectrum of velocity fluctuations cα 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=eun=e/2, with e(λm, λn) the total variance density, we find that

 
em=e2,en=98e,eθ=λm2e2.

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π)−1e(λ), and since e(λm, λn) = e(λ, Θ)/λ we find after substitution that

 
Dmm=σ¯g0e(λ)λdλ,
(27)
 
Dnn=94Dmm,and
(28)
 
Dθθ=σ¯g0λe(λ)dλ.
(29)

The relations (26)(28) express the wave ray diffusivities in terms of the submesoscale velocity variance spectrum and they represent the main results of this section.

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

 
[cxcy]=Cg[cosθsinθ]+[1+12cos2θ12cosθsinθ12sinθcosθ1+12sin2θ]u+U.

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 (Ux = Uy = 0), and assume θ¯=0. Further, deviations from the mean are assumed small, and approximating sinθ and cosθ by truncated second-order series approximations then yields,

 
[cxcy]=Cg[112(θ)2θ]+[3214(θ)212θ12θ1+12(θ)2]u+U.
(30)

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

 
[c¯xc¯y]=C¯g[112Vθθ0],
(31)

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

 
x=C¯g(t12Dθθt2),
(32)

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

 
[VxxVyy]=0t0t[cx(t1)cx(t2)cy(t1)cy(t2)]dt1dt22t[DnnDmm]+C¯g2[13Dθθ2t423Dθθt3].
(33)

Although initially spatial variance growth is dominated by Dmm, Dnn, 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 Vyy > Vxx, 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 ∂θcx ≈ 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¯(x,y,σ,θ,t) in the presence of random scattering. Assuming the action density is known at a time t, where NJ (with J = Cg/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

 
δVJN(X,t+τ)¯dxdydσdθJN(X,t+τ)¯δV.

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

 
J(X,t+τ)N(X,t+τ)¯δV=N[XΔ(τ),t]J[XΔ(τ),t]¯δV,
(34)

where the displacement Δ during τ is given by

 
Δ(t,τ)=tt+τ[c(t)+c(t)]dtc(t)τ+dc(t)dtτ22++X(τ).

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

 
τNJt¯+τ222NJt2¯+=ΔαNJXα¯+ΔαΔβ22NJXαXβ¯+.

Substitution of the definition of Δ, dropping terms proportional to τ2 and higher, and using that NJ¯=NJ then results in

 
NJt+c¯αNJxα=Dαβ2NJxαxβ.
(35)

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

 
Nt+c¯xNx+c¯yNy+c¯θNθ+c¯σNσ=Dxx2Nx2+2Dxy2Nxy+Dyy2Ny2+Dθθ2Nθ2,
(36)

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

 
Dxx=Dmm(1+54cos2θ),
 
Dyy=Dmm(1+54sin2θ),
 
Dxy=Dmm(54sinθcosθ),

with Dmm 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 Dxx 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

 
e(λ)={e0π(λλ0)n,λlowλλup,0,elsewhere,
(37)

where n is the power of spectral decay and e0 is the kinetic energy density3 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, e0 ranges from 1 to 10 m3 s−2. We use n = 2.5 and set e0 = 1 or 10 m3 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 and uy 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 Xj 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α=Vαα are in excellent agreement with estimates from Monte Carlo simulations, for a range of frequencies and different values of e0. Here we specifically consider results for waves with periods of 10 and 20 s and e0 = 1 or 10 m3 s−2 (Fig. 3). Spatial spreading in m and n coordinates (Sm and Sn) typically grows to O(1) km over 5 days (Fig. 3, left and center panels). Although distances increase for large e0 and f (since Dmmfe0 the standard deviations Sm and Sn are proportional to fe0t), in general Sm, Sn 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.

Fig. 3.

Evolution of standard deviation (left) Sm, (center) Sn, and (right) Sθ as a function of time. Considered waves have periods of 10 and 20 s and the variance spectrum e(λ) of the flow is parameterized with e0 = 10 m3 s−2 (upper curves) and e0 = 1 m3 s−2 (lower curves).

Fig. 3.

Evolution of standard deviation (left) Sm, (center) Sn, and (right) Sθ as a function of time. Considered waves have periods of 10 and 20 s and the variance spectrum e(λ) of the flow is parameterized with e0 = 10 m3 s−2 (upper curves) and e0 = 1 m3 s−2 (lower curves).

The directional spreading Sθ is also fe0t, 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¯gtx¯ in the x coordinate along the mean wave direction. This spatial lag is frequency independent (results with different periods coincide), and scales as ∝ e0t2. 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).

Fig. 4.

Evolution of standard deviation (left) Sx, (center) Sy, and (right) the diffusion induced lag C¯gtx¯ as a function of time. Considered waves have periods of 10 and 20 s and the variance spectrum e(λ) of the flow is parameterized with e0 = 10 m3 s−2 (upper curves) and e0 = 1 m3 s−2 (lower curves).

Fig. 4.

Evolution of standard deviation (left) Sx, (center) Sy, and (right) the diffusion induced lag C¯gtx¯ as a function of time. Considered waves have periods of 10 and 20 s and the variance spectrum e(λ) of the flow is parameterized with e0 = 10 m3 s−2 (upper curves) and e0 = 1 m3 s−2 (lower curves).

Spreading along the mean wave direction Sx for small t is determined by the fluctuating flow (dependent on frequency). However, as directionality grows, shear-dispersion dominates variance in Eq. (33) and Sx becomes frequency independent and ∝ e0t2 (Fig. 4, left panel). As a consequence Sx values are typically an order of magnitude larger, and comparable to the spatial lag in magnitude [in fact, Sx2(C¯gtx¯)/3 for large t]. However, spatial spreading on the order of days is typically dominated by lateral spreading. The Sy is te0ft and consequently increases with increased variance of the flow or increased frequency, and after 1 day Sy ~ 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 = t0t is

 
Δt=12Dθθx2C¯ge0f2x2,

where t0 is the travel time for a great-circle route (assuming that Δt/t0 ≪ 1) and x¯ 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.

Fig. 5.

Contours of the expected time lag (h) for swell propagating through a turbulent flow field. Flow field parameterized as in Eq. (37) with n = 2.5 and e0 = 10 m3 s2.

Fig. 5.

Contours of the expected time lag (h) for swell propagating through a turbulent flow field. Flow field parameterized as in Eq. (37) with n = 2.5 and e0 = 10 m3 s2.

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

 
Nt+c¯xNx+c¯yNy+c¯θNθ+c¯σNσ=Dθθ2Nθ2+S,
(38)

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).

Fig. 6.

Evolution of an ocean wave field with 15-s period according to Eq. (38) propagating at the group velocity in the mean direction of the initial swell field. For the initial wave field, we assume a spatially Gaussian distribution (500-km width) and a cosine squared directional distribution. (a) Directional spreading with (black) and without (gray) the inclusion of diffusive scattering for increasingly energetic submesoscale motions [parameterized with Eq. (37)]. (b) Comparison between directional spreading after 5 days calculated with and without diffusive term in. (c) Wave heights as a function of time with (black) and without (gray) diffusive term (the latter is indistinguishable from least energetic diffusive solution).

Fig. 6.

Evolution of an ocean wave field with 15-s period according to Eq. (38) propagating at the group velocity in the mean direction of the initial swell field. For the initial wave field, we assume a spatially Gaussian distribution (500-km width) and a cosine squared directional distribution. (a) Directional spreading with (black) and without (gray) the inclusion of diffusive scattering for increasingly energetic submesoscale motions [parameterized with Eq. (37)]. (b) Comparison between directional spreading after 5 days calculated with and without diffusive term in. (c) Wave heights as a function of time with (black) and without (gray) diffusive term (the latter is indistinguishable from least energetic diffusive solution).

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

 
dMdt=c,dωdt=σ2g(cosθtux+sinθtuy),
(A1)

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

 
ω=σ+σ2gun.
(A2)

The wave field is characterized by characteristic period Tw, wavelength Lw, (absolute) frequency σ0 = 2π/Tw, wavenumber k0 = σ0/c0, celerity c0 = g/σ0, and we introduce the scaled variables

 
t=t˜Tw,m=m˜Lw,n=n˜Lw,θ=θ˜,σ=σ0σ˜,ω=σ0ω˜,
(A3)

where ˜ denote scaled variables. Further, the flow and flow gradients are scaled as

 
Un=u0Un˜,Um=u0Um˜,mUn=u0LUmUn˜,nUn=u0LUnUn˜,
(A4)
 
un=u0un˜,um=u0um˜,mun=u0Lumun˜,nun=u0Lunun˜,
(A5)
 
tUx=u0TUtUx˜,tUy=u0TUtUy˜,tux=u0Tutux˜,tux=u0Tutux˜,
(A6)

with u0 the characteristic flow scale; LU the spatial scale and TU the time scale of the large scale flow; and Lu and Tu the spatial scale of the small scale flow. We define the scale parameters

 
ε=u0c0,λu=LwLu,λU=LwLU,μu=TwTu,μU=TwTU,
(A7)

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

 
TuTU=αT,LuLU=αL,
(A8)

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

 
dM˜dt˜=c˜dω˜dt˜=εμu(σ˜)2(cosθtux˜+sinθtuy˜)+αTεμu(σ˜)2(cosθtUx˜+sinθtUy˜),
(A9)

with

 
c˜=[0(2σ˜)100]+ε[U˜mU˜nαLλumUn˜12αLλuσ˜nUn˜]+ε[u˜mu˜nλumun˜12λuσ˜nun˜],
(A10)

and

 
ω˜=σ˜+εσ˜2(U˜+u˜).
(A11)

At the leading order we thus have ω˜=σ˜ and to O(ε) we have

 
σ˜=ω˜εω˜2(U˜+u˜).
(A12)

We now decompose absolute and intrinsic frequencies into mean and fluctuating components as ω˜=ω˜¯+ω˜ and σ˜=σ˜¯+σ˜. We further introduce the scaling assumption ω˜=O(εμu) so that ω˜ can be neglected to O(ε) (we motivate this assumption below).

The average σ¯ 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 σ˜¯=ω˜, and since accurate up to O(ε) we have ω˜=ω˜¯, we then obtain

 
σ˜¯=ω˜εσ˜¯2U˜σ˜=εσ˜¯2u˜.
(A13)

To the order considered changes in absolute frequency are thus associated with changes in mean intrinsic frequency, whereas σ˜ balances with the Doppler shift due to local velocity fluctuations.

From Eq. (A13) we find that (σ˜)1σ˜¯=O(ε) and we can thus approximate the propagation velocity (2σ˜)1 with a Taylor series approximation, which at the first order reduces to

 
12σ˜12σ˜¯σ˜2σ2˜¯12σ˜¯+εu˜2.

As a consequence, in the leading order approximation particle displacement only occurs in the normal direction, with cn=(2σ˜¯)1+O(ε).

a. Order of ω˜

To obtain the evolution equation for ω˜ we take the mean of Eq. (A9), and subsequently subtract it from Eq. (A9). To the leading order this results in

 
dω˜dt˜=εμuσ˜¯2(cosθtux˜+sinθtuy˜).
(A14)

As a consequence ω˜=O(εμu) 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˜, growth is not monotonic; instead σ˜ oscillates and a mean drift only accumulates slowly. To estimate the time scale when σ˜ becomes O(ε), we assume (analogous with the other variables in the main text) random walk for the perturbed variable so that

 
Sω˜=2Dω˜ω˜tt,
(A15)

with Sω˜ the standard deviation and Dω˜ω˜ a coefficient characterizing the diffusive strength. If the initial scale Sω is order O(εμ) this implies that Sω becomes O(ε) after t˜O(μ2). Setting Tw = O(10) s and Tu = O(104) we find that μ−2 represents O(106) wave periods, or O(100) days. As a consequence, we may neglect the ω˜ contribution on scales of O(10) days that are of interest for swell propagation.

APPENDIX B

Estimation of Vxx and Vyy

To estimate Vxx and Vyy in Eq. (33) we have to approximate

 
[VxxVyy]=2C¯g20t0t1[cx(t1)cx(t2)cy(t1)cy(t2)]dt2dt1,

where we made use of symmetry of the integrands about t1 = t2 so that now t2t1. To approximate Vyy (which we will make use of in approximating Vxx below) we note that cy(t1)cy(t2)=C¯g2θ(t1)θ(t2)+(9/4)ux(t1)ux(t2) [neglecting in Eq. (30) assumed small contributions originating from products between u′ and θ′], and assume quasi-stationary statistics so that

 
θ(t1)θ(t2)=0t10t2cθ(t3)cθ(t4)dt3dt4=0t1(t1τ)Γcθcθ(t,τ)dτ+0t2(t2τ)Γcθcθ(t,τ)dτ0t1t2(t1t2)Γcθcθ(t,τ)dτ,

where we introduced the coordinate transformation τ = t3t4 and t¯=t4 and integrated over t¯ 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 Vyy when |t1t2| ≤ τc, t1τc and t2τc). And since we then have τt1 and τt2 we find

 
θ(t1)θ(t2)2t2Dθθ.
(B1)

Substitution into Vyy of this approximation we then find after integration that

 
Vyy=23Cg2Dθθt3+2tDnn.
(B2)

To approximate Vxx, we substitute cx=cxc¯x so that

 
cx(t1)cx(t2)=14Cg¯2[θ(t1)]2[θ(t2)]2C¯g4Vθθ(t1)Vθθ(t2)+uy(t1)uy(t2).

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

 
[θ(t1)]2[θ(t2)]2Vθθ(t1)Vθθ(t2)+2[θ(t1)θ(t2)]2.

Consequently, with Eq. (B1) we have cx(t1)cx(t2)=2Dθθ2t22, so that after substitution and integration of the results we find

 
Vxx=13Cg2Dθθ2t4+2tDmm.
(B3)

REFERENCES

REFERENCES
Andrews
,
D. G.
, and
M. E.
Mcintyre
,
1978
:
An exact theory of nonlinear waves on a Lagrangian-mean flow
.
J. Fluid Mech.
,
89
,
609
646
, https://doi.org/10.1017/S0022112078002773.
Ardhuin
,
F.
,
N.
Rascle
, and
K. A.
Belibassakis
,
2008
:
Explicit wave-averaged primitive equations using a generalized Lagrangian mean
.
Ocean Modell.
,
20
,
35
60
, https://doi.org/10.1016/j.ocemod.2007.07.001.
Ardhuin
,
F.
,
B.
Chapron
, and
F.
Collard
,
2009
:
Observation of swell dissipation across oceans
.
Geophys. Res. Lett.
,
36
,
L06607
, https://doi.org/10.1029/2008GL037030.
Ardhuin
,
F.
,
S. T.
Gille
,
D.
Menemenlis
,
C. B.
Rocha
,
N.
Rascle
,
B.
Chapron
,
J.
Gula
, and
J.
Molemaker
,
2017
:
Small-scale open ocean currents have large effects on wind wave heights
.
J. Geophys. Res. Oceans
,
122
,
4500
4517
, https://doi.org/10.1002/2016JC012413.
Barber
,
N. F.
, and
F.
Ursell
,
1948
:
The generation and propagation of ocean waves and swell. I. Wave periods and velocities
.
Philos. Trans. Roy. Soc. London
,
240A
,
527
560
, https://doi.org/10.1098/rsta.1948.0005.
Battjes
,
J.
,
1974
:
Computation of set-up, longshore currents, run-up and overtopping due to wind-generated waves. Ph.D. thesis, Delft University of Technology, 251 pp
.
Bretherton
,
F. P.
, and
C. J. R.
Garrett
,
1968
:
Wavetrains in inhomogeneous moving media
.
Proc. Royal Soc. London
,
302A
,
529
554
, https://doi.org/10.1098/rspa.1968.0034.
Callies
,
J.
,
R.
Ferrari
,
J.
Callies
, and
R.
Ferrari
,
2013
:
Interpreting energy and tracer spectra of upper-ocean turbulence in the submesoscale range (1–200 km)
.
J. Phys. Oceanogr.
,
43
,
2456
2474
, https://doi.org/10.1175/JPO-D-13-063.1.
Callies
,
J.
,
R.
Ferrari
,
J. M.
Klymak
, and
J.
Gula
,
2015
:
Seasonality in submesoscale turbulence
.
Nat. Commun.
,
6
,
6862
, https://doi.org/10.1038/ncomms7862.
Cavaleri
,
L.
,
B.
Fox-Kemper
, and
M.
Hemer
,
2012
:
Wind waves in the coupled climate system
.
Bull. Amer. Meteor. Soc.
,
93
,
1651
1661
, https://doi.org/10.1175/BAMS-D-11-00170.1.
Einstein
,
A.
,
1905
:
On the motion of small particles suspended in a stationary liquid, as required by the molecular kinetic theory of heat
.
Ann. Phys.
,
322
,
549
560
, https://doi.org/10.1002/andp.19053220806.
Ewans
,
K. C.
,
2002
:
Directional spreading in ocean swell
.
Fourth Int. Symp. on Ocean Wave Measurement and Analysis
,
San Francisco, CA, American Society of Civil Engineers, 517–529
, https://doi.org/10.1061/40604(273)54.
Fabrikant
,
A. L.
, and
M. A.
Raevsky
,
1994
:
The influence of drift flow turbulence on surface gravity wave propagation
.
J. Fluid Mech.
,
262
,
141
156
, https://doi.org/10.1017/S0022112094000455.
Gallet
,
B.
, and
W. R.
Young
,
2014
:
Refraction of swell by surface currents
.
J. Mar. Res.
,
72
,
105
126
, https://doi.org/10.1357/002224014813758959.
Heller
,
E. J.
,
L.
Kaplan
, and
A.
Dahlen
,
2008
:
Refraction of a Gaussian seaway
.
J. Geophys. Res.
,
113
,
C09023
, https://doi.org/10.1029/2008JC004748.
Janssen
,
T. T.
, and
T. H. C.
Herbers
,
2009
:
Nonlinear wave statistics in a focal zone
.
J. Phys. Oceanogr.
,
39
,
1948
1964
, https://doi.org/10.1175/2009JPO4124.1.
Jiang
,
H.
,
A. V.
Babanin
,
G.
Chen
,
H.
Jiang
,
A. V.
Babanin
, and
G.
Chen
,
2016a
:
Event-based validation of swell arrival time
.
J. Phys. Oceanogr.
,
46
,
3563
3569
, https://doi.org/10.1175/JPO-D-16-0208.1.
Jiang
,
H.
,
J. E.
Stopa
,
H.
Wang
,
R.
Husson
,
A.
Mouche
,
B.
Chapron
, and
G.
Chen
,
2016b
:
Tracking the attenuation and nonbreaking dissipation of swells using altimeters
.
J. Geophys. Res. Oceans
,
12
,
1446
1458
, https://doi.org/10.1002/2015JC011536.
Kuik
,
A. J.
,
G. P.
Van Vledder
, and
L. H.
Holthuijsen
,
1988
:
A method for the routine analysis of pitch-and-roll buoy wave data
.
J. Phys. Oceanogr.
,
18
,
1020
1034
, https://doi.org/10.1175/1520-0485(1988)018<1020:AMFTRA>2.0.CO;2.
Lentz
,
S. J.
,
M.
Fewings
,
P.
Howd
,
J.
Fredericks
, and
K.
Hathaway
,
2008
:
Observations and a model of undertow over the inner continental shelf
.
J. Phys. Oceanogr.
,
38
,
2341
2357
, https://doi.org/10.1175/2008JPO3986.1.
Longuet-Higgins
,
M. S.
,
1970
:
Longshore currents generated by obliquely incident sea waves, 1
.
J. Geophys. Res.
,
75
,
6778
6789
, https://doi.org/10.1029/JC075i033p06778.
MacMahan
,
J. H.
,
E. B.
Thornton
, and
A. J. H. M.
Reniers
,
2006
:
Rip current review
.
Coastal Eng.
,
53
,
191
208
, https://doi.org/10.1016/j.coastaleng.2005.10.009.
McWilliams
,
J. C.
,
J. M.
Restrepo
, and
E. M.
Lane
,
2004
:
An asymptotic theory for the interaction of waves and currents in coastal waters
.
J. Fluid Mech.
,
511
,
135
178
, https://doi.org/10.1017/S0022112004009358.
Mei
,
C. C.
,
M.
Stiassnie
, and
D. K. P.
Yue
,
2005
:
Theory and Applications of Ocean Surface Waves
.
Advanced Series on Ocean Engineering, Vol. 23, World Scientific, 1136 pp
.
Metzger
,
J. J.
,
R.
Fleischmann
, and
T.
Geisel
,
2014
:
Statistics of extreme waves in random media
.
Phys. Rev. Lett.
,
112
, 203903, https://doi.org/10.1103/PhysRevLett.112.203903.
Munk
,
W. H.
, and
M. A.
Traylor
,
1947
:
Refraction of ocean waves: A process linking underwater topography to beach erosion
.
J. Geol.
,
55
,
1
26
, https://doi.org/10.1086/625388.
Munk
,
W. H.
,
G. R.
Miller
,
F. E.
Snodgrass
, and
N. F.
Barber
,
1963
:
Directional recording of swell from distant storms
.
Philos. Trans. Roy. Soc. London
,
255A
,
505
584
, https://doi.org/10.1098/rsta.1963.0011.
O’Reilly
,
W. C.
,
T. H. C.
Herbers
,
R. J.
Seymour
, and
R. T.
Guza
,
1996
:
A comparison of directional buoy and fixed platform measurements of Pacific swell
.
J. Atmos. Oceanic Technol.
,
13
,
231
238
, https://doi.org/10.1175/1520-0426(1996)013<0231:ACODBA>2.0.CO;2.
Smit
,
P. B.
, and
T. T.
Janssen
,
2013
:
The evolution of inhomogeneous wave statistics through a variable medium
.
J. Phys. Oceanogr.
,
43
,
1741
1758
, https://doi.org/10.1175/JPO-D-13-046.1.
Snodgrass
,
F. E.
,
G. W.
Groves
,
K. F.
Hasselmann
,
G. R.
Miller
,
W. H.
Munk
, and
W. H.
Powers
,
1966
:
Propagation of ocean swell across the Pacific
.
Philos. Trans. Roy. Soc. London
,
259A
,
431
497
, https://doi.org/10.1098/rsta.1966.0022.
Stopa
,
J. E.
,
F.
Ardhuin
,
R.
Husson
,
H.
Jiang
,
B.
Chapron
, and
F.
Collard
,
2016
:
Swell dissipation from 10 years of Envisat advanced synthetic aperture radar in wave mode
.
Geophys. Res. Lett.
,
43
,
3423
3430
, https://doi.org/10.1002/2015GL067566.
Sullivan
,
P. P.
, and
J. C.
McWilliams
,
2010
:
Dynamics of winds and currents coupled to surface waves
.
Annu. Rev. Fluid Mech.
,
42
,
19
42
, https://doi.org/10.1146/annurev-fluid-121108-145541.
Tolman
,
H. L.
, and
N.
Booij
,
1998
:
Modeling wind waves using wavenumber-direction spectra and a variable wavenumber grid
.
The Global Atmosphere and Ocean System
,
Vol. 6, Overseas Publishers Association, 295–309
, http://polar.ncep.noaa.gov/mmab/papers/tn162/OMB162.pdf.
Xu
,
Z.
, and
A. J.
Bowen
,
1994
:
Wave- and wind-driven flow in water of finite depth
.
J. Phys. Oceanogr.
,
24
,
1850
1866
, https://doi.org/10.1175/1520-0485(1994)024<1850:WAWDFI>2.0.CO;2.
Young
,
I. R.
,
A. V.
Babanin
, and
S.
Zieger
,
2013
:
The decay rate of ocean swell observed by altimeter
.
J. Phys. Oceanogr.
,
43
,
2322
2333
, https://doi.org/10.1175/JPO-D-13-083.1.

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.

2

Note that for Γσσ we would find that the first term in Eq. (22) vanishes and all contributions to the variance originate from the latter terms. Specifically, Eq. (26) would integrate to 0 in that case.

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.