# Momentum fluxes in wind-forced breaking waves

Nicolò Scapin, Jiarong Wu, J. Thomas Farrar, Bertrand Chapron, Stéphane Popinet, Luc Deike

► To cite this version:

Nicolò Scapin, Jiarong Wu, J. Thomas Farrar, Bertrand Chapron, Stéphane Popinet, et al. Momentum fluxes in wind-forced breaking waves. Journal of Fluid Mechanics, 2025, 1009, A20 (33p.). 10.1017/jfm.2025.77. hal-05105692

HAL Id: hal-05105692

https://hal.science/hal-05105692v1

Submitted on 11 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Momentum fluxes in wind-forced breaking waves

Nicolò Scapin<fnref n="1" />,<fnref n="2" /> ID, Jiarong Wu<fnref n="1" />,<fnref n="3" /> ID, J. Thomas Farrar<fnref n="4" /> ID,
Bertrand Chapron<fnref n="5" /> ID, Stéphane Popinet<fnref n="6" /> ID and Luc Deike<fnref n="1" />,<fnref n="2" /> ID

**Corresponding author:** Luc Deike, ldeike@princeton.edu

(Received 23 September 2024; revised 27 December 2024; accepted 8 January 2025)

We investigate the momentum fluxes between a turbulent air boundary layer and a growing-breaking wave field by solving the air-water two-phase Navier-Stokes equations through direct numerical simulations. A fully developed turbulent airflow drives the growth of a narrow-banded wave field, whose amplitude increases until reaching breaking conditions. The breaking events result in a loss of wave energy, transferred to the water column, followed by renewed growth under wind forcing. We revisit the momentum flux analysis in a high-wind-speed regime, characterized by the ratio of the friction velocity to wave speed $u_*/c$ in the range [0.3 – 0.9], through the lens of growing-breaking cycles. The total momentum flux across the interface is dominated by pressure, which increases with $u_*/c$ during growth and reduces sharply during breaking. Drag reduction during breaking is linked to airflow separation, a sudden acceleration of the flow, an upward shift of the mean streamwise velocity profile and a reduction in Reynolds shear stress. We characterize the reduction of pressure stress and flow acceleration through an aerodynamic drag coefficient by splitting the analysis between growing and breaking stages, treating them as separate subprocesses. While drag increases with $u_*/c$ during growth, it decreases during breaking. Averaging over both stages leads to a saturation of the drag coefficient at high $u_*/c$, comparable to what is observed at high wind speeds in laboratory and field conditions. Our analysis suggests that this saturation is controlled by breaking dynamics.

**Key words:** air/sea interactions, gas/liquid flows

<footnote n="1">¹Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ 08544, USA,</footnote>

<footnote n="2">²High Meadows Environmental Institute, Princeton University, Princeton, NJ 08544, USA,</footnote>

<footnote n="3">³Courant Institute of Mathematical Sciences, New York University, New York, NY, USA</footnote>

<footnote n="4">⁴Woods Hole Oceanographic Institution (WHOI), Woods Hole, MA, USA</footnote>

<footnote n="5">⁵IFREMER, Univ. Brest, CNRS, IRD, Laboratoire d'Océanographie Physique et Spatiale (LOPS), Plouzané, France</footnote>

<footnote n="6">⁶Institut Jean Le Rond d'Alembert, CNRS UMR 7190, Sorbonne Université, Paris 75005, France</footnote>
# 1. Introduction

Ocean waves modulate the exchanges of mass, momentum and energy at the ocean-atmosphere interface. Waves continuously modify the ocean surface roughness and affect the momentum transfer between the atmospheric and oceanic boundary layers (Sullivan & McWilliams 2010). When waves break, sea spray droplets are produced (Veron 2015), together with the generation of bubbles due to air entrainment (Deike & Melville 2018; Deike 2022), enhancing heat and mass transfer (Veron 2015; Deike 2022). Wind forcing and ocean waves mutually influence each other. On the air side, the irregular and time-varying surface wave topography alters the wind velocity, air temperature and humidity in space and time. On the water side, the wind drives wave formation, promotes growth and steepening up to breaking and drives upper-ocean turbulence.

Wind-induced wave breaking plays a crucial role in the interaction between the ocean and the atmosphere. Wave breaking limits wave steepness, influences the momentum exchange between the atmospheric boundary layer and the upper ocean and locally alters wind velocity profiles and those of any other tracers (Melville & Rapp 1985; Melville 1996). Wave breaking is a highly dissipative process that controls the energy transfer from the wind to water currents and the transition to turbulence in the upper ocean (Lamarre & Melville 1991; Melville 1996; Veron & Melville 2001) and influences the wave-induced Langmuir turbulence (McWilliams et al. 1997; Melville et al. 1998; McWilliams 2016).

Understanding the momentum and the energy exchanges in wind-forced breaking waves remains an active area of research, especially in the high-wind-speed regime above a wind speed of 20–25 m s⁻¹ (evaluated at 10 m height). Accurate evaluation of the momentum fluxes (or wind stress) between the wind and the wave field is necessary to properly represent the turbulent boundary layers in the lower atmosphere and in the upper ocean. Without waves, the momentum flux at the ocean-atmosphere interface would be solely due to viscous effects. Waves introduce a pressure component in the flux, whose significance in the stress partitioning increases as wind speed and local wave steepness increase (Edson et al. 2013), so that open ocean measurements report an increase in the momentum flux measured at the top of the wave boundary layer (conventionally taken at a 10 m height), compared to a flat wall (Edson et al. 2013; Ayet & Chapron 2022).

In coupled ocean-atmosphere numerical models, the momentum flux $\tau$ is related to the wind speed at a given height using the drag coefficient, $C_D$:

$$ \tau = \rho_a u_*^2 = \rho_a C_D U_{10}^2, \quad (1.1) $$

where $u_*$ is the friction velocity, $\rho_a$ the air density and the reference velocity $U_{10}$ taken at a conventional reference height of 10 m (Edson et al. 2013). Assuming a logarithmic wind velocity profile (Monin & Obukhov 1954; Janssen 2004) in (1.1), the calculation of the momentum flux $\tau$ is reduced to the estimation of the drag coefficient $C_D$.

Field observations have shown that $C_D$ is strongly dependent on the wind velocity $U_{10}$ (Edson et al. 2013; Ayet & Chapron 2022). However, significant scatter in the data is often observed at a given wind speed due to difficulty in making momentum flux measurements in the open ocean and the role of variables other than wind speed (waves, current) in controlling the drag (Ayet & Chapron 2022). The relationship between $C_D$ and $U_{10}$ is well constrained when the wind and the waves are at equilibrium at moderate wind speed, i.e. $U_{10} \le 7-15$ m s⁻¹), but the scatter is larger at low wind speed due to misalignment between wind and waves (Ayet & Chapron 2022; Manzella et al. 2024). For wind speeds exceeding $U_{10} > 20$ m s⁻¹, the $C_D-U_{10}$ relationship exhibits a saturation in the field and laboratory experiments (Sroka & Emanuel 2021), with considerable scatter (Donelan et al.
2004; Sullivan & McWilliams 2010; Sraj et al. 2013; Curcic & Haus 2020; Sroka & Emanuel 2021). The precise saturation value and the underlying physical principle remain uncertain (Chen et al. 2007; Takagaki et al. 2016; Komori et al. 2018) while being a critical parameter to understanding the intensification of tropical (Sroka & Emanuel 2021) and extra-tropical (Gentile et al. 2021, 2022) cyclones. The main hypotheses that have been proposed to explain the saturation of $C_D$ are: sea spray generation at high wind speed (Bye & Jenkins 2006; Veron 2015), and airflow separation and flattening wave crests, and a marked reduction in surface roughness (Donelan et al. 2004). In this latter scenario, the role of wave breaking in flow separation and modulating pressure stress remains to be quantified in order to better constrain momentum flux formulation (Kudryavtsev et al. 2014).

Several studies have investigated momentum flux and the partitioning of interfacial stress into pressure and viscous contributions at moderate wind speeds, both experimentally (e.g. Buckley et al. 2020; Yousefi et al. 2020) and numerically for idealized waves using three-dimensional large-eddy simulations (e.g. Sullivan et al. 2000; Yang & Shen 2010). Iafrati et al. (2019) and Wu & Deike (2021) considered two-phase Navier-Stokes equations in a two-dimensional configuration with idealized airflow, noting that the two-dimensional configuration cannot capture the realistic turbulent boundary layer (linear-logarithmic velocity profile and three-dimensional fluctuations). Yang et al. (2018) and Lu et al. (2024) employed three-dimensional direct numerical simulations (DNS) to study the disturbance of the air boundary layer induced by breaking waves and to explore the associated modulation of heat transfer, while Wu et al. (2022) employed three-dimensional DNS to investigate momentum fluxes and wave growth for non-breaking wind waves.

So far, no fully resolved simulations have been conducted for strongly forced wind wave breaking up to $u_*/c \approx 1$, including both the wind-wave growth stage up to breaking and the following breaking event. Such configurations, where breaking conditions are reached through wind forcing, would be ideal for studying the role of the breaking and growing wave stages on the turbulent airflow and momentum flux, permitting the separation of the contribution of both phenomena, which are usually time- and ensemble-averaged in field and laboratory conditions.

Here, we consider the case of waves forced by a turbulent wind up to $u_*/c = 0.9$, including breaking events. The turbulent airflow and the water waves are fully resolved without any model for the wave shape, which can grow and break, and without any subgrid turbulence model. This is accomplished by using a two-phase Navier-Stokes solver, with a geometric volume-of-fluid method to reconstruct the wave interface. For different values of $u_*/c$, we estimate the momentum fluxes by analysing the growing, breaking and post-breaking stages of the wave field separately. Using this novel approach, we clarify the relation between the stress partition at the interface and the modulation of the velocity profile in the airflow region during the different stages of the wave dynamics. We extract an equivalent of the drag coefficient during the wave growth and breaking stages and discuss how their average values reproduce trends observed in the laboratory and field observations.

This paper is organized as follows. In § 2, we introduce the overall methodology, including (i) the governing equations that we solve and a brief summary of the numerical algorithm (§ 2.1), (ii) the numerical set-up (§ 2.2) and (iii) the physical dimensionless parameters (§ 2.3). In § 3, we provide a qualitative description of the evolution of the coupled system composed of the wind and the growing-breaking wave field. In § 4, we discuss the momentum flux over breaking waves and how the wave growth, steepening and eventual breaking modulate the velocity profile and the Reynolds stress. In § 5, we
connect our analysis to formulations of the drag coefficient over growing and breaking waves and compare them with experimental results. Conclusions are presented in § 6.

