Turbulence-induced anti-Stokes flow: experiments and theory

Abstract

We report experimental evidence of an Eulerian-mean flow, $\overline {u}(z)$, created by the interaction of surface waves and tailored ambient sub-surface turbulence, which partly cancels the Stokes drift, $u_s(z)$, and present supporting theory. Water-side turbulent velocity fields and Eulerian-mean flows were measured with particle image velocimetry before vs after the passage of a wave group, and with vs without the presence of regular waves. We compare different wavelengths, steepnesses and turbulent intensities. In all cases, a significant change in the Eulerian-mean current is observed, strongly focused near the surface, where it opposes the Stokes drift. The observations support the picture that, when waves encounter ambient sub-surface turbulence, the flow undergoes a transition during which Eulerian-mean momentum is redistributed vertically (without changing the depth-integrated mass transport) until a new equilibrium state is reached, wherein the near-surface ratio between $|{\rm d}\overline {u}/{\rm d}z|$ and $|{\rm d}u_s/{\rm d} z|$ approximately equals the ratio between the streamwise and vertical Reynolds normal stresses. This accords with a simple statistical theory derived here and holds regardless of the absolute turbulence level, whereas stronger turbulence means faster growth of the Eulerian-mean current. We present a model based on Rapid distortion theory which describes the generation of the Eulerian-mean flow as a consequence of the action of the Stokes drift on the background turbulence. Predictions are in qualitative, and reasonable quantitative, agreement with experiments on wave groups, where equilibrium has not yet been reached. Our results could have substantial consequences for predicting the transport of water-borne material in the oceans.

1. Introduction

The phenomenon of Stokes drift implies that periodic water waves in irrotational flow induce a net Lagrangian-mean transport along their direction of propagation (van den Bremer & Breivik 2017). Although discovered theoretically a long time ago by Stokes (1847), Stokes drift has been elusive in laboratory experiments. A particular difficulty is in separating it from Eulerian-mean flows whose properties are determined by the boundary conditions of the flume itself (Monismith 2020), although recent investigations in wave flumes where waves (or wave groups) propagate on initially quiescent water have provided convincing evidence (Grue & Kolaas 2017; van den Bremer et al. 2019).

Floating and suspended matter of small enough size, such as microplastics (van Sebille 2020), oil spills (Boufadel et al. 2021), plankton (Hernandez-Carrasco et al. 2018), larvae and nutrients (Röhrs et al. 2014), are transported in the oceans with the Lagrangian-mean current $\bar {u}_L$ , equal to the Eulerian-mean current, $\bar {\boldsymbol{u}}(z)$ , plus the Stokes drift $\boldsymbol{u}_s$ , i.e. $ \bar {\boldsymbol{u}}_L(z) = \bar {\boldsymbol{u}}(z) + \boldsymbol{u}_s(z)$ . The importance of correctly modelling ocean transport makes prediction of the Lagrangian current a pivotal, but still open, question. For a monochromatic (nearly) linear wave of wavenumber $k_0$ and constant amplitude $a$ propagating in deep water the Stokes drift velocity along the direction of wave propagation is

(1.1) \begin{equation} u_s(z) = (ak_0)^2c(k_0){e}^{2k_0 z}, \end{equation}

where the intrinsic phase velocity (i.e. in the reference system where the surface is at rest) is $c$ . The phase velocity and group velocity, $c_g$ , are given by

(1.2) \begin{equation} c(k_0) = \sqrt {g/k_0} \text{ and } c_g(k_0) = \frac {1}{2}\sqrt {g/k_0}, \end{equation}

where $g$ is the gravitational acceleration. A contribution to oceanic transport velocities of magnitude $u_s$ can be highly significant and must be taken into account when predicting, e.g. the development of oil spills or the fate of microplastic particles (Hackett et al. 2006; Onink et al. 2019; Cunningham, Higgins & van den Bremer 2022).

The naïve approach for the modeller is to obtain the Lagrangian-mean current by adding the Stokes drift calculated from wave spectra to the Eulerian flow, which is assumed to be unaffected by waves. However, several field studies have observed that waves appear to have no discernible effect on the Lagrangian-mean current, contrary to theory (Smith 2006; Lentz et al. 2008). Smith (2006) found that even short wave groups experience an Eulerian-mean current, which acts to entirely cancel the Stokes drift at the surface, and that the counter-current is strongly correlated with the presence of Stokes drift – it appears only when the wave group is present and disappears once it has passed. Other observations find that the inclusion of Stokes drift does improve results, however, e.g. Röhrs et al. (2012), who used drifters in coastal waters and employ a coupled wave–ocean model.

Experiments of waves propagating on currents have also yielded results which are inconsistent with a simple addition of Stokes drift. In a careful laboratory study, Monismith et al. (2007) found no change in Lagrangian-mean flow when waves were added, i.e. the Stokes drift is locally cancelled by an equal and opposite Eulerian flow. Moreover, reanalysis of previous experiments by Swan (1990a , b ), Jiang & Street (1991) and Thais (1994) supported the same conclusion (the observation seems to have been made also in the PhD work of Cowen 1996). These results have remained something of a puzzle, as ‘[n]o existing theory of wave–current interactions explains this behaviour’ as Monismith et al. (2007) put it. Here, we demonstrate a mechanism that could resolve this conundrum at least in part.

While being careful not to draw definite conclusions concerning the above-mentioned results, one might remark that some level of turbulence was likely to have been present alongside the waves in all these experiments. The flow of Swan (1990b ) passed through a honeycomb flow straightener which will have generated significant turbulence levels, and the mean shear of the flow would provide further turbulence production. Jiang & Street (1991) discuss wave–turbulence interactions in their experiments in detail (further detailed by Cheung & Street 1988), as do Thais & Magnaudet (1996) for the experiment of Thais (1994). The experiment of Swan (1990a ) did not have imposed wind or current. Since measurements were made after the wavemaker had run for many hours, one may speculate that turbulence might well have been present, created, e.g. in the oscillating bottom boundary layer (the depth was less than a third of a wavelength) or thermal convection, with which the waves would have had ample time to interact, the flume being $18$ m long. Naturally, this can be no more than speculation so long afterwards.

There have been indications that the interaction between waves and pre-existing turbulence will result in an alteration of the Eulerian current. Waves have the effect of reorienting and intensifying the turbulence beneath, as predicted theoretically (Magnaudet & Masbernat 1990; Teixeira & Belcher 2002), and confirmed numerically (e.g. Guo & Shen 2013; Tsai et al. 2017; Xuan et al. 2019, 2020, 2024), experimentally (e.g. Bliven et al. 1984; Cheung & Street 1988; Thais & Magnaudet 1996; Savelyev, Maxeiner & Chalikov 2012; Smeltzer et al. 2023) and in field studies (Cavaleri & Zecchetto 1987; Qiao et al. 2016). Langmuir turbulence, the disordered pattern of long rolls approximately aligned with wind and waves due to Langmuir circulation formation at sea (McWilliams, Sullivan & Moeng 1997), was observed by Plueddemann et al. (1996) to persist for up to a day after the wind had stopped, sustained by the surface waves the wind had created.

Pearson (2018) predicted with a theoretical argument that the interaction between waves and ambient turbulence would, on average, produce an Eulerian-mean current which opposes the Stokes drift near the surface. His paper has received little attention up until now, but in our theoretical work later, we shall draw heavily on his work and apply it to our own settings. Pearson’s prediction shows that the Eulerian-mean current, $\bar {u}(z)$ , is opposite and similar to $u_s(z)$ near the surface but changes sign beneath the wave-influenced surface layer and integrates to zero as a function of depth. Although the wave–turbulence-induced Eulerian-mean current incurs no net change in mass transport, it will partly cancel the Stokes drift at the surface, and we therefore refer to it as an ‘anti-Stokes’ current.

The picture which emerges is that, when waves and turbulence first meet, the combined flow goes through a transient ‘spin-up’ period before a new quasi-equilibrium is reached, which includes the Eulerian anti-Stokes current. In our set-up, turbulence is created at the inlet of the mean flow and as it moves downstream towards the point of measurement it encounters waves and develops together with the velocity profile. When the time of effective interaction is short, the result of partial spin-up is captured, while turbulence moving in the presence of waves for sufficiently long could reach a quasi-steady state before it is measured. The conditions for ‘short’ and ‘long’ interaction time are satisfied in our experiments, for wave groups (Experiment 1) and regular waves (Experiments 2 and 3), respectively.

We review and extend statistical theory for both the transient and steady stages, and derive a theory based on Rapid distortion theory (RDT) which captures the underlying physics of the ‘spin-up’ period and shows that the Eulerian acceleration of the current will exhibit approximately the same depth-dependent behaviour as the resulting current that we observe. The process is intimately related to the so-called ‘CL2’ mechanism which creates Langmuir circulation (Craik & Leibovich 1976). It is worth noting that the CL2 mechanism, as lucidly reviewed by Leibovich (1983), requires the presence of an Eulerian-mean flow with slope (i.e. ${\rm d}\bar {u}/{\rm d} z$ ) of the same sign as that of the Stokes drift profile, whereas the anti-Stokes current induced by wave–turbulence interaction, $\bar {u}(z)$ , has opposite slope. Thus, $({\rm d}\bar {u}/{\rm d} z)\boldsymbol{\cdot }({\rm d} u_s/{\rm d} z) \lt 0$ , which implies that the induced current tends to stabilise the combined system with respect to the CL2 instability.

1.1. Possible previous observations

Indications are that the same physical phenomenon has been at play in wave–current experiments performed in the context of studying the bottom boundary layer in shallow wave–current flows, motivated by understanding sediment transport (thus not highlighting the significance for ocean modelling). A string of independent measurements of the mean Eulerian flow in the presence and absence of waves by van Hoften & Karaki (1977), Bakker & van Doorn (1978), van Doorn (1981), Kemp & Simons (1982, 1983), Rashidi, Hetsroni & Banerjee (1992), Klopman (1994), Mathisen & Madsen (1996) and Singh & Debnath (2017), all showed that the waves caused a significant alteration of the mean flow near the surface, adding a contribution in the direction opposite to wave propagation in the near-surface region. More recent experiments also report the same (Zhang & Simons 2019; Peruzzi et al. 2021). In addition to making the same observation, Umeyama (2005, 2009) found in his experiments that the vertical structure of the streamwise–vertical Reynolds shear stress depends strongly on the wave-propagation direction in the near-surface region.

Since these studies considered shallow currents where waves affect the bottom boundary layer, a direct comparison with our experiment in deep water is dubious, yet subsequent theoretical analysis gives reason to suspect a close connection. Nielsen & You (1996) and Dingemans et al. (1996) propose two early explanations for the difference in Eulerian-mean current; the former relies on a force balance on average including the mean stress from waves and turbulence represented by eddy viscosity, the latter on the creation of streamwise rolls due to the Craik–Leibovich vortex force due to the sidewall boundary layers, whose presence was observed by Klopman (1994). Groeneweg & Klopman (1998) developed a more sophisticated theory based on Generalised Lagrangian Mean (GLM) theory, similar in spirit to the physical process we consider herein, but their analysis is not easy to compare with ours since it involves a set of coupled nonlinear differential equations and a complex turbulence model. Interestingly, Groeneweg & Battjes (2003) reconcile all three descriptions, at least qualitatively, by extending their GLM theory to three dimensions. Huang & Mei (2003) provided a careful theory, once more primarily concerned with the effect of waves on the bottom boundary layer, but also shedding light on the observed changes near the surface. They, too, use the simple eddy-viscosity model of turbulence. In particular, they conclude that the mean wave-induced shear stress near the free surface is opposite in direction to the wave propagation, and is ‘due largely to the distortion of eddy viscosity near the surface’. With a simple mixing-length model, Umeyama (2005, 2009) as well as Yang et al. (2006) find reasonable agreement with experiments. Olabarrieta, Medina & Castanedo (2010) devise a simplified numerical model to avoid the restriction of low-steepness waves, and the perturbation theory of Tambroni, Blondeaux & Vittori (2015) yields a model able to predict the Eulerian-mean velocity profile throughout the water column; both of the latter employ depth-dependent eddy-viscosity models. Crucially, all of these many model explanations depend on the simultaneous presence of waves and turbulence to explain the change in mean flow.

