2.1. Governing equation for wave-induced airflow

We first define what constitutes a wave-induced or wave-coherent airflow field by employing a triple decomposition of an instantaneous turbulence field, $f$ . This decomposition splits $f$ into three components, namely a mean field $\langle f\rangle$ , a wave-coherent component $\widetilde {f}$ and a wave-independent or turbulent component $f^\prime$ , following the triple decomposition of Hussain & Reynolds (Reference Hussain and Reynolds1970)

(2.6) \begin{align} f(x,y,z,t)=\left \langle f\right \rangle (\zeta ) + \widetilde {f}(\xi ,\zeta ) +f^\prime (x,y,z,t). \end{align}

The mean $\langle f \rangle (\zeta )$ and wave-induced $\widetilde {f}(\xi ,\zeta )$ components of the triple decomposition are obtained by phase-averaging instantaneous snapshots of the turbulent flow field. It is more convenient to express them in mapped coordinates ( $\xi ,\zeta$ ), which fit the wave surface geometry. Following studies such as Yang & Shen (Reference Yang and Shen2010, Reference Yang and Shen2011), Yang et al. (Reference Yang, Meneveau and Shen2013), Cao et al. (Reference Cao, Deng and Shen2020) and Cao & Shen (Reference Cao and Shen2021), the phase averaging of an instantaneous flow field $f$ is defined as

(2.7) \begin{align} \overline {f}(\xi , \zeta ) = \frac {1}{N_t} \frac {1}{N_y} \sum _{n=1}^{N_t} \sum _{j=1}^{N_y} f(\xi - c t_n, \psi _j, \zeta ). \end{align}

Here, $\overline {f}(\xi , \zeta )$ denotes the phase-averaged quantity, which is obtained by first performing a spanwise average along the $\psi$ -direction followed by time averaging of the phase-shifted field. The parameters $N_t$ and $N_y$ are the number of time snapshots and the number of grid points in the spanwise direction that are used for phase averaging, respectively. The mean $\langle f \rangle$ is obtained by averaging the phase-averaged field in the $\xi$ -direction, and the wave-induced component is obtained from $\widetilde {f}(\xi , \zeta ) = \overline {f}(\xi , \zeta ) - \langle f \rangle (\zeta )$ (Yang & Shen Reference Yang and Shen2010).

Considering the periodic nature of the wavy boundary and the definition of phase-averaged expressions, we analyse the flow field within a wavelength. The wave-induced fluctuations ( $\widetilde {f}$ ) can be represented as a Fourier series in the form of a complex exponential given by

(2.8) \begin{align} \widetilde {f} = \sum \limits _{n=-\infty }^{\infty } \hat {f}_n \mathrm{e}^{\text{i} n k \xi }, \end{align}

where $\hat {f}_n$ denotes the complex Fourier coefficients. Assuming that $\widetilde {f}$ is dominated by its fundamental mode ( $n=-1,1$ ) with wavenumber $k$ , then

(2.9) \begin{align} \widetilde {f} = \hat {f}_1 \mathrm{e}^{\text{i} k \xi } + \hat {f}_1^* \mathrm{e}^{-\text{i} k \xi } + \text{harmonics}. \end{align}

Here, we have applied the identity $\hat {f}_{-1}=\hat {f}_{1}^*$ , which holds for real-valued functions $f$ . The superscript ‘ $*$ ’ denotes the complex conjugate. The phase difference between the wave-coherent field $\widetilde {f}$ and phase-averaged monochromatic wave surface $\widetilde {\eta }=a\cos (k\xi )$ is then calculated via

(2.10) \begin{align} \widetilde {f} = 2 \textrm {Re} \big\{\hat {f} \big\} \cos (k\xi ) - 2 \textrm {Im} \big\{\hat {f} \big\} \sin (k\xi ) = -2 \big|\hat {f} \big| \sin \big(k\xi -{\phi _{\widetilde {f}\, \widetilde {\eta }}}\big). \end{align}

In (2.10), we assume that $\hat {f} \approx \hat {f}_1$ ; therefore, the phase difference is computed as

(2.11) \begin{align} {\phi _{\widetilde {f}\, \widetilde {\eta }}} = \arctan \left (\frac {\textrm {Re}\{\hat {f}\}}{\textrm {Im}\{\hat {f}\}}\right )\!. \end{align}

By applying the triple decomposition from (2.6) to the NS equations (2.5a ) and (2.5b ) and removing the horizontally averaged mean quantities, we can obtain, similar to Cao et al. (Reference Cao, Deng and Shen2020), the following governing equations for the wave-coherent motions:

(2.12) \begin{align} \underbrace {\frac {\partial }{\partial \xi _j} \big ( \widetilde {u}_i \langle U_{\!j} \rangle + \left \langle u_i\right \rangle \widetilde {U}_{\!j}\big )}_{\textbf {I}} + \underbrace {\frac {\partial \tau _{\textit{ij}}^w}{\partial \xi _j}}_{\textbf {II}} + \underbrace {\frac {\partial \widetilde {\tau }_{\textit{ij}}}{\partial \xi _j} }_{\textbf {III}} &= -\underbrace { \frac {\partial \widetilde {\tau }_{\textit{ij}}^p}{\partial \xi _j}}_{\textbf {IV}} -\underbrace { \frac {\partial \widetilde {\tau }_{\textit{ij}}^\nu }{\partial \xi _j}}_{\textbf {V}}. \end{align}

The enumerated terms of (2.12) represent either the advection of the wave-coherent component or the wave-coherent forcing terms. In particular, term $\textbf{I}$ represents advection by the mean flow, term $\textbf{II}$ represents a force caused by the divergence of wave-coherent stresses ( $\tau _{\textit{ij}}^w=\widetilde {u}_i \widetilde {U}_{\!j}$ ), term $\textbf{III}$ represents a force caused by the divergence of wave-coherent turbulence stresses ( $\widetilde {\tau }_{\textit{ij}}=\widetilde {{u\!}_i^\prime U_{\!j}^\prime }={\tau }_{\textit{ij}}-\langle {\tau }_{\textit{ij}}\rangle$ , where ${\tau }_{\textit{ij}}=\overline {u_i^\prime U_{\!j}^\prime }=J^{-1}\overline {u_i^\prime u_l^\prime }\partial \xi _j/\partial x_l$ ), term $\textbf{IV}$ represents a force caused by wave-coherent pressure stress and term $\textbf{V}$ represents viscous forces in the curvilinear coordinate system.

2.2. Turbulence curvilinear model

By combining the leading-order terms of the momentum equation (2.12), which include the turbulence stress $(\widetilde {\tau }_{\textit{ij}})$ terms but ignore the higher-order $\mathcal{O}((ak)^2)$ wave-induced stress $(\tau _{\textit{ij}}^w)$ terms (Cao et al. Reference Cao, Deng and Shen2020), with the wave-coherent continuity equation

(2.13) \begin{align} \frac {\partial \widetilde {U}_i}{\partial \xi _i} = 0, \end{align}

we establish a basis for the linearised turbulence curvilinear model for wave-induced airflow, which is presented in (2.14a – c ) below

(2.14a) \begin{align} \!\!( \langle u \rangle -c)\frac {\partial \tilde {u}}{\partial \xi }+ \big(\tilde {w}+(\langle u \rangle -c )g\tilde {\eta }_\xi \big)\frac {\partial \langle u \rangle }{\partial \zeta }&=-\frac {\partial \tilde {p}_*}{\partial \xi }+\nu \! \left ( \!\frac {\partial ^2\tilde {u}}{\partial \zeta ^2}-\frac {\partial ^2\tilde {w}}{\partial \xi \partial \zeta } \!\right ) \!+\nu \frac {\partial }{\partial \zeta } \!\left ( \! g_\zeta \frac {\partial \langle u \rangle }{\partial \zeta } \!\right ) \nonumber \\ &\quad-\frac {\partial \widetilde {\tau }^d_{11}}{\partial \xi }-\frac {\partial \widetilde {\tau }^d_{13}}{\partial \zeta }, \\[-12pt] \nonumber \end{align}

(2.14b) \begin{align} (\langle u \rangle -c)\frac {\partial \tilde {w}}{\partial \xi }+\frac {\partial \tilde {p}_*}{\partial \zeta }&=\nu\! \left (\frac {\partial ^2\tilde {w}}{\partial \xi ^2}+\frac {\partial ^2\tilde {w}}{\partial \zeta ^2}\right )-\frac {\partial \widetilde {\tau }^d_{31}}{\partial \xi }-\frac {\partial \widetilde {\tau }^d_{33}}{\partial \zeta }, \\[-12pt] \nonumber\end{align}

(2.14c) \begin{align} \frac {\partial \tilde {u}}{\partial \xi }+g\tilde {\eta }_\xi \frac {\partial \langle u \rangle }{\partial \zeta }+\frac {\partial \tilde {w}}{\partial \zeta }&=0. \\[9pt] \nonumber\end{align}

Here, $\widetilde {\tau }^d_{\textit{ij}}=\widetilde {\tau }_{\textit{ij}}-\delta _{\textit{ij}}\widetilde {\tau }_{kk}/3$ is the deviatoric part of the wave-induced Reynolds stress tensor, and $\widetilde {p}_*=\widetilde {p}/\rho +\widetilde {\tau }_{kk}/3$ is the modified wave-induced pressure. For clarity and brevity, we omit the asterisk notation ‘ $*$ ’ hereafter and refer to $p_*$ simply as the pressure $p$ . Additionally, $\widetilde {\eta }=\overline {\eta }=a\cos k\xi$ is the phase-averaged expression of the wave elevation $\eta (x,y,t)$ given by (2.2a ). Details regarding the derivation of these equations without the turbulence stress terms are provided in the supplementary materials, available at https://doi.org/10.1017/jfm.2025.11015, of Cao et al. (Reference Cao, Deng and Shen2020). In the present study, we extend the formulation presented by Cao et al. (Reference Cao, Deng and Shen2020) by retaining the turbulence stress terms. By combining (2.14a – c ) and applying a Fourier transform to the resulting equation, we obtain the following equation for $\hat {w}(\zeta )$ :

(2.15) \begin{align} \left (\langle u \rangle -c){\nabla} ^2-\frac {\partial ^2 \langle u \rangle }{\partial \zeta ^2}-\frac {\nu }{\mathrm{i}k}{\nabla} ^4\right )\hat {w}&=\hat {\eta }\nu \frac {\partial ^2}{\partial \zeta ^2}\left (g\frac {\partial ^2\langle u \rangle }{\partial \zeta ^2}\right )+\mathrm{i}k\frac {\partial }{\partial \zeta } \big[\hat {\tau }^d_{11}-\hat {\tau }^d_{33} \big]\nonumber \\& \quad +\frac {\partial ^2\hat {\tau }^d_{13}}{\partial \zeta ^2} +k^2\hat {\tau }^d_{31}, \end{align}

where ${\nabla} ^2=\partial ^2/\partial \zeta ^2-k^2$ is the Laplacian operator. The left-hand side of (2.15) has the same form as that of the Orr–Sommerfeld equation (Lin Reference Lin1955; Orszag Reference Orszag1971). We also note that retaining the wave-induced stresses, $\tau _{\textit{ij}}^{w} = \widetilde {u}_i \widetilde {U}_{\!j}$ , which are assumed to be negligible ( $\sim \!\mathcal{O}((ak)^2)$ ) in this study, would render (2.15) nonlinear, necessitating an iterative solution procedure to obtain the $\widetilde {u}$ and $\widetilde {w}$ disturbances. Such an iterative approach that includes the wave-induced stresses has been explored in Cao et al. (Reference Cao, Liu, Xu and Deng2023), although the wave-induced Reynolds stresses were still considered negligible in that formulation. Solving (2.15) is equivalent to solving the following coupled equations for the real and imaginary parts of $\hat {w}$ :

(2.16) \begin{align} \left (( \langle u \rangle -c){\nabla} ^2-\frac {\partial ^2\langle u \rangle }{\partial \zeta ^2}\right )\textrm {Re}\{\hat {w}\}-\frac {\nu }{k}{\nabla} ^4 \textrm {Im}\{\hat {w}\}&=\hat {\eta }\nu \frac {\partial ^2}{\partial \zeta ^2}\left (g\frac {\partial ^2 \langle u \rangle }{\partial \zeta ^2}\right )\nonumber \\& \quad -k\frac {\partial }{\partial \zeta }\big(\textrm {Im}\big\{\hat {\tau }^d_{11} \big\}-\textrm {Im}\big\{\hat {\tau }^d_{33} \big\} \big)\nonumber \\& \quad +\frac {\partial ^2}{\partial \zeta ^2}\textrm {Re}\big\{\hat {\tau }_{13}^d \big\}+k^2\textrm {Re}\big\{\hat {\tau }_{31}^d \big\}, \\[-12pt] \nonumber \end{align}

(2.17) \begin{align} \left ((\langle u \rangle -c){\nabla} ^2-\frac {\partial ^2\langle u \rangle }{\partial \zeta ^2}\right )\textrm {Im}\{\hat {w}\}+\frac {\nu }{k}{\nabla} ^4 \textrm {Re}\{\hat {w}\}&=k\frac {\partial }{\partial \zeta }\big(\textrm {Re}\big\{\hat {\tau }^d_{11}\big\}-\textrm {Re}\big\{\hat {\tau }^d_{33} \big\} \big)\nonumber \\& \quad +\frac {\partial ^2}{\partial \zeta ^2}\textrm {Im}\big\{\hat {\tau }_{13}^d \big\}+k^2 \textrm {Im} \big\{\hat {\tau }_{31}^d \big\}. \\[9pt] \nonumber \end{align}

Next, the deviatoric part of the wave-induced Reynolds stress terms in (2.15) can be approximated using the Boussinesq eddy viscosity model, which assumes the following form:

(2.18) \begin{align} \widetilde {\tau }^d_{\textit{ij}} &=\begin{bmatrix} 2\nu _T\dfrac {\partial \widetilde {w}}{\partial \zeta }&-\nu _T\!\left (\dfrac {\partial \widetilde {u}}{\partial \zeta }+\dfrac {\partial \widetilde {w}}{\partial \xi }\right )-\nu _T\dfrac {\partial \langle u \rangle }{\partial \zeta }g_\zeta \widetilde {\eta },\\ &\\ -\nu _T\!\left (\dfrac {\partial \widetilde {u}}{\partial \zeta }+\dfrac {\partial \widetilde {w}}{\partial \xi }\right )&\nu _T\!\left (\dfrac {\partial \widetilde {u}}{\partial \xi }-\dfrac {\partial \widetilde {w}}{\partial \zeta }\right ) \end{bmatrix}\!. \end{align}