## 2. Methodology

In this section, we discuss the governing equations, numerical set-up and dimensionless parameters used to characterize the wind-forced wave growth and breaking processes.

### 2.1. Governing equations and numerical model

We investigate wind-forced breaking waves as a two-phase problem. We solve the incompressible, two-phase Navier-Stokes equations with surface tension, implemented in the open-source Basilisk solver (http://basilisk.fr/) (Popinet 2009, 2015), following the approach developed in Wu & Deike (2021) and Wu et al. (2022) to study wind-wave growth. Briefly, to distinguish the two phases, an indicator function $\mathcal{F}$ is introduced and set equal to 1 in the water phase within the volume $\Omega_w$ and 0 in the air phase within the volume $\Omega_a$. The two domains are separated by a zero-thickness interface $\Gamma$. The indicator function $\mathcal{F}$ is governed by the transport equation (Tryggvason et al. 2011):

$$ \frac{\partial \mathcal{F}}{\partial t} + \mathbf{u} \cdot \nabla \mathcal{F} = 0, \qquad (2.1) $$

where $\mathbf{u} = (u, v, w)$ is the one-fluid velocity, assumed continuous in the whole domain $\Omega = \Omega_w \cup \Omega_a$. The indicator function is employed to define the one-fluid density $\rho = \rho_w \mathcal{F} + \rho_a(1 - \mathcal{F})$ and the one-fluid dynamic viscosity $\mu = \mu_w \mathcal{F} + \mu_a(1 - \mathcal{F})$, where $\rho_w, \rho_a$ and $\mu_w, \mu_a$ are the densities and the dynamic viscosities of water and air, respectively. The transport of $\mathcal{F}$ is coupled with the incompressibility constraint and the Navier-Stokes equations for a Newtonian fluid. These equations, expressed in the one-fluid formulation (Tryggvason et al. 2011), read as follows:

$$ \nabla \cdot \mathbf{u} = 0, \qquad (2.2) $$

$$ \frac{\partial(\rho\mathbf{u})}{\partial t} + \nabla \cdot (\rho\mathbf{u}\mathbf{u}) = -\nabla p + \nabla \cdot (2\mu\mathbf{D}) + \rho|\mathbf{g}|e_z + \sigma\varkappa\delta_{\Gamma}\mathbf{n}, \qquad (2.3) $$

where $p$ is the pressure, $\mathbf{D}$ is the strain rate tensor, i.e. $\mathbf{D} = (\nabla\mathbf{u} + \nabla\mathbf{u}^T)/2$, $|\mathbf{g}|$ is the modulus of the gravitational acceleration, $e_z = (0, 0, -1)$ is unit vector oriented like gravity, $\sigma$ is the surface tension coefficient, $\delta_{\Gamma}$ is a Dirac distribution function satisfying the identity $\int_{\Gamma} \delta_{\Gamma} dS = 1$, $\mathbf{n}$ is the interface normal vector pointing outward to the liquid domain and $\varkappa = \nabla \cdot \mathbf{n}$ is the interfacial curvature.

Equations (2.1), (2.2) and (2.3) are solved using an adaptive mesh refinement (AMR) strategy on an octree grid, as implemented in Basilisk and described in Popinet (2015) and Van Hooft et al. (2018). The use of AMR significantly reduces the computational cost while efficiently representing different multiscale processes (Mostert et al. 2022). We employ a conservative and diffusion-free geometric volume-of-fluid method (Tryggvason et al. 2011) to discretize equation (2.1). The equation for the indicator function (2.1) is then solved together with the momentum equation (2.3) using the Bell-Colella-Glaz method (Bell et al. 1989) for the advection part and a standard second-order finite-volume scheme for the viscous term. The viscous term is advanced semi-implicitly in time using a Crank-Nicolson scheme, whereas an explicit treatment is kept for the advection part.
<figure><img src="image_3.png" /><figcaption>Figure 1. Computational domain and physical configuration illustrating the initial condition for the airflow, wave and water field. The air and water mean heights are $(L_0 - h_W)/\lambda$ and $h_W/\lambda$. The airflow is a fully developed turbulent boundary layer (mean profile in light-blue line, while turbulent eddies are illustrated in black). The wave field and the water region are initialized using an irrotational third-order Stokes wave solution (Deike et al. 2015). The initial wave profile $\eta_0$ has zero spatial mean and a steepness $a_0k = 0.3$. In the surface contour, dark-blue regions denote wave troughs, while yellow regions indicate wave crests.</figcaption></figure>

The momentum equation is discretized in a conservative and consistent manner (Popinet 2018; Mostert et al. 2022) to ensure accurate and stable integration of (2.3) in the presence of a large difference in the density and viscosity of the two phases. The capillary and gravity forces are discretized using a well-balanced formulation (Popinet 2018), maintaining an exact equilibrium among pressure gradient, capillary and gravity forces under static conditions, and minimizing the generation of artificial parasitic currents. These features are important for studying phenomena such as wave breaking (Deike et al. 2015, 2016; Mostert et al. 2022), wind-wave growth (Wu & Deike 2021; Wu et al. 2022) and other classes of two-phase turbulent flows (see e.g. Riviere et al. 2021; Perrard et al. 2021; Farsoiya et al. 2023).

## 2.2. Configuration and initialization

The governing equations (2.1), (2.2) and (2.3) are solved in a cubic computational box, as illustrated in figure 1, of size $L_0 = 4\lambda$, i.e. $[-2\lambda, 2\lambda] \times [-2\lambda, 2\lambda] \times [-h_W, 4\lambda - h_W]$, where $\lambda$ is the wavelength and $h_W$ is the mean water depth. The streamwise, spanwise and vertical directions correspond to the $x$, $y$ and $z$ directions, respectively. Periodic boundary conditions are prescribed for all the variables, i.e. **u**, **p** and $\mathcal{F}$, along the $x-y$ horizontal directions. On the top and bottom domain boundaries, which correspond to the grey geometric planes at $z = (L_0 - h_W)/\lambda$ and $z = -h_W/\lambda$ in figure 1, a homogeneous Neumann boundary condition is prescribed for the pressure $p$, the indicator function $\mathcal{F}$ and the two horizontal velocity components $u$ and $v$, while a no-penetration boundary condition is imposed on the vertical velocity component $w$.

Figure 1 illustrates the three domains of interest: the wave field occupies the region $z = \eta(x, y)$ and separates the water and air flows. The water flow occupies the region where
$z < \eta(x, y)$, with a mean depth of $h_W$. The airflow occupies the region where $z > \eta(x, y)$ and has mean height ($L_0 - h_W$). The wave field, water and air velocity must be properly initialized at the beginning of the simulation.

We initialize the wave and the associated water velocity field using a nonlinear potential flow solution $\Phi_0$ with free surface $\eta_0$ (Lamb 1993), and consider a third-order expansion of the system (third-order Stokes wave), with a given amplitude $a_0$ and wavenumber $k = 2\pi/\lambda$. Once $\Phi_0$ is known, the initial velocity in the water $\mathbf{u}_{w,0}$ can be readily evaluated from the velocity potential as $\mathbf{u}_{w,0} = \nabla\Phi_0$; for the expressions of $\eta_0$ and $\Phi_0$, see Lamb (1993) and Deike et al. (2015). We note that this initialization prescribes an initial orbital velocity field in the water characterized by a zero Eulerian mean. It is worth mentioning that several previous works (Chen et al. 1999; Iafrati 2009; Deike et al. 2015, 2016; Yang et al. 2018; Mostert et al. 2022) have demonstrated the relevance of the Stokes waves solution for the simulation of breaking waves with numerical results that can be accurately compared with laboratory experiments. However, unlike the Stokes wave solution, which disregards surface tension effects, the current approach incorporates them. Thus, based on the dispersion relation, the wave speed reads as $c = \sqrt{|g|/k + \sigma k/\rho_w}$. We consider the linear phase speed for gravity-capillary waves, as the typical nonlinear Stokes wave correction remains negligible during the initial growth stage. During the breaking stage, the wave speed undergoes a slight slowdown, which is not accounted for by this correction (Banner et al. 2014). However, the nominal wave speed set by the initial conditions remains an excellent approximation of the effective wave speed throughout the simulation.

We initialize the airflow region using a fully developed turbulent flow field at the desired friction Reynolds number, as defined in § 2.3, following Wu et al. (2022). This initialization involves a precursor simulation conducted independently in a single-phase set-up and employing the same domain. During this precursor simulation, the wave field with profile $\eta_0$ remains at rest, and a no-slip/no-penetration boundary condition is enforced on the wave surface for the velocity field, using the embedded boundary method (Johansen & Colella 1998) available within the Basilisk framework (Ghigo et al. 2021). The precursor simulation is performed long enough until the turbulent airflow achieves a statistically steady state by adding an external body force per unit mass acting in the streamwise direction, i.e. $\partial p_0/\partial x(1 - \mathcal{F})\mathbf{e}_x$ with $\mathbf{e}_x = (1, 0, 0)$, on the right-hand side of the momentum equation (2.3). In the expression of the body force, $\partial p_0/\partial x$ is a uniform pressure gradient driving the flow and, here, reads as

$$ \frac{\partial p_0}{\partial x} = \frac{\rho_a u_*^2}{L_0 - h_W}. \qquad (2.4) $$

The imposed pressure gradient sets the nominal friction velocity $u_*$ and prescribes the total stress $\rho_a u_*^2$ on the wave field. Once a statistically steady state is achieved for the precursor, the resulting fully developed turbulent field is employed as an initial condition for the airflow region in the two-phase simulations.

Once the simulation begins, the airflow and wave field dynamically evolve without any prescribed conditions at the two-phase interface. This fully coupled system undergoes a short self-adjustment stage due to changes in boundary conditions at the interface and the motion of the water layer (Wu et al. 2022). As discussed by Wu et al. (2022), the self-adjustment period is short compared with other processes of interest. In this work, we verify that the self-adjustment stage lasts less than $2\omega t$ in all cases, where $\omega = 2\pi/T_0$ is the angular frequency and $T_0 = \lambda/c$ is the wave period; therefore, all analyses are conducted for $\omega t > 2$.
## 2.3. Non-dimensional parameters

The flow field in the air and water phases depends on several dimensional parameters. Here, the number of independent physical dimensions is three, i.e. mass, length and time. According to Buckingham's $\pi$ theorem, the 11 physical variables ($\rho_a$, $\rho_w$, $\mu_a$, $\mu_w$, $L_0 - h_W$, $h_W$, $\lambda$, $a_0$, $\sigma$, $|g|$, $u_*$) can be reduced to 8 dimensionless groups so that the problem can be described by the following dimensionless groups (note that other combinations would have been possible):

$$ \left( \mathrm{Re}_{*,\lambda}, \mathrm{Re}_W, \mathrm{Bo}, \frac{\rho_w}{\rho_a}, \frac{L_0 - h_W}{\lambda}, \frac{h_W}{\lambda}, a_0 k, \frac{u_*}{c} \right). \qquad (2.5) $$

In (2.5), $Re_{*,\lambda} = \rho_a u_* \lambda / \mu_a$ and $Re_W = \rho_w c \lambda / \mu_w$ represent the friction and wave Reynolds numbers, respectively, reflecting the balance between inertia and viscous effects in the regions adjacent to the wave field, for the air and water phases. A third dependent friction Reynolds number is that based on the size of air boundary layer $L_0 - h_W$, i.e. $Re_* = Re_{*,\lambda} (L_0 - h_W) / \lambda = \rho_a u_* (L_0 - h_W) / \mu_a$, which reflects the importance of inertia over viscous forces in the region well above the wave field. In Appendix E, we show that $Re_{*,\lambda}$ is the physically relevant group in order to compare cases at different Reynolds numbers and, therefore, $Re_{*,\lambda}$ is used in conjunction with $Re_W$. The Bond number $\mathrm{Bo} = (\rho_w - \rho_a) |g| / (\sigma k^2)$ provides the ratio between gravitational and restoring capillary forces and $\rho_w / \rho_a$ represents the density ratio, which is set to that of water and air. Next, we have a set of ratios of length scale, or geometric parameters: $(L_0 - h_W) / \lambda$ represents the ratio between the mean airflow height and the wavelength; $h_W / \lambda$ is the ratio between the mean height of the water depth and the wavelength $\lambda$. Sensitivity to these parameters is shown in Appendix E, and we verify that the physical conclusions do not depend on their specific values. Finally, $a_0 k$ represents the wave's initial steepness and $u_*/c$ defines the ratio between the friction velocity and wave phase speed.

Here, we focus on the interaction between turbulent wind and breaking waves at high wind speed and systematically vary the ratio $u_*/c$ from 0.3 to 0.9 while also performing sensitivity tests in terms of $Re_{*,\lambda}$ and geometrical parameters, as discussed in § 4. The range of $u_*/c$ from 0.3 to 0.9, is relevant to discuss small-scale physics of wind-wave fields during tropical cyclone conditions (Sroka & Emanuel 2021).

We consider four water waves and prescribe the airflow's mean height $(L_0 - h_W) / \lambda = 3.36$. This value ensures that the air boundary layer is more than three times larger than one single wavelength, avoiding any confinement effect due to the top boundary, as discussed in Wu et al. (2022). We set the mean water depth $h_W / \lambda = 0.64$, which is considered adequate to satisfy the deep-water assumption (Deike et al. 2015, 2016; Yang et al. 2018). We consider $Re_W = 2.55 \times 10^4$ and $\mathrm{Bo} = 200$, so that the waves are within the gravity regime. We consider the air-water density ratio, $\rho_w / \rho_a = 816$, while the dynamic viscosity ratio (fixed through the choices of the Reynolds numbers) remains well above unity, as reported in table 1. The initial wave steepness is $a_0 k = 0.3$, below breaking threshold, typically [0.32 - 0.33] (Deike et al. 2015) for the present configuration. This choice ensures that the wind drives the wave field to break.

In the present configuration, and unlike in laboratory experiments or open-ocean conditions, $u_*/c$ is independent of $a_0 k$, which allows us to isolate the effect of $u_*/c$ on the flow at the same $Re_{*,\lambda}$ (by varying the wave speed $c$ while keeping $u_*$ fixed from the same precursor). Furthermore, we work at constant $Re_W$ so that the growth rate of the wave field is controlled by the wind forcing $u_*/c$, while the viscous dissipation is the same (Wu et al. 2022), but note that the Reynolds number (air and water side) and Bond number we consider are significantly smaller than those for ocean waves with wavelength 1 m or larger
<table><thead><tr><th>u<sub>*</sub>/c</th><th>Re<sub>*,λ</sub></th><th>Re<sub>W</sub></th><th>μ<sub>w</sub>/μ<sub>a</sub></th><th>(L<sub>0</sub> - h<sub>W</sub>)/λ</th><th>Bo</th><th>Le</th></tr></thead><tbody><tr><td>0.30</td><td>214</td><td>2.55 × 10<sup>4</sup></td><td>22.84</td><td>3.36</td><td>200</td><td>10</td></tr><tr><td>0.40</td><td>214</td><td>2.55 × 10<sup>4</sup></td><td>17.13</td><td>3.36</td><td>200</td><td>10</td></tr><tr><td>0.50</td><td>214</td><td>2.55 × 10<sup>4</sup></td><td>13.71</td><td>3.36</td><td>200</td><td>10</td></tr><tr><td>0.70</td><td>214</td><td>2.55 × 10<sup>4</sup></td><td>9.79</td><td>3.36</td><td>200</td><td>10</td></tr><tr><td>0.90</td><td>214</td><td>2.55 × 10<sup>4</sup></td><td>7.61</td><td>3.36</td><td>200</td><td>10</td></tr><tr><td>0.90</td><td>53.5</td><td>2.55 × 10<sup>4</sup></td><td>1.90</td><td>3.36</td><td>200</td><td>10</td></tr><tr><td>0.90</td><td>107</td><td>2.55 × 10<sup>4</sup></td><td>3.81</td><td>3.36</td><td>200</td><td>10</td></tr><tr><td>0.90</td><td>214</td><td>2.55 × 10<sup>4</sup></td><td>7.61</td><td>3.36</td><td>200</td><td>11</td></tr><tr><td>0.90</td><td>107</td><td>2.55 × 10<sup>4</sup></td><td>3.81</td><td>6.72</td><td>200</td><td>11</td></tr></tbody><tfoot><tr><td colspan="7">Table 1. Summary of the simulated cases for different values of u<sub>*</sub>/c. In the table, Re<sub>*,λ</sub> = ρ<sub>a</sub>u<sub>*</sub>λ/μ<sub>a</sub>, Re<sub>W</sub> = ρ<sub>w</sub>cλ/μ<sub>w</sub>, μ<sub>w</sub>/μ<sub>a</sub>, (L<sub>0</sub> - h<sub>W</sub>)/λ and Bo = (ρ<sub>w</sub> - ρ<sub>a</sub>)|g|/(σk<sup>2</sup>) are as defined in (2.5). The case with (L<sub>0</sub> - h<sub>W</sub>)/λ = 3.36 corresponds to four waves per box size, whereas the case with (L<sub>0</sub> - h<sub>W</sub>)/λ = 6.72 to eight waves per box size. For the different cases, the initial steepness is set equal to a<sub>0</sub>k = 0.3, and the density ratio is taken as ρ<sub>w</sub>/ρ<sub>a</sub> = 816.</td></tr></tfoot></table>

due to computational limitations. The underlying assumption to compare our results with laboratory or field conditions is that we perform simulations in quasi-asymptotic regimes of high values of these parameters (Deike et al. 2016; Mostert et al. 2022; Wu et al. 2022). The simulations are summarized in table 1.

The numerical grid in our simulations is adaptive, featuring a minimum grid size $\Delta = L_0/(2^{2L_c})$, where Le represents the maximum level of refinement. The AMR technique significantly reduces computational costs by maintaining a highly refined grid near the interface and in the boundary layers while allowing coarser grids in the bulk airflow, provided that refinement criteria are met (Popinet 2015; Van Hooft et al. 2018). The AMR technique has been shown to be accurate for homogeneous and isotropic turbulent flow (Riviere et al. 2021; Farsoiya et al. 2023), and wall-bounded turbulent flows and the present configuration in Wu et al. (2022), where the grid is dynamically adapted with respect to the norm of the second derivative of the velocity (in the air and water phases) and of the volume fraction. Following Wu et al. (2022), the refinement criteria for the air and water velocity components and the volume fraction are set equal to $\epsilon_{ua} = 0.3u_*$, $\epsilon_{uw} = 10^{-3}c$ and $\epsilon_F = 10^{-4}$, respectively.

In most simulations, Le = 10, with the maximum number of grid points reaching up to $70 \times 10^6$, which corresponds to 7% of an equivalent $1024^3$ uniform grid. Additionally, we present in Appendix C a grid refinement study up to Le = 11 for $u_*/c = 0.9$. In this case, the maximum number of grid points is around $500 \times 10^6$, which corresponds to about 6% of an equivalent $2048^3$ uniform grid. The significant reduction in the number of grid cells with AMR at Le = 10 - 11, compared with a uniform Cartesian grid, makes AMR an attractive approach for investigating wind-forced breaking waves.

The computational costs of these simulations include generating a precursor simulation and running the two-phase simulations. Four precursors are used for varying air-side Reynolds number ($Re_{*,\lambda} = 53.5 - 107 - 214$) and the ratio $(L_0 - h_W)/\lambda = 3.36 - 6.72$, amounting to $\approx 1.5 \times 10^5$ CPU hours each. Each two-phase simulation at Le = 10 requires $\approx 4.2 \times 10^5$ CPU hours, while the two simulations at Le = 11 require $\approx 1.20 \times 10^6$ CPU hours each. We employ 384 processors for the precursor simulations and for the two-phase cases 480 at Le = 10, and 980 processors at Le = 11. The total cost of the simulation campaign is about $6 \times 10^6$ CPU hours.
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

Figure 2. Wind-forced breaking waves for the case $u_*/c = 0.9$ at $\omega t = 20 - 40 - 70 - 150$ (a) just before breaking, (b) during breaking, (c) during second growing stage and (d) final stage with a more three-dimensional field and microbreaking, respectively. The turbulent airflow, i.e. the wind, and the waves move parallel along the positive streamwise direction, $x$. In all the panels, the plane on the left contains the contour of the local spanwise vorticity $\Omega_y = (\partial w/\partial x - \partial u/\partial z)$ normalized by the wave angular velocity $\omega = 2\pi/T_0$ and the plane on the right contains the contour of the local streamwise velocity $u$ normalized by the friction velocity $u_*$.

### 3. Evolution of wind-forced breaking waves

This section discusses the evolution of the fully coupled system composed of the turbulent wind, the wave field and the water column.

#### 3.1. Wave interface evolution

Figure 2 shows the wind-wave system for the case $u_*/c = 0.9$ at four characteristic times. At $\omega t = 20$ (figure 2a), 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 2b), the wave field grows again, starting with a smaller steepness (figure 2c) taken at $\omega t = 70$ and finally breaks again at $\omega t = 150$ (figure 2d). 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,
<figure><img src="image_8.png" /><figcaption>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_*/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_*/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.</figcaption></figure>

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):

