## Abstract

Nonlinear effects in Lagrangian sea surface motions are important to understanding variability in wave-induced mass transport, wave-driven diffusion processes, and the interpretation of measurements obtained with moored or free-drifting buoys. This study evaluates the Lagrangian vertical and horizontal motions of a particle at the surface in a natural, random sea state using second-order, finite-depth wave theory. In deep water, the predicted low-frequency (infragravity) surface height fluctuations are much larger than Eulerian bound wave motions and of the opposite sign. Comparison to surface elevation bispectra observed with a moored buoy in steady, high-wind conditions shows good agreement and confirms that—in contrast to the Eulerian sea surface motion with predominant phase coupling between the spectral peak and double-frequency harmonic components—nonlinearity in Lagrangian wave observations is dominated by phase-coupled infragravity motions. Sea surface skewness estimates obtained from moored buoys in deep and shallow sites, over a wide range of wind–sea and swell conditions, are in good agreement with second-order theory predictions. Theory and field data analysis of surface drift motions in deep water reveal energetic [*O*(10) cm s^{−1}] infragravity velocity fluctuations that are several orders of magnitude larger and 180° out of phase with Eulerian infragravity motions. These large fluctuations in Stokes drift may be important in upper-ocean diffusion processes.

## 1. Introduction

The nonlinearity of ocean surface waves affects the geometrical properties of the sea surface and is important for understanding wave-induced transport and drift characteristics. Second-order nonlinear effects include the familiar enhanced steepness of wave crests (Stokes 1847) and mean water level variations on the scale of wave groups (Longuet-Higgins and Stewart 1962). The associated deviations from Gaussian sea surface statistics and variations in the wave-induced surface drift (commonly known as “Stokes drift”) are important in the interpretation of remote sensing data, in particular the precise measurement of sea level with satellite altimeters (e.g., Srokosz 1986; Rodriguez 1988) and radar observations of surface currents (e.g., Longuet-Higgins 1986). Whereas the weakly nonlinear theory for a two-dimensional, random sea surface is well established (e.g., Phillips 1960; Hasselmann 1962), it is not well understood how nonlinearity is manifested in Lagrangian measurement records, such as obtained by moored and free-drifting surface-following instruments. Moreover, accurate field observations are scarce owing to the difficulty of obtaining nonintrusive in situ measurements of wave motion at the sea surface and the cost and limited availability of high-resolution airborne topographic mappers.

The most widely available wave-resolved sea surface observations are from moored surface-following buoys that measure surface height fluctuations with an internal sensor package equipped with accelerometers or a global positioning system (GPS) receiver. Recent advances in compact and inexpensive sensor packages have enabled the development of small drifting buoys that measure both the surface wave and drift properties (Herbers et al. 2012; Thomson 2012; Pearman et al. 2014). Whereas the accuracy of the buoy sensors is reasonably well established, the interpretation of measurements is complicated by the fact that surface-following buoys do not collect measurements at a fixed location but instead provide Lagrangian time series of the orbital motion of a water parcel at the surface. Srokosz and Longuet-Higgins (1986) and Longuet-Higgins (1986) present a second-order theory of Lagrangian buoy motion in deep water and show that the high-frequency bound waves observed in an Eulerian reference frame are replaced by a change in mean sea level in the Lagrangian surface record. Interestingly, this change in sea surface properties does not affect the sea surface variance and skewness (Srokosz and Longuet-Higgins 1986).

In the present work, we revisit some of the results by Srokosz and Longuet-Higgins (1986), compare field observations to theoretical predictions, and discuss the implied low-frequency (infragravity) modulations of surface elevation and Stokes drift that are important to understanding, for example, satellite altimetry, surface dispersion of pollutants, and the interpretation of infragravity wave signals in buoy records. We extend the second-order theory of Srokosz and Longuet-Higgins (1986) to finite water depth (section 2), explicitly consider the Lagrangian infragravity motion, and compare theoretical predictions to field observations from moored and drifting buoys. The theory and data analysis show that the nonlinearity of wave orbital motion manifests itself in infragravity fluctuations of surface elevation (section 3) and Stokes drift (section 4) that are orders of magnitude larger than their Eulerian counterparts. The results are summarized in section 5.

## 2. Lagrangian sea surface height variations

