2.1 Mesh motion equation

The mesh motion equations are coupled to the conservation equations of fluid motion to simulate the moving wavy boundary. It is assumed that the disturbance of the propagating wave on the whole field decreases with the increase of the distance from the moving boundary, and the disturbance at infinity is zero. Therefore, for the entire computational domain, the mesh motion can be described by the diffusion equation governing the grid point displacement,

(2.4)\begin{equation} K\frac{\partial^2D_{w,j}}{\partial x_i\partial x_i} = 0,\quad\text{for } j = 1,2,3, \end{equation}

where $({{D}_{w,1}},\,{{D}_{w,2}},\,{{D}_{w,3}})$ represent the displacement components for the $w$th grid point in the streamwise, spanwise and vertical directions, respectively. Here, $K$ denotes the diffusive coefficient related to the position of the grid point that can be defined by $K={{l}^{-w}}$, where $l$ is the average distance from the wall to every grid point, namely the inverse distance.

2.2 The Navier–Stokes equations in a moving frame

The grid point displacement at the boundary can be determined according to (2.1). Hence, the displacement of each grid point can be calculated by (2.4). The displacement vector for the $w$th grid point is $D_{w,i}^{t+\Delta t}-D_{w,i}^{t}$ from $t$ to $t+\Delta t$, and the velocity vector of the moving grid point is thus ${{u}_{w,i}}=({D_{w,i}^{t+\Delta t}-D_{w,i}^{t}})/{\Delta t}$. The fluid motion coupled with the mesh motion can be described by the Navier–Stokes equations in a moving frame:

(2.5)$$\begin{gather} \frac{\partial ({{u}_{i}}-{{u}_{w,i}})}{\partial {{x}_{i}}} = 0, \end{gather}$$

(2.6)$$\begin{gather}\left.\frac{\partial {{u}_{i}}}{\partial t}\right|_{\varsigma } + ( {{u}_{j}}-{{u}_{w,j}} )\frac{\partial {{u}_{i}}}{\partial {{x}_{j}}} = {-}\frac{1}{\rho }\frac{\partial p}{\partial {{x}_{i}}} + \nu \frac{{{\partial }^{2}}{{u}_{i}}}{\partial {{x}_{j}}\partial {{x}_{j}}} + \varPi {{\delta }_{1i}}. \end{gather}$$

Compared with the equations in the stationary coordinate frame, the mesh motion would introduce the time change rate of the velocity in the new coordinate frame $\varsigma$, which is independent of the Lagrangian and Eulerian coordinate systems. The mesh motion also introduces the velocity vector of the moving coordinate frame ${\boldsymbol {u}_{w}}=( {{u}_{w,1}},\,{{u}_{w,2}},\,{{u}_{w,3}} )$ in the mass conservation equation and the nonlinear terms of the momentum conservation equations.

The integral conservation equation of a general tensorial property $\psi$ on a control volume $V$ bounded by a closed surface $S$ with arbitrary motion is

(2.7)\begin{equation} \frac{{\rm d}}{{\rm d}t}\int_{V} \rho \psi\,{\rm d}V + \oint_{S}{\rho \boldsymbol{n}\boldsymbol{\cdot} ( \boldsymbol{u}-{{\boldsymbol{u}}_{w}} )}\psi \,{\rm d}S = {-}\oint_{S}{\rho \boldsymbol{n}\boldsymbol{\cdot} {{\boldsymbol{q}}_{\psi }}\,{\rm d}S} + \int_{V}{{{S}_{\psi }}\,{\rm d}V}, \end{equation}

where $\boldsymbol {u}$ is the fluid velocity, ${{\boldsymbol {u}}_{w}}$ is the velocity vector of the control volume, ${{\boldsymbol {q}}_{\psi }}$ denotes the volume source, ${{S}_{\psi }}$ is the surface source and $\boldsymbol {n}$ is the unit normal vector pointing outward on the boundary surface of a grid (Reference Jasak and TukovićJasak & Tuković 2006). Replacing $\psi$ with the fluid velocity vector $\boldsymbol {u}$ can lead to the integral conservation momentum equations under the moving coordinate frame,

(2.8)\begin{equation} \frac{{\rm d}}{{\rm d}t}\int_{V}{\rho \boldsymbol{u}\,{\rm d}V} + \oint_{S}{\rho \boldsymbol{n}\boldsymbol{\cdot} ( \boldsymbol{u}-{{\boldsymbol{u}}_{w}} )}\boldsymbol{u}\,{\rm d}S = {-}\oint_{S}{\rho \boldsymbol{n}\boldsymbol{\cdot} {{\boldsymbol{q}}_{\boldsymbol{u}}}\,{\rm d}S} + \int_{V}{{{S}_{\boldsymbol{u}}}\,{\rm d}V}. \end{equation}

