2.1. Governing equations

We solve the incompressible Navier–Stokes equations:

(2.1)$$\begin{gather} \frac{\partial{\rho}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot}(\rho \boldsymbol{u})=0, \end{gather}$$

(2.2)$$\begin{gather}\rho\left(\frac{\partial\boldsymbol{u}}{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}\right) ={-}\boldsymbol{\nabla} p +\boldsymbol{\nabla}\boldsymbol{\cdot}(2\mu\boldsymbol{\mathsf{D}}) + \sigma\kappa\delta_{S}(\boldsymbol{x} - \boldsymbol{x_{\mathcal{F}}}) \boldsymbol{n}, \end{gather}$$

(2.3)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0. \end{gather}$$

An additional scalar field representing the volume fraction of one of the two phases is introduced as $\mathcal {F}(x,y,z,t)$. The physical properties (i.e. density and the viscosity) are the $\mathcal {F}$ weighted averaged of the densities and the viscosities of the water and air phases:

(2.4a,b)\begin{equation} \rho = \mathcal{F}\rho_w+(1-\mathcal{F})\rho_a, \quad \mu = \mathcal{F}\mu_w+(1-\mathcal{F})\mu_a. \end{equation}

This together with (2.1)–(2.3) constitute the governing two-phase Navier–Stokes equations we numerically solve for. The $\mathcal {F}$ field evolves based on the continuity equation. A momentum-conserving scheme is implemented, and mass is well conserved with an error typically below 0.01 % (Mostert et al. Reference Mostert, Popinet and Deike2022).

2.2. Numerical set-up

The computation domain is of size $L_0\times L_0 \times L_0$, with four waves in the $x$ direction of wavelength $\lambda = L_0/4$ (wavenumber $k=2{\rm \pi} /\lambda =8{\rm \pi} /L_0$). The depth of the resting water is $H_w = L_0/{2{\rm \pi} }$, while the height of the airflow is $H_a = L_0(1-1/{2{\rm \pi} })$ (see figure 2). The top and the bottom are both free-slip boundary conditions, while the front and back, left and right, are periodic boundary conditions.

Figure 2. Snapshots of the air-side turbulent boundary layer and the evolving waves, for the strongest forcing case ($c/u_*=2, ak=0.2$). There are four waves in the computational domain, and the height of the water and the half-channel height for the air are shown. The colours indicate the instantaneous horizontal wind velocity $u_a$, and the surface water velocity $u_s$, respectively. The waves grow in amplitude and become short crested, which is a characteristic of wind waves. At later stage, the waves also appear to be 3-D because of the development of an underwater turbulent boundary layer. Here (a) $\omega t= 6$; (b) $\omega t= 26$; (c) $\omega t= 66$; (d) $\omega t= 86$.

We initialise the wave shape with a given surface elevation function $h_w(x,y,z, t=0)$, in this case chosen to be a third-order Stokes wave shape similar to that used in Wu & Deike (Reference Wu and Deike2021). The initial steepness $ak$ ranges from 0.1 to 0.3. The $\mathcal {F}(x,y,z, t=0)$ field is then initialised on a discretised grid based on the sign of ${y-h_w(x,y,z, t=0)}$. Here $\mathcal {F} = 1$ for the water phase ($y-h_w(x,y,z)<0$) and $\mathcal {F} = 0$ for the air phase ($y-h_w(x,y,z)>0$).

During the turbulence precursor preparation stage, the waves are kept stationary by setting $\mathcal {F}\boldsymbol {u}=\boldsymbol {0}$ at each time step. This configuration is equivalent to a turbulent boundary layer over stationary bumps. We force the turbulence with a pressure gradient (similar to a canonical channel flow), which sets the nominal friction velocity $u_*$ (i.e. total wall stress $\tau _{0}$):

(2.5)\begin{equation} \tau_{0} = \rho_{a} u_*^{2} = H_a\frac{\partial{p}}{\partial{x}}. \end{equation}

The friction Reynolds number is defined as $Re_{*} = \rho _{a}u_*H_a/\mu _{a}$ and set to $720$ for all cases. Notice that the height of the airflow is set to more than three times the wavelength $\lambda$, so that the effect of the top boundary is minimised. The physically more relevant Reynolds number is the one based on the wavelength

(2.6)\begin{equation} Re_{\lambda} = \frac{\rho_{a}u_*\lambda}{\mu_{a}}, \end{equation}

which is 214 (the ratio of the wave and the boundary layer length scales, equivalently $k\delta _{\nu }=0.029$). We use an adaptive mesh with a maximum refinement level 10 (see Appendix C for a detailed description of the AMR feature), which means that the smallest cell size is $\varDelta =L_0/2^{10}$. There are around $1.8\times 10^{7}$ grid cells in a typical simulation case, which is less than 2 % of the uniform spaced grid of the same resolution ($1024^{3}\approx 1.07\times 10^{9}$). We have validated the solver against a canonical flat wall case with $Re_*=180$ (Kim, Moin & Moser Reference Kim, Moin and Moser1987) (see Appendix C for details). The mean wind velocity profile of such a channel flow follows the law of the wall, and is similar to that of laboratory wind wave experiments (e.g. Buckley et al. Reference Buckley, Veron and Yousefi2020).

After the turbulence precursor reaches a statistically steady state, the waves are released at $t=0$ (meaning that there is no manually setting $\mathcal {F}\boldsymbol {u}=\boldsymbol {0}$ anymore, and the initial orbital velocity is added), and travel with a phase speed given by the free surface dispersion relation

(2.7)\begin{equation} c = \sqrt{g/k+\sigma k/\rho_{w}}, \end{equation}

where $g$ is the gravitational acceleration and $\sigma$ is the surface tension. The orbital velocity is initialised with the corresponding velocity field of the third-order Stokes wave (see Wu & Deike Reference Wu and Deike2021).

Since we initialise the waves with a solution of the free surface gravity wave equation, we expect the flow field to self-adjust under wind forcing during the very early stage of the simulation. The turbulent boundary layer also goes through a relaxation period when the near-wall flow adjusts to the moving boundary. We define an eddy turnover time scale $T_e = 2\lambda / \langle u \rangle (z=\lambda )$, where $\lambda$ is the wavelength and $\langle u \rangle (z=\lambda )$ is the mean horizontal velocity at vertical height $z=\lambda$. Physically it is the time scale for an eddy of size comparable to the wavelength $\lambda$ to reach equilibrium with the changing flow boundary conditions. Based on both the evolution of the wind stress and the mean profiles, we observe that the relaxation period last for approximately $4T_e$, therefore, the data of the first $4T_e$ are not included in the physical analysis. We note that any choices of initialisation will present certain limitations, as there is no exact solution of the full two-phase turbulent problem that we can use to start the simulation. After the initialisation, the waves and the turbulence interact in a fully coupled way without any prescribed interfacial conditions. The whole simulation is transient by nature, meaning that the wave amplitude changes with time, despite over a much longer time scale than both the turbulence time scale and the wave period.

The non-dimensional numbers relevant for the waves are

(2.8a,b)\begin{equation} Bo = \frac{(\rho_{w} - \rho_{a})g}{\sigma k^{2}}, \quad Re_w = \frac{\rho_{w}c\lambda}{\mu_{w}}. \end{equation}

