4.1 Validation

Figure 2.

Simulation results for wave breaking with $\unicode[STIX]{x1D716}=0.35$ and $u_{\ast }=0$ . (a) Wave geometry from $t=0$ to $3.0T$ . The wave surface is shifted by $ct/\unicode[STIX]{x1D706}$ in the $+z$ -direction. (b) Time series of total energy $E$ .

Figure 2(a) shows the surfaces of breaking waves with $\unicode[STIX]{x1D716}=0.35$ without wind, from $t=0$ to $3.0T$ . The wave surface at $t$ is shifted by $ct/\unicode[STIX]{x1D706}$ in the $+z$ -direction. The time interval between two wave surfaces is $0.1T$ . As shown, the breaking is gentle, and the wave amplitude decreases slowly. This observation is consistent with the results of Iafrati (Reference Iafrati2009). Figure 2(b) further compares the total energy $E/E_{0}$ with the results of Iafrati (Reference Iafrati2009) and Deike et al. (Reference Deike, Popinet and Melville2015). Here, $E$ is the mechanical energy, defined as the summation of the potential energy $E_{p}$ and the kinetic energy $E_{k}$ of water, and $E_{0}$ is the value of $E$ at $t=0$ . The definitions of $E_{p}$ and $E_{k}$ are given respectively as

(4.1) $$\begin{eqnarray}\displaystyle E_{p}=\int \unicode[STIX]{x1D70C}(z+d/2)\,\text{d}V, & & \displaystyle\end{eqnarray}$$

and

(4.2) $$\begin{eqnarray}\displaystyle E_{k}=\frac{1}{2}\int \unicode[STIX]{x1D70C}(u^{2}+v^{2}+w^{2})\,\text{d}V. & & \displaystyle\end{eqnarray}$$

In (4.1), $d$ is the water depth (see § 2). The integrations in (4.1) and (4.2) are performed in the water. Figure 2(b) shows that the present results agree with those in the literature. The test results on the accuracy of our code in simulating non-breaking and breaking waves are also given in Hu et al. (Reference Hu, Guo, Lu, Liu, Dalrymple and Shen2012) and Liu (Reference Liu2012).

4.2 Resolution study

Figure 3.

Contours of the ratio between the grid resolution $\unicode[STIX]{x1D6E5}$ and Kolmogorov length scale $\unicode[STIX]{x1D702}_{K}$ at $t=0.8T$ in case A55. The thick dash-dotted line represents the upper bound of wave surfaces.

We have examined the grid resolution to confirm that it is sufficient for performing DNS of wind over breaking waves. Figure 3 displays the contours of the ratio between the grid size $\unicode[STIX]{x1D6E5}=(\unicode[STIX]{x1D6E5}_{x}\unicode[STIX]{x1D6E5}_{y}\unicode[STIX]{x1D6E5}_{z})^{1/3}$ and Kolmogorov length scale $\unicode[STIX]{x1D702}_{K}=(\unicode[STIX]{x1D708}_{a}^{3}/\unicode[STIX]{x1D716})^{1/4}$ at $t=0.8T$ in case A55. Here, $\unicode[STIX]{x1D716}=\unicode[STIX]{x1D708}\overline{\unicode[STIX]{x1D622}_{ij}^{\prime \prime }\unicode[STIX]{x1D622}_{ij}^{\prime \prime }}$ is the dissipation rate of turbulence kinetic energy (TKE), with $\unicode[STIX]{x1D622}_{ij}^{\prime \prime }=\unicode[STIX]{x2202}u_{i}^{\prime \prime }/\unicode[STIX]{x2202}x_{j}$ being the fluctuation of the velocity gradient tensor. It is demonstrated in figure 3 that the value of $\unicode[STIX]{x1D6E5}/\unicode[STIX]{x1D702}_{K}$ is smaller than $1.4$ in the bulk of air. We have also examined the grid resolution at different instants (from $t=0$ to $t=3.0T$ with an interval of $0.1T$ ) in all simulation cases listed in table 1. The maximum of $\unicode[STIX]{x1D6E5}/\unicode[STIX]{x1D702}_{K}$ is smaller than 1.6, indicating that the grid resolution of the DNS is satisfied in the present study (Eswaran & Pope Reference Eswaran and Pope1988). Because we focus on the airflow in the present study, we show only the resolution in the air phase. In fact, the grid resolution normalized by the Kolmogorov scale is more adequate in water than in air due to the lower Reynolds number in water.

We note here that the instantaneous wave surface varies with the spanwise coordinate $y$ . Moreover, it differs among different ensemble runs. Therefore, the air–water interfaces should not be represented deterministically by a single line in a 2-D plot of ensemble-averaged statistics. In this paper, we use a thick dash-dotted line to denote the upper bound of all interfaces (figure 3). There are different approaches for defining the upper bound. The most strict definition corresponds to $\bar{\unicode[STIX]{x1D713}}=0$ , above which no water is found in any spanwise plane and in any sample. However, as pointed out by Brocchini & Peregrine (Reference Brocchini and Peregrine2001), because some small water droplets can be convected far away from the wave surface, it is more practical to use $\bar{\unicode[STIX]{x1D713}}=0.05$ to define the upper bound, which is adopted in the present study. In the region between the upper and lower bounds of the interfaces, the mean density is the fractional-averaged value of the densities of air and water, which varies spatially. In this sense, the flow is statistically compressible in this region (Hong & Walker Reference Hong and Walker2000; Chassaing et al. Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002; Hendrickson et al. Reference Hendrickson, Weymouth, Banerjee and Yue2013). Because we focus on the airflow in the present study, we only consider the flow field above the upper bound to reduce the complexity of the transport equations of the mean velocity and TKE.

4.3 Instantaneous wave surface

