3.1. Wave interface evolution

Figure 2 shows the wind-wave system for the case $u_\ast /c=0.9$ at four characteristic times. At $\omega t = 20$ (figure 2 a), the wave field has grown due to the wind input and is close to breaking conditions (incipient breaking). While the initial wave field is nearly monochromatic, the waves become sharp-crested as they grow under wind forcing and approach breaking. Once the breaking stage is concluded, around $\omega t=41$ (figure 2 b), the wave field grows again, starting with a smaller steepness (figure 2 c) taken at $\omega t = 70$ and finally breaks again at $\omega t = 150$ (figure 2 d). In this last condition, the wave field becomes clearly three-dimensional, with energy redistribution to higher wavenumber following the breaking and growing cycles, and intermittent microbreaking events visible. Despite this, the ratio between the effective main wavelength and the initial wavelength remains nearly unchanged, with wave elevation analysis showing variations of less than 6 % throughout the simulations (not shown). Moreover, most of the wave energy remains concentrated around the peak wavelength set by the initial conditions, as constrained by the boundary conditions. The peak downshift is effectively limited by these constraints, with energy redistribution to higher wavenumbers primarily occurring following breaking events.

3.2. Time evolution of the potential wave energy

We consider the potential energy of the wave, $E_W(t)$ , to characterize the evolution of a wave field and distinguish between the growing (i.e. $E_W$ increases) and the breaking stages (i.e. $E_W$ decreases):

(3.1) \begin{equation} E_W(t) = \rho _w|\textbf {g}|\int _{\Omega _w} z \,{\rm d}V - E_{p,0}\textrm {,} \end{equation}

where the integration volume in the water $\Omega _w$ is done up to the wave surface $\eta$ in the z direction. Parameter $E_{p,0}$ is the potential energy when the wave field is undisturbed and can be evaluated as $E_{p,0}=\rho _w|\textbf {g}|L_0^2\int _0^{h_W} z\,{\rm d}z = \rho _w|\textbf {g}|(h_WL_0)^2/2$ . A contribution due to surface tension exists in $E_W(t)$ but is negligible here due to the large $Bo$ . Note that we calculate $E_W(t)$ by integration in the water volume over the vertical coordinate $z$ (Lamb Reference Lamb1993) since it is also valid for nonlinear and breaking waves contrary to the formula based on linear approximation $E_W(t)=\rho _w|\textbf {g}|\int _\Gamma \sqrt {2(\eta -\overline {\eta })^2} {\rm d}S$ (Lamb Reference Lamb1993) which is not well defined during the breaking stage.

Figure 3. Wind-wave growth and breaking life cycle. (a) Evolution of the normalized potential wave energy $E_W/E_{W,0}$ with $E_{W,0}=E_W(t=0)$ for increasing $u_\ast /c$ from $0.3$ to $0.9$ as a function of the dimensionless time $\omega t$ where $\omega =2\pi /T_0$ is the angular frequency and $T_0$ is the wave period. The associated instantaneous steepness $ak$ values are shown on the second $y$ axis. (b) Sketch of the characteristics of dynamical regimes observed in the simulations, illustrated for $u_\ast /c=0.9$ . Here $G_{1,2}$ represent the first and the second growing stages; $B_{1,2}$ represent the first and the second breaking stages; $F$ is the final stage; and $G_{2,a}$ and $G_{2,b}$ are the fractions of the second growing stages with equal time windows of stages $G_1$ and $F$ . These windows are used to compute averages for the momentum flux and drag coefficient during growth and breaking.

Figure 3(a) depicts the time evolution of the potential wave energy normalized by its initial value, $E_W/E_{W,0}$ (left axis), highlighting the growth and breaking stages of the wave field under strong wind forcing. For reference, figure 3(a) also presents the instantaneous wave steepness $ak$ (right axis), calculated from the wave amplitude as $a(t) = \sqrt {(2/\Gamma )\int _\Gamma (\eta - \overline {\eta })^2{\rm d}S}$ . Alternatively, the wave steepness can be determined using the surface elevation derivatives as $\mathcal {S} = \max (\sqrt {(\partial \eta / \partial x)^2 + (\partial \eta /\partial y)^2})$ . Note that under the assumption of linear wave theory, $\mathcal {S}$ reduces to $\mathcal {S} = a(t)k$ . In this study, both methods were evaluated, and despite the assumptions of linear theory being unmet during breaking, they produced very similar results.

The time evolution of wave energy $E_W/E_{W,0}$ and the wave steepness $ak$ is controlled by the competing effects of the wind forcing $u_\ast /c$ and the dissipation. In each curve, the potential wave energy displays small oscillations of period $T_0/2$ , already observed in Iafrati (Reference Iafrati2009) and Deike et al. (Reference Deike, Popinet and Melville2015) in similar configurations without wind.

During the growing stage, the dissipation is viscous and controlled by the water-wave Reynolds number $Re_W$ . For the smaller forcing $u_\ast /c=0.3$ , the wind energy input and the viscous dissipation are in balance. Note that in ocean water waves, the wave Reynolds number is significantly higher and, therefore, the wave field would grow for $u_\ast /c=0.3$ . Performing simulation at higher wave Reynolds number would also lead to wave energy growing in time for such value of $u_\ast /c=0.3$ (see Wu & Deike Reference Wu and Deike2021; Wu et al. Reference Wu, Popinet and Deike2022; Zhang et al. Reference Zhang, Hector, Rabaud and Moisy2023). The wave growth rate increases with $u_\ast /c$ and for all stronger forcing, $E_W(t)$ grows until it reaches a critical point (breaking conditions), when part of the wave energy is dissipated. Once the breaking stage concludes, $E_W(t)$ grows again, and a second breaking event can occur at later times. When $u_\ast /c$ increases, the wave field reaches a critical amplitude for breaking $(ak)_c$ earlier. The observed critical breaking steepness $(ak)_c$ varies slightly from $(ak)_c\in [0.28\,{-}\,0.34]$ between the different wind forcing and first and second cycles, close to the values reported without wind for similar initial conditions (Deike et al. Reference Deike, Popinet and Melville2015, Reference Deike, Melville and Popinet2016) and experiments on focusing wave packets (Drazen et al. Reference Drazen, Melville and Lenain2008). The magnitude of the breaking event can be quantified by the maximum energy loss, $(\max (E_W)-\min (E_W))$ , which increases with the amplitude of breaking.