$$E_W(t) = \rho_w |\mathbf{g}| \int_{\Omega_w} z \, dV - E_{p,0}, \quad (3.1)$$

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 |\mathbf{g}| L_0^2 \int_0^{h_w} z \, dz = \rho_w |\mathbf{g}| (h_w L_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 1993) since it is also valid for nonlinear and breaking waves contrary to the formula based on linear approximation $E_W(t) = \rho_w |\mathbf{g}| \int_{\Gamma} \sqrt{2(\eta - \bar{\eta})^2} dS$ (Lamb 1993) which is not well defined during the breaking stage.

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 - \bar{\eta})^2 dS}$. Alternatively, the wave steepness can be determined using
the surface elevation derivatives as $S = \max(\sqrt{(\partial\eta/\partial x)^2 + (\partial\eta/\partial y)^2})$. Note that under the assumption of linear wave theory, $S$ reduces to $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_*/c$ and the dissipation. In each curve, the potential wave energy displays small oscillations of period $T_0/2$, already observed in Iafrati (2009) and Deike et al. (2015) 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_*/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_*/c = 0.3$. Performing simulation at higher wave Reynolds number would also lead to wave energy growing in time for such value of $u_*/c = 0.3$ (see Wu & Deike 2021; Wu et al. 2022; Zhang et al. 2023). The wave growth rate increases with $u_*/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_*/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. 2015, 2016) and experiments on focusing wave packets (Drazen et al. 2008). 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 2d). 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. 2008; Deike et al. 2015).

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_*/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
<figure><img src="image_9.png" /><figcaption>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/λ = 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_*$, (d–f) Contours of the instantaneous spanwise vorticity in the airflow sampled at the middle plane $y/λ = 0$ and normalized by the angular velocity, i.e. $Ω_y/ω$. (g–i) Contours of the spanwise vorticity in the airflow, averaged along the spanwise direction and normalized by the angular velocity, i.e. $Ω_y/ω$. All the panels are plotted in the region $-0.2λ ≤ z ≤ 0.4λ$ for $ωt = [20, 35, 70]$ for $u_*/c = 0.9$.</figcaption></figure>

