## Abstract

During its periodic motion, a particle floating at the free surface of a water wave experiences a net drift velocity in the direction of wave propagation, known as the Stokes drift (Stokes 1847 *Trans. Camb. Philos. Soc.* **8**, 441–455). More generally, the Stokes drift velocity is the difference between the average Lagrangian flow velocity of a fluid parcel and the average Eulerian flow velocity of the fluid. This paper reviews progress in fundamental and applied research on the induced mean flow associated with surface gravity waves since the first description of the Stokes drift, now 170 years ago. After briefly reviewing the fundamental physical processes, most of which have been established for decades, the review addresses progress in laboratory and field observations of the Stokes drift. Despite more than a century of experimental studies, laboratory studies of the mean circulation set up by waves in a laboratory flume remain somewhat contentious. In the field, rapid advances are expected due to increasingly small and cheap sensors and transmitters, making widespread use of small surface-following drifters possible. We also discuss remote sensing of the Stokes drift from high-frequency radar. Finally, the paper discusses the three main areas of application of the Stokes drift: in the coastal zone, in Eulerian models of the upper ocean layer and in the modelling of tracer transport, such as oil and plastic pollution. Future climate models will probably involve full coupling of ocean and atmosphere systems, in which the wave model provides consistent forcing on the ocean surface boundary layer. Together with the advent of new space-borne instruments that can measure surface Stokes drift, such models hold the promise of quantifying the impact of wave effects on the global atmosphere–ocean system and hopefully contribute to improved climate projections.

This article is part of the theme issue ‘Nonlinear water waves’.

## 1. Introduction

Although their leading-order motion is periodic—in other words backwards and forwards—surface gravity waves induce a net drift in the direction of wave propagation known as the Stokes drift [1]. Following the motion of particles underneath a wave, ‘in addition to the motion of oscillation the particles are transferred forwards, that is, in the direction of propagation, with a constant velocity’ [1], p. 207. The linearized trajectories of Lagrangian particles underneath linear unidirectional surface gravity waves are formed by closed ellipses, tending to circles in the limit when the water depth is large relative to the wavelength, as is illustrated in figure 1. A fluid particle, which oscillates backwards and forwards due to the linear wave motion, spends more time in the forward-moving region underneath the crest than in the backward-moving region underneath the trough and undergoes its forward motion at greater height, where the velocities are larger. As a result of these two effects, the particle experiences a net forward drift, which is proportional to the square of the steepness of the waves. Although other types of waves can display Stokes drift, such as vertically confined internal wave modes [4–6], oceanic Kelvin and Rossby waves [7–10] or acoustic waves [11], the focus of this review is on surface gravity waves.

Broadly, three categories of applications of Stokes drift can be distinguished. First, Stokes drift plays a role in wave-induced sediment transport and sandbar migration in the coastal zone, where it drives an opposite return flux, often in the form of an undertow, when it meets the no-flow boundary condition imposed by a beach. Second, Stokes drift plays an important role in explaining Langmuir turbulence, namely the formation of a series of shallow, slow, counter-rotating vortices near the ocean's surface with characteristic bands of floating seaweed, foam and debris accumulating in the convergence zones between the vortices, at scales of 2 m–1 km. The magnitude of the Stokes drift velocity and its rate of shear with depth is important for the inclusion of Langmuir turbulence and Coriolis–Stokes forcing in Eulerian ocean models. Third, in combination with Eulerian currents, driven by winds, density gradients and tides, Stokes drift transports heat, salt and other natural or man-made tracers including micro-plastic pollution in the upper ocean layer.

This paper reviews progress in fundamental and applied research on the induced mean flow associated with surface gravity waves since the first description of Stokes drift, now 170 years ago. The authors are unaware of previous comprehensive reviews of the fluid mechanics of Stokes drift and it applications, but refer the reader to insightful (mathematical) treatments of waves and mean flows by Craik [12], ch. 4, and the broader context of Lagrangian-mean theory by Bühler [13], ch. 10. Although it is our aim herein to be comprehensive albeit brief, we emphasize that our review is by no means complete, as the applications of Stokes drift are diverse and the literature addressing each application significant. In particular, there is a rich and varied literature on other wave effects that require inclusion in Eulerian ocean and near shore hydrodynamic models. Most conspicuously, we do not review the rich topic of wave breaking and its associated injection of kinetic energy and momentum [14–16] and we point the interested reader to the review by Sullivan & McWilliams [17]. Nor do we attempt to review in detail how wave shoaling and breaking along a shore will modify water levels and generate along- and across-shore currents.

This paper is laid out as follows. First, §2 introduces the fundamental physical processes behind Stokes drift for surface gravity waves. Specifically, §2a reviews the basic derivation of Stokes drift for periodic waves, followed by a discussion of the additional Eulerian return flow for wave groups in §2b. Then, §2c reviews how rotation can modify this picture and how Stokes drift gives rise to the Coriolis–Stokes force. Particle diffusion by Stokes drift and Eulerian wave-induced flows in stochastic sea states is discussed in §2d. Second, §3 reviews laboratory studies of Stokes drift, followed by a discussion of field measurements in §4, distinguishing measurements from buoys and drifters (§4a) and radar (§4b). The subsequent three sections discuss three main areas of application of Stokes drift: its role in formation of undertow in the coastal zone (§5); its integration in non-wave-resolving Eulerian ocean models to account for Coriolis–Stokes forcing and the production of Langmuir turbulence (§6); and its ability to transport tracers, such as oil and plastic pollution and in search and rescue operations (§7). Finally, conclusions are drawn and future directions indicated in §8.

## 2. Fundamental physical processes

### (a) Periodic waves: Stokes drift

Stokes drift is conceptually most straightforwardly defined as the difference between Lagrangian and Eulerian averages of a flow field (cf. [13]):
2.1
Stokes drift thus corresponds to the difference in wave-averaged velocity following a particle (Lagrangian) and in a stationary reference frame (Eulerian). The term ‘Stokes drift’ is reserved to such differences in velocity, whereas the term ‘Stokes correction’ can be applied more generally to the difference between Lagrangian and Eulerian averages of other quantities (e.g. pressure). Stokes drift is a property of the wave only and, to leading-order in wave steepness, it can then be calculated as the net velocity that results from small displacements of a Lagrangian particle during its phase cycle:
2.2
where *u*^{(1)} is the velocity field of the linear wave and *ξ*^{(1)} is the corresponding linear displacement vector with components , which can be evaluated from the linear velocity field according to ∂*ξ*^{(1)}/∂*t*=*u*^{(1)}. The overline denotes averaging over the linear waves, and the superscripts (1) correspond to first order in steepness.

