Introduction

Surface gravity waves are an important route by which the ocean exchanges energy, momentum, heat and gases with the overlying atmosphere (Cavaleri, Fox-Kemper & Hemer 2012;Villas Bôas et al. 2019). Sea-surface currents modify the wavenumber, direction and amplitude of surface waves, and affect the spatial variability of the wave field. The effect of currents on waves under the Wentzel-Kramers-Brillouin (WKB) approximation has been well studied (Kenyon 1971;Peregrine 1976;White & Fornberg 1998;Henderson et al. 2006;Heller, Kaplan & Dahlen 2008;Gallet & Young 2014). But the sparseness of ocean current observations makes it difficult to explicitly account for wave-current interactions in numerical surface-wave models. Thus, instead of explicit resolving sea-surface currents, a statistical approach to the effect of currents on surface waves is required. † Email address for correspondence: avillasboas@ucsd.edu A. B. Villas Bôas and W. R. Young

Recent studies of surface wave-current interactions suggest that the sea-state variability at meso-and submesoscales, here referred to as macroturbulence, is dominated by the variability of the current field (Ardhuin et al. 2017;Quilfen et al. 2018;Quilfen & Chapron 2019). At these scales, horizontally divergent motions associated with tides, inertia-gravity waves, and fronts contribute significantly to the surface kinetic energy (Bühler, Callies & Ferrari 2014;Rocha et al. 2016a;D'Asaro et al. 2018). If surface gravity waves respond differently to divergent than to rotational flows -and we show here that they do -then changes in the dominant regime of surface currents can result in significant changes in the surface-wave field. This offers the possibility that observations of surface gravity waves might be used to probe the structure of submesoscale ocean turbulence.

In the context of internal gravity waves, McComas & Bretherton (1977) showed how scale-separated wave interactions can be analysed with the WKB approximation and understood as diffusion of wave action. The induced diffusion approximation of McComas and Bretherton has recently been developed and extended by Kafiabad, Savva & Vanneste (2019, KSV, hereafter) to obtain an action-diffusion equation for the scattering of internal gravity waves by mesoscale ocean turbulence. Here we apply the KSV method to surface gravity waves. Crucial to this development is that the parameter def = |U|/c

(1.1) is small; above, U is the horizontal current at the sea surface and c = √ g/4k is the deep-water group speed at wavenumber k.

In § 2, and in appendix A, we use the formalism of KSV to derive an expression for a diffusivity tensor of surface-wave action. In § 3 we consider the simplifications that result from assuming that the sea-surface velocity U has isotropic statistics. We show that the horizontally divergent component of U has no effect on action diffusivity: diffusivity results solely from the vortical (solenoidal) component of U and produces an angular diffusivity that is expressed as a weighted integral of the solenoidal part of the energy spectrum of U as in (3.20). Smit & Janssen (2019) have also examined the action diffusion of surface waves using a framework based on Lagrangian random walk theory. Section 3 discusses the differences between Smit and Janssen's expression for the action diffusivity and ours. In § 4 the analytic results are tested with Monte Carlo ray tracing through an ensemble of stochastic velocity fields.

The induced diffusion approximation

For linear deep-water surface waves, the Doppler-shifted dispersion relation is

where k = (k 1 , k 2 ) is the wavenumber, σ = √ gk is the intrinsic wave frequency, with k = |k| and g the gravitational acceleration. Also in (2.1), U(t, x) = (U 1 , U 2 ) is the horizontal current at the sea surface. Provided that U(t, x) is slowly varying with respect to the waves (i.e. the temporal scales of variations in the current field are longer and the spatial scales are larger than those of the waves), wave kinematics is described by the ray equations. Using index notation the ray equations are

where c n = ∂ k n σ (k) is the group velocity. Under the same assumptions, wave dynamics is governed by the conservation of wave action density A(x, k, t)

with ẋn and kn given by (2.2) (Phillips 1966;Mei 1989). We follow KSV and develop a multiple-scale solution, based on 1, that enables one to average (2.3) over the ensemble of velocity fields U (see appendix A). Assuming that the statistical properties of U are stationary and homogeneous, one finds that

where Ā denotes the ensemble average of A. The diffusivity tensor D jn in (2.4) is expressed in terms of the two-point velocity correlation tensor

(2.5)

Because of the assumption of spatial homogeneity, V im depends only on the separation r = xx of the two points. The most convenient formula for explicit calculation of D jn is the Fourier space result

where c = g/2σ is the magnitude of the group velocity and

is the Fourier transform of V im (r). (In (2.6) and (2.7) the integrals cover the entire two-dimensional planes (q 1 , q 2 ) and (r 1 , r 2 ) respectively.) The diffusivity in (2.6) is the two-dimensional equivalent of (A7) in KSV. Our appendix A derivation, however, assumes only spatial homogeneity and stationarity of the velocity U, and does not require incompressibility of U.

One can verify from (2.6) that D jn k n = 0 and therefore there is no diffusion of wave action in the radial direction in k-space. Fast surface-wave packets propagate through a frozen field of macroturbulent eddies and thus preserve the absolute frequency

Thus, at leading order, both σ and k, are constant. In other words, absolute frequency conservation, together with 1, implies that there is no radial k-diffusion in (2.4). Thus, scattering by weak surface currents results mainly in directional diffusion of surface gravity waves.

Diffusion of wave action density by isotropic velocity fields

The derivation of (2.6) makes essential use of the assumption that the spatial statistics of U are spatially homogeneous. We now make the further assumption that the statistical properties of U are also isotropic and investigate the contributions of vertical vorticity and horizontal divergence to D jn . We follow Bühler et al. (2014) and represent U with a two-dimensional Helmholtz decomposition into rotational (solenoidal) and irrotational (potential) components

(3.1)

A. B. Villas Bôas and W. R. Young

The streamfunction ψ and velocity potential φ have the two-point correlation functions

If the velocity ensemble is not mirror invariant under reflection with respect to an axis in the (x, y)-plane, then there might also be a 'cross-correlation' between ψ and φ:

Because of isotropy, the scalar correlation functions introduced in (3.2) and(3.3) depend only on the distance r = |r| between points x and x . Therefore,

Using the notation r = (r 1 , r 2 ), the V 11 component of the velocity autocorrelation tensor in (2.5) can be expressed in terms of the scalar correlation functions as

Similar calculations for the other components of V im result in

with

and

(3.9)

The Fourier transform of (3.8) and (3.9) follows with ∂ r i → iq i and is equal to

where q ⊥ = (-q 2 , q 1 ) is the perpendicular vector to q = (q 1 , q 2 ). Also in (3.10) and (3.11)

(3.12) with J 0 the Bessel function of order zero, is the Fourier transform of the axisymmetric function C ψ (r). The expressions for Cφ (q) and Cψφ (q) are analogous to (3.12). Substituting (3.10) and (3.11) into (2.6) we have

where • • • above indicates the three other terms that arise from contracting (3.10) and (3.11) with k i k m . Each of these three terms, however, contains a factor k • q. Courtesy Diffusion of surface waves of δ(k • q) in the integrand of (2.6), the . . . in (3.13) makes no contribution to D jn and the diffusivity tensor reduces to

where the integral covers the entire (q 1 , q 2 )-plane. The diffusion tensor in (3.14) does not depend on the velocity potential φ. Using δ(q • k) to evaluate one of the two integrals in (3.14) one obtains

It is remarkable that the compressible and irrotational component of the velocity field, produced by the velocity potential φ, makes no contribution to the actiondiffusion tensor in (3.15). Dysthe (2001) shows that in the weak-current limit, 1, the ray curvature is equal to ζ /c, where ζ = ψ xx + ψ yy is the vertical vorticity of the surface currents; see § 68 of Landau & Lifshitz (1987) and Gallet & Young (2014) for alternative derivations. These ray-tracing results rationalize the result in (3.15) that diffusion of surface-wave action by sea-surface currents is produced only by the vortical and horizontally incompressible component of the sea-surface velocity.

This effect is illustrated in figure 1, where we show ray trajectories obtained by numerical integration of the ray equations (2.2) for waves with a period of 10 s propagating through three different types of surface flows (purely solenoidal, purely potential, and combined solenoidal and potential). These synthetic surface currents were created from a scalar function with random phase and prescribed spectral slope (q -2.5 in this case). In panel (a) this function is used as a streamfunction ψ to generate an incompressible vortical flow. In panel (b) the same function is used as a velocity potential φ to generate an irrotational horizontally divergent flow. In panels (b) and (e), with pure potential flow, the ray trajectories are close to straight lines (i.e. there is almost no scattering). The flow in panel (c) is constructed by summing the velocity fields in panels (a) and (b). Even though the flow in panel (c) is twice as energetic as that in panel (a), the ray trajectories in panels (d) and ( f ) are very similar. This is a striking confirmation of (3.14): the diffusivity is not affected by φ.

Because D jn k n = 0, the diffusive flux of wave action, -D jn ∂ k n Ā, is in the direction of k ⊥ = k θ , where (k, θ ) are polar coordinates in the k-plane and θ is a unit vector in the θ-direction. Using these polar coordinates simplifies the ∂ k j and ∂ k n derivatives on the right of (2.4) so that the averaged action equation becomes

is the directional diffusivity.

To conclude this section we express α in (3.17) in terms of the energy spectrum of the solenoidal component of the velocity Ẽψ (q), related to Cψ (q) by Ẽψ (q) = q 3 4π Cψ (q).

(3.18) The spectrum is normalized so that the root-mean-square velocity of the solenoidal component, U ψ , is

Then α(k) can be written as

Taking the trace of the velocity correlation tensors in (3.10) and (3.11) shows that the total energy spectrum is Ẽ = Ẽψ (q) + Ẽφ (q), (3.21)

where Ẽφ (q) is obtained by ψ → φ in (3.19). As anticipated in figure 1, the diffusivity α(k) in (3.20) depends only on the spectrum of the solenoidal component, Ẽψ (q).

Smit & Janssen (2019) arrive at an expression for α(k) differing from (3.20) in two respects: (i) the coefficient in front of the integral on the right is (1/c); and (ii) the integrand is Ẽ(q). This expression agrees with (3.20) only for the special class of isotropic velocity fields considered by Smit and Janssen in which the integral of Ẽφ (q) is equal to the integral of Ẽψ (q) (i.e. isotropic flows in which kinetic energy is equipartitioned between the solenoidal, ψ, and the potential, φ, components). An example of an equipartitioned flow is shown in figure 1(c,f ) and discussed further in § 4 (see the + simulations). Equipartition, however, is not characteristic of ocean macroturbulence (for example, large scales are in geostrophic balance and are therefore solenoidal). In this pure solenoidal case the diffusivity in Smit & Janssen (2019) would be too small by a factor of two. where denotes the action-weighted average and cos θ 0 is the initial value of cos θ . At large times cos θ → 0 with an e-folding time α -1 : this is long-time isotropization of the wave field by eddy scattering. To investigate short-time and small-angle scattering, consider for simplicity an initial condition such as that in figure 1 with initial direction θ 0 = 0. Then with αt 1 and θ 1, it follows from (4.3) that θ 2 ≈ 2αt. To test our result for the diffusivity α, we verify θ 2 ≈ 2αt by numerical integration of the ray-tracing equations (2.2) for surfaces waves with an initial period of 10 s propagating through an ensemble of stochastic velocity fields. The ensemble is created by assigning random phases to each Fourier component of the stream function ψ and velocity potential φ.

The energy spectrum of the sea-surface velocity is modelled with power laws Ẽψ (q) and Ẽφ (q) ∝ q -n , with q 1 < q < q 2 and no energy outside the interval (q 1 , q 2 ). The spectra are normalized with prescribed mean square velocities U 2 ψ and U 2 φ as in (3.19). For n = (1, 2) the integral in (3.20) is evaluated as: ≈ 2αt (solid lines). Here we show the results for an energy spectrum with spectral slopes following a q -n power law where n = 5/3, 2, 2.5, or 3. Circles • are the result for solenoidal flows; diamonds , for potential flows; and crosses + for the combination of solenoidal and potential. The solenoidal and potential flows have mean square velocity 0.01 m 2 s -2 , whereas the combined flow + has mean square velocity 0.02 m 2 s -2 . The initial period and direction of the waves are 10 s and 0 • , respectively. and for n = 2 α = 2 c q 1 q 2 q 2 -q 1 ln(q 2 /q 1 )U 2 ψ .

(4.6)

We take q 1 = 2π/150 km and q 2 = 2π/1 km and spectral slopes n = (5/3, 2.0, 2.5, 3.0).

For each n we consider three cases corresponding to the three columns in figure 1:

Figure 2 summarizes the results by showing θ 2 as a function of time obtained by averaging 2000 rays. The results are in agreement with θ 2 ≈ 2αt using α obtained from (4.4) and (4.6). In particular, there is good agreement between θ 2 for case • and the analytic result (solid lines). As expected, the potential component of the velocity has no effect on the diffusion of wave action. Thus, in case -pure potential flow -there is no diffusion of action. In case + the flow has twice as much kinetic energy (and shear) as in cases • and ; this is also an example of a flow with kinetic energy equipartitioned between the solenoidal and potential components (Smit & Janssen 2019). Doubling the strength of the flow, by adding a φ component, does not significantly increase action diffusion above that of case •.

At the final time, two days, the n = 5/3 simulations shown in figure 2 have θ 2 of order 30 • and the other, steeper, spectral slopes result in smaller directional spreading. Thus, none of the Monte Carlo simulations shown in figure 2 have lasted long enough to result in isotropization of the wave field. We verified, however, that longer simulations are in agreement with (4.3) when αt ∼ 1 (not shown).

We conclude this section by noting that numerical and observational evidence supports the hypothesis that Ẽ(q) ∼ q -2 on submesoscales (very roughly, scales less than 50 km). Horizontally divergent motions contribute significantly the total surface kinetic energy in this range, but the solenoidal component is not negligible (Rocha et al. 2016b;Torres et al. 2018;Kafiabad, Savva & Vanneste 2019;Morrow et al. 2019). There is considerable geographic variation. For example, the Gulf Stream region is an exception, with Ẽ(q) ∼ q -3 and little indication of horizontally divergent motions (Bühler et al. 2014). These results indicate that spectral slope -2 is relevant to oceanic application of (3.20). But with -2, the integral on the right of (3.20) is sensitive to high-wavenumber solenoidal energy (i.e. to the value of the high-wavenumber cutoff q 2 ), as in (4.6). This problem is worse for spectral slopes shallower than -2, and less severe in the Gulf Stream region with the steeper slope -3.

In the absence of a high-wavenumber transition to a spectral fall-off steeper than -2, the cutoff q 2 might be determined by the failure of the WKB approximation once the horizontal scales of U are comparable to the hundred-metre wavelengths of surface gravity waves. These considerations complicate the practical application of (3.20) in some oceanic regimes, and indicate the necessity of better understanding the interaction of surface gravity waves with wave-scale currents.

Conclusions

Our expression for the action diffusivity in (2.6) assumes that the WKB approximation is valid and that = |U|/c 1. Typical sea-surface currents are of order 0.1 m s -1 , while the swell band has group velocities that exceed 5 m s -1 . Thus, 1 is not restrictive. Our analysis also neglects effects associated with vertical shear of the flow, which would modify the Doppler-shifted dispersion relationship (Kirby & Chen 1989).

We derived an expression for the diffusivity of surface-wave action in (2.6) and demonstrated that for isotropic surface currents the action diffusivity can be expressed in terms of the kinetic energy spectrum of the flow as in (3.20). This result shows that the potential component makes no contribution to action diffusion. Our results are illustrated both qualitatively (figure 1) and quantitatively (figure 2) by numerical solution of the ray equations. Although the numerical examples presented here were obtained for synthetic flows having random phase, the results are also valid in the presence of coherent structures, such as axisymmetric vortices, as long as the statistics remain isotropic (not shown). To leading order, there is no difference between the diffusivity obtained for rays propagating through a pure solenoidal flow and the same solenoidal flow with the addition of an equally strong potential component. Provided that 1, the horizontally divergent and irrotational component of the sea-surface velocity has no effect on the action diffusion of surface gravity waves.

Recent studies motivated by the upcoming Surface Water and Ocean Topography (SWOT) satellite mission have found that surface kinetic energy spectra in the ocean are marked by a transition scale from balanced geostrophic motions (horizontally nondivergent) to unbalanced horizontally divergent motions such as inertia-gravity waves (for example, Rocha et al. 2016a,b;Torres et al. 2018;Qiu et al. 2018;Morrow et al. 2019). At scales shorter than this 'transition' scale, the kinetic energy spectrum of the potential component of the currents has been observed to dominate over the solenoidal component. In this regime, only a small fraction of the total kinetic energy of the flow would be contributing to the diffusion of surface-wave action.

Perhaps the most important application of our results is in the realm of operational surface-wave models. Wave models, such as WaveWatch III, solve the action balance equation (2.3) with additional terms to account for wind forcing, nonlinear interactions and wave dissipation (Wavewatch III Development Group 2009). Explicitly solving for wave-current interactions in surface-wave models poses two main challenges: it is computationally costly and surface current observations at scales shorter than 100 km are rare (Ardhuin et al. 2012). The wave action diffusivity calculated here can be easily implemented as an additional term in operational wave models, allowing the effects of the currents to be accounted for based on statistical properties of the sea-surface velocity. Although not discussed in the present manuscript, it is also worth noting that refraction of surface waves by meso-and submesoscale flows will ultimately lead to deviations of the wave propagation from the great-circle route, impacting path lengths and, subsequently, arrival times (Smit & Janssen 2019). A statistical approach to account for these effects in numerical wave models could potentially improve arrival-time predictions.