2.1. Numerical simulation of turbulent flow in a moving frame

The key to understanding wind-driven waves is determining wind pressure and tangential viscous stress. To anticipate the wind forces, we first examine the turbulent flow over a prescribed moving wavy boundary. The wave profile is described as $\eta (x,y,t)=a\cos (k(x-ct))$ , where $a$ is the wave amplitude, $k={2\pi }/{{{\lambda }_{x}}}$ stands for the wavenumber ( ${\lambda }_{x}$ is the wavelength in the direction of wave propagation) and $c$ is the phase velocity. A moving frame is introduced by incorporating the mesh motion equation into the three-dimensional incompressible Navier–Stokes equations:

(2.1) \begin{align} \frac {\partial ({{u}_{i}}-{{u}_{w,i}})}{\partial {{x}_{i}}}=0, \\[-28pt] \nonumber \end{align}

(2.2) \begin{align} \frac {\partial {{u}_{i}}}{\partial t}\bigg | _{\vartheta }+({{u}_{\kern-1pt j}}-{{u}_{w,j}})\frac {\partial {{u}_{i}}}{\partial {{x}_{\kern-1pt j}}}=-\frac {1}{\rho }\frac {\partial p}{\partial {{x}_{i}}}+\nu \frac {{{\partial }^{2}}{{u}_{i}}}{\partial {{x}_{\kern-1pt j}}\partial {{x}_{\kern-1pt j}}}+\varPi {{\delta }_{1i}}, \\[0pt] \nonumber \end{align}

where $({{x}_{1}},{{x}_{2}},{{x}_{3}})=(x,y,z)$ (with tensor notation $i$ the $i$ th component) respectively denote the streamwise, spanwise and vertical coordinates, $({{u}_{1}},{{u}_{2}},{{u}_{3}})=(u,v,w)$ represent the corresponding velocity components and repeated indices are implicitly summed over. Here, $p$ is the pressure, $\varPi$ is the external force transformed as a pressure gradient driving the flow, ${\delta }_{1i}$ is the Kronecker delta with ${{\delta }_{11}}=1$ and ${{\delta }_{12}}={{\delta }_{13}}=0$ , $\nu$ is the kinematic viscosity and $\rho$ is the fluid density. In contrast to equations in a stationary coordinate frame, mesh motion introduces a time-dependent change in the velocity within the new coordinate frame $\vartheta$ . This change is independent of both the Lagrange and Euler coordinate systems. Furthermore, the mesh motion adds the velocity vector of the moving coordinate frame ${{\boldsymbol{u}}_{w}}=({{u}_{w,1}},{{u}_{w,2}},{{u}_{w,3}})$ to the nonlinear terms. The motion of the moving mesh can be described by the diffusion equation, which governs the grid point displacement:

(2.3) \begin{align} K\frac {{{\partial }^{2}}{{D}_{w,j}}}{\partial {{x}_{i}}\partial {{x}_{i}}}=0,\quad \text{for}\quad j=1,2,3, \end{align}

where $({{D}_{w,1}},{{D}_{w,2}},{{D}_{w,3}})$ represent the displacement components for the $w$ th grid point in the streamwise, spanwise and vertical directions, respectively. The diffusive coefficient, denoted as $K$ , is related to the position of the grid point and is defined as $K=l^{-w}$ , where $l$ is the average distance from the wall to each grid point, specifically the inverse distance. Given the moving wave motion at the wall boundary, the velocity vector of the moving grid point can be expressed as ${{u}_{w,i}}={(D_{w,i}^{t+\Delta t}-D_{w,i}^{t})}/{\Delta t}$ . The mesh velocity is also incorporated into (2.1). By rewriting (2.2) in the integral form, we obtain

(2.4) \begin{align} \frac {{\rm d}}{{\rm d}t}\int _{V}{\rho \boldsymbol{u}\,{\rm d}V}+\oint _{S}{\rho \boldsymbol{n}\boldsymbol{\cdot }(\boldsymbol{u}-{{\boldsymbol{u}}_{w}})\boldsymbol{u}\,{\rm d}S}=-\oint _{S}{\rho \boldsymbol{n}\boldsymbol{\cdot }{{\boldsymbol{Q}}_{\boldsymbol{u}}}\,{\rm d}S}+\int _{V}{{{S}_{\boldsymbol{u}}}\,{\rm d}V}, \end{align}

where $V$ is the control volume and $S$ is the closed surface bounding the volume, $\boldsymbol{u}$ is the fluid velocity vector and ${\boldsymbol{u}}_{w}$ is the velocity vector of the control volume, ${\boldsymbol{Q}}_{\boldsymbol{u}}$ is the volume source and ${S}_{\boldsymbol{u}}$ is the surface source and $\boldsymbol{n}$ is the unit normal vector pointing outward on the boundary surface of the volume (Jasak & Tukovic Reference Jasak and Tukovic2006). Applying a second-order finite-volume discretisation to (2.4) yields

(2.5) \begin{align} \frac {{{({{\rho }_{P}}{{\boldsymbol{u}}_{P}}{{V}_{P}})}^{t+\vartriangle t}}-{{({{\rho }_{P}}{{\boldsymbol{u}}_{P}}{{V}_{P}})}^{t}}}{\vartriangle t}+\sum \limits _{\!f}{{{\rho }_{\!f}}(F-{{F}_{S}}){{\boldsymbol{u}}_{\!f}}}=-\sum \limits _{\!f}{{{\boldsymbol{S}}_{\!f}}\boldsymbol{\cdot }\rho {{\boldsymbol{Q}}_{\boldsymbol{u}}}}+{{S}_{u}}{{V}_{P}}, \end{align}

where the subscript $P$ represents the grid point, the subscript $f$ denotes the value on the face, $V_P$ is the volume of grid $P$ , $F={{\boldsymbol{S}}_{\!f}}\boldsymbol{\cdot }{{\boldsymbol{u}}_{\!f}}$ stands for the fluid flux with ${{\boldsymbol{S}}_{\!f}}=\boldsymbol{n}{{S}_{\!f}}$ and $F_S$ represents the mesh motion flux. In contrast to the momentum equation, which does not account for mesh motion, a new term, the ‘mesh motion flux’, has been introduced. This term must adhere to the law of space conservation, which implies that the mesh motion flux is determined by the rate of change in grid volume:

(2.6) \begin{align} \frac {{\rm d}}{{{\rm d}}t}\int _{V}{{{\rm d}}V}-\oint _{S}{\boldsymbol{n}\boldsymbol{\cdot }{{\boldsymbol{u}}_{w}}\,{{\rm d}}S=0}\quad \Longrightarrow \quad \frac {V_{P}^{t+\Delta t}-V_{P}^{t}}{\Delta t}-\sum \limits _{\!f}{{{F}_{S}}}=0, \end{align}

where the meanings of the symbols are consistent with those in (2.4) and (2.5). It is important to ensure that (2.6) is satisfied at each time step. Additionally, the mesh motion flux $F_S$ is determined by the volume swept by the moving surface during the current time step rather than by the mesh velocity ${\boldsymbol{u}}_{w}$ . This approach helps to prevent errors that could arise from interpolating the mesh motion velocity from the grid centre to the face centre. For more details, see Appendix A. It is also worth mentioning that significant contributions to the numerical simulations of turbulent flow over a moving wavy boundary have been made by Lian Shen and his collaborators, as well as Peter Sullivan and his group, through the use of curvilinear transformations (Sullivan et al. Reference Sullivan, McWilliams and Moeng2000, Reference Sullivan, Edson, Hristov and McWilliams2008; Shen et al. Reference Shen, Zhang, Yue and Triantafyllou2003; Yang & Shen Reference Yang and Shen2010, Reference Yang and Shen2017; Husain et al. Reference Husain, Buckley, Yousefi, Veron and Sullivan2019; Cao & Shen Reference Cao and Shen2021). This transformation to a wave-following (curvilinear) coordinate system offers several advantages for problems involving wave surfaces. It ensures conservation properties, accuracy of fluxes and proper representation of boundary conditions. In this study, we utilise the finite-volume method with mesh motion, along with the OpenFOAM solver, to perform numerical simulations. The numerical framework incorporates an arbitrary Lagrangian–Eulerian description within its moving-mesh module. This description reformulates the governing equations into a conservative form relative to mesh motion. As noted by Jasak & Tukovic (Reference Jasak and Tukovic2006), the spatial conservation of fluxes due to mesh motion must be maintained. In the moving-mesh method, the computational mesh is adjusted to conform to the boundary shape, which is updated at each time step of the transient simulation. The prescribed motion of the boundaries dictates the movement of all points within the mesh. However, the primary challenge in handling variable-geometry cases is maintaining the quality of the mesh. Poor mesh quality may introduce errors, limiting the handling of complex, largely deformed mesh motions. As noted by Tukovic & Jasak (Reference Tukovic and Jasak2012), a significant drawback of moving-mesh methods is that they are usually restricted to interfaces that undergo only moderate deformation and maintain constant topology. In cases where large interface deformations or changes in topology are required, some form of global or local remeshing is necessary. Since this study focuses on small-amplitude moving waves, we provide a detailed description of the numerical scheme and the handling of the deformed moving mesh in Appendix A. Furthermore, we include a careful comparison with previous numerical results of turbulent flow over moving wave boundaries, using the curvilinear transformation, to ensure the validity of this method for the current problem. Specifically, we solve the full Navier–Stokes equations in a moving reference frame without any additional modelling. To accurately capture the flow process, we aim to enhance the grid resolution to an acceptable level. In the current solver, the convection terms in the momentum equations are discretised using a second-order upwind difference scheme to account for the velocity gradient (see Appendix A for more details). Furthermore, a second-order backwards implicit scheme is employed for time stepping. The finite-volume method, based on the standard pressure–velocity coupling approach and the PIMPLE algorithm, is used to solve the governing equations (Gatin et al. Reference Gatin, Vukcevic, Jasak and Rusche2017; Gimenez et al. Reference Gimenez, Aguerre, Idelsohn and Nigro2019). Before solving the pressure–velocity coupled equations, we correct the face fluxes. After correcting the pressure, the face velocity vector is updated accordingly. These two steps ensure flux consistency throughout the simulation.

2.2. Simulation configuration

We applied an external force to drive turbulent flow, which was effectively converted into a pressure gradient to maintain a bulk velocity of

(2.7) \begin{align} {{U}_{0}}=\frac {\iint {U{\rm d}x\,{\rm d}z}}{\iint {{\rm d}x\,{\rm d}z}}, \end{align}