When reviewing previous numerical studies of waves in the presence of turbulence, one can also find evidence of the same Eulerian-mean current creation that we observe experimentally, even though the authors themselves have not discussed its significance especially. Borue, Orszag & Staroselsky (1995) merely remark that the change in mean current near the surface should be studied further, while Kawamura (2000) notes that the changes in current are higher the larger the Stokes drift magnitude, and discusses consequences for (Langmuir) turbulence production. Also Fujiwara, Yoshikawa & Matsumura (2020) find Eulerian-mean velocity profiles under Langmuir turbulence which are qualitatively very similar to those we measure under regular waves (their figures 6a and 9a) and report in § 2.1.2, but do not discuss this point particularly.

1.2. Other wave-driven currents

Eulerian-mean flow driven by waves also occurs without pre-existing turbulence or vorticity from two separate mechanisms. First, to compensate for the divergence of Stokes drift on the group scale, Stokes drift in otherwise quiescent water must also be accompanied by an Eulerian return flow (e.g. Longuet-Higgins & Stewart 1962; van den Bremer & Taylor 2016) for wave groups. In deep water, the depth-integrated Stokes drift and the depth-integrated Eulerian return flow are equal and opposite, and mass is preserved globally, not locally. This phenomenon is reviewed further in § 4.1. The Stokes drift profile is highly concentrated near the surface whereas the return flow varies slowly with depth. The Eulerian-mean ‘anti-Stokes’ current we observe varies rapidly with depth and cannot be explained by this mass conservation mechanism. Second, there is also surface streaming driven by viscosity confined to a thin viscous boundary layer beneath the surface, resulting from the imparting of wave momentum to the fluid as the waves decay (Longuet-Higgins 1953; Craik 1982). Tsai et al. (2017) and, recently, Fujiwara (2024), studied how this current can also interact with waves to generate small-scale turbulence of Langmuir type via the CL2 mechanism. The viscous sub-surface layer is likely thinner than what our measurements can resolve, and, additionally, the resulting Eulerian-mean flow is directed along, not against, the direction of wave propagation. Third, the Earth’s rotation causes a wave-induced Eulerian-mean flow that can exactly cancel the Stokes drift (for periodic waves and in the absence of viscosity), which is also known as the anti-Stokes flow (Hasselmann 1970). While this rotation-induced anti-Stokes flow could explain some field observations (Lentz et al. 2008; Röhrs et al. 2012), it cannot explain experimental results discussed above and presented herein as our characteristic time scales are vastly smaller than Earth’s period of rotation (the Rossby number is large). Hence, neither of these three mechanisms can explain the observations just mentioned, nor those by, e.g. Monismith et al. (2007), or indeed those we report herein.

2. Experimental methods

We report on measurements performed during three experimental campaigns in the water channel facility at NTNU Trondheim, shown in figure 1(a). A pump system circulates water through the test section of dimensions $11.2$ m $\times$ $1.8$ m $\times$ $1.0$ m (length $\times$ width $\times$ height).

Figure 1. Experimental set-up: (a) side view of water channel with flow from left to right and field of view (FOV) indicated with a green rectangle, (b) top view of measurement region for the stereo particle image velocimetry (PIV) set-up in Experiment 1, including positions of wave probes (WPs) and laser-induced fluorescence (LIF) camera for surface detection; (c) longitudinal view of the planar PIV set-up in Experiments 2 and 3. Experiment 2 employed three stacked cameras as shown, whereas in Experiment 3, a single PIV camera was used.

An active grid at the inlet of the test section allowed the turbulence to be generated and varied. The grid consists of square wings measuring $10$ cm across the diagonal, attached to $18$ vertically and $10$ horizontally oriented bars, each controlled independently by a stepper motor. Several different active-grid actuation cases were investigated, listed in Appendix A. The grid wings were rotated with random rotational velocity, acceleration and period within set limits (Hearst & Lavoie 2015; Smeltzer et al. 2023), or in one case flapped back and forth between two positions at irregular time intervals. The instantaneous rotation frequency of the grid wings varied about a mean active-grid frequency $\overline {f_G}$ by $\pm 0.5\overline {f_G}$ with a top-hat distribution. In Experiments 1 and 2 (see §§ 2.1.1 and 2.1.2) a surface plate was mounted from the grid that extends downstream approximately $1$ m to dampen surface disturbances produced by the grid. A diagram of the set-up is shown in figure 1, and further details are given by Jooss et al. (2021).

The bottom and sidewall boundary layers are thin enough not to reach the central part of the flow. They reach a maximum of approximately $20$ cm at the point where velocity measurements were made, and are thinner upstream where waves and turbulence interact, hence essentially all turbulence which affects our experimental results has been generated by the grid and not boundary layers. The upstream distance within which waves and turbulence interact is denoted $L_{\textit{FOV}}$ and is tabulated in table 1.

Table 1. Overview of the three experiments. Here, $h$ is the mean water depth, $f_0$ is the wave frequency in the laboratory frame (wavemaker frequency), $f_{\textit{ac}}$ is the PIV acquisition frequency, $N_{\textit{ens}}$ is the number of repetitions per case, $T_{\textit{PIV}}$ is the duration of each acquisition interval and $L_{\textit{FOV}}$ is the distance the flow travels as a free stream upstream of the measurement field of view. $\dagger$ : except case 1.C.2 where $ N_{\textit{ens}}=20$ .

Table 2. Measured flow quantities: mean current $U_0$ , carrier-wave number $k_0$ , measured group temporal width in the lab frame $\tau$ , root-mean-square turbulent velocity components (subscript ‘rms’), turbulent kinetic energy $e$ (calculated for Experiments 2 and 3 as discussed in main text), and streamwise–streamwise integral scale $L_x^x$ . Steepness is given as $k_0a_{{p}}$ for Experiment 1 where $a_{{p}}$ is peak amplitude, and $ak_0$ in Experiments 2 and 3. Where several values of steepness are listed (cases 1.C and 2.A–3.B), these are referred to elsewhere as 1.C.1, 1.C.2, 2.A.1.1, 2.A.1.2, etc.

A plunger wavemaker at the downstream end of the test section ( $10.2$ m from the active grid) was used to generate waves propagating upstream on the current. Waves with a group velocity lower than the mean flow were unable to propagate upstream, thus preventing high-frequency wave noise and unwanted free harmonics (parasitic waves) from the wavemaker from entering the test section. Wave properties were obtained from surface elevation measured by a pair of resistive wave probes (HR Wallingford) near the measurement position.

As indicated by the coordinate system in figure 1(a), the waves propagate in the positive $x$ -direction, while the flow is in the negative $x$ -direction. The mean free-surface level is at $z=0$ .

2.1. Experimental campaigns

Three separate experimental campaigns were conducted between 2020 and 2023. Experiment 1 was also reported in Smeltzer et al. (2023), where the focus was on the change in turbulent enstrophy and wave scattering.

The three experiments investigate essentially the same phenomenon, but with fundamental differences in experimental design and the nature of the acquired measurement data. A summary is provided in table 1.

2.1.1. Experiment 1: wave groups

Cases 1.A-1.D in tables 1–3 are from Experiment 1. Wave groups were generated that propagated upstream atop the turbulent flows. The water depth $h$ was $0.4$ m. The velocity field was measured using stereoscopic particle image velocimetry (SPIV), measuring all three velocity components in a plane perpendicular to the mean flow located a distance $8.38$ m downstream from the active grid (i.e. $83.8$ grid units), and $L_{\textit{FOV}}=7.38$ m downstream of the trailing edge of the surface plate. Two 25-megapixel cameras were mounted on either side of the test section, viewing the field of view at $\pm 45^\circ$ to the $x$ -axis as shown in figure 1 b). The field of view was $0.12\times 0.14$ m $^2$ . Particle images recorded by the two cameras were processed using a final pass of $48\times 48$ pixel interrogation window and a 50 % overlap, resulting in a velocity vector spacing of about $0.8$ mm. The free-surface intersection with the SPIV plane was detected from laser-induced fluorescence (LIF) images recorded by a camera viewing the plane at an oblique angle from the air side. A small amount of rhodamine-6G was added to the water generating image contrast between the air and water regions, and the free surface was detected from the image intensity gradient. Further details can be found in Smeltzer et al. (2023).

Table 3. Wave quantities derived from measured values in table 2. Here, $\lambda _0=2\pi /k_0$ is the carrier wavelength, $c(k_0)=\sqrt {g/k_0}$ is the carrier-wave phase velocity, $u_s(0)=(ak_0)^2c(k_0)$ is the Stokes drift velocity of the carrier wave at the surface, $T_0=\lambda _0/c(k_0)$ is the intrinsic wave period, $\tau _0$ the intrinsic group width from (3.2) and $u_{{rf}}$ is the Eulerian return flow from (4.2).

For each wave group, SPIV/LIF measurements taken during three time intervals were used: well before the group arrived, and at the leading and trailing edges of the group envelope, referred to as intervals 1, 2 and 3, respectively, as shown in figure 2 b). We consider the difference in mean velocity between intervals 1 and 3 here, with interval 2 as a check to verify that the change is indeed due to waves. Values for all three intervals in all cases can be found in the Supplementary Materials available at https://doi.org/10.1017/jfm.2026.11163. The duration of the measurements for each interval, $T_{\textit{PIV}}$ , was $10$ s and PIV and LIF were sampled at $f_{\textit{ac}}=8$ Hz. After each group, residual waves from reflections were allowed to dissipate for approximately five minutes before the next wave group was generated. The above procedure was performed a total of $N_{\textit{ens}}=60$ times to produce ensemble statistics, except for the case 1.C.2 which was performed $20$ times. In case 1.B, only vertically orientated bars of the active grid were actuated. For case 1.A, the grid was stationary with the wings aligned with the flow in the position of least blockage (see also Appendix A). The experimental conditions are listed in table 1.

Figure 2. (a) Example surface elevation of a single wave group from Experiment 1 as a function of time $t$ , measured by a wave probe at the measurement location. (b) Example from Experiment 1 of an ensemble-average group surface elevation amplitude envelope as a function of time normalised by the measured group temporal width $\tau$ for case 1.D (see (3.1)). The time intervals for SPIV measurement (1–3) are shown with vertical dashed lines. (c) Surface elevation measurements of one ensemble from Experiment 3 (cases 3.A and 3.B), which shows the onset of a regular-wave train. The red box indicates the interval used for analysis.

The wave groups were generated with the wavemaker motion having carrier frequency $f_0 = 1.02$ Hz and a Gaussian amplitude envelope of the form

(2.1) \begin{equation} S(t) = S_0 \exp\left [-\frac {(t-T_{{wm}}/2)^2}{2\tau _{{wm}}^2}\right ]\!, \end{equation}

for $0\leqslant t\leqslant T_{{wm}}$ , where $S_0$ is the peak stroke, $\tau _{{wm}}= 6\,{\textrm{s}}$ is the group width in time as applied to the wavemaker signals and $T_{{wm}} =24\,{s}$ was the duration over which the wavemaker plungers were actuated. The surface elevation for one wave group measured at the SPIV measurement location is shown in figure 2(a).

2.1.2. Experiment 2: regular waves