In all cases, after one of two breaking events, the wave field reaches a final stage with limited variation in the wave energy. This condition can be attributed to the balance between the wind forcing and the wave field characterized by a broader wave spectrum and intermittent microbreaking events occurring at some of the wave crests (see figure 2 d). Such conditions reduce the effectiveness of wind input in promoting wave growth and lead to a quasi-steady condition with wind input approximately balanced by energy dissipation. Note that the loss of energy of these microbreaking events is smaller since the steepness at breaking is smaller (see e.g. Drazen et al. Reference Drazen, Melville and Lenain2008; Deike et al. Reference Deike, Popinet and Melville2015).

For the purposes of the present analysis, we define different time windows that characterize the physics at play made of the growing and breaking stages, which are illustrated in figure 3(b). We consider the first and second growing stages, $G_1$ and $G_2$ (when $E_W$ increases). Each growing stage is followed by a corresponding breaking stage, $B_1$ and $B_2$ (when $E_W$ decreases). Note that for the highest $u_\ast /c=0.7\,{-}\,0.9$ , we run the simulations long enough to observe a breaking stage, followed by a quasi-stationary final state, $F$ , where $E_W$ exhibits a limited variation. In the second growing stage, we define two further time windows $G_{2,a}$ and $G_{2,b}$ . These time windows, together with $G_{1,2}$ , $B_{1,2}$ and $F$ , are useful to define the mean velocity profiles, perform momentum flux and drag coefficient analysis during the growing and breaking stages and determine the respective role of each regime in the averaged quantities. Note that $G_1$ - $G_{2,a}$ and $G_{2,b}$ - $F$ are chosen to have the same duration in time, so that the averaged quantities defined therein are computed over comparable temporal windows.

3.3. Airflow above growing and breaking waves

The flow structure near the interface is heavily modified by the presence of the waves. Figure 4(a–c) shows the instantaneous velocity contours of the streamwise air velocity in a reference frame moving at the wave speed, i.e. $u-c$ , in the region close to the wave surface ( $-0.2\leqslant z\leqslant 0.4\lambda$ ) and sampled at the spanwise mid-plane ( $y/\lambda =0$ ), at several characteristic times (during growth and breaking) for $u_\ast /c=0.9$ .

Figure 4. Illustration of airflow separation in strongly forced steep waves, just before breaking (a,d,g), during breaking (b,e,h) and during the second growing stage (c,f,i). (a–c) Contours of the streamwise instantaneous velocity in the airflow sampled at the middle plane $y/\lambda =0$ . The streamwise velocity is plotted in a reference frame moving with the wave and is normalized by the friction velocity, i.e. $(u-c)/u_\ast$ . (d–f) Contours of the instantaneous spanwise vorticity in the airflow sampled at the middle plane $y/\lambda =0$ and normalized by the angular velocity, i.e. $\Omega _y/\omega$ . (g–i) Contours of the spanwise vorticity in the airflow, averaged along the spanwise direction and normalized by the angular velocity, i.e. $\Omega _y/\omega$ . All the panels are plotted in the region $-0.2\lambda \leqslant z \leqslant 0.4\lambda$ for $\omega t=[20,35,70]$ for $u_\ast /c=0.9$ .

As the wave field grows and approaches breaking, i.e. $\omega t=20$ in figure 4(a), large negative values appear and localize not only near the wave troughs but also close to the wave crests. In these regions, the flow reverses and recirculates, with negative velocity up to $1.5u_\ast$ . Note that flow reversal starts to occur before the breaking stage, as the wave becomes steep and short-crested. During the post-breaking stage for $\omega t\gt 35$ , when the wave amplitude is greatly reduced, the regions with negative $u-c$ reduce compared to the growing stage (see figue 4b,c).

As noted in (Veron et al. Reference Veron, Saxena and Misra2007), although flow reversal suggests airflow separation from the wave field, a more rigorous approach involves examining the spanwise vorticity $\Omega _y$ , which is a Galilean-invariant quantity. The spanwise vorticity normalized by the angular velocity, i.e. $\Omega _y/\omega$ , is reported in figure 4(d–f). For the different panels, a thin layer of vorticity is observed to form on the windward side of the wave, i.e. before the wave crest, where the turbulent airflow shears this layer, causing it to detach from the surface and mix with the flow. Lacking support from the wave field, the detached vorticity ‘packet’ destabilizes, breaking down before gradually reforming at the next wave crest. Although figure 4(d–f) refers to a specific region of the domain, i.e. $y/\lambda =0$ , flow separation occurs all over the wave crests, as shown in figure 4(g–i), where the spanwise-averaged vorticity is displayed. This observation suggests that, in the present configuration with steep waves under breaking conditions, flow separation occurs over all the wave crests, as also observed in Buckley et al. (Reference Buckley, Veron and Yousefi2020).

Flow separation is an important feature of wind-forced breaking waves and has been argued to enhance the momentum flux in breaking waves (Banner & Melville Reference Banner and Melville1976; Reul et al. Reference Reul, Branger and Giovanangeli1999; Buckley et al. Reference Buckley, Veron and Yousefi2020), and is analysed in the next section.

4.1. Momentum budget

The momentum flux, $\boldsymbol \tau _{t}=(\tau _{t,x},\tau _{t,y},\tau _{t,z})$ , is defined as the total force exerted from the wind to the waves per unit of air mass density and can be decomposed into a pressure and a viscous contribution (Belcher & Hunt Reference Belcher and Hunt1998; Sullivan et al. Reference Sullivan, Banner, Morison and Peirson2000). For water waves moving parallel to the wind in the streamwise $x$ direction, the first force component, i.e. $\tau _{t,x}$ , is typically dominant and reads

(4.1) \begin{equation} \tau _{t,x} = \tau _{p,x} + \tau _{\nu ,x} = - \int _\Gamma p\textbf {n}\cdot \textbf {e}_x{\rm d}S + 2\mu _a\int _\Gamma (\textbf {D}\textbf {n})\cdot \textbf {e}_x{\rm d}S\textrm {,} \end{equation}

where both the pressure $p$ and the strain rate tensor $\textbf {D}$ are evaluated in the air phase. The two contributions to the right-hand side of (4.1) are termed pressure (or form) drag and viscous drag.

