Study on the Nearshore Evolution of Regular Waves under Steady Wind

Abstract

We present a study on regular wave propagation on a sloping bed under the action of steady wind, which is of a great significance to complement and replenish the interaction mechanisms of nearshore wave and wind. Physical experiments were conducted in a wind-wave flume, and the corresponding numerical model was constructed based on the solver Waves2FOAM in OpenFOAM, with large-eddy simulation (LES) used to investigate the turbulent flow. The comparisons between the measured and calculated results of the free surface elevation and flow velocity indicated that the numerical model could predict the associated hydrodynamic characteristics of a nearshore wave regardless of the presence or absence of wind. The results showed that wind had a significant impact on nearshore wave evolution. It was found that under the same wind speed coverage constraint, wave breaking occurred ahead of time. The smaller the surf similarity ${\xi }_0$ was, the higher the dispersion degree of wave breaking locations would be, and the breaker index of $H_b/h_b$ increased with wind speed under the same incident wave height. The main components of analysis for turbulent flow were the results of the cross-spectrum, the TKE (turbulent kinetic energy), and TDR (turbulent dissipation rate). The cross-spectrum illustrated that wind enhanced the degree of coherence of the residual velocity components and aggravated turbulence. The TKE indicated that in regions near the water surface, wind speed made it considerably larger and the average level rapidly decreased with depth. The TDR exhibited that the significant effect of wind was merely imposed after breaking, wherein the turbulence penetrated the deeper flow and the average level generally rose. The velocity profile on the slope showed that the wind accelerated the undertow, and the moment statistics indicated that the velocity distribution deviated gradually from the Gaussian distribution to the right.

1. Introduction

Wind is one of the driving forces of surface waves. When wind blows over the water surface, the water surface will fluctuate due to wind pressure, and the fluctuation will further develop into wind wave in the processes of continuous momentum and energy transfer in air-sea interaction. Under extreme climate, e.g., extreme wind waves and storm surges, it endangers the safety of coastal residents and buildings over the long term. In a marine dynamic environment, the mechanism of wind wave formation and propagation is a research topic that concerns many researchers. Abundant research results have been obtained ranging from field observation [1,2,3] to theoretical analysis [4,5,6,7], as well as experimental [8,9,10] and numerical study [11,12,13,14,15,16]. These studies are mainly focused on the generation of wind wave and air-sea interaction, rather than the mechanism of wind effect on waves. However, wind waves are usually inseparable from wind field when they propagate from deep water to shallow water, and the law of wave propagation with the action of wind is more complex than that of pure wind wave or swell in coastal regions.

Researchers have realized the action of wind on nearshore wave propagation since the 1950s, and managed to obtain some results [17]. Until now, research on this topic has been mainly implemented by experimental and numerical methods that investigate the relationship among different wind speed, wind direction, together with the shape, height, and location of a breaking wave. For field observation, it is difficult to obtain the temporal and spatial evolution law of waves due to the large scale of strong winds such as typhoons, and thus only some qualitative investigations were conducted where the existence of wind can change the type of breaking coastal waves [18,19,20]. Galloway et al. (1989) [21] conducted field measurements and found that offshore wind and onshore wind increased the number of plunging and spilling, respectively. Some researchers have carried out laboratory experiments and verified observations and findings of Galloway et al. For example, Douglass (1990) [22], King and Baker (1996) [23] had experimentally studied the influence of local wind on a nearshore breaking wave, and found that wind had a significant impact on the location, geometry, and type of breaking wave. Onshore wind was found to cause waves to break ahead of time in offshore deep water, whereas the offshore wind was found to delay the wave breaking in nearshore shallow water. Meanwhile, Feddersen and Veron (2005) [24] stated that wind could enhance the wave energy in shallow water and cause changes in wave shape. With respect to the numerical simulation, Kharif (2008) [25], Tian and Choi (2013) [26], Liu (2016) [27], and Hasan (2018) [28] applied the Jeffery’s shelter theory and Miles’ shear instability theory [4,5] to simulate the impact of wind on waves. While the core of the two theories is to predefine the wind pressure terms through empirical parameterization, it is essentially the single-phase flow and cannot truly embody the interaction processes between wind and wave. Recently, the enhanced computer performance had enabled the numerical simulation of multiphase flow, in which the two-phases flow model had been widely employed in the coupling of wind and wave [29,30,31]. The two-phases flow model can simultaneously perform the coupling between two different fluids, which avoids empirical parameterization, making it more reasonable compared with the single-flow model by principle and applicability [32]. The existing two-phase models can be subdivided into two kinds, one of which is solved in the air and water respectively, and the predicted air and water flow are coupled by the velocity continuity and stress balance at the free surface [12,33]. Despite that, this kind of model is not suitable for the simulation of a free surface related to deformation, breaking, and air entrapment, etc. [34,35]. The other depends on the Navier–Stokes equation where the two phases are solved simultaneously with an assumption that treats the two phases as a kind of fluid in the fixed Eulerian grids [36]. As for the complex topological deformation of free surface, e.g., the geometric VOF (volume of fluid) method [37,38] and level set method [39,40], are usually used to fulfill free surface reconstruction and transportation.