The process of wave breaking has been described in many papers (e.g. Duncan Reference Duncan1981, Reference Duncan1983; Longuet-Higgins Reference Longuet-Higgins1982; Bonmarin Reference Bonmarin1989; Chen et al. Reference Chen, Kharif, Zaleski and Li1999; Tulin & Waseda Reference Tulin and Waseda1999; Iafrati Reference Iafrati2009; Deike et al. Reference Deike, Popinet and Melville2015). The geometries of plunging and spilling breakers obtained from the present simulations are consistent with those observed in previous numerical simulations (Chen et al. Reference Chen, Kharif, Zaleski and Li1999; Iafrati Reference Iafrati2009; Deike et al. Reference Deike, Popinet and Melville2015; Lubin & Glockner Reference Lubin and Glockner2015). We demonstrate representative 3-D wave surfaces at different stages of wave breaking in this section. The plunging and spilling breakers are described based on cases A55 and A35, respectively.

Figure 4.

Successive snapshots of instantaneous wave surface of plunging breakers in case A55. The waves propagate from left to right.

Figure 4 shows successive snapshots of instantaneous wave surfaces visualized using the isosurface of $\unicode[STIX]{x1D719}=0$ in case A55. The initial geometry of the wave surface is plotted in figure 4(a). The wave front is steepened and becomes almost vertical at $t=0.3T$ (figure 4 b). An overturning jet is then formed at the wave front at $t=0.8T$ (figure 4 c). The jet impinges onto the wave surface, causing splash-up events at $t=1.0T$ (figure 4 d). At $t=3.0T$ (figure 4 e), the waves keep breaking but the wave steepness is much smaller than that before the wave breaking (figure 4 a).

Figure 5.

Successive snapshots of instantaneous wave surface of spilling breakers in case A35. The waves propagate from left to right.

Figure 5 shows the characteristic snapshots of spilling breaking waves in case A35. Different from case A55, no overturning jet is generated in case A35. Instead, as displayed in figure 5(b), a bulge is formed when the wave front is steepened. Figure 5(c,d) shows that the bulge grows during wave breaking, and the wave surface is roughened behind (to the left in the figure) the wave crest.

4.4 Mean velocity

Figure 6.

Vertical profiles of plane-averaged mean streamwise velocity $\langle u\rangle /u_{\ast }$ in cases (a) A35, (b) A55, (c) B55 and (d) C55. The results of cases A35, B35 and C35 at $t=0.0T$ are plotted respectively in (b–d) for comparison.

Figure 6(a) compares the vertical profiles of the plane-averaged mean streamwise velocity $\langle u\rangle$ normalized by $u_{\ast }$ at different time instants in case A35. Figure 6(b–d) shows respectively the vertical profiles of $\langle u\rangle /u_{\ast }$ in cases A55, B55 and C55, with the results in the respective cases A35, B35 and C35 at $t=0.0T$ superposed for comparison. Because the plane averaging is undefined below the wave crest, the lowest point of the profile is away from $z=0$ , which is different from the analyses of turbulence statistics on wave-surface-fitted grid (Hara & Sullivan Reference Hara and Sullivan2015). As shown in figure 6(a), the profile of $\langle u\rangle$ changes little during the spilling breaking process in case A35. Similarly, spilling breakers do not alter the profiles of $\langle u\rangle$ much in cases B35 and C35, so that we only show results at $t=0.0T$ in these two cases to keep the paper concise.

Figure 6(b) shows that at $t=0.0T$ , the magnitude of $\langle u\rangle$ below $z/\unicode[STIX]{x1D706}=0.4$ is slightly larger in case A35 than in case A55, but it is less sensitive to the wave steepness above $z/\unicode[STIX]{x1D706}=0.4$ . In the intermediate and old wave cases, $\langle u\rangle$ is more sensitive to the wave steepness. As demonstrated in figure 6(c), at $t=0.0T$ , $\langle u\rangle$ below $z/\unicode[STIX]{x1D706}=0.8$ is larger in case B35 than in case B55. Figure 6(d) shows that the profiles of $\langle u\rangle$ in cases C35 and C55 are drastically different. The vertical gradient of $\langle u\rangle$ is positive in case C35, but it is negative below $z/\unicode[STIX]{x1D706}=0.3$ in case C55. The negative velocity gradient near the wave surface is consistent with previous simulation results of wind over non-breaking waves at large wave age (Jiang et al. Reference Jiang, Sullivan, Wang, Doyle and Vincent2016). A further discussion on the negative velocity gradient in case C55 is given in § 4.5.

Near the wave surface, the profiles of $\langle u\rangle$ in group I cases are sensitive to the violent plunging wave breaking. As demonstrated in figure 6(b–d), the magnitude of $\langle u\rangle$ near the wave surface is larger at $t=3.0T$ than at $t=0.0T$ in all of the three group I cases, indicating the overall transport of momentum from water to air during the wave plunging at various wave ages. However, it is observed in figure 6(b) that in case A55, $\langle u\rangle$ experiences a transient decrease above $z/\unicode[STIX]{x1D706}=0.2$ before $t=1.0T$ . Figure 6(c,d) shows respectively that in cases B55 and C55, the magnitude of $\langle u\rangle$ decreases slightly below $z/\unicode[STIX]{x1D706}=0.2$ after $t=2.0T$ . A further analysis on the response of $\langle u\rangle$ to the wave plunging is given in § 4.5 together with the total momentum flux $\langle u^{\prime }w^{\prime }\rangle$ .

Figure 7.

Phase-averaged flow field at (a) $t=0.0T$ , (b) $t=0.7T$ , (c) $t=0.8T$ , (d) $t=1.0T$ and (e) $t=3.0T$ in case A55. Mean streamwise velocity $\bar{u}$ , mean vertical velocity $\bar{w}$ , and mean spanwise vorticity $\unicode[STIX]{x1D714}_{y}$ are shown in panels 1, 2 and 3, respectively. The dotted line demonstrates an instantaneous air–water interface in a spanwise plane. The dash-dotted line represents the upper bound of all wave surfaces (see § 4.2 for the definition). The solid and dashed lines in panel 3 represent positive and negative isopleths of $\unicode[STIX]{x1D714}_{y}\unicode[STIX]{x1D706}/c$ , respectively. The interval of the isopleths is 5. The wave propagates in the $+x$ -direction.