Cases 2.A and 2.B in tables 1–3 are from Experiment 2. Regular waves were generated propagating upstream atop two different flows with comparable mean velocity but different levels of turbulence as controlled by the active grid. A planar PIV set-up with a light sheet, orientated in the streamwise–vertical ( $xz$ ) plane, measured the in-plane streamwise and vertical velocity components. The field of view was centred $L_{\textit{FOV}}=8.5$ m downstream of the surface plate’s trailing edge (see figure 1). Three cameras stacked vertically (from the top to bottom: two 16 megapixel cameras and one 5.5 megapixel camera) covered a field of view extending over the entire water depth ( $h = 0.80$ m) of the channel, roughly $21.7$ cm wide. The parameters for Experiment 2 are listed in table 1. For both flow cases, waves of two frequencies $0.94$ and $1.16$ Hz each with three different steepness values were generated. For each case, $2000$ PIV images were acquired at a sampling frequency of $0.86$ Hz; thus the total measurement time $T_{\textit{PIV}}$ was approximately $39$ min. Using a final pass with $48 \times 48$ pixel interrogation window and a 50 % overlap resulted in a velocity vector spacing of approximately $1.6$ mm. For both cases 2.A and 2.B, PIV measurements were also performed without wave generation to characterise ambient flow conditions, reported in table 2.

2.1.3. Experiment 3: onset of regular waves

Cases 3.A and 3.B in tables 1–3 are from Experiment 3. This set of measurements was taken at the highest acquisition frequency in our study ( $15$ Hz) during a $25$ s interval just after the first arrival of regular waves as indicated in figure 2(c). The water depth was $0.50$ m, and a slower mean flow of $U_0=0.19\, \textrm{m s}^{-1}$ was used compared with the other cases (see table 1). The two cases are separated by the different sequences for the active grid. Case 3.A is similar to case 1.B as only vertically orientated grid bars were rotated at a random, top-hat-distributed rotation speed. Similarly for case 3.B, only the vertical grid bars were actuated, made to flap over an angle of $\pm 60^\circ$ about the fully open position during random time intervals of a top-hat-distributed duration between $0.5$ and $1.0$ seconds. (A full list of experimental conditions is found in Appendix A.) No surface plate downstream of the active grid was used in this experiment, and the field of view in the streamwise–vertical plane (see figure 1 c) was centred a distance $L_{\textit{FOV}}=8.5$ m downstream of the grid.

The wavemaker was actuated at $f_0=1.4$ Hz. Here, two steepnesses for each flow case are investigated, yielding four different steepness–turbulence combinations. For each of the four combinations, $32$ ensembles are measured, where the wave probe measurements of one ensemble are shown in figure 2(c). In order to avoid wave breaking at the leading edge of the wave train, the wavemaker stroke is linearly ramped up to the set amplitude over a period of $6$ seconds.

For the PIV measurements, a single 25-megapixel camera was used to cover the entire water column. Using a $64\times 64$ pixel interrogation window and a 50 % overlap resulted in a final velocity vector spacing of approximately $3$ mm. The field of view is $0.45\, {\rm m}\times 0.5$ m. PIV images were acquired at $f_{\textit{ac}}=15$ Hz. Also here, PIV measurements were performed without wave generation to characterise the ambient flow conditions, reported in table 2. For these background measurements $2000$ snapshots were captured at $1.0$ Hz.

3. Experimental measurements

3.1. Flow and wave characteristics

The measured physical characteristics of mean flow, turbulence and waves are listed in table 2, and some derived quantities we will make use of are quoted in table 3 for convenience.

An assumption of deep water is well satisfied since in all our cases $k_0h\gtrsim 3.6$ hence $1\gt \sqrt {\tanh (k_0h)}\gt 0.999$ . Here, and henceforth, $U_0$ is the mean absolute surface speed in the absence of waves (i.e. the mean surface velocity is $\boldsymbol{U}_0=(-U_0,0,0)$ ) and in all our cases the velocity profile is sufficiently depth uniform that vertical shear affects wave dispersion negligibly (e.g. Ellingsen & Li 2017). The wavenumber $k_0$ and the wavemaker carrier frequency $f_0$ are related approximately by $2\pi f_0 = \sqrt {gk_0} {-} U_0k_0$ , or equivalently $f_0 = (c(k_0){-}U_0)/\lambda _0$ with $c(k_0)$ from (1.2) and the wavelength is $\lambda _0=2\pi /k_0$ . We will use the measured (as opposed to derived) value of $\lambda _0$ and the wavenumber $k_0$ accordingly. Values for these quantities are listed in tables 2 and 3, respectively. The time it takes for the waves to propagate from wavemaker to the edge of the surface plate is around $40$ s for Experiment 1, $25$ s for cases 2.A.1 and 2.B.1, $56$ s for cases 2.A.2 and 2.B.2 and $35$ s for Experiment 3.

The root-mean-square (r.m.s.) of the turbulent velocity fluctuation after subtracting the average is defined as $\boldsymbol{u}_{\textit{rms}}=[u_{\textit{rms}},v_{\textit{rms}}, w_{\textit{rms}}]$ representing the streamwise, spanwise, and vertical components, respectively. The turbulent kinetic energy is as $e=({1}/{2}) |\boldsymbol{u}_{\textit{rms}}|^2$ . In cases 2.A–3.B, only $u_{\textit{rms}}$ and $w_{\textit{rms}}$ are available, so we use instead $e=({1}/{2}) (u_{\textit{rms}}^2+2w_{\textit{rms}}^2)$ , an assumption of cross-plane isotropy, as is the standard convention for grid turbulence (e.g. Comte-Bellot & Corrsin 1966; Lavoie et al. 2005) and was observed for our facility by Jooss et al. (2021). The assumption is not quite satisfied for our flow as table 2 shows for Experiment 1. However, given the inevitable anisotropy of the flow as the free surface is approached (also with no waves) and consequent depth variation of $e$ closer to the surface than about one integral scale (the blockage layer, see, e.g. Aarnes et al. 2025), representing turbulent kinetic energy (TKE) by a single number is always a simplified picture. The precise value of $e$ does not significantly affect our conclusions, and given that no simple expression for $e$ in terms of $u_{\textit{rms}}$ and $w_{\textit{rms}}$ can reconcile all cases in Experiment 1, we choose to go with the conventional expression.

The mean velocity $U_0$ was calculated from averaging all streamwise velocities over the lower part of the field of view; averaging was performed over $-12$ cm $\lt z\lt -5$ cm for Experiment 1 and $-20$ cm $\lt z\lt -10$ cm for Experiments 2 and 3. Full profiles of the Eulerian-mean velocity profiles for all experimental cases may be found in the Supplementary Materials.

When comparing turbulent quantities for the various cases one should bear in mind that these are measured at a fixed position in space, whereas the turbulence becomes gradually weaker as it travels downstream because of dissipation. The change in Eulerian-mean current is an integrated effect of wave–turbulence interactions that occurred upstream, where the turbulence intensity is in general a little higher. This is true of all cases, but because of the lower mean velocity, the turbulence that reaches the field of view in cases 3.A and 3.B has decayed for longer. A detailed study of turbulent decay in our lab with similar flow conditions as Experiment 1 was reported by Jooss et al. (2021). They found that for all cases, once the turbulence is fully developed TKE decays approximately as $\sim 1/x$ downstream for all cases. Thus, although the absolute values of TKE and Reynolds stresses are different at different positions, their relative magnitudes are expected to be constant. We consider only trends with respect to TKE, and ratios between Reynolds stresses which it seems reasonable to assume are well captured.

In previous experiments by Nepf & Monismith (1991) and Klopman (1994, 1997), secondary motion in the form of a pair of streamwise rolls were measured when waves were propagated on a current, triggered by the CL2 mechanism. We could measure the cross-plane flow in Experiment 1, and while some weak in-plane mean flow was observed, of the order of $1$ cm s⁻¹ or less, the pattern was not consistent with a coherent pair of rolls. The resulting convection, in magnitude and direction, would be insufficient to convect boundary-layer turbulence into the near-surface regions where turbulence and waves interact to any noticeable degree. During the travel from grid to measurement position, the vortex force near the boundaries has far less time to accelerate a similar amount of water as in Klopman (1994) and would likely not have time to develop. It is not perhaps so surprising that secondary flows would be different given that geometrical size and aspect ratios are rather different, and the current’s travel time from grid to measurement point is less than half that in either of the mentioned flows even for our slowest flow in Experiment 3.

It was found by Smeltzer et al. (2023) that the turbulent integral length scale, characteristic of the largest turbulent structures, is important for wave–turbulence interactions, especially so for the scatting of waves by turbulence. An advantage of our active-grid turbulence generation is that the integral scale can be varied independently of TKE. The most pertinent integral scale is related to correlation of streamwise turbulent velocity for points separated in the streamwise direction, $L_x^x$ . Due to differences in the way the experimental data were acquired, no single method for estimating $L_x^x$ can be applied to all cases. Estimating the integral length scales from experimental data uniquely and quantitatively is notoriously difficult, and the various methods in common use produce quantitatively different results. Additional challenges pertain to active-grid turbulence due to the slow spatial decay of the autocorrelation function (Puga & LaRue 2017; Mora et al. 2019). In Experiment 1, the integral scale was estimated with a zero-crossing method as described by Mora & Obligado (2020), and we used the same method to calculate the integral scale for Experiment 3, as listed in table 2. Since the data in Experiment 2 have low time resolution the same method cannot be applied there, and using another, spatially based method would not give directly comparable numbers, so we provide no integral scale for Experiment 2. We note that integral scales are significantly shorter in Experiment 3 than in Experiment 1 at similar turbulence levels, which can be explained by the lower mean-flow velocity $U_0$ .

Several practical aspects in evaluation of the wave and flow characteristics reported in table 2 differed for the three experiments, and are described in further detail below.

3.1.1. Wave group measurements (Experiment 1)

The mean flow and turbulence statistics without waves presented in table 2 were evaluated over the first interval, where any influence from waves can be assumed to be negligible. The mean flow in the negative $x$ -direction had approximately constant (absolute) value $U_0$ , found from averaging as described above. A small spanwise and vertical mean velocity was found, everywhere below $1$ cm s⁻¹, mostly much less, which we attribute to the channel flow conditioning not being perfect (achieving a perfectly straight and uniform flow in a channel of our size is very challenging). Notably, unlike the streamwise velocity, the cross-plane velocities did not change after the passage of the wave group, which implies that there is negligible contribution secondary flow due to wave–current interactions near the sidewalls mentioned above.

The characteristic peak wave amplitude $a_{{p}}$ and the observed group temporal width $\tau$ were estimated from the ensemble-averaged amplitude envelope of the wave groups as measured by the probes near the SPIV/LIF laser sheet as seen in figure 2(b). The average envelope was fitted to a Gaussian function of the form

(3.1) \begin{equation} a(t) = a_{{p}} \exp\! \left [-\frac {(t-t_{{p}})^2}{2\tau ^2}\right ]\!, \end{equation}

with $t_{{p}}$ the temporal location of the group peak. The peak wave steepness is $k_0a_{{p}}$ .

We define an intrinsic temporal group width $\tau _0$ , listed in table 3 as

(3.2) \begin{equation} \tau _0 = \tau \frac {c_g(k_0)-U_0}{c_g(k_0)}, \end{equation}

with $c_g$ defined in (1.2). The intrinsic temporal group width is expressed in a reference frame without mean flow and reflects the time scale during which the ambient turbulence interacts with the wave groups.

3.1.2. Measurements with regular waves (Experiments 2 and 3)

The ambient flow statistics were evaluated from PIV measurements acquired without waves. Similarly to Experiment 1, the mean-flow profile varied only slightly across the measurement plane, and a representative absolute value $U_0$ is given in table 2.

The wave steepness $ak_0$ was calculated using the average wave amplitude from the wave probe measurements in the proximity of the PIV measurement location. The amplitude of each individual wave oscillation varied slightly during the experiments, especially in cases with the highest level of turbulence, likely due to wave–turbulence interactions (e.g. Smeltzer et al. 2023). The variation is not of central interest to the present study, and thus only the mean steepness value is reported.

The regular waves were present throughout the entire test section during the experiments. The grid-generated and measured turbulence thus interacted with the waves over a length $L_{\textit{FOV}}$ as given in §§ 2.1.2 and 2.1.3, with associated interaction time $T_{\textit{int}}=L_{\textit{FOV}}/U_0$ . The frequency at which turbulence encounters waves, i.e. as seen by the moving flow, equals the intrinsic wave frequency (in Hz)