For simplicity, we consider periodic deep water surface gravity waves, where the term deep water refers to the water depth being large relative to the wavelength, a suitable approximation for most waves in the ocean. Satisfying the Laplace equation ∇^{2}*ϕ*^{(1)}=0 for an incompressible and irrotational fluid and the two linearized free surface boundary conditions, we then have and , where *ϕ*^{(1)} and *η*^{(1)} are the (linear) velocity potential and the free surface displacement (), respectively. The parameters *a*, *ω* and *k* denote the wave amplitude, frequency and wavenumber, and *x* and *z* are the horizontal and vertical coordinates, respectively. After solving ∂*ξ*^{(1)}/∂*t*=*u*^{(1)} for the linear displacement vector *ξ*^{(1)}, substitution of *u*^{(1)} and *ξ*^{(1)} into (2.2) gives
2.3
where *c*=*ω*/*k* is the phase speed, as first shown by Stokes [1] for deep water (*kh*≫1). Ursell [18] extended (2.3) to general water depth and obtained , where *h* is the water depth. Of the two contributions to the Stokes drift for a unidirectional wave in (2.2), the contribution from the vertical displacement is unaffected by finite depth apart from a change in the depth structure ( versus ), whereas the contribution from the horizontal displacement from the increasingly elliptical orbits becomes more and more significant for shallower depth. Longuet-Higgins [19] showed that the error in (2.3) is small: *O*((*ak*)^{6}). The non-closed trajectories underneath linear Stokes waves are computed explicitly in Constantin & Villari [20] (and for shallow-water and deep-water linear waves by Ionescu-Kruse [21] and Constantin *et al*. [22], respectively) and underneath nonlinear Stokes wave in Constantin [23] (and for deep-water nonlinear waves by Henry [24]), with supporting numerical evidence from a boundary integral formulation provided in Nachbin & Ribeiro-Junior [25]. For a spectrum consisting of waves of arbitrary direction and wavenumber, Kenyon [26] showed that the Stokes drift is given by
2.4
where *k*=|** k**| is the magnitude of the wavenumber vector,

*g*is the gravitational acceleration and

*F*(

**) is the energy (wave variance) spectrum in wavenumber coordinates.**

*k*Stokes drift is closely related to mean wave (pseudo-) momentum [27]. Returning to periodic, deep-water, unidirectional waves, the mean wave momentum per unit area (in the *x*-direction) is defined as *ρM*, with *ρ* denoting constant density and *M* given by
2.5
where *T*=2*π*/*ω* is the wave period and we only retain leading-order terms. From (2.3), we have
2.6
so that Stokes drift can be interpreted as the vertical distribution of mean wave momentum [28,29]. The quantity *M* also corresponds to the depth-integrated Stokes drift [30], known as the Stokes transport, and is often estimated in oceanographic applications.

### (b) Wave groups: Stokes drift and Eulerian return flow

In reality, the wave field on the open sea often has a group-like structure (e.g. [31]), resulting in a Stokes transport flux that is divergent on the scale of the group and must induce another flow at second order in steepness: the return flow [32]. This return flow is Eulerian and, together with the Stokes drift, it makes up the total Lagrangian induced mean flow for wave groups (cf. (2.1)). Although the return flow is commonly explained as driven by a gradient in radiation stress [33], it can also be explained as the irrotational response to balance the Stokes transport, which is divergent on the scale of the group and acts to ‘pump’ fluid from its trailing edge to its leading edge (e.g. [34]). Provided the water is deep enough, a spatial separation of these two aspects of the induced mean flow takes place: the Stokes transport dominates near the free surface over the e-folding depth (2*k*_{0})^{−1}, with *k*_{0} now denoting the peak of the wavenumber spectrum; the magnitude of the return flow decreases much more slowly with depth, on the scale of the group, and consequently dominates far below the free surface, as illustrated in figure 2. Combining the (local) Stokes transport and the (non-local) return flow leads to zero vertically integrated mass transport at the centre of the group, and hence there is zero vertically integrated momentum associated with the centre of a surface gravity wave group, as emphasized by McIntyre [27]. The return flow is accompanied by a depression in the wave-averaged surface elevation on the scale of the wave group, often referred to as a set-down [32]. The set-down itself does not affect the return flow in deep water [27], but acts to enhance it for intermediate and small water depth (e.g. [34]).

For the spectrum of linear waves representing a wave group, the harmonic components interact to give both ‘frequency-sum’ and ‘frequency-difference’ terms, as first described by Longuet-Higgins & Stewart [32] for unidirectional waves and to second-order in steepness. These ‘frequency-difference’ terms describe the return flow. Phillips [35] and Hasselmann [36], and separately in the offshore engineering literature Sharma & Dean [37], Dalzell [38] and Forristall [39], extended these results to multidirectional seas, allowing for interactions between parent-wave components of different frequencies and travelling in different directions. In the narrow-bandwidth limit of a single carrier wave travelling in one direction, differential equations describing these second-order bound interactions can also be calculated using a multiple-scales approach. Dysthe [40], for infinite water depth, and Davey & Stewartson [41], for finite water depth, show that the nonlinear evolution equations for the wave group are accompanied by a second set of differential equations describing the mean flow and the wave-averaged free surface.

Very recently, Haney & Young [42] have shown that a surface gravity wave group on a stratified fluid generates a trailing wake of internal gravity waves in addition to a slightly compressed return flow. Using parameters typical of short period surface wave swell (8 s period and 1 m amplitude), these authors find that the energy flux between the two types of waves is small, and the coupling between surface waves and internal waves is not a significant sink of energy for the surface waves nor a source for internal waves. In a more extreme case (20 s period and 4 m amplitude), this coupling becomes a significant source of energy for internal waves with frequencies close to the buoyancy frequency.

### (c) Rotation: Gerstner waves and the Coriolis–Stokes force

Long before Stokes [1] published his approximate irrotational solutions to the water wave equations, Gerstner [43] identified an exact solution that was obtained through consideration of the Lagrangian equations of motion, a solution that was later rediscovered by Froude [44], Rankine [45] and Reech [46]. So-called Gerstner waves are rotational, leading Lamb [47] to conclude that they cannot be generated by a system of conservative forces. As a result, Gerstner waves have received little attention in the literature, although a number of authors have recently considered whether Gerstner waves may explain observations of mean drift in wave flume experiments [48,49]. We will return to this discussion in §3.