where $U$ represents the instantaneous velocity. The Reynolds number was determined to be ${Re}={{{U}_{0}}{{\lambda }_{x}}}/{\nu }\approx 7300$ , with the corresponding friction Reynolds numbers $ \textit{Re}_{\tau }={{{u}_{\tau }}{{\lambda }_{x}}}/{\nu }$ presented in table 1. To calculate the friction velocity $u_{\tau }$ , we followed the methodology outlined by Bopp (Reference Bopp2018), which involves extrapolating turbulent shear stress from the outer layer. This approach is based on the principle that total wind stress can be expressed as the sum of turbulent, wave-coherent and viscous stresses in the outer layer. It allows us to isolate the turbulent stress, assuming that the wave-coherent and viscous stresses in the outer layer are negligible. In this study, we considered five wave slopes, $ak=0.12$ , 0.13, 0.15, 0.2 and 0.25, which meet the criteria for small-amplitude water waves. The wave age ranges from ${c}/{{{u}_{\tau }}}=1.00$ to ${c}/{{{u}_{\tau }}}=3.89$ . The size of the computational domain was defined as $ (2{{\lambda }_{x}},{{\lambda }_{x}},0.5{{\lambda }_{x}} )$ , which has been verified to be sufficiently large based on the two-point spatial correlation. The total number of grid points is ${{N}_{x}}\times {{N}_{y}}\times {{N}_{z}}=251\times 126\times 76$ , with dimensionless grid scales using the wave-following curvilinear coordinates $(\xi ,\zeta )=(x,( {z-\eta }/{H_h}))$ , where $\eta$ represents the wave elevation and $H_h$ denotes the height of the physical domain, as shown in table 1. We note that Komen et al. (Reference Komen, Camilo, Shams, Geurts and Koren2017) conducted tests on turbulent channel flow. Their results, which employed a 9-grid resolution, were found to be consistent with previous direct numerical simulation (DNS) findings for mean velocity, root-mean-square velocity and the budget terms of mean turbulent kinetic energy. The dimensionless grid scales in the streamwise, spanwise and vertical directions are given by $(\Delta {{x}^{+}},\Delta {{y}^{+}},\Delta {{z}^{+}})=(9,4.5,4.5)$ , with the first layer in the vertical direction set at $\varDelta z_{w\textit{all}}^{+}=0.25$ , where the ‘+’ sign indicates normalisation by the viscous scale. The current grid scales, as shown in table 1, are defined by $(\Delta {{\xi }^{+}},\Delta {{y}^{+}},\Delta \zeta _{b}^{+})\sim (5,5,0.5)$ , with the subscript ‘ $b$ ’ referring to the first layer of the bottom boundary, similar to the grid resolution used by Komen et al. (Reference Komen, Camilo, Shams, Geurts and Koren2017). The Kolmogorov scale ${\eta }_{K}$ is defined by ${{\eta }_{K}}={{\nu }^{{3}/{4}}}{{\epsilon }^{{-1}/{4}}}$ , where $\epsilon$ is the dissipation of turbulent kinetic energy. Assuming the wavelength, $\lambda _x$ , as a characteristic length scale, we can approximate $\epsilon \approx {u_{\tau }^{3}}/{\lambda _x}$ . With the current settings of ${{\lambda }_{x}}=0.1$ , ${{u}_{\tau }}\approx 0.0059{-} 0.0078$ and $\nu =1.1\times {{10}^{-6}}$ , we find that ${{\eta }_{K}}\approx 0.0007{-}0.0009$ , leading to $\varDelta \xi =\Delta y=0.0008\approx 0.89{{\eta }_{K}}{-} 1.14{{\eta }_{K}}$ . Since the grid points in the vertical direction are concentrated near the bottom boundary, the maximum grid spacing occurs at the upper boundary, denoted as $\varDelta {{\zeta }_{u}}\approx 0.0018\approx 2{{\eta }_{K}}{-} 2.57{{\eta }_{K}}$ , where the subscript ‘ $u$ ’ refers to the grid at the upper boundary. However, by considering the wavelength $\lambda _x$ as a characteristic length scale, we acknowledge that the turbulent flow is mainly influenced by wave motion. In scenarios involving turbulent flow over a slowly moving wavy boundary, such as those driven by strong wind forces, we reference the experimental findings by Yousefi et al. (Reference Yousefi, Veron and Buckley2020b ), which indicate that typical turbulent structures are predominately found within a height of $kz={2\pi z}/{{{\lambda }_{x}}}\approx 1$ (or a height characterising the wind-wave boundary layer). Thus, we can define the characteristic length scale as ${{{\lambda }_{x}}}/{2\pi }$ . Upon re-evaluation using the same parameters, we find that ${{\eta }_{K}}\approx 0.00046{-}0.00056$ . Consequently, we have $\varDelta \xi =\Delta y=0.0008\approx 1.43{{\eta }_{K}}{-} 1.74{{\eta }_{K}}$ and the maximum grid space in the vertical direction is $\varDelta {{\zeta }_{u}}\approx 0.0018\approx 3.21{{\eta }_{K}}{-} 3.91{{\eta }_{K}}$ . The first layer along the vertical direction is $\varDelta {{\zeta }_{b}}\approx 0.0001\approx 0.179{{\eta }_{K}}{-} 0.218{{\eta }_{K}}$ . These approximations meet the typical second-order accuracy requirements for DNS of turbulent channel flow. However, we assert that the turbulent flow over a moving wavy boundary is more complicated than that of a turbulent channel flow. The estimation is based on the Kolmogorov scale approximation; therefore, we provide more detailed results in Appendix B, which outline the ratio of the grid size to the Kolmogorov scale. To validate the wind forces obtained from our numerical simulations, we compare our results with those from previous experiments and simulations, as discussed in §§ 2.3, 2.4 and 2.5. In addition to the grid convergence validation presented in § 2.3, we also include a comparison of results for the case where the wave slope is $ak=0.25$ across different mesh resolutions (refer to Appendix B). For our computations, periodic boundary conditions were imposed in both the streamwise and spanwise directions. The bottom wall was modelled with a wavy motion to represent a moving-wall boundary, while a slip-wall condition was utilised for the upper wall. Table 1 provides detailed parameters of the simulation configuration, along with some resolved quantities, including the dimensionless height of the wave boundary layer denoted as $\delta _{w}^{+}={{{u}_{\tau }}}/{k\nu }$ and the height of the Stokes layer given by ${{\delta }_{s}}=\sqrt {{2\nu }/{ck}}$ .

Table 1.

Parameter settings. We consider five wave slopes in this study: $ak=0.12, 0.13, 0.15, 0.2, 0.25$ , which fall within the slow-wave regime, corresponding to a wave age of $1.00{-} 3.89$ . The bulk Reynolds number is ${Re}={U_0}{\lambda _x}/\nu \approx 7300$ , where $U_0$ is the bulk velocity and $\lambda _x$ is the wavelength, and the corresponding friction Reynolds number is $ \textit{Re}_{\tau }={u_{\tau }}{\lambda _x}/{\nu }=455{-} 710$ , where $u_{\tau }$ is the friction velocity. The grid space normalised by the viscous scale in streamwise, spanwise and the first layer in the vertical direction is approximately $(\Delta {{\xi }^{+}},\Delta {{y}^{+}},\Delta \zeta _{b}^{+} )\approx (5, 5, 0.5)$ , where $(\xi ,\zeta )=(x,( {z-\eta }/{H_h}))$ represents the wave-following coordinate frame ( $\eta$ is the wave elevation and $H_h$ is the height of the physical domain). The normalised thickness of the wave boundary layer $\delta _{w}^{+}={{{u}_{\tau }}}/{k\nu }$ and the Stokes layer ${{{\delta }_{s}}{{u}_{\tau }}}/{\nu }$ are also shown, where ${{\delta }_{s}}=\sqrt {{2\nu }/{ck}}$ .

Table 1. Long description

A table comparing parameter settings for different cases. The table has 11 rows and 7 columns. The columns are labeled as case, ak, c/uτ, (Δξ⁺, Δy⁺, Δz⁺b), δ⁺w, δs uτ/v, and Reτ. The rows are labeled with case identifiers from S1 to S11. Each row provides specific values for the parameters listed in the columns. Row 1: case S1, ak 0.12, c/uτ 3.89, (Δξ⁺, Δy⁺, Δz⁺b) (4.30, 4.30, 0.54), δ⁺w 85.60, δs uτ/v 6.64, Reτ 538. Row 2: case S2, ak 0.12, c/uτ 3.82, (Δξ⁺, Δy⁺, Δz⁺b) (4.30, 4.30, 0.54), δ⁺w 85.60, δs uτ/v 6.70, Reτ 538. Row 3: case S3, ak 0.13, c/uτ 3.69, (Δξ⁺, Δy⁺, Δz⁺b) (4.45, 4.45, 0.55), δ⁺w 88.61, δs uτ/v 6.93, Reτ 557. Row 4: case S4, ak 0.13, c/uτ 2.37, (Δξ⁺, Δy⁺, Δz⁺b) (4.45, 4.45, 0.56), δ⁺w 88.61, δs uτ/v 8.65, Reτ 557. Row 5: case S5, ak 0.15, c/uτ 3.43, (Δξ⁺, Δy⁺, Δz⁺b) (4.79, 4.79, 0.59), δ⁺w 95.42, δs uτ/v 7.46, Reτ 600. Row 6: case S6, ak 0.15, c/uτ 2.00, (Δξ⁺, Δy⁺, Δz⁺b) (3.64, 3.64, 0.45), δ⁺w 72.34, δs uτ/v 8.51, Reτ 455. Row 7: case S7, ak 0.15, c/uτ 1.00, (Δξ⁺, Δy⁺, Δz⁺b) (3.64, 3.64, 0.45), δ⁺w 72.34, δs uτ/v 12.03, Reτ 455. Row 8: case S8, ak 0.20, c/uτ 3.28, (Δξ⁺, Δy⁺, Δz⁺b) (5.01, 5.01, 0.63), δ⁺w 99.72, δs uτ/v 7.80, Reτ 627. Row 9: case S9, ak 0.25, c/uτ 2.89, (Δξ⁺, Δy⁺, Δz⁺b) (5.68, 5.68, 0.71), δ⁺w 113.00, δs uτ/v 8.84, Reτ 710. Row 10: case S10, ak 0.25, c/uτ 2.01, (Δξ⁺, Δy⁺, Δz⁺b) (5.63, 5.63, 0.70), δ⁺w 112.07, δs uτ/v 10.55, Reτ 701. Row 11: case S11, ak 0.25, c/uτ 1.53, (Δξ⁺, Δy⁺, Δz⁺b) (4.77, 4.77, 0.60), δ⁺w 94.87, δs uτ/v 11.15, Reτ 596.

2.3. Numerical verification

Wave motion plays a significant role in complex momentum transport alongside turbulent winds. To better describe wave motion, a wave-phase averaging method can be used to capture the disturbances caused by waves. The instantaneous velocity can be decomposed as (Hussain & Reynolds Reference Hussain and Reynolds1970)

(2.8) \begin{align} u\left ( x,z,t \right )=\langle u\left ( \varphi ,z \right )\rangle +{u}'\left ( x,z,t \right )=\overline {u}(z)+\tilde {u}\left ( \varphi ,z \right )+{u}'\left ( x,z,t \right ), \end{align}

with

(2.9) \begin{align} \langle u\left ( \varphi ,z \right )\rangle =\underset {N\to \infty }{\mathop {\lim }}\,\frac {1}{N}\sum \limits _{n=0}^{N-1}{u\left ( \varphi +2n\pi ,z \right )}, \\[-28pt] \nonumber \end{align}

(2.10) \begin{align} \overline {u}(z)=\underset {M\to \infty }{\mathop {\lim }}\,\frac {1}{M}\sum \limits _{m=0}^{M-1}\left\langle {u\left ( \varphi +\frac {2\pi m}{M},z \right )}\right\rangle \! , \\[0pt] \nonumber \end{align}

where $\tilde {u} (\varphi , z )$ is the wave-induced velocity (here $\varphi$ is the phase), determined by the difference between phase-averaged velocity $\langle u (\varphi , z )\rangle$ and ensemble-averaged velocity $\overline {u} ( z )$ , and ${u}' ( x,z,t )$ represents the turbulent fluctuating velocity.

To ensure that the computational domain is sufficiently large, we evaluate the two-point spatial correlation coefficient. The time series of velocity fluctuations is collected after the turbulence has developed. Figure 2 illustrates that the two-point spatial correlation coefficients decrease to negligibly small values, which indicates that the domain is adequately sized. Additionally, in our previous work on turbulent channel flow over stationary two- and three-dimensional wavy walls under the same Reynolds number conditions (Zhang et al. Reference Zhang, Wang and Liu2021, Reference Zhang, Wang and Liu2022a ,Reference Zhang, Wu, Liu and Wang b , Reference Zhang, Wang, Wu, Wang and Liu2023), we established the same domain dimensions in the streamwise and spanwise directions. The correlation coefficients in those studies were also found to be negligibly small, further supporting the adequacy of the domain size used in this research. In addition, we expanded the computation domain from the original dimensions of $ ( x,y,z )= ( 2{{\lambda }_{x}},{{\lambda }_{x}},0.5{{\lambda }_{x}} )$ to the enlarged dimensions of $ ( x,y,z )= ( 4{{\lambda }_{x}},2{{\lambda }_{x}},0.5{{\lambda }_{x}} )$ . Figure 3 compares turbulent statistics between these two cases, including phase-averaged streamwise and vertical velocities, turbulent kinetic energy and Reynolds shear stress. Our analysis reveals no significant differences between the two domains. The averaged flow remains confined to the scale of one wavelength. We also plotted the vertical profiles of these statistics, as shown in figure 4. The mean velocity, turbulent kinetic energy and Reynolds shear stress for the original domain closely match those of the enlarged domain. Therefore, we conclude that the domain size is large enough for the current settings. It is noteworthy that elongated streaks and secondary motions are intrinsic features of wind-wave dynamics, significantly influencing the momentum and energy exchange at the air–sea interface. As emphasised in § 3, the present study focuses on developing a nonlinear wind-wave model within a two-dimensional framework. To bridge this dimensional gap, the generalised wind-force model is derived by spanwise-averaging the three-dimensional numerical results. A detailed discussion on the role of three-dimensional vortical structures in turbulent momentum transport is provided in Appendix C. Nevertheless, while our wind-forcing formulations implicitly account for the mean effects of three-dimensional turbulence, the subsequent wave evolution in §§ 3 and 5 remains constrained to a two-dimensional framework. Consequently, the model does not resolve three-dimensional wave dynamics, such as oblique modes and directional spreading, thus representing an idealised mechanistic scenario rather than fully developed three-dimensional wind-wave turbulence. Despite these inherent limitations, the validations in this section – along with those in §§ 2.4 and 2.5 – substantiate the rationality and reliability of the proposed wind-force model within the defined scope of this research.