This expression retains the same form as that of the viscous stress tensor, but the molecular viscosity is replaced by the turbulence eddy viscosity $\nu _T$ . By applying the Fourier transform to (2.18) and substituting it into (2.15), we obtain the following modified form of the turbulence curvilinear model equation with the additional eddy viscosity terms:

(2.19) \begin{align} \left (\begin{matrix} (\langle u \rangle -c){\nabla} ^2-\dfrac {\partial ^2\langle u \rangle }{\partial \zeta ^2}-\dfrac {(\nu +\nu _T)}{\mathrm{i} k}{\nabla} ^4\\ -2\dfrac {\nu _T^\prime }{\mathrm{i}k}\big(D^2-k^2 \big)D-\dfrac {\nu _T^{\prime \prime }}{\mathrm{i}k}\big(D^2+k^2 \big) \end{matrix}\right ) \hat {w} & = \begin{matrix} \hat {\eta }(\nu +\nu _{T})\big[g\langle u\rangle ^{\prime \prime }\big]^{\prime \prime }+\hat {\eta }\nu _{T}^{\prime \prime } g\langle u\rangle ^{\prime \prime }\\ +\hat {\eta }\nu _{T}^\prime 2g^\prime \langle u\rangle ^{\prime \prime }+g\hat {\eta }\nu _{T}^\prime \big[2D^2-k^2 \big]\langle u\rangle ^\prime \end{matrix}. \end{align}

A detailed derivation of (2.19) is provided in § 1 of the supplementary materials document.

Interestingly, the differential operator on the left-hand side of (2.19) with the additional eddy viscosity terms is similar to the Orr–Sommerfeld equation derived by Pujals et al. (Reference Pujals, García-Villalba, Cossu and Depardon2009), who extended the classical formulation to turbulent flows by modelling Reynolds stresses using the Boussinesq eddy viscosity parameterisation. Equation (2.19) is subject to the following boundary conditions:

(2.20a) \begin{align} \hat {w}|_{\zeta =0} = \hat {w}_s, \quad & \hat {w}|_{\zeta \rightarrow \infty } = 0, \\[-12pt] \nonumber \end{align}

(2.20b) \begin{align} \frac {\mathrm{d} \hat {w}}{\mathrm{d} \zeta }\bigg |_{\zeta =0} = -\mathrm{i} k \hat {u}_s -\mathrm{i} k \hat {\eta } \bigg ( g \frac {\mathrm{d} U }{\mathrm{d} \zeta } \bigg ) \bigg |_{\zeta =0} , \quad & \frac {\mathrm{d} \hat {w}}{\mathrm{d} \zeta }\bigg |_{\zeta \rightarrow \infty } = 0. \\[9pt] \nonumber \end{align}

Here, $\hat {u}_s$ and $\hat {w}_s$ represent the Fourier transform of the wave orbital velocity prescribed at the bottom wave boundary of the airflow domain. In the current study, we employ (2.2b ) and (2.2d ) for $u_s(x,t)$ and $w_s(x,t)$ , respectively. The corresponding phase-averaged expressions are $\overline {u}_s=\widetilde {u}_s=akc \cos k\xi$ and $\overline {w}_s=\widetilde {w}_s=akc \sin k\xi$ . After the Fourier transform, they become

(2.21) \begin{align} \hat {u}_s&=\frac {1}{2}akc, \ \hat {w}_s=-\mathrm{i}\frac {1}{2}akc. \end{align}

Once $\hat {w}(\zeta )$ is obtained by solving the turbulence curvilinear model (2.19), the corresponding predictions for $\hat {u}(\zeta )$ , $\hat {p}(\zeta )$ and the form drag $F_p$ can be obtained using the expressions derived in §§ 2.2.1, 2.2.2 and 2.2.3. Furthermore, for $\widetilde {\eta } = a \cos k\xi$ , (2.10) shows that $\textrm {Re}\{\hat {f}\}$ and $\textrm {Im}\{\hat {f}\}$ respectively denote the in-phase and out-of-phase components of the wave-induced airflow perturbation relative to the wave surface elevation. In contrast, Cao et al. (Reference Cao, Deng and Shen2020) adopted the definition $\widetilde {\eta } = a \sin k\xi$ , which introduces a phase shift such that $\textrm {Im}\{\hat {f}\}$ and $\textrm {Re}\{\hat {f}\}$ respectively denote the in-phase and out-of-phase components. This distinction should be noted when comparing results between the two formulations, as it affects the interpretation of phase relationships between the airflow and the underlying wave.

In addition to the above discussion, while the present turbulence curvilinear model (2.15) is formulated for monochromatic surface waves, extending it to broadband wave conditions poses additional challenges. In broadband seas, separating wave-coherent motions from background turbulence is non-trivial because of overlapping contributions from long and short waves. A possible extension would be to apply the model in a Fourier sense, solving (2.15) independently for each wavenumber $k$ in the broadband spectrum, similar to the approach of Liu et al. (Reference Liu, Yang, Guo and Shen2010). However, such a formulation neglects nonlinear interactions among wave components, and a key challenge lies in constructing an appropriate representation of the wave-induced turbulence stresses for each spectral mode. Alternatively, studies such as Hristov, Friehe & Miller (Reference Hristov, Friehe and Miller1998) and Lu et al. (Reference Lu, Xiang, Zaki and Katz2025) have employed the Hilbert transform to extract wave-coherent fluctuations from flows over broadband waves. Incorporating such techniques may provide a promising pathway toward extending the present model to broadband sea states and improving the representation of wind–wave interactions, which should be studied in the future.

2.2.1. Expression for wave-induced streamwise perturbation

By applying Fourier transform to the linearised continuity (2.14c ), we obtain wave-induced streamwise velocity perturbation $\hat {u}$ as follows:

(2.22) \begin{align} \hat {u}(\zeta )&=-g\hat {\eta }\frac {\partial \langle u \rangle }{\partial \zeta }-\frac {1}{\mathrm{i}k}\frac {\partial \hat {w}}{\partial \zeta }. \end{align}

The expressions for the real and imaginary parts of $\hat {u}$ are

(2.23) \begin{align} \textrm {Re}\{\hat {u}\}&=-g\hat {\eta }\frac {\partial \langle u \rangle }{\partial \zeta }-\frac {1}{k}\frac {\partial \textrm {Im}\{\hat {w}\}}{\partial \zeta }, \\[-12pt] \nonumber \end{align}

(2.24) \begin{align} \textrm {Im}\{\hat {u}\}&=\frac {1}{k}\frac {\partial \textrm {Re}\{\hat {w}\}}{\partial \zeta }. \\[9pt] \nonumber \end{align}

Equation (2.22) utilises $\hat {w}$ predicted from solving the turbulence curvilinear model (2.15). Finally, using (2.9), we obtain the wave-induced streamwise velocity $\widetilde {u}(\xi ,\zeta )$ as

(2.25) \begin{align} \widetilde {u}(\xi ,\zeta )=2\textrm {Re} [\hat {u}(\zeta )] \cos k\xi -2\textrm {Im} [\hat {u}(\zeta )] \sin k\xi . \end{align}

2.2.2. Expression for wave-induced pressure

Applying Fourier transform to (2.14b ), yields

(2.26) \begin{align} \frac {\partial \hat {p}}{\partial \zeta }&=-( \langle u \rangle -c)\mathrm{i}k\hat {w}+\nu \left (-k^2\hat {w}+\frac {\partial ^2\hat {w}}{\partial \zeta ^2}\right )-\mathrm{i}k\hat {\tau }_{13}-\frac {\partial \hat {\tau }_{33}}{\partial \zeta }. \end{align}

By integrating (2.26) in the vertical direction from an arbitrary location $\zeta$ near the wave surface up to $\lambda$ , which is sufficiently high such that $\tilde {p}(\zeta =\lambda )\approx 0$ and $\hat {\tau }_{33}(\zeta =\lambda )\approx 0$ , we obtain

(2.27) \begin{align} \hat {p}(\zeta )&=\int _{\zeta }^{\lambda } (\langle u \rangle -c)\mathrm{i}k\hat {w} \ \textrm {d}\zeta - \int _{\zeta }^{\lambda } \nu\! \left (-k^2\hat {w}+\frac {\partial ^2\hat {w}}{\partial \zeta ^2}\right ) \ \textrm {d}\zeta +\int _{\zeta }^{\lambda } \mathrm{i} k \hat {\tau }_{13} \ \textrm {d}\zeta -\hat {\tau }_{33}(\zeta ). \end{align}

The $\hat {w}$ predicted from solving (2.15) is used in (2.27) to obtain $\hat {p}(\zeta )$ . From (2.27), the real and imaginary parts of $\hat {p}(\zeta )$ are

(2.28) \begin{align} \textrm {Re}\{\hat {p}\}&=-\int _{\zeta }^{\lambda } (\langle u \rangle -c)k\,\textrm {Im}\{\hat {w}\} \ \textrm {d}\zeta - \int _{\zeta }^{\lambda } \nu\! \left (\frac {\partial ^2}{\partial \zeta ^2}-k^2\right )\textrm {Re}\{\hat {w}\} \ \textrm {d}\zeta \nonumber \\& \quad -\int _{\zeta }^{\lambda } k\, \textrm {Im}\{\hat {\tau }_{13}\} \ \textrm {d}\zeta -\textrm {Re}\{\hat {\tau }_{33}\}(\zeta ), \\[-12pt] \nonumber \end{align}

(2.29) \begin{align} \textrm {Im}\{\hat {p}\}&=\int _{\zeta }^{\lambda } (\langle u \rangle -c)k\,\textrm {Re}\{\hat {w}\} \ \textrm {d}\zeta - \int _{\zeta }^{\lambda } \nu\! \left (\frac {\partial ^2}{\partial \zeta ^2}-k^2\right )\textrm {Im}\{\hat {w}\} \ \textrm {d}\zeta \nonumber \\& \quad +\int _{\zeta }^{\lambda } k\, \textrm {Re}\{\hat {\tau }_{13}\} \ \textrm {d}\zeta -\textrm {Im}\{\hat {\tau }_{33}\}(\zeta ). \\[9pt] \nonumber \end{align}

By using (2.28) and (2.29) in (2.9), we can obtain the wave-induced pressure $\widetilde {p}(\xi ,\zeta )$ as

(2.30) \begin{align} \widetilde {p}(\xi ,\zeta )=2\textrm {Re}[\hat {p}(\zeta )] \cos k\xi -2\textrm {Im}[\hat {p}(\zeta )] \sin k\xi . \end{align}

2.2.3. Form drag

The form drag is defined as the integral of the wave-induced pressure on the wave surface ( $\widetilde {p}(\xi ^\prime ,\zeta =0)$ ) over the wavelength $\lambda =2\pi /k$ (Sullivan et al. Reference Sullivan, McWilliams and Moeng2000; Yang & Shen Reference Yang and Shen2010) and is given by

(2.31) \begin{align} F_p=\frac {1}{\lambda }\int _{0}^{\lambda }\frac {\widetilde {p}(\xi ^\prime ,\zeta =0)}{u_\tau ^2}\frac {\textrm {d}\widetilde {\eta }}{\textrm {d}\xi ^\prime } \ \textrm {d}\xi ^\prime . \end{align}

Note that $\widetilde {p}$ in (2.31) is the modified pressure, which incorporates the trace of the wave-induced turbulence stresses. Using $\widetilde {p}=2\textrm {Re}\{\hat {p}\}\cos (k\xi )-2\textrm {Im}\{\hat {p}\}\sin (k\xi )$ from (2.30) and $\widetilde {\eta }=a\cos (k\xi )$ , which is the phase-averaged expression for wave elevation from (2.2a ), we obtain

(2.32) \begin{align} F_p&=\frac {1}{ u_\tau ^2\lambda }\int _{0}^{\lambda }(2\textrm {Re}\{\hat {p}(\zeta =0)\}\cos (k\xi ^\prime )-2\textrm {Im}\{\hat {p}(\zeta =0)\}\sin (k\xi ^\prime ))(-ak \sin (k\xi ^\prime )) \ \textrm {d}\xi ^\prime , \nonumber \\ &=\frac {-ak}{u_\tau ^2\lambda }\Biggl [\int _{0}^{\lambda }2\textrm {Re}\{\hat {p}(\zeta =0)\}\cos (k\xi ^\prime )\sin (k\xi ^\prime )\, \textrm {d}\xi ^\prime \nonumber \\ &\mspace {60mu}-\int _{0}^{\lambda }2\textrm {Im}\{\hat {p}(\zeta =0)\}\sin (k\xi ^\prime )\sin (k\xi ^\prime ) \ \textrm {d}\xi ^\prime \Biggr ],\nonumber \\ &=\frac {ak\,\textrm {Im}\{\hat {p}(\zeta =0)\}}{u_\tau ^2}. \end{align}

Equation (2.32) indicates that the form drag arises from the out-of-phase component of $\widetilde {p}$ relative to the cosine wave surface. Substituting (2.29) into (2.32) yields

(2.33) \begin{align} F_p&=F_p^{\textit{ad}v}+F_p^{\nu }+F_p^{\textit{turb}}, \end{align}

where $F_p^{\textit{ad}v}$ is the form drag caused by the linear advection term, $F_p^{\nu }$ is the viscous drag component and $F_p^{\textit{turb}}$ is the turbulent drag component. They are expressed as

(2.34) \begin{align} F_p^{\textit{ad}v}&=\frac {ak}{ u_\tau ^2}\int _{0}^{\lambda } (\langle u \rangle -c) k\, \textrm {Re}\{\hat {w}\}\,\textrm {d}\zeta , \\[-12pt] \nonumber \end{align}

(2.35) \begin{align} F_p^{\nu }&=-\frac {ak}{ u_\tau ^2}\int _{0}^{\lambda }\nu\! \left (\frac {\partial ^2}{\partial \zeta ^2}-k^2\right )\textrm {Im}\{\hat {w}\} \ \textrm {d}\zeta , \\[-12pt] \nonumber \end{align}