In all the cases presented in this paper, the Bond number $Bo=200$ so that the waves are in the gravity wave regime, and we have verified that further increasing $Bo$ does not affect the results presented here (see Appendix C). The density ratio $\rho _a/\rho _w$ is set to air–water conditions $1/850$, while the viscosity ratio $\mu _a/\mu _w$ is always larger than 50 and is adjusted to set the air friction Reynolds number $Re_{\lambda }$ (2.6) and the wave Reynolds number $Re_{w}$ (2.8a,b) independently. The wave Reynolds number is kept at $Re_w \approx 10^{5}$. Note that the value of $Re_w$ gives the linear dissipation rate (per radian) due to viscosity $\gamma _d$ (Lamb Reference Lamb1993),

(2.9)\begin{equation} \gamma_d ={-}4\nu_wk^{2}/\omega = \frac{8{\rm \pi} ck}{Re_w}/\omega = 8{\rm \pi}/Re_w \end{equation}

and $D = \gamma _d \omega E$ (equivalent to (1.11)).

Notice that the velocity ratio (wave age) $c/u_*$ is varied by changing $c$, while keeping $u_*$ constant, independently of the steepness $ak$. This configuration allows us to resolve the turbulent air flow and capture the wave growth for $c/u_*$ ranging from 2 to 8 and $ak$ from 0.1 to 0.3. Table 1 summarises the simulation conditions, together with the characteristic length scales of the turbulence $\delta _\nu$ and the capillary length $l_c$, relative to wavenumber $k$ and to the smallest grid size $\varDelta$.

Table 1. A table of controlling parameters $ak$ and $c/u_*$, and relevant length scales. The third and fourth columns are the viscous wall unit $\delta _{\nu } = \nu _a/u_*$ and the capillary length scale $l_c = 2{\rm \pi} \sqrt {\sigma /(\rho _w - \rho _a)g}$ relative to $1/k$, respectively, showing the physical relevance of the parameters. They are controlled by ${Re}_{*}=720$ and ${Bo}=200$ and kept constant. The last three columns are $a$, $\delta _{\nu }$ and $k$ relative to the smallest grid size $\varDelta$, showing the numerical resolution. In the simulations, $\varDelta = L_0/2^{N}$, where $N=10$ is the maximum refinement level of the octree adaptive grid. $^{*}$For wall-modelled LES, the roughness length $kz_0$ is usually reported instead of $k\delta _{\nu }$. If we use the $z_0 = 0.11\nu _a/u_* = 0.11\delta _{\nu }$ conversion for flat smooth surface, $kz_0 = 0.003$. Also notice that these length scales are not changed when we change $c/u_*$ because $k$ is fixed, in contrast to the realistic situation, where wavenumber $k$ is smaller for fast-moving waves.

4.1. Directly observed wave growth

We quantify the growth of the waves through the time evolution of the water surface elevation $h_w (x,y,t)$, which we use to directly compute the wave energy (neglecting the surface tension energy),

(4.1)\begin{equation} E_{rms}(t) = \rho_w g \langle h_w^{2}(x,y,t) \rangle, \end{equation}

with the spatial wave averaging of a quantity $q$, in the $x$–$y$ plane, being defined as

(4.2)\begin{equation} \langle q \rangle = \frac{1}{L_0^{2}} \int_{{-}L_0/2}^{L_0/2}\int_{{-}L_0/2}^{L_0/2} q \,{\rm d}x\,{\rm d}y. \end{equation}

Figure 5 shows the time evolution of $E_{rms}(t)$ for three different $c/u_*$ cases, with initial wave steepness $ak = 0.2$. The smallest wave age case has the strongest wind forcing, and therefore the largest growth rate. The $c/u_* = 8$ case presents an almost exact balance between the wind input and viscous dissipation, resulting in a nearly constant wave energy with time. From this directly observed wave growth, we can measure a temporal rate of change of energy ${\rm d}E/{\rm d}t$ (here after we omit the subscript $rms$ for brevity). The wave growth is rather slow, and happens over O(10) wave periods and O(100) to O(1000) turbulent wall time scale $t_{\nu }=\delta _{\nu }/u_*$. This slow change in the wave energy is related to the small density ratio $\rho _a/\rho _w$, which implies weak air–water coupling (see (1.8)).

Figure 5. Wave energy normalised by initial energy $E_0 \equiv E_{rms}(t=0)$, as a function of time, directly computed from water surface height output $h_w(x,y,t)$, for three different wave ages $c/u_* = 2,4,8$ and initial steepness $ak = 0.2$. The solid curves are exponential fits to the points, although we caution that the growth rates are so small that for the exponential growth cannot be distinguished definitively from a linear growth. The $c/u_*=2$ case grows the fastest while the $c/u_*=8$ case is very slowly decaying. Note that both $E_0$ and $\omega$ change with $c/u_*$ because $g$ is changed in the numerical set-up (see § 2).

4.2. Wind surface stress

Apart from the direct surface elevation $h_w$, we extract the surface stress from the simulation. The wind stress at the surface consists of two parts, the pressure variation $\boldsymbol {\tau _p}$ (i.e. drag force) and the viscous stress $\boldsymbol {\tau _\nu }$ (see Grare et al. Reference Grare, Peirson, Branger, Walker, Giovanangeli and Makin2013; Peirson & Garcia Reference Peirson and Garcia2008),

(4.3a,b)\begin{equation} \boldsymbol{\tau_p} ={-}p_s\boldsymbol{n}, \quad \boldsymbol{\tau_\nu} = \mu_a(\boldsymbol{\nabla}\boldsymbol{u_a}+\boldsymbol{\nabla}\boldsymbol{u_a}^{\rm T})\boldsymbol{\cdot} \boldsymbol{n} = (\tau_{\nu x}, \tau_{\nu y}, \tau_{\nu z}), \end{equation}

where $\boldsymbol {u_a}$ is the air velocity vector, $p_s$ is the surface pressure, $\boldsymbol {n}$ is the normal vector of the water surface. The stresses are computed in the post-processing steps, which are independent of the computational steps. We first interpolate the velocity and pressure fields from an unstructured octree grid onto a Cartesian grid using a nearest interpolation method. The stress computation is based on the interpolated Cartesian grid. More specifically, the pressure is further interpolated onto a surface defined by $\eta + 4\varDelta$, and the mean pressure along the $x$ direction is subtracted. For the shear stress, the velocity gradients are interpolated the gradients onto the same plane, while the normal vector $\boldsymbol {n}$ of the surface is constructed by the VoF method.

Figure 6 shows the instantaneous stress fields projected onto the wave-following surface $\eta + 4\varDelta$. Since the plane is in the viscous sublayer, it is considered close enough to the actual surface that the turbulent stress can be ignored. Both the pressure and the shear stress present clear wave coherent patterns, while also having 3-D structures due to the turbulence. For example, the streaks shown in figure 6(b) are approximately $100\delta _{\nu }$ apart, which is consistent with the typical structure of wall bounded turbulent flows. There is an order of magnitude difference between the pressure and shear stress (but not their horizontal projection in (4.4)). The maximum of the pressure appears on the windward face, which is to the left of the grey line indicating the wave crest in figure 6; this gives rise to the non-zero correlation in (1.1). The viscous shear stress also reaches maximum near the wave crest due to the straining of the shear layer.

Figure 6. Instantaneous pressure (a) and the horizontal component of viscous stress (b) projected onto the wave-following surface $4\varDelta =0.1/k$ above the water surface, at $\omega t = 38$, for the case of $c/u_*=2$ and $a_0k=0.2$. Notice that there is one order of magnitude difference in the colour scale range. The grey lines show where the wave crests are. There are clearly wave coherent patterns.

