Material Transport under a Wave Train in Interaction with Constant Wind: A Eulerian RANS Approach Combined with a Lagrangian Particle Dispersion Model

The interaction of a developed train of gravity deep water waves with suddenly applied winds is investigated in this manuscript. The direction of the wind is the same as that of the wave train (i.e., following) and its imposed surface shear stress is constant and steady. The focus of this study is on a micro-scale water wave field where the time scale is on the order of ten wave periods and the length scale is on the order of ten wave lengths. Accurate 2D Reynolds-averaged Navier–Stokes (RANS) multi-phase simulations of Navier–Stokes equations are performed in a Eulerian framework to capture the flow features, induced by the wave field and the surface wind. The sufficiently large spatial domain in the horizontal direction, combined with the sufficiently long simulation time, permits the development of surface currents and the consequent formation of a near-surface shear layer. The interaction of surface currents with the wave orbital velocity field results in the generation of spilling breaking waves. Downstream of the domain, vertical turbulent structures are observed as the result of such breaking waves. Lagrangian particle tracking is performed, using the RANS simulation velocity and eddy diffusivity data. A second order random acceleration particle tracking method is applied with the vitally important spatial gradients of the eddy diffusivity (Journal of Marine Science and Engineering, 2018, 6, 10.3390/jmse6010007.) also included in calculations. The spatial gradients of the eddy diffusivity were proved to be a key factor in material transport simulations. Our particle tracking results exhibit strong vertical mixing downstream of the domain and by means of visualizing the spiral trajectory of neutrally buoyant particles. Such enhanced vertical mixing (caused by horizontal winds) is the result of the strong near-surface advection (induced by currents) and the turbulence (induced by breaking waves). The objectives of this paper are twofold. Firstly, a numerical approach to simulate wind breaking waves is proposed based on: using K − e RANS model to capture turbulence features, employing the Volume of Fluid Method (VOF) to model the free surface flow, and applying a constant shear stress body force at the interfacial cells to simulate the wind force. Such treatment of the winds eliminates the need for fully resolving the air phase. The computed eddy viscosity profiles are in good agreement with the experimental profiles reported in the literature (νt = −κu∗z, Journal of Physical Oceanography, 1977, 7, pp. 248–255; Journal of Physical Oceanography, 1984, 14, pp. 855–863). Secondly, effects of the horizontally applied wind on the vertical mixing and eddy viscosity profiles on the water column are studied. It is observed that, away from the surface and outside the shear layer, the negative horizontal gradient of eddy diffusivity (induced by the dampening effect of breaking surface waves), combined with the downwards advection velocities (induced by breaking waves), results in an enhanced vertical mixing and reduced horizontal drift of transported material. Fluids 2018, 3, 40; doi:10.3390/fluids3020040 www.mdpi.com/journal/fluids Fluids 2018, 3, 40 2 of 18