Figure 5. Surface distribution of the pressure stress $p_an_x$ (a,c,e) and viscous stress $\tau _{sx}=2\mu _a(\textbf {D}\textbf {n})\cdot \textbf {e}_x$ (b,d,f) for $u_\ast /c=0.9$ . Note that both quantities are normalized by the total imposed stress $\rho _au_\ast ^2$ at the interface. The dot-dashed lines in all the panels represent the position of the wave crest.

We start by qualitatively inspecting both contributions during the growing and breaking stages. Figure 5 reports the instantaneous distribution of $p_an_x$ and $2\mu _a (\textbf {D}\textbf {n} )\cdot \textbf {e}_x$ for the case $u_\ast /c=0.9$ at three physical times $\omega t=[20, 35, 70]$ , corresponding to the growing pre-breaking stage, to the breaking stage and to the second growth post-breaking. Hereinafter, we define $\tau _{sx}=2\mu _a (\textbf {D}\textbf {n} )\cdot \textbf {e}_x$ for brevity and, therefore, $\tau _{\nu ,x}=\int _\Gamma \tau _{sx}dS$ in (4.1). The distributions $p_an_x$ and $\tau _{sx}$ are both projected on a wave-following surface, very close to the surface (vertically shifted by $0.1/k$ from the wave surface; see Wu et al. (Reference Wu, Popinet and Deike2022) and Appendix B for details).

Figure 5 shows that both $p_an_x$ and $\tau _{sx}$ display clear wave-induced coherent patterns, together with the development of three-dimensional structures induced by turbulence. When the wave field grows, i.e. figures 5(a,b) and 5(e,f), the pressure assumes negative values in the wave troughs and positive values near the wave crest. Similarly, the viscous stress also assumes positive value near the wave crest, whereas it approaches zero near the wave troughs and intermittently becomes negative at $\omega t=20\,{-}\,70$ , consistent with the measurements in Veron et al. (Reference Veron, Saxena and Misra2007). Furthermore, the peaks in pressure appear on the windward face (on the left of the dotted black lines at the wave crest), whereas the peak in the viscous stress is localized near the wave crest, mainly due to turbulence-induced straining of the shear layer. Note that the peaks of $pn_x$ are one order of magnitude larger than the peak in $\tau _{sx}$ . This result is consistent with the large initial slope of the waves for which the pressure force is the leading term in (4.1), as also shown in experimental works (Buckley et al. Reference Buckley, Veron and Yousefi2020). When the wave breaks, i.e. figures 5(c,d), the distribution of $pn_x$ changes dramatically with a loss of the wave-coherent pattern composed of high-pressure and low-pressure regions in the windward and leeward sides, respectively. Conversely, the viscous stress distribution is much less affected by breaking.

We can quantify how the abrupt change in the form drag distribution reflects in the total momentum budget. We integrate the streamwise component of the momentum equation (2.3) over the air volume $\Omega _a$ :

(4.2) \begin{equation} \underbrace {\rho _a\dfrac {\partial }{\partial t}\int _{\Omega _a} u\,{\rm d}V}_{\rho _a\partial U/\partial t} + \underbrace {\rho _a\int _{\Gamma }u(\textbf {u}_r\cdot \textbf {n})\,{\rm d}S}_{\rho _a\phi _{c,x}} = \tau _{p,x} + \tau _{\nu ,x} + \underbrace {\rho _au_\ast ^2A_\Gamma }_{\Pi _f}\textrm {,} \end{equation}

where ${\rm d}V$ is the elementary volume in the airflow region and $\Pi _f$ is the volume-integrated imposed pressure gradient (defined in (2.4)), i.e. $\Pi _f=\int _{\Omega _a}\partial p_0/\partial x (1-\mathcal {F}){\rm d}V=\rho _au_\ast ^2A_\Gamma$ with $A_\Gamma =\Omega _a/(L_0-h_W)$ . In (4.2), the left-hand side is the total rate of change in the air velocity, accounting for the temporal variation $\rho _a\partial U/\partial t$ and the convective contribution $\rho _a\phi _{c,x}$ . The former term, $\rho _a\partial U/\partial t$ , accounts for the response of the instantaneous flow field to the variation in the momentum total flux $\tau _{t,x}=\tau _{p,x}+\tau _{\nu ,x}$ , and the latter, $\rho _a\phi _{c,x}$ , accounts for the momentum flux originated from a non-zero relative velocity $\textbf {u}_r$ between the airflow and the wave field. Since the magnitude of $\textbf {u}_r$ is enhanced in the presence of airflow separation events, the term $\rho _a\phi _{c,x}$ can be employed to quantify such events over growing and breaking waves. The right-hand side includes the uniform forcing term $\Pi _f$ and the two forces per unit of mass $\tau _{p,x}$ and $\tau _{\nu ,x}$ , as in (4.1). Note that in a statistically steady condition with negligible relative velocity between the airflow and the wave field, the two terms on the left-hand side of (4.2) vanish, i.e. $\langle \rho _a\partial U/\partial t\rangle _t=\langle \rho _a\phi _{c,x}\rangle _t=0$ (with $\langle \rangle _t$ a time-averaged operator). Accordingly, equation (4.2) reduces to the simplified form where the total imposed stress at the interface $\rho _au_\ast ^2A_\Gamma$ balances the sum of pressure and viscous drag (Janssen Reference Janssen2004; Buckley et al. Reference Buckley, Veron and Yousefi2020; Funke et al. Reference Funke, Buckley, Schultze, Veron, Timmermans and Carpenter2021). More details on the evaluation of the momentum fluxes are provided in Appendix B.

Figure 6. Contributions of the momentum budget in the streamwise direction, as in (4.2), for (a) $u_\ast /c=0.5$ and (b) $u_\ast /c=0.9$ . On the $y$ -axis label $\mathcal {T}$ represents the variation in the instantaneous flow $\rho _a\partial U/\partial t$ , the viscous stress $\tau _{\nu ,x}$ , the pressure stress $\tau _{p,x}$ , the convective term $\rho _a\phi _{c,x}$ or the driving force $\Pi _f$ (defined in (2.4)). Each budget component is normalized by the total stress $\rho _au_\ast ^2A_\Gamma$ . For both cases, the normalized variation in the wave energy $E_W/E_{W,0}$ is reported in the top panel.