Examining the role of rotation in the mean flow in a periodic wave train, Ursell & Deacon [50] argued that there can only be zero time-averaged Lagrangian drift in a rotating frame, because an unopposed Stokes drift would violate the conservation of circulation. The system therefore responds to rotation with a Eulerian anti-Stokes flow that, in a time-averaged sense, is equal and opposite to the Stokes drift. Hasselmann [51] found that in a steady state, the Coriolis force introduces a stress perpendicular to the direction of wave propagation (i.e. along the wave crest), causing a Lagrangian particle to rotate in inertial circles and eventually to anti-align with the Stokes drift after one quarter day: the Eulerian anti-Stokes flow takes the form of inertial oscillations. Pollard [52] considered Gerstner waves in the presence of rotation and found that, although the dispersion relation of such waves is only slightly affected by rotation, their time-averaged Lagrangian velocity is zero, and inertial oscillations are generated in agreement with Ursell & Deacon [50] and Hasselmann [51]. Very recently, Constantin & Monismith [53] have reconsidered the solution by Pollard [52], but in the presence of a depth-invariant mean current, identification of an additional type of wave is made possible by the presence of the current that Constantin & Monismith [53] term an ‘inertial Gerstner wave’.

In fact, the so-called Hasselmann force [51] corresponds to the strong rotation limit of the vortex force identified in the seminal paper by Craik & Leibovich [54] that arises from interaction of the wave velocity field with the mean vorticity field. In the Craik–Leibovich (CL) momentum equation for the wave-averaged velocity ** u**, derived through perturbation methods [54] or using the generalized Lagrangian mean (GLM) theory of Andrews & McIntyre [55,56] (also including Coriolis forces), Stokes drift is responsible for three effects, the Coriolis–Stokes force, the Stokes correction to pressure and the vortex force:
2.7
where

**is the angular velocity of the Earth's rotation,**

*Ω**π*is pressure normalized by the fluid density and

*ν*is the viscosity. A similar equation, albeit for vorticity, was derived by Huang [57] using the methods of Craik & Leibovich [54] to study the development of the Ekman layer. The Coriolis–Stokes force, as previously considered by Ursell & Deacon [50] and Hasselmann [51], is a forcing of the mean flow in a rotating fluid: different from a Coriolis force in the absence of waves, the total Lagrangian velocity must be considered in this forcing term, and an additional term 2

**×**

*Ω*

*u*_{SD}arises as a result, which depends crucially on the vertical shear of the Stokes drift. The vortex force

*u*_{SD}×(

**∇**×

**) gives rise to the overturning force responsible for Langmuir circulation (see §6).**

*u*### (d) Stochastic sea states: particle dispersion

Although much less studied, the random velocity field due to small-amplitude surface gravity waves leads to the dispersion of particles in the same way as in a turbulent flow. This dispersion can be captured, for example, by computing the one-particle diffusivity of Taylor [58], which captures the temporal growth in the variance of the particle displacement with time or, more physically, the growth rate of a cloud of particles. In a random surface gravity wave field, the diffusion is not actually a consequence of the random linear wave motions themselves, but arises because of the interaction between the random Lagrangian wave-induced velocities. Dispersion is therefore a fourth-order effect in wave steepness and requires the inclusion of both second-order wave-induced components of the Lagrangian velocity: the Stokes drift and the Eulerian mean flow (cf. [59]).

By considering all the relevant wave–wave interactions in Lagrangian multi-chromatic second-order theory, Herterich & Hasselmann [60] computed the effective diffusivity of a random deep-water surface gravity wave field and found typical values of 10^{−2} m^{2} s^{−1} for a fully developed Pierson–Moskovitz spectrum [61] with a wind speed of 10 m s^{−1} (cf. molecular diffusivity is of the order of 10^{−9} m^{2} s^{−1}). Based on field measurements by Okubo [62], Herterich & Hasselmann [60] conclude that wave-induced diffusion will generally be negligible in the ocean, apart from on scales of 10–100 m. Buick *et al*. [63] perform experiments in a three-dimensional random sea basin and find good agreement with the theory of Herterich & Hasselmann [60]. Nevertheless, based on the diffusivity computed by Herterich & Hasselmann [60] and for realistic sea states, Pugliese Carratelli *et al*. [64] estimate wave-induced diffusion to be very important for small oil spills and to affect significantly some smaller-scale aspects (e.g. the formation of filaments) of larger-scale spills such as the one at Deepwater Horizon in the Gulf of Mexico in 2010. Spydell *et al*. [65] examined dispersion in the surf zone using drifters and found that the diffusivity from the Stokes drift of unbroken irrotational surface gravity waves (sea swell and infragravity waves) was much smaller than their observations, suggesting that rotational motions not directly associated with waves and Stokes drift are important to surf-zone dispersion. Although motivated by the possible explanation of mixing by internal waves, Bühler & Holmes-Cerfon [59] have recently examined the particle dispersion by random waves in rotating shallow water and found that rotation ‘chokes’ the Lagrangian transport and the diffusion due to random waves in shallow water.

## 3. Laboratory studies