(3.3) \begin{equation} f_{\textit{enc}} = f_0 + U_0/\lambda _0 = \sqrt {gk_0}/2\pi. \end{equation}

3.2. Measured change in Eulerian-mean velocity

The measured changes in Eulerian-mean velocities presented below are of the order of millimetres per second, yet, while this is the same order of magnitude as our PIV measurement accuracy, one should bear in mind that the variance of an average from hundreds of independent measurements is far lower than that of single measurements. It is crucial that a careful analysis of errors and statistical convergence be performed in order to establish confidence that our results are accurate and reliable. In Appendix B, we report results of these tests, supported by further data in the Supplementary Materials.

3.2.1. Velocity change after the passage of wave groups (Experiment 1)

We now consider the measured Eulerian-mean velocity change for Experiment 1, i.e. cases 1.A–1.D. Vertically sheared flows in our laboratory have been found to be stable and therefore the wave-modified Eulerian-mean velocity profile remains close to unchanged throughout interval 3 (also found to be true for strongly sheared currents, see § 4 of Pizzo et al. 2023). In figure 3(a) we show as an example the mean streamwise velocity profiles for case 1.D during the three measurement intervals, denoted $U_I(z)$ for intervals $I = \{1,2,3\}$ (see figure 2 b). As can be seen, the mean-flow speed had increased in the direction opposite to wave propagation after the passage of the wave groups. Similar plots for all cases are provided in the Supplementary Materials. Error bars are omitted from the figure for reasons of visibility, but discussed in Appendix B. The net Eulerian-mean-flow change $\Delta U=U_3 {-} U_1$ is shown in figure 3(b) for flow cases 1.A–1.D as expressed in the legend. For all cases, $\Delta U\lt 0$ near the surface, and decays to small absolute values at depths $|k_0z|\gtrsim 1$ .

Figure 3. Change in Eulerian-mean current due to the passage of a wave group. The waves travelled in the positive $x$ -direction, against the current. (a) An example of mean streamwise velocity depth profile before the arrival of a wave group, $U_1(z)$ , and after the group has passed, $U_3(z)$ , here for case 1.D; (b) mean streamwise velocity difference $\Delta U = U_3-U_1$ as a function of depth for the flow cases of Experiment 1. Error bars are omitted for visibility – see analysis in Appendix B; (c) the slope of $\Delta U(z)$ relative to the Stokes drift gradient (a prime denotes derivation with respect to $z$ ). Light smoothing (moving average with window size $8$ mm) was applied to the curves in panel (c) for better visibility.

Interestingly, $\Delta U$ clearly depends on the level of ambient turbulence; while the waves have approximately equal properties in cases 1.A, 1.B, 1.C.1 and 1.D, the velocity change is far higher at the highest turbulence level (case 1.D) than at the lowest level (case 1.A), with the intermediate cases 1.B and 1.C.1 in between. Moreover, comparing cases 1.C.1 and 1.C.2 where two wave steepnesses are tested on the same turbulent current indicates a positive correlation between $k_0a_{{p}}$ and current change $\Delta U$ . Put together, the results in figure 3(b) provide a strong indication that the change in current is due to an interaction between waves and turbulence.

We can exclude the possibility that the measured $\Delta U$ is due to the Eulerian return flow under groups of waves, as measured by van den Bremer et al. (2019). The return flow follows the group, i.e. it is very weak in intervals 2 and 3 which lie outside the main group itself. It is also near-uniform in depth, whereas the current measured in figure 3(a) is strongly depth dependent.

It is instructive to plot the ratio between ${\rm d}(\Delta U)/{\rm d} z$ and $-{\rm d} u_s/{\rm d} z$ as a function of $z$ , shown in figure 3(c), because it gives some indication of how far the turbulent flow has transitioned towards a new, quasi-equilibrium state. We will later present theory and evidence that at the end of the ‘spin-up’ period, this ratio should, in the final state, be approximately equal to $u_{\textit{rms}}^2/w_{\textit{rms}}^2$ nearest the surface. Since our bulk turbulence is slightly anisotropic (in the cases reported in Jooss et al. (2021), $1.2\lesssim u_{\textit{rms}}^2/w_{\textit{rms}}^2\lesssim 1.4$ ), a final current change $|\Delta U|\gtrsim |u_s|$ is expected nearest the surface. We shall later see that in the regular-wave cases where the equilibrium is likely reached, this relation holds well. Since none of the changes in currents in cases 1.A–1.D are close to reaching these values, it seems that the passing of the wave group has not led to a wave–turbulence interaction of sufficient duration for a final state to be reached, and the flow is still relatively early in the ‘spin-up’ stage.

There are at least two striking observations to make in figure 3(c). First, with the exception of the low-turbulence case, 1.A, the ratio between the slopes is close to constant with depth, which illustrates that $\Delta U \sim \exp (2k_0z)$ near the surface in these cases. The scaling is not perfect, particularly at the shallowest depths; this should not be surprising since the depth dependence of wave–turbulence interaction should scale not only with the wavelength, but also the turbulent integral length scale, which delimits the vertical extent of the topmost layer where the kinematic boundary condition at the surface begins to limit the vertical extent of turbulent eddies (the blocking effect, see, e.g. Teixeira & Belcher 2002). Second, while the value of $\Delta U(z)$ after the passing of a group depends strongly on the turbulence level and steepness, there is no such trend for the relative slope ( $-\Delta U^{\prime }(z)/u_s^{\prime }(z)$ ), case 1.A excepted.

In § 4.4 we will develop a RDT model describing the early onset of wave–turbulence interaction, which we can compare with the measurements in figure 3(b), given this evidence that the combined wave/turbulence flow is still far from fully developed in Experiment 1.

In conclusion, the evidence suggests that the change in Eulerian-mean current observed in our experiments is due to wave–turbulence interaction, and increases with increasing turbulence and increasing wave steepness.

3.3. Velocity change under regular waves (Experiments 2 and 3)

Cases 2.A–3.B all consider turbulence interacting with regular (i.e. continuous and periodic) waves. For cases 2.A.1–2.B.2, we evaluate the mean streamwise velocity in the presence of waves, and subtract off the mean velocity profile from the ambient flow case without waves. This velocity difference we define as $\Delta U(z)$ . Although the flow is wavy, the time series is long enough for the Eulerian-mean current, averaged over all $2000$ PIV images, to be well converged in the sense that further measurements would affect it insignificantly. See details in Appendix B.

The turbulent current is affected by the waves during the time it takes it to traverse the test section, a distance $L_{\textit{FOV}}$ as defined in §§ 2.1.2 and 2.1.3. The interaction time is thus $T_{\textit{int}}=L_{\textit{FOV}}/U_0$ , approximately $29$ s for Experiment 2 and $46$ s for Experiment 3, the latter having a slower flow speed. This corresponds to a number of wave cycles, $T_{\textit{int}}/T_0,$ between $36$ (cases 2.A and 3.A) and $83$ listed in table 4. Due to the slower flow speed. Cases 3.A and 3.B interact for considerably longer than 2.A and 2.B, both in terms of absolute time (in seconds) and in terms of number of intrinsic wave periods $T_0$ , i.e. $T_{\textit{int}}/T_0$ is the number of wave cycles that the turbulence has encountered.

Table 4. Derived wave quantities for use in RDT analysis. For cases 1.A–1.D, $T_{\textit{int}}=\sqrt {\pi }\tau _0$ is used with $\tau _0$ from (3.2), while for cases 2.A–3.B, $T_{\textit{int}}=L_{\textit{FOV}}/U_0$ ; $\beta _{{f}}(0)$ is found from (4.28) and (4.29) for groups and regular waves, respectively. Note that it is related to $T_{\textit{int}}$ via (4.33).

In figure 4 we show $\Delta U(z)$ for cases 2.A–3.B, where the different wave steepness values $ak_0$ are labelled in the legend. Several characteristic behaviours can be observed. First, all graphs have similar dependence on depth, with the highest absolute value nearest to the surface, decreasing down to a level where $k_0z$ is roughly in the range between $-1$ and $-2$ , then turning and slowly becoming more negative again. (The non-monotonicity (turning) would not be visible in figure 3, where measurements were only made for $k_0z\gt -1.2$ .) Indeed, some cases see the induced Eulerian-mean current take positive values over a certain depth range. (One may note that the current profile is similar in shape and magnitude to those of Rashidi et al. 1992). Second, there is again a clear tendency that higher steepness leads to a larger vertical variation of $\Delta U(z)$ . Finally, there appears to be a near-constant offset between graphs (an addition to the mean flow), which increases with steepness, and we will find in § 4.1 that it can be explained, at least in part, as the approximately depth-uniform Eulerian return flow $u_{{rf}}$ .

Figure 4. The wave-induced current $\Delta U$ under regular waves as a function of depth $k_0z$ for cases 2.A.1, 2.A.2, 2.B.1, 2.B.2, 3.A and 3.B in panels (a)–( f), respectively. For each case, different wave steepness values $ak_0$ are shown as indicated in the legend. The dashed lines are the theoretical Stokes drift profiles at the same location for each case, shown as $-u_s(z)$ , that is, with opposite sign to the Stokes drift. The filled circles at $k_0z=-4$ and $0$ indicate the theoretical value of the Eulerian return flow, $u_{{rf}}$ .

The reverse of the theoretical Stokes drift $-u_s(z)$ from (1.1) is plotted in all panels of figure 4 with dashed lines of the same colour for comparison, calculated for each case. For all cases we note that $|\Delta U|\lt u_s$ nearest the surface, but that the depth variation ${\rm d} |\Delta U|/{\rm d} z\sim {\rm d} u_s/ {\rm d} z$ , is higher in some cases, smaller in others. We will return to these points later when comparing the results with theory.

To further illustrate how the Eulerian-mean-current change depends on wave steepness, we plot the value of $\Delta U$ at a set reference depth as a function of steepness $ak_0$ in figure 5. The value at the shallowest depth available for all cases is used, i.e. $k_0z=-0.27$ . For cases 2.A.1 and 2.B.1, $\Delta U$ appears to scale as $(ak_0)^2$ , indicated by the dashed line. Cases 2.A.2 and 2.B.2 do not adhere to the scaling, possibly due to the high wave steepness values involved, so that interactions of order $(ak_0)^3$ and higher may become significant. A further possible explanation is the lack of scale separation: the wavelength is likely only slightly larger than the integral length scale in these two cases (based on comparison with cases 1.A–1.D in table 2, which have similar $U_0$ – unfortunately, a direct comparison of $L_x^x$ is not possible, as explained). The comparatively large integral length scale and high turbulence levels mean that stronger angular scattering of waves on turbulent velocity changes is to be expected (Villas Bôas & Young 2020; Smeltzer et al. 2023), a process which does not scale with steepness. Figure 5 should be interpreted only qualitatively, since the value of $\Delta U$ at a constant value of $k_0z$ is not entirely comparable between cases with different turbulence properties. As discussed in connection with figure 3(c) and at length in § 4.4, $\Delta U$ depends not only on $k_0$ but also on the turbulent integral length scale and anisotropy, which varies between cases.

Figure 5. The wave-induced current under regular waves at the ‘reference’ depth $k_0z=-0.27$ as a function of wave steepness for cases 2.A–3.B as indicated in the legend. The dashed line is proportional to $(ak_0)^2$ .

4. Theoretical considerations