Moreover, the turbulence structure and undertow in the process of wave nearshore propagation are worth researching and some researchers have already conducted studies on the topic. Hansen (1984) [41] deduced the theoretical results for the undertow, and Ting (1994, 1995) [42,43] experimentally investigated the turbulence and undertow in the surf zone for spilling and plunging waves. It was found that the average turbulence levels were much higher and the vertical variations of undertow were much smaller for plunging waves than those for spilling waves, which were numerically confirmed by Bakhtyar (2009) [44]. In addition, research on the turbulent boundary layer at the bottom of a nearshore wave has received attention. Cavallaro (2011) [45], Kranenburg (2012) [46], and Blondeaux (2018) [47] compared the precision and suitable conditions with respect to the bottom shear stress and TKE (turbulent kinetic energy), among different turbulent models. However, these studies were all conducted in such a way that the wind was not taken into account. Consequently, given the lack of experimental and numerical research, no clear conclusions have been made with respect to the wind effects on wave propagation. Likewise, the interaction between wave and wind is associated with an extremely complex physical mechanism that has not been completely solved. Therefore, it is of great significance to complement and replenish relevant studies, which is the main purpose of this paper.

In this paper, the 2nd Stokes wave propagation on the slope under steady wind was investigated. The solver Waves2FOAM in OpenFOAM, by equipping it with the wave-making module and active sponge layer wave attenuation technology, was applied to simulate the wind and wave field. Additionally, the VOF method with an artificial compression term was used to capture air-water interface, with a large-eddy simulation (LES) for the turbulent flow. The remaining part of the paper is organized as follows. In Section 2, the governing equation, turbulence model, and boundary conditions are briefly introduced. In Section 3, the experimental and numerical setup are described. The accuracy of the numerical model is verified by experimental data in Section 4, and the hydrodynamic characteristics of wave under wind action are discussed in Section 5. A summary of the main findings are given in Section 6.

2. Numerical Methods

2.1. Governing Equations

On the basis of the mass and momentum conservation of fluids, the continuity equation and momentum equation (Navier–Stokes equation) for incompressible viscous fluid are expressed as follows:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathit{U}\right)=0$$

(1)

$$\frac{\partial \rho \mathit{U}}{\partial t}+\nabla \cdot (\rho \mathit{U}{\mathit{U}}^T)=-\nabla p+\nabla \cdot \left(\mu \nabla \mathit{U}+\rho \tau \right)-\mathit{g}\cdot \mathsf{{\rm X}}\nabla \rho +{\mathit{f}}_{\sigma },$$

(2)

where $t$ is the time, $\rho$ is the fluid density, $\mathit{U}=\left(u, w\right)$ is the fluid velocity vector, $\mathbf{X}=\left(x, z\right)$ is the Cartesian coordinates, $\tau =\frac{1}{\rho }{\mu }_t\left(\nabla \mathit{U}+{\left(\nabla \mathit{U}\right)}^T\right)-\frac{2}{3}k\mathit{I}$ is the stress tensor in which ${\mu }_t$ denotes the dynamic viscosity, $k$ is the turbulent kinetic energy and $\mathit{I}$ is the Kronecker delta function, $\mathit{g}$ is gravitational acceleration, and ${\mathit{f}}_{\sigma }$ is the surface tension.

The VOF model has been adopted to capture the interface of the two-phase flow with artificial compression term [37,48]. The transport equation is expressed as:

$$\frac{\partial {\alpha }_p}{\partial t}+\nabla \cdot \left({\alpha }_p\mathit{U}\right)+\nabla \cdot \left({\alpha }_p\left(1-{\alpha }_p\right){\mathit{U}}_c\right)=0$$

(3)

where ${\alpha }_p$ represents the volume fraction of different phase fluids as follows:

$${\alpha }_p=\{\begin{array}{l}0, air\\ 1, water\\ 0<{\alpha }_p<1, interface\end{array}$$

(4)

${\mathit{U}}_c$ is the velocity field only applicable to the compression interface, which can be expressed by the following formula:

$${\mathit{U}}_c=\min \left[C_{{\alpha }_p}\left|\mathit{U}\right|,\max \left(\left|\mathit{U}\right|\right)\right]\frac{{\nabla }_{{\alpha }_p}}{\left|{\nabla }_{{\alpha }_p}\right|}$$

(5)

where the compression rate $C_{{\alpha }_p}$ is expressed by a default of 1, and the curvature of the interface is expressed by $\frac{{\nabla }_{{\alpha }_p}}{\left|{\nabla }_{{\alpha }_p}\right|}$.

The VOF model treats the fluid of different phases as mixed phases, wherein the fluid characteristics, involving the fluid density and dynamic viscosity, are mixed characteristics that can be written in the form of the following formula:

$$\{\begin{array}{l}\rho ={\alpha }_p{\rho }_w+(1-{\alpha }_p){\rho }_a\\ \mu ={\alpha }_p{\mu }_w+(1-{\alpha }_p){\mu }_a\end{array}.$$

(6)

2.2. Turbulence Model

The turbulence flow is accompanied with the processes of nearshore wave propagation, especially in wave breaking. In this paper, LES was adopted to model the turbulent flow in the processes of wind-wave interaction, which can be employed to calculate the turbulent structures in various scales and possesses a certain superiority compared with other models such as direct numerical simulation (DNS) and the Reynolds Average Navier–Stokes simulation (RANS) in terms of time cost and computation load. For details, one can refer to Bogey (2006) [49], and thus they are not repeated here.

2.3. Boundary Conditions

The realization flow of the inlet for wave making and autonomous wave absorption was depicted in considerable detail in Higurea (2013) [50], which will not be further covered here. The velocity function and free surface function derived from the 2nd Stokes wave are as follows:

$$u=A\sigma \left[\frac{coshk\left(z+h\right)}{coshkh}cos\theta +\frac{3}{4}Ak\frac{cosh2k(z+h)}{sin{h}^{4}kh}cos2\theta \right]$$

(7)

$$w=A\sigma \left[\frac{sinhk(z+h)}{coshkh}sin\theta +\frac{3}{4}Ak\frac{sinh2k(z+h)}{sin{h}^{4}kh}sin2\theta \right]$$

(8)

$$\eta =Acos\theta +\frac{A^{2}k}{4}\frac{coshkh(cosh2kh+2)}{sin{h}^{3}kh}cos2\theta$$

(9)

where $A$ denotes wave amplitude. The angular frequency $\sigma =2\pi /T$ and $T$ is the wave period. For the wave number $k=2\pi /L$ and $L$ is the wave length. $\theta =kx-ct$ denotes phase and the wave speed is $c=\sqrt{g(H+h)}$. $h$ is the water depth and $\eta$ is the free surface elevation.

In addition, the wave absorption function is as follows:

$${\alpha }_R\left({\chi }_R\right)=1-\frac{e^{{\chi }_R^{3.5}}-1}{e-1}$$

(10)

$$\varphi ={\varphi }_{calculated}+\left(1-{\alpha }_R\right){\varphi }_{t\arg et}$$

(11)

where ${\chi }_R\in \left[0, 1\right]$ and $\varphi$ denotes $\mathit{U}$ or $\alpha$.

3. Experimental and Numerical Setup

3.1. Experimental Setup

The physical experiment was carried out in a wind wave flume located in the Hydraulics Modeling Laboratory of the Changsha University of Science and Technology. The flume is 45 m long, 0.8 m wide, and 1 m deep. One end of the flume was equipped with a flap-typed wave maker (Dalian University of Technology, China), capable of generating various kinds of waves covering a larger range of wave heights and periods, wherein the 2nd Stokes wave was used in this study. A GXF diagonal fan (Bejing Jifeng Shunda Modern Ventilation Equipment Co., Ltd., China) was installed at the other end of the flume, which could create the wind speed with a range of 0~10 m/s. The schematic layout of the physical model is illustrated in Figure 1a. The flume (see Figure 1b) was covered by the roof of galvanized iron sheet. To prevent the wind from escaping from the gaps between the flume walls and cover plate, sealing plasticine was used to fill the gaps. Referring to the natural beach slope angle and the setting of physical models from Sibul and Ticker (1956) [17], Battjes (1974) [51], and Inami (2016) [52], a slope beach of 1:10 (see Figure 1c) was built with the PVC (polyvinyl chloride) panel at 13.25 m from the wave maker.

Due to the fact that the coastal water depth and wind conditions are too complex to observe and measure in the field environment, the model scale was not strictly considered in this study in view of any prototype. To examine the rationality of the magnitude between the wind speed and wave elements, a scale factor of 1:20 was roughly chosen to calculate the relevant variables of the prototype and experiment. According to the Froud similarity law, the range of incident wave heights was 3.0~11.0 cm, which was comparable to the prototype of 1.0~2.2 m. The maximum experimental wind speed and constant experimental water depth were set to 5 m/s and 0.4 m, corresponding to the prototypes of 22.36 m/s and 8 m, respectively.

The wave elevations were measured at eleven cross-shore locations (S1-S8). To monitor the incident wave in the flat-bed zone, the resistance-typed water level sensor (AWP-24-3, Pivotal Scientific Ltd., Canada) was placed at S1, with an accuracy of 0.1 mm and a sampling frequency of 50 Hz was used. S1 was located at a position 10.8 m away from the wave maker, where the wave reflection was relatively minor. Therefore, it could be treated as the reference location of stable incident wave height. Meanwhile, to measure the nearshore wave shoaling and breaking, the non-intrusive water level sensors (ULS 80D, General Acoustics Ltd., Germany), with an accuracy of 0.1 mm and setting sampling frequency of 50 Hz, were equally spaced from S2 to S8 with a distance of 0.6 m on the slope, wherein S2 was located at the slope toe. A velocimeter (Vectrino, Notek Ltd., Norway) was arranged at the slope toe to measure the flow velocity. Moreover, to monitor the wind speed, wind anemometers with a pitot tubes type (DP1000-T1, Shanghai Lingsheng Electronic Instrument Co., Ltd., China) and an accuracy of 1 mm/s and setting sampling frequency of 8 Hz, were installed at W1 (7.6 m) and W2 (13.1 m) with a height of 10 cm above the still water surface. Notably, Wu (1969) [53] proposed 10 cm as the height of characteristic wind speed in the experiment scale. In the absence and presence of wind, 120 and 240 waves were separately generated within 3 min and 6 min, respectively. The data in the last 2 min were considered valid and used for verification.

3.2. Numerical Setup

In accordance with the experimental setup in Figure 1, the main aspects of the experimental model were reproduced by setting the relevant numerical boundary conditions. According to Figure 2, the wave making and wind speed conditions were set at the inlet. Meanwhile, in an attempt to weaken the effect of wave reflection, the autonomous wave absorption boundary, with the length of the wave absorption zone set to 3 m to satisfy the parameter configuration at the request of at least one wavelength, was adopted at the inlet. Boundary conditions for other faces were: A slip wall condition at the top boundary, the no-slip wall conditions at the bottom and slope boundary, as well as the pressure outflow condition at the outlet. As for turbulent eddy viscosity, all boundary conditions were zero gradient. Additionally, a two-dimensional cartesian coordinate system was defined as indicated in Figure 2, where the coordinate origin was set at the slope foot and the positive direction of axis $x$ and $z$ represent the horizontal direction of wave propagation and the vertical upward of the flume, respectively.

The computational domain was composed of two sections, which was discretized with structured grids made up of a total of 0.5 million cells. In the flat-bed zone, the coarser grids, with a constant size of 0.0065 m × 0.0065 m, were used to calculate the steady evolution of wind and wave fields, which could effectively reduce the computation load on the precondition of assuring the calculation precision. Meanwhile, the asymptotic grids with an average size of 0.005 m, which varied in the $x$ and $z$ direction, were used in the slope zone to capture the wave breaking. The grid size reduced gradually from 0.0065 m to 0.0035 m in the $x$ direction, and owing to the contraction of outlet it fell to 0.001 m and 0.0035 m at the lower and upper edge of the outlet, respectively. The physical features of air and water were quite different. To ensure that the global Courant number was suitable for the calculation convergence, the time step had to be small enough. Therefore the adaptive time-step was employed during the simulation, with a maximum value of 0.025. In the meantime, the maximum Courant number was set at 0.2.

4. Model Validation

Based on the experimental and numerical model discussed above, the validity of the numerical model for the hydrodynamic characteristics such as free surface elevations and flow velocity was verified with experimental data in the typical cases of W0H07 and W5H07.

The comparison shows good agreement with the average root mean square error (RMSE) below 0.6 between the measured and computed time series of free surface elevations at the seven measured locations (S2-S8) with and without wind, as given in Figure 3 and Figure 4, respectively. At the offshore locations (S2-S4), the wave surface was symmetrical with a relatively stable wave height. When the wave propagated from S5 to S7, it became rapidly steeper and close to breaking due to the shoaling effect. After that the wave breaking occurred at a location between S7 and S8, with a sharp decrease of the maximum free surface elevation at S8.

Moreover, the measured and computed results of the two-dimensional velocity at a depth of 0.2 m (anemometer depth) at the slope toe cross-section in the windless case of W0H07 are presented in Figure 5, in which the measured results were nearly horizontally and vertically identical to that of computed, and the values were basically equally distributed near zero. In contrast to the windy case of W5H07, as indicated in Figure 6, despite the fact that the measured amplitude of horizontal velocity was slightly smaller than the computed results, the average levels show the same overall downtrend, which reflects the action of wind to some extent. The vertical velocity distribution was less affected by wind with no significant differences.

Overall, the comparisons between the measured and calculated results of the free surface elevation and flow velocity indicates a good model prediction regardless of the presence or absence of wind. Hence, 12 scenarios were simulated for subsequent analysis based on a combination of wind speeds and wave heights in

Table 1. Simulated scenarios (2nd Stokes).

Case	W (m/s)	H0 (cm)	T (s)	h (m)	ξ0
W0H05~W0H11	0.0	5, 7, 9,11	1.5	0.4	plunging
W3H05~W3H11	3.0	5, 7, 9,11	1.5	0.4	plunging
W5H05~W5H11	5.0	5, 7, 9,11	1.5	0.4	plunging

, wherein the surf similarity parameter ${\xi }_0=tan\alpha /{\left(H_0/L_0\right)}^{1/2}$ ranging from 0.57~0.84 satisfies the plunging wave.

5. Results and Discussions

In this section, we will quantitatively analyze the influence of wind on the evolution of near-shore waves taking into account wave breaking, turbulent flow, and velocity profile, which are crucial to the studies in air-sea interaction, hydrodynamics in surf-swash zone, sediment transport in coastal areas, and wave-structure interaction.

5.1. Wave Breaking

As described in Section 4, the processes of wave propagation on the slope was accompanied with the wave distortion and wave breaking, and little is known about wave breaking under wind action up to now.

The similarity shared by different cases makes it a feasible practice to focus on typical cases here. As a result, we merely present the instantaneous snapshots of wave surface during wave breaking processes in the typical cases of W0H07~W5H07, as indicated in Figure 7, in which three moments are defined separately: Initial breaking, breaking, and overturning jet initiation. The figures indicate that the breaker type was almost unchanged among the three cases, yet there exists some differences, as reflected in the wave breaking ahead of time with the increase of wind speed. On the other hand, taking the initial breaking time in the absence of wind as a reference, we supplemented the compared results of all cases as evident in Figure 8. As an example, Figure 8a indicates that although waves with a wind speed of 3 m/s and 5 m/s were already broken at $t/\sqrt{h/g}=94$, it is just the beginning of initial breaking in the windless case. Such features are apparently also included in other cases (Figure 8b–d) and not repeated here. It suggests that the wave breaking location would migrate towards the offshore direction with the impact of wind, which probably caused changes of breaking height and depth, and even critical breaking conditions. This clearly addresses a fact that the role of wind on wave breaking in time is coincident with that in space throughout the breaking processes. As indicated in Figure 8a, the breaking courses approached each other and even overlapped for different wind speeds. It is especially true in the case of a smaller incident wave height. As the incident wave height increased, this behavior tended to disappear and the spatial distribution of each course seemed to disperse, as shown in Figure 8b–d. It can be inferred that the incident wave height was one of the main factors affecting the dispersion degree of wave breaking locations under the same wind speed coverage constraint.

Taking into account the critical breaking depth, the wave breaking location can be calculated by the following formula:

$$x_b=L_x-\frac{h_b}{tan\alpha }$$

(12)

where the wave breaking location and breaking depth are represented by $x_b$ and $h_b$ and $L_x$ is the horizontal distance from slope toe to the shoreline.

According to the proposal by Collins and Weir (1969) [54], the breaking wave height is proportional to the breaking depth, which can be expressed as:

$$r_b=H_b/h_b=0.72+5.6tan\alpha$$

(13)

as well as by Le Mehuate (1967) [55] where the relationship between the breaking wave height and the height in deep water can be expressed as:

$$H_b/H_0=0.76{\left(tan\alpha \right)}^{1/7}{\left(H_0/L_0\right)}^{-1/4}.$$

(14)

Then, on the combination of Equations (12)–(14), the wave breaking locations can be predicted with the below formula:

$$x_b=L_x-\frac{0.76{\left(tan\alpha \right)}^{1/7}{(H_0/L_0)}^{-1/4}{H}_0}{0.72+5.6tan\alpha }\frac{1}{tan\alpha }.$$

(15)

The correlation between the numerical dimensionless wave breaking locations $\left(L_x-x_b\right)/h$ and surf similarity parameter (${\xi }_0$) is given in Figure 9. In view of ${\xi }_0$ embodying the intrinsic qualities of the slope, the results of the added cases of $H_0=5 \mathrm{cm}$ with the slope of 1:16 and 1:20 where ${\xi }_0$ satisfies the spilling wave are included. Also shown are the distribution curves of breaking locations by Equation (15). Under the condition of the same wind speed and breaking parameter, the wave breaking location of the 2nd Stokes wave was farther from the shoreline than the calculated results, mainly caused by the discrepancy between numerical calculation and empirical parameterization. Likewise, the distribution of wave breaking locations was identical in essence with the above observation in Figure 8, which lends support to a more accurate deduction. It could be deduced that under the same wind speed coverage constraint, the smaller the ${\xi }_0$ is, the higher dispersion degree of wave breaking locations will be, wherein the ${\xi }_0$ involves the slope factor, rather than simply the incident wave height.

To find out the influences on wave properties along the cross-shore distance induced by the wind, we plotted the cross-shore variation of the average wave height and water level in the typical cases, as shown in Figure 10. It is evident that wind raised the average wave height and average water level at the same cross-shore location before wave breaking. For the critical breaking height and depth, the effect would be enhanced wherein the variation of the former was significantly greater than the latter. As a result, it is expected that the critical breaker indices like $H_b/h_b$ will present an uptrend. We accessed the wind effect on the breaker index of $H_b/h_b$, as shown in Figure 11. With the absence of wind, $H_b/h_b$ was basically maintained at a constant level with a value of 1 except at ${\xi }_0=0.42$ where it was a gentle slope of 1:20 and the results agreed well with that given by Galvin (1968) [56]. In the presence of wind, the overall trend of $H_b/h_b$ showed an ineligible increase with wind speed, implying that the critical breaker index without wind was no longer suitable for the windy cases any more.

5.2. Turbulent Flow

Nearshore wave breaking is the main mode of wave energy consumption, in which the turbulence is fierce. It is generally accepted that the existence of wind will further heighten the turbulence level [57]. Hence, it is essential to quantitatively investigate the impact of wind on the turbulent flow when wave propagates on the slope.

Firstly, the cross-spectrum method (for more details see Appendix A) was applied to examine the correlation property of horizontal and vertical residual velocity, which can effectively reflect the turbulent features. On the support of the low-pass filtering in LES, the filtered residual velocity is:

$$\{\begin{array}{l}{u}^{\prime }=u-\overline{u}\\ w^{\prime }=w-\overline{w}\end{array}$$

(16)

where $u$ and $w$ are filtered velocity, the overbar and comma indicate the time-averaged and residual quantities. As illustrated in Figure 12, in cases of W0H07~W5H07, the results of the cross-spectrum analysis of three levels below the trough at the location of $\left(x-x_b\right)/h=-1.25$ show that the spectral energy was mainly concentrated in the low-frequency region ($f<2.5 \mathrm{Hz}$) where the wind input was absorbed and transported. In this region the cross spectral density, coherence, and phase lag of $u^{\prime }$ and $w^{\prime }$ were basically the same for different levels and the wind speed increased the coherence peak at a depth of 0.37 cm, rising from 0.76 to 0.9, which means that the wind enhanced the degree of coherence of velocity components, and aggravated turbulence. Then, the absorbed wind input spreads further, to a high-frequency region ($f>2.5 \mathrm{Hz}$) by nonlinear interaction, which mainly reflected that the spectral density with wind was larger than that without wind. Remarkably, the phase lag tended to disperse at $f=2.5 \mathrm{Hz}$ for different levels, with the degree increasingly evident as the wind speed increased.

Subsequently, the generation and dissipation of TKE are worth analyzing within the processes of wave breaking. Combining the kinetic energy governing equation of LES, the following formula was used to calculate TKE in the solvable scale [58]:

$$k=\frac{1}{2}\left(u^{\prime 2}+w^{\prime 2}\right).$$

(17)

Part of TKE in a solvable scale will dissipate in small scale fluctuations, which is defined as the TDR (turbulent dissipation rate) that cannot be solved directly from the LES. Via the kinetic energy equation for solvable scale [58], the TDR is found to be similar to the following form:

$$\epsilon =2<{\tau }_{ij}\overline{S_{ij}}>,$$

(18)

where ${\tau }_{ij}$ is the stress term and $\overline{S_{ij}}=\frac{1}{2}\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)$ denotes strain-rate tensor. A modeling method is proposed by Smagorinsky (1963) [59] to calculate ${\tau }_{ij}$, namely the well-known sub-grid stress model:

$${\tau }_{ij}={\upsilon }_{SGS}\overline{S_{ij}}={(C_{SGS}\Delta )}^2\left|\overline{S_{ij}}\right|\overline{S_{ij}}$$

(19)

where $\left|\overline{S_{ij}}\right|=\sqrt{2\overline{S_{ij}}\overline{S_{ij}}}$ is the module of analytic deformation rate, ${\upsilon }_{SGS}$ is the Smagorinsky dynamic viscosity, $C_{SGS}$ is the Smagorinsky constant of 0.16, and $\Delta$ represents the Smagorinsky scale that is identical to the magnitude of grids in the computational domain.

In accord with the time in Figure 7, Figure 13 and Figure 14 show the relevant snapshots of distribution of TKE and TDR obtained by Equations (17) and (18) during wave breaking. As shown in Figure 13, the maximum TKE occurred under the wave crest, near the windward wave surface. It could also be found that the turbulence owing to the surface roller was rapidly transported to the depths of the flow and the rear regions of wave front by the turbulence transformation mechanisms, and the TKE was more distributed horizontally rather than vertically [42,44]. Meanwhile, the distribution of TDR was relatively uniform with a constant order of magnitude in most regions except for the wave front and roller where the mixing of air and water with notable values exists, as shown in Figure 14. The influence of wind on TKE and TDR can be observed as well in Figure 13 and Figure 14, and as a result, there exists some differences that will be discussed in detail below.

We plotted the time variations of TKE and TDR at a depth of 0.37 m below the trough level at two cross-sections before and after breaking based on the above cases, as shown in Figure 15 and Figure 16, respectively. The results show that the TKE and TDR displayed a periodic variation throughout the collection time, wherein the TKE was a bimodal distribution in each wave period. Before wave breaking ($\left(x-x_b\right)/h=-0.35$), as plotted in Figure 15a, the maximum TKE occurred at an initial breaking time and increased with wind speed. In sharp contrast with the TKE, the TDR barely felt the impact of wind, as displayed in Figure 15b. The maximum TDR occurred before initial breaking time, resulting in the phase advance. Therefore, the coherence was relatively weak as compared to the case in the TKE. Likewise, Figure 16 shows similar results after wave breaking ($\left(x-x_b\right)/h=0.375$), since the mixing of air and water in the inner part of surf zone enhanced the turbulence intensity. Particularly in the presence of wind, the average TKE and TDR levels generally rose and were significantly larger than that before breaking, which was likely to cause the TKE and TDR to diffuse downward with more complex turbulence transformation mechanisms [44,60]. Then, we plotted the variation of time-averaged TKE and TDR with depth in Figure 17. As expected, in Figure 17a,c, both before and after breaking, the time-averaged TKE in regions near the water surface was considerably larger and rapidly decreased with the depth. It is also noted that wind enhanced the average TKE level. For time-averaged TDR, Figure 17b presents that before breaking it was found to be negligible even in the windy cases, while the significant effect of wind was imposed after breaking as given in Figure 16b, the turbulence penetrated to the deeper flow and the average level of TDR increased with wind speed.

5.3. Undertow

For the regular wave, the periodic motion of uprush and backrush occurred in the course of wave propagation on the slope. Since the large volume of water was carried shorewards by wave breaking, the longshore currents were fed continuously, and eventually returned through the surf zone because of gravity, usually as undertow. Again, in Figure 18, we present the time-averaged velocity profiles below the trough level in cases of W0H07~W5H07. Given that they share similar behaviors in other cases, the numerical results of three cross-sections ranging from the slope toe to the location in the proximity of the breaking location were used for comparation. In Figure 18a, even under different wind speeds, the magnitude of horizontal velocity profile was actually constant over the depth between the trough level and upper bottom boundary layer, consistent with Bakhtyar et al. (2009), although the latter had no wind, and maintained a larger level in the cross-sections before breaking (e.g. $x=1 \mathrm{m}$ and $x=2 \mathrm{m}$), whereas it decreased sharply being nearly close to zero in the slope toe cross-section ($x=0 \mathrm{m}$), where the undertow was less significant. Meanwhile, Figure 18b indicates that we cannot ignore the effect of the vertical velocity profile on undertow. The level of that tended to shift in the negative direction and the undertow intensity rose gradually with the distance away from the slope toe. We noticed the dependence of undertow on the wind speed, which manifested specifically as: The wind affecting the undertow intensity in various regions to different degrees. In the vicinity of the slope toe cross-section, we found that the undertow intensity itself was small whether there existed wind or not. The closer it was to the surf zone, the more intense the undertow would be as the wind speed increased. Considering the formation mechanism of undertow as well as the wind-induced setup added to the wave setup, which is responsible for the intensified imbalance between the setup pressure gradient and the wave radiation stress, and thus contributing to the acceleration of undertow for compensating the loss volume of water carried shorewards incurred by the effects of wave breaking and wind. The mean flow rate per meter width below the trough level over $6T$ is given in

Table 2. The mean flow rate below the trough level.

x (m)	hx (m)	W (m/s)	Mean Flow Rate, q (m3/s/m)
0	0.4	0, 3, 5	−0.0011, −0.0014, −0.0016
1	0.3	0, 3, 5	−0.0051, −0.0054, −0.0056
2	0.2	0, 3, 5	−0.0024, −0.0029, −0.0033

, the maximum is generated at $x=1 \mathrm{m}$, where the depth is larger with a stable return flow velocity, and wind enhances the mean flow rate

The moment of skewness and kurtosis was used to quantitatively evaluate the undertow intensity. The skewness represents the deviation of velocity from the standard Gaussian distribution, while the kurtosis indicates the peak level. The empirical formulas of skewness $C_{ss}$ and kurtosis $C_{ks}$ are known as:

$$C_{ss}=\frac{\frac{1}{n}{\displaystyle \sum _{i=1}^n{(x_i-\overline{x})}^3}}{{\left(\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{(x_i-\overline{x})}^2}}\right)}^3}$$

