2.1. Numerical method

We consider the three-dimensional turbulent shear flow beneath non-breaking progressive surface waves driven by surface stresses. The water density $\rho$ and the kinematic viscosity $\nu$ are constants. The velocity $\boldsymbol {\upsilon }=(u,\upsilon,w)$ and the dynamic pressure $p$ of the flow field are governed by the solenoidal condition and the Navier–Stokes equations. The nonlinear boundary conditions of mass and momentum conservation, namely the conditions of material surface and continuity of tangential and normal stresses, are satisfied on the wavy surface $z=\eta (x,y,t)$, where $z=0$ is the mean surface. The advection–diffusion equation governing the evolution of the temperature field $\theta$ is also integrated in time with the flow simulation. The temperature is treated as a passive tracer; the buoyancy effect due to temperature fluctuations hence does not modify the vertical momentum equation.

The numerical method for solving the governing equations in a time-dependent domain subject to the fully nonlinear surface conditions was documented in Tsai & Hung (Reference Tsai and Hung2007) where details of the numerics and the implementation for simulations of various wavy flows were reported. The simulation set-up in this study is similar to that of Tsai et al. (Reference Tsai, Chen, Lu and Garbe2013). The important features of the numerical method are summarized here.

To accurately track the Lagrangian free-surface boundary and satisfy the boundary conditions, the time-dependent physical domain in $(x,y,z)$ is mapped to a rectangular computational domain in $(\xi,\psi,\zeta )$ by the algebraic transformation

(2.1a–c)\begin{equation} \xi = x, \quad \psi = y, \quad \zeta = \frac{z + h}{h + \eta(x,y,t)}, \end{equation}

where $h$ is the constant depth, the $x$ and $y$ axes are in the streamwise and spanwise directions, respectively. The bottom and surface of the water column, therefore, are at $\zeta =0$ and 1, respectively. The governing equations and the boundary conditions are transformed to, and solved in, the computational domain. Such a computation strategy employing algebraic mapping has been implemented in Fulgosi et al. (Reference Fulgosi, Lakehal, Banerjee and De Angelis2003), Tsai & Hung (Reference Tsai and Hung2007), Tsai et al. (Reference Tsai, Chen, Lu and Garbe2013), Guo & Shen (Reference Guo and Shen2009), Yang & Shen (Reference Yang and Shen2011), Xuan & Shen (Reference Xuan and Shen2019) and Fujiwara et al. (Reference Fujiwara, Yoshikawa and Matsumura2020) for simulating the three-dimensional, fully nonlinear free-surface flow of viscous fluids.

The numerical model employs spectral discretization for horizontal differentials and finite differencing in the vertical with a mesh fine enough to resolve gravity–capillary waves down to capillary scale and the viscous sublayer immediately beneath the wavy surface. Low-storage Runge–Kutta scheme (e.g. Williamson Reference Williamson1980) is adopted for temporal integration of the Navier–Stokes equations for the velocity $\boldsymbol {\upsilon }(x,y,z,t)$ and the kinematic free-surface boundary condition for the surface elevation $\eta (x,y,t)$. Taking the divergence of the temporally discretized Navier–Stokes equations and applying the solenoidal constraint to the velocity field at the new time interval results in a Poisson equation for the pressure. The pressure Poisson equation is solved before integrating the Navier–Stokes equations to the next time interval. The normal-stress condition is employed as the Dirichlet condition for the pressure at the free surface in solving the Poisson equation. The tangential-stress surface conditions are satisfied implicitly in integrating the horizontal momentum equations at the free surface by rearranging the conditions to evaluate the vertical derivatives of the horizontal velocities. A given uniform heat flux gives rise to a Neumann condition for the temperature field.

The validity and effectiveness of the numerical method have been demonstrated in the previous studies of wave–turbulence interaction (Tsai, Chen & Lu Reference Tsai, Chen and Lu2015; Tsai et al. Reference Tsai, Lu, Chen, Dai and Phillips2017) and the turbulent boundary layer beneath wind waves (Tsai et al. Reference Tsai, Chen, Lu and Garbe2013; Lu et al. Reference Lu, Tsai, Garbe and Jähne2021).

2.2. Simulation set-up

The simulations are initiated by superimposing an ambient fluctuation velocity field onto the velocity field of progressive, monochromatic surface waves from the nonlinear Stokes-wave solution. The surface wave is characterized by the wavelength $\lambda$ and steepness $ak$, where $a$ is half the wave height and $k=2{\rm \pi} /\lambda$ is the wavenumber. The fluctuation velocity field is solenoidal and homogeneous in both the along- and cross-wind directions. Therefore, the only flow structure of the initial velocity field is the two-dimensional wavy motion in the along-wind vertical plane associated with the surface waves. The initial surface wave is imposed to shorten the spin-up time for the wave and flow fields to reach quasi-steady states (discussed in § 2.4). The simulations are carried out for at least 20 surface wave periods before the flow fields are used for analyses.

Tangential and normal surface stresses are applied at the water surface to maintain the propagation of the surface waves and the evolution of the underlying turbulent shear flow. The stress distributions are guided by the measurements of Banner & Peirson (Reference Banner and Peirson1998) and Peirson (Reference Peirson1997), and the analyses of Fedorov & Melville (Reference Fedorov and Melville1998). Despite the idealizations imposed in the numerical simulation, the computed surface elevation and surface-tangential velocity compared well with the measurements of Banner & Peirson (Reference Banner and Peirson1998) as reported in Tsai et al. (Reference Tsai, Chen, Lu and Garbe2013).

The length $L_x$, width $L_y$ and depth $h$ of the computational domain are 4, 2 and 0.8 times the wavelength $\lambda$ of the surface waves and discretized by 512, 256 and 128 grids. A stretched grid system is employed such that the discretization grids cluster when approaching the surface to allow proper resolution of the viscous sublayer adjacent to the upper surface. For the simulations considered (see § 2.3), there are at least ten grids in the near-surface region $|z^+ |<10$, where $z^+=h(1-\zeta )u_*/\nu$ is the distance from the free surface in wall unit and $u_*$ is the mean friction velocity of water at the wavy surface.

2.3. Simulation parameters

In Lu et al. (Reference Lu, Tsai, Garbe and Jähne2021), simulations employing the above-mentioned numerical model were conducted with various initial surface-wave steepness $ak$ and mean surface shear stress $\tau _0=\rho u_*^2$. These simulations are guided by a series of experiments conducted in the annular wind-wave facility Aeolotron at Heidelberg University (Schnieders et al. Reference Schnieders, Garbe, Peirson, Smith and Zappa2013; Schnieders Reference Schnieders2015). The measurement set-ups of the experiments were summarized in Lu et al. (Reference Lu, Tsai, Garbe and Jähne2021). Thermal images at various wind-wave conditions were taken in these experiments with $u_*$ ranging from 0.1 to $1.4\ {\rm cm}\ {\rm s}^{-1}$ as indicated by the black solid and open symbols in figure 2. (Figure 2 is reproduced from the results of figures 4 and 8 in Lu et al. (Reference Lu, Tsai, Garbe and Jähne2021) showing the mean surface streak spacings obtained from experiments and numerical simulations at various wind speeds. These comparisons will be further discussed in § 3.) In the experiments, non-breaking gravity–capillary waves were observed for $0.4\ \mbox {cm}\ \mbox {s}^{-1}\lessapprox u_*< \lessapprox 0.74\ \mbox {cm}\ \mbox {s}^{-1}$. Based on the observations, five scenarios of non-breaking wind-wave flows at $u_*=0.25$, 0.4, 0.5, 0.707 and $1\ \mbox {cm}\ \mbox {s}^{-1}$ are considered in the simulations of Lu et al. (Reference Lu, Tsai, Garbe and Jähne2021), as indicated by the solid red circles in figure 2. The simulation of $u_*=1\ \mbox {cm}\ \mbox {s}^{-1}$, was conducted to demonstrate the impact of non-breaking periodic waves in the high shear condition under which the surface waves observed in the wind–wave flume break. The analyses of surface thermal images by Lu et al. (Reference Lu, Tsai, Garbe and Jähne2021) further reveal that, when $u_* \lessapprox 0.5\ \mbox {cm}\ \mbox {s}^{-1}$, wind shear is not strong enough to maintain turbulence production.