Surface-following wave buoys provide Lagrangian measurements of the wave orbital motion at the sea surface. In the linear approximation, the measured vertical buoy displacement record is equivalent to an Eulerian measurement of surface elevation at a fixed location, but the horizontal wave orbital excursions introduce a distortion at second order in wave steepness (see Srokosz and Longuet-Higgins 1986; Longuet-Higgins 1986). Notably, in the Lagrangian frame of reference, second-order high-frequency bound waves are exactly cancelled out so that the characteristic steepening of wave crests and broadening of troughs in deep-water Stokes waves is not observed in a buoy record.

In addition to wave nonlinearity, buoy measurements are also affected by mooring response (for a moored buoy) or surface currents (for a drifting buoy). The mooring response is difficult to quantify and is not considered here under the assumption that it affects buoy motions primarily at time scales longer than the periods of the dominant waves and associated (infragravity) group modulations.

To describe the motions recorded by a small surface-following buoy, we consider a surface particle (at *z* = *η*) that follows the Lagrangian wave orbital motion while being advected with a surface current **U**. This surface current may include the Stokes drift as well as ambient tidal and wind-driven contributions for a drifting buoy and can be set equal to zero for a moored buoy. For simplicity, we assume here that variations in **U** on the space and time scales of the dominant waves are small so that, in the local wave field description, **U** can be approximately considered steady and uniform in space. Furthermore, we assume that the ratio between the mean surface current |**U**| and a characteristic wave phase speed *c*_{p} is of the same order, or smaller than, the wave steepness so that the effects of the current on second-order nonlinear wave kinematics may be neglected. The assumption that |**U**| ≪ *c*_{p} is reasonable in most open-ocean environments where surface currents are generally weak and typically one or two orders of magnitude smaller than ocean swell phase speeds. In the analysis presented here, we thus neglect all effects of wave–current interactions on local wave kinematics, with the exception of the familiar leading-order Doppler shift effect on the wave frequency as may be observed in drifting buoy records.

To evaluate the horizontal [**x**_{b}(*t*)] and vertical [*z*_{b}(*t*)] buoy positions in a stationary and spatially homogeneous sea state, we use a fully two-dimensional spectral description of the Eulerian sea surface, including second-order bound waves (e.g., Phillips 1960; Hasselmann 1962):

where the wavenumber vector is defined as **k** = *s*_{±}(*k* cos*θ*, *k* sin*θ*) with the sign index *s*_{±} defined as +1 for positive *ω* and −1 for negative *ω*. The wavenumber magnitude *k* ≡ |**k**| obeys the linear gravity wave dispersion relation *ω*^{2} = *gk* tanh(*kh*), where *g* is gravity and *h* is the water depth. The complex amplitudes obey the symmetry relation , where * denotes the complex conjugate. The quadratic terms in Eq. (1) are bound wave components with the sum frequency *ω*_{3} = *ω*_{1} + *ω*_{2} and wavenumber **k**_{3} = **k**_{1} + **k**_{2} of a pair of primary wave components. The coupling coefficient

depends on the primary wave frequencies, spreading angle and water depth.

To evaluate the vertical displacement time series

we take **x**_{b} to be the sum of the initial position **x**_{0}, a steady drift with velocity **U**, and a fluctuation approximated with the leading-order local wave orbital displacement

with the phase function

In the weak current (|**U**| ≪ *c*_{p}) approximation employed here, the only effect of the surface drift on the observed wave record is a small Doppler shift −**k** ⋅ **U** in the wave frequency [Eq. (5)]. Substitution of Eqs. (3)–(5) into Eq. (1) and expanding for small wave orbital displacements yield to second order in wave steepness:

with a coupling coefficient ,

that is the sum of the Eulerian coefficient and an additional Lagrangian term

The Lagrangian contribution to the coupling coefficient accounts for the horizontal buoy displacements by the wave orbital motion. In deep water, Eq. (8) is in agreement with the result derived by Srokosz and Longuet-Higgins [1986, their Eqs. (6.10) and (6.3)].

Since the Lagrangian corrections obey the antisymmetry relation

they do not affect the second and third cumulants of the sea surface height time series (Srokosz and Longuet-Higgins 1986). However, the diagonal (*ω*_{1} + *ω*_{2} = 0) contributions to Eq. (8) result in a change in mean water level (affecting the first cumulant), and the distortion of spectral properties is important in the detailed analysis of surface-following buoy measurements. For example, consider the interactions of a pair of wave trains in deep water with frequencies *ω*_{1,2} = *ω* ± Δ/2 (where 0 < Δ < *ω*), traveling in the same direction *θ*_{1,2} = 0. The sum–frequency interaction yields perfect cancellation between the Eulerian and Lagrangian terms:

causing the well-known absence of the double-frequency harmonic components in a Lagrangian wave record. Interestingly, this cancellation does not occur for the difference interaction

For small values of the difference frequency Δ, the magnitude of the negative term is a factor Δ/*ω* smaller than , and thus the familiar second-order setdown effect under wave groups is replaced by a much larger apparent setup signal in a Lagrangian buoy record.

These effects are illustrated in Fig. 1 with a numerical simulation of an energetic narrowband wave field in deep water. The Eulerian [Eq. (1)] and Lagrangian [Eq. (6)] sea surface height variations were evaluated at an arbitrary location **x**_{0} in the absence of surface drift (**U** = 0) for a random Gaussian sea state with a significant wave height of 8 m and spectral peak frequency of 0.1 Hz. To simulate a narrow swell beam, a two-dimensional Gaussian-shaped spectral energy distribution was used with standard deviations of 0.007 Hz (in frequency) and 5° (in direction). As expected, the results show the double-frequency harmonic components disappear in the Lagrangian reference frame (Fig. 1, middle panel) and the occurrence of an infragravity modulation of comparable magnitude, with maximum setup in the center of the wave groups (Fig. 1, bottom panel).

These infragravity surface height modulations are closely related to the mean Lagrangian water level change discussed in Srokosz and Longuet-Higgins (1986). That is, for a single wave train in deep water with frequency *ω* and amplitude *a* = 2|*A*_{ω}|, the difference interaction in Eq. (6) yields a mean water level change of *ω*^{2}*a*^{2}/2*g*, consistent with Srokosz and Longuet-Higgins (1986). In a bichromatic wave field, consisting of two wave trains with slightly different frequencies and the same amplitudes, this water level change is split equally between a mean setup [the self–self-difference interaction terms in Eq. (6)] and an infragravity group modulation [the cross-difference interaction terms in Eq. (6)]. In a random wave field with an infinite number of frequency components, there are no longer distinct mean and oscillating contributions to the sea level change but instead a continuous spectrum of low-frequency variations.

The situation is rather different in shallow water where the bound wave forcing approaches resonance [], causing strong amplification of the Eulerian coupling coefficient [Eq. (2)] for both sum and difference interactions. In contrast, the corresponding Lagrangian correction terms are not as strongly amplified [Eq. (8)], and thus buoy measurements of nonlinear sea surface properties in shallow water are not expected to be significantly distorted by the Lagrangian displacements.

## 3. Third-order statistics of sea surface height

To examine how nonlinearity affects Lagrangian surface height variations in natural wind-generated ocean waves and verify theoretical predictions of these effects, we analyzed third-order statistics of moored buoy observations. First, we present a bispectral analysis of a nearly fully developed wind sea in strong winds with the objective to characterize and verify the Lagrangian bound wave contributions in the spectral domain. Next, in order to quantify the nonlinearity over a wider range of conditions, we evaluate the sea surface skewness from long-term buoy records in deep and shallow water.

### a. Bispectral analysis

Although the nonlinear distortion of surface wave profiles is generally subtle (e.g., Fig. 1), the second-order bound waves cause deviations from Gaussian statistics that are important for, for example, the sea state bias in satellite altimetry (e.g., Rodriguez 1988) and nearshore sediment transport (e.g., Bailard and Inman 1981; Hoefel and Elgar 2003). Bispectral analysis (Hasselmann et al. 1963) is the natural tool to explore the non-Gaussian properties of a natural random wave field and identify nonlinear coupling between wave components across the frequency spectrum. Here, we evaluate the bispectrum of an idealized, small, surface-following buoy. Since we will apply the analysis to moored buoys, we set the mean drift velocity **U** in Eq. (5) equal to zero. Choosing a coordinate system centered on the mean buoy position so that **x**_{0} = 0, the vertical buoy elevation *z*_{b}(*t*) [Eq. (6)] can be expressed as a simple Fourier sum:

with a transform

In the limit of small frequency bandwidth Δ*ω*, a continuous bispectrum *B*(*ω*_{1}, *ω*_{2}) can be defined as

where denotes an averaging operator. Substitution of Eq. (13) in Eq. (14), assuming the primary wave amplitudes *A*_{ω,θ} are statistically independent and Gaussian, yields to lowest order:

with *E*(*ω*, *θ*) the (double sided in frequency) frequency–directional spectrum of primary waves

Equation (15) is the theoretical expression for the Lagrangian surface height bispectrum. The Eulerian surface height bispectrum is given by the same equation with replaced with [see Hasselmann et al. (1963), for a more general and formal derivation of the bispectrum].

To explore the Lagrangian sea surface statistics in a natural ocean wave field and verify the theoretical bispectrum [Eq. (15)], we use data from a moored Datawell DWR-G7 Directional Waverider buoy that was deployed in 157-m depth off the California coast near Bodega Bay during June 2010, as part of the Office of Naval Research High-Resolution Air–Sea Interaction (HIRES) research initiative (Herbers et al. 2012). A wind sea near full development with a significant wave height of about 4 m was observed over several days in persistent strong (13–15 m s^{−1}) winds. The steady conditions of this event allow for a detailed bispectral analysis to quantify nonlinear effects in the spectral domain. A 40-h-long record (1400 UTC 14 June through 0600 UTC 16 June) at the peak of this event was selected for analysis. Both spectra and bispectra were computed from 26.7-min-long data segments and smoothed through ensemble averaging and the merging of 13 frequency bands to a resolution of 0.0081 Hz. The frequency–directional wave spectrum was estimated from cross spectra of the three-component buoy displacement data using the maximum entropy method (MEM) of Lygre and Krogstad (1986). The frequency and frequency–directional spectra (Fig. 2; for convenience transformed to single-sided *f* spectra with the cyclic frequency defined as *f* = *ω*/2*π*) show the familiar properties of a unimodal wind-wave spectrum that is narrow at the peak and broadens at higher frequencies with a frequency^{−4} energy rolloff (e.g., Komen et al. 1994).