We analyse the temporal variation of the different terms in the momentum flux (4.2) for two representative cases ( $u_\ast /c=0.5\,{-}\,0.9$ ) in figure 6. The top panels display the normalized wave energy variation $E_W(t)$ to remind the reader of the overall time evolution. In both cases, the wave field experiences a first growth up to the breaking event ( $\omega t\in [0-65]$ for $u_\ast /c=0.5$ and $[0-20]$ for $u_\ast /c=0.9$ , the shorter growing cycle for $u_\ast /c=0.9$ being due to the higher wind forcing).

When the wave breaks, the momentum flux due to pressure $\tau _{p,x}/(\rho _au_\ast ^2A_\Gamma )$ drops with a corresponding acceleration of the flow, i.e. $\rho _a\partial U/\partial t/(\rho _au_\ast ^2A_\Gamma )$ increases, the drop being larger for the strongest wind forcing. Once the breaking stage ends, around $\omega t=75$ for $u_\ast /c=0.5$ and around $\omega t=35$ for $u_\ast /c=0.9$ , the wave field experiences a second growing cycle. During the interval, the pressure force remains the dominant contributor to the momentum budget and keeps increasing with a progressive deceleration of the airflow. For the case $u_\ast /c=0.9$ , the wind forcing drives the wave field to a second breaking event, which occurs around $\omega t=120$ . Similarly to the first breaking event, the pressure force decreases while the airflow accelerates, with the magnitude of the momentum loss being smaller (since the breaking is weaker).

In all cases, the momentum flux originating from the non-zero velocity between the airflow and the waves, $\rho _a\phi _{c,x}$ , represents a negligible part of the momentum budget. This is consistent with Yang et al. (Reference Yang, Deng and Shen2018), who found a negligible role of the convective term in the budget. Next, the viscous contribution $\tau _{\nu ,x}$ represents a small but not negligible contribution in the forces $\tau _{\nu ,x}$ . This is consistent with the high initial wave steepness, i.e. $a_0k=0.3$ . Notably, upon wave breaking, the viscous contribution experiences a slight increase to partially compensate for the loss of the pressure force.

The sensitivity of the pressure force to the breaking event can be qualitatively understood by considering the variations of the instantaneous slope. When the wave breaks, the instantaneous wave slope, $\partial \eta /\partial x$ , suddenly reduces and directly influences the dominant pressure stress $\tau _{p,x}$ (which is directly visible in a simplified version of (4.1); see Funke et al. (Reference Funke, Buckley, Schultze, Veron, Timmermans and Carpenter2021)).

4.2. Velocity profiles

In § 4.1, we showed that when the wave breaks, the momentum flux associated with pressure is reduced in favour of a flow acceleration in the air. We now characterize the associated streamwise vertical velocity profiles $\langle u_a\rangle$ during the growing–breaking cycle.

Following Sullivan et al. (Reference Sullivan, McWilliams and Moeng2000) and Wu et al. (Reference Wu, Popinet and Deike2022), the mean velocity profile is computed in wave-following coordinates by using an implicit transformation, which maps the Cartesian coordinates $(x,y,z)$ to wave-following coordinates $(\xi ,\eta ,\zeta )$ . This step allows us to perform the space average of the velocity field along the periodic directions by including the region below the wave crests. More details about the procedure to transform a generic field defined in a Cartesian coordinate system into a wave-following one are given in Appendix A.

Figure 7. Streamwise velocity profile normalized by the nominal friction, $\langle u_a\rangle /u_\ast$ , as a function of vertical wave-following coordinate $\zeta /\lambda$ in the airflow for $(a)$ $u_\ast /c=0.5$ and $(b)$ $u_\ast /c=0.9$ . For large enough values, $\zeta =z$ . The instantaneous velocity profiles are averaged in time over the cycle $G_1$ (dark-blue lines) and a fraction of the second growing cycle, $G_{2,a}$ , as defined in § 3.2 (see figure 3 b). The dot-dashed curves represent the instantaneous values during the breaking stage, i.e. $\omega t\in [58\,{-}\,98]$ for $u_\ast /c=0.5$ and $\omega t\in [22\,{-}\,42]$ for $u_\ast /c=0.9$ .

Figure 7 shows the streamwise velocity profile spatially averaged along the horizontal directions as a function of the vertical wave-following coordinate $\zeta$ for $u_\ast /c=0.5\,{-}\,0.9$ . In both cases, the velocity profiles are time-averaged over three time windows. The two windows correspond to the two growing cycles, $G_1$ and $G_{2,a}$ , introduced in figure 3(b). Furthermore, note that the velocity profiles in the airflow are not subtracted by the streamwise component of the surface water velocity, as its maximum value in the water throughout the simulation is of the order of $0.33u_\ast$ , which is negligible compared with the maximum streamwise velocity in the airflow (see figure 7).

By comparing the mean velocity profile $\langle u_a\rangle$ during the stages $G_1$ and $G_{2,a}$ , we clearly see that the breaking event leads to a flow acceleration mainly confined in the region $\zeta /\lambda \lt 1$ and a shift of the velocity profiles to larger values, as also clearly shown from the instantaneous profiles during the breaking stage. The flow acceleration becomes more pronounced as we compare the case $u_\ast /c=0.5$ with the case $u_\ast /c=0.9$ . The profile is non-zero at the water surface due to the presence of a developing underwater current.

Note that even during the growing stage, the simulation remains transient in nature. However, the turbulent flow adjusts to the moving wave field on a time scale $t_t=\nu _a/u_\ast ^2$ that is $40$ to $200$ times smaller than the wave period $T_0 = \lambda /c$ for the different analysed cases, i.e. $T_0/t_t = (u_\ast /c)Re_{\ast ,\lambda } \gt 40\,{-}\,200$ . Therefore, except during the breaking stage, the configuration can be approximated as quasi-stationary, allowing us to analyse the modulation of the velocity profile by the growing waves and the breaking event.

This analysis demonstrates the specific contribution of breaking events to the average velocity profiles, which, in laboratory or open-ocean conditions, are typically averaged over long time periods without distinguishing between breaking and growing stages.

4.3. Reynolds shear stress