Many authors have considered Stokes drift for periodic waves in laboratory wave flumes. To this day, there remains considerable confusion in the literature as to whether a net drift is or should be observed in laboratory experiments. It is self-evident that, in a closed experimental wave flume, the net depth-integrated mass flux must be zero, and the Stokes drift of a periodic wave train must be accompanied by a Eulerian return current driven by a set-up in the direction downstream of wave propagation, so that the steady-state depth-integrated Lagrangian drift is zero. Longuet-Higgins [66] was the first to point out that the mass-transport velocity in laboratory measurements in wave flumes can be significantly different from predictions based on the irrotational theory of Stokes [1]. He derived two classes of analytical solutions to explain how vorticity could be transported into the interior of the fluid: a ‘conduction’ solution and a ‘convection’ solution. Depending on the ratio of the wave amplitude *a* to the thickness of the boundary layer *δ*, the transport of vorticity takes place by viscous ‘conduction’ (*a*^{2}/*δ*^{2} small) from the bottom and free surface boundary layers, or convection with the mass-transport velocity (*a*^{2}/*δ*^{2} large), from the wavemaker or the beach at the other end, where vorticity can be generated. Furthermore, Longuet-Higgins [66] explored net transport in the thin oscillating boundary layer near the channel bottom, showing that the total Eulerian velocity just outside this layer is positive (in the same direction as the waves), that the total Lagrangian transport at this location is greater than the Stokes drift by a factor of 5/2 and that this enhanced transport does not disappear in the limit of infinite Reynolds number.

It is instructive to consider some of the concluding remarks by Longuet-Higgins [66] before reviewing more recent progress in detail. If one considers a motion that is started from rest, the motion in the interior of the fluid will always initially be irrotational, and it will take (considerable) time for vorticity to be advected or diffused from the vertical or horizontal boundaries, respectively. Finally, as also pointed out by Longuet-Higgins [66], the convection solution may not be stable and the instabilities may themselves be comparable to those responsible for Langmuir circulation [67]. A comprehensive early review of experimental observation of drift is given by Craik [67], distinguishing waves decaying temporally and spatially under the influence of viscosity and the effects of surface contamination.

Most practical experiments will be in the convection regime. For a wave frequency of 1.0 Hz and an amplitude of 3.0 cm (steepness *ak*=0.12), which can be generated in a typical laboratory flume, we obtain a Stokes drift velocity *u*_{SD}(*z*=0)=2.0 cm s^{−1}, so that vorticity takes 70 s to be advected over the distance of one wavelength, emphasizing the need for long-duration experiments for the convection solution to be established. Swan [68] demonstrates that convection indeed plays an important part within a relatively deep experimental wave flume: vorticity generated at the end conditions is convected backwards with the mass transport velocity and the near shore region progressively influences the entire length of the wave flume, although the flow field is not always stable. By installing a plastic sheet at the toe of the beach, Swan & Sleath [69] could obtain long-time stable conditions that agreed better with their fourth-order finite-depth extension the irrotational solution for Lagrangian transport in a closed domain. Umeyama [70] performs a similar expansion, but focuses explicitly on particle trajectories and finds reasonable agreement with experimentally obtained trajectories. Monismith *et al*. [48], who consider the experimental results by several authors, find that Stokes drift is not only cancelled by a Eulerian current in the depth-integrated sense, but that a cancellation takes place at all levels. Yet, many of the cases considered include a constant or sheared current even in the absence of waves. Based on their observations, Monismith *et al*. [48] hypothesize the existence of (rotational) Gerstner waves in wave tank experiments, a hypothesis further discussed in Weber [49]. Although the prospect of Gestner waves is theoretically appealing, resolution of this debate will more probably be found in careful consideration of the time-varying nature of any mean flow field, noting that even if stable, the convection solution of Longuet-Higgins [66] will take a very long time to establish, and the boundary conditions imposed at both ends.

In fact, Paprota *et al*. [71], who take their measurements after a relatively short wave train of periodic waves in a relatively long flume, find good agreement with the irrotational theory of Stokes [1], supplemented by a uniform return flow reflecting volume conservation in a closed domain. For waves of intermediate water depth (*kh*=*O*(1)) and very large steepness, Grue & Kolaas [72] find good agreement with nonlinear irrotational theory in the interior of the fluid in a set of very high-quality experiments. Their experiments are stopped long before the first waves reach the end of the flume, but the length of their wave train remains long relative to the water depth. Despite the good agreement with irrotational theory in the interior of the fluid, Grue & Kolaas [72] observe significant additional streaming and associated shear in both the bottom and free surface boundary layers, more than can be predicted by the boundary layer streaming solution of Longuet-Higgins [66], which the authors note may be invalid due to the large amplitude of the waves.

In the conduction regime, Groeneweg & Klopman [73] compare their more generally applicable GLM model for wave-current interaction to the conduction solution of Longuet-Higgins [66], showing near perfect agreement, and to the laboratory measurements in very shallow and long closed flumes by Mei *et al*. [74], finding good agreement for intermediate water depth (*kh*=1.0), but less good agreement for deep water (*kh*=1.8).

## 4. Field measurements

### (a) Buoys and drifters

Most information about waves on the sea surface is obtained from moored floating buoys that measure surface height fluctuations with internal accelerometers or global positioning system (GPS) receivers. Webb & Fox-Kemper [75] emphasize that it is often overlooked that the depth-dependent and depth-integrated Stokes drift are not easily measured from wave data, including from buoys, and that a large part of the uncertainty associated with these estimates stems from how the effect of directional spreading is handled. In recent years, making use of advances in compact and inexpensive sensor packages, small drifting buoys have been deployed that measure both surface waves and drift properties, and the data they generate compare well with conventional buoys [76]. In particular, the ‘Surface Wave Instrument Float with Tracking’ (SWIFT) [77] has been shown to be capable of measuring dissipation under breaking waves [77,78], and it is reasonable to believe that the buoy, which has very low windage (air resistance under the influence of wind), will indirectly measure the Stokes drift through its downward-looking acoustic Doppler current meter (ADCP).

Floating buoys and drifters do not collect measurements at a fixed point (Eulerian), but instead provide Lagrangian time series of the orbital motion of a water parcel at the surface, provided the buoy mooring is flexible and the density is close to that of water. It has been known since Srokosz & Longuet-Higgins [79] and Longuet-Higgins [80] that for deep water the high-frequency (HF) bound waves observed in a Eulerian reference frame that are responsible for the steepening of crests and the broadening of troughs are not present in Lagrangian records, and the mean water level in these records is subject to a set-up. Evidently, the wave period in a Lagrangian reference frame moving with the wave *T*_{L} also exceeds that in a Eulerian reference frame *T*_{E} as a result of the Stokes drift; the two are related by (*T*_{L}−*T*_{E})/*T*_{E}=*u*_{SD}/(*c*−*u*_{SD}) [80]. Recently, Herbers & Janssen [81] have extended the second-order Lagrangian wave theory of Srokosz & Longuet-Higgins [79] to finite depth and have examined second-order wave group signals in detail. Crucially, the set-down of the wave-averaged free surface, which forms in response to the divergence of the Stokes drift on the scale of a group (see §2b), can appear as a (significant) set-up in Lagrangian buoy records. Furthermore, in shallow water, where the magnitude of the set-down is typically strongly amplified, the distortion in its Lagrangian records is only expected to be small. Analysing velocity data from buoys, Herbers & Janssen [81] find clear evidence of the second-order Eulerian ‘frequency-difference’ or ‘infragravity waves’ that form in response to fluctuations in the Stokes drift in ‘groupy’ signals, emphasizing that such infragravity waves have periods small compared to the Earth's inertial period and thus do not cancel the Stokes drift exactly as suggested for periodic waves in the presence of rotation (see §2c).

### (b) Radar

HF radars have been used to measure surface currents since the early 1970s [82,83]. Surface waves scatter electromagnetic waves in the HF band (3–30 MHz) through Bragg diffraction. By analysing the Doppler spectrum of the back-scattered electromagnetic waves at grazing angle, an estimate of the near-surface current is found. This estimate rests on the assumption that the phase speed of the resonant surface wave is known through the dispersion relation. By comparing the location of the first-order Doppler peaks to their theoretical positions according to linear theory, an estimate of the underlying current component in the radial direction towards or away from the antenna is found. The radial current estimate *v*_{r} is equal to the shift in phase speed, *Δc*, and is considered a weighted average over depth [83], namely
4.1
where is the unit vector in the radar look direction, *k*_{B} is the Bragg wavenumber of the resonant surface wave, and ** u**(

*z*) is the quasi-Eulerian current defined as the difference between the Lagrangian current and the Stokes drift [84,85]. It can be argued that the vertically sheared Stokes drift should affect HF radar measurements through a similar mechanism, and Broche

*et al*. [86] pointed out that by substituting the Stokes drift associated with a finite amplitude Bragg wave

*u*

_{SD}=

*c*

_{B}(

*ka*)

^{2}e

^{2kBz}for the Eulerian current in (4.1) gives a phase velocity shift to the measured signal of 4.2 which is identical to the phase velocity shift expected from the second-order correction (arising at third order in the expansion) to the linear wave dispersion relation [1], 4.3 where is the unmodified dispersion relationship for deep water (correct up to second order).

Ardhuin *et al*. [87] argue that the current measured by HF radar is actually influenced by the Stokes drift associated with waves longer than the resonant wave and the Eulerian current itself (a filtered Stokes drift). This is in contrast with observations reported by Röhrs *et al*. [88], who compared near-surface drifters with Eulerian current measurements from an ADCP and HF radar currents. They found that the current measured with a typical HF radar system was in close agreement with ADCP currents, and that, after subtracting the surface Stokes drift velocity calculated by a wave model, the drifter current estimates agreed better with the HF currents.

For wave groups, the authors are only aware of observations by Smith [89], who used a long-range phased-array Doppler sonar to study the Eulerian and Lagrangian velocity at the surface of short wave groups. As such short wave groups pass, Eulerian counterflows occur that oppose the Stokes drift velocity at the surface (see §2b). Smith [89] found that the magnitude of this Eulerian return flow at the surface exceeds predictions based on an irrotational response, a finding that cannot be understood using existing theory.

## 5. Coastal zone: Stokes drift and the undertow

Many factors influence near shore circulation, including wind, wave breaking and Stokes drift. With waves directed towards a beach, Stokes drift leads to a build-up of fluid near the beach, which in turn generates a pressure gradient driving a bulk offshore flow, so that volume is conserved. The undertow is one of the most important mechanisms for sediment transport in near shore regions [90]. Svendsen [91] provided one of the first rigorous explanations of the undertow as driven by the local difference between radiation stress and the set-up pressure gradient, which develops in response to wave breaking. Turbulent shear stresses evidently play a key role. More recently, Guannel & Özkan Haller [92] have shown that nearly all existing theoretical formulations of undertow can be reconciled under the confines of linear wave theory, provided the bottom shear stress is incorporated correctly.

Moving further offshore, Lentz *et al*. [93] have presented measurements of the wave-driven offshore flow (undertow) seawards of the surf-zone and find that a typical Eulerian velocity profile has a maximum near the surface and decreases towards the bottom: it is equal and opposite to the Stokes drift. This is in contrast with the surf-zone, where the offshore Eulerian velocity profile is typically parabolic with its maximum at mid-depth. Extending theoretical work by Xu & Bowen [94], Lentz *et al*. [93] showed that their observations are consistent with an inviscid balance between the Coriolis force associated with the offshore flow and the Hasselmann or Coriolis–Stokes force associated with the Stokes drift. Instead of in the form of a steady undertow, mass balance can also be obtained in the form of unsteady transient rip currents (e.g. [95]). We refer the reader to Lentz & Fewings [96] for a recent review of inner-shelf circulation.

## 6. The role of the Stokes drift in Eulerian ocean models

Traditionally, ocean waves and the interior ocean have been modelled independently, each forced by atmospheric fluxes without feedback to the atmospheric boundary layer or an exchange of fluxes between waves and the ocean surface boundary layer (OSBL). This is inconsistent because ocean models do not resolve surface waves. The equations of motion, the turbulence kinetic energy (TKE) budget and the other conservation equations should all be modified to account for the presence of surface waves. In principle, this entails calculating the full two-dimensional wave spectrum, because the Stokes drift at a given vertical level requires the integration of (2.4) from Kenyon [26]. This is too computationally demanding and impractical for most model systems, although Webb & Fox-Kemper [75] did employ the full two-dimensional wave spectrum by coupling the WaveWatch-III model [97] in the Community Earth System Model [98]. Instead, simplified wave fields (monochromatic waves) have often been used [99–104].

Representation by a single monochromatic wave is problematic because the vertical rate of shear of the Stokes drift in a broad-banded spectrum is stronger than that of a monochromatic wave. Furthermore, the deep Stokes drift profile will be stronger than that of a monochromatic wave because the longer waves, which are missing in a monochromatic representation, penetrate much deeper [105]. Although approximate profiles exist that better fit the Stokes drift profile under a broad-banded wind-generated wave spectrum [106], it is clear that in order to model adequately the Stokes profile in mixed wind-driven sea and swell conditions, the full two-dimensional wave spectrum must be represented by a spectral wave model (see [107–112] for examples and descriptions of third-generation spectral wave models). This can only be achieved by coupling the wave model to the ocean model.

The depth of the OSBL is maintained by a balance between different turbulent processes, such as buoyancy production through heating and cooling and shear production. It is widely found that mixing in ocean models exhibits the greatest deficiencies in the wave-rich extra-tropics [113,114], and to first order an enhancement factor with a longitudinal dependency can make a significant impact, as was reported by Breivik *et al*. [115]. This is borne out by the observation by D’Asaro *et al.* [116] that inclusion of CL wave forcing seems to have little effect at tropical or subtropical latitudes, but appears to increase boundary layer depths at high latitudes by 15–20% on average. This is also consistent with the observation by Belcher *et al*. [117] that the effect of waves in enhancing boundary layer turbulence varies with latitude, and helps explain model biases in mixed layer depth in the extra-tropics. This has consequences beyond the modelling of the upper ocean because sea surface temperature biases can affect the deep convection of the atmospheric circulation [118].

Two ocean models in particular have been used in integrated wave–ocean model experiments, namely the Regional Ocean Model System (ROMS; see [119]) and the Nucleus of European Modelling of the Ocean (NEMO; see [120]). ROMS has been extended to incorporate the vortex-force formalism by Uchiyama *et al*. [121] and, through a coupled atmosphere–wave–ocean (WRF-SWAN-ROMS) set-up, has been used extensively for near shore applications, where wave effects play an important role [122–124]. NEMO has been employed mostly on larger scales and is the ocean model component in the coupled forecast model of the European Centre for Medium-Range Weather Forecasts (ECMWF; see [115,125]). Recently, however, NEMO has also been tested on regional domains of the North Sea and the Baltic Sea on much higher resolution (3.5 km) with wave fields from a wave model on the same resolution [126,127].

As the Stokes drift enters the wave-averaged momentum equation, as well as the tracer advection equation and the TKE equation, it is clear that estimates must be made of its importance for the circulation and hydrography in Eulerian ocean models. McWilliams & Restrepo [128] used a wind climatology to assess the impact on the general circulation from adding the Coriolis–Stokes and vortex forces, as well as Stokes drift to the tracer advection equation. These authors found that, in the extra-tropics, the impact can indeed be quite significant, with wave effects amounting to up to 40% of the wind-driven Ekman transport. Below, we will review briefly how the Stokes drift is thought to affect the wave-averaged Navier–Stokes equations appropriate for ocean models, putting practical applications to (2.7) (see also [100,128]).

### (a) Stokes drift in the momentum and tracer advection equations

First, the Coriolis–Stokes 2** Ω**×

*u*_{SD}force is of significance for ocean modelling outside the tropics, where the wave field is mostly dominated by swell [128]. A number of workers have recently introduced the term in models of varying complexity, from a simple vertical column, horizontally homogeneous idealizations [101,129,130] to fully resolved three-dimensional models [115,121].

Second, the CL vortex force, *u*_{SD}×(**∇**×** u**), like the Coriolis–Stokes force, arises as an explicit term in the wave-averaged Navier–Stokes equations (2.7) (see [54,131]). The term is of greatest importance when modelling near shore, shallow-water conditions because of the strong shear in Eulerian currents found there. Several studies have incorporated the vortex force in coastal applications, notably Uchiyama

*et al*. [121] and Warner

*et al*. [123]. The role of the CL vortex force in Langmuir turbulence is however important throughout the world's oceans, as will be discussed below.

Third, the tracer advection equation must also be modified in order for Eulerian ocean models to account properly for unresolved wave effects. Following McWilliams & Sullivan [100], the conservation equation for a scalar quantity *c* becomes
6.1
where SGS stands for sub-grid scale effects.

### (b) Langmuir turbulence

Langmuir cells were first observed by Langmuir [132], who hypothesized that the narrow convergent bands of debris often observed on the sea surface were related to wave activity. Craik & Leibovich [54] offered the first model explaining how waves could set up instabilities in a homogeneous ocean. Two instability mechanisms were put forward by which the Stokes drift could destabilize the water column and form counter-rotating gyres in accordance with observations. The first was premised on the periodicity of the Stokes drift in the cross-wind direction, which would lead to an instability now known as CL1, after Faller & Caponi [133]. This ‘direct drive’ instability arises because the vortex force, if *u*_{SD} varied in the cross-wind direction, would not be balanced, thus rotation would ensue. The second mechanism envisioned by Craik & Leibovich [54] is caused by the vertical shear of the Stokes drift which tilts the vertical vorticity due to a, possibly small, disturbance of the mean surface current into the horizontal. This has become known as the CL2 instability mechanism. Here, no coherent wave structure is required, only a horizontal shear in the Eulerian current. Only CL2 is believed to be of importance in the open ocean, because there is nothing to suggest that the Stokes drift under a broad wind–sea spectrum would exhibit the required periodicity in the cross-wind direction. Recently, Suzuki & Fox-Kemper [134] revisited the canonical CL mechanisms and demonstrated how Stokes advection, the Coriolis–Stokes force and the Stokes shear force mediate energy transfers between the mean flow and the wave field.

Since the crucial large eddy simulations (LES) by Skyllingstad & Denbo [103] and McWilliams *et al*. [131], which demonstrated how Langmuir cells in the open ocean can create turbulence in the OSBL, termed Langmuir turbulence, several authors have investigated how Langmuir mixing can be parametrized in Eulerian ocean models. This is usually achieved by comparing the parametrizations for idealized cases with LES [135–137]. The turbulence models fall in two broad categories: turbulent profile parametrization and explicit modelling of the TKE equation. These are discussed separately below.

We will now briefly follow Ardhuin & Jenkins [138] and investigate the effect of non-breaking waves on turbulence. The shear-induced production of TKE, per unit volume, resulting from the organized wave and mean current motions, is, using indicial notation
6.2
where fluctuation quantities are denoted by primes. By assuming that the GLM quantities can approximate the fluctuation quantities in (6.2), Ardhuin & Jenkins [138] followed the derivation by Andrews & McIntyre [55] of the relation between a Eulerian average and the corresponding GLM quantity (their eqn (2.28)):
6.3
Here, *ϕ* can represent any oceanic quantity like temperature, salinity or a vector component of the flow field (as will be assumed in the following), the superscript L denotes GLM quantities, the over-bar Eulerian wave averaging and the quantities *ξ*_{i} the fluctuation-related displacements. By realizing that the wave-induced Stokes drift is a GLM quantity, the wave-induced contribution to the production term (6.2) may be expressed as
6.4
It is now clear (as was also assumed by [131,139]) that it is the *shear* of the Stokes drift that gives rise to Langmuir production. It should also be noted that (6.4) actually represents the primary mechanism by which all non-breaking wave–mean flow interactions act to generate turbulence, not just Langmuir circulation.

#### (i) *K*-profile parametrizations of Langmuir production

The *K*-profile parametrization by Large *et al*. [140] represents the most straightforward turbulence closure through which to introduce the impact of Langmuir turbulence. The first of these parametrizations was presented by McWilliams & Sullivan [100] in which the turbulent Langmuir number
6.5
is used to boost the turbulent velocity scale
6.6
Here, *ϕ* is the Monin–Obukhov stability function, *κ* is von Kármán's constant and *u*_{*} is the water-side friction velocity. McWilliams & Sullivan [100] proposed an enhancement factor of the form
6.7
Here, *α* and *L*_{w} are constants. This parametrization was followed shortly after by a more complicated formula by Smyth *et al*. [141], which allowed stratification in the upper ocean to limit the intensity of Langmuir turbulence. Similar parametrizations, also based on the turbulent Langmuir number were later presented by Harcourt & D’Asaro [136] and Takaya *et al*. [142].

In one of the first studies involving a fully coupled atmosphere–wave–ocean model, Fan & Griffies [114] compared the impact of these two parametrizations and found the latter to yield closer agreement with the wintertime mixed-layer depth. A problem with such parametrizations is that although the Stokes drift should have a bearing on the mixing, the direction of the Stokes drift must also be taken into account, because crossing seas and a directionally wide spectrum will reduce the Stokes drift [75,87,105]. Van Roekel *et al.* [143] argued that the Stokes drift should be projected into the direction of the Langmuir cells before the turbulent Langmuir number is estimated, as misalignment affects the intensity of Langmuir turbulence. McWilliams *et al*. [144] also showed that the presence of swell at large angles to the local wind sea would modify the Langmuir production. This raises the question of how Langmuir turbulence can be parametrized, since swell is non-local and thus not related to the local wind. This remains an open question, but Li *et al*. [145] showed that a climatology of the enhancement factor calculated from a global wave model integration has comparable performance to a run incorporating the wave model itself, suggesting that the impact of swell is small enough to be ignored, at least to first order. Li *et al*. [146] showed that parametrizing the Stokes drift from the wind alone and estimating a layer-averaged broad-band Stokes profile [105,106] did in fact also yield results comparable to those from the full wave model.

#### (ii) Langmuir production in second-moment turbulence closure models

Second-moment closure models in the vein of Mellor & Yamada [147] handle the TKE directly by parametrizing the terms in the TKE equation [148,149]. As Kantha & Clayson [150] note, this makes physical sense for the TKE equation because it can be derived from first principles, but the corresponding equation for the turbulent length scale, or equivalently for the dissipation rate or turbulent time scale, cannot. Although the results reported are often in reasonable agreement with large-scale properties such as the mixed-layer depth and surface temperature, it is nevertheless unsatisfactory to mould the evolution equation for the turbulent length scale on the TKE without any guiding physical principle. Nevertheless, several attempts have been made at extending the standard second-order moment closure for ocean models [147,148] to account also for Langmuir turbulence. Kantha & Clayson [150], Carniel *et al*. [99] and Janssen [16] added a production term proportional to the shear of the Stokes drift velocity ∂*u*_{SD}/∂*z* (cf. (6.4)), in vertical one-equation turbulence models with an algebraic expression for the relation between the energy dissipation *ε* and the mixing length, *ε*≃*q*^{3}/*l*, where *q* is the root-mean square turbulent velocity and *l* is the mixing length.

For wave-averaged quantities and small-amplitude (non-breaking) waves, which are irrotational to first order, the turbulent kinetic energy equation can be written as
6.8
where is the TKE per unit mass (with *q* the turbulent velocity) and *ε* the dissipation rate. Further, *w*′*e*′ and *w*′*p*′ are the turbulent transport and pressure correlation terms, respectively [151]. Here, the gradient hypothesis is adopted (e.g. [152], p. 204), and the shear-production and the buoyancy-production terms are proportional to the squared shear and the buoyancy frequency *N*^{2}=−(*g*/*ρ*_{w})d*ρ*_{w}/d*z*, respectively (*ν*_{h,m} are turbulent diffusion coefficients). The Langmuir production term corresponds to in (6.4), here written per unit mass. This quantity obviously depends heavily on the form of the Stokes drift velocity profile. Only a modest impact on the overall mixing was found in these studies, and Kantha & Clayson [150] cautioned that the results were within the error bars of the observed dissipation rates found observationally. Harcourt [153,154] went further and introduced Langmuir turbulence into a full second-moment closure model for the oceanic surface boundary layer and found that the effect was to deepen the mixed layer.