Figure 2.

The two-point velocity spatial correlation coefficient at ${\zeta {{u}_{\tau }}}/{\nu }\approx 80$ along the streamwise (a) and spanwise (b) directions.

Figure 3.

Comparisons of turbulent statistics, including (a,e) phase-averaged streamwise velocity, (b,f) phase-averaged vertical velocity, (c,g) turbulent kinetic energy and (d,h) Reynolds shear stress, are made between the enlargement domain case (a–d) and the original domain case (e–h).

Figure 3. Long description

Panel A: A heat map representing the phase-averaged streamwise velocity for the enlargement domain case. The x-axis is labeled x/λx and the y-axis is labeled kz. The color scale ranges from 0 to 1.2, with lighter colors indicating higher velocities. The heat map shows a gradient with higher velocities at the top and lower velocities at the bottom, following the contour of the wavy boundary. Panel B: A heat map representing the phase-averaged vertical velocity for the enlargement domain case. The x-axis is labeled x/λx and the y-axis is labeled kz. The color scale ranges from -0.1 to 0.1, with lighter colors indicating positive velocities and darker colors indicating negative velocities. The heat map shows alternating regions of positive and negative velocities along the wavy boundary. Panel C: A heat map representing the turbulent kinetic energy for the enlargement domain case. The x-axis is labeled x/λx and the y-axis is labeled kz. The color scale ranges from 0 to 0.05, with lighter colors indicating higher kinetic energy. The heat map shows regions of higher kinetic energy near the peaks and troughs of the wavy boundary. Panel D: A heat map representing the Reynolds shear stress for the enlargement domain case. The x-axis is labeled x/λx and the y-axis is labeled kz. The color scale ranges from 0 to 0.012, with lighter colors indicating higher shear stress. The heat map shows regions of higher shear stress near the peaks and troughs of the wavy boundary. Panel E: A heat map representing the phase-averaged streamwise velocity for the original domain case. The x-axis is labeled x/λx and the y-axis is labeled kz. The color scale ranges from 0 to 1.2, with lighter colors indicating higher velocities. The heat map shows a gradient with higher velocities at the top and lower velocities at the bottom, following the contour of the wavy boundary. Panel F: A heat map representing the phase-averaged vertical velocity for the original domain case. The x-axis is labeled x/λx and the y-axis is labeled kz. The color scale ranges from -0.1 to 0.1, with lighter colors indicating positive velocities and darker colors indicating negative velocities. The heat map shows alternating regions of positive and negative velocities along the wavy boundary. Panel G: A heat map representing the turbulent kinetic energy for the original domain case. The x-axis is labeled x/λx and the y-axis is labeled kz. The color scale ranges from 0 to 0.05, with lighter colors indicating higher kinetic energy. The heat map shows regions of higher kinetic energy near the peaks and troughs of the wavy boundary. Panel H: A heat map representing the Reynolds shear stress for the original domain case. The x-axis is labeled x/λx and the y-axis is labeled kz. The color scale ranges from 0 to 0.012, with lighter colors indicating higher shear stress. The heat map shows regions of higher shear stress near the peaks and troughs of the wavy boundary.

Figure 4.

Comparisons of (a) mean velocity, (b) turbulent kinetic energy and (c) Reynolds shear stress profiles are made between the enlargement domain case (black line with circle symbols) and the original domain case (red line).

Next, we verify grid convergence by comparing the mean streamwise velocity, turbulent kinetic energy (TKE) and Reynolds shear stress (RSS) (normalised by the friction velocity) between the high-grid-resolution (HGR) and super-high-grid-resolution (SHGR) cases. In both cases, we set $ak=0.13$ and ${c}/{{{u}_{\tau }}=2.37}$ . The number of grid points for HGR is ${{N}_{x}}\times {{N}_{y}}\times {{N}_{z}}=251\times 126\times 76$ , while for SHGR it is ${{N}_{x}}\times {{N}_{y}}\times {{N}_{z}}=401\times 201\times 151$ . The normalised grid scales by viscous scale are $(\Delta {{\xi }^{+}},\Delta {{y}^{+}},\Delta \zeta _{b}^{+})=(4.45,4.45,0.55)$ for HGR and $(\Delta {{\xi }^{+}},\Delta {{y}^{+}},\Delta \zeta _{b}^{+})=(2.78,2.78,0.55)$ for SHGR. Figure 5 demonstrates that the velocity and high-order turbulent statistics for HGR align well with SHGR, indicating that the grid resolution for HGR has already achieved convergence.

Figure 5.

Grid convergence tests: a comparison between cases HGR and SHGR for (a) the ensemble-averaged streamwise velocity and (b) the ensemble-averaged TKE and RSS. The black lines are the results of HGR, while the red lines are the results of SHGR.

Finally, to verify the numerical model, we conducted simulations for cases corresponding to the experimental research by Yousefi et al. (Reference Yousefi, Veron and Buckley2020a ,Reference Yousefi, Veron and Buckley b ), referred to as case S3, and to the numerical investigations by Yang & Shen (Reference Yang and Shen2010, Reference Yang and Shen2017), referred to as case S10. Additionally, we included numerical studies (Sullivan et al. Reference Sullivan, McWilliams and Moeng2000; Akervik & Vartdal Reference Akervik and Vartdal2019; Zhang et al. Reference Zhang, Huang and Xu2019; Cao & Shen Reference Cao and Shen2021; Hao et al. Reference Hao, Cao and Shen2021; Esteghamatian et al. Reference Esteghamatian, Katz and Zaki2022) to validate the second-order statistics. Figure 6(a) illustrates a typical logarithmic law of velocity normalised by friction velocity, which is observed in the turbulent boundary layer near the wall. The velocity results from our simulations are in good agreement with both the experimental and numerical findings. Figure 6(b) presents the vertical profile of viscous stress, normalised by $\rho u_{\tau }^{2}$ , which aligns well with previous results of Yousefi et al. (Reference Yousefi, Veron and Buckley2020a ,Reference Yousefi, Veron and Buckley b ). In this context, the viscous stress is calculated using ${{\tau }_{\nu }}=\mu [ ({{u}_{z}}+{{w}_{x}} )+2{{\eta }_{x}} ({{w}_{z}}-{{u}_{x}} ) ]$ , where $\mu$ is the dynamic viscosity, and the subscripts $x$ and $z$ denote the gradient in the corresponding direction. In later discussions, we clarify that the viscous stress is normalised by $\rho u_{\tau }^{2}$ , which allows us to substitute dynamic viscosity $\mu$ with kinematic viscosity $\nu$ . It is important to note that in the numerical simulations, the pressure calculated is actually ${p}/{\rho }$ . This means we only need to specify the kinematic viscosity parameter, while the density parameter remains unspecified.

Figure 6.

Comparisons between the current study (indicated by red symbols and lines) and experimental results of Yousefi et al. (Reference Yousefi, Veron and Buckley2020a ,Reference Yousefi, Veron and Buckley b ) (EXP, $ak=0.13$ , $c/u_{\tau }=3.69$ ) as well as numerical outcomes for (a) streamwise velocity profiles, (b) viscous stress, (c) RSS, (d) wave-induced shear stress and (e,f) RSS profiles. The numerical results are from Yang & Shen (Reference Yang and Shen2017) (DNS1, $ak=0.25$ , $c/u_{\tau }=2$ ), Yang & Shen (Reference Yang and Shen2010) (DNS2, $ak=0.1$ , $c/u_{\tau }=2$ ), Sullivan et al. (Reference Sullivan, McWilliams and Moeng2000) (DNS3, $ak=0.1$ , $c/u_{\tau }=3.9$ ), Akervik & Vartdal (Reference Akervik and Vartdal2019) ( $ak=0.1$ , $c/u_{\tau }=4$ , $ \textit{Re}_{\tau }=200$ – rectangles, and $ \textit{Re}_{\tau }=950$ – circles), Cao & Shen (Reference Cao and Shen2021) ( $ak=0.15$ , $c/u_{\tau }=3.46$ and $ \textit{Re}_{\tau }=867$ ), Zhang et al. (Reference Zhang, Huang and Xu2019) ( $ak=0.1$ , $c/u_{\tau }=7.25$ and $ \textit{Re}_{\tau }\approx 3456$ ) and Hao et al. (Reference Hao, Cao and Shen2021) ( $ak=0.14$ , $c/u_{\tau }=7.7$ and $ \textit{Re}_{\tau }=2000$ ). We also present the results of turbulent channel flow over a hyperelastic-compliant wall as reported by Esteghamatian et al. (Reference Esteghamatian, Katz and Zaki2022).

Figure 6. Long description

Panel A: A line graph shows streamwise velocity profiles. The x-axis is labeled with u+ and ranges from 0 to 15. The y-axis is labeled with ζ+ and ranges from 10^-1 to 10^3. The graph includes data from the current study (red symbols and lines) and various experimental and numerical studies. Panel B: A line graph depicts viscous stress. The x-axis is labeled with u+ and ranges from 0 to 1. The y-axis is labeled with ζ+ and ranges from 10^0 to 10^3. The graph includes data from the current study and other studies. Panel C: A line graph shows Reynolds shear stress (RSS). The x-axis is labeled with x/λx and ranges from 0 to 2. The y-axis is labeled with -u’w’+ and ranges from 0 to 2. The graph includes data from the current study. Panel D: A line graph depicts wave-induced shear stress. The x-axis is labeled with kζ and ranges from 0 to 1. The y-axis is labeled with u+τ/w+ and ranges from -1 to 2. The graph includes data from the current study. Panel E: A line graph shows RSS profiles. The x-axis is labeled with ζ+ and ranges from 10^0 to 10^2. The y-axis is labeled with -u’w’+ and ranges from 0 to 2. The graph includes data from the current study and other studies. Panel F: A line graph depicts RSS profiles. The x-axis is labeled with ζ+ and ranges from 10^0 to 10^2. The y-axis is labeled with -u’w’+ and ranges from 0 to 2. The graph includes data from the current study and other studies.

Figure 6(c) shows the phase-averaged RSS along the waves at a height of ${\zeta }/{{{\lambda }_{x}}\approx 0.1}$ . Our results are consistent with those of Yang & Shen (Reference Yang and Shen2010). We also examined the vertical profiles of wave-induced shear stress, as illustrated in figure 6(d). Our findings align with previous experimental and numerical results (Sullivan et al. Reference Sullivan, McWilliams and Moeng2000; Yang & Shen Reference Yang and Shen2010; Yousefi et al. Reference Yousefi, Veron and Buckley2020b ). In figure 6(e), the vertical profile of RSS is displayed, with numerical results from Akervik & Vartdal (Reference Akervik and Vartdal2019), Zhang et al. (Reference Zhang, Huang and Xu2019), Hao et al. (Reference Hao, Cao and Shen2021) and Cao & Shen (Reference Cao and Shen2021) included for comparison. Although the parameters in these studies differ from those in the present research (as noted in the caption of figure 6), the vertical characteristics of the RSS profiles show notable similarities. Akervik & Vartdal (Reference Akervik and Vartdal2019) investigated flow driven by a pressure gradient with friction Reynolds numbers $ \textit{Re}_{\tau }=200$ and $ \textit{Re}_{\tau }=950$ , along with a wave age of $c/u_{\tau }=4$ . Our RSS results are comparable to their findings at $ \textit{Re}_{\tau }=200$ , including the peak RSS position. Additionally, Esteghamatian et al. (Reference Esteghamatian, Katz and Zaki2022) recently explored turbulent channel flow over a hyperelastic-compliant wall. The RSS results are in close alignment with those of Esteghamatian et al. (Reference Esteghamatian, Katz and Zaki2022) within the viscous sublayer, as illustrated in figure 6(f).