(2.36) \begin{align} F_p^{\textit{turb}}&=\frac {ak}{u_\tau ^2}\left [\int _{0}^{\lambda }k\,\textrm {Re}\{\hat {\tau }_{13}\}\textrm {d}\zeta -\textrm {Re}\{\hat {\tau }_{33}(\zeta =0)\}\right ]\!. \\[9pt] \nonumber \end{align}

To summarise § 2, we have derived the turbulence curvilinear model from the governing wave-induced momentum equations. Owing to its reduced-order form, solving this model enables rapid predictions of wave-induced airflow perturbations and the associated form drag.

3.1. Numerical model

The LES solver utilised in this study to simulate turbulent airflows over water waves is based on the method developed by Yang & Shen (Reference Yang and Shen2011) and Xuan & Shen (Reference Xuan and Shen2019). The governing equations are the non-dimensional, filtered form of the NS (2.1), which are transformed into a boundary-fitted curvilinear coordinate system using the transformation defined in (2.3). The non-dimensionalisation uses $u_\tau$ and $\delta$ as the characteristic velocity and length scales, respectively. This normalisation results in the multiplication of the $1/\textit{Re}_\tau$ factor with the divergence of the viscous stress term, where $ \textit{Re}_\tau =u_\tau \delta /\nu$ is the friction Reynolds number.

The LES equations solve for the filtered velocity $u_i$ and pressure $p$ corresponding to a filter length scale $\varDelta =\sqrt [3]{\Delta x \Delta y \Delta z}$ . To model the effects of subgrid-scale (SGS) motions, a divergence of the SGS stress tensor $(\sigma ^{\textit{SGS}}_{\textit{ij}})$ term is added to the momentum equations. The deviatoric trace-free part of the SGS stress tensor, $\sigma ^{d,\textit{SGS}}_{\textit{ij}}=\sigma ^{\textit{SGS}}_{\textit{ij}}-\delta _{\textit{ij}} \sigma ^{\textit{SGS}}_{kk}{/3}$ , is evaluated as $\sigma _{\textit{ij}}^{d,\textit{SGS}}=-2\nu _{T,\textit{SGS}}S_{\textit{ij}}$ , where the SGS eddy viscosity $\nu _{T,\textit{SGS}}$ is obtained from the dynamic Smagorinsky SGS model (Smagorinsky Reference Smagorinsky1963; Germano et al. Reference Germano, Piomelli, Moin and Cabot1991; Lilly Reference Lilly1992).

To simulate the LES equations, a hybrid finite-difference and pseudo-spectral method accelerated by graphics processing units (Xuan, Ren & Shen Reference Xuan, Ren and Shen2025) is used for spatial discretisation. The pseudo-spectral method based on the fast Fourier transform is used for the horizontal directions $\xi$ and $\psi$ . A second-order finite-difference scheme is used in the vertical direction $\zeta$ with refinement near the top and bottom boundaries. A staggered grid is imposed in the vertical direction to avoid the checkerboard oscillations, where $u,v$ and $p$ are stored at the cell centre and $w$ is located at the cell face. The fractional-step method (Kim & Moin Reference Kim and Moin1985) is implemented, in which a predictor step is evaluated first by advancing the velocity field over a time step according to the momentum equations but omitting the pressure term. The time advancement uses the second-order Adams–Bashforth scheme. Following this predictor step, the pressure Poisson equation is solved by enforcing the incompressibility condition. Owing to the grid transformation, the resulting pressure Poisson equation contains nonlinear terms, and is iteratively solved to obtain the pressure field (Yang & Shen Reference Yang and Shen2011). Finally, we obtain the divergence-free velocity through a correction step using the pressure field obtained from the pressure Poisson equation. This solver has been extensively tested and validated to simulate a variety of wind over wave problems, including monochromatic waves (Yang & Shen Reference Yang and Shen2010, Reference Yang and Shen2017; Cao et al. Reference Cao, Deng and Shen2020; Cao & Shen Reference Cao and Shen2021) and broadband waves (Yang et al. Reference Yang, Meneveau and Shen2013, Reference Yang, Meneveau and Shen2014; Hao & Shen Reference Hao and Shen2019). In the following section, we describe the simulation set-up used in this study to generate LES data for various wind over wave scenarios. The LES data are used to understand the wind over wave flow physics, validate the formulation of the turbulence curvilinear model and test the validity of the eddy viscosity parameterisations for wave-induced turbulence stress components.

3.2. Simulation set-up

The length, width and height of the computational domain (figure 1) are set to $(L_\xi ,\, L_\psi ,\, L_\zeta )$ $=(8\pi\! H,\, 4\pi\! H,\,H)$ . Here, $H$ represents the mean height of the domain, which is also taken as the characteristic length scale denoted as $\delta =H$ . The friction Reynolds number $ \textit{Re}_\tau =u_\tau H/\nu$ for the simulations is set to 1000. The grid size is $(N_\xi \times N_\psi \times N_\zeta ) = (864 \times 576 \times 144)$ , which provides a resolution of $\Delta \xi ^+ \approx 30$ , $\Delta \psi ^+ \approx 22$ and $\Delta \zeta _{\textit{min}}^+\approx 0.27$ . Here, $\Delta \zeta _{\textit{min}}$ denotes the minimum grid spacing near the boundary in the $\zeta$ direction. The superscript ‘+’ denotes normalisation in the inner units, for example, $\zeta ^+=\zeta u_\tau /\nu =(\zeta /H)\textit{Re}_\tau , u_i^+=u_i/u_\tau , t^+=tu_\tau ^2/\nu$ . This grid resolution satisfies the requirements for resolving the near-wall flow features in wall-resolved LES (WRLES) according to Choi & Moin (Reference Choi and Moin2012).

We first conduct precursor WRLES of a turbulent half-channel flow with a flat bottom boundary. Once the precursor simulation reaches a statistically stationary state, a monochromatic surface wave with a prescribed wave age ( $c/u_\tau$ ) is introduced at the bottom boundary of the simulation domain. Various wind over wave scenarios are simulated by specifying different wave-age values, given by $c/u_\tau = 0, \pm 2, \pm 7, \pm 15, \pm 25$ . The ‘ $+$ ’ and ‘ $-$ ’ signs of the wave age represent wind-following and wind-opposing wave cases, respectively. Since we utilise a boundary-fitted grid method that cannot capture wave breaking, we limit the wave steepness $(ak)$ by setting its value to 0.1. According to Stokes (Reference Stokes1847), the theoretical breaking limit is $ak = 0.443$ , whereas experimental and numerical studies (Longuet-Higgins & Cokelet Reference Longuet-Higgins and Cokelet1976; Duncan Reference Duncan1981; Dommermuth & Yue Reference Dommermuth and Yue1987; Iafrati Reference Iafrati2009; Perlin, Choi & Tian Reference Perlin, Choi and Tian2013) suggest $ak \approx 0.3$ as a more realistic threshold. Furthermore, Cao et al. (Reference Cao, Deng and Shen2020) demonstrated that, for a given wave age, variations in wave steepness within the non-breaking range primarily scale the magnitude of the wave-induced perturbations without significantly altering their vertical structure. For this reason, the present study varies the wave age while keeping $ak$ fixed at 0.1. We set the wavelength to $\lambda = (\pi /2)H$ , corresponding to a wavenumber $k = 4/H$ . Consequently, the simulation domain of length $L_\xi = 8\pi\! H$ accommodates sixteen wavelengths. Once each wind over wave simulation reaches a statistically steady state, more than $300$ time snapshots that span 250 viscous time units ( $tu_{\tau }^2/\nu$ ) are stored to calculate the turbulence statistics. A balanced mean momentum budget, where the total stress for each wind over wave case follows the linear relation $ u_\tau ^2(1 - \zeta /H)$ , as shown in § 3 of the supplementary material, is consistent with the LES reaching a statistically steady state.

We note that, although the present study employs WRLES at a finite Reynolds number, the conclusions drawn in this study on wave-induced turbulence remain applicable to high-Reynolds-number atmospheric boundary layer flows that are typically simulated using wall-modelled LES (WMLES) approach. A WMLES framework can capture wave-coherent motions if the wave phase is explicitly resolved using a boundary-fitted grid. This has been demonstrated in studies by Sullivan et al. (Reference Sullivan, McWilliams and Patton2014) and Hara & Sullivan (Reference Hara and Sullivan2015), who reported the presence of wave-coherent structures above the surface in wave-phase-resolved WMLES set-ups. Alternatively, recent studies such as Aiyer, Deike & Mueller (Reference Aiyer, Deike and Mueller2023) and Ayala et al. (Reference Ayala, Sadek, Ferčák, Cal, Gayme and Meneveau2024) have proposed WMLES frameworks where the waves are resolved only in the horizontal direction but not vertically. These approaches incorporate wave effects through a wave-induced wall shear stress correction into the equilibrium log-law stress. While this simplified formulation captures the mean-flow variations associated with wave-induced form drag, the lack of vertical wave resolution limits its ability to represent the modulation of turbulence and the excitation of wave-coherent motions. Further refinement of wall models is therefore necessary to provide a more complete description of the coupled wave–turbulence dynamics.

Because the turbulence curvilinear model (2.15) relies on the mean velocity profile and wave-induced turbulence stresses, understanding how these quantities vary over the different wave ages is essential. Therefore, we first analyse the turbulence statistics of the mean velocity profile in § 3.3, and then discuss the properties of the wave-induced perturbations and the corresponding turbulence stresses in § 3.4.

3.3. Mean velocity profile for wind over wave flows

In this section, we discuss the properties of the mean velocity profile across various wave ages. The progressive waves induce varying levels of drag on the airflow field, thus altering the magnitude of the mean velocity depending on the wave age. As a validation of the simulation set-up, figure 2(a) shows the mean velocity profile over a flat wall, which exhibits the expected logarithmic region for $\zeta ^+\gt 30$ and a linear viscous sublayer for $\zeta ^+\lt 5$ . Building on this baseline, figure 2(b) highlights the modifications to the mean-flow profile induced by progressive waves at various wave ages ( $c/u_\tau = 0, \pm 2, \pm 25$ ). Qualitatively, all the profiles have similar characteristics, where there exist a linear viscous sublayer region near the wave surface and a logarithmic region at larger heights ( $\zeta ^+\gt 30$ ). Qualitatively, for the wind-following wave scenarios, as the wave age increases, the mean velocity magnitude initially decreases and then increases, approaching the flat-wall scenario at the wave age $c/u_\tau =25$ . Furthermore, the slope of the logarithmic region increases with increasing wave age, in accordance with the findings of previous studies, such as Sullivan et al. (Reference Sullivan, McWilliams and Moeng2000), Sullivan et al. (Reference Sullivan, Edson, Hristov and McWilliams2008), Yang & Shen (Reference Yang and Shen2010) and Sullivan et al. (Reference Sullivan, McWilliams and Patton2014). This increase in mean velocity with wave age can be attributed to the low wave-induced drag under wind-following wave conditions. On the other hand, for the wind-opposing wave scenarios, the mean velocity magnitude tends to decrease as the wave-age magnitude increases, and the slope of the logarithmic region also decreases with increasing wave age, which indicates that a greater drag is induced by the wind-opposing waves. Although omitted from the figure for clarity, the intermediate cases ( $c^+ = \pm 7, \pm 15$ ) fit these observed trends: For wind-following waves, the $c^+ = 7$ profile shows a greater velocity deficit than the $c^+=25$ profile, while the $c^+=15$ profile is nearly indistinguishable from the flat-wall profile. Conversely, for wind-opposing waves, the $c^+=-7$ and $c^+=-15$ profiles show a progressively larger velocity reduction than the $c^+=-2$ profile, confirming that the greatest drag is induced by the largest magnitude of opposing wave age.

A detailed discussion on the behaviour of the form drag is provided in § 3.5.

Figure 2. Mean streamwise velocity profile $\langle u\rangle ^+=\langle u\rangle /u_\tau$ for (a) a flat wall and its comparison with (b) wind over wave cases.

3.4. Characteristics of wave-induced perturbations

In this section, we examine the properties of wave-induced turbulence. Section 3.4.1 presents contour visualisations of the wave-induced airflow perturbations and discusses their dependence on the wave age. The spectral characteristics of the wave-induced Reynolds stresses are analysed in § 3.4.2.

3.4.1. Contours of wave-induced airflow perturbations

Waves propagating at the water surface induce disturbances to the airflow. This section explores these wave-induced perturbations for the various wave ages. Wave-coherent motions are identified from the LES results by subtracting the mean flow from the phase-averaged flow data (§ 2.1). We find that for the chosen wave steepness, $a k = 0.1$ , the fundamental wave-coherent mode accounts for more than 90 % of the total wave-induced kinetic energy across all the wave ages considered. This dominant contribution supports the validity of representing the wave-coherent motions using only the fundamental mode in (2.9), while neglecting higher-order components. The profound effects of wave age on the wave-induced streamwise velocity ( $\widetilde {u}$ ), vertical velocity ( $\widetilde {w}$ ) and pressure ( $\widetilde {p}$ ) for two distinct wave-age regimes, namely, fast-moving waves ( $c/u_\tau = \pm 25$ ) and slow-moving waves ( $c/u_\tau = \pm 2$ ), are displayed in figure 3. The flow fields exhibit a characteristic undulating pattern, with the influence of the wave predominantly localised near the surface and the perturbation amplitude diminishing rapidly in the vertical direction.

The intensity of the wave-induced airflow perturbations are indicated by the contour colour bars, which reveal a clear dependence on the wave age. The fast-moving, wind-opposing wave ( $c/u_\tau =-25$ ) generates the most pronounced flow perturbations, which exceed those induced by the slow-moving wave ( $c/u_\tau =\pm 2$ ) and the fast-moving, wind-following wave ( $c/u_\tau =25$ ). This phenomenon is intrinsically linked to wave forcing through both the wave phase speed and orbital velocity, which induce streamwise and vertical pressure gradients that drive the observed secondary flow motions. The counter-propagating wave with $c/u_\tau = -25$ excites stronger perturbations than the other wave-age cases do. These results demonstrate how the wave speed and direction relative to the wind directly influence the strength of wave-induced secondary flows through the modulation of the pressure field.