From the stress field we can compute the wave-averaged integral quantities: the momentum flux (total stress $\tau _{total}$) and the energy flux (input rate $S_{total}$). The total horizontal wind stress

(4.4)\begin{equation} \tau_{total} = \langle \boldsymbol{\tau_p} \boldsymbol{\cdot} \boldsymbol{e_x} \rangle + \langle \boldsymbol {\tau_{\nu}} \boldsymbol{\cdot} \boldsymbol{e_x} \rangle = \left\langle p_s\frac{\partial h_w}{\partial x} \right\rangle + \langle \tau_{\nu x} \rangle \equiv F_p + F_s, \end{equation}

is the sum of the form drag force per unit area $F_p$ and the averaged viscous stress in the horizontal direction $F_s$. Notice that the linear approximation (${\rm d}\eta /{{\rm d} x} \ll 1$) is considered.

This stress (momentum) partition is closely related to, but different from the energy input partition. The total energy input rate by the wind stress (into both waves and underwater drift layer) is a product of the stress and the surface water velocity:

(4.5)\begin{equation} S_{total} = \langle \boldsymbol{\tau_{total}}\boldsymbol{\cdot} \boldsymbol{u_s} \rangle = \langle -p\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{u_s} \rangle + \langle \boldsymbol{\tau_{\nu}}\boldsymbol{\cdot} \boldsymbol{u_s} \rangle \equiv S_p + S_s. \end{equation}

The part of $\boldsymbol {u_s}$ that correlates with the pressure is the vertical orbital velocity $w_{orbit}$, which gives (1.1); the part of $\boldsymbol {u_s}$ that correlates with the viscous stress, however, contains both the wave horizontal orbital velocity $u_{orbit}$ and the drift velocity $u_d$,

(4.6)\begin{equation} S_s = \langle \boldsymbol{\tau_{\nu}}\boldsymbol{\cdot} \boldsymbol{u_s} \rangle \approx \langle \tau_{\nu x}{u_{sx}} \rangle = \langle \tau_{\nu x}{u_{orbit}} \rangle + \langle \tau_{\nu x}{u_{d}} \rangle \equiv S_{s,w} + S_{s,d}, \end{equation}

where $S_{s,w}$ and $S_{s,d}$ denote the energy input by the viscous shear stress into the waves and the drift, respectively. The development of the drift is discussed in Wu & Deike (Reference Wu and Deike2021), and here we focus on the energy input into the waves

(4.7)\begin{equation} S_{in} = S_p + S_{s,w} = c \left\langle p_s\frac{\partial h_w}{\partial x} \right\rangle + \langle \tau_{\nu x}{u_{orbit}} \rangle. \end{equation}

Notice that the assumption $S_p = cF_p$ means that the pressure-induced form drag contributes solely to the wave growth, while only a small variation of the viscous shear stress is correlated with $u_{orbit}$ and can contribute to the wave growth (Peirson & Garcia Reference Peirson and Garcia2008). In other words, it is not the mean stresses but the correlated part of the stresses with the wave surface velocity that contributes to the wave energy growth. In reality, $F_p$ and $F_s$ are of the same order of magnitude, but $S_p$ is generally thought to play a dominant role over $S_{s,w}$ (i.e. $S_{in} \approx S_p$), as mentioned in the introduction. We will examine both the momentum and the energy partitions using the simulation data.

In this paper we refer to the form drag $F_p$ as the wave form drag, and drag coefficient as the ratio $F_p/\tau _{total}$. Note that the wave drag force in the literature sometimes refers to the effective stress that contributes to the wave growth (from the energy flux $S_{in}$, instead of the momentum flux partition), and includes the pressure and the wave coherent viscous stress (Peirson & Garcia Reference Peirson and Garcia2008; Grare et al. Reference Grare, Peirson, Branger, Walker, Giovanangeli and Makin2013; Melville & Fedorov Reference Melville and Fedorov2015; Buckley et al. Reference Buckley, Veron and Yousefi2020, etc.),

(4.8)\begin{equation} \tau_{w} = S_{in}/c = F_p + S_{w,s}/c. \end{equation}

4.3. Wave energy growth rate versus pressure input rate

The direct wave growth and surface stress extracted from the simulation and introduced in § 4 offer two ways of computing the energy input rate into the wave $S_{in}$. First, we compute ${\rm d}E/{\rm d}t$ from figure 5 and correct for the dissipation (1.10); and second, we extract the surface pressure $p_s$ and compute the correlation (1.1).

Figure 7(a) shows a comparison of the results obtained using the two methods. The wind input rate $S_{in}(t)$ computed with (1.10) is plotted with dotted lines, and the pressure input rate $S_{p}(t)$ computed with (1.1) is with crosses, for $c/u_* = 2$ and $4$. In both cases, the pressure input $S_p$ closely follows the wave energy growth rate $S_{in}$, although there is a small gap for the $c/u_* = 2$ case. A further demonstration of the dominant role of the pressure term is shown in figure 7(b), where we plot the ratio $S_p/S_{in}$ averaged over 10 wave periods for all the cases. The ratio is very close to 1 for most cases, indicating that the pressure input $S_p$ is the major energy input term in $S_{in}$. Again, the smallest wave age cases ($c/u_*=2$) present the largest difference ($S_p/S_{in}=0.8$) and indicate that the wave coherent viscous stress might start to play a role in the strongly forced cases.

Figure 7. (a) The instantaneous pressure energy input rate $S_p=cF_p=c\langle p_s \partial h_w / \partial x\rangle$ closely follows the instantaneous wave energy growth rate (corrected with dissipation) $S_{in} = {\rm d}E/{\rm d}t+D$ for $c/u_*=2$ (dark orange) and $c/u_*=4$ (light orange); both are of $ak=0.2$. The curves are smoothed out using a moving window averaging. The variation in $S_p$ is mostly due to turbulence fluctuation. (b) The ratio between the averaged pressure energy input rate $S_p$ and the total input rate $S_{in} = (E(t_1)-E(t_2))/(t_1-t_2) + D$ computed over 10 wave periods. The ratio stays close to 1 for all the simulation cases with some variations.

Note that at high $c/u_*$, uncertainties in the budget are related to uncertainties of the decay rate for finite amplitude waves, which get amplified by the large $E$ for the fast travelling waves, together with the very small decay rate which are also hard to accurately capture numerically. Furthermore, the viscous stress input $S_{s}$ could potentially be negative for these fast travelling cases.

We want to point out that the dissipation correction is necessary in our cases, as the dissipation is non-negligible due to the limited $Re_w$. Although the wave Reynolds number $Re_w$ is constant (and therefore $\gamma _d$ is the same for different cases by (2.9)), we still have different values of $D$ for different cases of different wave frequency $\omega$ and initial energy $E_0$. The faster travelling waves have higher $E_0$ and therefore higher $D$, and the relative change in energy is much smaller. This relative change in energy (per radian) is reflected by the parameter $\gamma$ (defined in (1.8)). The underlying assumption of (1.8) is that the wave growth is exponential, and $\gamma$ represents the exponential growth rate per radian. In our simulations, we find that the growth rates are so small that for most cases, this exponential growth cannot be distinguished from a linear growth, and the growth rate computed by $\gamma ^{\prime } = S_{in}/(\omega E_0)$ shows more directly the trend of $S_{in}$. There is an uptake of $S_{in}$ as the instantaneous amplitude slowly increases over the interval of approximately 10 wave periods for the $c/u_* = 2$ case; in contrast, $S_{in}$ stays almost constant for the $c/u_* = 4$ case, as the amplitude growth is so small that its effect on $S_{in}$ is negligible.