Notably, whereas Sullivan et al. (Reference Sullivan, McWilliams and Moeng2000) and Yang & Shen (Reference Yang and Shen2010, Reference Yang and Shen2017) performed numerical simulations of turbulent flow over moving wavy boundaries using curvilinear coordinate transformations, the present study utilises a dynamic mesh technique. To establish the fidelity of our approach against these distinct numerical formulations, rigorous cross-validations with previous experimental and numerical benchmarks are presented in §§ 2.4 and 2.5 (see figures 7–10). Furthermore, Appendix D details the wave-induced velocities and pressure for case S5, alongside the phase-averaged turbulent shear stress for case S10; these results capture spatial patterns and magnitudes consistent with those reported by Yang & Shen (Reference Yang and Shen2010) and Cao et al. (Reference Cao, Deng and Shen2020). Our recent work (Zhang et al. Reference Zhang, Wang and Liu2024) also explored the associated momentum flux and kinetic energy transfer, where the underlying numerical framework was further validated against canonical smooth-wall channel flow data (Lee & Moser Reference Lee and Moser2015, Reference Lee and Moser2018). Following these comprehensive assessments, encompassing domain-size independence (figures 2–4), grid convergence (figure 5 and Appendix B), algorithmic verification (Appendix A) and benchmark comparisons (figures 6–10 and Appendix D), the numerical model is demonstrated to be highly robust, providing a high-fidelity data asset for closed wind-force modelling.

Figure 7.

Mean form drags (a) and mean viscous stress (b) versus the wave slope. The data from Banner (Reference Banner1990), Sullivan et al. (Reference Sullivan, McWilliams and Moeng2000), Peirson et al. (Reference Peirson, Walker and Banner2014), Bopp (Reference Bopp2018), Sullivan et al. (Reference Sullivan, Banner, Morison and Peirson2018), Buckley et al. (Reference Buckley, Veron and Yousefi2020) and Funke et al. (Reference Funke, Buckley, Schultze, Veron, Timmermans and Carpenter2021) are also plotted.

2.4. Mean wind stress and wave growth

Experimental evidence from Zavadsky & Shemer (Reference Zavadsky and Shemer2012) suggests that wind stress decreases linearly with vertical height, a finding supported by Buckley et al. (Reference Buckley, Veron and Yousefi2020). According to Bopp (Reference Bopp2018), the total wind stress can be determined by extrapolating turbulent stress from the outer layer, as it matches the sum of turbulent, wave-coherent and viscous stresses in that region. However, Buckley et al. (Reference Buckley, Veron and Yousefi2020) pointed out that extrapolating turbulent stress profiles to the wave surface does not align with the wind stress estimated using mean velocity profiles based on the classical logarithmic law. To address this discrepancy, we present the mean form drag ( ${{\overline {\tau }}_{p}}=\overline {p{{\eta }_{x}}}$ ) and mean viscous stress ( ${\overline {\tau }}_{\nu }$ ), both normalised by $\rho u_{\tau }^{2}$ , in figure 7, with $u_{\tau }$ estimated from RSS in the outer layer. For comparison, we have also included data from previous studies (Banner Reference Banner1990; Sullivan et al. Reference Sullivan, McWilliams and Moeng2000; Peirson et al. Reference Peirson, Walker and Banner2014; Bopp Reference Bopp2018; Sullivan et al. Reference Sullivan, Banner, Morison and Peirson2018; Buckley et al. Reference Buckley, Veron and Yousefi2020; Funke et al. Reference Funke, Buckley, Schultze, Veron, Timmermans and Carpenter2021).

Figure 7(a) illustrates that form drag increases as the wave slope rises, which aligns with findings of Peirson et al. (Reference Peirson, Walker and Banner2014) and Funke et al. (Reference Funke, Buckley, Schultze, Veron, Timmermans and Carpenter2021). However, there is a slight discrepancy in the experimental results reported by Buckley et al. (Reference Buckley, Veron and Yousefi2020). Notably, Buckley et al. (Reference Buckley, Veron and Yousefi2020) calculated form drag by subtracting viscous stress from the total wind stress, while our study computes form drag directly by multiplying pressure by the wave slope. Funke et al. (Reference Funke, Buckley, Schultze, Veron, Timmermans and Carpenter2021) recently reconstructed the pressure field by solving a pressure Poisson equation in airflow over water waves. Their results, depicted as red plus signs in figure 7(a), differ from those of Buckley et al. (Reference Buckley, Veron and Yousefi2020) (shown as black plus signs) under similar wind-wave conditions. Nonetheless, the form drag for case S3 closely resembles that of Funke et al. (Reference Funke, Buckley, Schultze, Veron, Timmermans and Carpenter2021), further confirming the reliability of our numerical model. The various methods for calculating form drag lead to discrepancies in growth rate, as shown in figure 8(b), which we discuss later. Figure 7(b) indicates that viscous stress decreases as wave slope increases, remaining in a region very close to that observed in previous studies. Furthermore, it approaches the dashed line obtained by subtracting form drag from total wind stress. These results suggest that extrapolating RSS from the outer layer to the surface waves can be an effective approach for calculating total wind stress.

Figure 8.

(a) Wave growth rate $\beta$ versus wave slope and (b) wave growth rate $\gamma$ versus inverse wave age. The red dashed line represents the theoretical limit while the dash-dotted lines denote $\gamma ={2\pi (0.04\pm 0.02)u_{\tau }^{2}}/{{{c}^{2}}}$ summarised by Plant (Reference Plant1982). The black solid line is the theoretical result of Miles (Reference Miles1993).

Figure 8. Long description

Two scatter plots depict wave growth rates against wave slope and inverse wave age, comparing theoretical limits and empirical data. Panel A: A scatter plot shows the wave growth rate (beta) versus wave slope (ak). The x-axis is labeled ak, ranging from 0.05 to 0.40. The y-axis is labeled beta, ranging from 0 to 200. The red dashed line represents the theoretical limit, approximately following the equation (2/(ak)^2). The plot includes data points from various studies, including those by Pole (1967), Peirson & Garcia (2008), Yang & Shen (2010), Grare et al. (2013), Peirson et al. (2014), and the present study, each represented by different symbols. Panel B: Another scatter plot shows the wave growth rate (gamma) versus inverse wave age (u_tau/c). The x-axis is labeled u_tau/c, ranging from 0.1 to 10. The y-axis is labeled gamma, ranging from 10^-4 to 10^1. The red dashed line represents the theoretical limit, and the dash-dotted lines denote summaries by Plant (1982). The black solid line represents the theoretical result of Miles (1993). The plot includes data points from studies by Plant (1982), Miles (1993), Buckley et al. (2020), and the present study, each represented by different symbols.

We begin our analysis by examining the energy input, specifically the work rate resulting from wind stress, to evaluate wave growth. It is widely accepted that wave pressure distribution is the main factor influencing their growth. However, energy input can be decomposed into contributions from both pressure and viscous stress, each supporting wave growth in different proportions, as confirmed by Buckley et al. (Reference Buckley, Veron and Yousefi2020). We computed the total energy input, which is the sum of $\overline {\langle \tau _{p}^{+}\rangle {{c}^{+}}}$ and $\overline {\langle \tau _{\nu }^{+}\rangle {{{\tilde {u}}}^{+}}}$ (or $\overline {S}_{\textit{in}}^{+}$ ) associated with wind waves. Here, the wind forces are normalised by $\rho u_{\tau }^{2}$ , while the phase speed is normalised by the friction velocity, $u_{\tau }$ . Since the viscous scale is characterised as a length scale, the corresponding time scale can be determined using ${\nu }/{u_{\tau }^{2}}$ , which serves as the basis for normalisation. Consequently, the rate of wave energy normalised by the wave angular frequency $\omega =2\pi f$ and wave energy ${{E}_{w}}=( {1}/{2}){{\rho }_{w}}g{{a}^{2}}$ is expressed as

(2.11) \begin{align} \gamma =\frac {2\pi }{\omega {{E}_{w}}}\frac {{\rm d}{{E}_{w}}}{{\rm d}t}=\frac {2\pi }{\omega {{E}_{w}}}\left ({{{\overline {S}}}_{\textit{in}}}-{{D}_{d\nu }}\right )\!, \end{align}

where ${D}_{d\nu }$ is the energy dissipation resulting from wave viscosity and can be determined as ${{D}_{d\nu }}=2{{\nu }_{w}}{{\rho }_{w}}{{k}^{2}}g{{a}^{2}}$ (Lamb Reference Lamb1932; Phillips Reference Phillips1977; Buckley et al. Reference Buckley, Veron and Yousefi2020). Here, to accurately reflect the flow of wind waves, we ensure the ratio of air to water kinematic viscosity is approximately $\nu /\nu _w\approx 15$ . In our simulations, we set the viscosity of wind flow to ${\nu }=1.1\times {{10}^{-6}}\,{{{\text{m}}^{2}}}\,{\text{s}^{-1}}$ , allowing us to estimate $\nu _w$ as ${{\nu }_{w}}\approx 0.07\times {{10}^{-6}}\,{{{\text{m}}^{2}}}\,{\text{s}^{-1}}$ . We define the wave Reynolds number as ${{{Re}}_{w}}={c}/{k{{\nu }_{w}}}$ and use the dispersion relation of linear gravity waves. The deep-water assumption $kh\to \infty$ gives $\omega =kc=\sqrt {gk}$ , which is used to derive the second term in (2.11); specifically, this means ${2\pi {{D}_{d\nu }}}/{\omega {{E}_{w}}}=8\pi {Re}_{w}^{-1}$ . The energy input term can be represented as the product of wave surface stress ${\overline {\tau }}_{w}$ and phase velocity $c$ . Therefore, the first term of (2.11) is

(2.12) \begin{align} \frac {\text{2}\pi }{\omega {{E}_{w}}}{{\overline {S}}_{\textit{in}}}=\frac {\text{4}\pi {{{\overline {\tau }}}_{w}}c}{\omega {{\rho }_{w}}g{{a}^{2}}}=\frac {\text{4}\pi {{{\overline {\tau }}}_{w}}c}{kc{{\rho }_{w}}g{{a}^{2}}}=\frac {\text{4}\pi {{{\overline {\tau }}}_{w}}}{{{\rho }_{w}}{{a}^{2}}{{k}^{2}}{{c}^{2}}}=2\pi \frac {\rho }{{{\rho }_{w}}}\frac {2{{{\overline {\tau }}}_{w}}}{\rho u_{\tau }^{2}}\frac {u_{\tau }^{2}}{{{c}^{2}}}\frac {1}{{{(ak)}^{2}}}. \end{align}

We further rewrite the wave surface stress as $\beta ={2\overline {\tau }_{w}^{+}}/{{{(ak)}^{2}}}$ , where the ‘+’ sign denotes normalisation by $\rho u_{\tau }^{2}$ . This coefficient is associated with wave growth and incorporates both pressure- and viscous-related contributions to energy input. Therefore, the non-dimensional wave growth rate can be formulated as

(2.13) \begin{align} \gamma =\frac {\text{2}\pi \rho \beta u_{\tau }^{2}}{{{\rho }_{w}}{{c}^{2}}}-\frac {8\pi }{{Re}_{w}}. \end{align}

Figure 8 illustrates the variations of two wave growth parameters in relation to wave slope and inverse wave age. In figure 8(a), the parameter $\beta$ decreases as the wave slope increases, approaching the theoretical limit of $\beta ={2}/{{{(ak)}^{2}}}$ , which indicates that the wind stress is completely converted into shear stress. The current value of $\beta$ aligns closely with the experimental findings reported by Grare et al. (Reference Grare, Peirson, Branger, Walker, Giovanangeli and Makin2013), but it is higher than the results by Peirson & Garcia (Reference Peirson and Garcia2008) and Yang & Shen (Reference Yang and Shen2010). We caution that the $\beta$ value determined by Yang & Shen (Reference Yang and Shen2010) only accounted for pressure-induced form drag, neglecting the contribution from viscous stress.