The growing and breaking cycle also affects the Reynolds shear stress $\langle u'w'\rangle$ (or turbulent intensity) in the near-wave-field region. Figure 8 displays the Reynolds stress spatially averaged over the two periodic directions and normalized by $u_\ast ^2$ , i.e. $-\langle u'w'\rangle /u_\ast ^2$ , for $u_\ast /c=0.5\,{-}\,0.9$ (using the wave-following coordinate as detailed in Appendix A). Each instantaneous profile is also time-averaged during the first and a fraction of the second growing cycle, i.e. $G_1$ and $G_{2,a}$ . Note that we also report the different instantaneous values during the breaking stage of each case.

Figure 8. Reynolds shear stress normalized by the square of the nominal friction, $-\langle u'w'\rangle /u_\ast ^2$ , as a function of the vertical wave-following coordinate $\zeta$ for (a) $u_\ast /c=0.5$ and (b) $u_\ast /c=0.9$ . The Reynolds shear stress is averaged over the same time windows as in figure 7. The dot-dashed curves represent the instantaneous values during the breaking stage, i.e. $\omega t\in [58\,{-}\,98]$ for $u_\ast /c=0.5$ and $\omega t\in [22\,{-}\,42]$ for $u_\ast /c=0.9$ . The dashed green curve represents the Reynolds stress on a flat stationary surface at $Re_\ast =720$ .

For both cases, $-\langle u'w'\rangle /u_\ast ^2$ shows an increasing trend in the near-wave region, reaching a peak at $\zeta /\lambda \approx 0.26$ (i.e. $z/\lambda \approx 0.38$ ). Beyond this point, it decreases and approaches the dashed green curve in figure 8, which represents the Reynolds stress for a flat and stationary wall at $Re_\ast =720$ (Pope Reference Pope2000). However, it does not reach this solution exactly due to the non-stationarity of the flow field induced by the moving waves, i.e. growing and breaking stages.

During the growing stage, the peak in the Reynolds stress is $0.62$ for $u_\ast /c=0.5$ , whereas it is $0.76$ for $u_\ast /c=0.9$ . This larger peak is due to the increased drag, primarily from the pressure component, at higher $u_\ast /c$ . When the wave field breaks, as shown from the dot-dashed curves, both peaks decrease without a significant change in the peak location, with a more pronounced reduction at larger $u_\ast /c$ related to the stronger breaking event.

Thus, wind-induced wave breaking greatly modulates the Reynolds stress and confirms the discussion made when analysing the pressure field. During the growing stage, the turbulent drag, i.e. $-\langle u'w'\rangle /u_\ast ^2$ , is larger as $u_\ast /c$ increases, while when the wave breaks, the turbulent drag is reduced, and such loss increases with $u_\ast /c$ . We note that both cases are run at the same $Re_{\ast ,\lambda }$ , so the phenomenology is indeed controlled by the wind forcing and breaking dynamics and is much less affected by the Reynolds number.

4.4. Extraction of the surface roughness from the velocity profiles

The breaking-induced shift can be quantified considering the velocity profiles in semi-logarithmic form and rescaled in wall units, as shown in figure 9. The velocity profiles display a consistent upshift which extends for the entire inner layer (Pope Reference Pope2000; Cimarelli et al. Reference Cimarelli, Romoli and Stalio2023), composed of the viscous sublayer ( $0\lt \zeta ^+\lt 5$ ), the buffer layer ( $5\lt \zeta ^+\lt 30$ ) and log-law region ( $\zeta ^+\gt 30$ ), where $\zeta ^+$ is the vertical coordinate in wall units $\zeta ^+=\zeta u_\ast /\nu _a$ with $\nu _a=\mu _a/\rho _a$ . The upshift is associated with the change in the instantaneous steepness due to breaking, leading to an overall drag reduction in the flow. Note that this upshift in the logarithmic region is mainly confined near the wave field ( $30\lt \zeta ^+\lesssim 100$ ). For larger $\zeta ^+$ , the mean velocity profile in the post-breaking stage, i.e. $G_{2,a}$ , approaches that in the pre-breaking stage, i.e. $G_1$ , confirming that the airflow modulation induced by wave breaking is negligible in this region.

Figure 9. Streamwise velocity profile normalized by the friction velocity, $\langle u_a^+\rangle =\langle u_a\rangle /u_\ast$ , as a function of vertical wave-following coordinate (in wall units) $\zeta ^+=\zeta u_\ast /\nu _a$ with $\nu _a=\mu _a/\rho _a$ for (a) $u_\ast /c=0.5$ and (b) $u_\ast /c=0.9$ , both at $Re_{\ast ,\lambda }=214$ ( $Re_\ast =720$ ). The velocity profiles are averaged over the same time windows as in figure 7. The dotted black lines refer to the fitted logarithmic law employed to estimate the intercept for each case. The continuous black line represents the mean velocity profile at $Re_\ast = 720$ for a flat stationary surface.

It is also worth mentioning that the streamwise velocity profile at the same $Re_\ast$ for a flat stationary wall is well above that with moving waves. Indeed, flows over waves experience much larger drag due to the added contribution of the pressure component, which is zero for flows over a flat wall (Belcher & Hunt Reference Belcher and Hunt1998). This observation agrees with experimental works where a downshift of the velocity profile was observed as the wave steepness was increased (Buckley & Veron Reference Buckley and Veron2016; Buckley et al. Reference Buckley, Veron and Yousefi2020).

Using the velocity profile in logarithmic form, we can quantify this upshift induced by the wave breaking. For this purpose, we express $\langle u_a\rangle$ using the wall-normal coordinate:

(4.3) \begin{equation} \dfrac {\langle u_{a,f}\rangle - \langle u_{a,i}\rangle }{u_\ast } = \dfrac {1}{\kappa }\log \left (\dfrac {\zeta _f^+-\zeta _i^+}{z_{0}^+}\right )\textrm {,} \end{equation}