The observed dependence on wave age holds true for the intermediate cases, $c^+ = \pm 7$ and $c^+ = \pm 15$ , which are omitted from the contour plots for visual clarity. The wind-opposing cases ( $c^+ = -7$ and $c^+ = -15$ ) have stronger wave-induced perturbations ( $\widetilde {u}, \widetilde {w}, \widetilde {p}$ ) than their wind-following wave counterparts ( $c^+ = 7$ and $c^+ = 15$ ). Specifically, the strength of the perturbation fields intensifies as the magnitude of opposing wave age increases from $c^+ = -7$ to $c^+ = -15$ , confirming the link between stronger opposing wave motion and enhanced secondary motions. It is also worth noting that the $c^+=15$ case exhibits the weakest overall perturbation fields, whereas the magnitudes of perturbations are slightly enhanced for the faster $c^+=25$ case. This increase in perturbation strength for the $c^+=25$ case compared with the $c^+=15$ case indicates that for large wind-following wave ages, the increased wave phase speed dominates the perturbation magnitude due to enhanced orbital velocities.

Figure 3. Contours of wave-induced (a,d,g, j) streamwise velocity $\widetilde {u}/u_\tau$ , (b,e,h,k) vertical velocity $\widetilde {w}/u_\tau$ and (e, f,i,l) pressure $\widetilde {p}/u_\tau ^2$ from LES of turbulent flow over progressive waves with wave ages $c/u_\tau =\pm 2, \pm 25$ .

3.4.2. Spectral behaviour of wave-induced turbulence

Continuing the discussion on wave-induced flow perturbations, in this section, we discuss the properties of the wave-induced Reynolds stress components $\widetilde {\tau }_{\textit{ij}}$ computed from the LES data. The Reynolds stress components are defined as $\tau _{\textit{ij}}=\overline {u^\prime _i U^\prime _j}+\overline { \tau }_{\textit{ij}}^{d,\textit{SGS}}$ , where $\overline {u_i' U_{\!j}'}$ and $\overline {\tau }_{\textit{ij}}^{d,\textit{SGS}}$ represent the Reynolds stress contributions from the resolved and unresolved (SGS) flow fields, respectively. The resolved component is evaluated as $\overline {u_i' U_{\!j}'} = \overline {u_i U_{\!j}} - \overline {u}_i \overline {U}_{\!j}$ , where $\overline {(\boldsymbol{\cdot })}$ denotes the phase-averaging operation defined in (2.7). The SGS contribution $\overline {\tau }_{\textit{ij}}^{d,\textit{SGS}}$ corresponds to the stress tensor in the curvilinear coordinate system, which is given by $\overline {\tau }_{\textit{ij}}^{d,\textit{SGS}} = J^{-1} \overline {\sigma }^{d,\textit{SGS}}_{\textit{il}} \partial \xi _j / \partial x_l$ , where $\overline {\sigma }^{d,\textit{SGS}}_{\textit{il}} = -2 \nu _{T,\textit{SGS}} \overline {S}_{\textit{il}}$ is the trace-free part of the SGS stress tensor introduced in § 3.1. Finally, upon subtracting the streamwise-averaged stress $\langle \tau _{\textit{ij}}\rangle$ from the phase-averaged stress $\tau _{\textit{ij}}$ , we obtain the wave-induced stress, $\widetilde {\tau }_{\textit{ij}}=\tau _{\textit{ij}}-\langle \tau _{\textit{ij}}\rangle$ . Because the monochromatic wave is aligned with the wind direction in this study, the $\widetilde {\tau }_{12},\widetilde {\tau }_{21},\widetilde {\tau }_{23}$ and $\widetilde {\tau }_{32}$ stress components have a mean value of zero by definition. Therefore, only the stress components $\widetilde {\tau }_{11}, \widetilde {\tau }_{22}, \widetilde {\tau }_{33}, \widetilde {\tau }_{13}$ and $\widetilde {\tau }_{31}$ are discussed here. The gradients of these wave-induced turbulence stress components appear as a forcing term, in addition to the forcing caused by the wave elevation, in the turbulence curvilinear model (2.15).

Figure 4. Vertical profiles of the real (a,b,c,d, e) and imaginary ( f,g,h,i, j) components, magnitude (k,l,m,n,o) and the phase difference relative to wave surface (2.11) (p,q,r,s,t) of the $\hat {\tau }^+_{11}$ (a, f,k,p), $\hat {\tau }^+_{22}$ (b,g,l,q), $\hat {\tau }^+_{33}$ (c,h,m,r), $\hat {\tau }^+_{13}$ (d,i,n,s) and $\hat {\tau }^+_{31}$ (e, j,o,t) wave-induced Reynolds stresses from LES for the $c/u_\tau =-25$ (), $c/u_\tau =-2$ (), $c/u_\tau =0$ (), $c/u_\tau =2$ () and $c/u_\tau =25$ () cases.

The fundamental mode shapes for $\widetilde {\tau }^+_{11}=\widetilde {\tau }_{11}/u_\tau ^2, \widetilde {\tau }^+_{22}=\widetilde {\tau }_{22}/u_\tau ^2,\widetilde {\tau }^+_{33}=\widetilde {\tau }_{33}/u_\tau ^2,\widetilde {\tau }^+_{13}=\widetilde {\tau }_{13}/u_\tau ^2$ and $\widetilde {\tau }^+_{31}=\widetilde {\tau }_{31}/u_\tau ^2$ are shown in figure 4(a– j) for wave ages $c/u_\tau =0, \pm 2, \pm 25$ . The intermediate wave-age cases ( $c/u_\tau = \pm 7, \pm 15$ ) are not shown for brevity; however, their behaviours are qualitatively consistent with the adjacent cases and are briefly discussed in the text. The real parts of the stress components in figure 4(a–e) represent the in-phase (relative to the wave surface elevation) contribution to the momentum transfer, whereas the imaginary parts of the stresses in figure 4( f– j) represent the out-of-phase contribution. The magnitudes in figure 4(k–o) indicate the overall strength of the stress components, and the phase differences in figure 4(p–t) reveal whether the turbulence stress leads or lags the wave surface. The profiles are plotted against the non-dimensional vertical coordinate $\zeta ^+$ .

Notably, the wave age influences the wave-induced Reynolds stress profiles significantly. For example, the real and imaginary parts of all the stress components $\hat {\tau }^+_{11}$ (figures 4 a and 4 f), $\hat {\tau }^+_{22}$ (figures 4 b and 4 g), $\hat {\tau }^+_{33}$ (figures 4 c and 4 h), $\hat {\tau }^+_{13}$ (figures 4 d and 4 i) and $\hat {\tau }^+_{31}$ (figures 4 e and 4 j) show distinct variations with the wave age, thus indicating the strong dependence of momentum transfer on the relative speed between the wave and the wind. The magnitudes of the stress components in figure 4(k–o) also reveal that the overall strength of the wave-induced turbulence stresses varies with the wave age.

A key characteristic of the wave-induced Reynolds stresses is that the magnitudes of the normal stress components are greater than those of the shear stress components, which, along with the differences among the various normal stress components, reflects the anisotropic nature of the wave-coherent turbulence field. The profiles also illustrate the vertical structure of the wave-induced turbulence stresses. For the wave-age cases of $c^+ = -25$ to $c^+ = 7$ , the components, particularly the normal stress components $|\hat {\tau }_{11}|$ (figure 4 k) and $|\hat {\tau }_{33}|$ (figure 4 m), exhibit multiple distinct peaks. The presence of these multiple maxima indicates the existence of different layers, such as an inner layer (with a peak typically located just below $\zeta ^+ \sim 10$ ) and an outer layer (with a second, stronger peak located just above $\zeta ^+ \sim 10$ ) (Belcher & Hunt Reference Belcher and Hunt1993) within the wave boundary layer. In all the panels, the magnitudes of the stresses decay at larger heights due to the rapid advection of turbulence by the mean flow (Belcher & Hunt Reference Belcher and Hunt1993; Mastenbroek et al. Reference Mastenbroek, Makin, Garat and Giovanangeli1996). The phase differences in figure 4(p–t) suggest that the stress fluctuations are not always in-phase with the wave motion, which indicates the complexity of phase shifting in wave-induced turbulent flows. The rapid change in the phase angle for the $\hat {\tau }_{11}$ and $\hat {\tau }_{33}$ components at $\zeta ^+\sim 10$ also represents the change in sign of the real and imaginary components as shown in figure 4(a,c,f,h).

Furthermore, the variation in the magnitudes of the stresses is highly complex and does not exhibit a monotonic trend across the different components. The magnitudes of the normal stress components $|\hat {\tau }_{11}|$ , $|\hat {\tau }_{22}|$ and $|\hat {\tau }_{33}|$ (figure 4 k–m) increase monotonically from $c^+ = -25$ to $c^+ = 2$ and then decrease for the $c^+ = 7, 15, 25$ cases. In contrast, the shear stress component $|\hat {\tau }_{13}|$ (figure 4 n) remains at comparable magnitude for $c^+ = -25$ to $c^+ = 0$ , but its magnitude diminishes for the $c^+ = 7, 15, 25$ cases. The magnitude of $|\hat {\tau }_{31}|$ shows a different behaviour: it increases monotonically from $c^+ = -25$ to $c^+ = 2$ , and remains comparable to the $c^+ = 0$ case for $c^+ = 7, 15, 25$ .

To summarise § 3.4, the wave age significantly influences the wave-induced turbulence stress fields and affects the momentum transfer in the wave boundary layer. This result underscores the importance of considering wave-age effects in modelling turbulent flows over moving wavy surfaces, which is discussed in §§ 4, 5 and 6.

Figure 5. Variation of wave-induced pressure normalised by its peak value, $\widetilde {p}/|\widetilde {p}|_{\textit{max}}$ , for wave ages (a) $c/u_\tau =-25$ (), (b) $c/u_\tau =-2$ (), (c) $c/u_\tau =2$ (), (d) $c/u_\tau =25$ () over the $\widetilde {\eta }=a\cos (k\xi )$ wave surface (). Panel (e) shows the variation of $|\widetilde {p}|_{\textit{max}}/u_\tau ^2$ and phase difference $\phi _{\widetilde {p}\,\widetilde {\eta }}(\zeta =0)$ (deg.) across the various wave ages, and panel ( f) shows the equivalence between the form drag $F_p$ estimated from (2.32) and the mean wave-coherent pressure stress $-\langle \tau _{13}^p\rangle (\zeta =0)$ at the wave surface.

3.5. Form drag

In this section, we discuss the variation of form drag across the various wave-age cases. Physically, the form drag ( $F_p$ ) defined in (2.31) can be interpreted as the correlation between the wave-induced pressure and the slope of the wave surface, $\langle \widetilde {p}\,\widetilde {\eta }_\xi \rangle$ (Miles Reference Miles1957). This indicates that the form drag arises from the phase difference ( $\phi _{\widetilde {p}\,\widetilde {\eta }}$ ) between the pressure distribution on the wave surface and the wave elevation. To illustrate this phase difference, figure 5(a–d) presents the wave-induced pressure distribution, normalised by its peak magnitude $|\widetilde {p}|_{\textit{max}}$ , which is clearly out of phase relative to the wave elevation $\widetilde {\eta } = a \cos (k\xi )$ across the different wave ages. The actual values of the phase difference, $\phi _{\widetilde {p}\,\widetilde {\eta }}(\zeta =0)$ , calculated from (2.11) for different wave ages are shown in figure 5(e). The phase difference $\phi _{\widetilde {p}\,\widetilde {\eta }}(\zeta =0)$ varies between $-55^\circ$ and $-90^\circ$ . Consequently, the out-of-phase component $\textrm {Im}\{\hat {p}(\zeta =0)\}$ contributes to the form drag, as expressed in (2.32).

The dependence of $|\widetilde {p}|_{\textit{max}}$ on the wave age is shown in figure 5(e). The case $c/u_\tau =-25$ excites the strongest $|\widetilde {p}|_{\textit{max}}$ . As $c/u_\tau$ increases to 25, the pressure perturbation diminishes to nearly zero. As a result, the variation in the strength of $|\widetilde {p}|_{\textit{max}}$ gives rise to differences in the form drag $(F_p)$ (2.32) across the different wave ages, as shown in figure 5( f). The form drag is highest for the $c/u_\tau =-25$ case and decreases for the $c/u_\tau =-15$ and $-7$ cases as $c/u_\tau$ increases. It then increases from $c/u_\tau =-2$ , reaching a maximum in the $c/u_\tau =2$ case, before declining again from $c/u_\tau =7$ to $25$ . In addition, because (2.31) represents a $\langle \widetilde {p} \eta _\xi \rangle$ correlation, the mean of wave-coherent pressure stress at the wave surface $-\langle \tau _{13}^p\rangle (\zeta =0)$ , where $\tau _{\textit{ij}}^p=J^{-1} p\partial \xi _j/\partial x_i$ , is equal to the form drag at the wave surface, as illustrated in figure 5( f).

The above variation in form drag across different wave ages aligns with the behaviour of the mean-flow profiles discussed in § 3.3. In particular, in the $c/u_\tau =-25$ case, which has the highest form drag, the mean-flow velocity is significantly reduced. In contrast, the mean-flow profile for $c/u_\tau =25$ closely resembles that of the flat-wall case, owing to the negligible form drag.

In summary, § 3 describes the LES set-up for wind over waves spanning a range of wave ages. The LES results reveal that progressive waves generate airflow perturbations near the wave surface, with their intensity varying across different wave ages. These perturbations impart a form drag on the overlying airflow, and the resulting deviations in the mean velocity profile relative to the flat-wall case depend on the wave age. In the following section, § 4, we assess the validity of the turbulence curvilinear model in accurately describing wave-induced perturbations by solving (2.15), using the LES-derived mean velocity and Reynolds stresses as inputs.

4.1. Comparison of models of $\textrm {Re}{\{\hat {w}\}}$ and $\textrm {Im}{\{\hat {w}\}}$ with LES results