The experimental results give reason to hypothesise that the changes in Eulerian-mean flow are a consequence of the encounter between waves and pre-existing turbulence. Before considering the turbulent current at all, however, we discuss the well-known Eulerian return flow which is present under wave groups also in quiescent water (§ 4.1). We thereafter consider the model situation in which irrotational waves appear in the presence of pre-existing turbulence and the two begin to interact. In § 4.2, the argument made by Pearson (2018) is revisited and adjusted to our case, that the combined flow will undergo a transition – a ‘spin-up’ as McWilliams et al. (1997) call it – to a new statistically steady state in which a new, depth-dependent Eulerian-mean current must be present. Next, a theory based on RDT is used to make more detailed predictions of the process during the `spin-up’ period, predicting the depth profile of the ‘wave-Reynolds’ stress which drives a mean current.

Throughout the theory sections, we assume that there is a negligible mean pressure gradient in the flow (apart from the local pressure gradient associated with a wave group) and that the mean flow after averaging over wave and turbulent motion (both time and ensemble average) lies in the streamwise–vertical plane, i.e. $\bar {\boldsymbol{u}} = (U_0,0,0)$ , i.e. $\bar {v}=\bar {w}=0$ , and all averaged quantities are presumed independent of the spanwise coordinate $y$ . At the time and length scales we consider, the Coriolis force is irrelevant.

We assume that a triple decomposition of the velocity field can be made according to

(4.1) \begin{equation} \boldsymbol{u}(\boldsymbol{x}, t) = \bar {\boldsymbol{u}}(\boldsymbol{x}) + \tilde {\boldsymbol{u}}(\boldsymbol{x}, t) + \boldsymbol{u}'(\boldsymbol{x}, t) \equiv \check {\boldsymbol{u}}(\boldsymbol{x}, t) + \tilde {\boldsymbol{u}}(\boldsymbol{x}, t), \end{equation}

where $\tilde {\boldsymbol{u}}$ and $\boldsymbol{u}'$ are the contributions from waves and turbulence, respectively, and the Eulerian-mean flow $\bar {\boldsymbol{u}}$ changes slowly compared with the period and wavelength of an individual wave. The wave-filtered (Eulerian, turbulent) current is denoted $\check {\boldsymbol{u}}=\bar {\boldsymbol{u}}+\boldsymbol{u}'$ . Performing the decomposition in practice is non-trivial when the waves are not monochromatic and/or there is no clear scale separation in time or space. We tried and evaluated two different methods for separations – the widely used ‘phase-conditioned averaging’ (PhCA) and proper orthogonal decomposition (POD) – and found the latter to perform better. A detailed comparison is found in Appendix C.

4.1. Irrotational Eulerian return flow

As is well known, the passage of a deep-water wave group does not entail a net mass flux, because the Stokes drift is cancelled in a depth-integrated sense by an Eulerian current in the opposite direction (Longuet-Higgins & Stewart 1962) known as the return flow. The current has been observed in laboratory studies, e.g. van den Bremer et al. (2019).

The return flow is found beneath a passing wave group, following the group and tending rapidly to zero ahead of and behind the group. For this reason, velocity measurements in cases 1.A–1.D which are taken before and after the passage of the wave group, are expected to see negligible influence from such flow. Cases 2.A–3.B, on the other hand, occur under regular waves – in practice a very long wave group. Since the length of the ‘group’ is far longer than the depth of the channel, the current will be approximately depth-uniform and, in the system following the mean current without waves, satisfying Lagrangian mass conservation, i.e.

(4.2) \begin{equation} u_{\textit{rf}} =-\frac1h \int_{-h}^0 u_s(z;h){\rm d} z \approx -\frac1h \int_{-\infty}^0 u_s(z){\rm d} z = - \frac{k_0a^2}{2h}c(k_0), \end{equation}

where we have assumed infinite water depth (i.e. $k_0h\gg 1$ ) for the Stokes drift profile, as we argued above, and inserted $u_s$ from (1.1). The relation is valid when the depth is greater than about half a wavelength but much smaller than the group length ( $\sim L_{\textit{FOV}}$ in our experiment), both of which are well satisfied in our experiments. We note that (4.2) is only approximate, altered somewhat by the effect of the channel’s bottom and sidewalls (see e.g. van den Bremer & Breivik (2017) for a discussion).

4.2. Statistical theory of Eulerian flow generation

Pearson (2018) argued that, when waves appear in a turbulent field, the combined flow will undergo a transition to a new statistically steady state. We will tailor his argument to our current application and discuss consequences in the context of our experimental observations. Momentum conservation, expressed through the time-averaged Navier–Stokes equations or the Craik–Leibovich equations as appropriate, implies that an Eulerian-mean current will manifest during the ‘spin-up’ period, before a quasi-steady state is reached (in our case slowly decaying due to dissipation).

4.2.1. Turbulent statistics in the transition period

While quiescent-water irrotational waves and turbulence may each be steady in isolation (except for their decay due to dissipation), together they form a system in disequilibrium which will go through a transient change to a new quasi-steady state. This is similar to the model of ‘spin-up from rest’ employed in theory for Langmuir turbulence (e.g. McWilliams et al. 1997). Revisiting and, to some extent, adapting the arguments of Pearson (2018) sheds light both on the early-stage ‘spin-up’ stage and the final situation. We first assume that a triple decomposition of the turbulent and wavy flow can be performed according to (4.1) (by no means a trivial point, an intricacy we shall return to later) and that all averaged quantities vary slowly as a function of $x$ over a wavelength so that their derivative can be neglected to leading order. Moreover, we assume that any changes in the mean flow develop slowly compared with a wave period. This is reasonable because we have a strong indication that the interaction during passage of the wave groups in Experiment 1, which effectively lasts for approximately $T_{\textit{int}}\sim 7$ wave periods, is far from enough to reach the quasi-equilibrium state: the ratio $|\Delta U'(z)/u_s'(z)|$ only reaches approximately $0.5$ , as figure 3(c) shows, whereas we shall see in the next section (in turn backed by experiments) that in the final quasi-equilibrium state this fraction is $\gtrsim 1$ . The wave field, and consequently the Stokes drift, changes slowly, so $\partial _t\boldsymbol{u}_s$ is negligible. The flow is assumed to be unforced, without, e.g. wind stress, and there is no influence from buoyancy or the Coriolis force. The details of the dissipation are not important to the argument, so we shall not consider them beyond assuming that viscous decay is slow compared with the mean-flow effects of wave–turbulence interactions (which is reasonable in light of the observations by Jooss et al. 2021), so that a quasi-steady state can be reached. The final assumption is that $\bar {\boldsymbol{u}}$ and $\boldsymbol{u}_s$ are both oriented along the $x$ -axis and vary slowly except with coordinate $z$ .

The turbulent flow is observed to develop slowly at the relevant scales, and can be assumed to be little affected by a Lagrangian average over a wave period. Wave averaging the equations of motion yields the so-called Craik–Leibovich equation which under the above assumptions may be written in the form (e.g. Suzuki & Fox-Kemper 2016)

(4.3) \begin{equation} \partial _t \check {\boldsymbol{u}} + (\check {\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{\nabla })\check {\boldsymbol{u}} = \boldsymbol{u}_s\times \check {\boldsymbol{\omega }} -\boldsymbol{\nabla }\left(\check {p} + { \frac 12} u_s^2 + \boldsymbol{u}_s\boldsymbol{\cdot }\check {\boldsymbol{u}} \right) + \text{diffusion,} \end{equation}

where $\check {\boldsymbol{\omega }}=\boldsymbol{\nabla }\times \check {\boldsymbol{u}} =(0,\partial _z\bar {u},0) + \boldsymbol{\nabla }\times \boldsymbol{u}'$ denotes the (Eulerian) wave-averaged vorticity field. Taking the $x$ -component and ignoring diffusion gives

(4.4) \begin{equation} \partial _t \check {u} + \bar {u}\partial _xu'+ (\boldsymbol{u}'\boldsymbol{\cdot }\boldsymbol{\nabla })u'+ w'\partial _z\bar {u} \approx -\partial _x (p' + u_s u' ), \end{equation}

where we have used $\bar {w}=0$ . Performing a Reynolds average over turbulent motions (denoted with an overbar), while noting that $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}'=0$ and that averaged quantities are independent of $x$ , eventually yields

(4.5) \begin{equation} \partial _t \bar {u} \approx - \partial _z \overline {u'w'}. \end{equation}

Equation (4.5) has two important consequences. First, that an inhomogeneous field of turbulence will drive a mean current for as long as a vertical gradient of the shear stress exists (and dominates over viscous forces). Second, when a steady state has been reached, i.e. $\partial _t\bar {u}=0$ , then $\overline {u'w'}$ is approximately constant with respect to depth (see also Pearson 2018, for more discussion), or, more precisely, its vertical gradient is balanced by viscous diffusion. Physically, $\overline {u'w'}\neq 0$ represents a vertical redistribution of streamwise momentum, allowing $\boldsymbol{u}=\boldsymbol{u}(z)$ to transition from its initial value to a presumed final steady-state value, so when the transition period is ended, this shear stress should thus be small compared with TKE. Equation (4.5) paves the way for further analysis of the ‘spin-up’ of the Eulerian-mean current, because the right-hand side can be related to the underlying physical process using a model based on RDT.

4.3. Statistics in the quasi-equilibrium state

Following Harcourt (2013) and Pearson (2018), it can be readily argued that the development of a mean flow via the turbulent shear stress $\overline {u'w'}$ is due to interaction between pre-existing turbulence and Stokes drift.

Multiplying the $x$ -component of (4.3) by $w'$ and the $z$ -component by $u'$ , adding them together and averaging, yields

(4.6) \begin{equation} \frac {\partial \overline {u'w'}}{\partial t}+\overline {w'w'} \frac {{\rm d} \bar {u}}{{\rm d} z} =- \overline {u'u'}\frac {{\rm d} u_s}{{\rm d} z} - \frac {{\rm d} }{{\rm d} z}\overline {u'w'w'} + \textrm{diffusion}, \end{equation}

where we have employed the chain rule and the fact that $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}'=0$ . The diffusion term contains viscous dissipation and turbulent pressure fluctuations, both of which may reasonably be assumed to be small in an oceanographic setting (Harcourt 2015; Pearson 2018) (unfortunately we cannot directly ascertain how accurate this assumption is in the present context; see also Pearson, Grant & Polton 2019). We should bear in mind that this assumption is questionable for our experiment, where waves are shorter and turbulence levels considerably higher than in typical field settings. This is especially true for cases 2.A.2 and 2.B.2. For instance, angular diffusion of the waves can be highly significant when the waves are not long compared with the turbulent integral scale as found by Smeltzer et al. (2023) which would cause a spatially variable Stokes drift. The approximation is more justified the larger the Stokes drift ( $\sim (ak_0)^2c(k_0)$ ) is compared with turbulent velocities ( $\sim u_{\textit{rms}}$ ). In cases 2.A.2 and 2.B.2 the turbulence is strong and the waves short, and turbulent diffusion and pressure correlation terms may not be negligible compared with the Stokes drift contributions even at surface level. Unfortunately, we are not in a position to quantify the diffusive terms in (4.6) from our data.

We now assume that a quasi-steady state has been reached, so that the explicit time derivative of $\overline {u'w'}$ is negligible. The term $\overline {u'w'w'}$ represents net vertical transport, and is expected to become very small in a quasi-steady state. Assuming the diffusion is small, it follows that

(4.7) \begin{equation} \overline {w'w'} \frac {{\rm d} \bar {u}}{{\rm d} z} \approx - \overline {u'u'}\frac {{\rm d} u_s}{{\rm d} z}, \end{equation}

which is to say that an Eulerian-mean current must have been created with the opposite sign compared with Stokes drift.

Note carefully that the relation (4.7) will only hold as long as the right-hand side term is large enough to dominate over the terms in (4.6) we have neglected. Since the Stokes drift decreases exponentially with depth, our assumptions will surely be highly suspect deeper than $k_0z \sim -2$ , possibly above if the turbulence is strong.

In a situation with two horizontal directions (such as in an ocean wave model), the above argument easily generalises to

(4.8) \begin{equation} \overline {w'w'} \frac {{\rm d} \bar {\boldsymbol{u}}_h}{{\rm d} z} \approx - \overline {\bar {\boldsymbol{u}}_h^{\prime}\boldsymbol{\cdot }\bar {\boldsymbol{u}}_h^{\prime}}\frac {{\rm d} \boldsymbol{u}_s}{{\rm d} z} ,\end{equation}

with $\boldsymbol{u}_s$ and $\bar {\boldsymbol{u}}_h$ lying in the horizontal plane. We emphasise that if a constant, depth-uniform current is initially present (as in our experiment), the current $\bar {\boldsymbol{u}}$ in (4.7) and (4.8) is the change in current due to wave–current interaction, or alternatively the current measured in the reference system following the original Eulerian-mean flow.

Beneath a free surface, turbulence is not isotropic, yet $\overline {u'u'}$ and $\overline {w'w'}$ can be expected to be of the same order of magnitude except very near the mean surface level where $\overline {w'w'}$ tends to zero. It is worth noticing that (4.7) is not a perfect cancellation between Eulerian flow and Stokes drift as suggested by the experimental (re)analysis of Monismith et al. (2007), but when turbulence is close to isotropic, the remaining Lagrangian-mean current could be difficult to distinguish from zero in an experiment. On the basis of available evidence it seems a reasonable conjecture that the wave-turbulence-generated anti-Stokes flow explains, at least partly, this surprising cancellation.

A key aspect to notice in (4.7) for modelling purposes is that the degree of cancellation of the mean Lagrangian current does not depend on the overall turbulence level, only on the anisotropy of turbulent fluctuations within the near-surface layer of thickness $\sim 1/k_0$ where Stokes drift is non-negligible. On the other hand, the rate of growth of the Eulerian current is proportional to $\partial _z\overline {u'w'}$ and is higher the more intense the pre-existing turbulence is. In an ocean setting some level of turbulence is nearly always present, so the partial cancellation of Stokes drift according to (4.7) is to be expected.

4.3.1. Comparison with experiment

In order to test (4.7), it is necessary to extract the variances $\overline {u'u'}$ and $\overline {w'w'}$ from the wave-turbulence PIV data. Accomplishing this with the best possible accuracy is paramount in order to estimate turbulent second-order moments, because the wave velocities far exceed those of the turbulence near the surface, so even a small percentage of the wave motion erroneously identified as turbulence could lead to considerable overestimates of $\overline {u'u'}$ and $\overline {u'w'}$ (less so $\overline {u'w'}$ due to the phase difference between $\tilde {u}$ and $\tilde {w}$ ).

In cases 3.A and 3.B one could employ the now-standard technique of PhCA (Umeyama 2005; Buckley & Veron 2017) to average out the wave motion, because the wave phase at each point in time and space can be estimated. For cases 2.A and 2.B the task is more difficult because measurements are made with frequency lower than the wave’s, and no phase information is available (Experiment 2, when performed, was not planned with this task in mind). Estimates of $\overline {u'u'}$ and $\overline {w'w'}$ can be made in all cases regardless of acquisition rate using the method of POD which decomposes the measured velocity field into spatial modes ordered according to their ‘energy’ content. A more detailed discussion of the performance of this method and validation may be found in Appendix C. In cases 3.A and 3.B direct comparison with PhCA can be made; for the purpose of comparing with (4.7), the two methods give very similar, but not quite identical results.

Although PhCA is in common use while POD has not been employed for this purpose previously to our knowledge, there are fundamental reasons why a data-driven method such as POD, in general, is more appropriate. In particular, PhCA relies on assumptions which typically, including in our case, are not fully met. A further discussion of this point is found in Appendix C, where we argue that POD is the preferred method of decomposition also in cases 3.A and 3.B.

The two sides of (4.7) are compared in figure 6 for cases 2.A to 3.B, the cases whose measurements were taken in the streamwise-vertical plane so that waves and turbulence can be separated. The POD procedure was employed for triple decomposition to evaluate turbulent second moments. The dashed lines show the vertical derivative of the Stokes drift, ${\rm d} u_s/{\rm d} z=2\sqrt {gk_0}(ak_0)^2\exp (2k_0z)$ , while the solid line is the measured value of $-(\overline {w'w'}/\overline {u'u'}){\rm d}\bar {u}/{\rm d} z$ , which is hypothesised to equal the Stokes drift gradient in the quasi-steady state. For completeness, the Reynolds stresses themselves are plotted in Appendix D.

Figure 6. Test of (4.7) for cases 2.A to 3.B. The theoretical Stokes drift gradient ${\rm d} u_s/{\rm d} z$ (the derivative of (1.1)) is shown as dashed lines of corresponding colour. Separation of turbulence from waves was performed with POD, discussed in Appendix C.

Panels (e) and ( f) of figure 6 are from Experiment 3 where the separation of waves and turbulence could be performed and validated with confidence, as discussed in Appendix C, but for cases 2.A and 2.B in panels (a)–(d) the uncertainly is difficult to quantify. For this reason, considerable caution should be exercised when interpreting the results in figure 6(a)–(d), particularly nearest the surface where the wave motion is most energetic, so that even a small percentage of wave motion remaining in the calculation of $\overline {u'u'}, \overline {w'w'}$ could cause a significant overestimation. We have decided to include the figures nevertheless, to illustrate the overall similarity of trend.

Given the highly irregular and unsteady nature of our strongly turbulent flow, the agreement is striking in many of the cases. As argued above, the assumptions behind (4.7) are expected to be reasonable when Stokes drift is sufficiently high compared with turbulent velocities, i.e. $k_0z \gtrsim -2$ and sufficiently large $ak_0$ , and subject to corrections if diffusion levels are high. The relatively poor fit for cases 2.A.2 and 2.B.2 is therefore not very surprising for the latter reason. The very low values of $ak_0$ in cases 2.A.1.1 and 2.A.2.1 is a likely reason why (4.7) is unlikely to be reasonable for these cases – again the Stokes drift term in (4.19) may not dominate over terms which are neglected. Indeed, the good performance in case 3.A.1 is more striking, but tallies with the fact that TKE is only half that of case 2.B and less than a quarter of that of case 2.A, with correspondingly less diffusion.

Some further words of caution are warranted. As discussed previously, turbulence was measured only at one position near the downstream end of the channel, while the wave–current interaction occurred upstream. We therefore cannot guarantee that the ratio $\overline {u'u'}/\overline {w'w'}$ has been the same throughout the time of interaction – indeed the theory of Teixeira & Belcher (2002) indicates this ratio has increased in magnitude near the surface during the spin-up period due to wave–current interaction. We hesitate to offer an explanation for the surprising fact that the experimental curves in figure 6 appear to curve back near the surface, corresponding to a similar behaviour of the fraction $\overline {u'u'}/\overline {w'w'}$ (see also figure 12 i), yet we remind the reader of the difficulty of separating waves and turbulence near the surface for these cases; we cannot rule out that the effect is spurious.

4.4. A model based on RDT

In § 3.2.1, we measured the change in Eulerian current after the passage of a wave group. Unlike for the regular-wave cases, the shear of the resulting current is considerably smaller than the gradient of the Stokes drift, so the two sides of (4.7) are highly dissimilar, which signifies that the influence of the group of waves has not brought the flow near the quasi-steady state. Since evidence suggests that the flow is still ‘spinning up’ when the group has passed, a theory which describes the development soon after the onset of waves is expected to describe the physical process in further detail.

The early onset of the (change in the) Eulerian current is governed approximately by (4.5), and we will use RDT to study how its right-hand side – that is, the Reynolds stress $\overline {u'w'}$ – depends on depth $z$ and time while changes are still moderate. The theoretical predictions will be compared with the observations made for cases 1.A–1.D, shown in figure 3(b). Teixeira & Belcher (2002) showed using RDT that a shear stress that varies in the vertical is generated by the passage of a progressive monochromatic surface wave over isotropic turbulence. We will show next that something similar occurs for a finite wave group, such as considered in the experiments of the present study. Since the effect we consider is inviscid in nature, we will neglect viscosity in the following.

The RDT is a theory first proposed by Batchelor & Proudman (1954) where the straining of turbulence due to the distortions from the surrounding flow is assumed to dominate over that due to turbulence acting on itself, so that the term $(\boldsymbol{u}'\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}'$ is negligible in the Navier–Stokes equation, yielding a linearised theory. The leading cause of change for the turbulence is presumed to be the mean Lagrangian motion of turbulence-containing parcels of fluid due to the larger-scale surrounding motion. The RDT approach is formally valid whenever the distortions applied to turbulence are sudden. This amounts to assuming that the distortions (for example, mean-flow gradients) are applied over a time shorter than an eddy turn-over time. The spectral formulation of RDT also requires that there is a spatial scale separation between the turbulence and the mean flow. In the present wave-associated mean flow, this corresponds to $\lambda \gg L$ where $\lambda$ is the wavelength of the waves and $L$ is the integral length scale of the turbulence (cf. Teixeira & Belcher 2002). Both of these criteria are reasonably well satisfied in the Experiments 1.A–1.D where the relevant turbulent length scale, $L_x^x$ is at most half a wavelength. Visual inspection of the velocity field clearly shows that the typical coherent eddies are much smaller than $\lambda _0$ . The values of $L_x^x$ and details of its calculation with the method of Mora & Obligado (2020) were reported in Smeltzer et al. (2023), and are quoted in table 2. In practice, RDT is known to provide useful results even when these conditions are not strictly fulfilled (e.g. Hunt & Carruthers 1990; Mann 1994; Cambon & Scott 1999).

In the present case, the effect that is essential in order to explain the generation of a Eulerian-mean current is not related to the individual wave oscillations but rather to the systematic tilting and stretching of the vorticity of the turbulence by the Stokes drift of the wave. Therefore, we adopt a linearised version of the Craik–Leibovich equation, (4.3), whose key term is the ‘vortex force’ $\boldsymbol{u}_s\times \check {\boldsymbol{\omega }}$ . Applying a simplified phase-averaged formalism is a good approximation because individual wave cycles have little effect on the net stretching and tilting of vortices, and it is the integrated effect of the whole wave group which matters; this is illustrated, e.g. in figure 3 of Teixeira & Belcher (2010) (regular waves are considered there, but the arguments hold equally for the relatively long wave groups we consider here). Only a small uncertainty is thereby introduced which is far smaller than the quantitative effects of the simplifications underlying RDT itself. The corresponding vorticity equation, which will be used in RDT, is

(4.9) \begin{equation} \frac {\partial \boldsymbol{\omega }'}{\partial t} + (\boldsymbol{u}_s\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\omega }' = (\boldsymbol{\omega }'\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_s. \end{equation}

In (4.9), $\boldsymbol{\omega }' = \boldsymbol{\nabla }\times \boldsymbol{u}'$ is the turbulent vorticity. As before, $\boldsymbol{u}_s$ is orientated along the $x$ axis and is independent of $x,y$ and $t$ . The mean background current is assumed to be initially depth uniform; hence, it has only trivial effect on the flow system, and we set it to zero by a change of coordinate system. It should be noted that these equations represent the effect of the Stokes drift of an irrotational surface wave; the rotational correction to the wave motion due to the small change in Eulerian current (see Ellingsen 2016) is neglected.

Let us first consider the turbulence away from the air–water interface, which to a first approximation can be considered homogeneous and isotropic, being denoted by the superscript $(H)$ . In this section, it is useful to let the subscript $i=1,2,3$ denote the component of a vector or tensor along directions $x,y$ and $z$ , respectively. The turbulent, homogeneous (or bulk) vorticity $\boldsymbol{\omega }'^{(H)}$ is expressed as a three-dimensional Fourier integral, as

(4.10) \begin{equation} \omega '^{(H)}_i(\boldsymbol{x},t)=\iiint \hat {\omega }^{(H)}_i(\boldsymbol{\kappa },t) \, {e}^{{i} \boldsymbol{\kappa } \boldsymbol{\cdot }\boldsymbol{x}} \, {\rm d} \kappa _1 \, {\rm d} \kappa _2 \, {\rm d} \kappa _3, \end{equation}

where the hat denotes a Fourier transformed turbulent perturbation quantity and $\boldsymbol{\kappa }=(\kappa _1,\kappa _2,\kappa _3)$ is the wavenumber vector.

Two equations result from (4.9) together with (4.10)

(4.11) \begin{align} \frac {{\rm d} \hat {\omega }^{(H)}_i}{{\rm d} t}&= \frac {\partial u_{si}}{\partial x_j} \hat {\omega }^{(H)}_j , \end{align}

(4.12) \begin{align} \frac {{\rm d} \kappa _i}{{\rm d} t} &= - \frac {\partial u_{sj}}{\partial x_i}\kappa _j. \end{align}

Since $u_{si} = u_s \delta _{1i}$ , (4.11)–(4.12) reduce to

(4.13a) \begin{align} \frac {{\rm d} \hat {\omega }^{(H)}_1}{{\rm d} t} &= \frac {{\rm d} u_s}{{\rm d} z}\hat {\omega }^{(H)}_3, \quad \frac {{\rm d} \hat {\omega }^{(H)}_2}{{\rm d} t}=0,\quad \frac {{\rm d} \hat {\omega }^{(H)}_3}{{\rm d} t}=0, \end{align}

(4.13b) \begin{align} \frac {{\rm d} \kappa _1}{{\rm d} t}&=0, \quad \frac {{\rm d} \kappa _2}{{\rm d} t}=0, \quad \frac {{\rm d} \kappa _3}{{\rm d} t}=-\frac {{\rm d} u_s}{{\rm d} z}\kappa _1. \end{align}

These equations can be integrated in time to yield

(4.14a) \begin{align} \hat {\omega }_1^{(H)}(t)&=\hat {\omega }_{10}^{(H)}+\hat {\omega }_{30}^{(H)} \int _0^t \frac {{\rm d} u_s}{{\rm d} z}{\rm d} t', \quad \hat {\omega }_2^{(H)}(t)=\hat {\omega }_{20}^{(H)}, \quad \hat {\omega }_3^{(H)}(t)=\hat {\omega }_{30}^{(H)}, \end{align}

(4.14b) \begin{align} \kappa _1(t)&=\kappa _{10}, \quad \kappa _2(t)=\kappa _{20}, \quad \kappa _3(t)=\kappa _{30}-\kappa _{10} \int _0^t \frac {{\rm d} u_s}{{\rm d} z}{\rm d} t', \end{align}

where the subscript ‘ $0$ ’ applied to a variable denotes its value at the initial time, before the turbulence is distorted. It is convenient to define

(4.15) \begin{equation} \beta = \int _0^t \frac {{\rm d} u_s}{{\rm d} z}{\rm d} t' = \frac {{\rm d} \Delta x}{{\rm d} z}, \end{equation}

where $\Delta x(z,t)$ is the total fluid parcel displacement in the $x$ -direction associated with the Stokes drift.

In order to calculate statistics of the turbulent velocity, it is necessary to relate the turbulent velocity fluctuations before and after distortion. Continuity of $\boldsymbol{u}'$ implies

(4.16) \begin{equation} \boldsymbol{\nabla} ^2 \boldsymbol{u}' = - \boldsymbol{\nabla } \times \boldsymbol{\omega }' ,\end{equation}

which in spectral space becomes

(4.17) \begin{equation} \hat {u}_i^{(H)} = \textrm {i} \varepsilon _{\textit{ijl}} \frac {\kappa _j}{\kappa ^2} \hat {\omega }_l^{(H)}, \end{equation}

where $\varepsilon _{\textit{ijl}}$ is the Levi-Civita permutation symbol and $\kappa = |\boldsymbol{\kappa }|$ .

We wish to relate the turbulent velocity after distortion by the Stokes drift to the velocity of the initial undistorted turbulence. The Fourier transform of the velocity of the homogeneous turbulence after distortion by the Stokes drift can be related to the corresponding Fourier transform of the vorticity using (4.17). In turn, the vorticity of the homogeneous turbulence after distortion can be related to the initial turbulent vorticity using (4.14a )–(4.14b ). Finally, the initial turbulent vorticity can be related to the initial turbulent velocity through

(4.18) \begin{equation} \hat {\omega }_{i0}^{(H)} = \textrm {i} \varepsilon _{\textit{ijl}} \kappa _{j0} \hat {u}_{l0}^{(H)}, \end{equation}

which comes from the definition of vorticity. From all of the above one finds that the distorted velocity can finally be expressed as

(4.19a) \begin{align} &\hat {u}_1^{(H)}(t)=\left ( 1 + \frac {\kappa _1 \kappa _3 \beta }{\kappa ^2} \right ) \hat {u}_{10}^{(H)}+ \frac {\kappa _1^2 \beta }{\kappa ^2} \hat {u}_{30}^{(H)}, \end{align}

(4.19b) \begin{align} &\hat {u}_3^{(H)}(t) = \left ( \frac {\kappa _0^2}{\kappa ^2} - \frac {\kappa _1 \kappa _{30} \beta }{\kappa ^2} \right ) \hat {u}_{30}^{(H)} - \frac {\kappa _{12}^2 \beta }{\kappa ^2} \hat {u}_{10}^{(H)}, \end{align}

where $\kappa _{12}=(\kappa _1^2+\kappa _2^2)^{1/2}$ , $\boldsymbol{\kappa }_0=(\kappa _{10},\kappa _{20},\kappa _{30})$ and $\kappa _0=|\boldsymbol{\kappa }_0|$ , and we only focus on $\hat {u}_1$ and $\hat {u}_3$ because these are the velocity components necessary to calculate the shear stress.

The previous calculations only apply to the turbulence far away from the air–water interface, at least a distance $\sim L$ where $L$ is a characteristic integral length scale (but affected by the Stokes drift, because of the scale separation $L \ll \lambda$ ). In order to take into account the blocking effect of the air–water interface, where we assume that the interface affects the turbulence essentially as a frictionless wall for depths of $O(L)$ (because the air–water interface has a large density contrast), this effect can be accounted for by adding an irrotational correction to the homogeneous turbulent velocity field, as done before by Hunt & Graham (1978) and Teixeira & Belcher (2002). Note that if the waves that generate the Stokes drift are irrotational (which is true to a good degree of approximation), then this irrotational correction remains irrotational. Hence, the turbulent velocity components affected both by the Stokes drift and by blocking can be expressed as two-dimensional Fourier integrals along the horizontal directions

(4.20) \begin{equation} u'_i(\boldsymbol{x},t) = \iint \hat {u}_i(\kappa _1,\kappa _2,z) e^{i(\kappa _1 x+ \kappa _2 y)} {\rm d}\kappa _1 \, {\rm d}\kappa _2, \end{equation}

because the inhomogeneity imposed by blocking does not allow Fourier transformation in the vertical direction.

Based on Hunt & Graham (1978) and Teixeira & Belcher (2002), the Fourier transforms of $u'_1$ and $u'_3$ can be written

(4.21a) \begin{align} &\hat {u}_1= \int \left ( \hat {u}_1^{(H)} e^{i \kappa _3 z} - i \frac {\kappa _1}{\kappa _{12}} e^{\kappa _{12} z} \hat {u}_3^{(H)} \right ) {\rm d}\kappa _3, \end{align}

(4.21b) \begin{align} &\hat {u}_3 = \int \hat {u}_3^{(H)} \big ( e^{i \kappa _3 z} - e^{\kappa _{12} z} \big ) {\rm d}\kappa _3. \end{align}

For a flow associated with a surface wave, the blocking condition actually applies perpendicularly to the wavy air–water interface, which when adopting a model accounting for the individual wave cycles requires the use of a curvilinear coordinate system, as in Teixeira & Belcher (2002). However, since our present model only includes the Stokes drift effect, no explicit wavy deformations of the interface are accounted for, and the surface is assumed to be in its average state, i.e. flat and horizontal at $z=0$ . The solutions (4.21) take this into account.

We wish to evaluate $\overline {u'w'} = \overline {u'_1 u'_3}$ after some time of deformation. This can be done by using (4.14b), (4.19), (4.20), and also noting that, by definition

(4.22) \begin{equation} \overline {\hat {u}_{i0}^{(H)}(\boldsymbol{\kappa }_0)^* \hat {u}_{j0}^{(H)}(\boldsymbol{\kappa }_0')} = \varPhi _{ij}^{(H)}(\boldsymbol{\kappa }_0) \delta (\boldsymbol{\kappa }_0 - \boldsymbol{\kappa }_0'), \end{equation}

where an asterisk denotes the complex conjugate, $\varPhi _{ij}^{(H)} (\boldsymbol{\kappa }_0)$ is the three-dimensional spectrum of the initial homogeneous and isotropic turbulence and $\delta$ is the three-dimensional Dirac delta, to show that $\overline {u'_1 u'_3}$ is given by

(4.23) \begin{equation} \overline {u'_1 u'_3} = \iiint \Big ( M_{13} \varPhi _{13}^{(H)}+M_{11} \varPhi _{11}^{(H)} + M_{33} \varPhi _{33}^{(H)} \Big ) {\rm d}\kappa _1 \, {\rm d}\kappa _2 \, {\rm d}\kappa _{30}, \end{equation}

where

(4.24a) \begin{align} M_{13}=&\left ( 1+ \frac {\kappa _1 \kappa _3 \beta }{\kappa ^2} \right ) \left ( \frac {\kappa _0^2}{\kappa ^2} - \frac {\kappa _1 \kappa _{30} \beta }{\kappa ^2} \right ) - \frac {\kappa _1^2 \kappa _{12}^2 \beta ^2}{\kappa ^4} \notag \\ &-\left [ \left ( 1 + \frac {\kappa _1 \kappa _3 \beta }{\kappa ^2} \right ) \left ( \frac {\kappa _0^2}{\kappa ^2} - \frac {\kappa _1 \kappa _{30} \beta }{\kappa ^2} \right ) - \frac {\kappa _1^2 \kappa _{12}^2 \beta ^2}{\kappa ^4} \right ] e^{\kappa _{12} z} \cos ( \kappa _3 z) \nonumber \\ &+ 2 \frac {\kappa _1 \kappa _{12} \beta }{\kappa ^2} \left ( \frac {\kappa _0^2}{\kappa ^2} - \frac {\kappa _1 \kappa _{30} \beta }{\kappa ^2} \right ) e^{\kappa _{12} z} \sin (\kappa _3 z), \end{align}

(4.24b) \begin{align} M_{11}=&-\frac {\kappa _{12}^2 \beta }{\kappa ^2} \left [ 1 + \frac {\kappa _1 \kappa _3 \beta }{\kappa ^2} - \left ( 1 + \frac {\kappa _1 \kappa _3 \beta }{\kappa ^2} \right ) e^{\kappa _{12} z} \cos (\kappa _3 z) +\frac {\kappa _1 \kappa _{12} \beta } {\kappa ^2} e^{\kappa _{12} z} \sin (\kappa _3 z) \right ], \end{align}

(4.24c) \begin{align} M_{33}=&\left ( \frac {\kappa _0^2}{\kappa ^2} - \frac {\kappa _1 \kappa _{30} \beta }{\kappa ^2} \right )\! \left [ \!\frac {\kappa _1^2 \beta }{\kappa ^2} \left ( 1 - e^{\kappa _{12} z} \cos (\kappa _3 z) \right ) - \frac {\kappa _1}{\kappa _{12}} \left ( \!\frac {\kappa _0^2}{\kappa ^2} - \frac {\kappa _1 \kappa _{30} \beta }{\kappa ^2} \!\right )\! e^{\kappa _{12} z} \sin (\kappa _3 z) \!\right ]\!. \end{align}

To calculate the shear stress using (4.23), an expression for $\varPhi _{ij}^{(H)}$ must be assumed. We employ the model (Batchelor 1953)

(4.25) \begin{equation} \varPhi _{ij}^{(H)}(\boldsymbol{\kappa }_0) = \left ( \delta _{ij} - \frac {\kappa _{i0} \kappa _{j0}}{\kappa _0^2} \right ) \frac {E(\kappa _0)}{4 \pi \kappa _0^2}, \end{equation}

where $E(\kappa _0)$ is the energy spectrum of the homogeneous and isotropic turbulence, and since the calculations are inviscid, a suitable assumption for the energy spectrum is that first introduced by Von Kármán (1948), expressed as

(4.26) \begin{equation} E(\kappa _0)= \frac {g_2 q^2 L (\kappa _0 L)^4}{\left [ g_1 + (\kappa _0 L)^2 \right ]^{17/6}}, \end{equation}

where $q$ is a characteristic r.m.s. value of a turbulent velocity component and $L$ is the longitudinal integral length scale. The constants $g_1 \approx 0.558$ and $g_2 \approx 1.196$ ensure that the definitions of $q$ and $L$ are consistent.

Note that while a number of simplifying assumptions have been made which mean that full quantitative agreement with observations is not to be expected, the above theory involves no fitting or particular tailoring to our specific system. The turbulence is described by just two scalar parameters pertaining to the homogeneous bulk, $q$ and $L$ , a crude simplification which nevertheless is expected to capture the most essential features. Importantly, these are quantities obtainable from point measurements in experiments as well as in the field.

4.5. Input parameters

Three quantities must be evaluated in order to obtain numerical values for $\overline {u'w'}$ from (4.23): $L, q$ and $\beta$ .

The simplest to determine is $q$ , which according to RDT satisfies

(4.27) \begin{equation} q^2 = {\frac{2}{3}}e, \end{equation}

assuming isotropic turbulence, where the TKE $e$ is found in table 2.

The exact choice of turbulent integral length scale $L$ involves some amount of judgement. We shall take the most immediate available option and use the streamwise integral scale $L_x^x$ for $L$ , reported for cases 1.A–1.D by Smeltzer et al. (2023). One does well to note that estimating the integral scale from experimental data is not trivial and several methods exist, which tend to yield significantly different results. The difference between the experiments means that no single method of calculating $L_x^x$ can be used for all, so direct quantitative comparison is dubious; however, RDT is only expected to describe the ‘spin-up phase’ cases 1.A–1.D from Experiment 1, which are directly comparable.

The displacement $\beta$ defined in (4.15) is determined by the duration and extent of interaction between waves and turbulence in the time between generation by the active grid, and measurement a distance $L_{\textit{FOV}}$ downstream. For a Gaussian wave packet like those in cases 1.A–1.D, van den Bremer et al. (2019) showed that the final value of $\beta$ after the passage of the group, $\beta _{{f}}$ (subscript ‘ ${f}$ ’: final), takes the form

(4.28) \begin{equation} \beta _{{f}}(z) = 4 \sqrt {\pi } \sigma _0 k_0 (k_0a_{{p}})^2{e}^{2 k_0 z}, \end{equation}

where $\sigma _0$ is the spatial standard deviation of the Gaussian wave group and $k_0a_{{p}}$ the peak steepness. By definition of wave group, $\sigma _0= c_g(k_0) \tau _0$ (see (1.2) and table 2). The expression presumes that the turbulence moving downstream from the active grid has interacted with the whole wave group before reaching the field of view; this is true in our Experiment 1.

It might be interesting to also calculate the values of $\beta _{{f}}$ for the cases with regular waves, which gives an indication of the extent of interaction between waves and turbulence before the turbulent flow reaches the field of view. The cases 2.A–3.B consider regular waves so that the interaction occurs at a near-constant rate during the travel time $t_{\textit{travel}}$ , hence in this case we may approximate

(4.29) \begin{equation} \beta = t_{\textit{travel}}\frac {{\rm d} u_s}{{\rm d} z} = \frac {2k_0L_{\textit{FOV}}}{U_0}(ak_0)^2c(k_0){e}^{2 k_0 z}. \end{equation}

Values of $\beta _{{f}}(0)$ for all cases are given in table 4, and vary by more than an order of magnitude from the longest waves of lowest steepness to the short and steep waves.

4.6. The RDT results

Figure 7(a) shows the profile of the Reynolds shear stress provided by RDT after passage of the Gaussian wave group, using the input parameters specified above. Figure 7(b) shows the vertical derivative of this shear stress, which is proportional to the forcing of the mean current, according to (4.5). The results are expected to represent a similar situation to our Experiment 1 where wave groups were employed and the flow is still in a state of transition when the wave group has passed. One may also note the similarity in shape with the measured currents also for regular waves presented in figure 4, indicative that the same process has been at play.

Figure 7. Results from RDT: (a) profile of the normalised shear stress $\overline {u'w'}/q^2$ as a function of depth $k_0 z$ ; (b) profile of the normalised vertical derivative of the shear stress; (c) the normalised shear stress as a function of $\beta$ for $k_0 z=-0.4335$ , the depth where $\overline {u'w'}/q^2$ attains its maximum magnitude. The dashed line is the 1:1 line, illustrating a linear dependence. (d) The RDT estimates of the anti-Stokes velocity profile after the passage of the wave groups, corresponding to figure 3(b).

Next, figure 7(c) shows the dependence of the shear stress at its maximum on $\beta$ (which is proportional to the square of the wave slope $ak_0$ as (4.28) and (4.29) show). The vertical derivative of the shear stress at the surface, to which the current at the surface must be proportional, is clearly proportional to this maximum. This implies that the current at the surface (departing from a quiescent flow) must be proportional to $(ak_0)^2$ , which is consistent with the results of figure 5 for cases 2.A and 2.B.

A shear stress profile such as depicted in figure 7(a) would, according to (4.5), lead to a current intensifying continuously and indefinitely in time. However, the current will generate a shear stress of opposite sign to that associated with the Stokes drift, and when these shear stresses have evolved so that they cancel each other, the current will reach a steady state (see Pearson 2018). Moreover, dissipation, which has been neglected in the present RDT treatment, is always present, and would act in a similar way to limit the current growth.

4.6.1. Effective interaction time

While the primary goal of our RDT model is to elucidate the physical process behind our observations, we have also used (4.23) to obtain a rough estimate of the current resulting from the wave groups in Experiment 1 (cases 1.A–1.D) encountering the incoming turbulence. In order to obtain quantitative values, we use a simplified procedure: noting that $\partial _z \overline {u'w'} \propto (ak_0)^2$ , we estimate the current after the passage of the group as

(4.30) \begin{align} U_{\textit{RDT}}(z) &\approx -\!\int \partial _z \overline {u'w'}(t) {\rm d} t \approx -(\partial _z \overline {u'w'})_{\textit{RDT}}\int _{-\infty }^\infty \frac {\tilde {a}(t)^2}{a^2_{{p}}} {\rm d} t\notag \\ &\approx -(\partial _z \overline {u'w'})_{\textit{RDT}}T_{{\textit{int},\textit{grp}}}, \end{align}

where we have introduced an `effective interaction time' for a wave group,