Introduction
The interaction of surface gravity waves and surface winds results in some of the most important oceanographic processes such as: breaking waves, material transport and vertical mixing, Langmuir cells, momentum and scalar transfer at the air-sea interface and coastal upwelling [1][2][3]. Strong winds generate shear layers in the upper mixing layer, which interact with the wave orbital velocity field [4,5]. In shallow waters, where the orbital wave velocities are strong enough to affect the bottom shear layer, the bottom boundary can be affected by strong surface wind. There is a wealth of knowledge, gained by field measurements, laboratory experiments, or numerical simulations, about the effects of the wind induced surface shear current on the bottom boundary layer in shallow waters. In a broader context, the theory of near-surface response of wave fields to wind currents is quite well developed with the key feature being reported as the vertical variations of eddy viscosity. The earliest of such studies in the literature is the classical wind-driven currents by Eckman [6] in which the eddy viscosity was assumed constant across the entire vertical column. Thomas [7] proposed a more realistic model for eddy viscosity: ν t = −κu * b z b where κ is von Kármán's constant and u * b is the bottom surface friction velocity, given by u * b = (τ b /ρ) 1/2 , where τ b is bottom shear stress and ρ is water density. Madsen [8] proposed a modified equation for eddy viscosity, which was based on the surface friction velocity u * s given by ν t = κu * s (h − z b ), where h is the water depth. The model was based on the assumption that eddy viscosity increases linearly with the depth. Using acoustic bubble measurements, Thorpe [9] reported that eddy viscosity profile is of the form proposed by [8], except for a near-surface zone, where it is directly affected by breaking waves. Such near-surface behavior will later shown to be the case in this manuscript.
As far as the vertical material or properties transport in such flows are concerned, eddy viscosity was shown [10,11] to greatly affect the transport phenomenon and was reported as the key feature in such wave and wind interactions. Extensive studies have dealt with the problem of deep water wave and wind interactions [4,12,13]; however, the problem of air-sea interaction continues to be a key unresolved issue for over half a century [14]. This is particularly true in regard to the numerical simulation studies and given the two-stage, unsteady, and abrupt characteristics of the breaking phenomena. In an ideal setup and to fully represent the resulting breaking waves, the air phase must be fully resolved prior to being imposed on a train of waves. To eliminate such an expensive requirement, a new numerical method is proposed here. Such a method is based on applying the wind shear stress as a body force at the air-water interface.
Although recent studies have been successfully implementing Direct Numerical Simulation (DNS) to deep water breaking waves in the absence of winds [3,[15][16][17], there are very few DNS or Large Eddy Simulation (LES) studies in which wind breaking waves in deep water are simulated [18,19]. Some studies model air flows over idealized waves [20,21], while others perform very accurate Direct Numerical Simulations over very small domains and short periods [22]. Pizzo et al. [15], as the closet example to this current work, studied the Lagrangian drift of particles under dispersive focusing breaking waves (with no winds involved), using their accurate DNS data. The present study, on the contrary, uses the less accurate RANS method to model a more realistic wave flow situation and without the irrotational flow assumption. To assess the material transport in this work, Lagrangian particle tracking is performed on the previously resolved RANS data. Given their simplicity, Lagrangian approaches are very suitable for studying geophysical and environmental processes such as: material transport in oceans [23,24], fate of air pollutant particles [25], and oil slick movement during oil spills [26]. However, the effect of transport induced by smaller turbulent scales should be properly implemented when applying Lagrangian particle tracking methods. Otherwise, inaccurate results could be obtained, as was reported in [27,28], where the effect of the inclusion of the gradient of diffusivity in the Random Walk Method was studied. It was demonstrated [27,28] that, by not including the gradient of diffusivity in the calculation, particles accumulate in the areas with lower diffusivity. This problem should be handled by direct calculation of the spatial gradients of eddy diffusivity tensor components and its inclusion in the random process [29][30][31][32]. It is shown here that the negative horizontal gradient of eddy diffusivity, induced by the damping effect of the breaking waves, is the key factor in the enhanced vertical mixing under wind breaking waves. This study is the continuation of our previous work [27] where the effect of non-breaking waves on oil droplet transport was studied.

Eulerian RANS Framework
For an unsteady, viscous, incompressible, two-dimensional flow, the Reynolds averaged Navier-Stokes equations are given as: ∂ρK ∂t ∂ρ ∂t where Equations (1) and (2) represent the conservation of mass and conservation of momentum in two spatial directions, respectively. Equations (3) and (4) represent the turbulence k-epsilon model consisting of a transport equation for turbulent kinetic energy (TKE) (Equation (3)) and a transport equation for TKE dissipation (Equation (4)) [33]. In the equations above, i = 1, 2 with index 1 corresponding to streamwise and 2 to wall-normal directions, respectively. A repeated index indicates summation over the index. Here, t is time, u i is the Reynolds-averaged fluid velocity vector, P is the Reynolds-averaged pressure, f i is the body force term, ρ is constant density of the fluid, µ is dynamic viscosity, K is TKE, is TKE dissipation rate, P k is TKE production rate due to the mean velocity shear, α k and α are the inverse Prandtl numbers for K and , respectively, C 1 = 1.42 and C 2 = 1.68 are model constants, and µ t is eddy viscosity. Note that the eddy viscosity tensor is taken to be diagonal with equal diagonal entries, thus we consider an isotropic eddy diffusivity. The K − model's isotropic eddy viscosity is given by with C µ = 0.0845. The last term on the right-hand side of Equation (4) is specific to the Re-Normalisation Group (RNG) K − model [34] and is given as with η ≡ SK/ with η 0 = 4.38 and β = 0.012 being model constants. Note that S = (2S ij S ji ) 1/2 with To model the multi-phase free-surface flow, the Volume Of Fluid (VOF) method is used. The two-phase fluid is composed of air and water. The governing equations (described above) for both of the phases are the same. However, the density, dynamic viscosity, and eddy viscosity depend on the phase in the local cell and may vary between adjacent cells. At each cell, the dynamic viscosity and density are calculated by: where α w is a scalar value representing volume fraction of water with value of 1 corresponding to a full water cell and 0 corresponding to an air cell. An extra transport equation for α w is solved to track the interface during the simulation [35].

Lagrangian Particle Dispersion Framework
The governing advection diffusion in a Eulerian framework [36] can be written as where C is the Reynolds-averaged tracer concentration (e.g., oil with C being oil mass per unit volume of water), D i is the turbulent eddy diffusivity in the corresponding directions and, given the difficulties associated with numerical solutions of highly advective flows [37], the solution of Equation (8) can require artificial numerical diffusion leading to inaccurate concentration solutions. By the passive tracers assumption (i.e., the oil particle motion has no effect on the dynamics of the carrier fluid phase), the oil is treated as a discrete phase that is only dispersed and diffused by the flow. Such flow induced particle advection and diffusion are studied using the Random Walk Method. The method consists of time integration of the Lagrangian velocity and acceleration equations for individual particles (i.e., Random acceleration method). To track a particle located at the position x (n) i and at the starting time t (n) , the location at the time t (n) + ∆t is found by the following pair of stochastic equations [38] needing to be solved: where R is a random Gaussian distribution. The second term on the right-hand side of Equation (10) represents the advection induced by the carrier velocity field. The third term is the divergence of the eddy diffusivity whose significance will be discussed later in the manuscript. The last term is a stochastic model representing the fluctuating turbulent field constructed from the RANS simulation data. Using the Boussinesq hypothesis and assuming that the eddy diffusivity is isotropic [33,39,40], we drop its index and model the eddy diffusion coefficient as with ν = µ/ρ being the kinematic viscosity of fluid. Note that for a high Reynolds number flow, the eddy viscosity is orders of magnitude larger than dynamic viscosity and hence the latter could be neglected. Furthermore, in (9), u p and ρ p are particles velocity and density and the drag force acting on the particle by the surrounding fluid id simulated by F D (u i − u p ) given by: where d p is particles diameter, C D the drag coefficient for spherical objects, and R e the relative Reynolds number, defined as R e = ρd p |u i − u p |/µ.

Flow Configuration
The domain of simulation consists of a rectangular box of 20 m × 3 m in x 1 and x 2 directions, respectively (see Figure 1). A mesh comprising 835,229 nodes and 834,239 mixed quadrilateral cells was used for the simulation. Given the sensitivity of the VOF method to coarse meshes near the interface area, the mesh was refined significantly near the free surface. The refinement was performed within a rectangular region. The size of the rectangular region was 8 × 0.3 m in the x 1 and x 2 directions, respectively, with the middle of the region corresponding to the mean water level. Within this rectangular region, the mesh size varied between 0.001 m and 0.005 m. Outside of this region, the mesh size varied between 0.001 m and 0.01 m. Flow simulations were performed using the commercial package ANSYS FLUENT ( ANSYS R Academic Research, Release 15, Canonsburg, PA, USA) [35]. The numerical discretization was based on the Finite Volume Method employing the VOF technique to capture the air-water interface and the RNG K − turbulence model. The Finite Volume Method was used with second-order upwind schemes to discretize the equations in space. The implicit and first order SIMPLE method is employed for the time discretization, given its stability and fast convergence [35]. The particularly small chosen time step of 0.0005 s to capture the unsteadiness of the flow [41] and to fulfill the CFL number conditions and restrictions. The solution was started from the rest and with an initial water level equal to h. At the input boundary (the left edge in Figure 1), a Dirichlet velocity boundary condition is imposed. The velocity value at the input boundary is calculated using Equation (14). Standard wall models are applied at the right and bottom edges of the rectangular domain, while the pressure outlet boundary condition is applied on the top edge.

Wave Formulation
The simulation was started initially with the modeling of non-breaking deep water waves. The setup is depicted in Figure 2 with the characteristic wave length λ being the horizontal distance between the crests (or troughs), wave height H is the vertical distance between the crest and trough of the wave, and the wave period is T. The wave number is defined as k = 2π/λ and the angular frequency is given by σ = 2π/T. For deep water waves, the wave length is linked to the frequency using the dispersion relation [42] as where g is acceleration due to gravity. The analytical solution for velocities using the first order linear wave theory is given as where h is the mean water depth. The chosen wave for simulation has a period of T = 1 s and a height of H = 0.15 m. The wave length calculated by the dispersion relation in (13) is equal to λ = 1.56 m.
To ensure a deep water wave condition, the mean water depth was set equal to h = 1.2 m. Furthermore, water density was taken as 998.2 kg/m 3 , air density as 1.225 kg/m 3 , water dynamic viscosity as 1.003 × 10 −3 kg/ms, and air dynamic viscosity as 1.7894 × 10 −5 kg/m −s . Figure 2. Details of multi-phase wave simulation. The areas colored in red represent cells with water phase (i.e., α w = 1 in Equation (7)) and areas colored in blue represent cells with air phase (i.e., α w = 0). Note that the air-water interface is characterized by a thickness associated with the Volume Of Fluid (VOF) method.

Wind-Wave Interaction Simulation
The constant wind that blows in the positive x 1 direction is simulated by adding a volume body force to the momentum Equation (2) as given by the following equation: whereκ is the gradient of the normal vector at the interface, ρ is the cell density given by Equation (7) and u * is the friction velocity based on the surface shear stress. To assess the response of near-surface current profile to the wind strength and also to study the resulting vertical mixing, two different cases with different winds are simulated. For the weaker wind, the surface body force is 0.5 [N/m] and, for the stronger wind, the surface body force is 0.6 [N/m]. Friction Reynolds number can be calculated by Re τ = u * δ/ν with an average friction velocity for each wind case calculated using Equation (15) and δ the turbulent length scale. Choosing the shear layer thickness (≈0.1 as shall be seen) as the turbulent length scale, our calculated Reynolds numbers for each wind case are 2280 and 2500 for the developed flow.

Particle Tracking
After simulating the flow, the results were imported into a particle tracking code. The imported results included velocity components, eddy viscosity or diffusivity and water volume fraction α w . The code is capable of constructing a triangular unstructured mesh over any set of points and using linear interpolation to calculate the in-between values of the variables. The particle search algorithm locates the particle and links it to the corresponding triangle in which the particle has traveled to. Moreover, the transient input data are updated at each tracking time step. Using an unstructured linear mesh provides the capability of constructing the spatial gradients of eddy diffusivity. Often overlooked, the gradient of eddy diffusivity is particularly important in mass transport studies and must not be neglected [27,28]. Particle tracking in general can be performed for light and buoyant (e.g., oil droplets, bubbles) or denser and submerging (e.g., sediment ) particles. To avoid the complication brought about with buoyancy and inertia effects, given that the main goal is to compare two cases with different wind strengths, we choose to conduct the experiment on a group of 1000 of neutrally buoyant particles. Furthermore, the tracking was performed over a 10 wave period time span after the wind was activated. The particles were uniformly distributed and seeded at two different streamwise locations; the first particle group was distributed in the spatial domain of x 1 ∈ [4.5, 5] m and the fixed elevation of x 2 = 1.15 m, while the second group occupied the spatial span of x 1 ∈ [5.5, 6] m while sharing the same vertical location with the first group.

Numerical Simulation Validation
The solution procedure in this study comprises the application of a sudden constant following wind on a developed regular wave train. The simulation of the regular wave was started from the rest. The simulation was continued until a fully developed flow was attained. To collect the essential data to perform the test particle tracking, the simulation was continued for another 10 wave periods. These numerical results were compared with the analytical solutions, obtained by the linear wave theory [42]. Our previous paper [27] showed the validity of the simulations and therefore such results are not repeated here. To simulate wind breaking waves, the wind body force was added on and the simulation is conducted for another 10 wave periods. Note that, as two different winds scenarios are simulated, two different simulations are conducted with the same initial conditions and different wind body forces.

Non-Breaking Wave Simulation
In our previous work [27], the validity of non-breaking wave simulations were demonstrated by providing the plots for: comparison of the time series of velocity components with the analytical solution from second order linear wave theory, and the comparison of vertical profiles of velocity components with their analytical profiles. Furthermore, the calculated averaged Stokes drift velocity (based on the particle tracking results) showed good agreement with the theoretical value (based on Stokes drift theory). To present the following results, the wind force was applied to the previous non-breaking wave simulations results.

Wind-Wave Simulation
As was mentioned earlier, the following wind was suddenly applied to a train of developed deep water waves. The interaction of the surface winds and the orbital velocity of waves results in small spilling breaking waves as shown in Figure 3. By choosing the characteristic wave slope S = ak, our calculated initial wave slope is 0.302. As the result of wind-induced surface current, the wave steepens slightly until the crest becomes unstable. The slight increase of the steepness, however, is not significant enough to push the breaking wave type into the more extreme plunging or collapsing scenarios. Based on our observations, the fourth wave in the train exhibits the most severe cases of breaking and the fifth wave suffers height loss and energy dissipation (which results in its flattening). This is figuratively consistent with the observations of [17], where the total dissipation of the wave in the first four periods were reported (although the breaking wave generation mechanism in that study was completely different and was based on the Dispersive Focusing Method). The dissipation of the fifth wave is evident in Figure 3. The figure illustrates the water volume fraction α w at the time t = 7 [s] after the winds were started. The plot is zoomed in at the horizontal span of x 1 ∈ [6, 7.6] m (where the fifth wave is located and also where an enhanced vertical mixing will be observed in Figure 4). Furthermore, velocity vectors are superposed to the contours in order to illustrate the relative strength of advective currents for each wind case. The figure illustrates formation of the surface currents which are stronger near the crest (x 1 ∈ [7, 7.6] m, x 2 ∈ [1.1, 1.3] m) and with higher magnitude for the stronger wind (bottom panel). Furthermore, this figure shows an area with strong downward vertical velocity vectors (x 1 ∈ [6.6, 6.9] m, x 2 ∈ [0.9, 1.3] m). The downward velocity magnitude is higher for the case with weaker wind. The lower velocity magnitude for the case with stronger wind can be related to the stronger near-surface current that affects the vectors by means of energy dissipation. Moreover, the breaking occurs as the result of the collision of upwards orbital velocities, induced by the waves, and the positive horizontal current flow, induced by the wind. The aforementioned collision phenomenon is evident in the area x 1 ∈ [7.4, 7.6] m, x 2 ∈ [1.2, 1.3] m. Given the stronger near-surface current for the case with stronger wind, stronger breakers are observed for this case. Figure 5 shows the wave-length and time averaged profiles of wall-normal eddy viscosity. Time averaging is performed over time span of the last three wave periods and wave-length averaging is performed over the third (occupying the streamwise span of x 1 ∈ [3.5, 5]) and the fourth (occupying the streamwise span of x 1 ∈ [5, 6.5]) waves in the train. The plot corresponds to both of the cases with weaker (solid and dotted blue lines) and stronger (solid and dotted red lines) winds as well as the cases with no wind (regular waves). To demonstrate that the observed dissipation is indeed a direct result of wind-wave interactions and not a numerical induced dissipation, the profiles of eddy viscosity corresponding to both the third and fourth waves for the case with no wind are also included (solid and dotted green lines).   Based on Figure 5, downstream of the domain where the turbulence is developed, the wall-normal profiles of time averaged eddy viscosity show close agreement with Madson's relation [8] and for both the wind cases. Based on Madson's relation [8], the eddy viscosity is assumed to be proportional to the depth where κ is von Kármán's constant (≈0.4) andx 2 is the vertical component of a coordinate system which originates at the free surface and is positive downwards. By transforming such vertical coordinate to our chosen Cartesian system (which origins at the bottom wall and is positive upwards), we can write: Note that η is the free surface elevation measured from the mean water level. As the eddy viscosity profile is presented in a time averaged manner, it is fair to assume that the averaged value of < η > x 1 ,t should be approximately equal to 0. Given this, the averaged profiles should start at the vertical axis at x 2 = h −ζ whereζ is the thickness of the shear layer, the layer which is directly affected by breaking waves. Using acoustic bubble measurements, Thorpe [9] showed that eddy viscosity profile is of the form (16), except for a near-surface zone where the flow is directly affected by breaking waves, with the estimated thickness of this zoneζ to be approximately 10 times the (g/σ 2 )H value. Given our simulated wave data, the calculated shear layer thickness,ζ is 0.2236 [m], which is based on the location where the eddy viscosity reaches its maximum value and its vertical gradient changes sign. This value matches the "medium profile" of Jenkins [43] (ζ = 6(g/σ 2 )H). For the matching profiles (downstream and shown in solid red and blue lines) in Figure 5, the calculated shear layer thickness is approximately the same as the theoretical value. For illustration purposes and to show the slope match, we had to slightly shift the profiles in Figure 5. Furthermore, the gradual decrease of eddy viscosity value in streamwise x 1 direction is the direct result of the dissipation induced by the breaking waves. Such drastic decline in eddy viscosity values produces a negative horizontal gradient of eddy diffusivity, which, as will be shown, causes enhanced mixing and greatly affects material transport. Downstream of the domain (x 1 = 7 [m]), the eddy viscosity reaches the equilibrium value (solid red line) and hence matches a theoretical value given by (16).

Drift Current Formation
Our results attest to the formation and further development of near-surface drift currents as the direct result of applied winds. Figure 6 presents the instantaneous snapshots of streamwise velocity component (u 1 [m/s]) at t = 5 [s] (left-hand side panels) t = 10 [s] (right-hand side panels). The top panels correspond to the case with the weaker wind and bottom panels correspond to the case with the stronger wind. To illustrate details of the surface current, the plot is zoomed in a domain of x 1 ∈ [5, 6.5] m and x 2 ∈ [0.8, 1.4] m. As it can be seen, the evolution of drift currents in time results in the negative velocity areas (x 1 ∈ [6, 6.5] m, x 2 ∈ [0.8, 1] m) shrinking and the velocity magnitudes decreasing. On the other hand, the positive velocity areas (x 1 ∈ [5, 6.5] m, x 2 ∈ [1, 1.2] m) start to grow and the velocity magnitudes start to increase. The effect is more pronounced for the case with the stronger wind (bottom panels). Furthermore, by comparing the wave heights between the cases with different winds, it can be concluded that the stronger wind blows the waves flat [44].  To assess the drift current velocity, Figure 9 is presented. The figure shows the horizontal structure of the drift currents during a wave period and in three streamwise locations. Based on this figure, upstream of the domain (the top panel corresponding to x 1 = 3 [m]), the surface current has the maximum local velocity (approximately 1 [m/s] for the weaker wind and 1.2 [m/s] for the stronger wind). Towards downstream, the streamwise velocities at the surface start to decrease (as the result of breaking waves and the ensuing dissipation) and the negative wall-normal velocity components start to increase. Note that the high local values at the interface are the result of lower density at the interface and hence the average current velocity should be calculated by integrating the local values in wall normal direction and only after filtering the wave velocity. Finally, the observed high horizontal velocities at the surface greatly affect the light material that can float on water surface (e.g., oil) and greatly increase their drift motion, as was shown in [15] where added drifts as large as 10 times the Stokes drift were reported for the particles at the surface. As it will be seen in the upcoming section, our particle tracking experiments take place deeper in the water column and the particles therefore are seeded and tracked outside the surface shear layer. The primary focus in our studies is on the vertical mixing and to demonstrate the drastically different effect of turbulence on the material transport in such regions.

Particle Tracking and Vertical Mixing
In this section, Lagrangian ensemble-averaged particle trajectories are presented in order to better understand the combined effect of advective velocities (induced by the breaking wave motion) and the diffusion (induced by turbulence). Results are obtained by tracking 1000 individual neutrally buoyant particles and calculating the plume trajectory by averaging the trajectories over all of the particles. The chosen number of particles was sufficiently large to ensure that the results are not affected by a further increase of the number of particles. Figure 4 presents the ensemble averaged trajectories for particles released at two different streamwise locations, x 1 ∈ [4.5, 5] m, x 1 ∈ [6.5, 7] m while all the particle sets had the identical wall-normal location of x 2 = 1.15 m. The chosen location is outside the surface shear layer as was observed in Figure 5. As was mentioned, the primary focus is on the vertical mixing. The dashed-dotted lines in the graphs correspond to the trajectories due to the pure advection and their values are calculated by setting eddy diffusivity and its gradients equal to 0 (i.e., the last two terms in the right-hand side of Equation (10)). Comparing the wind-wave induced particle trajectories in Figure 4 with the classical wave induced drift motion trajectories spirals, one realizes that the horizontally blown wind induces enhanced vertical mixing. Such mixing is induced through the downward strong advection in areas where wave breaking occurs (upstream of the domain) and through the strong turbulence diffusion prior to breaking (downstream of the domain). Downstream of the domain the turbulence diffusion is not affected by the breaking waves. Thereby a thicker near-surface shear layer is formed. Such layer contains higher eddy viscosity/diffusivity values (evident in Figure 5, blue lines). Comparing the particle trajectories in upstream (left top and bottom panels) and downstream (right top and bottom panels) of the domain, one realizes that in the upstream region both the drift motion and vertical mixing affect the trajectories. Downstream of the domain, however, vertical mixing is the dominating force that results in enhanced mixing and negative drift motion for the stronger wind case (bottom right panel). For the case with the stronger wind, sudden drops in trajectories are associated with the stronger breakers on the surface. As a result of wind breaking waves, turbulent diffusion plays a more significant role in the material transport. This is evident in all the bottom panels of Figure 4 where there are clear different trends of particle movements when the cases with pure advection (green lines) are compared with the cases with advection/diffusion (blue lines). To better understand the effect of turbulence diffusion, we next present the averaged wall-normal profile of eddy diffusivity. Figure 10 presents the vertical profile of streamwise-time averaged eddy diffusivity profiles for the case with wind body force f = 0.5 [N/M] (blue labeled as "Weaker wind"), the case with the wind body force f = 0.6 [N/M] (red labeled as "Stronger wind") and the case with no wind stress and with non-breaking waves. Interestingly, the streamwise-time averaged eddy diffusivity matches a theoretical line similar to that of eddy viscosity reported in ( [8]).
The surface shear layer is significantly thinner for the case with pure waves and no wind. The thinner boundary layer results in much sharper gradients of eddy viscosity/diffusivity at the surface. Sharp gradients of the aforementioned quantities result in vertical diffusion mixing near the surface for the case with no wind. However, as it was mentioned before, a combination of two factors results in stronger mixing for the cases with blowing winds-firstly, the strong downward advection downstream of the domain, which is induced by the breaking waves; secondly, the thicker near-surface shear layer with stronger turbulence activities upstream of the domain, observed in Figure 5. The negative drift downstream of the domain for the case with stronger wind (right bottom panel of Figure 4) is to be explained as the result of the combined effect of the weakened downward wave orbital motion and the negative stream-wise gradient of eddy diffusion ∂D/∂x 1 (as the result of wind induced dissipation that is higher upstream) manifesting as a negative streamwise advective velocity.

Conclusions
Material transport under the wind breaking waves was investigated. The present study was focused on micro-scale spatial and temporal scales compared to our previous work [27], which focused on material transport under non-breaking waves without the presence of wind. Based on our current results, in the water column and outside the shear layer, mass transport is significantly affected by the wind and through the breaking waves on the surface. A combination of two different mechanisms, namely, the enhanced downwards advection induced by the interaction of the breaking waves with the surface current, and the diffusion that is induced by turbulence, results in enhanced mixing in the water column. A key issue explored in our previous work [27,[45][46][47][48] was vertical variations of eddy viscosity and the effect of the inclusion of the gradient of eddy diffusivity in the Random Walk particle dispersion practices. While studying the vertical distribution of particles in the water column induced solely by turbulence diffusion, Visser [28] described the mechanism by which the center of mass of the particles move in the direction of increasing diffusivity. As a result, material is transported from areas of lower eddy diffusivity to areas of higher diffusivity.This mechanism was once again proved to be a very important factor in the material transport studies and therefore should not be overlooked.
The results of the current study also agree with our previous findings and confirm the necessity of the inclusion of the gradient of eddy diffusivity in the Random Walk Method. Another important result presented in this paper was the use of the Volume of Fluid method to simulate the wind flow by applying a wind body force in interfacial cells (e.g., Finite Volumes, Finite Elements, etc.) The proposed approach results in a significant decrease in terms of degrees of freedom by eradicating the need to fully model the wind in the air phase. The proposed method led to computation of vertical eddy viscosity profiles, which show good agreement with those from the literature that are based on experimental studies. Our result also manifested the evolution of a surface current, which, by its interaction with the wave motion, results in breaking surface waves. Breaking waves result in the development of an upper mixing layer with the estimated thickness reported in the literature [9,43]. Our calculated thickness shows agreement with those from the literature. To relate the wind stress to the wind velocity for practical applications, given the well-developed theory of response of near-surface current to wind, the appropriate values of surface shear velocity u * could be calculated based on the desired wind speed. Finally, since some of the flow features can not be captured in a 2D simulation, it is vitally important to extend this work to a 3D framework. The proposed method can be used in a 3D framework to study a variety of oceanographic phenomena such as generation of Langmuir cells and wind-induced coastal upwelling. The approach could be applied to shallow water waves to assess the effect of wind on increasing the bottom shear stress and linking it to disruption of the log layer. Turbulence greatly affects material transport; recently, Ghodke et al. [49,50] in their groundbreaking and fully resolved numerical study reported the effects of turbulence on incipient material transport mechanisms. Such high-fidelity simulations can further provide time-accurate eddy viscosity data that otherwise are extremely challenging to obtain using experiments.
Author Contributions: R.G. and N.S. conceived and designed the experiments follwing the ideas from R.G.'s previous paper, which can serve as part 1 to this paper; Nityanand Sinha performed the experiments and post processed the data. R.G. wrote the paper.