Figure 7 shows the ensemble-averaged mean streamwise velocity $\bar{u}$ , mean vertical velocity $\bar{w}$ , and mean spanwise vorticity $\bar{\unicode[STIX]{x1D714}}_{y}$ near the wave surface at successive time instants in case A55. Note that the flow physics far away from the wave surface are dominated by wind, so that we use the characteristic velocity scale of wind $u_{\ast }$ to non-dimensionalize $\langle u\rangle$ in figure 6. In contrast, the effects of waves are dominant near the wave surface, and as a result, we use $c$ as the characteristic velocity scale in figure 7. As introduced in § 4.2, the thick dash-dotted line represents the upper bound of wave surfaces in different spanwise planes and simulation samples, which is different from instantaneous wave surface. For comparison, we also plot an instantaneous air–water interface at an $x$ – $z$ plane using the dotted line in figure 7. At $t=0.0T$ , the instantaneous wave surface collapses to the upper bound of interfaces, such that only the upper bound is shown in figure 7(a). It can be observed from figure 7(b–e) that at later times, the dotted lines are below the dash-dotted lines, which is expected based on the definition of the upper bound of interfaces given in § 4.2. Also note that as explained in § 3, ensemble-averaged quantities are defined in the region $[0,\unicode[STIX]{x1D706}]$ in the streamwise direction. Due to the use of the periodic boundary condition in the streamwise direction, we show the results in the region $[0,2\unicode[STIX]{x1D706}]$ by duplicating the flow field in $[0,\unicode[STIX]{x1D706}]$ to $[\unicode[STIX]{x1D706},2\unicode[STIX]{x1D706}]$ . This ensures that all flow structures can be demonstrated in a continuous frame.

As demonstrated in figure 7(a2), when the wave is onset to break, the vertical velocity is negative and positive on the windward and leeward of the wave crest, respectively. This can be explained by the fact that the vertical component of the orbital velocity at the wave surface is positive on the leeward, and is negative on the windward. Therefore, the wave pushes the air on the leeward to move upward, while the air on the windward moves downward to follow the motion of the wave surface.

Figure 7(a3) shows that at $t=0.0T$ , the vorticity near the leeward of the wave surface is of small magnitude, and a shear layer detached from the wave surface appears in the downstream of the wave crest. This indicates that air particles do not move downward along the wave surface when they flow past the wave crest. Instead, they move upward, and a separation region is developed on the leeward of the wave surface. The airflow separation is consistent with observations from laboratory experiments (Banner & Melville Reference Banner and Melville1976; Reul et al. Reference Reul, Branger and Giovanangeli1999; Veron et al. Reference Veron, Saxena and Misra2007; Reul et al. Reference Reul, Branger and Giovanangeli2008; Buckley & Veron Reference Buckley and Veron2016).

The overturning jet (figure 7 b,c) formed during the first plunging event alters the flow structures significantly. When the jet impinges onto the wave surface, it pushes the air to move forward. This process accelerates the airflow in the streamwise direction, and as a result, a high-speed region is developed ahead of the wave front (figure 7 b1,c1). Meanwhile, the jet brings the air particles around it to impinge onto the wave surface. Therefore, the vertical velocity of airflow displayed in figure 7(b2,c2) becomes negative on the leeward of the wave crest, and the airflow separation over the wave crest is weakened. This is evident from figure 7(b3,c3), which shows that the detached shear layer downstream of the wave crest is significantly weakened. Furthermore, figure 7(b2) shows that the upward-moving air on the leeward is pushed away from the wave surface, and a counterclockwise-rotating vortex is generated between the downward and upward airflow. This vortex corresponds to the negative vorticity in figure 7(b3,c3). The airflow above the vortex is decelerated to compensate the air moving downward with the jet, such that a low-speed region appears above the vortex (figure 7 b1,c1).

From the comparison between figure 7(a1,e1), it can be seen that the streamwise velocity above the wave crest increases during the wave breaking in case A55. This observation is consistent with the plane-averaged mean velocity $\langle u\rangle$ shown in figure 6(b). Two physical processes responsible for the acceleration of wind can be observed from figure 7. The first process corresponds to the air motion around the large-scale vortex. As illustrated in figure 7(c), when the airflow goes along the cycling routine ABCDA, the low-speed air moves downward with the jet from A to B, and is then accelerated from B to C. The high-speed air moves upward from C to D, losing momentum to the flow above the wave crest. Finally, the low-speed air moves back to point A (note that the mean streamwise velocity is positive in the earth frame plotted, but is negative in the frame travelling with the wave) and starts a new cycle. In this process, the vortex transports the high-speed air from the region ahead of the wave front to that above the wave crest. The second process lies in the fact that the velocity of water near the wave crest increases during wave breaking. Correspondingly, the air velocity above the wave crest is also accelerated. This is evident from figure 7(b1,c1), which shows that a high-speed region occurs over the wave crest. The two mechanisms are both associated with the first plunging event. As a result, the plane-averaged velocity near the wave crest increases more significantly before $t=1.0T$ than after that. For example, figure 6(b) shows that the value of $\langle u\rangle /u_{\ast }$ at $z/\unicode[STIX]{x1D706}=0.1$ increases by 48 % before $t=1.0T$ , whereas it increases by only 12 % after that.

Figure 8.

Phase-averaged flow field at (a) $t=0.0T$ , (b) $t=0.8T$ and (c) $t=3.0T$ in case B55. See figure 7 for the detailed caption.

Figure 8 shows the ensemble-averaged mean streamwise velocity $\bar{u}$ , mean vertical velocity $\bar{w}$ and mean spanwise vorticity $\bar{\unicode[STIX]{x1D714}}_{y}$ in case B55 for the wave with intermediate wave age (table 1). Figure 8(a1) does not show a shear layer over the wave crest in case B55 before the wave breaking, indicating that airflow separation is not evident at larger wave age. This is consistent with the observation in the laboratory experiment of Buckley & Veron (Reference Buckley and Veron2016). Similar to case A55, a spanwise vortex with negative $\bar{\unicode[STIX]{x1D714}}_{y}$ is generated during the first plunging event (figure 8 b3). From the comparison between figure 8(a1,c1), it is known that the mean velocity at the late stage of wave breaking ( $t=3.0T$ ) is slightly larger than that before wave breaking ( $t=0.0T$ ). Such an acceleration effect of wave breaking on the mean streamwise velocity is similar to that at small wave age in case A55.

Figure 9.