In this section, we discuss the properties of $\textrm {Re}\{{\hat {w}}\}$ and $\textrm {Im}\{{\hat {w}}\}$ , which are the modes of the wave-induced vertical airflow disturbances, and compare the viscous and turbulence curvilinear models with LES data. The quantity Re{ $\hat {w}$ } represents the in-phase (relative to the cosine wave surface profile) component of the vertical velocity perturbation induced by the wave. The value of $\textrm {Re}\{\hat {w}\}$ obtained from solving the turbulence curvilinear model (2.15) is equivalent to solving (2.16). In contrast, the viscous curvilinear model neglects the forcing due to turbulence stress terms in (2.15), resulting in accounting for only the forcing contributions from the wave elevation ( $\hat {\eta }$ ) and $\textrm {Im}\{\hat {w}\}$ terms in (2.16).

The effect of retaining the turbulence stress forcing terms is illustrated in figure 6(a–i), which shows the vertical profile of Re{ $\hat {w}$ } for a range of wave ages and compares the viscous and turbulence curvilinear models with the LES data. For the fast-propagating, wind-opposing waves in figure 6(a–c) that correspond to the $c/u_\tau =-25$ , $-15$ and $-7$ cases, both the viscous and turbulence curvilinear models capture the correct trend of the vertical variation in Re{ $\hat {w}$ }. However, for the $c/u_\tau =-7$ case, the turbulence curvilinear model agrees with the LES data much better than the viscous curvilinear model does, in that the former captures the peak magnitude of $\textrm {Re}\{\hat {w}\}$ correctly. This comparison indicates that as the magnitude of the wind-opposing wave age decreases, turbulence stresses, which are neglected in the viscous model, play a crucial role in determining the shape of the vertical velocity profile.

Figure 6. Vertical profiles of $\textrm {Re}\{\hat {w}\}/u_\tau$ obtained from the viscous () and turbulence () curvilinear models compared with the LES data () for various wave ages: (a,i) $c/u_\tau =\pm 25$ , (b,h) $c/u_\tau =\pm 15$ , (c,g) $c/u_\tau =\pm 7$ , (d, f) $c/u_\tau =\pm 2$ , (e) $c/u_\tau =0$ .

As the wave age increases further in figure 6(d– f), which corresponds to the cases of $c/u_\tau =-2$ , $0$ (stationary wavy wall) and $2$ , the discrepancy between the viscous curvilinear model and the LES data becomes even larger. The turbulence curvilinear model, however, is in better agreement with the LES data, which highlights the increasing importance of wave-induced turbulence stresses at small magnitudes of wave age. The peak magnitude of Re{ $\widetilde {w}$ } and its vertical decay with height are captured by the turbulence curvilinear model more accurately, which suggests a better representation of the physical processes involved.

For the wind-following wave cases of $c/u_\tau =7,15$ and 25 in figure 6(g–i), both models predict similar peak values. Nevertheless, the turbulence curvilinear model provides better agreement with the LES data than does the viscous model. The above results suggest that when the wave-induced turbulence stresses discussed in § 3.4.2 are incorporated into the turbulence curvilinear modelling framework (2.16), the resulting description of the $\textrm {Re}\{\hat {w}\}$ mode shapes is improved.

Figure 7. Vertical profiles of $\textrm {Im}\{\hat {w}\}/u_\tau$ obtained from the viscous () and turbulence () curvilinear model compared with the LES data () for various wave ages: (a,i) $c/u_\tau =\pm 25$ , (b,h) $c/u_\tau =\pm 15$ , (c,g) $c/u_\tau =\pm 7$ , (d, f) $c/u_\tau =\pm 2$ , (e) $c/u_\tau =0$ .

Next, we examine the vertical profiles of the imaginary component of the wave-induced vertical velocity (Im{ $\hat {w}$ }) displayed in figure 7 across a range of wave ages, and compare the viscous and turbulence curvilinear models with the LES data. From the triple-decomposition expression (2.6), Im{ $\hat {w}$ } represents the out-of-phase (with respect to the wave surface elevation $\eta$ ) component of the vertical velocity perturbation induced by the wave. The magnitude of $\textrm {Im}\{\hat {w}\}$ is greater than that of $\textrm {Re}\{\hat {w}\}$ , because the airflow vertical motion is mostly induced by the wave surface vertical motion, which has a $90^\circ$ phase shift relative to the wave surface elevation due to the free surface kinematic boundary condition. This phase shift relative to the wave surface is qualitatively illustrated by the $\widetilde {w}$ contours plotted in figure 3.

Notably, Re{ $\hat {w}$ } and Im{ $\hat {w}$ } are coupled. The linearised governing equation for $\textrm {Im}\{\hat {w}\}$ is given by (2.17), where a viscous effect on $\textrm {Im}\{\hat {w}\}$ arises through the term $-(\nu /k) {\nabla} ^4 \textrm {Re}\{\hat {w}\}$ , i.e. through $\textrm {Re}\{\hat {w}\}$ . Interestingly, owing to the coupling between the equations for $\textrm {Re}\{\hat {w}\}$ (2.16) and $\textrm {Im}\{\hat {w}\}$ (2.17), the forcing from the wave elevation and turbulence stresses exerts a coupled influence on $\textrm {Re}\{\hat {w}\}$ and $\textrm {Im}\{\hat {w}\}$ . In the context of turbulence modelling for wave-induced turbulence stresses, this phenomenon indicates the necessity of capturing both the magnitude and phase of the wave-induced turbulence stresses to accurately predict the characteristics of $\hat {w}$ . In figure 7(a–i), the $\textrm {Im}\{\hat {w}\}$ predictions from both the viscous and turbulence curvilinear models are similar and compare well with the LES results, except for the slight under-prediction of the peak magnitude by the viscous model for the $c/u_\tau =2, 7$ and $15$ cases. For the other wave-age cases, both models capture the peak magnitude as well as the decay of $\textrm {Im}\{\hat {w}\}$ at larger heights, in close agreement with the LES data. This comparison suggests that the viscous terms dominate the turbulence stress terms in (2.15) for predicting $\textrm {Im}\{\hat {w}\}$ . Nevertheless, although the turbulence stress terms may seem less significant for predicting $\textrm {Im}\{\hat {w}\}$ in the turbulence curvilinear model, the coupling between (2.16) and (2.17) highlights the fact that the turbulence stress terms in (2.16) for $\textrm {Re}\{\hat {w}\}$ indirectly influence the properties of $\textrm {Im}\{\hat {w}\}$ in (2.17). Consequently, the turbulence curvilinear model provides a better prediction of wave-induced vertical fluctuations than does the viscous curvilinear model.

In summary, the results in § 4.1 demonstrate that incorporating turbulence stresses is essential for accurately modelling wave-induced vertical velocity perturbations at small wave ages. By accounting for these stresses, the turbulence curvilinear model provides a more comprehensive representation of the phase relationship between the wave motion and airflow across a wide range of wave ages, which underscores the importance of considering turbulence stress effects at small magnitudes of wave age where turbulent mixing significantly influences the flow dynamics.

4.2. Comparing models of $\textrm {Re}{\{\hat {u}\}}$ and $\textrm {Im}{\{\hat {u}\}}$ with LES results

In this section, we further illustrate the effectiveness of the turbulence curvilinear model in predicting the wave-induced streamwise perturbation $(\widetilde {u})$ . We compare the viscous and turbulence curvilinear model predictions of the real component of the wave-induced streamwise velocity (Re{ $\hat {u}$ }) with the LES data in figure 8. Given the solution for $\textrm {Re}\{\hat {w}\}$ from the viscous or turbulence curvilinear models, $\widetilde {u}$ is determined from the continuity (2.14c ), and its real component is given by (2.23).

Equation (2.23) indicates that $\textrm {Re}\{\hat {u}\}$ , which denotes the in-phase component of $\hat {u}$ with respect to the phase of the wave surface elevation, is governed by the vertical gradient of the out-of-phase component of $\hat {w}$ . As a result, the accuracy of the $\textrm {Re}\{\hat {u}\}$ prediction depends on the accuracy of the $\textrm {Im}\{\hat {w}\}$ prediction. In § 4.1, it is shown that the viscous and turbulence curvilinear models yield similar predictions for $\textrm {Im}\{\hat {w}\}$ across all wave ages, except that the viscous model slightly under-predicts the values in the $c/u_\tau =2,7$ and $15$ cases. This underestimation leads to deviations in the predictions of $\textrm {Re}\{\hat {u}\}$ made by the viscous curvilinear model (figure 8 f,e,h). In contrast, the turbulence curvilinear model predictions align well with the LES results for these wave ages.

Figure 8. Vertical profiles of $\textrm {Re}\{\hat {u}\}/u_\tau$ obtained from the viscous () and turbulence () curvilinear model compared with the LES data () for various wave ages: (a,i) $c/u_\tau =\pm 25$ , (b,h) $c/u_\tau =\pm 15$ , (c,g) $c/u_\tau =\pm 7$ , (d, f) $c/u_\tau =\pm 2$ , (e) $c/u_\tau =0$ .

Next, we discuss the characteristics and model predictions of Im{ $\hat {u}$ }, which represents the out-of-phase component of the wave-induced streamwise velocity perturbation relative to the wave elevation. The efficacy of the viscous and turbulence curvilinear models in predicting Im{ $\hat {u}$ } is demonstrated through the comparison of the model predictions with LES results in figure 9. The quantity $\textrm {Im}\{\hat {u}\}$ is obtained from (2.24), which represents the imaginary component of the continuity (2.14c ) and depends on the vertical gradient of $\textrm {Re}\{\hat {w}\}$ . Therefore, similar to the performance of the predictions of $\textrm {Re}\{\hat {w}\}$ for the $c/u_\tau =-25,-15$ and $-7$ cases, figure 9(a–c) shows that both the viscous and turbulence curvilinear models capture the general trend of Im{ $\hat {u}$ }. However, the turbulence curvilinear model reproduces the LES data more accurately, particularly in the near-surface region. This improvement suggests that turbulence stresses, which are neglected in the viscous model, play an important role in determining the phase relationship between the wave motion and the wave-induced streamwise velocity for these rapid, opposing waves. The vertical gradient of $\textrm {Re}\{\hat {w}\}$ approaches zero at the wave surface as shown in figure 6. As a result, $\textrm {Im}\{u\}(\zeta =0)=\textrm {Im}\{\hat {u}_s\}=0$ , which aligns with the boundary condition $\hat {u}_s=akc{/2} { \ +\ \rm i \ 0}$ from (2.21).

Figure 9. Vertical profiles of $\textrm {Im}\{\hat {u}\}/u_\tau$ obtained from the viscous () and turbulence () curvilinear model compared with the LES data () for various wave ages: (a,i) $c/u_\tau =\pm 25$ , (b,h) $c/u_\tau =\pm 15$ , (c,g) $c/u_\tau =\pm 7$ , (d, f) $c/u_\tau =\pm 2$ , (e) $c/u_\tau =0$ .

When the wave age is nearly zero as in the $c/u_\tau =-2$ , 0 (stationary wavy wall) and 2 cases, figure 9(d– f) demonstrates a growing discrepancy between the viscous model’s predictions and the LES data. The turbulence curvilinear model, in contrast, maintains better agreement, which highlights the increasing importance of turbulent mixing in accurately representing the phase lag between the wave and the streamwise airflow. The turbulence curvilinear model’s ability to capture the peak magnitude and the vertical decay of Im{ $\hat {u}$ } indicates that turbulent diffusion becomes increasingly dominant as the wave speed magnitude decreases.

For the wind-following wave cases with $c/u_\tau =7$ and 15, figure 9(g,h) shows that the turbulence curvilinear model provides a better representation of the LES data than does the viscous curvilinear model. The viscous model prediction for Im $\{\hat {u}\}$ deviates from the LES data, which leads to differences in the observed peak and the vertical decay characteristics. In figure 9(i), for the $c/u_\tau =25$ case, the viscous curvilinear model provides a good prediction of the peak magnitude and its location near the wave surface, whereas the turbulence curvilinear model captures the decay trend at larger heights more accurately.

To summarise § 4.2, we again confirm that the turbulence curvilinear model provides improved predictions of both $\textrm {Re}\{\hat {u}\}$ and $\textrm {Im}\{\hat {u}\}$ relative to those of the viscous curvilinear model across the various wave ages.

4.3. Comparison of models of $\textrm {Re}{\{\hat {p}\}}$ and $\textrm {Im}{\{\hat {p}\}}$ with LES results

In this section, we examine the properties of Re{ $\hat {p}$ } and Im{ $\hat {p}$ }, which are the Fourier modes of the wave-induced pressure perturbation, and compare the viscous and turbulence curvilinear models with the LES data. The quantity Re{ $\hat {p}$ } represents the in-phase component of the pressure fluctuations induced by the wave and reflects the direct effect of wave-induced pressure gradients on the flow field. The expression for $\textrm {Re}\{\hat {p}\}$ is obtained by integrating the vertical momentum (2.14b ) and is given by (2.28). We rewrite (2.28) as

(4.1) \begin{align} \textrm {Re}\{\hat {p}\}&=\textrm {Re}\{\hat {p}_{\textit{ad}v}\}+\textrm {Re}\{\hat {p}_{\mathrm{\nu }}\}+\textrm {Re}\{\hat {p}_{\textit{turb}}\}, \end{align}

where $\textrm {Re}\{\hat {p}_{\textit{ad}v}\}, \textrm {Re}\{\hat {p}_{\mathrm{\nu }}\}$ and $\textrm {Re}\{\hat {p}_{\textit{turb}}\}$ are the real parts of the wave-induced pressure due to advection, viscous and turbulence stress effects, respectively. These components have the following expressions:

(4.2a) \begin{align} \textrm {Re}\{\hat {p}_{\textit{ad}v}\}&=-\int _{\zeta }^{\lambda } (\langle u \rangle -c)k\,\textrm {Im}\{\hat {w}\} \textrm {d}\zeta , \\[-12pt] \nonumber \end{align}

(4.2b) \begin{align} \textrm {Re}\{\hat {p}_{\mathrm{\nu }}\}&=- \int _{\zeta }^{\lambda } \nu\! \left (\frac {\partial ^2}{\partial \zeta ^2}-k^2\right )\textrm {Re}\{\hat {w}\} \ \textrm {d}\zeta , \\[-12pt] \nonumber \end{align}

(4.2c) \begin{align} \textrm {Re}\{\hat {p}_{\textit{turb}}\}&=-\int _{\zeta }^{\lambda } k\, \textrm {Im}\{\hat {\tau }_{13}\} \ \textrm {d}\zeta -\textrm {Re}\{\hat {\tau }_{33}\}(\zeta ). \\[9pt] \nonumber \end{align}