The second-order finite volume discretization of (2.8) for a structural grid, achieved using the midpoint rule, can convert the surface integral into the sum of area parts. After discretization, the momentum equations can be expressed as

(2.9)\begin{equation} \frac{{{( {{\rho }_{P}}{{\boldsymbol{u}}_{P}}{{V}_{P}} )}^{t + \Delta t}}-{{( {{\rho }_{P}}{{\boldsymbol{u}}_{P}}{{V}_{P}} )}^{t}}}{\Delta t} + \sum_{f}{{{\rho }_{f}}( F-{{F}_{S}} ){{\boldsymbol{u}}_{f}}} = {-}\sum_{f}{{{\boldsymbol{S}}_{f}}\boldsymbol{\cdot} \rho {{\boldsymbol{q}}_{\boldsymbol{u}}}} + {{S}_{\boldsymbol{u}}}{{V}_{P}}, \end{equation}

where the subscript $P$ denotes the grid point, $f$ the face value, ${{V}_{P}}$ the grid volume, $F={{\boldsymbol {S}}_{f}}\boldsymbol {\cdot } {{\boldsymbol {u}}_{f}}$ the fluid flux with ${{\boldsymbol {S}}_{f}}=\boldsymbol {n}{{S}_{f}}$, and the mesh motion flux is ${{F}_{S}}$. As described in (2.9), the introduced time change rate for the time term is converted to the mesh motion flux, which should meet the need for space conservation,

(2.10)\begin{equation} \frac{{\rm d}}{{\rm d}t}\int_{V}{\,{\rm d}V}-\oint_{S}{\boldsymbol{n}\boldsymbol{\cdot} {{\boldsymbol{u}}_{w}}}\,{\rm d}S = 0. \end{equation}

Equation (2.10) should be satisfied for every time step, with its discretizing form expressed as

(2.11)\begin{equation} \frac{V_{P}^{t + \Delta t}-V_{P}^{t}}{\Delta t}-\sum_{f}{{{F}_{S}}} = 0. \end{equation}

It is noted that the mesh motion flux ${{F}_{S}}$ is determined by the volume swept by the moving surface $f$ during the current time step, not by the mesh velocity ${{\boldsymbol {u}}_{w}}$.

The present paper uses the finite volume method to solve the governing equations. The OpenFOAM (open field operation and manipulation) solver is used to accomplish the simulation. It is noted that we here directly solve the full Navier–Stokes equations without a turbulent model, which can be regarded as an approximation of direct numerical simulation (DNS), termed the quasi-DNS. The DNS requires fine mesh, resolving the Kolmogorov scale. However, quasi-DNS loosens this constraint, as verified by Reference Komen, Camilo, Shams, Geurts and KorenKomen et al. (2017) (figures 3 and 5 in their paper). Therefore, to accomplish the simulation accurately, we use a relatively fine mesh (discussed in the supplementary material for grid convergence). In the current solver, the convection terms in the momentum equations are spatially discretized by the second-order upwind difference scheme, and the second-order backward implicit scheme is adopted for time stepping. The finite volume method based on the standard pressure velocity coupling method and PIMPLE algorithm is a variant of the PISO (pressure-implicit with splitting of operators) method to solve the governing equations. The PIMPLE algorithm combines the SIMPLE (semi-implicit method for pressure-linked equations) strategy and the PISO algorithm and regards each time step as a steady-state flow (Reference Gatin, Vukčević, Jasak and RuscheGatin et al. 2017; Reference Gimenez, Aguerre, Idelsohn and NigroGimenez et al. 2019). The standard PISO algorithm is used for the last step when the solution is obtained according to the steady-state algorithm to a certain extent.

2.3 Simulation configuration