Phase-averaged flow field at (a) $t=0.0T$ , (b) $t=0.8T$ and (c) $t=3.0T$ in case C55. See figure 7 for the detailed caption.

Figure 9 shows the ensemble-averaged flow field in case C55 with large wave age (table 1). Similar to case B55 (figure 8 a3), airflow separation is not evident in case C55 (figure 9 a3). As shown in figure 9(b3), the wave breaking also generates a spanwise vortex in case C55. From the comparison among figures 7(c1), 8(b1) and 9(b1), it can be seen that the streamwise scale of the high-speed region is larger in case C55 than in cases A55 and B55. The vortex also expands in the streamwise direction in case C55. It is known from the comparison between figure 9(a1,c1) that a plunging beaker with large wave age also causes the increase of the air velocity. This is consistent with the observation from figures 7 and 8 for cases A55 and B55, respectively, and is also consistent with the result of the plane-averaged mean velocity $\langle u\rangle$ shown in figure 6(d).

To further study the effect of wave breaking on the mean velocity, we analyse the transport equation of the integral mean velocity, which reads

(4.3) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\int _{\unicode[STIX]{x1D6FA}}\overline{u}\,\text{d}V & = & \displaystyle \underbrace{-\int _{TB}\overline{uu_{ri}}n_{i}\,\text{d}s}_{C_{TB}}\underbrace{-\int _{TB}\frac{\bar{p}}{\unicode[STIX]{x1D70C}_{a}}n_{1}\,\text{d}s}_{\unicode[STIX]{x1D6F1}_{TB}}+\underbrace{\int _{TB}\unicode[STIX]{x1D708}_{a}\frac{\unicode[STIX]{x2202}\bar{u}}{\unicode[STIX]{x2202}x_{i}}n_{i}\,\text{d}s}_{D_{TB}}\nonumber\\ \displaystyle & & \displaystyle +\,\underbrace{\unicode[STIX]{x1D706}\langle uw\rangle |_{z_{1}}}_{C_{1}}+\underbrace{\unicode[STIX]{x1D706}\unicode[STIX]{x1D708}_{a}\left.\frac{\unicode[STIX]{x2202}\langle u\rangle }{\unicode[STIX]{x2202}z}\right|_{z_{1}}}_{D_{1}}.\end{eqnarray}$$

The equation above is derived by performing ensemble averaging on (2.2), and integrating the resultant equation in a control volume $\unicode[STIX]{x1D6FA}$ between the top bound (denoted by TB in the equation) of the wave surfaces and the horizontal plane at $z_{1}=0.14\unicode[STIX]{x1D706}$ , which is chosen to ensure that the wave crest is close but below it. It is shown in the following analyses that the evolution of the mean streamwise velocity near the wave surface is dominated by the momentum fluxes across the wave surface. Changing the value of $z_{1}$ does not alter our conclusions. Note that according to the definition of the ensemble averaging (3.2), $\unicode[STIX]{x1D6FA}$ is a 2-D geometry, while the top bound of wave surfaces and horizontal plane reduce to lines in (4.3). The left-hand side of (4.3) is the material derivative of the integral mean velocity. The first, second and third terms on the right-hand side of (4.3) are respectively the convection flux $C_{TB}$ , pressure flux $\unicode[STIX]{x1D6F1}_{TB}$ and viscous flux $D_{TB}$ across the top bound of wave surfaces, where $\boldsymbol{u}_{r}=[u_{r1},u_{r2},u_{r3}]$ is the relative velocity of airflow with respect to the moving speed of the top bound of wave surfaces, with $\boldsymbol{n}=[n_{1},n_{2},n_{3}]$ being its normal direction, pointing from air into water. Because the top bound of wave surfaces is not a material surface, fluid particles can cross it, resulting in a relative velocity between the airflow and top bound. However, this relative velocity is small, and as a result, the convection flux does not play a major role at the top bound of the wave surface, which is evident from the results shown below. The fourth and fifth terms on the right-hand side of (4.3) are respectively the convection flux $C_{1}$ and viscous flux $D_{1}$ across plane $z_{1}$ . Note that the definition of the plane averaging (3.5) is used to derive the last two terms.

Figure 10.

Time history of momentum fluxes across the top bound of wave surfaces and horizontal plane at $z_{1}$ in cases (a) A35, (b) A55, (c) B55 and (d) C55.

Figure 11.

Phase-averaged pressure field near the wave surface in case A55 at (a) $t=0.0T$ and (b) $t=0.6T$ . The pressure is normalized by $\unicode[STIX]{x1D70C}_{a}c^{2}$ . The thick dashed-dotted line represents the top bound of wave surfaces (see § 4.2 for the definition). The waves propagate in the $+x$ -direction.

Figure 10 compares the time series of the flux terms on the right-hand side of equation (4.3) in cases A35, A55, B55 and C55. The total flux is also shown in the figure. All fluxes in figure 10 are normalized using $\unicode[STIX]{x1D706}c^{2}$ . At $t=0.0T$ , the initial field is a fully developed turbulent flow, so that the total flux is negligibly small, i.e. the flux terms are balanced, in all cases. In cases A35 and A55 (figure 10 a,b) for young waves, the pressure flux at the top bound of wave surface $\unicode[STIX]{x1D6F1}_{TB}$ is negative. In case B55 for the intermediate waves (figure 10 c), the magnitude of $\unicode[STIX]{x1D6F1}_{TB}$ becomes significantly smaller. In case C55 (figure 10 d) for the old waves, $\unicode[STIX]{x1D6F1}_{TB}$ reverses its sign to be positive. It has also been reported in the literature of wind over non-breaking waves (Sullivan et al. Reference Sullivan, McWilliams and Moeng2000; Yang & Shen Reference Yang and Shen2010) that the value of pressure flux at the wave surface increases as the wave age increases.