Figure 8(b) plots the non-dimensional wave growth rate, primarily based on the study conducted by Plant (Reference Plant1982). This research specifically addresses turbulent flow over a prescribed travelling-wave boundary, focusing on single-phase flow, which differs from the two-phase air–water flow examined in earlier studies. Consequently, we determine the wave parameter by matching the actual ratio of air density to water density. The wave Reynolds number is found to be of the order of $\sim\!{{10}^{4}}$ , resulting in relatively low dissipation; however, this factor has been included in the calculation of the non-dimensional wave growth rate. The red dash-dotted lines in the graph, derived from Plant (Reference Plant1982) based on multiple experimental studies, are represented as $\gamma ={2\pi (0.04\pm 0.02)u_{\tau }^{2}}/{{{c}^{2}}}$ . The theoretical outcome from Miles (Reference Miles1993), which accounts for wave-induced turbulent effects, is depicted by the solid black line. Our results agree well with previous data, fall within the range of empirical evidence reported by Plant (Reference Plant1982) and are consistent with the theoretical prediction of Miles (Reference Miles1993). It is important to note that the rate reported by Buckley et al. (Reference Buckley, Veron and Yousefi2020) is slightly higher than our findings. We attribute this discrepancy to the different methods used to calculate form drag, as mentioned above. The pressure reconstruction employed by Funke et al. (Reference Funke, Buckley, Schultze, Veron, Timmermans and Carpenter2021) significantly reduces form drag, thereby decreasing the energy input and resulting in a lower growth rate. Recently, Carpenter et al. (Reference Carpenter, Buckley and Veron2022) demonstrated that the Miles critical layer mechanism leads to significant wave growth in young, wind-forced waves (with a wave age of approximately 6.3) by comparing the growth rate determined through pressure reconstruction with that obtained using alternative methods. We caution that our numerical model simplifies the growth of wind waves into a turbulent boundary layer over a prescribed, moving wavy boundary, consistent with the Miles theory. Consequently, our numerical results closely match theoretical expectations. Nonetheless, the consistency with experimental results highlights the potential applicability of the Miles theory for predicting wave growth under specific wave-age conditions.

2.5. The general form of wind stress: phase shift and magnitude comparison

The main forces affecting wave dynamics are wind pressure and tangential viscous stress; their specific forms are essential for predicting wave behaviour. We begin by examining the general form of wind pressure. In a recent study, Wu et al. (Reference Wu, Popinet and Deike2022) revisited the growth of wind waves using simulations of two-phase air–water interfacial flow. They described the surface pressure as ${{p}_{s}}=\sum \nolimits _{n=1}^{\infty }{{{{\hat {p}}}_{n}}\cos (nkx+2\pi n{{\varphi }_{p\eta }})}$ , where ${\hat {p}}_{n}$ is the pressure amplitude of mode $n$ and ${\varphi }_{pn}$ signifies the pressure phase shift. We start by plotting the pressure-induced form drag profiles, which are determined by $\langle \tau _{p}^{+}\rangle =\langle {{p}^{+}}{{\eta }_{x}}\rangle$ , along the air–water interface, as shown in figure 9(a). The results align closely with those found by Wu et al. (Reference Wu, Popinet and Deike2022), although there is a slight discrepancy. Figure 9(c) displays the approximate sinusoidal pattern of the pressure profiles, showing a similar trend but slightly higher values compared with those reported by Loft et al. (Reference Loft, Kuhl, Buckley, Carpenter, Hinze, Veron and Rung2023). These discrepancies primarily occur on the windward side of the waves, where strong positive pressure or form drag is present, and can be attributable to differences in coupling methods. The two-phase air–water interfacial flow method employed by Wu et al. (Reference Wu, Popinet and Deike2022) and Loft et al. (Reference Loft, Kuhl, Buckley, Carpenter, Hinze, Veron and Rung2023) is a fully coupled approach designed to replicate wind–wave interactions. In contrast, our current method utilises a one-way coupling. It is well established that the strong pressure on the windward side transfers energy into the waves, driving their evolution. Consequently, the pressure or form drag in our simplified numerical model is noticeably lower. As a result, the wave growth rate we observe is notably lower than the predictions made by the Miles theory and the results presented in figure 15 of Wu et al. (Reference Wu, Popinet and Deike2022). Nevertheless, our findings provide a reasonable basis for developing general forms of wind stress, as they yield a rational wave growth rate that compares favourably with various experimental results. As illustrated in figure 9(b), the mean pressure varies linearly with wave slope. Furthermore, the phase shift of $2\pi {{\varphi }_{p\eta }}\approx \text{120}{^{\circ }}$ in figure 9(d) is consistent with the findings of Kihara et al. (Reference Kihara, Hanazaki, Mizuya and Ueda2007), Druzhinin et al. (Reference Druzhinin, Troitskaya and Zilitinkevich2012) and Husain et al. (Reference Husain, Buckley, Yousefi, Veron and Sullivan2019), supporting the experimental range of $( {\pi }/{2},\pi )$ reported by Donelan et al. (Reference Donelan, Babanin, Young and Banner2006) and Grare (Reference Grare2009).

Figure 9.

(a) Pressure-induced form drag profiles along waves; (b) mean pressure versus wave slope; (c) pressure profiles along waves; (d) phase shift between pressure and wave profiles versus wave age. For comparison, the results of the two-phase air–water interfacial flow simulations by Wu et al. (Reference Wu, Popinet and Deike2022) ( $ak=0.15$ , $c/u_{\tau }=2$ ) and Loft et al. (Reference Loft, Kuhl, Buckley, Carpenter, Hinze, Veron and Rung2023) ( $ak=0.09$ , $c/u_{\tau }=6.4$ ) are also included. The phase shift results by Kihara et al. (Reference Kihara, Hanazaki, Mizuya and Ueda2007) (red plus), Druzhinin et al. (Reference Druzhinin, Troitskaya and Zilitinkevich2012) (red cross) and Husain et al. (Reference Husain, Buckley, Yousefi, Veron and Sullivan2019) (red solid circle) are also plotted.

Belcher & Hunt (Reference Belcher and Hunt1998) pointed out that wave-induced pressure leads to form drag, which is the dominant mechanism for momentum exchange in wind waves (Cao & Shen Reference Cao and Shen2021). The present study examines turbulent flow over a moving wavy boundary, characterised by one-way coupling, in contrast to the strong coupling typically observed in realistic wind waves. Consequently, the pressure distribution induced by the moving wave can be approximated as equivalent to the pressure at the interface, which can be utilised to assess variations in wave height caused by the wind. The wave-induced pressure can be determined by ${{\tilde {p}}^{+}}=\langle {{p}^{+}}\rangle -{{\overline {p}}^{+}}$ . Since the pressure shown in figure 9(c) varies approximately sinusoidally, it is reasonable to represent the wave-induced pressure in the form ${{\tilde {p}}^{+}}={{\hat {p}}^{+}}\cos (kx+\text{2}\pi {{\varphi }_{p\eta }})$ . In this analysis, we focus solely on the dominant pressure mode related to the linear wave wall boundary. In fact, even when considering the effects of wind on nonlinear waves, it remains valid to use the pressure from linear wave wall turbulence as a reference. Yang & Shen (Reference Yang and Shen2010) conducted numerical investigations of turbulent flow over various wavy boundaries and found that, in the case of slow waves, the Stokes wave scenario exhibits nearly identical values for both maximum and minimum surface pressures when compared with the Airy wave case. Based on the results illustrated in figure 9, the profiles of wave-induced pressure can be articulated in the form of ${{\tilde {p}}^{+}}\sim ak\cos (kx+2\pi {{\varphi }_{p\eta }})$ . Therefore, for a surface wave $\eta =ak\cos (kx)$ , we have ${{\tilde {p}}^{+}}\sim \eta (kx+2\pi {{\varphi }_{p\eta }})$ .

Next, we focus on the viscous stress along the waves, as shown in figure 10. We also include experimental data from Buckley et al. (Reference Buckley, Veron and Yousefi2020) and Yousefi et al. (Reference Yousefi, Veron and Buckley2020a ) for comparison. It is noted that our current wave age, ${c}/{{{u}_{\tau }}}=3.89$ , slightly differs from the value of ${c}/{{{u}_{\tau }}}=3.91$ reported by Buckley et al. (Reference Buckley, Veron and Yousefi2020). This discrepancy arises because the friction velocity was not predetermined; instead, we set different phase-to-bulk velocity ratios to approximate the wave age. Moreover, in the experimental study, profiles were obtained by averaging stress measurements taken between 284 and 664 μm above the air–water interface, corresponding to the region of ${{\zeta }^{+}}\approx 3.38{-}7.5$ . Figure 10(a) suggests that the viscous stress profile for case S1 displays a slight deviation from the experimental results; however, the overall wavy pattern is consistent with the findings of Buckley et al. (Reference Buckley, Veron and Yousefi2020). The peak value of viscous stress occurs on the windward side near the crests, but this peak position subtly shifts downstream with increasing wave slope, likely due to flow sheltering effects (Buckley et al. Reference Buckley, Veron and Yousefi2020). Figure 10(b) indicates that our numerical results agree well with the experimental results reported by Yousefi et al. (Reference Yousefi, Veron and Buckley2020a ), albeit with a minor phase shift.

Figure 10.

Profiles along the wave for (a) phase-averaged tangential viscous stress and (b) wave-induced tangential viscous stress. The circles denote the experimental results (EXP) of Buckley et al. (Reference Buckley, Veron and Yousefi2020). The triangles represent the experimental results (EXP) of Yousefi et al. (Reference Yousefi, Veron and Buckley2020a ).

Figure 10. Long description

The image contains two line graphs labeled Panel A and Panel B. Panel A shows the phase-averaged tangential viscous stress along the wave, while Panel B shows the wave-induced tangential viscous stress along the wave. Both graphs have the x-axis labeled as x/λx, ranging from 0 to 2. The y-axis in Panel A is labeled as <τθ>+, ranging from -0.5 to 1.5, and in Panel B is labeled as τθ’, ranging from -1 to 1. Each graph includes multiple lines representing different data series labeled as S1, S2, S3, S5, S8, S9, and S10. The circles and triangles denote experimental results from Buckley et al. (2020) and Yousefi et al. (2020a), respectively. Below each graph, there is a visual representation of the wave profile labeled as η. The graphs illustrate the variations in tangential viscous stress along the wave, with different data series showing distinct trends and patterns. The experimental data points are scattered along the lines, indicating the correlation between the theoretical models and experimental results.

Buckley et al. (Reference Buckley, Veron and Yousefi2020) employed a phase-averaged decomposition to effectively quantify the energy transfer from wind-tangential viscous stress into wave currents and wave motion. The ensemble-averaged viscous stress was found to contribute to wave currents. Therefore, it is reasonable to assume that $\overline {\tau }_{\nu }^{+}$ provides a mean shear to stimulate wave currents, while $\tilde {\tau }_{\nu }^{+}$ contributes to wave motion. Consequently, the results shown in figure 10(b) can serve as a valuable reference for understanding wave motion. The tangential viscous stress related to wave motion exhibits a very slight phase shift relative to wave elevation. As a consequence, the tangential viscous stress associated with wave motion can be expressed as $\tilde {\tau }_{\nu }^{+}\sim ak\cos (kx)$ . It is worth noting that wind pressure is typically about one order of magnitude larger than tangential viscous stress. Overall, both wind pressure and tangential viscous stress are related to surface displacement. These approximations will serve as references in the subsequent theoretical derivations.

3.1. Theoretical derivation of wind-induced waves

Traditionally, water waves are treated as an incompressible and inviscid fluid that flows irrotationally, leading to the application of potential theory. In this section, we build on the explicit forms of wind pressure and tangential viscous stress to develop an improved model for water waves, inspired by quasi-linear potential theory under the assumption of weak viscosity. Based on the Helmholtz decomposition of the velocity field, the essence of the quasi-potential theory is to retain complete information about the irrotational component while linearly approximating the solenoidal part. Although the previous section involved computations of turbulent wind in three dimensions, we focus on a two-dimensional water-wave problem and reserve the three-dimensional wind-wave problem for future studies. To avoid any confusion, we will redefine all variables used in the subsequent derivations. First, the kinematic boundary condition can be stated as

(3.1) \begin{align} \frac {\partial \eta }{\partial t}+u(x,\eta ,t)\frac {\partial \eta }{\partial x}=w(x,\eta ,t),\quad \text{at}\quad z=\eta (x,t), \end{align}

where $\eta (x,t)$ represents the free-surface displacement and $(u,w)$ is the velocity field of the water. Moreover, the dynamic boundary condition requires a force balance at the wind–wave interface, expressed as

(3.2) \begin{align} -({{p}_{w}}-p)\boldsymbol{n}+\tau \boldsymbol{\cdot }\boldsymbol{n}={{\tau }_{s}}\boldsymbol{\cdot }\boldsymbol{n}, \end{align}

where $p_w$ is the water pressure, $p$ the wind pressure, $\tau$ the water viscous stress, $\tau _s$ the wind viscous stress and $\boldsymbol{n}$ the unit normal vector pointing upward at the free surface. Following Dias et al. (Reference Dias, Dyachenko and Zakharov2008), we describe the continuity of tangential and normal stresses at $z=\eta (x,t)$ as