The viscous and turbulence curvilinear model predictions of $\textrm {Re}\{\hat {p}\}$ are compared with the LES data in figure 10(a–i) across different wave ages. The contributions of the individual components from (4.2a – c ) to the turbulence curvilinear model are also illustrated in figure 10. It is clear from figure 10(a–i) that both the viscous and turbulence curvilinear models capture the overall vertical variation of $\textrm {Re}\{\hat {p}\}$ , closely matching the LES results, except for the $c/u_\tau =7$ and $15$ cases. It is evident that, for all cases, the advection component $\textrm {Re}\{\hat {p}_{\textit{ad}v}\}$ contributes the most to $\textrm {Re}\{\hat {p}\}$ . Since $\textrm {Re}\{\hat {p}_{\textit{ad}v}\}$ is determined from $\textrm {Im}\{\hat {w}\}$ (4.2a ), which exhibits minimal difference between the viscous and turbulence curvilinear models (see figure 7), the resulting predictions for $\textrm {Re}\{\hat {p}\}$ are similar between the two models, with a slight under-prediction relative to the LES results. For the $c/u_\tau =7$ and $15$ cases, the turbulence curvilinear model provides a better prediction of $\textrm {Re}\{\hat {p}\}$ than does the viscous curvilinear model. For these two cases, figure 10(g,h) shows that the finite contribution from $\textrm {Re}\{\hat {p}\}_{\textit{turb}}$ enhances the turbulence curvilinear model prediction. Because the viscous model neglects the turbulence stress effects, its predictions do not agree with LES at $c/u_\tau =7$ and $15$ .

Figure 10. Vertical profiles of $\textrm {Re}\{\hat {p}\}/u_\tau ^2$ obtained from the viscous () and turbulence () curvilinear model compared with the LES data () for various wave ages: (a,i) $c/u_\tau =\pm 25$ , (b,h) $c/u_\tau =\pm 15$ , (c,g) $c/u_\tau =\pm 7$ , (d, f) $c/u_\tau =\pm 2$ , (e) $c/u_\tau =0$ . The contributions of $\textrm {Re}\{\hat {p}_{\textit{ad}v}\}$ (), $\textrm {Re}\{\hat {p}_{\nu }\}$ (), $\textrm {Re}\{\hat {p}_{\textit{turb}}\}$ () from the turbulence curvilinear model are also plotted in each panel.

Notably, although the turbulence curvilinear model correctly captures the trend in the vertical variation of $\textrm {Re}\{\hat {p}\}$ for the $c/u_\tau =7$ case, it does not predict the magnitude accurately compared with LES. In contrast, the turbulence model prediction for the $c/u_\tau =15$ case agrees well with LES. The larger error for $c/u_\tau =7$ could be attributed to the neglect of wave-induced stresses ( $\widetilde {\tau }_{\textit{ij}}^w$ ) in the turbulence curvilinear model (2.15). In § 6, from the discussion of the wave-coherent energy budget, we find that the magnitude of wave-coherent energy production due to the correlation between the wave-coherent pressure stress and the velocity gradients is greater for the $c/u_\tau =7$ case than for other wave-age cases. This difference results in enhanced wave-induced stress for the $c/u_\tau =7$ case, which leads to the under-prediction of the $\textrm {Re}\{\hat {p}\}$ magnitude when the wave-induced stress is neglected.

Next, we discuss the properties of Im{ $\hat {p}$ }, which represents the out-of-phase component (relative to the cosine profile of the wave surface) of the pressure fluctuations induced by the wave and reflects the phase delay between the wave motion and the resulting pressure response. Accurately predicting this component is particularly important, because it is the out-of-phase part of the wave-induced pressure perturbation that contributes to the form drag, as shown in (2.32). The expression for $\textrm {Im}\{\hat {p}\}$ is obtained by integrating the vertical momentum (2.26) and is given by (2.29), which yields

(4.3) \begin{align} \textrm {Im}\{\hat {p}\}&=\textrm {Im}\{\hat {p}_{\textit{ad}v}\}+\textrm {Im}\{\hat {p}_{\nu }\}+\textrm {Im}\{\hat {p}_{\textit{turb}}\}, \end{align}

where $\textrm {Im}\{\hat {p}_{\textit{ad}v}\}$ , $\textrm {Im}\{\hat {p}_{\nu }\}$ and $\textrm {Im}\{\hat {p}_{\textit{turb}}\}$ are the advection, viscous and turbulence stress components that contribute to $\textrm {Im}\{\hat {p}\}$ . From (2.29), these components are expressed as

(4.4a) \begin{align} \textrm {Im}\{\hat {p}_{\textit{ad}v}\}&=\int _{\zeta }^{\lambda } (\langle u \rangle -c)k\,\textrm {Re}\{\hat {w}\} \ \textrm {d}\zeta , \\[-12pt] \nonumber \end{align}

(4.4b) \begin{align} \textrm {Im}\{\hat {p}_\nu \}&=- \int _{\zeta }^{\lambda } \nu\! \left (\frac {\partial ^2}{\partial \zeta ^2}-k^2\right )\textrm {Im}\{\hat {w}\} \ \textrm {d}\zeta , \\[-12pt] \nonumber \end{align}

(4.4c) \begin{align} \textrm {Im}\{\hat {p}_{\textit{turb}}\}&=\int _{\zeta }^{\lambda } k\, \textrm {Re}\{\hat {\tau }_{13}\} \ \textrm {d}\zeta -\textrm {Im}\{\hat {\tau }_{33}\}(\zeta ). \\[9pt] \nonumber \end{align}

The profiles of Im{ $\hat {p}$ } from the LES data for various wave ages are compared with the respective viscous and turbulence curvilinear model predictions in figure 11(a–i). In addition to showing $\textrm {Im}\{\hat {p}\}$ , figure 11 also presents its individual components (4.4a – c ) predicted by the turbulence curvilinear model. Only the resulting $\textrm {Im}\{\hat {p}\}$ profile is shown for the viscous curvilinear model.

From figure 11, it is evident that the predictions of the turbulence curvilinear model $\textrm {Im}\{\hat {p}\}$ are in better agreement with the LES results than those of the viscous model are. It is also evident that the contribution of $\textrm {Im}\{\hat {p}_{\textit{turb}}\}$ to $\textrm {Im}\{\hat {p}\}$ is noticeable, which, when neglected as in the viscous curvilinear model, results in significant deviations in the results. The error associated with the neglect of the turbulence stress component in the viscous curvilinear model predictions of $\textrm {Re}\{\hat {w}\}$ is also reflected in the corresponding $\textrm {Im}\{\hat {p}_{\textit{ad}v}\}$ predictions, resulting in additional errors in the predictions for $\textrm {Im}\{\hat {p}\}$ .

Figure 11. Vertical profiles of $\textrm {Im}\{\hat {p}\}/u_\tau ^2$ obtained from the viscous () and turbulence () curvilinear model compared with the LES data () for various wave ages: (a,i) $c/u_\tau =\pm 25$ , (b,h) $c/u_\tau =\pm 15$ , (c,g) $c/u_\tau =\pm 7$ , (d, f) $c/u_\tau =\pm 2$ , (e) $c/u_\tau =0$ . The contributions of $\textrm {Im}\{\hat {p}_{\textit{ad}v}\}$ (), $\textrm {Im}\{\hat {p}_{\nu }\}$ (), $\textrm {Im}\{\hat {p}_{\textit{turb}}\}$ () from the turbulence curvilinear model are also plotted in each panel.

To summarise § 4.3, the comparison in this section demonstrates that incorporating turbulence stresses is essential for accurately modelling wave-induced pressure perturbation across a wide range of wave ages. By accounting for these stresses, the turbulence curvilinear model provides more physically faithful representations of both the magnitude and phase relationship between the wave motion and the pressure field. The limitations of the viscous curvilinear model in capturing the observed trends underscore the importance of considering turbulence effects, particularly in scenarios where turbulent mixing significantly influences the magnitude and phase variations of the pressure fluctuations and the resulting form drag, which is discussed next.

4.4. Comparison of models of form drag with LES data

This section compares the accuracy of form drag predictions from the viscous and turbulence curvilinear models. The form drag is evaluated using (2.32), where the $\textrm {Im}\{\hat {p}(\zeta =0)\}$ component of the wave-induced pressure perturbation contributes to the form drag. The accuracy of the model predictions is assessed by comparing the model results with the LES data, as shown in figure 12. Under higher wave age and wind-following wave conditions ( $c/u_\tau =7$ to $25$ ), as well as for the wind-opposing wave case with wave age $c/u_\tau =-25$ , the two models produce similar results, which indicates that the effects of turbulence stresses are less significant in these regimes. However, for the cases with the wave age $c/u_\tau$ ranging from $-15$ to 2, the turbulence curvilinear model provides better agreement with LES than does the viscous curvilinear model. Overall, the mean absolute percentage error (MAPE) of the turbulence curvilinear model relative to LES is only 4 %, whereas the corresponding error for the viscous curvilinear model is 29 %.

Because $\textrm {Im}\{\hat {p}_{\textit{ad}v}\}$ is the dominant contributor to $\textrm {Im}\{\hat {p}\}$ , as discussed in § 4.3, the resulting form drag induced by the linear advection term, $F_p^{\textit{ad}v}$ in (2.34), predominates the total drag, $F_p$ given by (2.33). This is evident from the comparison of the $F_p^{\textit{ad}v}$ predictions of the turbulence and viscous curvilinear models with the LES results shown in figure 12. In the viscous curvilinear model, the behaviour of $F_p^{\textit{ad}v}$ closely resembles that of $F_p$ . This result is consistent with the findings of Åkervik & Vartdal (Reference Åkervik and Vartdal2019), Cao et al. (Reference Cao, Deng and Shen2020) and Cao & Shen (Reference Cao and Shen2021). The viscous curvilinear model under-predicts the form drag for cases with $c/u_\tau$ ranging from $-7$ to $2$ . The MAPE obtained for $F_p^{\textit{ad}v}$ in the viscous curvilinear model is 29 %. In contrast, the MAPE obtained for $F_p^{\textit{ad}v}$ in the turbulence curvilinear model decreases to 16 %, which indicates that $F_p^{\textit{ad}v}$ of the turbulence curvilinear model aligns more closely with the LES results. Moreover, when both $F_p^{\textit{ad}v}$ and $F_p^{\textit{turb}}$ are considered, the error decreases to 4 %, which underscores the importance of both the linear advection and turbulence stress contributions for accurately capturing the total drag exerted by waves on the airflow.

Figure 12. Comparison of the form drag $F_p$ (2.32) predictions from the viscous and turbulence curvilinear models with LES results. The contributions of the advection-induced form drag $(F_p^{\textit{ad}v})$ are also plotted.

In summarise § 4, a priori tests demonstrate that incorporating turbulence stress terms improves the prediction of wave-induced airflow perturbations and form drag. The turbulence curvilinear model performs better than the viscous curvilinear model does. In the following section, we discuss solving the turbulence curvilinear model using eddy viscosity-based parameterisations for the wave-induced turbulence stresses.

5.1. Properties of eddy viscosity estimated from mean turbulence stress and velocity gradient

In this section, the characteristics of the eddy viscosity $\nu _T(\zeta )$ that corresponds to the background turbulence across the various wave ages are discussed. This $\nu _T(\zeta )$ for the background turbulence is computed from the mean shear stress $\langle \tau _{13}\rangle$ and the mean velocity gradient $\partial \langle u\rangle /\partial \zeta$ using $\nu _T(\zeta ) = \langle \tau _{13}\rangle / (\partial \langle u\rangle / \partial \zeta )$ from LES data for each wind over wave case. The corresponding vertical profiles of the eddy viscosity are shown in figure 13(a–b). The coloured markers indicate $\nu _T(\zeta )$ for different wave ages. The $\nu _T(\zeta )$ is normalised by the molecular viscosity $\nu$ to delineate the regions influenced by the mean turbulent flow. The vertical axis in figure 13(a) is normalised by the viscous length scale ( $\nu /u_\tau$ ), which emphasises the inner region, whereas in figure 13(b), it is normalised by the wavenumber $(k)$ to highlight the outer region.

Figure 13(a) shows that $\nu _T(\zeta )/\nu \lt 1$ within the viscous sublayer region ( $\zeta ^+ \lt 10$ ), which indicates that viscous stresses dominate the turbulence stresses in this zone. Furthermore, the $\nu _T(\zeta )$ profiles exhibit self-similarity within the viscous sublayer, which suggests that wave age has a minimal influence in this region. The mean-flow profiles in figure 2(a) also demonstrate a self-similar behaviour in the viscous sublayer across all the wave-age cases, which follows the linear relationship $\langle u\rangle ^+ = \zeta ^+$ .

Figure 13. Panels (a) and (b) present the vertical profiles of the eddy viscosity that corresponds to the mean turbulent flow in the background. The vertical axis is shown on a logarithmic scale in panel (a) and on a linear scale in panel (b). The coloured circular markers indicate the eddy viscosity profiles calculated from mean Reynolds stress and velocity gradient using $\nu _T(\zeta ) = \langle \tau _{13} \rangle / ( \partial \langle u \rangle / \partial \zeta )$ for the wave-age cases $c/u_\tau = -25$ (), $-15$ (), $-7$ (), $-2$ (), $0$ (), $2$ (), $7$ (), $15$ () and $25$ (). The eddy viscosity profile for the flat-wall case is shown with black circular markers. The dashed black line represents the model based on $\nu _{T,{I}}$ (5.1), whereas the coloured solid lines, following the same colour scheme as that used for the circular markers, correspond to $\nu _{T,\textit{II}}$ (5.2). The dashed cyan line corresponds to $\nu _{T,\textit{III}}$ (5.3), the Cess-type eddy viscosity profile. Panel (c) shows the variation of the parameter $C_1$ in (5.2) for $\nu _{T,\textit{II}}$ capturing the different magnitudes of the eddy viscosity across the various wave-age cases. The dotted black line in panel (c) is a polynomial fit to $C_1$ given by $C_1(c/u_\tau )=\alpha _1 {(c/u_\tau )}^3+\alpha _2{(c/u_\tau )}^2 + \alpha _3 {(c/u_\tau )} +\alpha _4$ , where $|c/u_\tau |\leqslant 25$ and the coefficients are $\alpha _1=-1.3\times 10^{-5}, \alpha _2=3.95\times 10^{-4}, \alpha _3=-1.11\times 10^{-2}, \alpha _4=0.964$ .