Overall, we are able to show directly that the pressure energy input plays a dominant role in wave growth for gravity waves of realistic wave age, especially when there is a finite amplitude established ($ak\geq 0.1$). This is a different picture in terms of forcing mechanism from our previous 2-D study (Wu & Deike Reference Wu and Deike2021), where the waves are gravity-capillary waves with $ak=0.05$, and the laminar wind has a linear velocity profile with much stronger shearing ($c/u_*$ around 1).

4.4. Transient effect and microbreaking of the strongly forced cases

As we have already seen in figures 5 and 7, both the wave amplitude and the related wave-averaged quantities ($F_p$, $S_{p}$, etc.) are not stationary, especially for the small $c/u_*$ cases. Before we discuss the scaling of these quantities, we show the typical time evolution of wave form drag $F_p$ with two representative cases ($c/u_*=2$ and 8).

In figure 8 we plot the time evolution of wave form drag $F_p$ (the blue curves) as a fraction of the total wind stress, together with a few wave characteristics: the orange curves show the wave amplitude, the solid ones for the root mean square (r.m.s.) wave amplitude, defined as $a_{rms} = \sqrt {2}\langle h_w \rangle ^{1/2}$, and the dashed ones for the peak-to-peak wave amplitude, defined as half of the difference between the peak and the trough $a_{pp} = [\max (h_w)-\min (h_w)]/2$; the green curves show the curvature around the wave crest, which is taken as the minimum value of the curvature $\kappa _{min}$. These quantities are sampled at a higher frequency than that shown in figures 5 and 7.

Figure 8. Time evolution of the wave form drag and wave characteristics, namely steepness and curvature around the crest for (a) $c/u_*=2$ and (b) $c/u_*=8$. The solid orange curves and the dotted orange curves represent the steepness that corresponds to r.m.s. wave amplitude $a_{rms}$ and the peak-to-peak wave amplitude $a_{pp}$, respectively. The green curves represent the curvature around the wave crest $\kappa _{min}$ normalised by wavenumber $k$.

The $t<0$ part of the curves are from the turbulence precursor where the waves are artificially kept stationary as described in § 2. After the waves are released, there is a transitory phase where the wave form drag $F_p$ jumps up, but it soon falls back and reaches a stationary level, not far from the precursor one. As we have mentioned in § 2, we find that the transitory period lasts approximately 4$T_e$ regardless of the wave frequency $\omega$, which corresponds to the flow adjusting to the initial conditions. Consequently, we do not consider the data for $t/T_e<4$ in our analysis. The ratio of time scale $\omega T_e$ is different for different $c/u_*$, which is why the extent of the transitory period looks different. In general, $T_e$ is much smaller than the wave period.

After $4T_e$, in both cases the wave form drag $F_p$ value fluctuates due to the presence of the turbulence. What is clearly different is that in the $c/u_*=8$ case, the mean value is relatively stable while in the $c/u_*=2$ case, the wave form drag $F_p$ value keeps increasing. The significant increase in $F_p$ is related to the relatively fast wave growth, associated with an increase in the r.m.s. amplitude, as well as an increase in the non-dimensional curvature. The curvature $\kappa _{min}/k$ is taken as a measurement of how ‘sharp’ the wave crest is (but does not carry information on how 3-D the flow is). For the slowest wave case, it significantly increases above the value of the initial third-order Stokes wave and later saturates, as shown in figure 8(a). The curvature metric could be as important when determining the occurrence of airflow separation (see more in § 5.4). Around $\omega t = 80$, the underwater drift transits into turbulence, and the surface develops more prominent 3-D structure (ripples), which could also affect the wave form drag.

At later time (after around $\omega t = 90$), the r.m.s. amplitude is still increasing even though the peak-to-peak amplitude starts to plateau. The saturation is due to the microbreaking of the waves, which is shown in figure 9. This microbreaking behaviour is characterised by a confined collapse of the water surface near the crest. This coincides with a sharp increase in the wave drag $F_p$. We have only run one case for long enough time to observe the whole history of wave growing until the point of breaking. In Appendix B, we include another case with initial $ak=0.3$ and $c/u_*=2$, which exhibits a quite different behaviour, in terms of associated form drag. It does not take too long before reaching the breaking point, and the breaking is much more perceivable (spilling breakers) with droplets ejection. There is a reduction of wave form drag $F_p$ instead of an increase because of the sharp decrease of wave amplitude. A systematic study of the effect of microbreaking on the form drag within this fully coupled approach is left for future studies.

Figure 9. Microbreaking event around $\omega t=113$. A close-up view shows the microbreaking features. Initial $ak=0.2$, $c/u_*=2$. The instantaneous r.m.s. steepness $a_{rms}k=0.285$, and the instantaneous peak-to-peak steepness $a_{pp}k=0.339$.

This highlights the importance of a fully coupled approach, especially for the strongly forced condition. For the discussion in § 5, however, we focus on the effect of the initial conditions $ak$ and $c/u_*$ on the surface pressure distribution. This is done by taking a small enough averaging window after $t/T_e>4$ so that the transient effect is not prominent, and $a_{rms}(t)k$ is close enough to $ak$. For example, for the case shown in figure 8(a) we take the window of time $\omega t \in [10,30]$, and for the case shown in figure 8(b) we take the window of time $\omega t \in [25,130]$. Our results can then be compared with previous numerical studies where the motion and shape of the waves are prescribed, as well as to the experimental results. In § 6 when the $F_p$ and $S_{in}$ scalings are concerned, we bring some of the transient effect back into discussion.

5.1. Definitions

To understand better the dynamics of the wind–wave interaction and to compare with theoretical formulations introduced in § 1.3, we proceed by analysing the detailed structure of the surface pressure distribution $p_s$. This is done by extracting the principal mode of the Fourier transformation, a method that was also used in previous numerical work, see e.g. Kihara et al. (Reference Kihara, Hanazaki, Mizuya and Ueda2007) and Druzhinin et al. (Reference Druzhinin, Troitskaya and Zilitinkevich2012). The structure of the pressure field is shown in figure 6(a) and clearly contains wave-induced signals, while also being influenced by the instantaneous turbulence. To distinguish the wave-induced effect from the turbulent fluctuation, we introduce phase averaging. For any quantity $q(x,y,z)$, the phase average is

(5.1)\begin{equation} \bar{q}(\theta,z) = \frac{1}{N_{{w}} L_0}\sum_{n=1}^{N_{{w}}-1} \int_{{-}L_0/2}^{L_0/2}\,{{\rm d}y} \; q(x = \lambda(n+\theta/2{\rm \pi}),y,z), \end{equation}

where $\lambda = 2{\rm \pi} /k = L_0/4$ is the wavelength of the initial waves, and $N_{{w}}=4$ is the number of waves in the $x$ direction. The phase $\theta$ can be extracted from the surface elevation $h_w(x,y,t)$ and is therefore generalisable to cases which are not strictly sinusoidal.

The surface pressure can be generally described as the sum of different frequency modes,

(5.2)\begin{equation} p_s(\theta, t) = \sum_{n=1}^{\infty} \hat{p}_n\cos(n\theta + \phi_{pn}), \end{equation}

where $\phi _{pn}$ is the pressure phase shift and $\hat {p}_{n}$ is the pressure amplitude of mode $n$. Meanwhile, the surface elevation can be written as