## 7. Tracer transport: oil, micro-plastics and search and rescue

Stokes drift is important for the drift of objects and matter at or near the sea surface. It has long been standard practice to incorporate an estimate of the surface Stokes drift when calculating drifter trajectories [155] and the dispersion of oil [100]. Moderately accurate empirical relations exist between the wind speed and the surface Stokes drift, which admit estimation of the Stokes drift from the local wind, assuming an average wind fetch and duration [26,87]. Yet, the impact of swell cannot be assessed without resorting to a wave model. Even if the surface Stokes drift can be estimated reasonably well from the local wind, the rapid decrease with depth of the Stokes drift complicates the picture for submerged objects and suspended matter. The strong vertical shear of the Stokes drift has been found to have a profound impact on the fate of oil. Drivdal *et al*. [156], using the General Ocean Turbulence Model (GOTM; see [157]) with modifications to account for an impulse formulation of the Coriolis–Stokes effect and enhanced mixing due to breaking waves, investigated how the increased mixing by wave breaking and Stokes shear production, as well as the stronger veering by the Coriolis–Stokes force would affect the drift of oil particles from an underwater release of oil spilled from a North Sea platform. They found the net drift to be both slower and more deflected away from the wind direction than predicted by oil drift models that did not include wave effects. Jones *et al*. [158] similarly found that the continuous mixing of oil particles by breaking waves can lead to large discrepancies in the trajectories of oil spills, as the subsurface waters experience a weaker Stokes drift.