(3.3) \begin{align} \left \{\begin{aligned} &{{\rho }_{w}}{{\nu }_{w}}\left (\frac {\partial u}{\partial z}+\frac {\partial w}{\partial x}\right )+{{\rho }_{w}}{{\nu }_{w}}{{\mathcal{N}}_{t}}={{\tau }_{a}},\\ &{{p}_{w}}-p=2{{\rho }_{w}}{{\nu }_{w}}\frac {\partial w}{\partial z}+{{\rho }_{w}}{{\nu }_{w}}{{\mathcal{N}}_{n}}, \end{aligned}\right . \end{align}

where ${\rho }_{w}$ and ${\nu }_{w}$ are the density and kinematic viscosity of water, respectively, ${\mathcal{N}}_{t}$ and ${\mathcal{N}}_{n}$ denote the nonlinear terms and $\tau _a$ refers specifically to the wind tangential viscous stress.

The general form of tangential viscous stress is ${{\tau }_{a}}={{\tau }_{a0}}+{{\hat {\tau }}_{a}}{{\text{e}}^{{\rm i}(kx-\omega t)}}$ , where ${\tau }_{a0}$ can be determined through ensemble averaging along the wave. This formulation implies that the tangential viscous stress can be decomposed into two components: a mean shear component, ${\tau }_{a0}$ , and a wave-related contribution, ${{\hat {\tau }}_{a}}{{\text{e}}^{{\rm i}(kx-\omega t)}}$ , as discussed in § 2.5. Since this study focuses on the effects of wind on wave motion, we will not consider the mean shear stress related to wave currents. Additionally, we can further decompose the velocity into irrotational and solenoidal components. Introducing the potential function $\phi$ and the stream function $\psi$ , the Helmholtz decomposition yields

(3.4) \begin{align} u=\frac {\partial \phi }{\partial x}-\frac {\partial {{\psi }}}{\partial z},\quad w=\frac {\partial \phi }{\partial z}+\frac {\partial {{\psi }}}{\partial x}. \end{align}

Substituting the expressions (3.4) into the linearised two-dimensional Navier–Stokes equations (the linearisation aims at expressing ${\partial \psi }/{\partial x}$ in terms of $\eta$ , as $ {\partial \eta }/{\partial t}={\partial \phi }/{\partial z}+{\partial \psi }/{\partial x}$ at $z=0$ ), yields the following equations for the potential and stream functions:

(3.5) \begin{align} \frac {{{\partial }^{\text{2}}}\phi }{\partial {{x}^{2}}}+\frac {{{\partial }^{\text{2}}}\phi }{\partial {{z}^{2}}}=0,\quad \frac {\partial \phi }{\partial t}=-\frac {{{p}_{w}}}{{{\rho }_{w}}}-gz+\frac {p}{{{\rho }_{w}}}, \\[-28pt] \nonumber \end{align}

(3.6) \begin{align} \frac {\partial {{\psi }}}{\partial t}={{\nu }_{w}}\left (\frac {{{\partial }^{2}}{{\psi }}}{\partial {{x}^{2}}}+\frac {{{\partial }^{2}}{{\psi }}}{\partial {{z}^{2}}}\right ), \\[0pt] \nonumber \end{align}

where $g$ is the acceleration due to gravity. Following Lamb (Reference Lamb1932) and Dias et al. (Reference Dias, Dyachenko and Zakharov2008), the solutions for the potential and stream functions can be expressed as

(3.7) \begin{align} \phi ={{\phi }_{0}}{{\text{e}}^{{\rm i}(kx-\omega t)}}{{\text{e}}^{\left | k \right |z}}\quad \text{and}\quad {{\psi }}={{\psi }_{0}}{{\text{e}}^{{\rm i}(kx-\omega t)}}{{\text{e}}^{mz}}, \end{align}

where ${{m}^{2}}={{k}^{2}}-\text{i}( {\omega }/{{{\nu }_{w}}})$ . The linearisation of the kinematic boundary condition reads $ {\partial \eta }/{\partial t}=w$ . The velocity is evaluated at $z=0$ in the linear approximation, and the wave profile can be written as

(3.8) \begin{align} \eta (x,t)=\hat {\eta }{{\text{e}}^{{\rm i}(kx-\omega t)}}=\frac {1}{\omega }(\text{i}\left | k\right |{{\phi }_{0}}-k{{\psi }_{0}}){{\text{e}}^{{\rm i}(kx-\omega t)}}. \end{align}

Linearising the dynamic boundary conditions yields

(3.9) \begin{align} \left \{ \begin{aligned} &{{\rho }_{w}}{{\nu }_{w}}\left (\frac {\partial u}{\partial z}+\frac {\partial w}{\partial x}\right )={{{\hat {\tau }}}_{a}}{{\text{e}}^{{\rm i}(kx-\omega t)}},\\[6pt]&{{p}_{w}}-p=2{{\rho }_{w}}{{\nu }_{w}}\frac {\partial w}{\partial z}. \end{aligned} \right . \end{align}

Since we ignore the effect of wave currents (or mean shear), the tangential viscous stress related to wave motion is assumed to be

(3.10) \begin{align} {{\hat {\tau }}_{a}}{{\text{e}}^{{\rm i}(kx-\omega t)}}=\varGamma \hat {\eta }\mathrm{e}^{{\rm i}(kx-\omega t)}. \end{align}

After some algebra, the relation between ${\phi }_{0}$ and ${\psi }_{0}$ and the dispersion relation $\omega (k)$ can be respectively written as

(3.11) \begin{align} {{\psi }_{0}}=\frac {\frac {\varGamma \left |k\right |}{{{\rho }_{w}}{{\omega }^{2}}}-\frac {2{{\nu }_{w}}\left |k\right |k}{\omega }}{1+\frac {2\text{i}{{\nu }_{w}}{{k}^{2}}}{\omega }-\frac {\text{i}\varGamma k}{{{\rho }_{w}}{{\omega }^{2}}}}{{\phi }_{0}} \end{align}

and

(3.12) \begin{align} \begin{aligned} &{{\left (2-\frac {\text{i}\omega }{{{\nu }_{w}}{{k}^{2}}}\right )}^{2}}+\frac {g}{\nu _{w}^{2}{{\left |k\right |}^{3}}}-4\sqrt {1-\frac {\text{i}\omega }{{{\nu }_{w}}{{k}^{2}}}}\\ &=\left (\frac {2\omega }{{{\nu }_{w}}{{k}^{2}}}\sqrt {\frac {\text{i}\omega }{{{\nu }_{w}}{{k}^{2}}}-1}-\frac {g}{{{\nu }_{w}}{{k}^{3}}}\right )\left (\frac {\frac {2\varGamma {{\nu }_{w}}{{k}^{2}}\left |k\right |}{{{\rho }_{w}}{{\omega }^{3}}}-\frac {2\varGamma {{\nu }_{w}}{{k}^{2}}\left |k\right |}{{{\rho }_{w}}{{\omega }^{2}}}-\frac {\text{i}\varGamma \left |k\right |}{{{\rho }_{w}}\omega }}{\frac {\text{i}\varGamma k}{{{\rho }_{w}}\omega }-\frac {1}{\omega }-\frac {2\text{i}{{\nu }_{w}}{{k}^{2}}}{\omega }}\right ). \end{aligned} \end{align}

According to Lamb (Reference Lamb1932), by introducing the dimensionless parameter ${{\nu }^{+}}= {{{\nu }_{w}}{{k}^{2}}}/{\left |\omega \right |}$ and considering small viscosity, namely ${{\nu }_{w}}{{k}^{2}}\ll \left |\omega \right |$ , equation (3.12) can be expanded in powers of ${\nu }^{+}$ . Ignoring higher-order terms of ${\nu }^{+}$ yields the following relation:

(3.13) \begin{align} \begin{aligned} &{{\rho }_{w}}\left |\omega \right |\big (1+2\text{i}{{\nu }_{w}}{{k}^{2}}\big )\big [{{\left |\omega \right |}^{2}}+4\text{i}{{\nu }_{w}}{{k}^{2}}\left |\omega \right |-g\left |k\right |\big ]\\[4pt]=\,&\varGamma \left |\omega \right |\big [\text{i}k{{\left |\omega \right |}^{2}}+2{{\nu }_{w}}k\left (\left |k\right |-2k\right )\left |\omega \right |+\text{i}gk\left |k\right |\left ({{\nu }_{w}}-1\right )\big ]+2\varGamma g\nu _{w}^{2}{{k}^{3}}\left |k\right |. \end{aligned} \end{align}

It is recalled that the tangential viscous stress is smaller than the pressure stress. If we ignore the terms related to $\varGamma$ , then an approximation of $\omega (k)$ can be obtained, $\omega =\pm \sqrt {g|k|}-2\text{i}{{\nu }_{w}}{{k}^{2}}$ , which aligns with the finding of Dias et al. (Reference Dias, Dyachenko and Zakharov2008). Following the derivation of Dias et al. (Reference Dias, Dyachenko and Zakharov2008), we can describe the flow using only the potential part of the velocity. First, $ {\partial {{\psi }}}/{\partial x}$ can be written as

(3.14) \begin{align} \frac {\partial {{\psi }}}{\partial x}\bigg |_{z=0} =\text{i}k{{\psi }_{0}}{{\text{e}}^{{\rm i}(kx-\omega t)}}=\left (\frac {\varGamma k}{{{\rho }_{w}}\omega }-2{{\nu }_{w}}{{k}^{2}}\right )\eta \end{align}

or

(3.15) \begin{align} \frac {\partial {{\psi }}}{\partial x}\bigg |_{z=0} =2{{\nu }_{w}}\frac {{{\partial }^{2}}\eta }{\partial {{x}^{2}}}-\frac {\varGamma }{{{\rho }_{w}}}\partial _{t}^{-1}{{\partial }_{x}}\eta =2{{\nu }_{w}}\frac {{{\partial }^{2}}\eta }{\partial {{x}^{2}}}-\int _{0}^{t}{\frac {\varGamma }{{{\rho }_{w}}}\frac {\partial \eta }{\partial x}\,{\rm d}t}, \end{align}

where the pseudo-differential operator $\partial _{t}^{-1}$ is defined through its symbol $\mathrm{i}/\omega$ . The kinematic boundary condition is thus determined only by $\phi$ and $\eta$ , which is expressed as

(3.16) \begin{align} \frac {\partial \eta }{\partial t}=\frac {\partial \phi }{\partial z}+2{{\nu }_{w}}\frac {{{\partial }^{2}}\eta }{\partial {{x}^{2}}}-\int _{0}^{t}{\frac {\varGamma }{{{\rho }_{w}}}\frac {\partial \eta }{\partial x}\,{\rm d}t}\quad \text{at}\;z=0. \end{align}

The linearised Bernoulli equation with wind pressure and tangential viscous stress then reads

(3.17) \begin{align} \frac {\partial \phi }{\partial t}+g\eta +\frac {p}{{{\rho }_{w}}}=-2{{\nu }_{w}}\frac {{{\partial }^{2}}\phi }{\partial {{z}^{2}}}\quad \text{at}\;z=0. \end{align}

So, these are the linearised kinematic and dynamic boundary conditions for water waves driven by both wind pressure and tangential viscous stress. Next, we generalise these equations to the nonlinear situation. According to the results of Ruvinsky et al. (Reference Ruvinsky, Feldstein and Freidman1991) and Dias et al. (Reference Dias, Dyachenko and Zakharov2008), in the limit of low viscosity and small tangential viscous stress, the nonlinear Bernoulli equation reads

(3.18) \begin{align} \frac {\partial \phi }{\partial t}+\frac {1}{2}\left [{{\left (\frac {\partial \phi }{\partial x}\right )}^{2}}+{{\left (\frac {\partial \phi }{\partial z}\right )}^{2}}\right ]+g\eta +\frac {p}{{{\rho }_{w}}}=-2{{\nu }_{w}}\frac {{{\partial }^{2}}\phi }{\partial {{z}^{2}}}\quad \text{at}\;z=\eta , \end{align}

and the nonlinear kinematic boundary condition can be expressed as

(3.19) \begin{align} \frac {\partial \eta }{\partial t}+\frac {\partial \eta }{\partial x}\frac {\partial \phi }{\partial x}=\frac {\partial \phi }{\partial z}+2{{\nu }_{w}}\frac {{{\partial }^{2}}\eta }{\partial {{x}^{2}}}-\int _{0}^{t}{\frac {\varGamma }{{{\rho }_{w}}}\frac {\partial \eta }{\partial x}\,{\rm d}t}\quad \text{at}\;z=\eta . \end{align}

Thus, the governing equations for nonlinear waves driven by both wind pressure and tangential viscous stress (the component leading to wave-related motion) propagating on water of depth $h$ have been obtained. In summary, equations (3.5), (3.18) and (3.19), together with the impermeability boundary condition ( $\phi _z=0$ ) at the bottom, form a closed system for the unknowns $\phi$ and $\eta$ . This system is referred to as the quasi-potential theory.

It must be noted that, from a formal perspective, our proposed model closely resembles that developed by Maleewong & Grimshaw (Reference Maleewong and Grimshaw2024), who explored the evolution of wave groups in water of finite depth, accounting for various factors such as wind pressure, tangential viscous stress, bottom stress and water viscosity. However, we have modified our model to better align with physical reality, particularly concerning turbulent winds. For instance, Maleewong & Grimshaw (Reference Maleewong and Grimshaw2024) used a simplified version of the Miles model by setting $p=\rho {W_{r}^2}\overline {\beta }\eta _x$ (which involves neglecting $\overline {\alpha }$ in (3.26)), where $W_r$ is a reference velocity and $\overline {\beta }$ is an empirical, dimensionless constant. This treatment indicates that the wind pressure aligns with the wave slope or has a phase shift of $\pi /2$ with the wave elevation. However, in our simulations of the verified numerical model from the previous section, we found that the phase shift of wind pressure exceeds $\pi /2$ , a result consistent with earlier studies (Kihara et al. Reference Kihara, Hanazaki, Mizuya and Ueda2007; Druzhinin et al. Reference Druzhinin, Troitskaya and Zilitinkevich2012; Husain et al. Reference Husain, Buckley, Yousefi, Veron and Sullivan2019). Experimental research conducted by Donelan et al. (Reference Donelan, Babanin, Young and Banner2006) and Grare (Reference Grare2009) also reported a phase shift range of ${(\pi }/{2},\pi )$ for wind pressure. Additionally, a recent two-phase air–water numerical model developed by Wu et al. (Reference Wu, Popinet and Deike2022) indicated a phase shift greater than $\pi /2$ for $ak=0.15$ and $c/u_{\tau }=2$ , which characterises a strong wind force. To account for these findings in the modelling, we employed a generalised formulation of wind pressure.

3.2. Asymptotic expansion

In this part, we use asymptotic expansion to derive the wave profile influenced by wind pressure and tangential viscous stress. To facilitate comparison with Zdyrski & Feddersen (Reference Zdyrski and Feddersen2020), we make the deep-water assumption by sending the water depth, $h$ , to infinity. It is noted that this asymptotic procedure can also be applied to scenarios with finite depth. Following Mei et al. (Reference Mei, Stiassnie and Yue2005) and Zdyrski & Feddersen (Reference Zdyrski and Feddersen2020), the variables are normalised as follows:

(3.20) \begin{gather} {{x}^{+}}=kx,\quad {{z}^{+}}=kz,\quad {{\eta }^{+}}=k\eta ,\quad {{t}^{+}}=t\sqrt {gk},\quad {{\phi }^{+}}=\sqrt {\frac {{{k}^{3}}}{g}}\phi ,\nonumber \\ {{p}^{+}}=\frac {pk}{{{\rho }_{w}}g},\quad {{\varGamma }^{+}}=\frac {\varGamma }{{{\rho }_{w}}g},\quad \nu _{w}^{+}=\frac {{{\nu }_{w}}{{k}^{2}}}{\sqrt {gk}}. \end{gather}

Suppressing the superscript ‘+’ in all variables for ease of notation, the dimensionless system reads

(3.21) \begin{align} \begin{aligned} \frac {{{\partial }^{2}}\phi }{\partial {{x}^{2}}}+\frac {{{\partial }^{2}}\phi }{\partial {{z}^{2}}}=0\quad &\text{for}\quad z\lt \eta ,\\ \frac {\partial \phi }{\partial z}\rightarrow 0\quad &\text{at}\quad z\rightarrow -\infty ,\\ \frac {\partial \eta }{\partial t}+\frac {\partial \eta }{\partial x}\frac {\partial \phi }{\partial x}=\frac {\partial \phi }{\partial z}+2{{\nu }_{w}}\frac {{{\partial }^{2}}\eta }{\partial {{x}^{2}}}-\int _{0}^{t}{\varGamma \frac {\partial \eta }{\partial x}\,{\rm d}t}\quad &\text{at}\quad z=\eta ,\\ \frac {\partial \phi }{\partial t}+\frac {1}{2}\left [{{\left (\frac {\partial \phi }{\partial x}\right )}^{2}}+{{\left (\frac {\partial \phi }{\partial z}\right )}^{2}}\right ]+\eta +p=-2{{\nu }_{w}}\frac {{{\partial }^{2}}\phi }{\partial {{z}^{2}}}\quad &\text{at}\quad z=\eta . \end{aligned} \end{align}

Motivated by Zdyrski & Feddersen (Reference Zdyrski and Feddersen2020), the method of multiple scales, which is built around the recognition that the solution changes in different ways over different time scales, will be used in the subsequent analyses. The idea of the method is to introduce new independent variables, $t_n=\varepsilon ^nt$ ( $n=0,1,2,\dots$ and $\varepsilon \ll 1$ ), one to describe each time scale, and to regard the solution as a function of this extended set of independent variables:

(3.22) \begin{align} \begin{aligned} &\eta =\sum \limits _{n=1}^{\infty }{{{\varepsilon }^{n}}{{\eta }_{n}}(x,{{t}_{0}},{{t}_{1}},\dots )},\quad \phi =\sum \limits _{n=1}^{\infty }{{{\varepsilon }^{n}}{{\phi }_{n}}(x,z,{{t}_{0}},{{t}_{1}},\dots )},\\ &{{p}}=\sum \limits _{n=1}^{\infty }{{{\varepsilon }^{n}}{{p}_{n}}(x,{{t}_{0}},{{t}_{1}},\dots )},\quad \varGamma =\sum \limits _{n=1}^{\infty }{{{\varepsilon }^{n}}{{\varGamma }_{n}}(x,{{t}_{0}},{{t}_{1}},\dots )}. \end{aligned} \end{align}

Here, to emphasise the effect of tangential viscous stress, we assume it is relatively large. This assumption can be established because we found the ratio of tangential viscous stress to pressure approximates $O(\varepsilon )\sim O(1)$ . We substitute the above asymptotic series into the governing system (3.21) and expand the kinematic and dynamic boundary conditions at $z=0$ . Similarly, we first approximate the magnitude of viscosity as $\sim\!O(\varepsilon )$ , namely ${{\nu }_{w}}\approx \varepsilon {{\nu }_{1}}$ , where ${{\nu }_{\text{1}}}\sim O(1)$ , to emphasise its effect. By collecting terms by order, we derive the following equations at $z=0$ :

(3.23) \begin{align} o(\varepsilon ){:}\left \{\begin{aligned} &\frac {\partial {{\eta }_{1}}}{\partial {{t}_{0}}}-\frac {\partial {{\phi }_{1}}}{\partial z}=0,\\[5pt] &{{\eta }_{1}}+\frac {\partial {{\phi }_{1}}}{\partial {{t}_{0}}}+{{p}_{1}}=0, \end{aligned}\right . \end{align}

and

(3.24) \begin{align} o({{\varepsilon }^{2}}){:} \left \{\begin{aligned} &\frac {\partial {{\eta }_{1}}}{\partial {{t}_{1}}}-\frac {\partial {{\phi }_{2}}}{\partial z}+\frac {\partial {{\eta }_{2}}}{\partial {{t}_{0}}}+\frac {\partial {{\eta }_{1}}}{\partial x}\frac {\partial {{\phi }_{1}}}{\partial x}={{\eta }_{1}}\frac {{{\partial }^{2}}{{\phi }_{1}}}{\partial {{z}^{2}}}+2{{\nu }_{1}}\frac {{{\partial }^{2}}{{\eta }_{1}}}{\partial {{x}^{2}}}-\int _{0}^{{{t}_{0}}}{{{\varGamma }_{1}}\frac {\partial {{\eta }_{1}}}{\partial x}{\rm d}{{t}_{0}}},\\[5pt]&{{\eta }_{2}}+\frac {\partial {{\phi }_{1}}}{\partial {{t}_{1}}}+{{p}_{2}}+2{{\nu }_{1}}\frac {{{\partial }^{2}}{{\phi }_{1}}}{\partial {{z}^{2}}}+\frac {\partial {{\phi }_{2}}}{\partial {{t}_{0}}}+{{\eta }_{1}}\frac {{{\partial }^{2}}{{\phi }_{1}}}{\partial z\partial {{t}_{0}}}+\frac {1}{2}\left [{{\left (\frac {\partial {{\phi }_{1}}}{\partial x}\right )}^{2}}+{{\left (\frac {\partial {{\phi }_{1}}}{\partial z}\right )}^{2}}\right ]=0. \end{aligned}\right . \end{align}

Further expanding the variables into the Fourier series yields the component equations; here, the Fourier series are

(3.25) \begin{align} \begin{aligned} &{{\phi }_{n}}={{{\hat {\phi }}}_{n,0}}({{t}_{0}},{{t}_{1}},\dots )+\sum \limits _{m=1}^{n}{{{{\hat {\phi }}}_{n,m}}({{t}_{0}},{{t}_{1}},\dots ){{\text{e}}^{{\rm i}mx+mz}}},\\ &{{\eta }_{n}}=\sum \limits _{m=0}^{n}{{{{\hat {\eta }}}_{n,m}}({{t}_{0}},{{t}_{1}},\dots ){{\text{e}}^{{\rm i}mx}}},\quad {{p}_{n}}=\sum \limits _{m=0}^{n}{{{{\hat {p}}}_{n,m}}({{t}_{0}},{{t}_{1}},\dots ){{\text{e}}^{{\rm i}mx}}}. \end{aligned} \end{align}

We now turn our attention to ${\hat {p}}_{n,m}$ . According to the Miles theory, the pressure profile along the waves can be described as

(3.26) \begin{align} {{p}_{M}}=\left (\overline {\alpha }+\text{i}\overline {\beta }\right )\rho U_{\!f}^{2}k{{\eta }_{a}}, \end{align}

where $\overline {\alpha }$ and $\overline {\beta }$ are two empirical, dimensionless constants, $\eta _a$ is the analytic representation of $\eta$ and $U_{\!f}$ refers to the wind velocity at a reference height. Generally, pressure can be expressed as $p=r\rho {U_{\!f}^{2}}k\eta (kx+2\pi {{\varphi }_{p\eta }},t)$ , as discussed in § 2.5. This expression indicates a phase shift with wave profile, where $r$ is a dimensionless constant and $2\pi {{\varphi }_{p\eta }}$ is the phase shift between pressure and wave profile. Based on this assumption, we have ${{\hat {p}}_{n,m}}={{\hat {P}}_{m}}{{\hat {\eta }}_{n,m}}$ , where ${{\hat {P}}_{m}}=kP{{\text{e}}^{2\pi{\rm i}m{{\varphi }_{p\eta }}}}$ , and $P$ , characterised as the modulus of pressure, is assumed to be $O(1)$ . For the $O(\varepsilon )$ equations, we can directly obtain the results:

(3.27) \begin{align} {{\phi }_{1}}=-\text{i}{{\omega }_{0}}{{A}_{1}}({{t}_{1}}){{\text{e}}^{{\rm i}(x-{{\omega }_{0}}{{t}_{0}})+z}}\quad \text{and}\quad {{\eta }_{1}}={{A}_{1}}({{t}_{1}}){{\text{e}}^{{\rm i}(x-{{\omega }_{0}}{{t}_{0}})}}, \end{align}

where ${{\omega }_{0}}=\pm \sqrt {1+{{{\hat {P}}}_{1}}}$ , the plus sign denoting that the waves are propagating to the right. This result is consistent with the findings of Zdyrski & Feddersen (Reference Zdyrski and Feddersen2020), who examined the effects of wind pressure on the surface of an irrotational flow. According to their study, the wave profile can be expressed as

(3.28) \begin{align} k\eta ={{a}_{1}}k{{\text{e}}^{{\rm i}(kx-{{\omega }_{0}}{{t}_{0}})}}+{{({{a}_{1}}k)}^{2}}\frac {{{a}_{2}}}{a_{1}^{2}k}{{\text{e}}^{{\rm i}\left [ 2(kx-{{\omega }_{0}}{{t}_{0}})+\varTheta \right ]}}, \end{align}

where ${{a}_{1}}=\varepsilon {{A}_{1}}$ . The wave shape parameters, referred to as the harmonic phase $\varTheta$ and amplitude ${{{a}_{2}}}/{a_{1}^{2}k}$ , are specified in their study. The inclusion of tangential viscous stress indicates its contribution to the equations of $O({{\varepsilon }^{2}})$ . When comparing the magnitudes of pressure and tangential viscous stress, it is found that only $\varGamma _1$ is significant in the present system when expanded to $O({{\varepsilon }^{2}})$ . For the $O({{\varepsilon }^{2}})$ equations, ${{\hat {\eta }}_{2,0}}=0$ leads to the dominance of the first and second Fourier components, which can be expressed, respectively, as

(3.29) \begin{align} m=1{:} \left \{\begin{aligned} &{{{\hat {\phi }}}_{2,1}}-\frac {\partial {{{\hat {\eta }}}_{2,1}}}{\partial {{t}_{0}}}=\left (\frac {\partial {{A}_{1}}}{\partial {{t}_{1}}}-{{\varGamma }_{1}}\frac {{{A}_{1}}}{{{\omega }_{0}}}+2{{\nu }_{1}}{{A}_{1}}\right ){{\text{e}}^{-{\rm i}{{\omega }_{0}}{{t}_{0}}}},\\[5pt]&\frac {\partial {{{\hat {\phi }}}_{2,1}}}{\partial {{t}_{0}}}+\big (1+{{{\hat {P}}}_{1}}\big ){{{\hat {\eta }}}_{2,1}}=\left (\text{i}{{\omega }_{0}}\frac {\partial {{A}_{1}}}{\partial {{t}_{1}}}+2\text{i}{{\omega }_{0}}{{\nu }_{1}}{{A}_{1}}\right ){{\text{e}}^{-{\rm i}{{\omega }_{0}}{{t}_{0}}}} \end{aligned}\right . \end{align}

and

(3.30) \begin{align} m=2{:} \left \{\begin{aligned} &2{{{\hat {\phi }}}_{2,2}}-\frac {\partial {{{\hat {\eta }}}_{2,2}}}{\partial {{t}_{0}}}=2\text{i}{{\omega }_{0}}A_{1}^{2}({{t}_{1}}){{\text{e}}^{-2{\rm i}{{\omega }_{0}}{{t}_{0}}}},\\[5pt]&\frac {\partial {{{\hat {\phi }}}_{2,2}}}{\partial {{t}_{0}}}+\big (1+{{{\hat {P}}}_{2}}\big ){{{\hat {\eta }}}_{2,2}}=\omega _{0}^{2}A_{1}^{2}({{t}_{1}}){{\text{e}}^{-2{\rm i}{{\omega }_{0}}{{t}_{0}}}}. \end{aligned}\right . \end{align}

Eliminating ${\hat {\eta }}_{2,1}$ and ${\hat {\eta }}_{2,2}$ gives

(3.31) \begin{align} \frac {{{\partial }^{2}}{{{\hat {\phi }}}_{2,1}}}{\partial t_{0}^{2}}+\big (1+{{\hat {P}}_{1}}\big ){{\hat {\phi }}_{2,1}}=\omega _{\text{0}}^{\text{2}}\left (2\frac {\partial {{A}_{1}}}{\partial {{t}_{1}}}-{{\varGamma }_{1}}\frac {{{A}_{1}}}{{{\omega }_{0}}}+4{{\nu }_{1}}{{A}_{1}}\right ){{\text{e}}^{-{\rm i}{{\omega }_{0}}{{t}_{0}}}} \end{align}

and

(3.32) \begin{align} \frac {{{\partial }^{2}}{{{\hat {\phi }}}_{2,2}}}{\partial t_{0}^{2}}+2\big (1+{{\hat {P}}_{2}}\big ){{\hat {\phi }}_{2,2}}=\left [2(1+{{{\hat {P}}}_{2}})-2\omega _{0}^{2}\right ]\text{i}{{\omega }_{0}}A_{1}^{2}({{t}_{1}}){{\text{e}}^{-2{\rm i}{{\omega }_{0}}{{t}_{0}}}}. \end{align}

Equation (3.31) requires $2({\partial {{A}_{1}}}/{\partial {{t}_{1}}})-{{\varGamma }_{1}}( {{{A}_{1}}}/{{{\omega }_{0}}})+4{{\nu }_{1}}{{A}_{1}}=0$ to avoid the secular term. Therefore,

(3.33) \begin{align} {{A}_{1}}({{t}_{1}})=-\frac {1}{2}\mathcal{A}\exp \left [ \left (\frac {{{\varGamma }_{1}}}{2{{\omega }_{0}}}-2{{\nu }_{1}}\right ){{t}_{1}}\right ]\!, \end{align}

where $\mathcal{A}$ is a parameter related to wave amplitude. Assuming the initial wave profile resembles a classical Stokes wave, it follows that ${{A}_{1}}(0)=-({1}/{2})\mathcal{A}= {{{a}_{1}}}/{\varepsilon }$ ; therefore, ${{A}_{1}}({{t}_{1}})=({{{a}_{1}}}/{\varepsilon })\exp [ (({{{\varGamma }_{1}}}/{2{{\omega }_{0}}})-2{{\nu }_{1}} ){{t}_{1}} ]$ . On the other hand, ${\hat {\phi }}_{2,2}$ in (3.32) can be expressed as

(3.34) \begin{align} {{\hat {\phi }}_{2,2}}=\text{i}{{\omega }_{0}}\left (\frac {{{{\hat {P}}}_{1}}-{{{\hat {P}}}_{2}}}{1+2{{{\hat {P}}}_{1}}-{{{\hat {P}}}_{2}}}\right )A_{1}^{2}{{\text{e}}^{-2{\rm i}{{\omega }_{0}}{{t}_{0}}}}. \end{align}

Substituting (3.34) into (3.30) yields

(3.35) \begin{align} {{\hat {\eta }}_{2,2}}=\left ( \frac {1+{{{\hat {P}}}_{1}}}{1+2{{{\hat {P}}}_{1}}-{{{\hat {P}}}_{2}}} \right )A_{1}^{2}{{\text{e}}^{-2{\rm i}{{\omega }_{0}}{{t}_{0}}}}. \end{align}

It then follows that the wave profile can be expressed as

(3.36) \begin{align} \eta ={{a}_{1}}\exp \left [\left (\frac {{{\varGamma }_{1}}}{2{{\omega }_{0}}}-2{{\nu }_{1}}\right ){{t}_{1}}\right ]{{\text{e}}^{{\rm i}(x-{{\omega }_{0}}{{t}_{0}})}}+{{{a}_{2}}\exp \left [2\left (\frac {{{\varGamma }_{1}}}{2{{\omega }_{0}}}-2{{\nu }_{1}}\right ){{t}_{1}}\right ]}{{\text{e}}^{{\rm i}\left [2(x-{{\omega }_{0}}{{t}_{0}})\right ]}}, \end{align}

where $a_2=({(1+{{{\hat {P}}}_{1}})a_1^2})/({1+2{{{\hat {P}}}_{1}}-{{{\hat {P}}}_{2}}})$ . It is noted that (3.36) differs from (3.28) as it incorporates tangential viscous stress and water viscosity. We next focus on the wind effect on the wave profile based on (3.36). The real part of $\varGamma _1/(2\omega _0)-2\nu _1$ characterises the amplitude amplification factor within $t_1$ scale, while its imaginary part can be incorporated into the phase shift:

(3.37) \begin{align} {{\text{e}}^{{\rm i}\left \{(x-{{\omega }_{0}}{{t}_{0}})+\operatorname {Im}\left [\left (\frac {{{\varGamma }_{1}}}{2{{\omega }_{0}}}-2{{\nu }_{1}} \right ){{t}_{1}}\right ]\right \}}}\quad \text{and}\quad {{\text{e}}^{{2{\rm i}}\left \{(x-{{\omega }_{0}}{{t}_{0}})+\operatorname {Im}\left [\left (\frac {{{\varGamma }_{1}}}{2{{\omega }_{0}}}-2{{\nu }_{1}} \right ){{t}_{1}}\right ]\right \}}}\,; \end{align}

thus, the wave profile can be recast to

(3.38) \begin{align} \begin{aligned} \eta =&\,{{a}_{1}}\exp \left \{\operatorname {Re}\left [\left (\frac {{{\varGamma }_{1}}}{2{{\omega }_{0}}}-2{{\nu }_{1}}\right ){{t}_{1}}\right ]\right \}{{\text{e}}^{{\rm i}\left \{x-{{\omega }_{0}}{{t}_{0}}+\operatorname {Im}\left [\left (\frac {{{\varGamma }_{1}}}{2{{\omega }_{0}}}-2{{\nu }_{1}}\right ){{t}_{1}}\right ]\right \}}}\\[5pt]&\,+{{{{a}_{2}}\exp \left \{2\operatorname {Re}\left [\left (\frac {{{\varGamma }_{1}}}{2{{\omega }_{0}}}-2{{\nu }_{1}}\right ){{t}_{1}}\right ]\right \}}}{{\text{e}}^{2{\rm i}\left \{x-{{\omega }_{0}}{{t}_{0}}+\operatorname {Im}\left [\left (\frac {{{\varGamma }_{1}}}{2{{\omega }_{0}}}-2{{\nu }_{1}}\right ){{t}_{1}}\right ]\right \}}}\\[6pt]=&\,{{a}_{1}}\tilde {\gamma }\cos (\varXi )+{{({{a}_{1}}\tilde {\gamma })}^{2}}\operatorname {Re}[\tilde {\beta }]\cos (\text{2}\varXi )-{{({{a}_{1}}\tilde {\gamma })}^{2}}\operatorname {Im}[\tilde {\beta }]\sin (\text{2}\varXi ), \end{aligned} \end{align}

where $\operatorname {Re}[\boldsymbol{\cdot }]$ and $\operatorname {Im}[\boldsymbol{\cdot }]$ represent the real and imaginary part of the quantity, respectively, and

(3.39) \begin{align} \begin{aligned} & \tilde {\beta }=\frac {{{a}_{2}}}{a_{1}^{2}},\quad \varXi =x-\mathrm{Re}\left [\omega _0 t_0\right ]+\operatorname {Im}\left [\left (\frac {{{\varGamma }_{1}}}{2{{\omega }_{0}}}-2{{\nu }_{1}}\right ){{t}_{1}}\right ]\!,\\[5pt]& \tilde {\gamma }=\exp \left \{\mathrm{Im}\left [\omega _0 t_0 \right ]+\mathrm{Re}\left [\left (\frac {{{\varGamma }_{1}}}{2{{\omega }_{0}}}-2{{\nu }_{1}}\right ){{t}_{1}} \right ]\right \}\!. \end{aligned} \end{align}

Here, $\tilde {\gamma }$ denotes the amplitude amplification factor. It is also noted that all variables are dimensionless.

We caution that the term $ (( {{{\varGamma }_{1}}}/{2{{\omega }_{0}}})-2{{\nu }_{1}} ){{t}_{1}}$ becomes significant only over a long period. The wave variation is determined by the combined effects of wind pressure, tangential viscous stress and phase shift. This implies that both wind pressure and tangential viscous stress must be considered together; neglecting viscous stress would lead to the omission of this term. It is worth noting that we approximate $\nu _{w}\approx \varepsilon \nu _1\approx O(\varepsilon )$ in scenarios involving large viscosity. However, in low-viscosity scenarios, where $\nu _{w}\approx \varepsilon ^3\nu _1\approx O(\varepsilon ^3)$ or $\nu _{w}\approx \varepsilon ^4\nu _1\approx O(\varepsilon ^4)$ , the effect of viscosity on the amplitude amplification factor can be neglected, resulting in $\tilde {\gamma }=\exp \left \{\mathrm{Im} [\omega _0 t_0 ]+\mathrm{Re} [({{{\varGamma }_{1}}}/{2{{\omega }_{0}})}{{t}_{1}} ]\right \}$ . If we ignore tangential viscous stress and viscosity, the wave profile simplifies to a form similar to that derived by Zdyrski & Feddersen (Reference Zdyrski and Feddersen2020). Additionally, equation (3.38) highlights the presence of a second harmonic term, $-{{({{a}_{1}}\tilde {\gamma })}^{2}}\operatorname {Im}\{\tilde {\beta }\}\sin (2\varXi )$ . When $\operatorname {Im}\{\tilde {\beta }\}\lt 0$ , this term causes a slight steepening of the wave’s forward face while flattening the rear slope, which is consistent with the findings of Longuet-Higgins (Reference Longuet-Higgins1982).