Figure 10(a) shows that all fluxes change little during wave breaking in case A35 with spilling breakers. As a result, the mean streamwise velocity $\langle u\rangle$ shown in figure 6(a) does not change much. In contrast, with plunging breakers, the total flux increases rapidly before $t=0.7T$ in cases A55, B55 and C55. The pressure flux is the dominant term that causes the increase of the total flux. As shown in figure 10(b–d), the pressure flux increases significantly after $t=0.3T$ . To further explore the physical process responsible for the increase of the pressure flux across the wave surface, we compare the ensemble-averaged pressure at $t=0.0T$ and $t=0.6T$ in case A55 in figure 11. It can be observed that at $t=0.0T$ , the pressure on the windward is higher than that on the leeward, which results in a form drag on the wind over steep waves. When the overturning jet is formed at $t=0.6T$ , as shown in figure 11(b), the pressure on the leeward increases significantly. This can be explained by the fact that the overturning jet increases the air pressure significantly to push the downstream air to move forward (see also figure 7 b1). The pressure field demonstrated in figure 11(b) explains why the pressure flux changes its sign to become positive during wave plunging.

At the horizontal plane $z_{1}$ , both the convection flux $C_{1}$ and viscous flux $D_{1}$ are positive in case A55 (figure 10 b), indicating that the air near the wave surface gains momentum from the air above this plane. The magnitude of $D_{1}$ decreases after the wave breaking, consistent with the observation from figure 6(b) that the vertical gradient of mean velocity $\unicode[STIX]{x2202}\langle u\rangle /\unicode[STIX]{x2202}z$ decreases during wave breaking in case A55. The magnitude of $C_{1}$ also decreases, indicating that the turbulence is suppressed there. Figure 10(c,d) shows respectively that in cases B55 and C55, $D_{1}$ is relatively small compared with the other fluxes. The value of $C_{1}$ is positive before wave breaking, characterizing a downward momentum flux. However, $C_{1}$ becomes negative at $t=0.8T$ and $t=0.3T$ in cases B55 and C55, respectively, indicating that the mean momentum gained from the pressure flux at the top bound of wave surfaces is transported to higher altitudes.

4.5 Momentum flux

Figure 12.

Vertical profiles of (a) plane-averaged total momentum flux $-\langle u^{\prime }w^{\prime }\rangle /u_{\ast }^{2}$ and (b) wave-coherent momentum flux $-\langle \tilde{u} \tilde{w}\rangle /u_{\ast }^{2}$ in case A55 at $t=0.0T$ , $1.0T$ , $2.0T$ and $3.0T$ . The result of case A35 at $t=0.0T$ is plotted for comparison.

The plane-averaged total momentum flux $-\langle u^{\prime }w^{\prime }\rangle$ plays an important role in the evolution of the plane-averaged mean streamwise velocity, of which the governing equation reads

(4.4) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\langle u\rangle }{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}\langle u\rangle }{\unicode[STIX]{x2202}z^{2}}-\frac{\unicode[STIX]{x2202}\langle u^{\prime }w^{\prime }\rangle }{\unicode[STIX]{x2202}z}. & & \displaystyle\end{eqnarray}$$

Figure 12 compares the profiles of $-\langle u^{\prime }w^{\prime }\rangle$ and its wave-coherent part $-\langle \tilde{u} \tilde{w}\rangle$ in case A55. The result of case A35 at $t=0.0T$ is superimposed for comparison. At $t=0.0T$ , the profiles of $-\langle u^{\prime }w^{\prime }\rangle$ in cases A35 and A55 are close to each other, especially in the region away from the wave surface. As introduced in § 3, the initial field of the simulation is the fully developed wind turbulence over prescribed waves. In this statistically steady flow field, the left-hand side of (4.4) is trivial, and as a result, equation (4.4) can be further simplified to

(4.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}\langle u\rangle }{\unicode[STIX]{x2202}z}-\langle u^{\prime }w^{\prime }\rangle =u_{\ast }^{2}. & & \displaystyle\end{eqnarray}$$

Figure 6(b) shows that $\unicode[STIX]{x2202}\langle u\rangle /\unicode[STIX]{x2202}z$ decreases as $z/\unicode[STIX]{x1D706}$ increases and becomes negligibly small far away from the wave surface. Therefore, as expected, it is observed from figure 12(a) that the value of $-\langle u^{\prime }w^{\prime }\rangle /u_{\ast }^{2}$ at $t=0.0T$ increases and approaches 1.0 as $z/\unicode[STIX]{x1D706}$ increases. This trend is consistent with the results of wind over non-breaking young waves in the literature (e.g. Sullivan et al. Reference Sullivan, McWilliams and Moeng2000; Yang & Shen Reference Yang and Shen2010; Druzhinin et al. Reference Druzhinin, Troitskaya and Zilitinkevich2012; Sullivan et al. Reference Sullivan, McWilliams and Patton2014).

The first plunging event disturbs the air near the wave surface, in a peak in the profile of $-\langle u^{\prime }w^{\prime }\rangle$ near $z/\unicode[STIX]{x1D706}=0.27$ at $t=1.0T$ (figure 12 a). Below and above this peak, the value of $-\unicode[STIX]{x2202}\langle u^{\prime }w^{\prime }\rangle /\unicode[STIX]{x2202}z$ is positive and negative, respectively. Because the viscous shear stress is less important than $-\langle u^{\prime }w^{\prime }\rangle$ in this region, it is understood from (4.4) that during wave breaking, the wind speed $\langle u\rangle$ is accelerated and decelerated below and above the peak of $-\langle u^{\prime }w^{\prime }\rangle$ , respectively. The peak of $-\langle u^{\prime }w^{\prime }\rangle$ moves upward with its magnitude decreasing after $t=1.0T$ , leading to an expansion of the acceleration region. As a result, the mean velocity increases after $t=1.0T$ (figure 6 b).

Figure 12(b) shows that at $t=0.0T$ , the value of $-\langle \tilde{u} \tilde{w}\rangle$ is positive near the wave crest height in case A55. It becomes negative above $z/\unicode[STIX]{x1D706}=0.15$ , and diminishes away from the wave surface ( $z/\unicode[STIX]{x1D706}>0.3$ ), indicating that the wave-coherent motion is confined in the near-surface region. At $t=1.0T$ , the value of $-\langle \tilde{u} \tilde{w}\rangle$ above the wave crest becomes positive. After $t=2.0T$ , the magnitude of $-\langle \tilde{u} \tilde{w}\rangle$ becomes negligibly small. Although the profile of $-\langle \tilde{u} \tilde{w}\rangle$ experiences a drastic change during wave breaking in case A55, it should not be regarded as the main reason for the occurrence of the peak in the profile of $-\langle u^{\prime }w^{\prime }\rangle$ , because the magnitude of $-\langle \tilde{u} \tilde{w}\rangle$ is much smaller than that of $-\langle u^{\prime }w^{\prime }\rangle$ . In the presence of wave plunging, the turbulence momentum flux (or, Reynolds stress) plays a much more important role than the wave-coherent momentum flux.