(4.31) \begin{equation} T_{{\textit{int},\textit{grp}}} = \frac 1{a_{{p}}^{2}}\int _{-\infty }^\infty \tilde {a}(t)^2{\rm d} t, \end{equation}

an estimate for the duration for which the turbulent current has been influenced by the waves before it reaches the point where it is measured, taking into account that the wave amplitude varies throughout the group. As seen by the turbulence, $\tilde {a}(t)$ is a Gaussian with $\tau _0$ from (3.2) replacing $\tau$ in (3.1), inserting which into (4.31) gives

(4.32) \begin{equation} T_{{\textit{int},\textit{grp}}} = \sqrt {\pi }\tau _0. \end{equation}

We next insert into (4.30) the calculated Reynolds stress vertical gradient shown in figure 7(b) (multiplied by $k_0q^2$ to obtain a dimensional value), using values for $k_0, q=2e/3,L=L_x^x$ and $T_{{int,grp}}$ from table 4. The resulting estimated velocity profiles are shown in figure 7(d).

It is worth noting the relation between the effective interaction time and the fluid displacement $\beta$ . For regular waves, the turbulence interacts constantly with waves throughout its passage from grid to field of view, so the interaction time in the presence of regular waves is $T_{{int,reg}} = L_{\textit{FOV}}/U_0$ . Using (4.28) and noting that for deep-water waves, $\sigma _0=c_g\tau _0=({1}/{2}) c\tau _0$ , we find that for regular waves as well as wave groups it holds that