Figure 2. Variations of the dimensional mean streak spacing $\bar {d}$ with the friction velocity $u_*$ obtained from experiments and numerical simulations. The comparisons are reproduced from figures 4 and 8 in Lu et al. (Reference Lu, Tsai, Garbe and Jähne2021), and will be further discussed in § 3. The red solid circles denote the results obtained from the thermal surface images of the numerical simulations. The black solid symbols denote the results obtained from the infrared images taken in the experiments of Schnieders (Reference Schnieders2015) and analysed by Lu et al. (Reference Lu, Tsai, Garbe and Jähne2021). The open symbols are the results reported in Schnieders et al. (Reference Schnieders, Garbe, Peirson, Smith and Zappa2013). The dashed line depicts the scaling $\bar {d}u_*/\nu = \overline {d^+}=100$.

Accordingly, in the following analyses, we focus on the case of $u_*=0.707\ \mbox {cm}\ \mbox {s}^{-1}$ in which both the non-breaking surface waves and the shear turbulence are pronounced, and their mutual interactions are significant. The final steepness of the surface waves $ak \cong 0.22$. The friction Reynolds number $Re_\tau =u_* \lambda /\nu \cong 530$, where $\lambda$ is the wavelength of the surface wave. The ratio of the friction velocity to wave linear phase velocity $u_*/c_0\cong 0.02$, which is of $O((ak)^2 )$. Surface waves with and without surface tension are considered to assess the impact of capillary waves. These simulations are denoted by I1 and I0, meaning intermediate-amplitude waves (denoted by I) with and without surface tension (denoted by 1 and 0, respectively). To examine the effect of carrier-wave steepness, the initial wave amplitude is reduced such that the final $ak \cong 0.135$. This simulation is referred to as L0, meaning low-amplitude wave (denoted by L) without surface tension. The simulation of turbulent flow bounded by a flat surface and driven by the same shear stress is also conducted for comparison with the shear turbulence beneath a wavy surface. This special case is called NW, which stands for simulation with no surface waves.

The flow parameters considered in the present simulations realize the condition in the laboratory experiments. The turbulent Langmuir numbers, $La_t=(u_*/U_s)^{1/2}$ (McWilliams et al. Reference McWilliams, Sullivan and Moeng1997), are approximately 0.65 and 1 for cases I0 and L0, respectively, where $U_s=(ak)^2c_0$, corresponding to the condition of weak wave forcing considered in the typical LESs of oceanic Langmuir turbulence solving the CL equations. However, as shown in the simulation and analysis results of the following sections, the surface waves of cases I0 and L0 are strong enough to generate the large-scale vortical structures of LCs.

2.4. Decomposition and averaging of the flow field

The mobile water surface poses the difficulty of evaluating statistical properties in the fixed coordinate system of the physical domain. As such, analyses of the flow field are carried out in the wave-following coordinate defined by (2.1a–c). Since the simulated flow consists of a turbulent shear layer undulated by the oscillatory motions of the surface waves, a flow property of the instantaneous flow field, $f(x,y,\zeta,t)$, is then decomposed into mean, wave-correlated, and fluctuation components (e.g. Sullivan, McWilliams & Moeng Reference Sullivan, McWilliams and Moeng2000; Tsai et al. Reference Tsai, Chen, Lu and Garbe2013; Xuan et al. Reference Xuan, Deng and Shen2020)

(2.2)\begin{equation} f(x,y,\zeta,t) = \bar{f}(\zeta,t) + \tilde{f}(x,\zeta,t) + f'(x,y,\zeta,t) , \end{equation}

where $f'(x,y,\zeta,t)$ is the fluctuation component. The mean component $\bar {f}(\zeta,t)$ is obtained by taking the spatial average of $f$ over the wave-following constant $\zeta$ plane

(2.3)\begin{equation} \bar{f}(\zeta,t) = \frac{1}{N_{\lambda} \lambda L_y} \sum_{m=1}^{N_{\lambda}} \int_{0}^{L_y} \int_{0}^{\lambda} f(x + m\lambda, y, \zeta, t)\,{{\rm d}\kern0.06em x}\, {{\rm d}y}, \end{equation}

where $N_{\lambda }$ is the number of surface waves in the computation domain. The wave-correlated component

(2.4)\begin{equation} \tilde{f}(x,\zeta,t)=\left\langle \,f \right\rangle_{y}(x,\zeta,t) -\bar{f}(\zeta,t), \end{equation}

where

(2.5)\begin{equation} \left\langle \,f \right\rangle_{y}(x,\zeta,t) = \frac{1}{N_{\lambda} L_y} \sum_{m=1}^{N_{\lambda}} \int_{0}^{L_y} f(x + m\lambda, y, \zeta, t) \,{\rm d}y, \end{equation}

is the phase average of $f$. In some analyses, the ensemble average in time is also employed. They will be noted in the following discussions.

Despite the wind forcings, the wave-correlated velocities of the simulated flow field resemble closely that of a nonlinear Stokes wave, indicating the effectiveness of (2.4) in decomposing the wave-correlated components, $\tilde {u}$ and $\tilde {w}$. The same finding was also reported in the wave-phase-resolved LES of Xuan et al. (Reference Xuan, Deng and Shen2019). For the simulated waves, the wave-correlated velocities decrease to 10 % of their surface values at the depth $z^+ \cong 200$. This depth gives a vertical length scale of the penetration depth of the LCs.

Figure 3 presents the vertical variations of the root-mean-square (r.m.s.) velocities, $\sigma _u$, $\sigma _\upsilon$ and $\sigma _w$, and the r.m.s. fluctuation velocities, $\sigma _{u'}$, $\sigma _{\upsilon '}$ and $\sigma _{w'}$, of case I0 at various time instances from $t=25T_0$ to $30T_0$, where the r.m.s. of quantity $f$ is defined as $\sigma _f = (\kern0.7pt\overline {f^2})^{1/2}$. The profiles of the r.m.s. velocities exhibit strong anisotropy, with $\sigma _u > \sigma _w > \sigma _\upsilon$. In contrast, the fluctuation intensities, $\sigma _{u'}$, $\sigma _{\upsilon '}$ and $\sigma _{w'}$, are compatible with one another; this is distinct from that in turbulent wall flow, in which $\sigma _{u'} > \sigma _{\upsilon '} > \sigma _{w'}$. For Langmuir turbulence considered in the present simulation, the circulatory motion of LCs contributes to the decomposed spanwise and vertical fluctuation components employing decomposition (2.2), resulting in comparable fluctuation velocity partitions $\sigma _{u'} \sim \sigma _{\upsilon '} \sim \sigma _{w'}$. The notable fluctuating velocities, attributed to both shear turbulence and LCs, are confined within the surface layer $z^+ \lessapprox 200$ where the wave motions dominate the flow. The vertical distributions of these averaged fluctuation velocities change insignificantly in time, indicating the Langmuir turbulence reaches the stationary state. In the following discussions, simulated flows within this time interval are analysed.