Figure 13.

Vertical profiles of (a) plane-averaged total momentum flux $-\langle u^{\prime }w^{\prime }\rangle /u_{\ast }^{2}$ and (b) wave-coherent momentum flux $-\langle \tilde{u} \tilde{w}\rangle /u_{\ast }^{2}$ in case B55 at $t=0.0T$ , $1.0T$ , $2.0T$ and $3.0T$ . The result of case B35 at $t=0.0T$ is plotted for comparison.

Figure 14.

Vertical profiles of (a) plane-averaged total momentum flux $-\langle u^{\prime }w^{\prime }\rangle /u_{\ast }^{2}$ and (b) wave-coherent momentum flux $-\langle \tilde{u} \tilde{w}\rangle /u_{\ast }^{2}$ in case C55 at $t=0.0T$ , $1.0T$ , $2.0T$ and $3.0T$ . The result of case C35 at $t=0.0T$ is plotted for comparison.

Figures 13 and 14 show the vertical profiles of $-\langle u^{\prime }w^{\prime }\rangle$ and $-\langle \tilde{u} \tilde{w}\rangle$ in cases B55 and C55, respectively. At $t=0.0T$ , the value of $-\langle u^{\prime }w^{\prime }\rangle$ near the wave surface is larger in cases B55 and C55 than in cases B35 and C35, respectively. In case C55, the value of $-\langle u^{\prime }w^{\prime }\rangle /u_{\ast }^{2}$ is larger than 1.0 below $z/\unicode[STIX]{x1D706}=0.30$ (figure 14 a). According to equation (4.5), this indicates that the value of $\unicode[STIX]{x2202}\langle u\rangle /\unicode[STIX]{x2202}z$ is negative, which is consistent with the observation from figure 6(d). At $t=1.0T$ , the value of $-\langle u^{\prime }w^{\prime }\rangle$ becomes negative in both cases B55 and C55. As discussed in § 4.4, this corresponds to the transport of momentum from the vicinity of the wave surface to the region away from the wave surface. At $t=3.0T$ , the value of $-\unicode[STIX]{x2202}\langle u^{\prime }w^{\prime }\rangle /\unicode[STIX]{x2202}z$ is negative below $z/\unicode[STIX]{x1D706}=0.2$ in cases B55 and C55. According to (4.4), this tends to decelerate the wind speed $\langle u\rangle$ near the wave surface, which is consistent with the results shown in figure 6(c,d).

A comparison among figures 12–14 indicates that the wave-coherent momentum flux $-\langle \tilde{u} \tilde{w}\rangle$ plays a more important role in cases B55 and C55 than in case A55. As shown in figure 13(b), the profile of $-\langle u^{\prime }w^{\prime }\rangle$ features a similar shape to that of $-\langle \tilde{u} \tilde{w}\rangle$ at $t=1.0T$ and $2.0T$ in case B55. In case C55, the contribution of $-\langle \tilde{u} \tilde{w}\rangle$ is more dominant. As shown in figure 14, the value of $-\langle \tilde{u} \tilde{w}\rangle /u_{\ast }^{2}$ is larger than 1.0 below $z/\unicode[STIX]{x1D706}=0.3$ . In response, the value of $-\langle u^{\prime }w^{\prime }\rangle /u_{\ast }^{2}$ also becomes larger than 1.0 in this region, which in turn leads to the negative value of $\unicode[STIX]{x2202}\langle u\rangle /\unicode[STIX]{x2202}z$ near the wave surface as shown in figure 6(d). At $t=1.0T$ and $2.0T$ , the magnitude of $-\langle \tilde{u} \tilde{w}\rangle$ is comparable to that of $-\langle u^{\prime }w^{\prime }\rangle$ in case C55, indicating that the evolution of $-\langle u^{\prime }w^{\prime }\rangle$ is dominated by the wave-coherent motion, when the wave age is large.

Figure 15.

Contours of ensemble-averaged Reynolds shear stress $-\overline{u^{\prime \prime }w^{\prime \prime }}/c^{2}$ at (a) $t=0.0T$ , (b) $t=0.8T$ and (c) $t=3.0T$ in case A55. The solid and dashed lines represent isopleths of positive and negative values of $-\overline{u^{\prime \prime }w^{\prime \prime }}/c^{2}$ , respectively. The interval of the isopleths is 0.02. The thick dashed-dotted line represents the top bound of the wave surfaces (see § 4.2 for the definition). The waves propagate in $+x$ direction.

Figure 15 displays the contours of the ensemble-averaged Reynolds shear stress $-\overline{u^{\prime \prime }w^{\prime \prime }}$ near the wave surface in case A55. We recall here that the thick dashed-dotted line is the top bound of the wave surfaces, which is different from instantaneous wave surface. The definition of top bound of wave surface is given in § 4.2, and the top bound of the wave surfaces and instantaneous wave surface are compared in figures 7–9. Figure 15(a) shows that at $t=0.0T$ , the value of $-\overline{u^{\prime \prime }w^{\prime \prime }}$ is negative in a thin layer attached to the windward of the wave surface. Above this layer, $-\overline{u^{\prime \prime }w^{\prime \prime }}$ is positive, and a peak is located to the upper left of the windward of the wave surface. At $t=0.8T$ , the coherent vortex to the right of the wave crest enhances the momentum flux around it, and as a result, a large magnitude of $-\overline{u^{\prime \prime }w^{\prime \prime }}$ is observed around the vortex in figure 15(b). From the comparison between figure 15(a,c), it is known that the magnitude of $-\overline{u^{\prime \prime }w^{\prime \prime }}$ is smaller at $t=3.0T$ than at $t=0.0T$ , indicating that the turbulence is reduced after wave plunging. This is consistent with the results of the plane-averaged total momentum flux demonstrated in figure 12(a).