where $\kappa$ is the Von Kármán constant, $z_{0}^+$ is the surface roughness and $\zeta _{i,f}^+$ are the initial and final extents of the log-law region, i.e. $\zeta _i^+$ and $\zeta _f^+$ . The corresponding velocities at these positions are termed $\langle u_{a,i}\rangle$ and $\langle u_{a,f}\rangle$ , respectively. Using the velocity profile in the logarithmic region, both $\kappa$ and $z_0^+$ can be estimated using a least-squares fit procedure (Sullivan et al. Reference Sullivan, McWilliams and Moeng2000; Lin et al. Reference Lin, Moeng, Tsai, Sullivan and Belcher2008) applied to the growing and breaking regimes. To evaluate $z_{0}^+$ and $\kappa$ , we consider $\zeta _i^+=30$ for all the stages, which corresponds, in physical units, to $\zeta _i=0.56\lambda$ , $0.28\lambda$ and $0.14\lambda$ for $Re_{\ast ,\lambda }=53.5$ , $107$ and $214$ , respectively. The final extent $\zeta _f^+$ is taken equal to a common value of $\zeta _f^+\approx Re_{\ast ,\lambda }/2$ to capture the breaking-induced drag reduction, which is confined in the near-wave region, as clearly shown in figure 9. Note that such $\zeta _f^+$ corresponds to a physical reference height of $\zeta _f=\lambda /2$ , similar to the reference height employed to compute the wind-wave growth in Jeffreys (Reference Jeffreys1925) and Donelan et al. (Reference Donelan, Babanin, Young and Banner2006) and to the height to compute the aerodynamic drag coefficient in § 5.1.

Results are reported in table 2. We note that the surface roughness slightly increases with $u_\ast /c$ , whereas when the wave field breaks, its value drops and decreases with $u_\ast /c$ . This trend is also confirmed during the second breaking cycle for the cases where it is available, although the change in $z_0^+$ during breaking is smaller due to the reduced breaking strength. Note that $z_0^+$ immediately after the breaking stage, i.e. over $G_{2,a}$ , is smaller than $z_0^+$ evaluated before the second breaking, i.e. $G_{2,b}$ (when available), as the wave field has grown between the two windows. It is worth remarking that irrespective of the cases, we found a common value for the Von Kármán constant $\kappa \approx 0.40$ , in agreement with the employed value in ocean and atmosphere models.

Table 2. Surface roughness for the stages $G_1$ , $G_{2,a}$ , $G_{2,b}$ and $F$ , as defined in § 3.2 (see figure 3 b). The surface roughness is expressed in viscous (or plus) units, $z_0^+=z_0u_\ast /\nu _a$ with $\nu _a=\mu _a/\rho _a$ . The values are extracted from the velocity profiles in logarithmic form using a best-fit procedure reported in figure 9. For the case $u_\ast /c=0.3$ , only one surface roughness value is reported since the wave field, in this case, is in equilibrium with the flow. For the cases $u_\ast /c=[0.7\,{-}\,0.9]$ at $Re_{\ast ,\lambda }=214$ and for the case $u_\ast /c=0.9$ at $Re_{\ast ,\lambda }=107$ with $(L_0-h_W)/\lambda =6.72$ , $z_0^+$ is also reported for the second breaking cycle.

We note that the drag reduction mentioned here is studied at a fixed $Re_{\ast ,\lambda }$ and with $(L_0-h_w)/\lambda =3.36$ . To assess the sensitivity of the results to these parameters, we perform three additional simulations at $u_\ast /c=0.9$ , two at $Re_{\ast ,\lambda }=53.5\,{-}\,107$ with $(L_0-h_w)/\lambda =3.36$ and one at $Re_{\ast ,\lambda }=107$ with $(L_0-h_w)/\lambda =6.72$ . Results are discussed in Appendices D and E, and the surface roughness of these cases is reported in table 2. We observe a breaking-induced drag reduction of comparable magnitude to that measured at the largest $Re_{\ast ,\lambda }$ . This supports that the drag modulation induced by wave breaking is primarily controlled by $u_\ast /c$ and is less sensitive to $Re_{\ast ,\lambda }$ and $(L_0-h_W)/\lambda$ .

5.1. Estimation of the aerodynamic drag coefficient

In § 4.1, we showed that the pressure drag force $\tau _{p,x}$ increases in time as long as the wave field grows (associated with increased wave slope) and suddenly decreases when the wave field breaks. In § 4.2, we observed that the breaking event is associated with a flow acceleration, i.e. $\langle u_a\rangle /u_\ast$ shifts towards larger values, and airflow separation occurs over the wave crests and troughs. To discuss the relative strength of the variation of the pressure force with the relative strength of mean flow variation, we introduce a dimensionless aerodynamic drag coefficient $C_{D,a}$ :

(5.1) \begin{equation} C_{D,a} = \dfrac {2\overline {\tau }_{p,x}}{\rho _a\overline {U}_{ref}^2(\overline {\zeta }=\zeta _{ref})A_\Gamma }, \end{equation}

where $\overline {\tau }_{p,x}$ is the time-averaged mean momentum flux associated with pressure:

(5.2) \begin{equation} \overline {\tau }_{p,x} = \dfrac {1}{\triangle T_{G,B}}\int _{\triangle T_{G,B}} \tau _{p,x}\, {\rm d}t\textrm {,} \end{equation}

where $\triangle T_{G,B}$ is the time interval when the wave grows or breaks, respectively, and the integral is taken over this time interval. An equivalent definition to (5.2) is employed for $\overline {U}_{ref}$ . We split the dynamics into the growing and breaking stages to understand their respective contribution to the total aerodynamic drag, with the same convention as before, so that the length of the growth and breaking intervals varies with $u_\ast /c$ .

The choice of the reference far-field boundary layer velocity $\overline {U}_{ref}$ evaluated at a reference height $\overline {\zeta }=\zeta _{ref}$ is based on two requirements: (i) being sufficiently far from the wave surface to always reside in the air during growth and breaking and (ii) being sufficiently close to the wave field to be affected by its dynamics (growth and breaking). It follows that $\overline {U}_{ref}$ increases as the flow accelerates in the breaking stage. Here, the requirements are satisfied typically if $\lambda /4 \lt \zeta _{ref} \lt \lambda$ , and we consider $\zeta = \lambda /2$ , similar to wind-wave growth sheltering discussion (Jeffreys Reference Jeffreys1925; Donelan et al. Reference Donelan, Babanin, Young and Banner2006). Different $\zeta _{ref}$ within the interval only changes the magnitude $C_{D,a}$ , without altering the discussion.