Figure 3. Vertical distributions of the r.m.s. velocities, $\sigma _u$, $\sigma _\upsilon$ and $\sigma _w$ (a–c), and the r.m.s. turbulent velocities, $\sigma _{u'}$, $\sigma _{\upsilon '}$ and $\sigma _{w'}$ (d–f), of case I0 at six equal-spanning time instances from $t=25T_0$ to $30T_0$. The profiles at $t=30T_0$ are shown in black; the rest are in grey.

4.1. Conditional averaging of LCs

As shown in figure 4(d), the predominant streaks meander, and so do the axes of the large-scale vortical cells that accompany these streaks. Accordingly, streamwise averaging to educe the large-scale vortical structures can be improved by utilizing the decomposed image characterized by predominant streaks, $\theta _{\mathcal {V}_{LC}}$. In the decomposed image $\theta _{\mathcal {V}_{LC}}$, the local temperature minima along various cross-wind transects are searched for as shown in figure 7(b); their horizontal coordinates are denoted by $(x_i,y_i)$. Connecting these coordinates forms the skeleton of the streaks. Skeletons with non-dimensional lengths less than 300 wall units are discarded. A conditional average of a physical quantity $q$ is then obtained by shifting the detected $y_i$ to the cross-wind origin of a new coordinate $\varUpsilon =y-y_i$ and taking the streamwise average

(4.1)\begin{equation} \left[ q \right] (\varUpsilon, \zeta) = \frac{1}{N_t N_s} \sum_{n=1}^{N_t} \sum_{i=1}^{N_s} q(x_i, y-y_i, \zeta, t_n), \end{equation}

where $N_s$ is the total number of detected sections, and $N_t$ is the total number of ensemble time instances. Here, $\varUpsilon = 0$ therefore becomes the symmetric plane of the counter-rotating pair of averaged LCs.

Figure 7. (a) Surface temperature distribution $\theta$ of case I0 and (b) the decomposed distribution characterized by predominant streaks $\theta _{{\mathcal {V}}_{LC}}$. The identified skeleton of the predominant streaks is marked by black dots.

The averaged distributions of the fluctuation velocities on the cross-wind vertical plane ($[ \upsilon ' ], [ w' ]$), and the corresponding streamwise vorticity $[ \omega _{x}' ]$ of the simulations I0 and L0 are shown in figure 8. The result reveals a more orderly cellular structure, consistent with the counter-rotating circulatory flow of LCs: the averaged flow converges near the surface and diverges in the submerged water, resulting in a downwelling beneath the surface streak ($\varUpsilon = 0$) and upwellings on both sides of the structure. The surface converging velocities are more significant than the submerged diverging velocities; the downwelling is stronger than the upwellings, and both decay rapidly with depth.

Figure 8. The averaged distributions of the fluctuation velocities on the cross-wind vertical plane ($[ \upsilon ' ], [ w' ]$), and the corresponding streamwise vorticity $[ \omega _{x}' ]$ of cases I0 (a–c) and L0 (d–f). The ensemble averaging is taken from $t=25T_0$ to $30T_0$ for case I0 and from $t=30T_0$ to $35T_0$ for case L0.

The averaged structure depicted in figure 8 also reveals the characteristic cross-wind length scale of the predominant LCs. The transverse span of a pair of Langmuir cells, $d_s$, can be estimated as the width between the two zero spanwise velocity contours at the water surface. For both cases I0 and L0, the estimated width of an averaged LC pair $d_s\cong 0.34\lambda$ (see figure 8a,d), resulting in a non-dimensional wavenumber $\ell _s=\lambda /d_s \cong 2.9$. Despite the comparable transverse length scales, the maximum velocities and the core vorticity of the averaged LCs of I0 are stronger than those of L0; the maximum $[ \upsilon ' ]$ and downwelling $[ w' ]$ of I0 are approximately 1.3 times those of L0, the maximum $[ \omega _{x}' ]$ of I0 is approximately 1.5 times that of L0.

The transverse length scale of the large-scale vortical structure can also be revealed from the velocity spectrum of the flow field. Figure 9 presents the streamwise averages of the spectral density distributions of the fluctuation velocities, $\langle {\lvert \widehat {u'}\rvert }\rangle _x$, $\langle {\lvert \widehat {\upsilon '} \rvert }\rangle _x$ and $\langle {\lvert \widehat {w'} \rvert }\rangle _x$, and the turbulent kinetic energy, $\langle {(\lvert \widehat {u'} \rvert ^2 + \lvert \widehat {\upsilon '} \rvert ^2 + \lvert \widehat {w'} \rvert ^2)^{1/2}}\rangle _x$, at the surface, where $\langle {\cdot }\rangle _x$ is the streamwise average, $\widehat {\boldsymbol {\upsilon }'}(x,k_y)$ is the Fourier transform of $\boldsymbol {\upsilon }'$ and $k_y$ is the non-dimensional spanwise angular wavenumber. The distributions clearly show peaks at $k_y/2{\rm \pi} =3$ for cases I0 and L0, which is very close to the transverse wavenumber $\ell _s\cong 2.9$ of the averaged vortical structure depicted in figure 8. The amplitudes of the $\upsilon '$ and $w'$ spectra of case I0 are higher than those of case L0, consistent with the averaged structure. In contrast, the velocity spectra are broad-banded for case NW with no surface waves; no peak spectra of velocities are observed.

Figure 9. The streamwise averages of the spectral density distributions of the fluctuation velocities, $\langle {\lvert \widehat {u'}\rvert }\rangle _x$, $\langle {\lvert \widehat {\upsilon '} \rvert }\rangle _x$ and $\langle {\lvert \widehat {w'} \rvert }\rangle _x$, and the turbulent kinetic energy, $\langle {(\lvert \widehat {u'} \rvert ^2 + \lvert \widehat {\upsilon '} \rvert ^2 + \lvert \widehat {w'} \rvert ^2)^{1/2}}\rangle _x$, of cases I0 (a), L0 (b) and NW (c).

4.2. Stability analysis of the CL equations

The results presented so far clearly reveal the possible formation of Langmuir cells beneath the surface waves. It is widely accepted that LC arises through the interaction of the Lagrangian drift of the surface waves with the shear current. Perturbations to the shear current result in distortion of the vortex lines initially comprising spanwise vorticity $\bar {\omega }_y(z)$. Despite the distortion, each vortex line advects at its original velocity plus Stokes drift $U_s$. Since ${\rm d}U_s/{\rm d}z\ne 0$, the distorted vortex lines are then exposed to different values of $U_s$ at different $z$, the vortex line is tilted to realize streamwise vorticity. The equations that govern this process beneath irrotational surface waves in the presence of weak shear were first derived by Craik & Leibovich (Reference Craik and Leibovich1976), and are known as the CL equations. The CL equations are wave averaged in the sense that they exploit the fact that LC evolves over a longer time scale compared with the wave period and so take an average over the wave field. The rectified effect of the wave field is then exposed as the Stokes drift. But for the streamwise vorticity to realize the spanwise and vertical velocities and thus circulatory rolls requires an instability mechanism.

To proceed, following Leibovich & Paolucci (Reference Leibovich and Paolucci1981) and Phillips (Reference Phillips2001), we consider the linearized CL equations, which can be written as

(4.2) \begin{equation} \frac{\partial \mathfrak{\textbf{v}}}{\partial t} + \mathfrak{w} \frac{\partial U_e}{\partial z} \boldsymbol{i} = U_s \boldsymbol{\nabla}_p \mathfrak{u} - \boldsymbol{\nabla}_p \mathfrak{p} + La \boldsymbol{\nabla}_{p}^2 \mathfrak{\textbf{v}}, \end{equation}

where $\mathfrak{\textbf{v}}(y,z,t) = \mathfrak {u}\boldsymbol {i} + \mathfrak {v}\boldsymbol {j} + \mathfrak {w}\boldsymbol {k}$ is the non-dimensional wave-averaged perturbed rotational velocity, $U_e$ is the Eulerian mean velocity, $\mathfrak {p}$ is the wave-averaged perturbed pressure and $\boldsymbol {\nabla }_p = \partial / \partial y \boldsymbol {j} + \partial / \partial z \boldsymbol {k}$. The above process to realize the streamwise vorticity of LCs is captured in the vortex force term $U_s \boldsymbol {\nabla }_p \mathfrak {u}$. In the non-dimensional equation (4.2), the characteristic length scale $\mathcal {L}=k^{-1}$; the streamwise velocity is non-dimensionalized by the characteristic velocity scale $\mathcal {U} = u_*^2\nu _T^{-1} k^{-1}$, where $\nu _T$ is the eddy viscosity; the spanwise and vertical velocities are non-dimensionalized by $\mathcal {V} = \mathcal {U}^{1/2} U_s^{1/2}$; the time is non-dimensionalized by $\mathcal {L/V}$; and the pressure is non-dimensionalized by $\rho \mathcal {V}^2$. This form then exposes the Langmuir number defined as $La = \nu _T \mathcal {V}^{-1} \mathcal {L}^{-1}$.

In the CL equation (4.2), the basic state consists of surface waves of $ak\sim O(\epsilon )$ and a weak Eulerian mean shear flow (Craik Reference Craik1982). The level of shear is quantified by $\mathcal {U}/c\sim O(\epsilon ^s)$, where $c$ is the phase velocity of the waves and the exponent $s\ge 0$ (Phillips Reference Phillips1998, Reference Phillips2001). For a weak mean shear flow, $s=2$, the characteristic scales of the perturbed velocities are $\mathcal {U}\sim O(\epsilon ^2)$ and $\mathcal {V}\sim O(\epsilon ^2)$. For stronger but still weak shear, $s=1$, the CL instability mechanism, i.e. (4.2), is essentially unchanged, apart from scaling factors (Craik Reference Craik1982); the characteristic scales of the perturbed velocities become $\mathcal {U}\sim O(\epsilon )$ and $\mathcal {V}\sim O(\epsilon ^{3/2})$.

The CL equation (4.2) admits two instability mechanisms, CL types 1 and 2 (CL1 and CL2), according to whether the wave drift $U_s$ varies laterally. Since our simulations impose unidirectional surface waves, the averaged large-scale vortical structures observed in figure 8 are likely to be formed by the CL2 instability; and the spanwise width of the averaged LC pair would be close to the wavelength of the most unstable spanwise disturbance.

Transverse instability of the rectified flow is then sought by assuming the disturbance to be spanwise periodic with wavenumber $\ell$ and represented by the normal mode

(4.3a,b) \begin{equation} \mathfrak{\textbf{v}} = \hat{\mathfrak{\textbf{v}}}(z) {\rm Re} \left\{ {\exp} \left( {\rm i}\ell y +\sigma t \right) \right\} \quad \mbox{and}\quad{\mathfrak{p}} = {\hat{\mathfrak{p}}}(z) {\rm Re} \left\{ {\exp} \left( {\rm i}\ell y +\sigma t \right) \right\}, \end{equation}

where the spanwise wavenumber $\ell$ is equivalent to the ratio of the wavelength of surface wave to the wavelength of spanwise disturbance, and $\sigma$ is complex. Substituting (4.3a,b) into the linearized CL equation (4.2), representing $\hat {\mathfrak{\textbf{v}}}$ and $\hat {\mathfrak {p}}$ by Chebyshev polynomials and using the collocation technique results in a generalized complex eigensystem for the truncated polynomial coefficients as the eigenvector and $\sigma$ as the eigenvalue. The transverse disturbance is unstable when the growth rate $\mbox {Re} \{ \sigma \} \equiv \sigma _r > 0$. For the analysis of asymptotic stability considering constant mean velocity shear as $t\rightarrow \infty$, we concur with the results of Tsai et al. (Reference Tsai, Lu, Chen, Dai and Phillips2017) employing the Galerkin technique of Phillips (Reference Phillips2001) up to two decimal places. The critical reciprocal Langmuir number $La_c^{-1}$ and its corresponding critical wavenumber $\ell _c$ computed using the Galerkin and collocation techniques are $(La_c^{-1},\ell _c )=(2.1435,0.3182)$ and $(2.1383,0.3247)$, respectively. Note that these values of $La_c^{-1}$ are higher than the 1.46 and 1.4395 globally stable values reported, respectively, by Leibovich & Paolucci (Reference Leibovich and Paolucci1981) and Phillips (Reference Phillips2001) because of a factor of 2 difference in the magnitude of the Stokes drift.

For stability analysis of the simulated flow, the mean shear velocity $U_e$ is computed from the ensemble average of the mean streamwise velocity $\bar {u}$ when the flow reaches a stationary state. The wave drift profile $U_s$ is calculated by the second-order Stokes drift since the wave-correlated velocities $(\tilde {u},\tilde {w})$ of the simulated flow field closely resemble those of a third-order Stokes wave. The eddy viscosity $\nu _T$ is evaluated by taking the average of $-\overline {u'w'} ( {\rm d}\bar {u} / {\rm d}z )^{-1}$ over the entire computation water column. For both cases I0 and L0, the exponent of the shear index $s$ ranges from 0.8 to 1.3, indicating the basic state consists of a weak Eulerian mean shear flow, and the CL equation (4.2) can be used to explain the formation mechanism of observed large-scale vortical structures.

Stability diagrams of the CL equation showing the range of unstable wavenumber $\ell$ for varying reciprocal Langmuir number $La^{-1}$ of cases I0 and L0 are presented in figures 10(a) and 10(c). The thick, solid lines mark the margin of the neutral stability boundary. Variation of the most unstable wavenumber $\ell$ with $La^{-1}$ is depicted by the thick, red dashed line. Contour lines of various fractions of the maximum growth rate $\sigma _r^{max}$ are superimposed on the instability region by the dashed lines. Also presented in figures 10(b) and 10(d) are the corresponding evolutions of the reciprocal Langmuir number $La^{-1}$ with time. The reciprocal Langmuir number decays rapidly from the unrealistic values of $O(10^{6\sim 8})$ in the early stage of the simulation to 85 at $t\cong 15T_0$ for case I0 and to 80 at $t\cong 23T_0$ for case L0. After this simulation spin-up period, $La^{-1}$ approaches asymptotically to $\sim$40 at $t\cong 30T_0$ for case I0 and to $\sim$35 at $t\cong 40T_0$ for case L0. Therefore, the subsequent evolution of the Langmuir turbulence can be considered as being dominated by the most unstable transverse instabilities incepted at $t\cong 15T_0\sim 16T_0$ with $La^{-1}\cong 80\sim 85$ for case I0, and at $t\cong 24T_0\sim 25T_0$ with $La^{-1}\cong 49.5\sim 52$ for case L0. From the stability diagram, the corresponding wavenumber of the unstable mode with the maximum growth rate is $\ell _{max}\cong 2.7\sim 2.9$ for case I0 and $\ell _{max}\cong 2.65\sim 2.77$ for case L0. These $\ell _{max}$ of the instabilities with the maximum growth rates are close to the wavenumber $\ell _s \cong 2.9$ of the averaged vortical structure depicted in figure 8. The close wavenumbers of the most unstable modes of cases I0 and L0 also explain that the two instabilities grow into vortical structures of similar transverse length scales shown in figure 8.

Figure 10. Stability diagrams of the CL equation showing the range of unstable wavenumber $\ell$ for varying reciprocal Langmuir number $La^{-1}$ of cases I0 (a) and L0 (c). The thick, black solid lines mark the margin of neutral stability. Variation of the most unstable wavenumber with $La^{-1}$ is depicted by the thick, red dashed line; variations of unstable wavenumber with fractions of maximum growth rate are depicted by the thin, black dashed lines ($0.9\sigma ^{max}_r$, $0.7\sigma ^{max}_r$ and $0.5\sigma ^{max}_r$). Temporal evolutions of the reciprocal Langmuir number $La^{-1}$ of cases I0 (b) and L0 (d).

Note that the instability of the linearized CL equations (4.2) depends on the Eulerian mean velocity $U_e$ and the Stokes drift of the waves $U_s$. These two velocity profiles are different for cases I0 and L0. However, the ranges of $\ell _{max}$ of the corresponding $La^{-1}$ are similar for the two waves of different steepnesses. Whether such similarity is due to the coincidence of the spin-up processes in the simulations is not conclusive. In contrast, the two instabilities of similar transverse length scales develop into vortical structures of different strengths for surface waves of different steepnesses. In the CL equations, the cumulative effect of the waves appears as the vortex force, $\boldsymbol{\mathcal {f}}_{\boldsymbol{\!v}}=\boldsymbol {u_s} \times (\boldsymbol {\nabla } \times \mathfrak{v})$, driving the counter-rotating Langmuir cells. Therefore, steeper surface waves accompanying larger wave drift $\boldsymbol {u_s}$ result in a stronger vortex force that intensifies the vortical structures of LC, as revealed by the averaged structures of cases I0 and L0 in figure 8.

It should also be emphasized that the less unstable modes of other wavelengths can also develop into LCs with a wide range of transverse length scales. They are filtered out in the averaging process of § 4.1.

5.1. Eduction of CVSs

There exist various methods to identify CVSs within turbulent shear flows. In these methods, the swirling core of a CVS is defined by mathematical criteria based on the local analysis of the velocity-gradient tensor $\boldsymbol {\nabla }\boldsymbol {\upsilon }$ (e.g. Hunt, Wray & Moin Reference Hunt, Wray and Moin1988; Chong, Perry & Cantwell Reference Chong, Perry and Cantwell1990; Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999) or on the Hessian of the pressure (Jeong & Hussain Reference Jeong and Hussain1995). Despite the differences in the definitions, Chakraborty, Balachandar & Adrian (Reference Chakraborty, Balachandar and Adrian2005) showed inter-relationships among these criteria for various canonical turbulent flows. Chen et al. (Reference Chen, Tsai, Druzhinin and Troitskaya2020) further confirmed such relationships are also valid for undulated turbulent shear flow over progressive surface waves. With no particular preference, we choose the $\lambda _{ci}$ criterion of Zhou et al. (Reference Zhou, Adrian, Balachandar and Kendall1999) to extract the CVSs in this study: a CVS is defined as the region with conjugate complex eigenvalues $\lambda _{cr} \pm {\rm i} \lambda _{ci}$ of the velocity-gradient tensor $\boldsymbol {\nabla }\boldsymbol {\upsilon }$, where the imaginary part of the eigenvalue $\lambda _{ci}$ represents the local swirling strength rate of the CVS. The swirling strength of the small-scale CVSs is higher than that of the large-scale vortical structure of LCs. Employing a low value of $\lambda _{ci}$ would identify large-scale structures but include too many structures of broad-range length scales. Therefore, by choosing the proper value of $\lambda _{ci}$, the small-scale CVSs are retained, and the large-scale structures are excluded. We have experimented with various values of $\lambda _{ci}$ to choose the one that is the least sensitive to the following analysis results.

Figure 11 depicts instantaneous distributions of the $\lambda _{ci}^2 = 0.06$ isosurfaces for the simulated flows of I0, L0 and NW. The red and blue colours of the isosurfaces denote the vortical structures with positive and negative streamwise vorticities, respectively. The distributions reveal the following notable features of CVS distributions: the population of CVSs beneath steeper surface waves (case I0) is higher than that beneath less steep waves (L0) and a flat surface (case NW). Most of the educed CVSs are elongated in the streamwise direction; some horseshoe-shaped vortices, with both upward and downward heads, are also observed. The CVSs tend to accumulate and align along elongated regions, and the distributions depend on wave phase for the flows with waves (cases I0 and L0). In contrast, the CVSs are distributed more randomly without spatial coherence in flow beneath a flat surface (case NW).

Figure 11. Instantaneous distributions of the $\lambda _{ci}^2 = 0.06$ isosurfaces for the flows of cases I0 (a), L0 (b) and NW (c). The red and blue colours of the isosurfaces represent the vortical structures with positive and negative streamwise vorticities, respectively. (d–f) The corresponding skeleton distributions of QSV (blue), forward horseshoe vortex (FHV; red) and reversed horseshoe vortex (RHV; green).

The complex variabilities of different vortical structures result in difficulty in revealing their population partitions and the spatial correlations with the surface waves and Langmuir cells. This calls for a formal scheme to classify the vortical structures and to quantify their distributions.

5.2. Spatial distribution and partition of CVSs

The method developed by Chen et al. (Reference Chen, Tsai, Druzhinin and Troitskaya2020) is adopted here to identify and classify the CVSs of different geometries. Details of the scheme are described in Chen et al. (Reference Chen, Tsai, Druzhinin and Troitskaya2020) and summarized in the following.

Given an instantaneous spatial distribution of $\lambda _{ci}$, an isolated region with $\lambda _{ci}^2$ higher than a given threshold is identified as a vortical structure. The geometry of the vortical structure is then characterized by the skeleton of the isolated $\lambda _{ci}^2$ region. The skeleton of an identified vortical structure is defined by connecting the grid points with the minimum increment in $\lambda _{ci}^2$ of the neighbouring grids. Each identified structure is then classified according to the topological characteristics of its skeleton. Three types of vortical structures, including QSV, FHV and RHV, are identified as schematized in figure 12. The QSV possesses a geometric structure similar to that of the horseshoe vortex but with one dominant leg and one degraded leg; therefore, it can be further classified into four types based on the orientation and geometry properties, as depicted in figure 12. Type-1 and 2 QSVs bend in at the downstream end and can be considered as a degenerate or incomplete FHV, whereas type-3 and 4 QSVs bend in at the upstream end and thus are geometrically similar to an RHV.

Figure 12. Schematics of FHV, RHV and four types of QSV. The surface waves propagate in the $x$ direction.

With the identification and classification scheme, the population partitions of different types of CVSs and their spatial variations can then be quantified by computing the cumulative occurrence count that a spatial coordinate $(x,z^+)$ is occupied by one of the three types of CVSs, as shown in figure 13 for the flows of cases I0 and L0. The distributions are ensemble results computed from the flows beneath four surface waves at 21 independent time instances. The corresponding histogram variations in wave phase $\phi$, which are calculated by summing up the occurrence counts over $z^+$, are also shown in figure 13. Comparing the distributions of cases I0 and L0 reveals that more CVSs form beneath steeper waves, consistent with the qualitative observation from figure 11.

Figure 13. Spatial distributions of the CVSs in the simulated flows of I0 (a–d) and L0 (e–h). The two panels in the top two rows are distributions of accumulative occurrence count of characteristic vortices of FHV (a,e), RHV (b,f) and QSV (c,g), respectively. The bottom panels (d,h) are the corresponding histogram variations in wave phase $\phi$. The wave trough and the wave crest are defined to be $0^{\circ }$ and $180^{\circ }$, respectively.

Similar to the turbulent airflow above the wavy surface (Chen et al. Reference Chen, Tsai, Druzhinin and Troitskaya2020), the distributions of the CVSs in the turbulent aqueous flows beneath surface waves also strongly depend on the wave phase. The QSVs dominate the distribution, but in contrast to airflow, the populations of FHV and RHV are insignificant. The QSVs are observed beneath the entire wavy surface; the distribution initiates near the wave crest ($x / \lambda \cong 0.5$ in figure 13(c,g), and extends upstream along the windward face $(0 < x/\lambda < 0.5)$, subducts below the wave trough $(x / \lambda =0)$ and continues toward the next forward face ${(0.5 < x/\lambda < 1)}$. The peak of the occurrence-count histogram is located near the middle of the windward face ($\theta \cong 90^{\circ }$ in figure 13d,h).

The corresponding distributions of the four types of QSVs are shown in figure 14. Again, type-3 and 4 QSVs, which can be considered as incomplete or degenerate RHVs with elongated vortices bending inward at the submerged upstream end, significantly outgrow type-1 and 2 QSVs.

Figure 14. Spatial distributions of accumulative occurrence count of four types of QSVs for the simulated flows of I0 (a–d) and L0 (e–h).

The above analyses reveal that type-3 and type-4 QSVs accumulate beneath the backward surface of the carrier wave and dominate the formation of CVSs. Figures 11(d)–11(e) depict the instantaneous skeleton distributions of QSVs for the corresponding distributions of $\lambda _{ci}^2$ isosurfaces shown in figures 11(a)–11(b), indicating the QSVs are more populous along several elongated regions; and the QSVs with positive and negative streamwise vorticities tend to align alternatively.

Comparing the QSVs distributions in figures 11(a) and 11(d) with the corresponding surface temperature distributions in figure 6(b) further reveals that the accumulation regions of QSVs roughly coincide with predominant streaks. This can also be demonstrated by examining the transverse distance $d_k$ between the skeleton of the QSV and the nearest skeleton of the predominant streak depicted in figure 7(b). Figure 15 presents the probability density histograms $\mathcal {P}$ and their corresponding cumulative probability distributions $\mathcal {C}_{\mathcal {P}}$ of the non-dimensional distance $d^+_k={d_k}u_*/\nu$ for the identified QSV skeletons in the simulated flows of cases I0 and L0. The results reveal that approximately $50\,\%$ of the QSVs are within the interval of 30 wall units from the predominant streaks.

Figure 15. Probability density histograms (red vertical bars) and cumulative probability distributions (black lines and rectangular symbols) of the non-dimensional distance $d^+_k$ from the skeleton of the QSV to the nearest skeleton of the predominant streak for the flows of cases I0 (a) and L0 (b).

The above observations suggest that the spanwise heterogeneity of the QSV distribution could be attributed to the larger-scale circulatory structure of LCs. In § 6, the averaged spatial relation between QSV and LC, and the interplay between these two vortical structures are examined.

6.1. Variable-interval space averaging of QSVs

To explore the relative spatial locations of the fine-scale QSVs and the large-scale Langmuir cells, the method of variable-interval space averaging (VISA, e.g. Kim & Moin Reference Kim and Moin1986) is combined with the CVS eduction and classification scheme described in § 5. The modified VISA procedure consists of two major steps.

First, the QSVs of the designated type whose skeletons intersect with the plane of $z_0^+ = h(1-\zeta _0)u_*/\nu \cong 13.23$ are identified. The vertical coordinate $\zeta _0\cong 1-13.23\nu /(h u_*)$, where a layer of discretized points is located in the computational domain, is around the skeleton midpoints of most QSVs; $\zeta _0$ is chosen by experimenting with various values to produce the best averaging result. The intersection-point coordinate of the $i$th QSV skeleton is denoted by $(x_{0_i}, y_{0_i}, \zeta _0)$. A realization interval encompassing the identified QSV is defined by choosing a rectangular domain centred horizontally at $(x_{0_i}, y_{0_i})$.

Aligning $(x_{0_i}, y_{0_i}, \zeta _0)$ of the identified intervals, averaging of the velocity field of the realized QSV events is then taken as

(6.1)\begin{equation} \{\boldsymbol{\upsilon}\}({X},\varUpsilon,\zeta) = \frac{1}{N_t N_s} \sum_{n=1}^{N_t} \sum_{i=1}^{N_v} \boldsymbol{\upsilon} \left( x-x_{0_i}, y-y_{0_i}, \zeta, t_n \right), \end{equation}

where $N_v$ is the total number of realizations of QSVs, $N_t$ is the total number of ensemble time instances, $({X}=0,\varUpsilon =0,\zeta =0)$ is the origin of the averaging interval. The VISA velocity field $\{\boldsymbol {\upsilon }\}$ is analysed to educe the flow structures.

Since type-3 and 4 QSVs dominate the population of educed CVSs, as shown in figure 14, the flow intervals encompassing these two types of CVSs are realized. Figure 16 depicts the plane distributions of the averaged velocities ($\{u'\},\{\upsilon '\},\{w'\}$) and the isosurface $\lambda _{ci}^2=0.06$ from different perspective views of case I0, delineating the dominant flow pattern and vortical structure around QSVs. The results deduced from the realizations of type-3 and 4 QSVs are presented in the right and left columns, respectively. The ensemble averaging is taken from flows at 21 time instances from $t=25T_0$ to $30T_0$ of the simulation.

Figure 16. The distributions of the velocities from different perspective views for the VISA flow of case I0 ($ak=0.22$). The averaged QSVs are revealed by the isosurface $\lambda _{ci}^2=0.06$. The red and blue colours of the isosurfaces represent the vortical structures with positive and negative streamwise vorticities, respectively. The results deduced from the realizations of type-3 and 4 QSVs are shown in (e–h) and (a–d), respectively. The coordinates are made non-dimensional by the friction length. The ensemble averaging is taken from flows at 21 time instances from $t=25T_0$ to $30T_0$ of the simulation.

The $\lambda _{ci}^2$ isosurfaces highlight the characteristic geometries of type-3 and 4 QSVs: both averaged vortical structures elongate in the streamwise direction from the submerged $(z^+ \cong 25)$ upstream head to the near-surface downstream end. The averaged structures do not have the submerged spanwise heads; the lengths are shorter than the observed QSVs of the instantaneous flow field in figure 11. These are attributed to the smearing effect of averaging, particularly near the two ends of the QSVs. The averaged structure deduced from the realizations of type-4 QSV (see figure 16a–d) is dominated by an elongated structure of $\{\omega _x\} < 0$ that forms on the right edge (looking downstream, $\varUpsilon < 0$) of a streamwise high-speed streak. A degenerated structure of $\{\omega _x\} > 0$ is attached to the downstream end of the dominant $\{\omega _x\} < 0$ structure and extends toward the left side of the streak. This is due to the coexistence of some type-4 QSVs on the other side of the dominant streaks ($\varUpsilon > 0$). A similar but opposite-orientation pattern of vortical structure is observed in the averaged flow based on realizations of the type-3 QSV (figure 16e–h). These results reveal two distinct features of the QSV structures: first, the type-3 and 4 QSVs form at different but preferable spanwise and vertical locations relative to the Langmuir cells. Second, the type-3 and 4 QSVs do not appear simultaneously in the averaged flow, suggesting the intermittent appearances of the two types of QSVs along both sides of the predominant streaks in the instantaneous flow field.

Both averaged velocity fields based on realizations of type-3 and 4 QSVs reveal flow patterns similar to that of averaged LCs depicted in figure 8: a streak with higher streamwise velocity forms on the water surface, the near-surface flows converge toward the surface streak, downward beneath the surface streak, and diverge in the submerged depth. However, in contrast to the averaged flow of LC, the averaged flow based on QSV realization exhibits spanwise unbalance. For averaging based on type-4 QSV realizations (figure 16d), the flow structure on the opposite side of the QSV $(\varUpsilon > 0)$ remains analogous to that of LCs in figure 8; the near-surface convergence $(\{\upsilon '\} > 0)$ above the QSV and the submerged divergence $(\{\upsilon '\} < 0)$ beneath the QSV, however, are more intensified and skewed toward the type-4 QSV. An opposite unbalanced flow pattern is observed for averaging based on type-3 QSV realization (figure 16h). Such averaged flow structures again imply that in the instantaneous flow field, the vortical structures of $\omega _x > 0$, primarily type-3 QSVs, tend to form on the port side of the dominant thermal streaks. Likewise, the major vortical structures that appear on the starboard side of the dominant streaks are mainly type-4 QSVs with $\omega _x < 0$. Therefore, the resulting averaged flow beneath a dominant streak is formed by a localized swirling flow of either type-3 or 4 QSV superimposed on the larger-scale flow of LCs, as illustrated by the schematics in figures 17(a) and 17(b). Of note is that QSVs can also be observed in the region away from the predominant streaks, as revealed in figures 11 and 15. These QSVs are perhaps formed by the same mechanism as that in wall turbulence.

Figure 17. (a,b) Schematics of vortical structures near Langmuir cells viewing downstream. The counter-rotating circulatory pair denotes the Langmuir cells. The dark-red colour beneath the origin of the $\varUpsilon -z$ axes represents the high-speed surface jet; the wind is in the direction out of the paper. The type-4 and 3 QSVs are depicted by blue and red circles in (a,b), respectively. (c) Schematic of a perturbed streak between the counter-rotating Langmuir cells.

Figure 18 presents the plane distributions of the velocities and the isosurface $\lambda _{ci}^2 = 0.06$ for the averaged flow of case L0 ($ak=0.135$). The averaged velocity distributions and vortical structure of L0 are similar to those of I0 ($ak=0.22$). However, decreasing the steepness of the carrier surface waves weakens the strength of the QSVs; the length of the averaged QSV and the maximum averaged velocities decrease. This observation indicates that the attenuation of LCs attributed to decreasing wave steepness weakens the formation of QSVs in the vicinity of LCs, inferring the role of LCs in enhancing the nearby QSVs.

Figure 18. The distributions of the velocities from different perspective views for the VISA flow of case L0 ($ak=0.135$). The averaged QSVs are revealed by the isosurface $\lambda _{ci}^2=0.06$. The red and blue colours of the isosurfaces represent the vortical structures with positive and negative streamwise vorticities, respectively. The results deduced from the realizations of type-3 and 4 QSVs are shown in (e–h) and (a–d), respectively. The coordinates are made non-dimensional by the friction length. The ensemble averaging is taken from simulated flows at 21 time instances from $t=30T_0$ to $35T_0$ of the simulation.

The above observations validate the proposition raised in § 3 regarding how the steepness of surface waves affects the distribution of surface streak spacings in Langmuir turbulence: the interaction between the surface waves and the shear layer gives rise to the LCs that enhance the formation of QSVs. The QSVs are more populous in the vicinity of LCs, resulting in finer streak spacing. Increasing the steepness of the surface waves strengthens the LCs, thus intensifying the enhancement of LCs on the QSVs. Such a trend is consistent with the finding revealed in the probability distributions of streak spacing shown in figure 5 and the distributions of CVSs presented in § 5.2. The observation sheds light on understanding the wind wave–turbulence interaction by revealing another route of energy transport from surface waves to turbulence. The recent study by Xuan et al. (Reference Xuan, Deng and Shen2020) also highlighted such a wave–turbulence energy transport route.

6.2. Enhanced formation of QSVs by LCs

The spatial correlation between the Langmuir cells and the nearby type-3 and 4 QSVs observed in the previous section indicates the possible interaction between the LCs and the QSVs. The various production budgets of streamwise enstrophy $\varOmega _x = 0.5\omega _x^2$ in the vorticity transport attributed to the flow field is then examined to elucidate such an interaction.

Applying the decomposition (2.2) to the velocity and vorticity fields, and noting that the flow set-up of long-crest waves renders the decomposed components of velocity and vorticity as $u=\bar {u}(z,t)+\tilde {u}(x,z,t)+u'(x,y,z,t)$, $\upsilon \approx \upsilon '(x,y,z,t)$, $w\approx \tilde {w}(x,z,t)+w'(x,y,z,t)$, $\omega _x\approx \omega '_x(x,y,z,t)$, $\omega _y=\bar {\omega }_y(z,t)+\tilde {\omega }_y(x,z,t)+\omega '_y(x,y,z,t)$ and $\omega _z\approx \omega '_z(x,y,z,t)$, the transport equation of streamwise enstrophy $\varOmega _x$ is expressed as

(6.2) \begin{align} \frac{{\rm D}\varOmega_x}{{\rm D}t} &= \underbrace{\omega_x \omega_x \overline{\left( \frac{\partial u}{\partial x} \right)}}_{\bar{S}_{11}} + \underbrace{\omega_x \omega_x \!\widetilde{\,\left( \frac{\partial u}{\partial x}\right)\,}}_{\tilde{S}_{11}}\! +\underbrace{\omega_x \omega_x {\left( \frac{\partial u}{\partial x} \right)'}}_{S'_{11}} +\underbrace{\omega_x \bar{\omega}^{I\!I}_y {\left(\frac{\partial u}{\partial y}\right)'}}_{\bar{T}_{12}} +\underbrace{\omega_x \tilde{\omega}^{I\!I}_y {\left(\frac{\partial u}{\partial y}\right)'}}_{\tilde{T}_{12}} \nonumber\\ &\quad +\underbrace{\omega_x \omega'^{I\!I}_y {\left(\frac{\partial u}{\partial y}\right)'}}_{T'_{12}} +\underbrace{\omega_x \omega_z^I \overline{\left( \frac{\partial u}{\partial z} \right)}}_{\bar{T}_{13}} +\underbrace{\omega_x \omega_z^I \widetilde{\left( \frac{\partial u}{\partial z} \right)}}_{\tilde{T}_{13}} +\underbrace{\omega_x \omega_z^I {\left( \frac{\partial u}{\partial z} \right)'}}_{T'_{13}} \nonumber\\ &\quad +\nu {\nabla}^2 \varOmega_x - \nu \left[ {\left( \frac{\partial \omega_x}{\partial x} \right)}^2 + {\left( \frac{\partial \omega_x}{\partial y} \right)}^2 + {\left( \frac{\partial \omega_x}{\partial z} \right)}^2 \right], \end{align}

where the spanwise and vertical vorticities are expressed as $\omega _y = \partial u/\partial z - \partial w/\partial x \equiv \omega _y^I + \omega _y^{I\!I}$ and $\omega _z = \partial \upsilon /\partial x-\partial u/\partial y \equiv \omega _z^I + \omega _z^{I\!I}$. In deriving (6.2), the contributions of vorticity $\omega ^{I}_y$ to the budget, namely the terms associated with $(\partial u/\partial z)(\partial u/\partial y)$, cancel out the contributions of $\omega ^{I\!I}_z$, namely the terms associated with $-(\partial u/\partial y)(\partial u/\partial z)$ (Brooke & Hanratty Reference Brooke and Hanratty1993). Therefore the component of vertical vorticity $\omega _z^{I\!I}=-\partial u/\partial y$ makes no net contribution to the production budget, highlighting the impact of the flow process that induces the vertical vorticity component $\omega ^I_z=\partial \upsilon /\partial x$. To present the actual contribution of the wave motion, only the variation along the wave-following coordinate $\xi$ is considered in evaluating the derivatives of $x$ and $y$ of streamwise velocity. In (6.2), the first three terms $\bar {S}_{11}$, $\tilde {S}_{11}$ and $S'_{11}$ represent the production of $\varOmega _x$ due to the stretching of $\omega _x$ by mean, wave straining and fluctuating streamwise variation of the streamwise velocity, respectively; $\bar {T}_{12}$, $\tilde {T}_{12}$ and $T'_{12}$ represent the production of $\varOmega _x$ attributed to the turning of spanwise vorticity, $\bar {\omega }^{I\!I}_y$, $\tilde {\omega }^{I\!I}_y$ and $\omega '^{I\!I}_y$, by fluctuating spanwise variation of the streamwise velocity, respectively; $\bar {T}_{13}$, $\tilde {T}_{13}$ and $T'_{13}$ are the productions of $\varOmega _x$ due to the tilting of vertical vorticity $\omega ^{I}_z=\partial \upsilon /\partial x$ by mean, wave straining and fluctuating vertical variation of the streamwise velocity, respectively; the last two terms are viscous diffusion and dissipation, respectively.

Note that Lu et al. (Reference Lu, Tsai, Garbe and Jähne2021) explained the impact of surface waves on the scaling of streak spacing by comparing the difference in enstrophy transport beneath the gravity–capillary waves and a flat surface without knowing the interaction mechanism between the LC vortical structures and QSVs in Langmuir turbulence. We reappraise the analysis of streamwise enstrophy transport in Lu et al. (Reference Lu, Tsai, Garbe and Jähne2021) to validate the interaction process depicted in figure 17 and described above. The analysis is similar to that in Lu et al. (Reference Lu, Tsai, Garbe and Jähne2021), but the interpretations are new.

Figure 19 presents the phased-averaged distributions of the production terms in the enstrophy transport equation (6.2) for the simulated flow of case I0. For comparison, the distributions of the total production budget $S_{11}+T_{12}+T_{13}$, the production budget by turning $T_{12}$ and the production budget by tilting $T_{13}$ are also depicted in figures 19(j), 19(k) and 19(l), respectively, where $S_{11}=\bar {S}_{11}+\tilde {S}_{11}+S'_{11}$, $T_{12}=\bar {T}_{12}+\tilde {T}_{12}+T'_{12}$ and $T_{13}=\bar {T}_{13}+\tilde {T}_{13}+T'_{13}$. The distributions on the $x-z$ plane are obtained by taking the phase average (2.5). The corresponding vertical variations of the gross contributions by the various terms are shown in figures 20(a)–20(c). The gross contribution is computed by taking the spatial average (2.3).

Figure 19. Phase-average distributions of various production terms in the transport equation of enstrophy $\varOmega _x$ for the flow of case I0; (a) $\bar {S}_{11}$, (b) $\bar {T}_{12}$, (c) $\bar {T}_{13}$, (d) $\tilde {S}_{11}$, (e) $\tilde {T}_{12}$, (f) $\tilde {T}_{13}$, (g) $S'_{11}$, (h) $T'_{12}$, (i) $T'_{13}$, (j) $S_{11}+T_{12}+T_{13}$, (k) $T_{12}$ and (l) $T_{13}$. The distributions are ensemble results computed from the flows under four carrier waves at 21 independent time instances from $t=25T_0$ to $30T_0$.

Figure 20. Vertical variations of various enstrophy $\varOmega _x$ production terms averaged over the $\xi$–$\psi$ plane for flows of cases I0 (upper row, a–c) and L0 (lower row, d–f); $S_{11}=\bar {S}_{11}+\tilde {S}_{11}+S_{11}'$, $T_{12}=\bar {T}_{12}+\tilde {T}_{12}+T_{12}'$ and $T_{13}=\bar {T}_{13}+\tilde {T}_{13}+T_{13}'$.

The profiles of the gross contributions in figure 20 indicate that the most significant production to the streamwise enstrophy is contributed by the $\bar {T}_{13}$ term. This production occurs within the layer beneath the windward surface, as shown in figure 19(c). This is also the region where the QSVs accumulate, as depicted in figures 13 and 14. The VISA results of figures 16 and 18 further reveal that the QSVs tend to form in the edge vicinity of the surface streaks induced by the LCs, with small tilting and turning angles relative to the streaks. These observations, therefore, suggest the following interaction process between the QSVs and the streaks in enhancing streamwise enstrophy production. A schematic helping explain the process is depicted in figure 17(c).

Suppose that the background flow is turbulent such that the fluctuations act over time scales that are small with respect to the LCs, and the QSVs have somehow formed to realize the respective velocity fields of LCs and mean shear flow. The counter-rotating cells of LCs induce high-speed streamwise streaks in the near-surface region and introduce spanwise gradient of the streamwise velocity $\partial u/\partial y<0$ and $\partial u/\partial y>0$ on the port and starboard sides of the streaks, respectively, as illustrated in the schematics of figures 17(a) and 17(b). The QSVs in the edge vicinity of the streak, of which the downstream near-surface portions tilt upward and turn toward the streak, perturb the nearby flows. Specifically, the type-3 QSV, which tends to form on the port side of the streak, induces $\partial \upsilon /\partial x=\omega _z^I>0$ in the flow region above the upper part of the QSV, as depicted in figure 17(c). Similarly, the type-4 QSV on the starboard side of the streak induces $\partial \upsilon /\partial x<0$. The vertical vorticity in these localized regions is tilted by vertical shear through $\omega _x\omega _z^I \overline {(\partial u/\partial z)}$ and therefore contribute to the production of streamwise enstrophy.

To validate the above proposition regarding the production process of streamwise enstrophy arising from the interaction between the QSVs and the LC streak, we present in figure 21 the distributions of related flow properties in the VISA flows around the QSVs. The distributions of the spanwise gradient of streamwise velocity clearly show that type-3 and 4 QSVs form on the port ($\{\partial u'/\partial y\}>0$) and starboard ($\{\partial u'/\partial y\}<0$) sides of the streaks, respectively; they also reveal the streaks being undulated to be parallel with the QSVs. Consequently, strong streamwise gradient of spanwise velocity occurs in a localized region around the near-surface part of the QSVs, as shown in figures 21(b) and 21(f). The resulting contribution to the production $\omega _x (\partial \upsilon /\partial x)\overline {(\partial u/\partial z)}$ is depicted in figures 21(d) and 21(h). Note that such a production process is not observed in the flow regions near the streaks without the QSVs, therefore underlining the interaction between the QSVs and LCs.

Figure 21. The distributions of $\{-\partial u'/\partial y\}$ (a,e), $\{\partial \upsilon /\partial x\}$ (b,f), $\{\omega _x\}$ (c,g) and $\{\omega _x\}\{\partial \upsilon /\partial x\}$ (d,h) at $z^+=3.29$ for the VISA flow of case I0. The averaged QSVs are revealed by the isosurfaces $\lambda _{ci}^2=0.06$. The red and blue colours of the isosurfaces represent the vortical structures with positive and negative streamwise vorticities, respectively. The results deduced from the realizations of type-3 and 4 QSVs are shown in the right and left columns, respectively. The coordinates are made non-dimensional by the friction length. The ensemble averaging is taken from simulated flows at 21 time instances from $t=25T_0$ to $30T_0$ of the simulation.

The distribution of $\bar {T}_{13}$ production initiates from the wave crest and extends upstream in the layer beneath the windward surface. In the same region, stretching production due to wave straining, $\tilde {S}_{11}$, is also significant, as shown in figure 19(d). The combined effect of these two mechanisms intensifies the production of streamwise enstrophy and enhances the formation of QSVs there. These QSVs are convected upstream in the $-x$-direction relative to the wave (Xuan et al. Reference Xuan, Deng and Shen2019), and accumulate in the surface layer extending from the windward surface to the wave trough, as revealed by the spatial distributions of accumulative occurrence count of QSVs depicted in figure 13(c), and figures 14(c) and 14(d). Beneath the leeward surface, compression straining of wave results in negative production of streamwise enstrophy, $\tilde {S}_{11}<0$. However, the gross contribution of $\tilde {S}_{11}$ is positive, as shown in in figure 20(a), indicating the wave stretching is more significant than the wave compression in the production of streamwise enstrophy.

The vertical profiles of the gross contributions to streamwise enstrophy by various mechanisms for the flow of case L0 are presented in figures 20(d)–20(f). Comparing the results of I0 and L0 reveals that decreasing the steepness of the surface wave reduces all the contributions to the production of streamwise enstrophy, in particular, the $\tilde {S}_{11}$ and $\bar {T}_{13}$ terms. Reduction in the contribution attributed to stretching/compression of wave straining ($\tilde {S}_{11}$ term) is evident since the streamwise velocity associated with wave motions ($\tilde {u}$) decreases with the steepness $ak$. Reducing the wave steepness also decreases the Stokes drift and attenuates the LCs and the accompanying streaks. This weakens the interaction between the streaks and the nearby QSVs, consequently decreasing the production owing to the interaction ($\bar {T}_{13}$ term).