The surface elevation bispectrum *B*(*ω*_{1}, *ω*_{2}) was predicted in both the Lagrangian reference frame of a surface-following buoy and the Eulerian reference frame of a fixed horizontal location by substituting the observed *E*(*ω*, *θ*) in Eq. (15), using the appropriate coupling coefficients and , respectively. These predicted bispectra (transformed to cyclic frequencies) are compared with the observed bispectrum, estimated from the vertical buoy displacement time series in Fig. 3 [only the real part of *B*(*ω*_{1}, *ω*_{2}) is shown; the imaginary part theoretically vanishes]. As expected from the properties of the coupling coefficients, discussed earlier, the Eulerian and Lagrangian predictions differ dramatically. The predicted Eulerian bispectrum (bottom panel) shows the familiar pattern (Hasselmann et al. 1963; Elgar and Guza 1985) of positive values at frequencies *f*_{1}, *f*_{2} > 0.08 Hz along ridges with *f*_{1} or *f*_{2} close to the peak frequency (0.09 Hz), indicating the coupling between the primary spectral peak components and in-phase harmonic components. At lower frequencies these ridges change sign to smaller negative values owing to the coupling of primary waves and 180° out-of-phase bound infragravity components. In contrast, in the predicted Lagrangian bispectrum, the negative values at infragravity frequencies are replaced with much larger positive values, and bispectral levels are weak at higher frequencies, consistent with the absence of harmonics in these relatively deep-water conditions (*kh* ≈ 5.1 at the peak frequency) and more energetic in-phase infragravity components (Fig. 1).

The main features of the observed bispectrum are clearly in agreement with the Lagrangian theory prediction of nonlinear interactions dominated by energetic phase-coupled infragravity motions. Small differences between the Lagrangian theory prediction and the observed bispectrum may be due to several possible sources of errors, including the mooring response of the buoy (not accounted for in the prediction for an idealized free-floating buoy), the lack of a high-resolution directional spectrum estimate, and the finite record length (bispectral estimates have much greater statistical uncertainty than ordinary spectral estimates). The high degree of similarity in the observed bispectrum and Lagrangian prediction is in sharp contrast with the completely different structure of the Eulerian prediction, and thus this comparison demonstrates the large differences between Eulerian and Lagrangian sea surface elevation records in a natural wind sea.

The Lagrangian bispectral analysis applied here to a moored buoy can readily be extended to account for the slow advection by a constant surface current **U** of a drifting buoy. In general, for every triad consisting of two primary waves and a secondary wave, the triad frequencies and wavenumbers obey the interaction rules *ω*_{1} + *ω*_{2} + *ω*_{3} = 0 and **k**_{1} + **k**_{2} + **k**_{3} = 0. In the case of a drifting buoy, the relative (Doppler shifted) frequencies *ω*_{i,r} = *ω*_{i} − **k**_{i} ⋅ **U**(*i* = 1, 2, 3) obey the same interaction rule *ω*_{1,r} + *ω*_{2,r} + *ω*_{3,r} = 0. Hence, a bispectral analysis applied to records collected with a drifting buoy is expected to yield results that are nearly identical to those of a moored buoy but projected in a slightly distorted relative frequency frame of reference.