a reference frame moving at the wave speed, i.e. $u-c$, in the region close to the wave surface ($-0.2 ≤ z ≤ 0.4λ$) and sampled at the spanwise mid-plane ($y/λ = 0$), at several characteristic times (during growth and breaking) for $u_*/c = 0.9$.

As the wave field grows and approaches breaking, i.e. $ω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_*$. 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 $ωt > 35$, when the wave amplitude is greatly reduced, the regions with negative $u-c$ reduce compared to the growing stage (see figure 4b,c).

As noted in (Veron et al. 2007), although flow reversal suggests airflow separation from the wave field, a more rigorous approach involves examining the spanwise vorticity $Ω_y$, which is a Galilean-invariant quantity. The spanwise vorticity normalized by the angular velocity, i.e. $Ω_y/ω$, 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/λ = 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. (2020).

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 1976; Reul et al. 1999; Buckley et al. 2020), and is analysed in the next section.

## 4. Momentum fluxes over strongly forced breaking waves

We now analyse the momentum fluxes between the turbulent airflow and the growing and breaking waves. In § 4.1, we examine the partition of the total momentum flux into pressure and viscous forces and their temporal variation during a breaking cycle. In § 4.2, we investigate how the velocity profile and the Reynolds stress of the airflow are modulated by waves during the growing and breaking stages.

### 4.1. Momentum budget