(4.33) \begin{equation} \beta (z) = 2k_0 T_{\textit{int}} u_s(z). \end{equation}

4.6.2. Comparison with experiment

There is little reason to expect close quantitative agreement between our measured $\Delta U(z)$ in Experiment 1 and $U_{\textit{RDT}}(z)$ because the assumptions behind (4.30) can only hold very early in the ‘spin-up’ period; the relation assumes that $\overline {u'w'}$ retains the vertical shape shown in figure 7(a) throughout the duration of the interaction between wave group and turbulence. Gradually, however, the turbulent flow field begins to ‘push back’ via the growth of other terms in (4.6) until $\partial _z\overline {u'w'}$ eventually becomes very small. Investigating this transient process in detail is a highly relevant question for future study, experimentally and numerically (e.g. in the vein of Guo & Shen 2013).

Despite the naïveté of (4.30) and the numerous assumptions made when applying RDT, the agreement with the measurements of $\Delta U$ in figure 3(b) is not only qualitatively but also quantitatively reasonable. The predictions vary with depth notably faster than the measured currents, which we conjecture is a consequence of the smoothing effect of turbulence as it develops beyond the linear-growth regime described by (4.30), a turbulent mixing which RDT by construction cannot account for. Agreement is comparatively poor in case 1.B where $U_{\textit{RDT}}$ is far smaller than the corresponding $\Delta U$ , which we cannot explain at present (note that this same case displays a surprisingly strong modification of the underlying turbulence considering that the TKE of case 1.B is intermediate, as seen in figure 4 of Smeltzer et al. (2023)).