### b. Sea surface skewness

Bispectra can provide detailed insight in nonlinear coupling of wave components in the spectral domain, but the interpretation is often complex and these higher-order spectra are cumbersome to use in the analysis of large datasets. A useful bulk measure of the deviation from Gaussian statistics in a surface elevation record *z*_{b}(*t*) is the third moment , which—by definition—equals the integral of the bispectrum over all *ω*_{1}, *ω*_{2} frequency pairs:

and, using Eq. (15), can be expressed in terms of the primary ocean wave spectrum

Since the Lagrangian contribution *L*^{η} [Eq. (8)] to [Eq. (7)] obeys the antisymmetry relation [Eq. (9)], the Lagrangian contribution to the bulk integral [Eq. (18)] is antisymmetric about both the *ω*_{1} and *ω*_{2} axes and vanishes altogether in the integration across the entire *ω*_{1}–*ω*_{2} plane. Hence, although a Lagrangian surface elevation record *z*_{b}(*t*) may look very different from the Eulerian surface *η*(*t*) at a fixed horizontal position, the corresponding third moments and are in theory equal and given by Eq. (18) or alternatively the same integral with the Eulerian coupling coefficient .

The third moment is conveniently normalized to a skewness measure

a parameter that is often used to quantify wave nonlinearity. To evaluate the sea surface skewness in natural wind sea and swell conditions and to examine the fidelity of estimates obtained with operational moored waverider buoys, we used long-term observations from two buoys in the Coastal Data Information Program (CDIP) network (https://cdip.ucsd.edu/). The selected buoys are the Point Reyes Buoy, located in deep water (575 m) off Point Reyes, California, and the nearby San Francisco Bar Buoy, located in shallow water (15 m) on the ebb tidal shoal in the entrance to San Francisco Bay. Both buoys are Datawell Mark III Directional Waverider buoys. A 3-month-long period from 1 November 2009 to 31 January 2010, with a representative range of Pacific swell and local wind-sea conditions, was selected for analysis. On several occasions the significant wave height recorded by the Point Reyes Buoy exceeded 6 m (Fig. 4, upper panel).

The buoy data were processed in 4-h-long records, discarding any data with gaps or anomalous spikes. Spectra and bispectra were computed using the same procedures discussed earlier for the buoy off Bodega Bay. However, since the CDIP buoys are equipped with an internal filter that removes signals with periods longer than 30 s, the skewness estimates presented here are restricted to a frequency range of 0.03–0.64 Hz that excludes the lower part of the infragravity range. To predict the skewness for each data record, first an estimate of *E*(*ω*,*θ*) was obtained with the MEM method from the buoy measurements, and this estimate was used with Eq. (15) to predict the bispectrum *B*(*ω*_{1}, *ω*_{2}). Finally, the observed and predicted skewness values were obtained by integrating the corresponding bispectra [Eq. (17)] within the same restricted frequency range and normalizing Eq. (19) with the measured variance, also in the same frequency range.

At both sites the predicted skewness values are positive, ranging from 0 to 0.1 at the deep site and from 0 to 0.7 at the shallow site (Fig. 4). These differences are indicative of the much stronger nonlinearity in shallow water (e.g., Elgar and Guza 1985). At both sites the observed and predicted skewness values are generally in good agreement. It should be noted, however, that these skewness values (both observed and predicted) do not include coupling to lower infragravity (<0.03 Hz) frequencies. Predicted skewness values (not shown) that include these lower frequencies are on average 90% higher at the deep site and 25% higher at the shallow site. Thus, while the encouraging agreement between observations and predictions suggests that surface-following buoy measurements of ocean surface waves can provide quantitative estimates of sea surface skewness, care must be taken to resolve the infragravity band, especially in deep water where the dominant contributions to the skewness come from the coupling between the primary sea–swell waves and infragravity components (e.g., Fig. 3). To accurately account for these low-frequency contributions to the bispectrum may require some refinement in the design and operation of moored buoys. In particular, attention should be paid to minimizing the effects of the mooring and the implications of built-in filters on the low-frequency response, and the collection of longer continuous data records than the customary 20–30-min records used in routine wave observations should be considered. Since the low-frequency response is critical for estimating the sea surface skewness from buoy records, free-drifting buoys potentially provide an attractive platform for this application by eliminating any possible distortion from mooring forces.

## 4. Stokes drift fluctuations

Whereas moored buoys are widely used to collect surface wave measurements, small drifting buoys can provide measurements of both waves and surface currents (e.g., Herbers et al. 2012; Thomson 2012; Pearman et al. 2014). The concurrent observation of waves and surface drift is of particular interest in surface dispersion and mixing studies because traditional Eulerian current measurements are difficult to make at the sea surface, and the wave-induced Stokes drift has a completely different profile in the more natural Lagrangian reference frame (e.g., Phillips 1977). Whereas the mean Stokes drift has been the topic of numerous studies (e.g., Hasselmann 1970; Xu and Bowen 1994; Polton et al. 2005; Lentz et al. 2008; Aiki and Greatbatch 2012), infragravity fluctuations in the wave-induced surface drift on the scale of wave groups have received almost no attention (e.g., Smith 2006). Here, we examine the fluctuating surface drift observed with an idealized drifting buoy in a random sea state using second-order wave theory. For simplicity, we consider a steady and homogenous wave field in the absence of any external forcing.

An expression for the Lagrangian (horizontal) surface velocity **u** at the drifter location **x**_{b}(*t*) can be derived by evaluating the horizontal momentum equation at the sea surface. Neglecting vertical shear and Coriolis effects, the flow is driven by the horizontal pressure gradient at the surface that includes a nonhydrostatic contribution:

Equation (20), derived from the nonlinear surface boundary conditions and the momentum equations, is a fully nonlinear relation between the surface elevation *η* and surface drift **u**.

Substitution of the surface elevation function [Eq. (1)] in Eq. (20) and integrating with respect to time, yields to second order

The integration constant **U**_{E} accounts for the ambient Eulerian surface current (e.g., the return flow driven by Stokes–Coriolis forcing and tidal- and wind-driven flow; see Lentz et al. 2008). The Lagrangian Stokes drift contribution to the mean surface current is implicitly included in the nonlinear interaction term (i.e., the zero-frequency contribution of self–self-difference interactions). The coupling coefficient is given by

The first term on the right-hand side of Eq. (22) describes the Lagrangian effect of the horizontal wave orbital excursions on the surface drift, whereas the second and third terms are quasi-Eulerian second-order contributions driven by nonhydrostatic pressure contributions and pressure gradients resulting from the second-order surface slopes, respectively. Our results [Eqs. (21) and (22)] are obtained from an expansion around the moving sea surface *z* = *η*, and although expressed different algebraically, they are in exact agreement with Eqs. (2.12) and (2.13) of Herterich and Hasselmann (1982), who expanded around the mean sea surface *z* = 0.

To illustrate the low-frequency drift fluctuations and connect this analysis with the classical steady Stokes drift solution for a plane wave, consider a pair of unidirectional wave components (*θ*_{1} = *θ*_{2} = 0) in deep water with frequencies *ω*_{1,2} = *ω* ± Δ/2 and the same amplitude = *a*/2. The interaction coefficients [Eq. (22), using Eqs. (10) and (11)] simplify to

Again, as for the surface elevation, in deep water the double-frequency harmonics vanish in the surface drift, and when the difference frequency Δ is small, the negative Eulerian contribution (−2*ω*^{2}Δ/*g*) to the infragravity response [Eq. (23b)] is *O*(Δ/*ω*) smaller than the positive Lagrangian contribution.

For this idealized wave field with modulated wave groups, Eq. (21) [using Eqs. (23a) and (23b)] reduces to

where, for simplicity, we set the arbitrary initial drifter position **x**_{0} = 0 and neglected the small Doppler shift induced by the Stokes drift in the phase function Eq. (5).

The difference frequency variations in the surface drift predicted by Eq. (24) are in agreement with the classical Stokes drift theory. In the limit of small Δ, the bichromatic wave field can be approximated as a sinusoidal wave train of frequency *ω* with an amplitude slowly varying between a maximum value 2*a* in the center of the wave groups and vanishing amplitude in between the groups. The classical steady Stokes drift prediction locally associated with this modulated wave field (not considering the return flow driven by the divergent mass flux) yields a surface drift velocity varying between 0 (in between the groups) to 4*ω*^{3}*a*^{2}/*g* (center of the groups), which—in the limit of small Δ—is in exact agreement with Eq. (24).

Dynamically, the mean contribution to the Stokes drift [the first term in Eq. (24)] is affected by the Coriolis force, resulting in Eulerian counterflows that in theory—if enough time is available for a stationary solution to develop—cancel the mean Stokes drift (see, e.g., Hasselmann 1970; Polton et al. 2005). However, the infragravity fluctuations, with periods small compared to the inertial period, are not affected by Earth’s rotation and will thus not be balanced by Eulerian counterflows. Hence, the Lagrangian second-order theory suggests the presence of energetic fluctuations in the surface drift on the scale of wave groups that are much larger than the Eulerian bound wave velocities and should be readily detectable in drifter records.

To determine whether such Stokes drift fluctuations are indeed present in a natural sea state, we analyzed data from free-drifting buoys deployed in deep water about 60 km offshore of Monterey Bay, California. The sea state was a mix of swell and locally generated wind sea from the northwest with a significant wave height of about 3.3 m. The buoys include a Datawell DWR-G7 Directional Waverider buoy and three small wave-resolving drifters (WRD) equipped with GPS and accelerometers [see Pearman et al. (2014) for a detailed description of these drifters]. The drifters were deployed in a cluster and allowed to drift for about 8 h before retrieval. Drift velocity measurements, based on the Doppler shift in the GPS signal, were recorded on board the drifters. An example time series of measured velocities from one of the drifters is shown in Fig. 5 (upper panel). The observed orbital velocities (green curve, bandpassed in the 0.05–0.5-Hz swell–sea frequency range) show the expected modulations on (infragravity) time scales of a few minutes, characteristic of wave groups. The corresponding record of infragravity velocity fluctuations (blue curve, bandpassed in the 0.002–0.02-Hz range) shows a correlation with the wave groups with the maxima in the center of the wave groups. This pattern is qualitatively consistent with the in-phase infragravity fluctuations predicted by the Lagrangian theory, in contrast to the 180° phase difference of the Eulerian infragravity bound waves. The observed infragravity drift variations are as large as 20 cm s^{−1} and comparable to the order of magnitude of the theoretical mean Stokes drift.

The in-phase coupling between infragravity drift fluctuations and surface wave groups is also clear in the observed bispectrum of the velocity component in the dominant wave direction (Fig. 5, bottom panel). The bispectrum was estimated from the entire 8-h-long record based on 27.3-min-long FFT segments and merging 13 × 13 bands. To further reduce the statistical uncertainty, the resulting bispectrum (resolution 0.0079 Hz) was ensemble averaged over the three drifters. Similar to the surface height bispectrum presented earlier (Fig. 3), the real part of the velocity bispectrum shows positive values for pairs of frequencies (*f*_{1}, *f*_{2}) with one component in the swell–sea band and the other in the infragravity band. The imaginary part of the bispectrum shows no detectable coupling. This observed near-zero biphase is consistent with the real and positive coupling coefficient [Eq. (23b)] and clearly demonstrates the in-phase coupling between infragravity drift fluctuations and surface wave groups.

Although the observed phase-coupled infragravity motions indeed exhibit the phase characteristics of the predicted Stokes drift fluctuations, other motions (e.g., free infragravity waves or turbulence) may contribute to the infragravity velocity field that are not coupled to the local swell–sea waves and thus do not contribute to the bispectrum. To compare the spectral energy levels of observed infragravity velocities and predicted Stokes drift fluctuations, high-resolution spectra of the velocity components *u* (in the dominant wave direction) and *υ* (in the transverse direction) were computed from a 7.28-h-long data record for all three drifters, using 1.82-h-long Hamming-windowed FFT segments with 50% overlap and merging three bands to yield a resolution of 0.0005 Hz. The resulting *u* and *υ* spectra, averaged over the three drifters, are shown in Fig. 6, together with the total velocity spectrum. The observed spectra are relatively flat at infragravity frequencies with some polarization along the dominant wave direction. To compare these observed spectra with the theoretically expected infragravity Stokes drift fluctuations, we simplify the general expression of the Lagrangian velocity field [Eq. (21)] by neglecting the effects of directional spreading and the Doppler shift associated with the mean current [Eq. (5)] to form a spectrum of the second-order velocity fluctuations:

where *E*^{η}(*ω*) is the frequency spectrum of primary waves and the coupling coefficient can be approximated for a narrow primary wave spectrum with Eq. (23). Similarly, the Eulerian velocity spectrum of second-order bound waves can be approximated by replacing in Eq. (25) with the corresponding Eulerian coefficient . Predictions of and at infragravity frequencies, using the Datawell buoy surface height record to estimate the primary wave spectrum *E*^{η}(*ω*), are included in Fig. 6 (the displayed spectra are single-sided *f*-spectra). Whereas the spectral levels predicted by the Lagrangian theory (blue curve) are similar to the observed spectral levels (albeit slightly lower), the Eulerian bound wave spectral levels (red curve) are several orders of magnitude smaller. Moreover, both the observed and theoretical Lagrangian infragravity drift spectra are nearly white, whereas the Eulerian spectrum drops off sharply toward low frequencies.

Although a precise quantitative comparison will need to account for directional spreading effects and also take into consideration the presence of free infragravity waves radiated from shore (e.g., Herbers et al. 1995), these preliminary results indicate that the infragravity Stokes drift fluctuations are indeed important in natural ocean waves, and these motions are much more energetic than the Eulerian bound infragravity wave field. Whereas recent field studies show the mean Stokes drift to be approximately canceled by an Eulerian return flow (Lentz et al. 2008), consistent with a balance between the Stokes–Coriolis wave stress and the Coriolis force acting on the Eulerian return flow (Hasselmann 1970), such a balance is not expected for the infragravity Stokes drift fluctuations that have small periods compared to the (inertial period) time scale of the Coriolis adjustment. The large observed drift fluctuations (Fig. 5) are consistent with the absence of a balancing Eulerian flow and may be important in upper-ocean mixing and dispersion of pollutants.

## 5. Conclusions

The Lagrangian properties of ocean surface waves are important to studies of upper-ocean mixing, the interpretation of remote sensing data, and the in situ sensing of waves and currents with drifting buoys. Although the theory for second-order Lagrangian wave properties is well established from earlier studies (notably, Srokosz and Longuet-Higgins 1986; Longuet-Higgins 1986; Herterich and Hasselmann 1982), their dynamics in a natural sea state are not well understood. In this study, we evaluate the spectral properties of Lagrangian surface height and drift measurements from second-order wave theory, consider skewness estimates from Lagrangian records, identify energetic infragravity modulations in surface Stokes drift, and compare predictions with field observations of moored and free-drifting buoys.

In a Lagrangian reference frame, high-frequency second-order bound waves are effectively shifted to infragravity frequencies, and a Lagrangian wave record that resolves these lower-frequency bound wave contributions has in theory the same skewness as an Eulerian wave record. Skewness estimates obtained from long-term moored buoy observations in deep and shallow water are in good agreement with the theoretical predictions. These results suggest that moored buoy networks, which are widely used to collect routine wave observations, can also provide reliable estimates of sea surface skewness, a parameter that plays an important role in the sea state bias of remote sensing systems and nearshore sediment transport.

Whereas the Lagrangian motion of a surface-following buoy in deep water exactly cancels the Eulerian sum frequency bound waves, no such cancellation occurs at difference (infragravity) frequencies. Instead, the setdown under wave groups is replaced with a much larger setup signal in a Lagrangian wave record. In deep water, the Lagrangian contributions dominate the infragravity wave signal, whereas in shallow water, where Eulerian bound waves approaching resonance are amplified, Lagrangian distortions are generally relatively small. Bispectral analysis of moored buoy observations in relatively deep water confirms both the suppression of double-frequency harmonics and strong in-phase coupling at infragravity frequencies predicted by the theory.

Lagrangian effects manifest themselves in a similar fashion in surface drift fluctuations, with infragravity variations that are in phase with the wave groups and much larger than the (180° out of phase) Eulerian infragravity bound wave velocities. Drifter observations confirm the presence of such energetic infragravity fluctuations [*O*(10) cm s^{−1}] that are an order of magnitude larger than the predicted Eulerian velocities.

These infragravity modulations in the wave-induced surface drift, which to our knowledge have not been explicitly identified before, may be important to upper-ocean mixing and diffusion processes.

## Acknowledgments

This research is supported by the Office of Naval Research (Littoral Geosciences and Optics Program and the Physical Oceanography Program) and the National Science Foundation (Physical Oceanography Program). We greatly appreciate the contributions of Paul Jessen, LCDR Douglas Pearman, LCDR Steve McIntyre, and Erik van Ettinger to the field experiments. The captains and crews of the R/V *Gordon Sproul* and the R/V *Point Sur* provided excellent support. The Point Reyes and San Francisco Bar buoy data were generously provided by the Coastal Data Information Program (https://cdip.ucsd.edu/).