The momentum flux, $τ_t = (τ_{t,x}, τ_{t,y}, τ_{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 1998; Sullivan et al. 2000). For water waves moving parallel to the wind in the streamwise $x$ direction, the first force component, i.e. $τ_{t,x}$, is typically dominant and reads

$$τ_{t,x} = τ_{p,x} + τ_{v,x} = - \int_{Γ} p \mathbf{n} \cdot \mathbf{e}_x dS + 2μ_a \int_{Γ} (\mathbf{D}n) \cdot \mathbf{e}_x dS, \quad (4.1)$$

where both the pressure $p$ and the strain rate tensor $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.

We start by qualitatively inspecting both contributions during the growing and breaking stages. Figure 5 reports the instantaneous distribution of $p_a n_x$ and $2μ_a (\mathbf{D}n) \cdot \mathbf{e}_x$ for the case $u_*/c = 0.9$ at three physical times $ω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 $τ_{sx} = 2μ_a (\mathbf{D}n) \cdot \mathbf{e}_x$ for brevity and, therefore, $τ_{v,x} = \int_{Γ} τ_{sx} dS$ in (4.1). The distributions $p_a n_x$ and $τ_{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. (2022) and Appendix B for details).

Figure 5 shows that both $p_a n_x$ and $τ_{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 $ωt = 20-70$, consistent with the measurements in Veron et al. (2007). 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 $p n_x$ are one order of magnitude larger than the peak in $τ_{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
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>Figure 5. Surface distribution of the pressure stress $p_a n_x$ (a,c,e) and viscous stress $\tau_{sx} = 2\mu_a (\mathbf{D}n) \cdot \mathbf{e}_x$ (b,d,f) for $u_*/c = 0.9$. Note that both quantities are normalized by the total imposed stress $\rho_a u_*^2$ at the interface. The dot-dashed lines in all the panels represent the position of the wave crest.</figcaption></figure>

(Buckley et al. 2020). 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$:

$$ \underbrace{\rho_a \frac{\partial}{\partial t} \int_{\Omega_a} u \, dV}_{\rho_a \partial U / \partial t} + \underbrace{\rho_a \int_{\Gamma} u(\mathbf{u}_r \cdot \mathbf{n}) \, dS}_{\rho_a \phi_{c,x}} = \tau_{p,x} + \tau_{v,x} + \underbrace{\rho_a u_*^2 A_{\Gamma}}_{\Pi_f}, \quad (4.2) $$

where $dV$ 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}) dV = \rho_a u_*^2 A_{\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_{v,x}$, and the latter, $\rho_a \phi_{c,x}$, accounts for the momentum flux originated from a non-zero relative velocity $\mathbf{u}_r$ between the airflow and the wave field. Since the magnitude of $\mathbf{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_{v,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_a u_*^2 A_{\Gamma}$ balances the sum of pressure and viscous drag (Janssen 2004; Buckley et al. 2020; Funke et al. 2021). More details on the evaluation of the momentum fluxes are provided in Appendix B.

We analyse the temporal variation of the different terms in the momentum flux (4.2) for two representative cases ($u_*/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_*/c = 0.5$ and $[0 - 20]$ for $u_*/c = 0.9$, the shorter growing cycle for $u_*/c = 0.9$ being due to the higher wind forcing).

When the wave breaks, the momentum flux due to pressure $\tau_{p,x}/(\rho_a u_*^2 A_{\Gamma})$ drops with a corresponding acceleration of the flow, i.e. $\rho_a \partial U / \partial t / (\rho_a u_*^2 A_{\Gamma})$ increases, the drop being larger for the strongest wind forcing. Once the breaking stage ends, around $\omega t = 75$ for $u_*/c = 0.5$ and around $\omega t = 35$ for $u_*/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_*/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. (2018), who found a negligible role of the convective term in the budget. Next, the viscous contribution $\tau_{v,x}$ represents a small but not negligible contribution in the forces $\tau_{v,x}$. This is consistent with the high initial wave steepness, i.e. $a_0 k = 0.3$. Notably, upon wave breaking, the viscous contribution experiences a slight increase to partially compensate for the loss of the pressure force.
<figure><img src="image_16.png" /><figcaption>Figure 6. Contributions of the momentum budget in the streamwise direction, as in (4.2), for (a) $u_*/c = 0.5$ and (b) $u_*/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_{v,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_a u_\tau^2 A_\Gamma$. For both cases, the normalized variation in the wave energy $E_W / E_{W,0}$ is reported in the top panel.</figcaption></figure>

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. (2021)).

## 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. (2000) and Wu et al. (2022), 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 shows the streamwise velocity profile spatially averaged along the horizontal directions as a function of the vertical wave-following coordinate $\zeta$ for $u_*/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_*$, which is negligible compared with the maximum streamwise velocity in the airflow (see figure 7).
<figure><img src="image_17.png" /><figcaption>Figure 7. Streamwise velocity profile normalized by the nominal friction, $\langle u_a \rangle/u_*$, as a function of vertical wave-following coordinate $\zeta/\lambda$ in the airflow for (a) $u_*/c = 0.5$ and (b) $u_*/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 3b). The dot-dashed curves represent the instantaneous values during the breaking stage, i.e. $\omega t \in [58-98]$ for $u_*/c = 0.5$ and $\omega t \in [22-42]$ for $u_*/c = 0.9$.</figcaption></figure>

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 < 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_*/c = 0.5$ with the case $u_*/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 = v_a/u_*^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_*/c)Re_{*,\lambda} > 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_*^2$, i.e. $-(\langle u'w' \rangle/u_*^2)$, for $u_*/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.

For both cases, $-(\langle u'w' \rangle/u_*^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_* = 720$ (Pope 2000). 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.
<figure><img src="image_18.png" /><figcaption>Figure 8. Reynolds shear stress normalized by the square of the nominal friction, $-\langle u'w' \rangle/u_*^2$, as a function of the vertical wave-following coordinate $\zeta$ for (a) $u_*/c = 0.5$ and (b) $u_*/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_*/c = 0.5$ and $\omega t \in [22-42]$ for $u_*/c = 0.9$. The dashed green curve represents the Reynolds stress on a flat stationary surface at $Re_* = 720$.</figcaption></figure>

During the growing stage, the peak in the Reynolds stress is 0.62 for $u_*/c = 0.5$, whereas it is 0.76 for $u_*/c = 0.9$. This larger peak is due to the increased drag, primarily from the pressure component, at higher $u_*/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_*/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_*^2$, is larger as $u_*/c$ increases, while when the wave breaks, the turbulent drag is reduced, and such loss increases with $u_*/c$. We note that both cases are run at the same $Re_{*,\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 2000; Cimarelli et al. 2023), composed of the viscous sublayer ($0 < \zeta^+ < 5$), the buffer layer ($5 < \zeta^+ < 30$) and log-law region ($\zeta^+ > 30$), where $\zeta^+$ is the vertical coordinate in wall units $\zeta^+ = \zeta u_*/v_a$ with $v_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 < \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.

It is also worth mentioning that the streamwise velocity profile at the same $Re_*$ 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 1998). This observation agrees with experimental
<figure><img src="image_19.png" /><figcaption>Figure 9. Streamwise velocity profile normalized by the friction velocity, $\langle u_a^+ \rangle = \langle u_a \rangle / u_*$, as a function of vertical wave-following coordinate (in wall units) $\zeta^+ = \zeta u_*/v_a$ with $v_a = \mu_a / \rho_a$ for (a) $u_*/c = 0.5$ and (b) $u_*/c = 0.9$, both at $Re_{*,\lambda} = 214$ ($Re_* = 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_* = 720$ for a flat stationary surface.</figcaption></figure>

works where a downshift of the velocity profile was observed as the wave steepness was increased (Buckley & Veron 2016; Buckley et al. 2020).

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:

$$ \frac{\langle u_{a,f} \rangle - \langle u_{a,i} \rangle}{u_*} = \frac{1}{\kappa} \log \left( \frac{\zeta_f^+ - \zeta_i^+}{z_0^+} \right), \quad (4.3) $$

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. 2000; Lin et al. 2008) 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_{*,\lambda} = 53.5$, $107$ and $214$, respectively. The final extent $\zeta_f^+$ is taken equal to a common value of $\zeta_f^+ \approx Re_{*,\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 (1925) and Donelan et al. (2006) 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_*/c$, whereas when the wave field breaks, its value drops and decreases with $u_*/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,
<table><thead><tr><th>Re<sub>*,λ</sub></th><th>(L<sub>0</sub> - h<sub>W</sub>)/λ</th><th>u<sub>*</sub>/c</th><th>z<sup>+</sup><sub>0,G<sub>1</sub></sub></th><th>z<sup>+</sup><sub>0,G<sub>2,a</sub></sub></th><th>z<sup>+</sup><sub>0,G<sub>2,b</sub></sub></th><th>z<sup>+</sup><sub>0,F</sub></th></tr></thead><tbody><tr><td>214</td><td>3.36</td><td>0.3</td><td>1.05</td><td>—</td><td>—</td><td>—</td></tr><tr><td>214</td><td>3.36</td><td>0.4</td><td>1.10</td><td>0.87</td><td>—</td><td>—</td></tr><tr><td>214</td><td>3.36</td><td>0.5</td><td>1.14</td><td>0.83</td><td>—</td><td>—</td></tr><tr><td>214</td><td>3.36</td><td>0.7</td><td>1.18</td><td>0.80</td><td>1.07</td><td>0.89</td></tr><tr><td>214</td><td>3.36</td><td>0.9</td><td>1.21</td><td>0.77</td><td>1.15</td><td>0.85</td></tr><tr><td>53.5</td><td>3.36</td><td>0.9</td><td>0.26</td><td>0.17</td><td>—</td><td>—</td></tr><tr><td>107</td><td>3.36</td><td>0.9</td><td>0.55</td><td>0.36</td><td>0.50</td><td>0.40</td></tr><tr><td>107</td><td>6.72</td><td>0.9</td><td>0.57</td><td>0.38</td><td>0.52</td><td>0.42</td></tr></tbody></table>

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.

We note that the drag reduction mentioned here is studied at a fixed $Re_{*,\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_*/c = 0.9$, two at $Re_{*,\lambda} = 53.5 - 107$ with $(L_0 - h_w)/\lambda = 3.36$ and one at $Re_{*,\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_{*,\lambda}$. This supports that the drag modulation induced by wave breaking is primarily controlled by $u_*/c$ and is less sensitive to $Re_{*,\lambda}$ and $(L_0 - h_w)/\lambda$.

## 5. Drag coefficients over growing and breaking waves

In this section, we connect our analysis of the momentum flux under growing and breaking waves at high wind speeds to formulations of the drag coefficient. We split our analysis between the growing waves and the breaking waves, for different $u_*/c$ and $Re_{*,\lambda}$. Two formulations are discussed: one based on first principles and directly integrating the momentum flux, and the one used in the laboratory and field observations.

### 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. $(u_a)/u_*$ 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}$:

$$C_{D,a} = \frac{2\tau_{p,x}}{\rho_a \overline{U}_{ref}^2 (\bar{\zeta} = \zeta_{ref}) A_{\Gamma}} \quad (5.1)$$
<figure><img src="image_20.png" /><figcaption>Figure 10. Aerodynamic drag coefficient $C_{D,a}$ (defined by (5.1)) for different $u_*/c$ in the growing (blue colours) and breaking (red colours) time intervals. For $u_*/c < 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. 2022). For $u_*/c > 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_1, G_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_*/c$, while breaking induces a decrease in drag with increasing $u_*/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.</figcaption></figure>

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

$$ \bar{\tau}_{p,x} = \frac{1}{\Delta T_{G,B}} \int_{\Delta T_{G,B}} \tau_{p,x} \, dt, \qquad (5.2) $$

where $\Delta 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 $\bar{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_*/c$.

The choice of the reference far-field boundary layer velocity $\bar{U}_{ref}$ evaluated at a reference height $\bar{\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 $\bar{U}_{ref}$ increases as the flow accelerates in the breaking stage. Here, the requirements are satisfied typically if $\lambda/4 < \zeta_{ref} < \lambda$, and we consider $\zeta = \lambda/2$, similar to wind-wave growth sheltering discussion (Jeffreys 1925; Donelan et al. 2006). Different $\zeta_{ref}$ within the interval only changes the magnitude $C_{D,a}$, without altering the discussion.

Figure 10 shows the estimated $C_{D,a}$ as a function of $u_*/c$ and for both growing and breaking stages. We also re-analyse the non-breaking growing simulations from Wu et al. (2022) at lower $u_*/c$ and varying initial slope. In the non-breaking growing regime ($u_*/c < 0.35$), $C_{D,a}$ increases with $u_*/c$, with different values of $C_{D,a}$ at the same $u_*/c$ corresponding to different initial wave steepness (see also Wu et al. 2022).
For the range of $u_*/c$ that includes growth and breaking ($u_*/c > 0.35$), $C_{D,a}$ follows an increasing trend with $u_*/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_*/c = [0.7 - 0.9]$ at $Re_{*,\lambda} = 214$ and $u_*/c = 0.9$ at $Re_{*,\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_{*,\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. 2004; Janssen 2004; Bye & Jenkins 2006; Buckley et al. 2020; Curcic & Haus 2020; Sroka & Emanuel 2021), by considering the classic oceanographic definition of the drag coefficient at 10 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_a C_D U_{10}^2$ with $\tau_{t,x} = \rho_a u_*^2$ since at 10 m height the airflow and the waves are assumed to be in equilibrium. For neutral atmospheric conditions, the velocity profile at $\bar{z}$ can be assumed to follow the logarithmic law of the wall (Edson et al. 2013; Ayet & Chapron 2022), i.e. $U_{10} = (u_*/\kappa) \log(\bar{z}/z_0)$ as in (4.3), and $C_D$ becomes

$$C_D(\bar{z} = 10 \text{ m}) = \frac{\kappa^2}{\log^2(\bar{z}/z_0)} \qquad (5.3)$$

Using the DNS data generated for this work and those reported in Wu et al. (2022), 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_*/c$ and $Re_{*,\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_*/c$ exhibit a similar qualitative trend to that of $C_{D,a}$ with $u_*/c$. We observe a systematic increase of the drag in the growing stage (dark triangles and blue symbols) with $u_*/c$, while the drag in the breaking stage (red
symbols) decreases with $u_*/c$. When averaging over both growing and breaking regimes, we observe a saturation at high wind speed ($u_*/c > 0.4$). As for the aerodynamic drag coefficient $C_{D,a}$, the results are not very sensitive to $Re_{*,\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_*/c < 0.3$ and $C_D$ between 1.5 and $2.0 \times 10^{-3}$ for $u_*/c > 0.4$, closely match those reported in laboratory experiments (Buckley et al. 2020; Curcic & Haus 2020) for similar values of $u_*/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_*/c$ and $a_{rmsk}$ 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_*/c$ may differ from our simulations. Similarly, directly comparing realistic $u_*/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.

The present analysis can be related to the discussion from Sullivan et al. (2018) and Wu et al. (2022) on the control of the form drag by the wave slope, $a_{rmsk}$, 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_*/c$.

Moreover, the drag coefficient $C_D$ is shown as a function of $a_{rmsk}$ in figure 11. During the growth phase of the wave field, $C_D$ increases with $a_{rmsk}$ up to approximately $a_{rmsk} \approx 0.15$ for non-breaking waves with finite steepness (Wu et al. 2022). As $a_{rmsk}$ reaches larger values, wave breaking initiates around $a_{rmsk} \approx [0.20 - 0.23]$. Focusing on the growing phase of the wave field, $C_D$ increases with $a_{rmsk}$ for the different $u_*/c$; however, as we approach the largest values, $a_{rmsk}$ saturates, leading to a corresponding saturation in $C_D$. Conversely, during the breaking stage, $a_{rmsk}$ decreases as $u_*/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_{rmsk}$ 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. (2023) 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. 2014).

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 2006;
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

Figure 11. Drag coefficient $C_D$ evaluated at $\bar{z} = 10$ m using (5.3) as a function of $u_*/c$ (a) and $a_{rmsk}$ (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.* (2022). 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_*/c = [0.4 - 0.5 - 0.7 - 0.9]$ at $Re_{*,\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.* (2020) up to $u_*/c \approx 0.71$ and from Curcic & Haus (2020) up to $u_*/c \approx 2.25$.
Veron 2015), 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_*/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_*/c$.

## 6. Conclusions

We used state-of-the-art two-phase-fluid DNS of wind-forced breaking waves to get new insights into the processes controlling momentum flux at high wind speed. We consider the ratio of wind friction velocity to wave phase speed, $u_*/c$, ranging from 0.3 to 0.9. The wave field is initialized with an amplitude below the breaking threshold, allowing the wind to be the primary driver of wave growth until the breaking conditions are reached and simulate the growing-breaking wave life cycle. We analyse the momentum flux, separating the growth and breaking stages, and treating them as two distinct subprocesses, and revisit drag coefficient analysis.

The momentum flux analysis underscores the dominance of the pressure force over the viscous force throughout all stages of wave evolution. When breaking occurs, the pressure force sharply decreases, and this reduction is compensated by a sudden acceleration of the flow field to conserve momentum. The remaining terms in the momentum budget, including a convective term accounting for the relative velocity between the airflow and waves, and a viscous force contribution, show smaller magnitudes and are less affected. The modulation of airflow caused by wave breaking results in an upshift in the streamwise velocity profile, corresponding to the drop in pressure force.

The relative strengths of these effects are quantified using the aerodynamic drag coefficient $C_{D,a}$, defined as the ratio of mean pressure force to mean velocity at a reference height within the wave-boundary layer ($z = \lambda/2$). During wave growth, $C_{D,a}$ increases with $u_*/c$, while during breaking, it decreases with $u_*/c$. When averaging over both growing and breaking regimes, we obtain a saturation of the drag coefficient $C_{D,a}$ at high enough $u_*/c$. A similar behaviour is obtained if one computes from the DNS results the drag coefficient commonly employed in air-sea interaction, $C_D$, which is based on wind friction velocity $u_*$ and the reference velocity at 10 m height, $U_{10}$. The trend in the DNS closely mirrors field observations and laboratory experiments at high wind speeds for $C_D$, and the magnitudes of the drag in both DNS and laboratory are in remarkable agreement. We observe a difference of the critical $u_*/c$ for which the drag starts to saturate, which might be attributed to the differences in set-up, narrow- or broad-banded wave field as well as averaging procedure.

Our simulations demonstrate that the drag coefficient, whether defined as $C_{D,a}$ or $C_D$, saturates when wave-breaking events are included. Notably, the production of droplets remains negligible due to the moderate $Bo$, indicating that breaking-induced dynamics – such as reduced pressure stress, loss of coherence between the airflow and wave field and airflow separation – is sufficient to cause drag saturation.

The present framework can now be used to explore the associated energy fluxes going into the water column, including wind-wave growth, energy dissipation by breaking and oceanic turbulence development, further quantifying the growing-breaking life cycle of ocean waves under strong wind forcing.

**Funding.** This work is supported by the National Science Foundation under grant 2318816 to L.D. (Physical Oceanography Program), the NASA Wind Vector Science Team, grant 80NSSC23K0983 to L.D. N.S. acknowledges the support of the fellowship High Meadows Environmental Institute Postdoctoral Teaching
Program. Computations were performed using the Stellar machine, granted by the Cooperative Institute for Earth System Modeling (CIMES), and managed by Princeton Research Computing. This includes the Princeton Institute for Computational Science and Engineering, the Office of Information Technology's High-Performance Computing Center and the Visualization Laboratory at Princeton University.

**Declaration of interests.** The authors report no conflict of interest.

## Appendix A. Procedure for converting Cartesian to wave-following coordinates

In this appendix, we detail the procedure to transform a generic field defined in a Cartesian coordinate system into a wave-following one. This transformation involves three steps:

(i) The simulation outputs different instantaneous two-dimensional fields, hereinafter termed $q = f(x, z)$, which has been already averaged along the spanwise $y$ direction. Next, this generic two-dimensional field $q$ is transformed from a Cartesian coordinate system into a wave-follower coordinate using a one-dimensional mapping function (Sullivan et al. 2000; Wu et al. 2022):

$$ \begin{bmatrix} x \\ z \end{bmatrix} = \begin{bmatrix} x(\xi, \zeta) \\ z(\xi, \zeta) \end{bmatrix} = \begin{bmatrix} \xi \\ \zeta + \bar{\eta}(\xi) \exp(-k|\zeta|) \end{bmatrix}, \qquad (A1) $$

where $\bar{\eta}(\xi)$ is the spanwise-averaged surface elevation function. For a sufficiently large value of $k|\zeta|$, $\exp(-k|\zeta|) \approx 0$ and the vertical Cartesian coordinate coincides with the vertical wave-following coordinate. Note that alternative mapping functions to (A1) have been proposed, but the sensitivity of the results to the choice of mapping is typically negligible (Sullivan et al. 2000).

(ii) Each wave-following two-dimensional field is averaged along the periodic direction $x$, including also the regions near the wave troughs.

(iii) The obtained vertical profile, i.e. $q = g(\zeta)$, is time-averaged over the cycles $G_1$ and $G_{2,a}$ as defined in figure 3(b). Note that during $G_1$ and $G_{2,a}$ ($G_{2,b}$ and $F$), the airflow is quasi-stationary, and employing time-averaging is an adequate approach to define turbulence statistics.

## Appendix B. Calculation of the momentum fluxes

We provide details on how we evaluate the momentum fluxes in (4.2), i.e. $\rho_a \phi_{c,x}$, $\tau_{p,x}$ and $\tau_{v,x}$. These fluxes at the boundaries are not zero only at the two-phase interface, since at the other faces of the computational box (shown in figure 1), they are zero due to periodicity along the horizontal directions and at the top boundary due to the imposed no-penetration/free-slip condition.

The calculation of the momentum fluxes is conducted in two steps. First, the interface $\Gamma$ is slightly shifted vertically by $\Delta_s = 0.1/k$, corresponding approximately to four grid points at a resolution $L_e = 10$. As noted by Wu et al. (2022), this step is a numerical strategy to ensure that the pressure and the velocity gradients are evaluated in the air domain only. As long as $\Delta_s$ stays within the viscous sublayer, its numerical value has a negligible impact on the evaluation of the different terms of the momentum budget. Next, the surface integrals in (4.1) and (4.2) are transformed into volume integrals using the Gauss theorem, e.g. $\int_{\Gamma} p \mathbf{n} \cdot \mathbf{e}_x dS = \int_{\Omega_a} \nabla \cdot (p \mathbf{e}_x) dV = \int_{\Omega_a} (\partial p / \partial x) dV$. This surface-to-volume integral allows evaluating the momentum fluxes over the air volume $\Omega_a$ rather than on $\Gamma$. Owing to the volume-of-fluid reconstruction, $\Omega_a$ remains a conserved quantity in our simulations, even when the wave field breaks, whereas $\Gamma$ does
<figure><img src="image_23.png" /><figcaption>Figure 12. Evolution of the normalized wave energy $E_W/E_{W,0}$ as a function of the dimensionless time $\omega t$ with $E_{W,0} = E_W(t=0)$ for $Re_{*,\lambda} = 107 - 214$ at $u_*/c = 0.9$. Note that these two cases share the same $Re_* = 720$ but different $(L_0 - h_W)/\lambda = 3.36 - 6.72$. The continuous coloured lines refer to $Le = 10$, the black symbols to $Le = 11$.</figcaption></figure>

not. Importantly, the volume and surface integral formulations yield the same numerical results when the wave field grows.

## Appendix C. Convergence studies between $Le = 10$ and $Le = 11$

We perform a grid convergence study for wave energy and momentum flux for two cases at $u_*/c = 0.9$ and $Re_{*,\lambda} = 214$ with $(L_0 - h_W)/\lambda = 3.36$ and $(L_0 - h_W)/\lambda = 6.72$ (four and eight waves per box size, respectively). Two levels of refinement are considered, $Le = 10$ and $Le = 11$. At $Le = 10$, we have 256 (128) and 512 (256) grid points per wavelength in the case of $(L_0 - h_W)/\lambda = 3.36$ ($(L_0 - h_W)/\lambda = 6.72$) (and twice as much at $Le = 11$). As we will see, the number of grid points per wave boundary layer (as defined in Mostert et al. (2022)) is enough to obtain grid convergence of energy dissipation due to breaking.

Figure 12 displays the normalized wave energy $E_W/E_{W,0}$ as a function of the dimensionless time $\omega t$. For $(L_0 - h_W)/\lambda = 3.36$, employing a maximum refinement level of $Le = 10$ guarantees grid-converged results across all wave field growing and breaking stages (results are identical at $Le = 11$). Moreover, this level of refinement proves adequate to achieve grid-independent results also for the momentum budget components, as shown in figure 13(a).

Increasing $(L_0 - h_W)/\lambda$ from 3.36 to 6.72, while keeping the other dimensionless parameters fixed, makes the grid requirement more demanding, as grid convergence still requires the same number of grid points per wavelength (or per wave boundary layer, which scales with the wavelength). In particular, during the first growing stage, the wave energy and the components of the momentum budgets are well captured at a resolution of $Le = 10$, as shown in figures 12 and 13(b). Conversely, a deviation with respect to
<figure><img src="image_24.png" /><figcaption>Figure 13. Grid-convergence study for the momentum budget in the streamwise direction as in (4.2) for (a) $(L_0 - h_W)/\lambda = 3.36$ and $Re_{*,\lambda} = 214$ and (b) $(L_0 - h_W)/\lambda = 6.72$ and $Re_{*,\lambda} = 107$ at $u_*/c = 0.9$. Both cases share the same $Re_* = 720$. The coloured curves refer to $Le = 10$; the symbols refer to resolution $Le = 11$. On the y-axis label $\mathcal{T}$ represents the variation in the instantaneous flow $\rho_a \partial U / \partial t$, the viscous stress $\tau_{v,x}$, the pressure stress $\tau_{p,x}$, the convective term $\rho_a \phi_{c,x}$ or the driving force $\Pi_f$. Each budget component is normalized by the total stress $\rho_a u_*^2 A_\Gamma$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 14. (a) Contributions of the momentum budget in the streamwise direction, as in (4.2), for $Re_{*,\lambda} = 53.5 - 107 - 214$ with $(L_0 - h_W)/\lambda = 3.36$. On the y-axis label of both panels, $\mathcal{T}$ represents the variation in the instantaneous flow $\rho_a \partial U / \partial t$, the viscous stress $\tau_{v,x}$ and the pressure stress $\tau_{p,x}$. Each budget component is normalized by the total stress $\rho_a u_*^2 A_\Gamma$. (b) Streamwise velocity profile normalized by the nominal friction, $\langle u_a^+ \rangle = \langle u_a \rangle / u_*$, as a function of vertical wave-following coordinate (in wall units) $\zeta^+ = \zeta u_*/v_a$ for $Re_{*,\lambda} = 53.5 - 107 - 214$ with $(L_0 - h_W)/\lambda = 3.36$ in coloured lines (continuous and dashed lines for the first and second growing cycle, $G_1$ and $G_{2,a}$, respectively). The dotted black lines are the linear fit. In blue colour, we report the corresponding cases pertaining to the stationary and flat wall at the same $Re_*$ (the equivalent $Re_{*,\lambda}$ is displayed in parenthesis). Note that the velocity profiles are almost indistinguishable in the case of a flat stationary wall at $Re_* = 360 - 720$, while that at $Re_* = 180$ (blue dotted lines) is slightly upshifted.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 15. Contributions of the momentum budget in the streamwise direction, as in (4.2), for the following cases: $Re_{*,\lambda} = 214$ with $(L_0 - h_W)/\lambda = 3.36$, $Re_{*,\lambda} = 107$ with $(L_0 - h_W)/\lambda = 6.72$ and the case with $Re_{*,\lambda} = 107$ with $(L_0 - h_W)/\lambda = 6.72$. On the y-axis label $\mathcal{T}$ represents the variation in the instantaneous flow $\rho_a \partial U / \partial t$, the viscous stress $\tau_{v,x}$ and the pressure stress $\tau_{p,x}$. The convective term $\rho_a \phi_{c,x}$ or the driving force $\Pi_f$ are omitted for clarity. Each budget component is normalized by the total stress $\rho_a u_*^2 A_\Gamma$.</figcaption></figure>

$Le = 11$ occurs during the breaking stage and the second growing stages of the wave field. Accordingly, fixing all the other parameters, cases with $(L_0 - h_W)/\lambda > 6.72$ requires $Le = 11$ to obtain grid-independent results.

## Appendix D. Sensitivity of the momentum flux to $Re_{*,\lambda}$

We now present the results of a sensitivity study to the turbulent air-side friction Reynolds number, $Re_{*,\lambda}$. We consider three cases at a fixed $u_*/c = 0.9$: $Re_{*,\lambda} = 53.5$ and $Re_{*,\lambda} = 107$ with $(L_0 - h_W)/\lambda = 3.36$, and $Re_{*,\lambda} = 107$ with $(L_0 - h_W)/\lambda = 6.72$, which correspond to a friction Reynolds number $Re_* = 180 - 360 - 720$, respectively.

We start comparing the three cases by inspecting the momentum budget in the streamwise direction, as displayed in figure 14(a). The budget clearly shows that the time-varying pressure flux $\tau_{p,x}$ has a smaller magnitude over the entire growing and breaking cycles as $Re_{*,\lambda}$ decreases. A larger viscous flux, $\tau_{v,x}$, compensates for this decrease.

The overall reduced pressure drag experienced by the wave field at lower $Re_{*,\lambda}$ modifies the streamwise velocity profile, as displayed in figure 14(b). As we increase $Re_{*,\lambda}$, the velocity profile has a consistent downshift compared with the flat and stationary wall scenario (blue lines in figure 14b), which is caused by the increased drag due to the pressure component (absent for a flat wall). When the wave breaks, in all the cases, there is an upshift of the velocity profile, and the associated drag reduction is solely induced by wave breaking. As visible in the figure, the upshift is largest for $Re_{*,\lambda} = 214$, and it slightly decreases for $Re_{*,\lambda} = 53.5 - 107$. Despite an unavoidable $Re_{*,\lambda}$ effect,
<figure><img src="image_27.png" /><figcaption>Figure 16. Streamwise velocity profile normalized by the nominal friction $\langle u_a \rangle/u_*$ in (a) physical or (b) wall units for the cases at $Re_{*,\lambda} \approx 107$ and $(L_0 - h_W)/\lambda = 3.36 - 6.72$. The solid black line in (b) represents the velocity profile of the flow over a solid wall at the same $Re_{*,\lambda} = 107$. Note that differently from the profiles in (a), the velocity profiles as a function of the viscous unit do not collapse since the normalization factor for the viscous unit, i.e. $v_a/u_*$, is independent of the number of waves per unit of the box size.</figcaption></figure>

this analysis confirms that the breaking-induced flow acceleration and consequent drag reduction is associated with the change in the wave steepness of the wave field and surface roughness and is not significantly affected by $Re_{*,\lambda}$.

## Appendix E. Sensitivity of the momentum flux on $(L_0 - h_W)/\lambda$

We present a sensitivity study on the ratio of the box size to the wavelength, i.e. $(L_0 - h_W)/\lambda$, or equivalently, the number of waves per box size. To this end, we run an additional case at $Re_{*,\lambda} = 107$ with $(L_0 - h_W)/\lambda = 6.72$ and $u_*/c = 0.9$. Given the larger $(L_0 - h_W)/\lambda = 6.72$, we set $Le = 11$, as discussed in Appendix C.

This additional case is compared with: (i) the case at different $(L_0 - h_W)/\lambda = 3.36$ but equal $Re_{*,\lambda} = 107$ and (ii) the case at different $Re_{*,\lambda} = 214$ and $(L_0 - h_W)/\lambda = 6.72$, but equal $Re_* = 720$.

Figure 15 displays the different components of the momentum flux as in (4.2). We can see that despite the different values of $(L_0 - h_W)/\lambda$, the cases with the same $Re_{*,\lambda} = 107$ display very similar momentum fluxes. Conversely, the case at $Re_{*,\lambda} = 214$ shows larger pressure and smaller viscous fluxes. These observations suggest two conclusions. First, while the ratio $(L_0 - h_W)/\lambda$ imposes stricter resolution criteria, its effect on the results is limited. Second, when comparing momentum fluxes across different cases, the relevant dimensionless parameter is $Re_{*,\lambda}$.

Figure 16 reports the velocity profiles for the two cases at $Re_{*,\lambda} \approx 107$ with $(L_0 - h_W)/\lambda = 3.36 - 6.72$. The velocity profiles in figure 16(a) of the two cases display a good collapse during the two considered growing cycles in the region immediate to the wave field up to $z/\lambda \approx 1$. Here, the profiles are similar since the two cases share the same $Re_{*,\lambda}$, which is the controlling parameter for the turbulent intensity near the wave field. Away from the wave-field region, i.e. $z/\lambda > 1$, the velocity profile pertaining to the same growing cycle diverges and the case with $(L_0 - h_W)/\lambda = 6.72$ is downshifted. Indeed, in
the airflow regions well above the wave field, the turbulent intensity is controlled by $Re_*$,
which is larger for the case at larger ($L_0 - h_W)/\lambda$.

**Figure 16(b)** reports the velocity profiles in semi-logarithmic form, which both display an upshift induced by wave breaking. The upshift is of similar magnitude for both cases, which further confirms that the breaking-induced drag reduction is mainly controlled by $u_*/c$, whereas ($L_0 - h_W)/\lambda$ has a limited effect.

REFERENCES

AYET, A. & CHAPRON, B. 2022 The dynamical coupling of wind-waves and atmospheric turbulence: a review of theoretical and phenomenological models. *Boundary-Layer Meteorol.* **183** (1), 1–33.

BANNER, M., BARTHELEMY, X., FEDELE, F., ALLIS, M., BENETAZZO, A., DIAS, F. & PEIRSON, W. 2014 Linking reduced breaking crest speeds to unsteady nonlinear water wave group behavior. *Phys. Rev. Lett.* **112** (11), 114502.

BANNER, M. & MELVILLE, W. 1976 On the separation of air flow over water waves. *J. Fluid Mech.* **77** (4), 825–842.

BELCHER, S. & HUNT, J. 1998 Turbulent flow over hills and waves. *Annu. Rev. Fluid Mech.* **30** (1), 507–538.

BELL, J., COLELLA, P. & GLAZ, H. 1989 A second-order projection method for the incompressible Navier–Stokes equations. *J. Comput. Phys.* **85** (2), 257–283.

BUCKLEY, M. & VERON, F. 2016 Structure of the airflow above surface waves. *J. Phys. Oceanogr.* **46** (5), 1377–1397.

BUCKLEY, M., VERON, F. & YOUSEFI, K. 2020 Surface viscous stress over wind-driven waves with intermittent airflow separation. *J. Fluid Mech.* **905**, A31.

BYE, J. & JENKINS, A. 2006 Drag coefficient reduction at very high wind speeds. *J. Geophys. Res. Oceans* **111** (C3).

CHEN, G., KHARIF, C., ZALESKI, S. & LI, J. 1999 Two-dimensional Navier–Stokes simulation of breaking waves. *Phys. Fluids* **11** (1), 121–133.

CHEN, S., PRICE, J., ZHAO, W., DONELAN, M. & WALSH, E. 2007 The CBLAST-hurricane program and the next-generation fully coupled atmosphere-wave-ocean models for hurricane research and prediction. *Bull. Am. Meteorol. Soc.* **88** (3), 311–317.

CIMARELLI, A., ROMOLI, F. & STALIO, E. 2023 On wind-wave interaction phenomena at low reynolds numbers. *J. Fluid Mech.* **956**, A13.

CURCIC, M. & HAUS, B. 2020 Revised estimates of ocean surface drag in strong winds. *Geophys. Res. Lett.* **47** (10), e2020GL087647.

DAVIS, J.R., THOMSON, J., HOUGHTON, I.A., DOYLE, J.D., KOMAROMI, W.A., FAIRALL, C.W., THOMPSON, E.J. & MOSKAITIS, J.R. 2023 Saturation of ocean surface wave slopes observed during hurricanes. *Geophys. Res. Lett.* **50** (16), e2023GL104139.

DEIKE, L. 2022 Mass transfer at the ocean–atmosphere interface: the role of wave breaking, droplets, and bubbles. *Annu. Rev. Fluid Mech.* **54** (1), 191–224.

DEIKE, L. & MELVILLE, W. 2018 Gas transfer by breaking waves. *Geophys. Res. Lett.* **45** (19), 10–482.

DEIKE, L., MELVILLE, W. & POPINET, S. 2016 Air entrainment and bubble statistics in breaking waves. *J. Fluid Mech.* **801**, 91–129.

DEIKE, L., POPINET, S. & MELVILLE, W. 2015 Capillary effects on wave breaking. *J. Fluid Mech.* **769**, 541–569.

DONELAN, M., BABANIN, A., YOUNG, I. & BANNER, M. 2006 Wave-follower field measurements of the wind-input spectral function. part ii: parameterization of the wind input. *J. Phys. Oceanogr.* **36** (8), 1672–1689.

DONELAN, M., HAUS, B., REUL, N., PLANT, W., STIASSNIE, M., GRABER, H., BROWN, O. & SALTZMAN, E. 2004 On the limiting aerodynamic roughness of the ocean in very strong winds. *Geophys. Res. Lett.* **31** (18), L18306.

DRAZEN, D., MELVILLE, W. & LENAIN, L. 2008 Inertial scaling of dissipation in unsteady breaking waves. *J. Fluid Mech.* **611**, 307–332.

EDSON, J., JAMPANA, V., WELLER, R., BIGORRE, S., PLUEDDEMANN, A., FAIRALL, C., MILLER, S., MAHRT, L., VICKERS, D. & HERSBACH, H. 2013 On the exchange of momentum over the open ocean. *J. Phys. Oceanogr.* **43** (8), 1589–1610.

FARSOIYA, P., LIU, Z., DAISS, A., FOX, R. & DEIKE, L. 2023 Role of viscosity in turbulent drop break-up. *J. Fluid Mech.* **972**, A11.

FUNKE, C., BUCKLEY, M., SCHULTZE, L., VERON, F., TIMMERMANS, M. & CARPENTER, J. 2021 Pressure fields in the airflow over wind-generated surface waves. *J. Phys. Oceanogr.* **51** (11), 3449–3460.
GENTILE, E., GRAY, S., BARLOW, J., LEWIS, H. & EDWARDS, J. 2021 The impact of atmosphere-ocean-wave coupling on the near-surface wind speed in forecasts of extratropical cyclones. *Boundary-Layer Meteorol.* **180** (1), 105–129.

GENTILE, E., GRAY, S. & LEWIS, H. 2022 The sensitivity of probabilistic convective-scale forecasts of an extratropical cyclone to atmosphere-ocean-wave coupling. *Q. J. R. Meteorol.* **148** (743), 685–710.

GHIGO, A., POPINET, S. & WACHS, A. 2021 A conservative finite volume cut-cell method on an adaptive cartesian tree grid for moving rigid bodies in incompressible flows.

IAFRATI, A. 2009 Numerical study of the effects of the breaking intensity on wave breaking flows. *J. Fluid Mech.* **622**, 371–411.

IAFRATI, A., DE VITA, F. & R., VERZICCO 2019 Effects of the wind on the breaking of modulated wave trains. *Eur. J. Mech. (B/Fluids)* **73**, 6–23.

JANSSEN, P. 2004 *The Interaction of Ocean Waves and Wind*. Cambridge University Press.

JEFFREYS, H. 1925 On the formation of water waves by wind. *Proc. R. Soc. Lond.* A **107** (742), 189–206.

JOHANSEN, H. & COLELLA, P. 1998 A cartesian grid embedded boundary method for poisson's equation on irregular domains. *J. Comput. Phys.* **147** (1), 60–85.

KOMORI, S., IWANO, K., TAKAGAKI, N., ONISHI, R., KUROSE, R., TAKAHASHI, K. & SUZUKI, N. 2018 Laboratory measurements of heat transfer and drag coefficients at extremely high wind speeds. *J. Phys. Oceanogr.* **48** (4), 959–974.

KUDRYAVTSEV, V., CHAPRON, B. & MAKIN, V. 2014 Impact of wind waves on the air-sea fluxes: a coupled model. *J. Geophys. Res.: Oceans* **119** (2), 1217–1236.

LAMARRE, E. & MELVILLE, W. 1991 Air entrainment and dissipation in breaking waves. *Nature* **351** (6326), 469–472.

LAMB, S. 1993 *Hydrodynamics*. 6th edn. Cambridge University Press.

LIN, M., MOENG, C., TSAI, W., SULLIVAN, P. & BELCHER, S. 2008 Direct numerical simulation of wind-wave generation processes. *J. Fluid Mech.* **616**, 1–30.

LU, M., YANG, Z., HE, G. & SHEN, L. 2024 Numerical investigation on the heat transfer in wind turbulence over breaking waves. *Phys. Rev. Fluids* **9** (8), 084606.

MANZELLA, E., HARA, T. & SULLIVAN, P. 2024 Reduction of drag coefficient due to misaligned wind-waves. *J. Geophys. Res.: Oceans* **129** (5), e2023JC020593.

MCWILLIAMS, J. 2016 Submesoscale currents in the ocean. *Proc. R. Soc. Lond. A: Math. Phys. Engng Sci.* **472**(2189), 20160117.

MCWILLIAMS, J., SULLIVAN, P. & MOENG, C.-H. 1997 Langmuir turbulence in the ocean. *J. Fluid Mech.* **334**, 1–30.

MELVILLE, W. 1996 The role of surface-wave breaking in air-sea interaction. *Annu. Rev. Fluid Mech.* **28** (1), 279–321.

MELVILLE, W. & RAPP, R. 1985 Momentum flux in breaking waves. *Nature* **317** (6037), 514–516.

MELVILLE, W., SHEAR, R. & VERON, F. 1998 Laboratory measurements of the generation and evolution of langmuir circulations. *J. Fluid Mech.* **364**, 31–58.

MONIN, A. & OBUKHOV, A. 1954 Basic laws of turbulent mixing in the surface layer of the atmosphere. *Contrib. Geophys. Inst. Acad. Sci. USSR* **151** (163), e187.

MOSTERT, W., POPINET, S. & DEIKE, L. 2022 High-resolution direct simulation of deep water breaking waves: transition to turbulence, bubbles and droplets production. *J. Fluid Mech.* **942**, A27.

PERRARD, S., RIVIÈRE, A., MOSTERT, W. & DEIKE, L. 2021 Bubble deformation by a turbulent flow. *J. Fluid Mech.* **920**, A15.

POPE, S. 2000 *Turbulent Flows*. Cambridge University Press.

POPINET, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. *J. Comput. Phys.* **228** (16), 5838–5866.

POPINET, S. 2015 A quadtree-adaptive multigrid solver for the serre-green-naghdi equations. *J. Comput. Phys.* **302**, 336–358.

POPINET, S. 2018 Numerical models of surface tension. *Annu. Rev. Fluid Mech.* **50** (1), 49–75.

REUL, N., BRANGER, H. & GIOVANANGELI, J.-P. 1999 Air-flow separation above unsteady breaking waves. *Phys. Fluids* **648**, 1959–1961.

RIVIÈRE, A., MOSTERT, W., PERRARD, S. & DEIKE, L. 2021 Sub-hinze scale bubble production in turbulent bubble break-up. *J. Fluid Mech.* **917**, A40.

SRAJ, I., ISKANDARANI, M., SRINIVASAN, A., THACKER, W., WINOKUR, J., ALEXANDERIAN, A., LEE, C.-Y., CHEN, S. & KNIO, O. 2013 Bayesian inference of drag parameters using axbt data from typhoon fanapi. *Mon. Weath. Rev.* **141** (7), 2347–2367.

SROKA, S. & EMANUEL, K. 2021 A review of parameterizations for enthalpy and momentum fluxes from sea spray in tropical cyclones. *J. Phys. Oceanogr.* **51** (10), 3053–3069.
SULLIVAN, P. & MCWILLIAMS, J. 2010 Dynamics of winds and currents coupled to surface waves. *Annu. Rev. Fluid Mech.* **42** (1), 19–42.

SULLIVAN, P., MCWILLIAMS, J. & MOENG, C. 2000 Simulation of turbulent flow over idealized water waves. *J. Fluid Mech.* **404**, 47–85.

SULLIVAN, P.P., BANNER, M.L., MORISON, R.P. & PEIRSON, W.L. 2018 Turbulent flow over steep steady and unsteady waves under strong wind forcing. *J. Phys. Oceanogr.* **48** (1), 3–27.

TAKAGAKI, N., KOMORI, S., SUZUKI, N., IWANO, K. & KUROSE, R. 2016 Mechanism of drag coefficient saturation at strong wind speeds. *Geophys. Res. Lett.* **43** (18), 9829–9835.

TRYGGVASON, G., SCARDOVELLI, R. & ZALESKI, S. 2011 *Direct Numerical Simulations of Gas–liquid Multiphase Flows*. Cambridge University Press.

VAN HOOFT, J., POPINET, S., VAN HEERWAARDEN, C., VAN DER LINDEN, S., DE ROODE, S. & VAN DE WIEL, B. 2018 Towards adaptive grids for atmospheric boundary-layer simulations. *Boundary-Layer Meteorol.* **167** (3), 421–443.

VERON, F. 2015 Ocean spray. *Annu. Rev. Fluid Mech.* **47** (1), 507–538.

VERON, F. & MELVILLE, W. 2001 Experiments on the stability and transition of wind-driven water surfaces. *J. Fluid Mech.* **446**, 25–65.

VERON, F., SAXENA, G. & MISRA, S. 2007 Measurements of the viscous tangential stress in the airflow above wind waves. *Geophys. Res. Lett.* **34** (19), L19603.

WU, J. & DEIKE, L. 2021 Wind wave growth in the viscous regime. *Phys. Rev. Fluids* **6** (9), 094801.

WU, J., POPINET, S. & DEIKE, L. 2022 Revisiting wind wave growth with fully coupled direct numerical simulations. *J. Fluid Mech.* **951**, A18.

YANG, D. & SHEN, L. 2010 Direct-simulation-based study of turbulent flow over various waving boundaries. *J. Fluid Mech.* **650**, 131–180.

YANG, Z., DENG, B.-Q. & SHEN, L. 2018 Direct numerical simulation of wind turbulence over breaking waves. *J. Fluid Mech.* **850**, 120–155.

YOUSEFI, K., VERON, F. & BUCKLEY, M. 2020 Momentum flux measurements in the airflow over wind-generated surface waves. *J. Fluid Mech.* **895**, A15.

ZHANG, J., HECTOR, A., RABAUD, M. & MOISY, F. 2023 Wind-wave growth over a viscous liquid. *Phys. Rev. Fluids* **8** (10), 104801.