The situation is more complex still for objects that partially protrude from the water. These will be subject to forces by the wind, currents and Stokes drift. In most cases, their drift properties must be quantified empirically [159–161], and at best one can hope to find a relation between the wind speed and direction by measuring these simultaneously with their leeway, that is the motion of the object relative to the ambient current and the wind speed and direction [161]. In such cases, quantification of the impact of the Stokes drift velocity on the overall leeway is difficult, and it is commonly assumed that the relation between the wind speed and the leeway also accounts for the Stokes drift. This implicitly assumes that the Stokes drift is in the direction of the local wind [162]. In light of the large uncertainties associated with trajectories of drifting objects, this may in many cases be a reasonable assumption, but with the notable exception of swell-dominated regions. Trinanes *et al*. [163] report that in the search for Malaysian Airline MH370, the low-windage debris found in the western Indian Ocean had drifted for months in an area with strong swell at large angles to relatively weak winds. In this case, the explicit inclusion of the Stokes drift altered the trajectory of the debris field dramatically.

The relatively new activity of mapping the pathways of plastic pollution including micro-plastic in the world's oceans [164,165] and its impact on biology [166] is a field in which, like oil fate modelling, the effect of waves is felt through both Stokes drift and enhanced mixing. The drift pattern of near-surface suspended matter will be determined by the total current *u*_{L}=*u*_{E}+*u*_{SD}, while wave-induced mixing determines the vertical distribution of particles [167]. It is not self-evident that tracers of all sizes and densities behave as ideal Lagrangian tracers [168,169].