(5.3)\begin{equation} h_w (\theta, t) = \sum_{n=1}^{\infty} a_n \cos (n\theta +\phi_n), \end{equation}

with $h_w (\theta, t) \approx a\cos (\theta )$ since the surface elevation $h_w$ is largely monochromatic in our simulation (and we can always shift the reference point so that the phase $\phi _1$ is zero).

Once given the surface pressure distribution $p_s$ (5.2), the wave form drag $F_p$ (4.4) becomes

(5.4)\begin{align} F_p = \left\langle p_s\frac{\partial h_w}{\partial x} \right\rangle &\approx \sum_{n=1}^{\infty} \hat{p}_{n} \,a_n nk \langle \cos(n\theta+\phi_{pn})\sin(n\theta + \phi_n) \rangle \end{align}

(5.5)\begin{align} &= \hat{p}_{1} \,ak \langle \cos(\theta+\phi_{p1})\sin(\theta) \rangle = \frac{1}{2}\,ak\,\hat{p}_1\sin(\phi_{p1}), \end{align}

and $S_p$ follows as $S_p=cF_p$. Finding the drag force and the pressure input rate now simplifies to finding the pressure perturbation amplitude $\hat {p}_1$ and the phase shift $\phi _{p1}$ that correspond to wavenumber $k$. Notice how a non-zero phase shift $\phi _p$ is necessary for a non-zero $F_p$ and $S_p$. Since (5.4) shows that only the principal mode ($n=1$) contributes to the wave growth, we then focus on how $\hat {p}_1$ and $\phi _{p1}$ depend on $c/u_*$ and $ak$ qualitatively.

5.2. Streamline and asymmetric pressure patterns

Figure 10(a–c) shows the phase-averaged vertical velocity $\bar {w}$, for three flow conditions ($c/u_*=2,4,8; ak=0.1$). The alternating patterns demonstrate the perturbation by the waves, as opposed to uniform zero for a flat surface. In the slowest wave cases (i.e. $c/u_* =2$), the alternating sign mostly comes from the straining and relaxing of the shear layer (because the airflow follows the boundary shape). In the intermediate wave speed cases ($c/u_* = 4$ and 8), the wave orbital velocity becomes significant and it leaves an imprint on the airflow (because the airflow follows the vertical motion of the boundary). Here we are plotting below $kz=3$, however, we noticed that the wave-induced perturbation in $\bar {w}$ extends higher up with increasing $c/u_*$, to almost $kz=2{\rm \pi}$, in the $c/u_*=8$ case.

Figure 10. Vertical velocity field, streamline and one-dimensional (1-D) stress distribution for three different wave ages $c/u_*=2,4,8$. In panels (a–c) (initial $ak=0.1$), the solid black lines are streamlines in the moving wave frame of reference (i.e. plotted with $\bar {w}$ and $\bar {u}-c$), and the colour shows the phase-averaged vertical velocity $\bar {w}$. Notice that the higher the $c/u_*$, the farther the wave-induced perturbation extends above the waves. Panels (d–f) (initial $ak=0.1$) are the asymmetric pressure distributions (green lines) that result from the distorted streamlines. The purple line is the shear stress. The phase shift $\phi _p$ between the pressure $p_s$ and the water surface elevation $h_w$ gives rise to the drag force and energy input. In panels (g–i) the shape of the $p_s$ distribution is consistent across different steepness, shown by different colours. The amplitude, however, seems to increase from low ($ak=0.1$) to moderate ($ak=0.15$) initial steepness, but not change much from moderate to high initial steepness ($ak=0.2,0.25$). The grey lines in all plots indicate the wave surface position, with exaggerated steepness.

In figure 10(a–c) we also plot the streamlines in the wave frame of reference (i.e. with $\bar {w}$ and $\bar {u} - c$). There are recirculation cells because the vertical velocity is of alternating signs, and the horizontal velocity changes sign at some height. This height is often called the critical height, and it depends on the value of $c/u_*$. The higher $c/u_*$ is, the farther away the critical height is from the water surface. These recirculating cells influence the pressure variation $p_s$ at the water surface in a complicated way, which is plotted in figure 10(d–f) with green lines.

Figure 10(d–f) is the averaged stress distribution for both the pressure and the viscous shear stress (shown in figure 6). We see clearly that the pressure maximum is on the windward face, and the phase shift is marked by $\phi _{p1}$. Notice that even for the smallest steepness case ($ak=0.1$), the shapes of the pressure distribution are not sinusoidal. For example, at $c/u_*=2$, the trough of the pressure signal is rather flat, which is a sign of sheltering (with or without a certain level of airflow separation). For the $c/u_*=8$ cases, the pressure distribution is tilted forward. Only in the $c/u_* = 4$ case does the pressure distribution roughly resemble a sinusoidal wave. The pressure structures are in qualitative agreement with those found in simulations (Kihara et al. Reference Kihara, Hanazaki, Mizuya and Ueda2007; Yang & Shen Reference Yang and Shen2010) and experiments (Mastenbroek et al. Reference Mastenbroek, Makin, Garat and Giovanangeli1996) with corresponding $c/u_*$. The non-sinusoidal pressure shape is the signature of higher frequency modes and would contribute to the growth of corresponding wave frequencies.

The bottom row shows how the 1-D pressure distribution changes with different initial $ak$, ranging from $0.1$ to $0.25$ (colour coded). The shapes are similar for the same $c/u_*$, with the amplitude of the pressure variation increasing with wave steepness $ak$. The largest difference is between $ak=0.1$ and the other three $ak$ values, where the amplitude of the pressure seems to saturate at high $ak$.

5.3. Pressure amplitude and phase shift

Figure 11 shows the pressure amplitude $\hat {p}_1$ and phase shift $\phi _{p1}$ as a function of $c/u_*$ and $ak$. These quantities are computed by Fourier transform of the phase-averaged surface pressure $p_s$. The ‘surface’ is defined as the wave-following surface $4\varDelta = 0.1/k$ away from the air–water interface. We have tested the sensitivity to the location within the first eight grid points and it does not present much difference (as long as we are in the viscous layer).

Figure 11. (a) Pressure amplitude $\hat {p}_1$ normalised by the nominal wall stress $\tau _0 = \rho _a u_*^{2}$, and in addition $ak$, plotted against $c/u_*$. (b) Pressure phase shift $\phi _{p1}$ as a function of $c/u_*$. Notice that because of the $F_p=(1/2)\hat {p}_1\,ak\sin (\phi _{p1})$ relation, the drag force is the largest when $\phi _p = 90^{\circ }$, and zero when $\phi _p = 0^{\circ }$ or $180^{\circ }$. The results from Kihara et al. (Reference Kihara, Hanazaki, Mizuya and Ueda2007) of $Re_{\lambda } = 161$ and $ak=0.1$ are plotted with black crosses. The results from Druzhinin et al. (Reference Druzhinin, Troitskaya and Zilitinkevich2012) of $Re = 10\,000$ and $ak=0.2$ are plotted with black plus signs. (c) The wave form drag $F_p$ is not a strong function of $c/u_*$ for all values of the steepness $ak$. We also show the full integral value $\langle p_s \partial {h_w}/\partial {x}\rangle$ in comparison with the single mode representation $(1/2)\hat {p}_1ak\sin (\phi _{p1})$. The markers and colours are the same with those in figure 7(b) and 10.