(20)

$$C_{ks}=\frac{\frac{1}{n}{\displaystyle \sum _{i=1}^n{(x_i-\overline{x})}^4}}{{\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{(x_i-\overline{x})}^2}\right)}^2}$$

(21)

where $\overline{x}$ is the mean of sample $x_i$. Then, the skewness and kurtosis with depth at three cross-sections ($x=0, 1, 2 \mathrm{m}$) associated with the horizontal and vertical velocity are plotted in Figure 19. The horizontal skewness was above the zero line, and it gradually tended to be right skew with the increase of wind speed except at the slope toe where it shifted left skew. Meanwhile, the kurtosis among the three cross-sections was lower than 3, while the significant difference was mainly reflected in the vertical direction, the horizontal kurtosis of various locations was at nearly the same level, it can also be observed that the influence of wind on the kurtosis was negligible.

6. Conclusions

In this paper, the effect of steady wind on the 2nd Stokes wave propagation on the slope with Waves2FOAM was investigated, wherein specific aspects were classified into three parts: Wave breaking, turbulent flow, and undertow. The main findings are as follows:

Under the same incident wave height, wind induced earlier wave breaking while the breaking location migrated towards the offshore direction. Moreover, within the same wind speed coverage constraint, the smaller the ${\xi }_0$ was, the higher the dispersion of wave breaking locations would be. As the critical breaker index of $H_b/h_b$ increased with wind speed, the critical breaker index without wind was no longer suitable for the windy cases any more.

We examined the wind effect on turbulent flow in the wave breaking processes. Wind enhanced the coherence peak and thus the turbulence intensity of wave in the low-frequency range ($f<2.5 \mathrm{Hz}$). The spectral density was larger than that with the absence of wind in the high-frequency range ($f>2.5 \mathrm{Hz}$), which could serve to account for the energy transfer processes of wind input. Meanwhile, both the TKE and TDR displayed periodic variation before and after wave breaking, in which the TKE showed a bimodal distribution. The TKE in regions near the water surface was considerably larger with the wind speed and average level rapidly decreasing with the depth. The TDR with depth was barely impacted by the wind before breaking, whereas a significant effect of wind was found after breaking. The turbulence penetrated to the deeper flow and the average level generally rose.

The dependence of undertow on the wind speed was found in the study. The closer it was to the surf zone, the more intense the undertow would be as the wind speed increased. The wind-induced setup added to the wave setup, which could aggravate the imbalance between the setup pressure gradient and the wave radiation stress, and contribute to the acceleration of undertow for compensating the loss volume of water carried shorewards incurred by the effects of wave breaking and wind.