The Stokes drift will also play a role in the nascent field of Lagrangian coherent structures (LCSs), which deals with the trapping and expulsion of passive tracers and drifters by material surfaces advected by the flow field [170]. LCSs have been shown to be able to trap and concentrate particles in large-scale flows and, likewise, to repel particles from certain areas. Such behaviour can be greatly altered by the presence of the Stokes drift for two reasons. First, as Stokes drift is determined by the wave field, which in turn is driven by the wind, the Stokes drift varies on much larger spatial scales than the Eulerian current field, which is dictated by baroclinic processes as well as wind forcing. Second, as the Stokes drift follows the direction of the wave field, it often advects fluid onshore, something the Eulerian current cannot. Allshouse *et al*. [171] show that in the case of an oil field off the western coast of Australia, the wind-induced drift could potentially take an oil spill all the way to the vulnerable Ningaloo coral reef, whereas studies of Eulerian currents alone suggest that LCSs would capture the oil in large gyres far offshore.

## 8. Conclusion and future directions

Surface Stokes drift can been measured directly through the drift of passive objects with minimum windage. Using photographic techniques, it is possible to observe the full vertical extent of the Stokes drift directly in laboratory settings (figure 1), but such measurements in laboratory wave flumes are associated with difficulties of their own, which must be resolved by systematic study. Future experimental studies must be conducted keeping in mind clearly the role of the different experimental parameters: (i) the depth of the experimental flume, relative to both the wavelength and the group length; (ii) the duration of the experiment, relative to both the individual wave period and the time scale associated with the advection or diffusion of vorticity from the boundaries; (iii) the steepness of the waves, noting that nonlinear contributions beyond second order only play a role in finite depth; (iv) the ratio of wave amplitude to boundary layer thickness to examine the effect of boundary layer streaming beyond second order; and (v) the role of (sheared) currents that are present even in the absence of waves.