Figure 11(a) shows that the amplitude $\hat {p}_1$ first increases with $c/u_*$ until $c/u_*\approx 6$ and then decreases, for all steepness $ak$. Figure 11(b) shows that the phase shift $\phi _{p1}$ follows the opposite trend. The net result is that the drag force shown in figure 11(c) is not a strong function of $c/u_*$, which is in agreement with previous studies in the slow wave regime. Figure 11(c) also confirms (5.4): the dotted and solid lines show the single mode representation and the integral representation of wave form drag $F_p$, respectively, which agree very well, even when the pressure distribution is not necessarily sinusoidal.

Taking a closer look at the phase shift, it is around $90^{\circ }$ for the strongest forcing cases $c/u_*=2$, and then goes under $90^{\circ }$ between $c/u_*=2$ and $6$, and then slightly above $90^{\circ }$ at $c/u_*=8$. This indicates that the sheltering mechanism is dominant in the strong forcing conditions, and that the theories based on linear stability analysis might be at work in the higher wave age cases (more on this in § 7.1). Good agreement was found with results from Kihara et al. (Reference Kihara, Hanazaki, Mizuya and Ueda2007) (marked with black crosses) at $ak=0.1$. The configuration in Kihara et al. (Reference Kihara, Hanazaki, Mizuya and Ueda2007) is similarly a pressure-driven channel flow, with $Re_{\lambda }=161$. Data from Druzhinin et al. (Reference Druzhinin, Troitskaya and Zilitinkevich2012) show larger $\phi _{p1}$ for all wave age although the trend with wave age is similar. They used a bulk Reynolds number of $10\,000$ with a Couette flow configuration. We could not infer the exact friction Reynolds $Re_{\lambda }$ value from the information provided in the paper, but they should be of the same order of magnitude as in our set-up.

The pressure amplitude $\hat {p}_1$ is normalised by $ak\tau _0$ in figure 11(a), and this choice is made by a commonly adopted scaling argument. Intuitively, and also used in the theoretical studies mentioned in § 1.3, the pressure variation amplitude should scale with $\rho _a (ak) U_{ref}^{2}$, with $U_{ref}$ being some characteristic wind velocity (not necessarily the friction velocity $u_*$). From figure 11(a) we see that this scaling does not collapse $\hat {p}_1$ with respect to $ak$, at least not when $u_*$ is used. Now defining the ratio between $\hat {p}_1$ and $ak \rho _a u_*^{2}$ as $P$, i.e.

(5.6)\begin{equation} P = \hat{p}_1/ak \rho_a u_*^{2}. \end{equation}

This ratio $P$ represents $(U_{ref}/u_*)^{2}$, the ratio between the should-be characteristic velocity $U_{ref}$ and the friction velocity $u_*$. From figure 11(a) we see that $P$ ranges from around 15 to 45, indicating that $U_{ref}/u_*$ is around 4 to 7.

5.4. A note on airflow separation and microbreaking for steep waves

We observe intermittent airflow separation when the wave steepness $a_{rms}k$ reaches a value between around 0.23 to 0.27, for the $c/u_*=2$ and 4 cases. In figure 12(a–c) we show examples of the instantaneous horizontal velocity at the centre slice ($y=0$), for three wave ages $c/u_*=2,4,8$, when the $a_{rms}k$ steepness is around 0.24. There is airflow separation and recirculation for the $c/u_*=2$ case, indicated by the confined negative $u$ zone. In contrast, there is no such airflow separation for the $c/u_*=4, 8$ cases (or much rarer throughout the whole domain), suggesting that the increased regular wave induced motion (which scales with $akc$) near the bottom boundary could suppress the airflow separation. The airflow separation in the $c/u_*=2$ case is, however, highly intermittent and localised. In fact, for the phase-averaged horizontal velocity shown in the second row, the airflow separation is not distinguishable.

Figure 12. The instantaneous horizontal velocity (a–c) at $y=0$ and the phase averaged horizontal velocity (d–f) in laboratory frame of reference, for three different wave ages at comparable instantaneous steepness $a_{rms}k$. There is airflow separation for the $c/u_*$ case but the effect is intermittent localised. For the $c/u_*=4$ and $8$ cases there is no separation at similar steepness. Panels (d–f) show that the phase-averaged flow field $\bar {u}_a$ is similar in effect to that of a non-separated case.

Notice that the microbreaking mentioned in § 4.4 does not occur until $a_{rms}k \approx 0.3$. This means that airflow separation can occur before the waves break, due to the sharp directional change of the lower boundary, when the waves are steep and short-crested. This finding is consistent with the previous experimental and numerical findings (Donelan et al. Reference Donelan, Babanin, Young and Banner2006; Yang & Shen Reference Yang and Shen2010; Druzhinin et al. Reference Druzhinin, Troitskaya and Zilitinkevich2012; Sullivan et al. Reference Sullivan, Banner, Morison and Peirson2018b; Buckley et al. Reference Buckley, Veron and Yousefi2020). We comment that it is, however, not practical to determine an exact steepness value of $a_{rms}k$ at which separation starts to occur. It is likely also dependent on other geometric quantities such as $\kappa _{min}/k$, as they are a more local measure of the change of direction in the boundary, and therefore closely related to vorticity generation in the boundary layer (Batchelor Reference Batchelor2000). In fact in Buckley et al. (Reference Buckley, Veron and Yousefi2020) figure 16, the likelihood of airflow separation is reported experimentally, and it increases with steepness but decreases with wave age, which is what we observe as well. We caution that the occurrence of separation and the exact separation point can also depend on the Reynolds number of the flow, which is much lower in the DNS than in realistic wind wave airflow.

Nonetheless, the onset of airflow separation does not significantly affect the discussion in § 5.3. In fact, if we consider the phase-averaged velocity and surface pressure, the separated and non-separated sheltering cases exhibit similar features. That is to say, even the separating cases can be readily incorporated into the current framework of principal $p_s$ mode analysis. Although it is possible that the separation point might shift the phase $\phi _{p1}$, and this explains why for $c/u_*=2$, the $ak=0.2$ and $0.25$ cases have different $\phi _{p1}$ from the $ak=0.15$ and $0.1$ cases (see figure 11b).

6.1. Wave drag $F_p/\tau _0$

In § 5.3, we have shown that the drag is not a strong function of $c/u_*$ in the slow wave regime. However, it is strongly dependent on the steepness. Instead of showing the wave form drag $F_p$ as a function of initial $ak$, figure 13 shows the drag coefficient $F_p/\tau _0$ as a function of the r.m.s. steepness $a_{rms}k$. Since for the small wave age cases (the green dots), there is a significant increase of $a_{rms}k$ due to the wind forcing, we take multiple averaging windows. The points that belong to the same case are connected with a line. The bar in the $x$ axis is the range of $a_{rms}k$ in the averaging time window. The bar in the $y$ axis is the standard deviation of $F_p$ fluctuation, which is mostly due to the turbulent fluctuation (see figure 8).