Although the $\nu _T/\nu$ profiles may initially appear self-similar in the outer region ( $\zeta ^+ \gg 10$ ) in figure 13(a), linear vertical scaling reveals the influence of wave age on the eddy viscosity profiles, as shown in figure 13(b). While the $\nu _T/\nu$ profile in the $c/u_\tau =25$ case closely resembles that of a flat wall, the peak magnitude of eddy viscosity increases progressively when $c/u_\tau$ decreases from 25 to $-25$ . This trend indicates that all the wind over wave cases with $c/u_\tau$ ranging from $-25$ to $15$ induce stronger turbulent fluctuations than the $c/u_\tau = 25$ and flat-wall cases do.

The black dashed line in figure 13(a,b) represents the eddy viscosity model with van Driest damping for capturing the decay of Reynolds stresses in the viscous region (van Driest Reference van Driest1956; Piomelli & Balaras Reference Piomelli and Balaras2002), which is given by

(5.1) \begin{align} \nu _{T,{I}}=u_\tau \kappa \zeta \left [1-\exp \left (-{\frac {\zeta }{H}\frac {\textit{Re}_\tau }{25}}\right )\right ]\!. \end{align}

Similarly, the coloured solid lines in figure 13(a,b) represent a modified version of the eddy viscosity $\nu _{T,{I}}$ given in (5.1). This modification includes a power-law-based damping function (Lars, Hermann & Klaus Reference Lars, Hermann and Klaus2003) to account for the decay of Reynolds stress at greater heights, along with a constant coefficient $C_1$ that varies with wave age to model the peak magnitude differences across the various wind over wave scenarios. The resulting expression for this wave-age-dependent eddy viscosity is as follows:

(5.2) \begin{align} \nu _{T,\textit{II}}=C_1 u_\tau\! \left (1-\frac {\zeta }{H}\right )^{0.8} \kappa \zeta\! \left [1-\exp \left (-{\frac {\zeta }{H}\frac {\textit{Re}_\tau }{25}}\right )\right ]\!. \end{align}

The above models $\nu _{T,{I}}$ and $\nu _{T,\textit{II}}$ , describe eddy viscosity as the product of a characteristic velocity scale ( $u$ ) and a characteristic length scale ( $\ell$ ), following Prandtl’s mixing-length theory, where $\nu _T \sim u \ell$ .

Equation (5.1) models eddy viscosity as $\nu _T \sim u_\tau \kappa \zeta$ , with the friction velocity $u_\tau$ serving as the characteristic velocity scale and $\ell \sim k\zeta$ as the characteristic length scale. Similarly, (5.2) expresses eddy viscosity as $\nu _T \sim C_1 u_\tau (1 - \zeta /H)^{0.8} \kappa \zeta$ , where the characteristic velocity scale is given by $u \sim C_1 u_\tau (1 - \zeta /H)^{0.8}$ and the characteristic length scale remains $\ell \sim \kappa \zeta$ . Both model expressions incorporate the van Driest damping function $(1 - \exp [-{(\zeta /H) \textit{Re}_\tau / 25}])$ to account for the reduction in eddy viscosity within the viscous sublayer.

It is clear from figure 13(a,b) that the effect of wave age and the decrease in the eddy viscosity at larger heights are not captured by (5.1). On the other hand, (5.2) models the variation in eddy viscosity across different wave ages through the model parameter $C_1(c/u_\tau )$ , and the factor $(1 - \zeta /H)^{0.8}$ models the decay of characteristic velocity at larger heights (Lars et al. Reference Lars, Hermann and Klaus2003). Both parameters, $C_1(c/u_\tau )$ (figure 13 c) and the exponent 0.8 in (5.2) are determined by fitting the model expression to the eddy viscosity profiles from the LES data. Furthermore, the dotted line in figure 13(c) represents a polynomial fit for $C_1(c/u_\tau )$ , which is given by $C_1(c/u_\tau )=\alpha _1 {(c/u_\tau )}^3+\alpha _2 {(c/u_\tau )}^2 + \alpha _3 {(c/u_\tau )} +\alpha _4$ with fitting coefficients $\alpha _1 = -1.3 \times 10^{-5}$ , $\alpha _2 = 3.95 \times 10^{-4}$ , $\alpha _3 = -1.11 \times 10^{-2}$ and $\alpha _4 = 0.964$ . This polynomial fit is valid for wave ages within the range of $|c/u_\tau | \leqslant 25$ .

In addition to the $\nu _{T,{I}}$ and $\nu _{T,\textit{II}}$ profiles, figure 13(a,b) also presents the Cess-type eddy viscosity profile (Reynolds & Hussain Reference Reynolds and Hussain1972; Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009) given by

(5.3) \begin{align} \nu _{T,\textit{III}}&=\frac {\nu }{2}\!\left [1+\frac {\kappa ^2 \textit{Re}_\tau ^2}{9} (1-\zeta _{\textit{Cess}}^2 )^2 (1+2\zeta _{\textit{Cess}}^2 )^2 (1-\exp [ (|\zeta _{\textit{Cess}}|-1 )\textit{Re}_\tau /25 ] )^2\right ]^{1/2} \nonumber \\ & \quad -\frac {\nu }{2}, \end{align}

where $\zeta _{\textit{Cess}}=\zeta -1$ . This analytical expression, originally developed for fully developed turbulent pipe flow (Reynolds & Tiederman Reference Reynolds and Tiederman1967), provides a smooth transition of eddy viscosity from the viscous sublayer through the buffer and logarithmic layers. As shown in figure 13(b), the eddy viscosity profile for the wind-following wave case at $c^+ = 25$ closely follows the $\nu _{T,\textit{III}}$ profile, consistent with the corresponding mean velocity profile in figure 2(b), which approaches the flat-wall behaviour.

We next evaluate the suitability of these eddy viscosity profiles for representing wave-induced turbulence stresses and predicting the associated wave-coherent airflow perturbations.

5.2. Suitability of mean-flow-based eddy viscosity model

Figure 14. Comparison of $\mathrm{Re}\{\hat {w}\}$ (a,b,c,d) and $\mathrm{Im}\{\hat {w}\}$ (e, f,g,h) predicted from the viscous curvilinear model () and the turbulence curvilinear models with $\nu _{T,{I}}$ (), $\nu _{T,\textit{II}}$ (), $\nu _{T,\textit{III}}$ () profiles and the a priori test of the turbulence curvilinear model () with the LES () results.

In this section, we evaluate the suitability of the eddy viscosity models given by (5.1), (5.2) and (5.3) for solving the turbulence curvilinear model (2.19). We compare the turbulence curvilinear model predictions of the vertical profiles of $\textrm {Re}\{\hat {w}\}$ and $\textrm {Im}\{\hat {w}\}$ using (5.1), (5.2) and (5.3) across various wave ages in figure 14. The figure presents the predictions obtained using the eddy viscosity models $\nu _{T,{I}}$ , $\nu _{T,\textit{II}}$ and $\nu _{T,\textit{III}}$ alongside the LES data. Additionally, results from the viscous curvilinear model and the a priori test of the turbulence curvilinear model in § 4 are also included for comparison purposes.

Although the use of the eddy viscosity models (5.1), (5.2) and (5.3) to solve the turbulence curvilinear model equation provides good predictions for $\textrm {Im}\{\hat {w}\}$ (figure 14 e–h), the predictions for $\textrm {Re}\{\hat {w}\}$ (figure 14 a–d) exhibit significant deviations from the LES data. Furthermore, the predictions of $\textrm {Re}\{\hat {w}\}$ from the eddy viscosity-based turbulence curvilinear model underperform compared with both the viscous model and the turbulence curvilinear model a priori test results. In addition, in figure 15, the predictions of the form drag from the turbulence curvilinear model obtained by using the eddy viscosity models deviate significantly from the LES data. This discrepancy indicates that the mean-flow-based eddy viscosity models do not satisfactorily capture the wave-induced stress characteristics.

Figure 15. Predictions of form drag ( $F_p$ ) from the turbulence curvilinear model using the $\nu _{T,{I}}$ (5.1), $\nu _{T,\textit{II}}$ (5.2) and $\nu _{T,\textit{III}}$ (5.3) eddy viscosity parameterisations along with the a priori test result (see § 4), and their comparison with the LES data of $F_p$ .

Consequently, the eddy viscosity derived from the mean turbulence shear stress component ( $\langle \tau _{13}\rangle$ ) and the mean velocity gradient ( $\partial \langle u\rangle /\partial \zeta$ ) fails to accurately represent the wave-induced turbulence stresses ( $\widetilde {\tau }_{\textit{ij}}$ ). This incorrect prediction of $\widetilde {\tau }_{\textit{ij}}$ using the mean-flow-based eddy viscosity is the primary reason for the discrepancy in predicting $\hat {w}$ as shown above when the turbulence curvilinear model is employed. In the following section, we explore the appropriate eddy viscosity profile required for capturing the wave-induced turbulence stresses.

5.3. Eddy viscosity characteristics from wave-coherent turbulent motions

Since the mean-flow-based eddy viscosity profiles (§ 5.1) do not adequately capture wave-induced Reynolds stresses (§ 5.2), we now examine the eddy viscosity profiles associated with the latter and contrast them with their mean-flow-based counterparts. In § 3.4.2, we showed that wave-induced turbulence stresses are anisotropic. As a result, we relax the isotropic eddy viscosity assumption in (2.18) and instead construct a tensorial form of the eddy viscosity as follows:

(5.4) \begin{align} {\hat {\nu }_{T,\textit{ij}}}=\begin{bmatrix} \dfrac {\hat {\tau }_{11}}{2\dfrac {\partial \hat {w}}{\partial \zeta }}& \dfrac {-\hat {\tau }_{13}}{\dfrac {\partial \hat {u}}{\partial \zeta }+\mathrm{i}k\hat {w}+\dfrac {\partial \langle u \rangle } {\partial \zeta }g_\zeta \hat {\eta }}\\ &\\ \dfrac {\hat {\tau }_{31}}{\dfrac {\partial \hat {u}}{\partial \zeta }+\mathrm{i}k\hat {w}} & \dfrac {\hat {\tau }_{33}}{\mathrm{i}k\hat {u}-\dfrac {\partial \hat {w}}{\partial \zeta }} \end{bmatrix}\!, \text{for $i,j=1,3$} . \end{align}

Note that (5.4) does not have the $\widetilde {\tau }_{22}$ stress component, as it does not contribute directly to the turbulence forcing terms in the turbulence curvilinear model (2.15).

Because the Fourier modes $\hat {\tau }_{\textit{ij}}$ and $\hat {u}_i$ in (5.4) are complex-valued numbers, the resulting eddy viscosity profile, ${\hat {\nu }_{T,\textit{ij}}}=|{\hat {\nu }_{T,\textit{ij}}}(\zeta )|\mathrm{e}^{\text{i}\phi _{\textit{ij}}(\zeta )}$ , is also complex, where $|{\hat {\nu }_{T,\textit{ij}}}|$ is the magnitude of the eddy viscosity and $\phi _{\textit{ij}}(\zeta )=\text{arctan}(\textrm {Im}\{{\hat {\nu }_{T,\textit{ij}}}\}/\textrm {Re}\{{\hat {\nu }_{T,\textit{ij}}}\})$ is the phase angle. The $\textrm {Re}\{{\hat {\nu }_{T,\textit{ij}}}\}$ and $\textrm {Im}\{{\hat {\nu }_{T,\textit{ij}}}\}$ components both contribute to determining $\textrm {Re}\{\hat {w}\}$ and $\textrm {Im}\{\hat {w}\}$ when substituted into (2.19).

Figure 16. Vertical profiles of the magnitude (a–h) and phase (i–l) of the complex-valued eddy viscosity profile ${\hat {\nu }_{T,\textit{ij}}}=|{\hat {\nu }_{T,\textit{ij}}}(\zeta )|\mathrm{e}^{\text{i}\phi _{\textit{ij}}(\zeta )}$ constructed from (2.18) for the cases with $c/u_\tau =-25$ (), $-15$ (), $-7$ (), $-2$ (), $0$ (), $2$ (), $7$ (), $15$ (), $25$ (). In panels (a–d) and (i–l), the vertical axis is presented on a logarithmic scale with inner-unit normalisation, whereas in panels (e–h), it is shown on a linear scale and normalised by the wavenumber.

Next, to understand the characteristics of these complex-valued eddy viscosities from (5.4), we construct the vertical profiles of the magnitude and phase of $\hat {\nu }_{T,ij}$ using LES data. The results are shown in figure 16(a–l). In figure 16(a–h), the magnitude of the eddy viscosity ( $|{\hat {\nu }_{T,\textit{ij}}}|$ ) is normalised by the molecular viscosity to illustrate the enhancement of mixing by wave-induced turbulence stresses relative to molecular diffusion. In figure 16(a–d), the vertical coordinate is normalised by the viscous length scale ( $\nu /u_\tau$ ) and plotted on a logarithmic scale to highlight the inner region. On the other hand, in figure 16(e–h), the vertical coordinate is normalised by the wavenumber ( $k$ ) and a linear scale is used to emphasise the outer region behaviour.

From figure 16(a–d), we observe a viscous sublayer in the region $\zeta \lt 10(\nu /u_\tau )$ , where $|{\hat {\nu }_{T,\textit{ij}}}| \lt \nu$ . For $\zeta \gg 10(\nu /u_\tau )$ or $k\zeta \gt 1$ , $|{\hat {\nu }_{T,\textit{ij}}}|$ increases above $\nu$ and approaches a constant value at larger heights, so that $|{\hat {\nu }_{T,\textit{ij}}}|/\nu$ spans $\mathcal{O}(10)$ – $\mathcal{O}(10^2)$ across various wave-age cases, as shown in figure 16(e–h). This increased eddy viscosity magnitude relative to $\nu$ reflects the strong influence of wave-induced turbulence. Additionally, the wave-induced eddy viscosity exhibits a complicated vertical variation in its phase angle (figure 16 i–l). This non-trivial height dependence of $\phi _{\textit{ij}}$ underscores its role in modulating momentum transfer by influencing both the in-phase and out-of-phase components of the wave-induced fluctuations. Such behaviour of $\hat {\nu }_{T,ij}$ further emphasises the limitations of simple mean-flow-based eddy viscosity closures discussed in § 5.1, which cannot adequately capture the phase-dependent structure of turbulence in wave-induced flows.