Figure 10. Aerodynamic drag coefficient $C_{D,a}$ (defined by (5.1)) for different $u_\ast /c$ in the growing (blue colours) and breaking (red colours) time intervals. For $u_\ast /c\lt 0.35$ , the simulated wave field is only growing (one would have to run the simulations longer to obtain breaking), so that all data are in the growing regime (and include increasing $a_0k$ ; see Wu et al. Reference Wu, Popinet and Deike2022). For $u_\ast /c\gt 0.35$ , both growing and breaking are present, and both ranges are separated by the red dash-dotted line. Whenever available, the data pertaining to the second growing and breaking cycles $G_2$ , $B_2$ are displayed. The growing and breaking stages $G_1$ , $G_2$ , $B_1$ and $B_2$ are defined in figure 3(b). Growing dynamics displays a systematic increase in drag with $u_\ast /c$ , while breaking induces a decrease in drag with increasing $u_\ast /c$ . The averages of the breaking and growing cycles are indicated in magenta, and we observe a saturation of the averaged drag at high wind speed.

Figure 10 shows the estimated $C_{D,a}$ as a function of $u_\ast /c$ and for both growing and breaking stages. We also re-analyse the non-breaking growing simulations from Wu et al. (Reference Wu, Popinet and Deike2022) at lower $u_\ast /c$ and varying initial slope. In the non-breaking growing regime ( $u_\ast /c \lt 0.35$ ), $C_{D,a}$ increases with $u_\ast /c$ , with different values of $C_{D,a}$ at the same $u_\ast /c$ corresponding to different initial wave steepness (see also Wu et al. Reference Wu, Popinet and Deike2022).

For the range of $u_\ast /c$ that includes growth and breaking ( $u_\ast /c \gt 0.35$ ), $C_{D,a}$ follows an increasing trend with $u_\ast /c$ in the growing (pre-breaking) stage (blue symbols) and a decreasing trend in the breaking stage (red symbols). Averaging $C_{D,a}$ over both the growing and breaking stages, the drag coefficient exhibits a saturation with a slight decrease at the highest wind speed (magenta symbols). Similar results are observed during the second breaking cycle, as shown in figure 10 for the available cases ( $u_\ast /c=[0.7\,{-}\,0.9]$ at $Re_{\ast ,\lambda }=214$ and $u_\ast /c=0.9$ at $Re_{\ast ,\lambda }=107$ ) with a slightly smaller drag coefficient, compared with the first breaking cycle. Thus, we can expect that averaging over multiple growing–breaking cycles would lead to a similar overall behaviour. Finally, the results are not very sensitive to the friction Reynolds number $Re_{\ast ,\lambda }$ and the number of waves in the simulation box $(L_0-h_W)/\lambda$ .

From the analysis presented in this section, we can conclude that the aerodynamic drag saturation is driven by breaking events, which cause a marked reduction in the pressure drag, airflow separation and a loss of coherence between the wave and the pressure fields.

5.2. Comparison of the drag coefficient with laboratory and field measurements

The above discussion suggests that in laboratory and field measurement analysis, the drag saturation at high wind speed is the consequence of the spatial and temporal averages over breaking events and wave growth. In this section, we can qualitatively compare our results, obtained for narrowbanded idealized conditions, with the more complex multiscale observations at high wind speed (Donelan et al. Reference Donelan, Haus, Reul, Plant, Stiassnie, Graber, Brown and Saltzman2004; Janssen Reference Janssen2004; Bye & Jenkins Reference Bye and Jenkins2006; Buckley et al. Reference Buckley, Veron and Yousefi2020; Curcic & Haus Reference Curcic and Haus2020; Sroka & Emanuel Reference Sroka and Emanuel2021), by considering the classic oceanographic definition of the drag coefficient at $10\rm\, m$ height, termed $C_D$ . As mentioned in § 1, equation (1.1), $C_D$ is related to the definition of the total momentum flux $\tau _{t,x} = \rho _aC_D U_{10}^2$ with $\tau _{t,x}=\rho _au_\ast ^2$ since at $10\rm \, m$ height the airflow and the waves are assumed to be in equilibrium. For neutral atmospheric conditions, the velocity profile at $\overline {z}$ can be assumed to follow the logarithmic law of the wall (Edson et al. Reference Edson, Jampana, Weller, Bigorre, Plueddemann, Fairall, Miller, Mahrt, Vickers and Hersbach2013; Ayet & Chapron Reference Ayet and Chapron2022), i.e. $U_{10} = (u_\ast /\kappa )\log (\overline {z}/z_0 )$ as in (4.3), and $C_D$ becomes

(5.3) \begin{equation} C_D(\overline {z}=10\,\ \textrm {m}) = \dfrac {\kappa ^2}{\log ^2\left (\overline {z}/z_0\right )}\textrm {.} \end{equation}

Using the DNS data generated for this work and those reported in Wu et al. (Reference Wu, Popinet and Deike2022), we compute $C_D$ from (5.3) by using the velocity profiles evaluated during the first and a fraction of the second growing cycles for the different $u_\ast /c$ and $Re_{\ast ,\lambda }$ (and in the second growing–breaking cycle when available). We convert the velocity profile in physical units by considering the length and time scale of the numerical system defined based on the air–water properties and the wave scales given by the Bond number. Using a best-fit procedure on the equation of the logarithmic profile ( $\kappa =0.40$ , as discussed in § 4.2), a dimensional friction velocity and roughness $z_0$ can be retrieved, and $C_D$ can be evaluated using (5.3). Note that we checked a posteriori that this analysis conserves the ratio of the friction velocity over the wave phase speed. The calculation of $C_D$ is done on two averaging windows representative of the flow during the first growing cycle ( $G_1$ ) and the disturbed flow just after breaking ( $G_{2,a}$ ), since defining a logarithmic velocity profile requires a quasi-stationary profile. Such an approach differs slightly from the one employed for $C_{D,a}$ , but we will see that the physical conclusions are identical, while the 10-m-based drag coefficient allows us to compare with existing laboratory data.

The values of $C_D$ for different $u_\ast /c$ exhibit a similar qualitative trend to that of $C_{D,a}$ with $u_\ast /c$ . We observe a systematic increase of the drag in the growing stage (dark triangles and blue symbols) with $u_\ast /c$ , while the drag in the breaking stage (red symbols) decreases with $u_\ast /c$ . When averaging over both growing and breaking regimes, we observe a saturation at high wind speed ( $u_\ast /c\gt 0.4$ ). As for the aerodynamic drag coefficient $C_{D,a}$ , the results are not very sensitive to $Re_{\ast ,\lambda }$ and $(L_0-h_W)/\lambda$ , as discussed in Appendix D and Appendix E.