Figure 13. The wave form drag $F_p$ (or in some cited works wave drag $\tau _w$ defined by (4.8)) as fraction of the total stress $\tau _0$, plotted as a function of r.m.s. steepness $a_{rms}k$. For the $c/u_*=2$ cases (green points), we take multiple averaging windows because of the transient evolution of $F_p$. The bar in the $x$ axis is the range of $a_{rms}k$ in the averaging time window. The bar in the $y$ axis is the standard deviation of $F_p$ fluctuation (mostly due to turbulent fluctuation). Points that belong to the same initial $ak$ case are connected with a line. Other numerical data: stars Kihara et al. (Reference Kihara, Hanazaki, Mizuya and Ueda2007), $c/u_*=2,4,8$, mostly overlapping with the $ak=0.1$ results; pentagons, Yang & Shen (Reference Yang and Shen2010), $c/u_*=2$. Experimental results: solid circles, Peirson & Garcia (Reference Peirson and Garcia2008); solid crosses, experimental observation from Mastenbroek et al. (Reference Mastenbroek, Makin, Garat and Giovanangeli1996); plus signs, numerical prediction from Mastenbroek et al. (Reference Mastenbroek, Makin, Garat and Giovanangeli1996); solid diamonds, Banner (Reference Banner1990); open circles, Banner & Peirson (Reference Banner and Peirson1998); light crosses, Grare et al. (Reference Grare, Peirson, Branger, Walker, Giovanangeli and Makin2013); open diamond, Buckley et al. (Reference Buckley, Veron and Yousefi2020); open squares, Funke et al. (Reference Funke, Buckley, Schultze, Veron, Timmermans and Carpenter2021). The last three data sets denoted with open marks are purely wind generated waves, and the Grare et al. (Reference Grare, Peirson, Branger, Walker, Giovanangeli and Makin2013) data set has mixed types, while the others are all mechanically generated waves (or similar numerical set-ups). The Banner (Reference Banner1990) and the Banner & Peirson (Reference Banner and Peirson1998) datasets include waves with microbreaking. Dashed line: the quadratic representation $F_p = 1/2\beta (a_{rms}k)^{2}$ with a constant $\beta$; solid line: the Belcher correction (6.2).

Again, if we analyse the first data point of each green dots group and the blue and purple dots of the same initial $ak$, they are close to each other, as we have found in § 5.3, meaning that the ratio $c/u*$ has little effect on the wave form drag. For the small steepness regime ($ak<0.2$), the data roughly scales with $(ak)^{2}$, with some small variation in different $c/u_*$,

(6.1)\begin{equation} F_p \sim (ak)^{2} \tau_0. \end{equation}

More specifically the prefactor is $1/2P\sin (\phi _p)$, with $P$ defined in (5.6). For higher steepness $ak=0.2, 0.25$, we see a plateau in $F_p/\tau _0$ and a departure from the $(ak)^{2}$ scaling, and slightly larger variation with $c/u_*$.

However, if we trace an individual case of $c/u_*=0.2$, the picture is further complicated. For example, for the $ak=0.2$ case, the $F_p$ value undergoes stages of growth and increases from $0.25\tau _0$ to around $0.7\tau _0$ over the course of $a_{rms}k$ from 0.2 to 0.27. The time evolution of the strongly forced cases overall seems to better fits the $(a_{rms}k)^{2}$ scaling than the ensemble of cases, which falls shorts of the $(ak)^{2}$ scaling. There also seems to be a wave history effect: for example, the initial $ak=0.15$ case shows higher value of $F_p/\tau$ when $a_{rms}(t)k$ reaches 0.2, when compared with the case with initial $ak=0.2$. This is probably related to the wave geometry and its short-crestness that evolves with the amplitude growth, that is not captured by only varying the initial amplitude for the Stokes waves.

Figure 13 also shows numerical and experimental data from the literature. If only considering the initial $ak$ effect, our results agree very well with previous numerical studies across different $ak$. The $ak=0.1$ result is very close to those from Kihara et al. (Reference Kihara, Hanazaki, Mizuya and Ueda2007), and the $ak=0.25$ results are within the range of those reported by Yang & Shen (Reference Yang and Shen2010). Note that these simulations are performed with prescribed wave boundary shape and motion. This agreement serves as a further validation for the current numerical method. On the other hand, it suggests that the one-way coupled approach could suffice for predicting the wave form drag of weakly coupled cases where the waves’ growth is very slow. The necessity of the fully coupled approach comes when the waves are strongly forced, grow relatively fast and exhibits strong nonlinear behaviour such as short-crestness and microbreaking.

For comparison with experimental studies, we note that some of the data plotted in figure 13 are actually $\tau _w$ defined by (4.8) instead of $F_p$. Since we have already verified that the pressure is responsible for over 80 % of the energy flux, the $F_p$ and $\tau _w$ values do not differ by much for the cases discussed here; at least the small difference does not affect the general trend of $F_p/\tau _0$ with increasing $ak$.

Peirson & Garcia (Reference Peirson and Garcia2008) (solid circles) measured $\tau _w$ by the spatial wave energy growth, and their data match with ours quite well. They also suggested a correction to the $(ak)^{2}$ relation inspired by Belcher (Reference Belcher1999), with two fitted parameters $\beta _f$ and $\beta _t$,

(6.2)\begin{equation} \tau_w/\tau_0 = (\beta_f+\beta_t)(ak)^{2}/[2+\beta_f (ak)^{2}], \end{equation}

that seem to fit a compilation of the data sets well (see their figure 5). This correction is plotted with the solid line in figure 13. The other experimental studies have reported a wave drag coefficient somewhat higher. Mastenbroek et al. (Reference Mastenbroek, Makin, Garat and Giovanangeli1996) measured the wave drag coefficient by using a fixed pressure probe at a fixed height $kh={\rm \pi}$, and Grare et al. (Reference Grare, Peirson, Branger, Walker, Giovanangeli and Makin2013) used PIV viscous stress measurement, pressure with fixed or wave-following probe for different subset data. Buckley et al. (Reference Buckley, Veron and Yousefi2020) and Funke et al. (Reference Funke, Buckley, Schultze, Veron, Timmermans and Carpenter2021) were obtained from the same data set; Buckley et al. (Reference Buckley, Veron and Yousefi2020) used PIV viscous stress measurement and computed pressure as a stress residual, while Funke et al. (Reference Funke, Buckley, Schultze, Veron, Timmermans and Carpenter2021) reconstructed the pressure field by solving the Poisson equation.

From the synthesis of data, we can see that the numerical estimations of $F_p/\tau _0$ are in general lower than the experimental measured ones. Since we have seen that there is a wave history effect which is related to the evolving wave shape, it explains why numerical simulations with more idealised wave shape (Airy waves or Stokes waves) might be missing that effect and therefore predicting lower wave form drag $F_p$. Yang & Shen (Reference Yang and Shen2010) noticed that nonlinearity can play an appreciable effect by comparing their Airy waves and Stokes waves results. The current work further shows that the steep wave shape can deviate even more from the Stokes waves and increase the wave form drag. However, there remains significant scatters within the experimental data using different methods to measure the stress. The ones inferred from the wave growth seem to be consistently lower than the ones measured from the air stress in experiments, and the differences are beyond the scatters that might be introduced by different wave ages. Although we have verified using our simulation that the measurement of $F_p$ directly from the air stress or indirectly from the wave growth should be consistent, there remain a few possible reasons for the scatters in the experimental data: one is the existence of 3-D smaller scale waves (roughness elements) that increases the drag; the other is the uncertainty caused by the air-side measurement, especially the pressure extrapolation error from a finite height to the surface, discussed in Donelan et al. (Reference Donelan, Babanin, Young and Banner2006) and Grare et al. (Reference Grare, Peirson, Branger, Walker, Giovanangeli and Makin2013). A further examination of the extrapolation error will require a study of the vertical pressure structure.

6.2. Growth rate $\gamma$

The energy input by pressure is closely linked to the wave form drag by $S_p = cF_p$; or considering the more general definition of wave drag (4.8), the total wind input is $S_{in} = c\tau _w$. The two are used interchangeably in the present discussion, i.e. $S_{in} = c\tau _w \approx cF_p$.