Figure 17. Variations in the critical-layer thickness $\zeta _{\textit{critical}}$ and inner-layer thickness $\zeta _{\textit{inner}}$ across various wave ages.

Assuming the $|{\hat {\nu }_{T,\textit{ij}}}|\sim u_\tau \ell$ scaling, we can also interpret figure 16(a–h) in terms of the variation of the mixing length, $\ell ^+ =\ell /(\nu /u_\tau )(\equiv |{\hat {\nu }_{T,\textit{ij}}}|/\nu )$ , which corresponds to the sizes of the wave-induced turbulent eddies. Therefore, the vertical profiles of $|{\hat {\nu }_{T,\textit{ij}}}|/\nu$ in figure 16(a–h) also signify that, when moving away from the wave surface, the eddy length scale increases first and then asymptotically approaches a constant value at greater heights. We find that this limiting length scale is determined by the inner-layer thickness ( $\zeta _{\textit{inner}}$ ) of the wave boundary layer, such that $|{\hat {\nu }_{T,\textit{ij}}}|\sim u_\tau \zeta _{\textit{inner}}$ . The quantity $\zeta _{\textit{inner}}$ corresponds to the height where the advection time scale $T_{\textit{ad}v}\sim \lambda /|\langle u\rangle (\zeta )-c|$ equals the eddy-turnover time scale $T_{\textit{eddy}}\sim \kappa \zeta /u_\tau$ (Belcher & Hunt Reference Belcher and Hunt1993). For $\zeta \lt \zeta _{\textit{inner}}$ , the eddy-turnover time scale is smaller than the advection time scale ( $T_{\textit{eddy}}\lt T_{\textit{ad}v}$ ), whereas for $\zeta \gt \zeta _{\textit{inner}}$ , the eddy-turnover time scale is larger ( $T_{\textit{eddy}}\gt T_{\textit{ad}v}$ ). The physical meaning of this result is that when moving away from the wave surface, the sizes of wave-induced turbulent eddies increase until $\zeta \sim \zeta _{\textit{inner}}$ and they are subsequently advected rapidly for $\zeta \gt \zeta _{\textit{inner}}$ . Following Belcher & Hunt (Reference Belcher and Hunt1993) and Mastenbroek et al. (Reference Mastenbroek, Makin, Garat and Giovanangeli1996), we determine the inner-layer thickness from

(5.5) \begin{align} k\zeta _{\textit{inner}} & = \frac {2\kappa u_\tau }{|\langle u \rangle (\zeta _{\textit{inner}})-c|}. \end{align}

Because the mean wind speed in (5.5) depends on the wave age (see figure 2) and $c$ also explicitly appears in the denominator, the resulting inner-layer thickness varies across the different wind over wave cases, as illustrated in figure 17. The $\zeta _{\textit{inner}}$ shown in figure 17 is obtained through a graphical procedure. We define

(5.6) \begin{align} f\!(\zeta ) = \biggl |k\zeta - \frac {2\kappa u_\tau }{|\langle u\rangle (\zeta ) - c|}\biggr |. \end{align}

The inner-layer height $\zeta _{\textit{inner}}$ corresponds to the value of $\zeta$ for which $f(\zeta _{\textit{inner}})=0$ , i.e. where the line $k\zeta$ intersects the curve $2\kappa u_\tau / |\langle u\rangle (\zeta ) - c|$ . This graphical procedure is repeated for each wave-age case to determine the corresponding $\zeta _{\textit{inner}}$ , as illustrated in figure 17. We note that if an explicit model for the mean velocity profile were available for the different wave-age cases, then $\zeta _{\textit{inner}}$ could be obtained iteratively from (5.5). The $c/u_\tau =-25$ case has the smallest $k\zeta _{\textit{inner}}$ , $\mathcal{O}(10^{-2})$ ( $\sim\! 2.5$ viscous units), which increases to $k\zeta _{\textit{inner}}\sim \mathcal{O}(1)$ ( $\sim\! 250$ viscous units) for the $c/u_\tau =15$ case and decreases again for the $c/u_\tau =25$ case with $k\zeta _{\textit{inner}}\sim 0.05$ ( $\sim\! 12.5$ viscous units). Correspondingly, the magnitude of $|{\hat {\nu }_{T,\textit{ij}}}|$ in figure 16(a–h) follows this trend for $\zeta _{\textit{inner}}$ . The magnitude of the velocity difference $|\langle u\rangle (\zeta _{\textit{inner}})-c|$ is larger for $c/u_\tau =-25$ , which results in a smaller $\zeta _{\textit{inner}}$ . The velocity difference becomes progressively smaller from the $c/u_\tau =-25$ to $c/u_\tau =15$ cases, leading to an increase in the inner-layer thickness. When $c/u_\tau$ further increases from 15 to 25, the wave phase speed surpasses the wind speed and the velocity difference $|\langle u\rangle (\zeta _{\textit{inner}})-c|$ increases. As a result, the $\zeta _{\textit{inner}}$ value for $c/u_\tau = 25$ is reduced relative to $c/u_\tau = 15$ and becomes comparable to that for $c/u_\tau = -7$ .

Figure 18. (a–h) Vertical profiles of the normalised eddy viscosity magnitude $|{\hat {\nu }_{T,\textit{ij}}}|/(m_{\textit{ij}} u_\tau \kappa \zeta _{\textit{inner}})$ , where $m_{\textit{ij}}$ is the proportionality constant for the different turbulence stress components to ensure that $|{\hat {\nu }_{T,\textit{ij}}}|/(m_{\textit{ij}} u_\tau \kappa \zeta _{\textit{inner}})\sim \mathcal{O}(1)$ . (i) Dependence of the proportionality constant $m_{\textit{ij}}$ on the wave age.

The variation in the critical-layer thickness ( $\zeta _{\textit{critical}}$ ), which is defined as the height where the mean-flow velocity equals the wave phase speed ( $|\langle u\rangle (\zeta _{\textit{critical}})-c|=0$ ) (Miles Reference Miles1957) is also depicted in figure 17. A critical layer is present only in the $c/u_\tau =2$ , $7$ and $15$ cases and it resides within the inner layer. Importantly, the existence of this critical layer does not alter the upper limit of the size of wave-induced turbulent eddies, as the condition $T_{\textit{eddy}} \lt T_{\textit{ad}v}$ remains valid at $\zeta _{\textit{critical}}$ . Consequently, for these wind-following wave cases, the wave-induced turbulent eddies are expected to grow to a size of the order of $\zeta \sim \zeta _{\textit{inner}}$ .

The validity of the $|{\hat {\nu }_{T,\textit{ij}}}|\sim u_\tau \zeta _{\textit{inner}}$ scaling is demonstrated in figure 18(a–h), which presents $|{\hat {\nu }_{T,\textit{ij}}}|$ normalised by $m_{\textit{ij}}u_\tau \zeta _{\textit{inner}}$ using an appropriate choice of proportionality constant $m_{\textit{ij}}$ . Focusing on the inner region, figure 18(a–d) displays plots of normalised $|{\hat {\nu }_{T,\textit{ij}}}|$ against $\zeta ^+=\zeta u_\tau /\nu$ on the vertical axis. Similarly, figure 18(e–h) shows the same normalised eddy viscosity magnitudes, but with $k\zeta$ on the vertical axis, which highlights the behaviour of these eddy viscosities in the outer region. From figure 18(a–h), we observe that the resulting normalised eddy viscosity profiles are approximately self-similar. This self-similarity in eddy viscosity profiles confirms that $u_\tau \zeta _{\textit{inner}}$ is indeed an appropriate scaling for the wave-induced turbulence eddy viscosity. The variation of $m_{\textit{ij}}$ across different wave ages is shown in figure 18(i). The proportionality constant $m_{\textit{ij}}$ is evaluated from the LES data as $m_{\textit{ij}} = \langle |{\hat {\nu }_{T,\textit{ij}}}| / (u_\tau \zeta _{\textit{inner}}) \rangle _{\zeta ^+ \gt 100}$ . Here, $\langle \boldsymbol{\cdot }\rangle _{\zeta ^+ \gt 100}$ denotes the mean value taken in the region $\zeta ^+ \gt 100$ . Since $\hat {\nu }_{T,ij}$ is approximately constant above $\zeta ^+ \sim 100$ (as shown in figure 17 a–d, this averaging yields a single representative value of the proportionality constant $m_{\textit{ij}}$ such that $|{\hat {\nu }_{T,\textit{ij}}}|/(m_{\textit{ij}}u_\tau \zeta _{\textit{inner}})\sim 1$ for $\zeta ^+\gt 100$ . The values of $m_{11}, m_{33}, m_{13}, m_{31}$ remain constant for wind over wave flows with wave ages ranging from $c/u_\tau =-25$ to $7$ , but differ for the $c/u_\tau =15$ and $25$ cases. The difference observed in the last two cases can be attributed to the attenuation of wave-induced perturbations at $c/u_\tau =15$ and $25$ , which is due to a decrease in wave-induced turbulent energy production that is discussed below in § 6. The reduction in wave-coherent energy production alters the behaviour of the wave-induced turbulence for the $c/u_\tau =15$ and $25$ cases compared with other wave-age scenarios, which results in changes to the proportionality constant $m_{\textit{ij}}$ for these two cases.

To summarise § 5, the wave-induced turbulence eddy viscosity differs significantly from the mean-flow-based eddy viscosity. This distinction is the primary reason for the underperformance of the turbulence curvilinear model equations when the mean-flow-based eddy viscosity is used. The eddy viscosity model for wave-induced Reynolds stresses must capture the complex characteristics of the ${\hat {\nu }_{T,\textit{ij}}}=|{\hat {\nu }_{T,\textit{ij}}}|\mathrm{e}^{\text{i}\phi _{\textit{ij}}(\zeta )}$ profiles presented in figure 16. Therefore, we conclude that it is inappropriate to use (5.1), (5.2) and (5.3) to solve (2.19), as $\nu _{T,{I}}$ , $\nu _{T,\textit{II}}$ and $\nu _{T,\textit{III}}$ do not correctly represent the behaviour of wave-induced Reynolds stresses. The $|{\hat {\nu }_{T,\textit{ij}}}|\sim u_\tau \zeta _{\textit{inner}}$ scaling provides valuable information for developing an appropriate eddy viscosity model. Nevertheless, the complex vertical variation of the phase angle $\phi _{\textit{ij}}(\zeta )$ , as illustrated in figure 2.18(i–l), requires further investigation.

Advanced turbulence modelling approaches, such as the Reynolds stress transport model (Launder, Reece & Rodi Reference Launder, Reece and Rodi1975), the macroscopic forcing method (MFM)-based turbulence model (Mani & Park Reference Mani and Park2021) and Pope’s nonlinear eddy viscosity model (Pope Reference Pope1975) implemented within a machine learning framework (Ling, Kurzawski & Templeton Reference Ling, Kurzawski and Templeton2016), may provide powerful means of capturing the complex, anisotropic nature of turbulence stresses in wind–wave interactions. The Reynolds stress transport model directly solves the transport equations for the Reynolds stresses, whereas the MFM-based model yields a tensorial eddy viscosity linking the mean turbulent fluxes to the local strain-rate tensor. Pope’s nonlinear eddy viscosity model expresses the Reynolds stress as a linear expansion of ten tensorial basis functions that depend on the strain- and rotation-rate tensors, for which the coefficients may be evaluated via a tensor-based neural network (TBNN) framework (Ling et al. Reference Ling, Kurzawski and Templeton2016). All these models account for turbulence anisotropy, which is critical for accurately capturing wave-induced fluctuations.

Previous studies (Harris et al. Reference Harris, Belcher and Street1996; Mastenbroek et al. Reference Mastenbroek, Makin, Garat and Giovanangeli1996; Meirink & Makin Reference Meirink and Makin2000) have demonstrated that solving the coupled continuity, momentum and Reynolds stress transport equations with mixing-length or two-equation-type (e.g. $k$ – $\epsilon$ , $k$ – $\omega$ ) turbulence closures can reproduce wave-induced perturbations reasonably well in comparison with experimental observations (e.g. Stewart Reference Stewart1970). While such models may provide a comprehensive representation of wave-induced turbulence, their implementation lies beyond the scope of the present study. The focus here is on assessing the effectiveness of Prandtl-type mixing-length-based eddy viscosity models in predicting wave-induced fluctuations. With suitable turbulence parameterisation, these simplified closures can serve as computationally efficient tools for quantifying wind–wave interactions without solving the full nonlinear Reynolds stress transport equations. From this modelling perspective, MFM- or TBNN-based turbulence parameterisation frameworks may offer cost-effective alternatives to the Reynolds transport approach. Extending the present turbulence curvilinear model framework to incorporate such advanced closure strategies will be pursued in our future work.

Next, if one compares wind over wave flows with turbulent boundary layer flows over rough walls, wave-induced fluctuations correspond to dispersive fluctuations arising from surface heterogeneity. Likewise, wave-induced turbulence stresses are conceptually similar to dispersive Reynolds stresses. Studies by Meyers, Ganapathisubramani & Cal (Reference Meyers, Ganapathisubramani and Cal2019) and Abdeldayem et al. (Reference Abdeldayem, Bon, Cal and Meyers2025) have derived linearised models which are similar in spirit to the turbulence curvilinear model to predict dispersive fluctuations in flows over heterogeneous surfaces. In these studies, the dispersive Reynolds stresses were represented using a mean-flow-based eddy viscosity, which captured weak dispersive fluctuations with good accuracy but achieved only approximately 55 % accuracy for strong dispersive fluctuations compared with DNS results (Abdeldayem et al. Reference Abdeldayem, Bon, Cal and Meyers2025). This result highlights the need for improved turbulence closures.

In summary, the turbulence associated with the mean flow and that arising from wave-coherent motions exhibit fundamentally different characteristics. As such, the wave-coherent turbulence cannot be accurately represented using mean-flow-based eddy viscosity closures, and advanced turbulence modelling strategies are required to capture its flow physics.