Figure 16.

Contours of ensemble-averaged Reynolds shear stress $-\overline{u^{\prime \prime }w^{\prime \prime }}/c^{2}$ at (a) $t=0.0T$ , (b) $t=0.8T$ and (c) $t=3.0T$ in case B55. See figure 15 for the detailed caption.

Figure 17.

Contours of ensemble-averaged Reynolds shear stress $-\overline{u^{\prime \prime }w^{\prime \prime }}/c^{2}$ at (a) $t=0.0T$ , (b) $t=0.8T$ and (c) $t=3.0T$ in case C55. See figure 15 for the detailed caption.

Figures 16 and 17 display the contours of the ensemble-averaged Reynolds shear stress $-\overline{u^{\prime \prime }w^{\prime \prime }}$ near the wave surface in cases B55 and C55, respectively. From the comparison among figures 15(a), 16(a) and 17(a), it can be seen that at $t=0.0T$ , the region with negative $-\overline{u^{\prime \prime }w^{\prime \prime }}$ on the windward of the wave expands as $c/u_{\ast }$ increases, and the positive and negative contours become more symmetric about the wave crest. The same trend of the ensemble-averaged Reynolds stress contours in response to the increase of $c/u_{\ast }$ has been reported by Yang & Shen (Reference Yang and Shen2010) for wind over non-breaking waves. Similar to that in case A55 (figure 15 b), the spanwise vortex induces large magnitude of $-\overline{u^{\prime \prime }w^{\prime \prime }}$ around it at $t=0.8T$ in cases B55 and C55 (figures 16 b and 17 b). The comparisons between figure 16(a,c) and between 17(a,c) show that the magnitude of $-\overline{u^{\prime \prime }w^{\prime \prime }}$ is slightly smaller at $t=3.0T$ than at $t=0.0T$ in both cases B55 and C55, indicating the reduction of turbulence after wave breaking. This is consistent with case A55 for young waves (figure 15 a,c).

4.6 Kinetic energy of velocity fluctuations

Figure 18.

Vertical profiles of (a) plane-averaged KE-F $k_{pl}/u_{\ast }^{2}$ and (b) KE-W $k_{w}/u_{\ast }^{2}$ in case A55 at $t=0.0T$ , $1.0T$ , $2.0T$ and $3.0T$ . The result of case A35 at $t=0.0T$ is plotted for comparison.

Figure 19.

Vertical profiles of (a) plane-averaged KE-F $k_{pl}/u_{\ast }^{2}$ and (b) KE-W $k_{w}/u_{\ast }^{2}$ in case B55 at $t=0.0T$ , $1.0T$ , $2.0T$ and $3.0T$ . The result of case B35 at $t=0.0T$ is plotted for comparison.

Figure 20.

Vertical profiles of (a) plane-averaged KE-F $k_{pl}/u_{\ast }^{2}$ and (b) KE-W $k_{w}/u_{\ast }^{2}$ in case C55 at $t=0.0T$ , $1.0T$ , $2.0T$ and $3.0T$ . The result of case C35 at $t=0.0T$ is plotted for comparison.

In this section, we further analyse the effect of wave breaking on the kinetic energy of velocity fluctuations (KE-F). Figures 18–20 compare the vertical profiles of the plane-averaged KE-F, $k_{pl}=\langle u_{i}^{\prime }u_{i}^{\prime }\rangle /2$ , and the kinetic energy of wave-coherent velocity (KE-W), $k_{w}=\langle \tilde{u} _{i}\tilde{u} _{i}\rangle /2$ , in cases with young waves (A35 and A55), intermediate waves (B35 and B55) and old waves (C35 and C55), respectively.

From the comparison between (a,b) in figures 18–20, it is known that the contribution from $k_{w}$ to $k_{pl}$ becomes more important as the wave age increases, consistent with the results of $-\langle u^{\prime }w^{\prime }\rangle$ and $-\langle \tilde{u} \tilde{w}\rangle$ displayed in figures 12–14. The wave steepness influences both $k_{pl}$ and $k_{w}$ near the wave surface. As shown in figures 18–20, at $t=0.0T$ , the magnitudes of $k_{pl}$ and $k_{w}$ in cases A55, B55 and C55 are respectively larger than those in cases A35, B35 and C35 below $z/\unicode[STIX]{x1D706}=0.4$ . However, the wave steepness does not influence the KE-F far away from the wave surface. This is consistent with the results of the mean velocity and total momentum flux demonstrated in figures 6, 12–14.

The wave plunging influences the KE-F near the wave surface significantly. As shown in figure 18(a), in case A55, the value of $k_{pl}$ below $z/\unicode[STIX]{x1D706}=0.55$ is higher at $t=1.0T$ than at $t=0.0T$ , indicating that the turbulence is transiently enhanced due to the wave breaking. The KE-F also experiences a transient growth in case C55. However, we do not observe such a transient growth of KE-F in case B55. The turbulence is eventually reduced after wave plunging ceases in cases A55, B55 and C55.

Figure 21.

Vertical profile of plane-averaged KE-F production $P_{k}/(u_{\ast }^{3}/\unicode[STIX]{x1D706})$ in cases (a) A55, (b) B55 and (c) C55.

To further study why the transient growth of KE-F takes place in cases A55 and C55 but is absent in case B55, we compare the production of KE-F, $P_{k}$ , in these three cases in figure 21. Here, $P_{k}$ is defined as

(4.6) $$\begin{eqnarray}\displaystyle P_{k}=-\langle u^{\prime }w^{\prime }\rangle \frac{\unicode[STIX]{x2202}\langle u\rangle }{\unicode[STIX]{x2202}z}. & & \displaystyle\end{eqnarray}$$

The production term represents a rapid effect on the KE-F (Townsend Reference Townsend1980), and the behaviour of KE-F is often dominated by its production term (Mansour, Kim & Moin Reference Mansour, Kim and Moin1988; Smits, McKeon & Marusic Reference Smits, McKeon and Marusic2011). Note that the full production term of KE-F reads