We have seen that the drag force $F_p$ is not a strong function of $c/u_*$, so the pressure energy input rate $S_{p} = cF_p$ increases with $c/u_*$ as shown in the inset of figure 14, i.e. in the slow wave regime, the energy flux is higher for waves travelling faster (at a fixed $u_*$). This could appear in contradiction to the observation that the slowest travelling waves have the fastest growing energy curve in figure 5. This is, however, not self-contradicting, because the curves in figure 5 reflect the relative rate of change of energy, which is $S_{in}$ further normalised by the total energy $E$, and $E$ is larger for faster waves. (Note that in figure 5, since we consider the net energy growth, another factor that is the viscous decay is also larger for the faster waves, since $\gamma _d$ is constant in our simulation.)

Figure 14. Non-dimensional growth rate scaling computed with the points in figure 11. Inset figure shows the energy input rate $S_p$ increasing with increasing $c$, while the main plot shows the non-dimensional growth rate ($\gamma$) decreasing with increasing $c/u_*$. The grey line is the $(u_*/c)^{2}$ fitting. It demonstrates that the non-dimensional growth rate scaling is dominated by the $\omega E$ normalisation.

This normalisation by the total energy $E$ and angular frequency $\omega$, i.e. the definition of growth rate per radian $\gamma = S_{in}/(\omega E)$, was introduced by Miles (Reference Miles1957), and is based on the assumption that the growth is exponential. Considering the definitions of wave energy and the gravity wave dispersion relation,

(6.3a,b)\begin{equation} E = \frac{1}{2}\rho_w g a^{2}, \quad \omega = kc = \sqrt{gk} \end{equation}

and using the assumption that $F_p \sim (ak)^{2}$ (which we have seen to be questionable at high $ak$), and by introducing the prefactor $\beta$ (Miles Reference Miles1957), we obtain

(6.4)\begin{equation} F_p = \frac{1}{2} \beta (ak)^{2} \tau_0 = \frac{1}{2} \beta (ak)^{2} \rho_a u_*^{2}, \end{equation}

which becomes

(6.5)\begin{equation} \gamma = \frac{S_{in}}{\omega E} = \frac{cF_p}{\omega E} = \beta \frac{\rho_a}{\rho_w}\left(\frac{u_*}{c}\right)^{2}. \end{equation}

It is worth noticing that this relationship, widely used in the literature, presents some strong self-correlation between the normalisation of $S_{in}$ by $\omega$ in the left-hand side and the phase speed $c=\omega /k$ on the right-hand side. The resulting $(u_*/c)^{2}$ scaling is reflected in figure 14.

The representation of (6.5) in figure 15 is often taken as an indirect proof of Miles’ theory. Plant (Reference Plant1982) compiled laboratory and field measurements known to that date (plotted in grey symbols in figure 15), which became the benchmark and established the $(u_*/c)^{2}$ scaling, although the empirical range of $\beta$ (indicated in grey dotted lines) is higher than the original prediction from Miles (Reference Miles1957).

Figure 15. Growth rate parameter $\gamma$ as function of inverse wave age $u_*/c$. The value of $a_{rms}k$ is denoted with the colour scale. Notice that the added averaging windows in figure 13 result in more points for the $c/u_*=2$ cases, but the $\gamma$ values are very close to each other, due to the fact that the time evolving $F_p(t)$ scales with $(a_{rms}(t)k)^{2}$ relatively well. The points from cited works are also colour-coded whenever the steepness value can be identified. Numerical works: blue crosses, Yang et al. (Reference Yang, Meneveau and Shen2013) with JONSWAP spectrum; open triangles, Kihara et al. (Reference Kihara, Hanazaki, Mizuya and Ueda2007), $ak=0.1$. Experimental works: open diamonds, Buckley et al. (Reference Buckley, Veron and Yousefi2020), $ak$ values as the colours indicate; grey symbols: data compiled by Plant (Reference Plant1982) with no steepness information. Dotted lines, the range of $\beta$ proposed by Plant (Reference Plant1982) based on empirical evidence.

We caution that while the $(u_*/c)^{2}$ scaling seems to hold, there is a wide scatter in the $\beta$ value at a given value of $u_*/c$, with sometimes is over an order of magnitude variation. We also note that alternatives for the reference velocity have been proposed (e.g. the sheltering coefficient at half-wavelength by Donelan et al. (Reference Donelan, Babanin, Young and Banner2006) or the middle-layer velocity from Belcher (Reference Belcher1999)), and the reported values of the $\beta$ parameter by experimental and numerical studies could be presented in terms of another reference velocity, leading to estimations of the sheltering coefficient (see Peirson & Garcia (Reference Peirson and Garcia2008) and Yang et al. (Reference Yang, Meneveau and Shen2013), for example).

A large contributing factor to the scatter is the role of the wave steepness at a given wave age, as already discussed by Peirson & Garcia (Reference Peirson and Garcia2008) and Buckley et al. (Reference Buckley, Veron and Yousefi2020). The steepness is indicated in figure 15 with different shades of red for the data sets where the wave steepness can be identified. As we have mentioned, the assumption that the wave form drag scales with the steepness $(ak)^{2}$ does not hold for moderate to high steepness ($ak>0.15$).

The other factor is again the uncertainty in the pressure-slope correlation (1.1) measurements. The data sets compiled by Plant (Reference Plant1982) were all obtained by measuring the aerodynamic pressure, with either fixed or wave-following probes. This is to some extent due to the difficulty in directly measuring the wave growth as an alternative: for the fast-moving waves, measuring the extremely small growth in amplitude is prone to errors; and for the less controlled field campaigns, it is hard to single out the wind input from the nonlinear interactions and dissipation. It is of crucial importance, therefore, that we find ways to quantity the uncertainties in these pressure measurements.

In summary, the $(u_*/c)^{2}$ scaling in figure 15, despite being robust because of the normalisation, inherits the uncertainty reflected in figure 13. The normalisation of $S_{in}$ by $\omega$ and $E$ following (1.8) is questionable with the growth rate being very small due to the small density ratio $\rho _a/\rho _w$ so that the exponential growth cannot be verified in a convincing way; and the normalisation makes the $\gamma$ parameter too skewed by the wave characteristics.

We want to mention that it remains to be studied how the results from the current study and the other lab experiments with nearly monochromatic wave trains can be extended to broadband ocean wave spectra. The method to date (Snyder et al. Reference Snyder, Dobson, Elliott and Long1981; Donelan et al. Reference Donelan, Babanin, Young and Banner2006; Yang et al. Reference Yang, Meneveau and Shen2013) is to keep the linear assumption, and the correlation term (1.1) becomes the cross-spectrum

(6.6)\begin{equation} Q(\omega) = \langle p_s(\omega)h_w(\omega)^* \rangle. \end{equation}

Interestingly, the numerical study of a broad-spectrum wave field from Yang et al. (Reference Yang, Meneveau and Shen2013) (blue crosses) reported growth rate of very similar magnitude to our study. The numerical methods are very different: the points from Yang et al. (Reference Yang, Meneveau and Shen2013) are from computing (6.6) in one run for different wave frequency $\omega$; while the points in our study are from different runs with different initial $c/u_*$ and $ak$. The steepness $a(\omega )k$ is not reported in Yang et al. (Reference Yang, Meneveau and Shen2013), therefore it is hard to draw a definite conclusion.