Owing to its Lagrangian character, the Stokes drift remains an elusive quantity which is not measurable with Eulerian field instruments. Worldwide Stokes drift climatologies so far have therefore been calculated from wave models [104,172]. Space-borne remote sensing appears, however, to hold some promise for mapping Stokes drift on a global scale. Chapron *et al*. [173] reported the first direct current measurements using the Doppler centroid technique with synthetic aperture radar ASAR. Chapron *et al*. [173] considered their radial current estimate *U*_{D} to correspond to a mean motion of the radar-scattering sea surface elements. As such, the estimate should also contain a mixture of the Stokes drift and the quasi-Eulerian current [174]. Presently, an effort is being made to launch a dedicated space-borne instrument for ocean current monitoring, the ‘Sea surface KInematics Multiscale’ instrument (SKIM) [175]. This instrument, if constructed, will consist of a pencil-beam rotating Ka-band altimeter, which goes from the nadir (regular altimeter) to a 12° incidence angle. Although the instrument will measure the quasi-Eulerian current contaminated by a wave bias [173] related to the Stokes drift, the ability to measure also the surface wave spectrum (down to an estimated 20 m wavelength) will allow for a direct estimate of the (radial) Stokes drift.

Separately, it is worth noting that the number of wave buoys in the ocean, from which we derive most of our knowledge about waves, is extremely limited: Ardhuin *et al*. [175] estimate one buoy per 1000 km of European shoreline. Cheaper sensor and GPS transmitter technology may enable future observation on an unprecedented scale, including by large numbers of small free drifting buoys, and we note new commercial technologies such as the Spoondrift Spotter and the Autonaut Unmanned Surface Vessel, among others.

The past decades have shown that waves play a significant role in the upper ocean layer through enhanced mixing and significant alteration of the momentum budget of the ocean. The only way to model consistently the coupled atmosphere–ocean system is, ultimately, to incorporate a wave model and thus provide fluxes and fields to the interior ocean. As the resolution of ocean models and indeed coupled systems continues to increase and coastal areas become better resolved, this will become increasingly important. It is also clear that wave effects can now be modelled in areas where such were not even considered two decades ago. The interaction between sea ice in the marginal ice zone and the oceanic wave field is one such topic, which is actively investigated [176]. Ardhuin *et al*. [177] found that icebergs are associated with anomalies in the climatology of wave heights in the Southern Ocean. Coupled climate and forecast systems contain thermodynamic ice models that interact with the ocean and atmosphere components. As wave models become incorporated into these coupled model systems, it seems likely that the Stokes drift will be found to play a significant role in the distribution and climatology of sea ice as well as being an essential mixing process in the the wider climate system.

## Data accessibility

This article has no additional data.

## Authors' contributions

T.S.v.d.B. coordinated the writing of §§1, 2, 3, 4a and 5; Ø.B. that of §§4b, 6 and 7. Both authors read, edited and approved the manuscript.

## Competing interests

The authors declare that they have no competing interests.

## Funding

Ø.B. acknowledges support from CMEMS COPERNICUS through the WAVE2NEMO contract as well as from the Research Council of Norway through the projects RETROSPECT (grant no. 244262) and CIRFA (grant no. 237906).

## Footnotes

One contribution of 19 to a theme issue ‘Nonlinear water waves’.

- Accepted August 31, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.