The turbulent flow over a propagating wave is simplified as channel flow with a translating wavy bottom, driven by the external force transformed as a pressure gradient to ensure the fixed bulk velocity ${{U}_{0}}={\iint {U\,{\rm d}x\,{\rm d}z}}/{\iint {{\rm d}x\,{\rm d}z}}$, where $U$ is the instantaneous velocity, which is the same as the numerical study by Reference Åkervik and VartdalÅkervik & Vartdal (2019). The Reynolds number based on the bulk velocity and wavelength is $Re={{{U}_{0}}{{\lambda }_{x}}}/{\nu }\approx 7300$ to provide fully developed turbulence. The corresponding wave friction Reynolds numbers ${Re}_{\tau }={{{u}_{\tau }}{{\lambda }_{x}}}/{\nu }$ are shown in table 1. In the wind–wave turbulent boundary layer flow, the wave slope and wave age determine the dynamics of wind–wave interactions. The former is defined by the ratio of amplitude $a$ to the wavelength ${{\lambda }_{x}}$ (or the product of amplitude $a$ and wavenumber $k$). It is noted that the present paper considers the Airy waves with the phase speed $c$, as described in (2.1). The Airy wave is the linear water wave solution and requires a small wave slope (Reference Yang and ShenYang & Shen 2010). The wave slope in the present study is thus $ak=0.1$, 0.13, 0.15, 0.20, 0.25. Three regimes of wind–wave conditions are considered: slow (${c}/{{{u}_{\tau }}}=2.01$, 3.28, 3.43, 3.69); intermediate (${c}/{{{u}_{\tau }}}=16.36$, 16.77, 17.25); and fast (${c}/{{{u}_{\tau }}}=35.80$, 42.83, 58.93, 63.65) waves. The size of the computational domain is $( x,y,z )=( 2{{\lambda }_{x}},{{\lambda }_{x}},0.5{{\lambda }_{x}})$, which is verified to be large enough (see the supplementary material). The grid points are evenly spaced in both streamwise and spanwise directions. In the vertical direction, the grid points are clustered at the wavy wall boundary through exponential transformation to enhance the accuracy for the boundary layer, in which the first layer of the grid meets the need of $\Delta \zeta _{w}^{+}<1$, where + denotes the normalization scaled by viscous unit, namely $\Delta \zeta _{w}^{+}={\Delta {{\zeta }_{w}}{{u}_{\tau }}}/{\nu }$. Here, we used wave-following curvilinear coordinate $(\xi,\zeta )=(x,{(z-\eta )}/{H})$ to show the dimensionless grid scales, where $\eta$ is the wave elevation and $H$ denotes the height of the physical domain (Reference Yang and ShenYang & Shen 2017; Reference Cao, Deng and ShenCao, Deng & Shen 2020). The high grid resolution (HGR) and super high grid resolution (SHGR) cases are used for grid convergence, and cases S4 and S5 are used to compare with the numerical results by Reference Yang and ShenYang & Shen (2010, Reference Yang and Shen2017), discussed in the supplementary material. The total number of grid points for case SHGR is ${{N}_{x}}\times {{N}_{y}}\times {{N}_{z}}=401\times 201\times 151$, while the other cases use ${{N}_{x}}\times {{N}_{y}}\times {{N}_{z}}=251\times 126\times 76$, with the dimensionless grid scale being shown in table 1. It is seen that the grid scale reaches a quasi-DNS (Reference Komen, Camilo, Shams, Geurts and KorenKomen et al. 2017) but cannot reach the Kolmogorov scale. Still, the grid convergence in the supplementary material demonstrates the suitability of the chosen grid scale. The periodic conditions are applied along the streamwise and spanwise directions. The propagating wave boundary is applied to the bottom wall, while the upper wall is set to be a slip wall. After the full development of the turbulent flow, the statistical averaging begins.

Table 1. Parameter settings. The HGR and SHGR cases are used for the verification of grid convergence. The present study considers five wave slopes: $ak$ = 0.1, 0.13, 0.15, 0.2, 0.25. The cases are divided into slow wave (S1–S5), intermediate wave (I1–I3) and fast wave (F1–F4).

2.4 Phase-averaging decomposition

To separate the effects of wave motion and turbulent motion, we use the phase-averaging method to obtain the disturbance feature caused by the waves, which had been widely used in previous investigations of wind waves (Reference Sullivan, McWilliams and MoengSullivan et al. 2000; Reference Kihara, Hanazaki, Mizuya and UedaKihara et al. 2007; Reference Yang and ShenYang & Shen 2010, Reference Yang and Shen2017; Reference Yousefi, Veron and BuckleyYousefi et al. 2020a, Reference Yousefi, Veron and Buckley2021; Reference Wu, Popinet and DeikeWu et al. 2022). The instantaneous quantity can be decomposed into phase-averaging and fluctuating components, and the phase-averaging quantity (also called the wave-coherent quantity) can further be decomposed into ensemble-averaged and wave-induced quantities. The wave-induced quantity represents the disturbance of waves to the above macroscopic flow. Hence, the instantaneous quantity can be decomposed into three parts (triple decomposition):