The key conclusion to be drawn is perhaps that the change in Eulerian current can indeed be ascribed to the interaction between waves and turbulence, to wit, the gradual tilting and stretching of vortices by the Stokes drift, as previously studied by Teixeira & Belcher (2002).

5. Conclusions

We have presented experimental evidence that an Eulerian-mean flow directed opposite to the waves’ propagation direction is created near the water surface when waves propagate atop a turbulent flow, and argued via two different theoretical approaches how the current is the result of waves and turbulence interacting.

Three experiments were conducted, all including waves propagating upstream on initially depth-uniform flows with different turbulent properties. Bespoke turbulence was created with an active grid and velocity fields were measured with PIV. One experiment compared flow conditions before and after the passage of groups of waves, while two studied the mean flow under regular waves, one with low acquisition frequency over a long time, the other at higher frequency in repeated intervals.

Experiments as well as theory show how, when irrotational waves and a turbulent current encounter each other, the combined flow goes through a period of transition until a new quasi-steady state is reached. The fundamental mechanism involved is the rearrangement of horizontal momentum driven by the Reynolds stress which arises when the Stokes drift acts on turbulent eddies, as studied, e.g. by Teixeira & Belcher (2002).

Our experiment with groups of waves studies both the transition period and the final equilibrium state. The former is investigated by allowing the turbulence a limited time to interact with passing wave group. In the two experiments involving regular waves, on the other hand, the final quasi-steady state appears to have been reached by the time the flow is measured.

We present two separate theoretical models which can describe the quasi-steady state and the transition period, respectively. For the latter situation an approximate relation between the mean current shear ${\rm d} \boldsymbol{u}/{\rm d} z$ and the Stokes drift gradient ${\rm d} u_s/{\rm d} z$ is derived following Pearson (2018), valid nearest the surface where $u_s$ is significant. The relation involves the turbulent variances $\overline {u'u'}$ and $\overline {w'w'}$ and shows good enough agreement with experiments to inspire confidence in its use, in the cases where underlying assumptions are satisfied.

A further approximate relation is adapted from Pearson (2018) and relates the rate of change of the Eulerian-mean current to the Reynolds stress $\partial _z\overline {u'w'}$ . Hypothesising that the underlying mechanism is the action of Stokes drift on turbulent eddies, RDT is used to estimate $\overline {u'w'}(z)$ and hence the time-varying $\boldsymbol{u}(z,t)$ in the early phase of interaction. The predictions of RDT are compared with the measurements of $\boldsymbol{u}(z)$ due to passing wave groups, by using measured values of TKE, integral length scale and group envelope width. Qualitative and quantitative agreement with measurements is found, despite a series of simplifying assumptions.