(4.7) $$\begin{eqnarray}\displaystyle P_{k}=-\langle u_{i}^{\prime }u_{j}^{\prime }\rangle \frac{\unicode[STIX]{x2202}\langle u_{i}\rangle }{\unicode[STIX]{x2202}x_{j}}. & & \displaystyle\end{eqnarray}$$

However, according to the definition of plane averaging, $\unicode[STIX]{x2202}\langle \;\rangle /\unicode[STIX]{x2202}x=\unicode[STIX]{x2202}\langle \;\rangle /\unicode[STIX]{x2202}y=0$ holds theoretically. Furthermore, the magnitudes of the mean spanwise velocity $\langle v\rangle$ and mean vertical velocity $\langle w\rangle$ are negligibly small in comparison with the mean streamwise velocity $\langle u\rangle$ . As a result, the production of KE-F is simplified to (4.6).

It is evident from figure 21(a,c) that in cases A55 and C55, the magnitude of $P_{k}$ is larger at $t=1.0T$ than at $t=0.0T$ . However, the physical process leading to the increase of $P_{k}$ is different in these two cases. From equation (4.6), it is known that the value of $P_{k}$ is determined by the vertical gradient of mean velocity $\unicode[STIX]{x2202}\langle u\rangle /\unicode[STIX]{x2202}z$ and the total momentum flux $-\langle u^{\prime }w^{\prime }\rangle$ . In case A55, $\unicode[STIX]{x2202}\langle u\rangle /\unicode[STIX]{x2202}z$ decreases during the first plunging event (as can be seen from the $\langle u\rangle$ profile in figure 6 b), which tends to reduce the value of $P_{k}$ . On the opposite, the magnitude of $-\langle u^{\prime }w^{\prime }\rangle$ increases and a positive peak occurs in its profile (figure 12 a) due to the disturbance effect of the wave breaking on the turbulent flow, which tends to cause an increase of $P_{k}$ . Figure 21(a) indicates that the latter effect is stronger in the young wave case, and as a result $k_{pl}$ experiences a transient growth during the first plunging event (figure 18 a). In case C55, the transient growth of KE-F is correlated to the negative value of $\unicode[STIX]{x2202}\langle u\rangle /\unicode[STIX]{x2202}z$ near the wave surface (figure 6 d), which leads to a negative value of $P_{k}$ at $t=0.0T$ as shown in figure 21(c). During the first plunging event, because the wave breaking induces a large negative value of $-\langle \tilde{u} \tilde{w}\rangle$ (figure 14 b), the value of $-\langle u^{\prime }w^{\prime }\rangle$ becomes negative correspondingly (figure 14 a). As a result, $P_{k}$ becomes positive at $t=1.0T$ near the wave surface. Consequently, the magnitude of KE-F also increases transiently in case C55 during the first plunging event. In case B55, however, none of the two processes described above is present. As demonstrated in figure 13(a), no positive peak occurs in the profile of $-\langle u^{\prime }w^{\prime }\rangle$ during the first plunging event, while figure 6(c) shows that $\unicode[STIX]{x2202}\langle u\rangle /\unicode[STIX]{x2202}z$ is positive in case B55. Therefore, the turbulence does not experience a transient growth in this intermediate wave case. The production of KE-F is smaller at $t=3.0T$ than at $t=1.0T$ in all cases, which leads to the magnitude decrease of KE-F after wave plunging (figures 18 a, 19 a and 20 a).

Figure 22.

Contours of ensemble-averaged TKE, $k_{ph}/c^{2}$ , at (a) $t=0.0T$ , (b) $t=0.8T$ and (c) $t=3.0T$ in case A55. The interval of the isopleths is 0.1. The thick dashed-dotted line represents the top bound of wave surfaces (see § 4.2 for the definition). The waves propagate in the $+x$ -direction.

Figure 22 displays the contours of the ensemble-averaged TKE, $k_{ph}=\overline{u_{i}^{\prime \prime }u_{i}^{\prime \prime }}/2$ , near the wave surface in case A55. As shown in figure 22(a), at $t=0.0T$ , the value of $k_{ph}$ is small near the leeward of the surface. Figure 22(b) shows that at $t=0.8T$ , $k_{ph}$ increases near the large spanwise vortex (cf. figure 7 c3), indicating that the turbulence is enhanced by the vortex. This is consistent with the result of $-\overline{u^{\prime \prime }w^{\prime \prime }}$ shown in figure 15(b). At $t=3.0T$ , as the turbulence is reduced, the magnitude of $k_{ph}$ is smaller than that at $t=0.0T$ .

Figure 23.

Contours of ensemble-averaged TKE, $k_{ph}/c^{2}$ , at (a) $t=0.0T$ , (b) $t=0.8T$ and (c) $t=3.0T$ in case C55. The interval of isopleths is 0.02 in (a,c), while it is 0.1 in (b). The thick dashed-dotted line represents the top bound of wave surfaces (see § 4.2 for the definition). The waves propagate in the $+x$ -direction.

Figure 23 illustrates the contours of $k_{ph}$ near the wave surface in case C55. As shown in figure 23(a), the peak of $k_{ph}$ occurs near the wave surface at $t=0.0T$ . Similar to that in case A55, the magnitude of $k_{ph}$ increases significantly at $t=0.8T$ near the spanwise vortex (cf. figure 9 b3). The value of $k_{ph}$ near the vortex at $t=0.8T$ scales with $c^{2}$ in cases A55 and C55, indicating that the high TKE in this region is dominated by the wave breaking. At $t=3.0T$ , the magnitude of $k_{ph}$ near the wave surface becomes smaller than that at $t=0.0T$ . This is similar to the results of case A55 shown in figure 22.

The figures of $k_{ph}$ of case B55 are not shown here to keep the paper concise. Similar to that in cases A55 and C55, the magnitude of TKE in case B55 experiences an increase during the wave plunging near the wave surface, and is eventually reduced after the wave breaking. However, the TKE far away from the wave surface does not increase transiently in case B55 during wave plunging (figure 19 a). This indicates that although TKE is generated in the vicinity of the wave surface during wave plunging, little is transported to high altitude before it is dissipated. As shown above, the transient growth of KE-F away from the wave surface during wave plunging in cases A55 and C55 is caused by the large magnitude of local KE-F production.