It is remarkable that the numerical values of $C_D$ in figure 11, with $C_D$ between 0.2 and $0.8 \times 10^{-3}$ for $u_\ast /c\lt 0.3$ and $C_D$ between 1.5 and $2.0\times 10^{-3}$ for $u_\ast /c\gt 0.4$ , closely match those reported in laboratory experiments (Buckley et al. Reference Buckley, Veron and Yousefi2020; Curcic & Haus Reference Curcic and Haus2020) for similar values of $u_\ast /c$ (shown in green symbols). In laboratory experiments at high wind speeds, the drag coefficient is determined by long-time averaging over multiple growing and breaking cycles, leading to drag saturation.

We note that the saturation of $C_D$ occurs at lower $u_\ast /c$ and $a_{{rms}}k$ in our simulations than in laboratory experiments. This discrepancy may be attributed to differences in set-up: our forced narrowbanded wave field contrasts with the laboratory’s multiscale finite fetch conditions at high wind speeds, where averaging spans multiple breaking–growing cycles, and the partitioning between regimes at the same $u_\ast /c$ may differ from our simulations. Similarly, directly comparing realistic $u_\ast /c$ values in field conditions is challenging due to the broad-banded nature of the wave field in the open ocean at high-wind-speed conditions. Additionally, the lower wave Reynolds number in our simulations affects the effective growth rate and time to reach breaking conditions and makes direct comparison in physical units, e.g. actual wind speed and wavelength, more challenging. Nevertheless, the ability of the simulations to capture drag saturation at high wind speed suggests that we are indeed approaching the correct asymptotic limit for high Reynolds numbers.

Figure 11. Drag coefficient $C_D$ evaluated at $\overline {z}=10$ $\textrm {m}$ using (5.3) as a function of $u_\ast /c$ (a) and $a_{{rms}}k$ (b) with $a_{{rms}}=a/\sqrt {2}$ . The figure includes the calculation of $C_D$ with the data generated in this work and those retrieved from Wu et al. (Reference Wu, Popinet and Deike2022). For all the cases, the reported $C_D$ is an average value between the first growing cycle, $G_1$ , and a fraction of the second growing cycle, $G_{2,a}$ (immediately after the breaking event). Whenever available, the data pertaining to the second growing cycle $G_{2,b}$ and the final stage $F$ are displayed. The employed time window to define $C_D$ follows the convention given in figure 3(b). For the cases at $u_\ast /c=[0.4\,{-}\,0.5\,{-}\,0.7\,{-}\,0.9]$ at $Re_{\ast ,\lambda }=214$ with $(L_0-h_W)/\lambda =3.36$ , we separate between these two stages (blue and red circles). The green symbols display the experimental datasets from Buckley et al. (Reference Buckley, Veron and Yousefi2020) up to $u_\ast /c\approx 0.71$ and from Curcic & Haus (Reference Curcic and Haus2020) up to $u_\ast /c\approx 2.25$ .

The present analysis can be related to the discussion from Sullivan et al. (Reference Sullivan, Banner, Morison and Peirson2018) and Wu et al. (Reference Wu, Popinet and Deike2022) on the control of the form drag by the wave slope, $a_{{rms}}k$ , for steep non-breaking waves. In the present work, we showed that the breaking event, which is responsible for wave slope saturation controls the saturation of the drag coefficient when presented as a function of $u_\ast /c$ .

Moreover, the drag coefficient $C_D$ is shown as a function of $a_{{rms}}k$ in figure 11. During the growth phase of the wave field, $C_D$ increases with $a_{{rms}}k$ up to approximately $a_{{rms}}k \approx 0.15$ for non-breaking waves with finite steepness (Wu et al. Reference Wu, Popinet and Deike2022). As $a_{{rms}}k$ reaches larger values, wave breaking initiates around $a_{{rms}}k \approx [0.20\,{-}\,0.23]$ . Focusing on the growing phase of the wave field, $C_D$ increases with $a_{{rms}}k$ for the different $u_\ast /c$ ; however, as we approach the largest values, $a_{{rms}}$ saturates, leading to a corresponding saturation in $C_D$ . Conversely, during the breaking stage, $a_{{rms}}k$ decreases as $u_\ast /c$ increases, resulting in a reduction in $C_D$ . When considering the mean value of $C_D$ over a complete growing and breaking cycle, the variations in both $C_D$ and $a_{{rms}}$ become significantly smaller, converging to similar values, so that all data points in the saturated drag coefficient also have saturated slope and appear clustered on the graph. This behaviour highlights that wave breaking dynamics not only governs the saturation of $C_D$ but also controls the saturation of the wave slope. This discussion presents similarities to Davis et al. (Reference Davis, Thomson, Houghton, Doyle, Komaromi, Fairall, Thompson and Moskaitis2023) showing saturation of the wave mean square slope during tropical cyclones observed from drifting buoys, suggesting that simple parameterizations of the drag coefficient based on a wave-slope metric could be promising.

We suggest that experimental and field analysis of the drag coefficient reporting a statistical analysis of the partitioning between breaking and growing (temporally and/or spatially) could help better understand, quantify and eventually parameterize drag saturation at high wind speed and constrain momentum flux models leveraging breaking occurrence statistics (Kudryavtsev et al. Reference Kudryavtsev, Chapron and Makin2014).

Finally, it is very important to remark that while the saturation and reduction of $C_D$ has sometimes been attributed to the production of sea sprays (Bye & Jenkins Reference Bye and Jenkins2006; Veron Reference Veron2015), in the current set-up, the production of droplets is negligible during the breaking stage of the wave field. Yet, we observe that the drag coefficient, whether defined as $C_{D,a}$ or $C_D$ , displays a saturation at high $u_\ast /c$ when wave breaking is considered. This aspect shows that wave-breaking events, which cause the marked reduction of the pressure stress and the associated flow separation and determine the loss of coherence between the pressure and the wave fields, are sufficient to saturate the drag coefficient at high $u_\ast /c$ .