(2.12)$$\begin{gather} {{f}_{i}}( x,z,t ) = \langle\,{{f}_{i}}( \varphi ,z )\rangle + {{{f}'_i}}( x,z,t ) = \overline{{{f}_{i}}}( z ) + \widetilde{{{f}_{i}}}( \varphi ,z ) + {{f}'_{i}}( x,z,t ), \end{gather}$$

(2.13)$$\begin{gather}\langle\, {{f}_{i}}( \varphi ,z )\rangle = \underset{N\to \infty }{\mathop{\lim }}\,\frac{1}{N}\sum_{n = 0}^{N-1}{{{f}_{i}}( \varphi + 2n\pi ,z )}, \end{gather}$$

(2.14)$$\begin{gather}\overline{{{f}_{i}}}( z ) = \underset{M\to \infty }{\mathop{\lim }}\,\frac{1}{M}\sum_{m = 0}^{M-1}\left\langle\, {{f}_{i}}\left( \varphi + \frac{2\pi m}{M},z \right)\right\rangle. \end{gather}$$

Here $\widetilde {{{f}_{i}}}( \varphi,z )$ is the wave-induced quantity, determined by the difference between the phase-averaging quantity $\langle\, {{f}_{i}}( \varphi,z )\rangle$ and the ensemble-averaged quantity $\overline {{{f}_{i}}}( z )$, and ${{f}'_{i}}( x,z,t )$ represents the turbulent fluctuation. It is worth noting that the present study only considers the single wave. The triple decomposition is a suitable way to separate the flow fields. In the real marine environment, the broad-banded wave fields are common, which means the failure of this decomposition. Reference Hristov, Friehe and MillerHristov, Friehe & Miller (1998) pointed out that the waves in the open ocean are continuously spread throughout spectral scales. Taking a period of the most energetic (peak) mode in the wave spectrum as the characterized period and conducting phase-averaging leads to a strong attenuation due to the destructive interference of multiple modes. They proposed a novel approach, employing an eikonal-like representation of the wave field based on the concept of an analytic signal and the Hilbert transform to identify the wave-coherent component. For monochromatic signals of wave profiles, the results filtered by this method are equivalent to those obtained by phase-averaging.

2.5 Verification of the numerical model

To verify the numerical model, figure 1 shows the comparisons of momentum statistics with the experimental results by Reference Yousefi, Veron and BuckleyYousefi et al. (2020a, Reference Yousefi, Veron and Buckleyb) and numerical results by Reference Yang and ShenYang & Shen (2010, Reference Yang and Shen2017). All quantities are scaled by the friction velocity. We have used the friction velocity determined by extrapolating the ensemble-averaged RSS from the outer layer based on the physical essence of the friction velocity. The downward shifted logarithmic and linear laws compared with the flat-wall boundary layer flow (FWBL) in figure 1(a) agree well with the experimental and other simulation results. In the supplementary material, we also compare the results of channel flow with the direct numerical simulation results by Reference Lee and MoserLee & Moser (2015, Reference Lee and Moser2018). Figure 1(b) shows the balance of momentum flux for case S1. The RSS and wave-induced shear stress (WSS) trends agree qualitatively with Reference Åkervik and VartdalÅkervik & Vartdal (2019) (pressure-driven turbulent flow), although the wave friction Reynolds number ($Re_{\tau }=557$) for case S1 is different from that in Reference Åkervik and VartdalÅkervik & Vartdal (2019) ($Re_{\tau }=200$, 260, 395, 950). It is noted that the Reynolds number is an essential parameter for the scenarios of wind–wave interactions. Reference Åkervik and VartdalÅkervik & Vartdal (2019) provides insights into the friction Reynolds number effect, which covers a range of wave ages and Reynolds numbers, that the budget of momentum flux shows similar pattern for low- and high-friction Reynolds numbers for different wave ages. In fact, the wave age can be written as a ratio of Reynolds number, namely ${c}/{{{u}_{\tau }}}=2\pi ({c}/{k\nu })({\nu }/{{{u}_{\tau }}\lambda _{x}})=2\pi ({Re_{w}}/{Re_{\tau }})$, where $Re_{w}$ is the wave Reynolds number. Since the present results are normalized by the friction velocity, we do not pay much attention to the friction Reynolds number effect. It is also found that the viscous shear stress (VSS) agrees well with Reference Yousefi, Veron and BuckleyYousefi et al. (2020b). Figure 1(c) indicates a slight discrepancy compared with experimental results, which is attributed to the different coordinate frame, as reported in figure 13 by Reference Yousefi, Veron and BuckleyYousefi et al. (2020a). Figure 1(d) shows consistency with the numerical results by Reference Yang and ShenYang & Shen (2010). It is noted that Reference Yang and ShenYang & Shen (2010) considered the stress-driven turbulent Couette flow over a progressive surface wave train, which differs from the pressure-driven turbulent flow in the current study. Therefore, the RSS is close to the wall stress of the upper plane, whereas the present RSS is nearly reduced to zero. Still, the RSS variation in the near wave region is approximately consistent for both cases. Furthermore, we plot the RSS and VSS along the wave propagation direction at the same height as Reference Yang and ShenYang & Shen (2010) and Reference Yousefi, Veron and BuckleyYousefi et al. (2020a) (the RSS locates at the height ${\zeta }/{{{\lambda }_{x}}\approx 0.1}$, and the VSS is obtained by averaging stress measurements between 284 and 664 $\mathrm {\mu }$m above the air–water interface corresponding to the region of ${{\zeta }^{+}}\approx 3.38$–7.5 in the present study), as shown in figure 1(e,f). We see both RSS and VSS basically agree with the previous investigations. It is noted that we use a common method to determine the friction velocity adopted in previous studies (Reference Hamed, Kamdar, Castillo and ChamorroHamed et al. 2015; Reference BoppBopp 2018), that the total wind stress can be determined by extrapolating turbulent stress from the outer layer since it matches the sum of turbulent, wave-coherent and viscous stress in that region. To address this, we simulated additional slow wave cases to highlight the wind stress partition. Figure 2 shows the form drag and viscous stress normalized by $\rho u_{\tau }^2$, where $u_{\tau }$ is determined by this method. For comparison, the data from many previous studies are also plotted (Reference BannerBanner 1990; Reference Banner and PeirsonBanner & Peirson 1998; Reference Sullivan, McWilliams and MoengSullivan et al. 2000; Reference Caulliez, Makin and KudryavtsevCaulliez et al. 2008; Reference Peirson and GarciaPeirson & Garcia 2008; Reference SavelyevSavelyev 2009; Reference Grare, Peirson, Branger, Walker, Giovanangeli and MakinGrare et al. 2013b; Reference Peirson, Walker and BannerPeirson et al. 2014; Reference BoppBopp 2018; Reference Sullivan, Banner, Morison and PeirsonSullivan et al. 2018; Reference Buckley, Veron and YousefiBuckley et al. 2020; Reference Funke, Buckley, Schultze, Veron, Timmermans and CarpenterFunke et al. 2021). Figure 2 suggests that the mean form drags fall within the region overlapping with findings from Reference Peirson, Walker and BannerPeirson et al. (2014) and Reference Funke, Buckley, Schultze, Veron, Timmermans and CarpenterFunke et al. (2021). The viscous stress also falls into the region very close to that observed in previous studies. Moreover, it approaches the dashed line determined by subtracting form drag from the total wind stress. Therefore, the extrapolating RSS from the outer layer to the surface waves can be a viable approach for calculating total wind stress. By evaluating these turbulent high-order statistics, form drag and viscous stress, it is inferred that the present numerical model can provide rational and reliable data.

Figure 1. Comparisons of the momentum statistics between the present study and other numerical and experimental results: EXP, experimental results by Reference Yousefi, Veron and BuckleyYousefi et al. (2020a) and Reference Yousefi, Veron and BuckleyYousefi, Veron & Buckley (2020b) for $ak=0.13$ and ${c}/{{{u}_{\tau }}}=3.69$; DNS1, numerical results by Reference Yang and ShenYang & Shen (2017) for $ak=0.25$ and ${c}/{{{u}_{\tau }}}=2$; DNS2, numerical results by Reference Yang and ShenYang & Shen (2010) for $ak=0.1$ and ${c}/{{{u}_{\tau }}}=2$. (a) Velocity profiles (the law of flat wall boundary layer flow is also shown by the black doubled dotted dashed line). (b) Vertical profiles of the shear stress (momentum flux) components (the dashed line denotes the extrapolation of RSS from the outer layer). Comparisons of RSS and WSS with (c) EXP and (d) DNS results. The discrepancy with EXP is ascribed to the different coordinate frames (it can be seen in figure 13 of Reference Yousefi, Veron and BuckleyYousefi et al. (2020a)). (e) The comparison of RSS with DNS2 along the wave propagation direction, located at the height ${\zeta }/{{{\lambda }_{x}}\approx 0.1}$. (f) The comparison of VSS with EXP along the wave propagation direction.

Figure 2. (a) Mean form drag and (b) viscous stress. The data from previous studies are also plotted (Reference BannerBanner 1990; Reference Banner and PeirsonBanner & Peirson 1998; Reference Sullivan, McWilliams and MoengSullivan et al. 2000; Reference Caulliez, Makin and KudryavtsevCaulliez, Makin & Kudryavtsev 2008; Reference Peirson and GarciaPeirson & Garcia 2008; Reference SavelyevSavelyev 2009; Reference Grare, Peirson, Branger, Walker, Giovanangeli and MakinGrare et al. 2013b; Reference Peirson, Walker and BannerPeirson, Walker & Banner 2014; Reference BoppBopp 2018; Reference Sullivan, Banner, Morison and PeirsonSullivan et al. 2018; Reference Buckley, Veron and YousefiBuckley, Veron & Yousefi 2020; Reference Funke, Buckley, Schultze, Veron, Timmermans and CarpenterFunke et al. 2021).

4.1 Energy transfer between wave-coherent, and turbulent and wave-induced motions (TKE or WKE production)

This section discusses the energy transfer between phase-averaged (wave-coherent), wave-induced and turbulent motions. We here use the two-dimensional statistical field; consequently, $T_{t}$ (denotes the energy transfer between wave-coherent and turbulent motions) and $T_{w}$ (denotes the energy transfer between wave-coherent and wave-induced motions) can be expanded as (see the supplementary material)

(4.1)$$\begin{gather} {{T}_{t}} = \underbrace{-\langle {u}'{u}'\rangle \frac{\partial \langle u\rangle }{\partial x}}_{{{T}_{t,11}}}\underbrace{-\langle {u}'{w}'\rangle \frac{\partial \langle u\rangle }{\partial z}-\langle {u}'{w}'\rangle \frac{\partial \langle w\rangle }{\partial x}}_{{{T}_{t,13}}}\underbrace{-\langle {w}'{w}'\rangle \frac{\partial \langle w\rangle }{\partial z}}_{{{T}_{t,33}}}, \end{gather}$$

(4.2)$$\begin{gather}{{T}_{w}} = \underbrace{-\tilde{u}\tilde{u}\frac{\partial \langle u\rangle }{\partial x}}_{{{T}_{w,11}}}\underbrace{-\tilde{u}\tilde{w}\frac{\partial \langle u\rangle }{\partial z}-\tilde{u}\tilde{w}\frac{\partial \langle w\rangle }{\partial x}}_{{{T}_{w,13}}}\underbrace{-\tilde{w}\tilde{w}\frac{\partial \langle w\rangle }{\partial x}}_{{{T}_{w,33}}}, \end{gather}$$

where ${{T}_{t,11}}$ (${{T}_{t,13}}$ or ${{T}_{t,33}}$) is TKE production through SRNS (RSS or VRNS) performing work on the velocity gradient, and ${{T}_{w,11}}$ (${{T}_{w,13}}$ or ${{T}_{w,33}}$) is WKE production through wave-induced streamwise normal stress (WSS or wave-induced vertical normal stress) performing work on velocity gradient.

Figure 6 shows the total and components of TKE production caused by the interaction between turbulent stresses and phase-averaged shear (also known as the energy transfer between wave-coherent and turbulent motions), with all components normalized by ${u_{\tau }^{4}}/{\nu }$. A slow wave in figure 6(a) generates locally enhanced turbulent production on the leeward side. Conversely, slightly negative production occurs above the windward side, where acceleration induces a favourable pressure gradient, which agrees well with Reference Yousefi, Veron and BuckleyYousefi et al. (2021). However, an intermediate wave in figure 6(e) sharply decreases this production on the leeward side while intensifying it on the windward side. A fast wave in figure 6(i) further strengthens the turbulent production on the windward side but decreases it into a negative value on the leeward side, which denotes turbulent consumption, where energy is transferred from turbulence to wave-coherent motion.

Figure 6. Energy transfer between wave-coherent and turbulent motions. The total and components of TKE production ($T_{t}^{+},T_{t,11}^{+},T_{t,13}^{+},T_{t,33}^{+}$) are plotted for (a–d) slow wave (S1 with $ak=0.13$ and $c/u_{\tau }=3.69$), (e–h) intermediate wave (I1 with $ak=0.13$ and $c/u_{\tau }=17.25$) and (i–l) fast wave (F1 with $ak=0.13$ and $c/u_{\tau }=35.80$). The red dashed line is the critical height.

For the wave ages considered in the present study, $T_{t,11}^{+}$ and $T_{t,13}^{+}$ dominate the total production. Specifically, a negative (positive) $T_{t,11}^{+}$ on the leeward (windward) side and a positive (negative) $T_{t,13}^{+}$ on the leeward (windward) side can be observed across different wave regimes. It is noted that Reference Yousefi, Veron and BuckleyYousefi et al. (2021), by orthogonal curvilinear coordinate frame, showed that $T_{t,13}^{+}$ is always positive. But our simulation results agree with those of Reference Hudson, Dykhno and HanrattyHudson et al. (1996), Reference Yang and ShenYang & Shen (2010, Reference Yang and Shen2017) and Reference Buckley and VeronBuckley & Veron (2019) under the Cartesian coordinate frame. By adjusting the relative magnitude of $T_{t,11}^{+}$ and $T_{t,13}^{+}$, the waves determine whether turbulence is generated or consumed. In figure 6(a–c), $T_{t,11}^{+}$ competes with $T_{t,13}^{+}$, with $T_{t,11}^{+}< T_{t,13}^{+}$ ($T_{t,11}^{+}>T_{t,13}^{+}$) prevailing on the leeward (windward) side in a slow wave. But when the wave becomes fast, as shown in figure 6(j,k), $T_{t,11}^{+}$ is always stronger than $T_{t,13}^{+}$ on both sides. In other words, $T_{t,13}^{+}$ negatively contributes to $T_{t}^{+}$, with normal stress playing a vital role in the energy transfer between wave-coherent and turbulent motions. The intermediate wave in figure 6(f,g) can be regarded as the transitional condition, where the competition starts to change. Overall, a fast wave enhances the streamwise Reynolds normal stress, the windward side's negative RSS and the velocity gradient, as shown in the supplementary material, to produce additional energy transfer between wave-coherent and turbulent motions.

To further investigate the energy transfer between wave-coherent and wave-induced motions, we present the contours of total and components of WKE production (non-dimensionalized by ${u_{\tau }^{4}}/{\nu }$) in figure 7. The positive (negative) total WKE production $T_{w}^{+}$ denotes that the wave-induced motion produces (disrupts) kinetic energy. Figure 7(a) shows that wave-induced motion generates (disrupts) kinetic energy on the windward (leeward) side for a slow wave. However, the rising wave age in figure 7(e,i) results in a relatively symmetrical distribution and strengthened WKE production.

Figure 7. Energy transfer between wave-coherent and wave-induced motions. The total and components of WKE production ($T_{w}^{+},T_{w,11}^{+},T_{w,13}^{+},T_{w,33}^{+}$) are plotted for (a–d) slow wave (S1 with $ak=0.13$ and $c/u_{\tau }=3.69$), (e–h) intermediate wave (I1 with $ak=0.13$ and $c/u_{\tau }=17.25$), and (i–l) fast wave (F1 with $ak=0.13$ and $c/u_{\tau }=35.80$). The red dashed line is the critical height.

The components of WKE production in figure 7 reveal an intriguing energy production mechanism. A slow wave leads to a typical interaction between $T_{w,11}^{+}$ and $T_{w,13}^{+}$, as shown in figure 7(b,c), where $T_{w}^{+}$ is determined through the near wave's WSS and the upper wave-induced streamwise normal stress, with the rare $T_{w,33}^{+}$ contribution in figure 7(d). This agrees with Reference Yousefi, Veron and BuckleyYousefi et al. (2021), who found that $T_{w,11}^{+}$ and $T_{w,13}^{+}$ contribute to the total wave production for high winds (corresponding to slow waves). Nevertheless, a fast wave transforms another interaction mechanism shown in figure 7(j–l): both $T_{w,13}^{+}$ and $T_{w,33}^{+}$ contribute to WKE production. Similarly, the intermediate wave in figure 7(f–h) is treated as the transitional condition, where this interaction mechanism starts to change.

4.2 Energy transfer between turbulence and waves (wave–turbulence exchange)

As described in the supplementary material, ${{W}_{t}}$ represents the energy transfer between wave and turbulence (Reference Reynolds and HussainReynolds & Hussain 1972; Reference Hara and BelcherHara & Belcher 2004; Reference Yousefi, Veron and BuckleyYousefi et al. 2021), which means waves produce or disrupt the TKE, called wave–turbulence exchange. While ${{W}_{w}}$ denotes the WKE transportation, indicating the wave motion generated by the wave-induced shear. As emphasized by Reference Yousefi, Veron and BuckleyYousefi et al. (2021), ${{\bar {W}}_{w}}$ redistributes wave-induced motions. It can be incorporated into the transport term of the kinetic energy density of the wave-induced flow. Some previous studies suggest that ${{\bar {W}}_{w}}$ is negligible compared with other transport terms (Reference Makin and KudryavtsevMakin & Kudryavtsev 1999; Reference Hara and BelcherHara & Belcher 2004; Reference Hara and SullivanHara & Sullivan 2015). Therefore, we only show ${{W}_{t}}$ here and the WKE transportation ${{W}_{w}}$ is discussed in the supplementary material. Following § 4.1 and expanding $W_t$ as

(4.3)\begin{equation} {{W}_{t}} = \underbrace{-\langle {u}'{u}'\rangle \frac{\partial \tilde{u}}{\partial x}}_{{{W}_{t,11}}}\underbrace{-\langle {u}'{w}'\rangle \frac{\partial \tilde{u}}{\partial z}-\langle {u}'{w}'\rangle \frac{\partial \tilde{w}}{\partial x}}_{{{W}_{t,13}}}\underbrace{-\langle {w}'{w}'\rangle \frac{\partial \tilde{w}}{\partial z}}_{{{W}_{t,33}}}, \end{equation}

where ${{W}_{t,11}}$ (${{W}_{t,13}}$ or ${{W}_{t,33}}$) is wave–turbulence exchange through SRNS (RSS or VRNS) performing work on the wave-induced velocity gradient.

Figure 8 shows $W_{t}^{+}$ and its components rescaled by ${u_{\tau }^{4}}/{\nu }$. A positive value represents that the energy is transferred from waves to turbulence, whereas a negative value means the other way around. It is typical for all wave regimes that waves inject (extract) energy into (from) turbulence on the windward (leeward) side, as shown in figure 8(a,e,i). Notably, the enhanced positive $W_{t}^{+}$ on the windward side for a slow wave extends downstream, resulting from the combination of turbulent stress and wave-induced strain rate. Figure 8(a–d) present that a slow wave determines $W_{t}^{+}$ primarily through $W_{t,11}^{+}$ and $W_{t,13}^{+}$. However, as the wave age increases in figures 8(e–h) and 8(i–l), $W_{t}^{+}$ is mainly governed by $W_{t,11}^{+}$. This arises due to the strengthening of SRNS and wave-induced velocity gradient by a fast wave (see the supplementary material), thus facilitating wave–turbulence exchange.

Figure 8. Energy transfer between turbulence and waves (wave–turbulence exchange). The total and components ($W_{t}^{+},W_{t,11}^{+},W_{t,13}^{+},W_{t,33}^{+}$) are plotted for (a–d) slow wave (S1 with $ak=0.13$ and $c/u_{\tau }=3.69$), (e–h) intermediate wave (I1 with $ak=0.13$ and $c/u_{\tau }=17.25$) and (i–l) fast wave (F1 with $ak=0.13$ and $c/u_{\tau }=35.80$). The red dashed line is the critical height.

4.3 Correlation between energy transfer and wave regimes

We further correlate the energy transfer with the wave regimes. Firstly, the averaging of $W_{t}^{+}$ leads to $\bar {W}_{t}^{+}=-{{\overline {\langle {{{{u}'_i}}}{{{{u}'_j}}}\rangle {{{\tilde {S}}}_{ij}}}}^{\,+}}=-{{\overline {\widetilde {{{{{u}'_i}}}{{{{u}'_j}}}}{{{\tilde {S}}}_{ij}}}}^{\,+}}$ and $\bar {T}_{t}^{+}=\bar {M}_{t}^{+}+\bar {W}_{t}^{+}$ (detailed information can be seen in the supplementary material), which means the energy transferred into turbulence can be divided into ensemble-averaged and wave-induced parts. Similarly, $\bar {T}_{w}^{+}=\bar {M}_{w}^{+}+\bar {W}_{w}^{+}$ denotes that the energy transferred into wave motion includes the terms of ensemble-averaged production and transportation.

We integrated the profiles of $\bar {M}_{t}^{+}$, $\bar {M}_{w}^{+}$ and $\bar {W}_{t}^{+}$ along the vertical direction to obtain the local-averaged energy transfer. We also defined the total energy as $TE=\bar {M}_{t}^{+}+\bar {M}_{w}^{+}+\bar {W}_{t}^{+}$. The normalized energy components can be expressed as $L_{M}^{t}={\int {\bar {M}_{t}^{+}\,{\rm d}\zeta }}/{\int {TE\,{\rm d}\zeta }}$, $L_{M}^{w}={\int {\bar {M}_{w}^{+}\,{\rm d}\zeta }}/{\int {TE\,{\rm d}\zeta }}$ and $L_{W}^{t}={\int {\bar {W}_{t}^{+}\,{\rm d}\zeta }}/{\int {TE\,{\rm d}\zeta }}$, respectively. Figure 9(a) shows the normalized energy components vary as a function of the wave age. As the wave age increases, the normalized mean TKE production decreases, while the normalized WKE production shows the reverse trend. The normalized wave–turbulence exchange is slightly enhanced with the rising wave age. It is seen that the WKE production contributes negatively to the total energy for the slow wave condition. This energy deficit is balanced through the mean TKE production enhancement. However, an intermediate wave partitions the total energy into nearly equal WKE and mean TKE productions, potentially demonstrating the intermediate wave age as a transition, repartitioning kinetic energy. Then, the WKE contributes more to the total energy than the mean TKE for a fast wave. Notably, the wave–turbulence exchange in the present study is positive and thus can be added into the mean TKE production term and merged as total TKE production. We show $L_{M}^{t}+L_{W}^{t}$ and $L_{M}^{w}$ in figure 9(b). A transitional variation of TKE and WKE productions at intermediate wave age can be clearly observed.

Figure 9. (a) Normalized energy components. (b) Normalized total TKE (TKE produced by wave-coherent and wave-induced motions) and WKE productions.

4.4 Wave growth rate

Following Reference Cao and ShenCao & Shen (2021), the wave-induced pressure can lead to a form drag on the wave surface. This form drag would result in the growth of wave, usually characterized as wave growth rate parameter $\beta$, which is defined as

(4.4)\begin{equation} \beta = \frac{2{{F}_{p}}}{{{(ak)}^{2}}},\quad \text{where}\ {{F}_{p}} = \int_{{0}}^{{{\lambda }_{x}}}{{{{\tilde{p}}}^{ + }}\frac{\partial \eta }{\partial x}\,{{\rm d} x}}. \end{equation}

Here the $+$ denotes normalization by $\rho u_{\tau }^2$. Figure 10 shows the wave growth rate parameter $\beta$ varying with wave age. For comparison, we plot some experimental and numerical results from Reference PlantPlant (1982), Reference Mastenbroek, Makin, Garat and GiovanangeliMastenbroek et al. (1996), Reference CohenCohen (1997) and Reference Cao and ShenCao & Shen (2021). It is seen that the wave growth rate falls into the region overlapping with previous studies, which further verifies the present numerical model. Moreover, a transition of $\beta$ from positive to negative can be found near the intermediate waves.

Figure 10. Wave growth rate varying with the wave age.