LARGE-WAVE SIMULATION OF TURBULENT FLOW INDUCED BY WAVE PROPAGATION AND BREAKING OVER CONSTANT SLOPE BED

In the present study, the three-dimensional, incompressible, turbulent, free-surface flow, developing by the propagation of nonlinear breaking waves over a rigid bed of constant slope, is numerically simulated. The main objective is to investigate the process of spilling wave breaking and the characteristics of the developing undertow current employing the large-wave simulation (LWS) method. According to LWS methodology, large velocity and freesurface scales are fully resolved, and subgrid scales are treated by an eddy viscosity model, similar to large-eddy simulation (LES) methodology. The simulations are based on the numerical solution of the unsteady, threedimensional, Navier-Stokes equations subject to the fully-nonlinear free-surface boundary conditions and the appropriate bottom, inflow and outflow boundary conditions. The case of incoming second-order Stokes waves, normal to the shore, with wavelength to inflow depth ratio λ/dΙ  6.6, wave steepness H/λ  0.025, bed slope tanβ = 1/35 and Reynolds number (based on inflow water depth) Red = 250,000 is investigated. The predictions of the LWS model for the incipient wave breaking parameters breaking depth and height are in very good agreement with published experimental measurements. Profiles of the time-averaged horizontal velocity in the surf zone are also in good agreement with the corresponding measured ones, verifying the ability of the model to capture adequately the undertow current.


INTRODUCTION
Wave breaking, one of the major near-shore processes, is responsible for the development of wavegenerated currents that result into the initiation and development of sediment transport.Therefore, the investigation of wave breaking, as a phenomenon interconnected with a series of significant coastal processes, is of great interest.Wave breaking takes place when wave height and steepness become very large as the water depth becomes shallower.For the case of spilling breaking, a vortex structure, usually called "surface roller", is formed under the collapsing wavefront just after breaking.
As mentioned above, closely related to wave breaking is the generation of wave-induced currents, i.e., the longshore current, which is developed only when the direction of the breaking wave is oblique to the shoreline, and the cross-shore current, which is known as the undertow current.Both of them are developing in the surf zone, i.e., the coastal area where wave energy dissipation occurs after breaking.The undertow current owes its existence to the mean shear stress field, developing in the surf zone, in order to balance the pressure gradient and the momentum fluxes due to the wave set-up and the wave height dissipation.Close to the bottom, the current is offshore directed, while near the free surface is oriented towards the shore, ensuring that the total cross-shore water flux is zero.
The wave propagation and breaking in the coastal zone has been investigated by several researchers for many years, by use of numerical and physical models.As for the numerical simulation of spilling wave breaking, there are two broad categories of models based on the treatment of the surface roller: those that incorporate empirical models, often called surface roller (SR) models, and those that simulate the surface roller as part of the turbulent flow, induced by the wave propagation, utilizing one of the turbulence modeling methods, i.e., the Reynolds-Averaged Navier-Stokes (RANS) equations and the large-eddy simulation (LES) method.
The SR model, in which incipient breaking is defined by empirical criteria, may be coupled with Boussinesq (Briganti et al. 2004;Madsen et al. 1997;Schäffer et al. 1993;Veeramony and Svedsen 2000) or Euler (Dimas and Dimakopoulos 2009) equations giving very good results for the normal to the shore wave breaking, but it cannot be easily extended to the case of oblique to the shore wave breaking.RANS models, where all turbulent scales are treated by a closure model, have been applied for the case of two-dimensional turbulent flow during spilling breaking (Bradford 2000;Lin and Liu 1998;Torres-Freyermuth et al. 2007).LES models, where only the large scales of turbulence are resolved while the small ones are modeled, have been used for two-dimensional (Hieu et al. 2004;Zhao et al. 2004) and three-dimensional (Christensen and Deigaard 2001;Christensen 2006) flows.LES method requires more computational resources than RANS, but generally achieves better results for the dynamics of spilling breakers.Both methods require the use of a scheme for capturing the free-surface elevation, such as the Volume of Fluid (VOF) method and the Marker and Cell (MAC) method.
Recently, Dimakopoulos and Dimas (2011), in order to investigate the turbulent flow induced close to the free surface by the oblique wave breaking over a constant slope bed, employed the Large-Wave Simulation (LWS) method coupled with their inviscid flow solver (Dimas and Dimakopoulos 2009).Their model was validated by comparison of their numerical results to corresponding experimental measurements presented in Ting and Kirby (1994;1996), which was also used for the calibration of the LWS model.
The novelty of the present study is the coupling of the LWS methodology (Dimas and Fialkowski 2000) with a numerical solver of a three-dimensional viscous flow, in order to calculate the waveinduced currents in the surf zone, the generation of which is substantially affected by the presence of the bed shear stress.Specifically, numerical simulations of the three-dimensional, free-surface flow, induced by the propagation of nonlinear breaking waves, normal to the shoreline, over a constant slope, rigid bed, are presented.The LWS methodology is based on the decomposition of the flow variable scales (velocity, pressure and free-surface elevation), to resolved (large) and subgrid (small) scales.One of the main objectives of this work, is to overcome some of the limitations of the aforementioned studies, mainly arising by the use of inviscid (Euler) and depth-averaged (Boussinesq) models.For example, inviscid models are not able to capture correctly the profile of the undertow current, since they ignore viscous effects close to the bed, while Boussinesq models are incapable of calculating flow variable profiles over the depth.In the following sections, the flow equations, the main features of the LWS methodology, the numerical method, the simulation results and the main conclusions, are presented.

FLOW EQUATIONS
The three-dimensional, incompressible free-surface flow, for a fluid of constant viscosity, is governed by the continuity and the Navier-Stokes equations 2 1 Re where i, j = 1,2,3, t is time, x 1 , x 2 are the horizontal coordinates, x 3 is the vertical coordinate, positive in the direction opposite to gravity, u 1 , u 2 and u 3 are the corresponding velocity components, p is the dynamic pressure and Re d is the Reynolds number.Equations ( 1) and (2) are expressed in dimensionless form with respect to the inflow depth d I , the gravity acceleration g and the water density π, therefore Re d = (gd I ) 1/2 d I /ν, where v is the kinematic water viscosity.For viscous flow, the kinematic and the normal stress dynamic boundary conditions at the free surface are, respectively, where η is the free-surface elevation, Fr is the Froude number, which under the present dimensionless formulation is equal to one, while in Eq. (4) the atmospheric pressure is considered equal to zero.The shear stress dynamic boundary condition at the free surface, is expressed for each of the horizontal coordinates x 1 and x respectively.In addition, the no-slip and non-penetration boundary conditions at the bottom are respectively, where d is the bottom depth measured from the still free-surface level.
Given that the free surface is time-dependent, the Cartesian coordinates are transformed, in order for the computational domain to become time-independent, according to the expressions and where -1 ≤ s 3 ≤ 1.In the transformed domain, s 3 = 1 corresponds to the free surface and s 3 = -1 to the bottom.By application of Eq. ( 9), the continuity and Navier-Stokes equations ( 1) and ( 2) are transformed, respectively, where k= 1,2, hereafter, and

LWS methodology
As aforementioned, the large-wave simulation (LWS) method is based on the application of a volume filter to the velocity components and pressure, as in LES, and an exclusive surface filtering operation for the free-surface elevation.Therefore, each flow variable, f, is decomposed into resolved, f , (large) and subgrid, f ', (small) scales, in a manner illustrated in Fig. 1 for the decomposition of the free-surface elevation, η.The filtering operation is applied on Eqs. ( 10) and ( 11), resulting into the continuity and Navier-Stokes equations for the resolved scales of the flow, which, respectively, are where includes the subgrid scale (SGS) terms, i.e., the eddy SGS stresses and the wave SGS stresses, which, respectively, are The eddy SGS stresses appear both in the LES and the LWS method, while the wave SGS stresses appear exclusively in the LWS method.The transformation, given by Eq. ( 9), and the filtering procedure are also applied successively to the boundary conditions (3) -( 8) at the free surface (s 3 = 1) and the bottom (s 3 = -1), resulting into the transformed boundary conditions for the resolved scales.
It must be noted that the filtering procedure takes into to account some simplifying assumptions, which are specified analytically in Dimakopoulos and Dimas (2011).In addition to this, the SGS terms that result from the filtering of the viscous terms of Eq. ( 11), are considered to be negligible compared to the rest of the terms of the same equation.
In the present study, the eddy and wave SGS stresses, appearing in Eq. ( 16), are computed by use of Smagorinsky eddy-viscosity models (Rogallo and Moin 1984).Specifically, the model for the eddy SGS stresses is where C = 0.1 is the model parameter, set according to the usual practice in LES method, Δ = (Δ 1 Δ 2 Δ 3 ) 1/3 is the smallest resolved scale based on the grid size, S ij is the strain-rate tensor of resolved scales 1 2 is its magnitude.The model for the wave SGS stresses, based on the one presented in Dimas and Fialkowski (2000), is given as where C η is the model parameter and S ij η is a modified strain-rate tensor of resolved scales where δ ij is the Kronecker delta.

NUMERICAL METHOD
The flow simulations are based on the numerical solution of the transformed Navier-Stokes equations, which is achieved by use of a fractional time-step scheme for the temporal discretization, and a hybrid scheme for the spatial discretization.The hybrid scheme includes central finite differences, on a uniform grid with size Δs 1 , for the discretization along the streamwise direction s 1 , a pseudo-spectral approximation method with Fourier modes along the spanwise direction s 2 , and a pseudo-spectral approximation method with Chebyshev polynomials along the vertical direction s 3 .
The transformed Navier-Stokes equations ( 16) can be written in the form 1 Re where is the transformed pressure head and T j  is the transformed Laplacian operator.The time-splitting scheme for the temporal discretization, in which each time-step consists of three stages, achieves the calculation of the velocity field at the next time step successively, by adding the corresponding corrections of each of the three stages to the field of the previous time-step.Moreover, the dynamic pressure field is obtained in the second stage of each time step, while the free-surface elevation is calculated from the kinematic boundary condition at the end of each time step.
At the first stage of each time-step, the nonlinear term, i A , the SGS term, i T , and the viscous term, i V , of the transformed equations of motion ( 24) are treated explicitly by an Euler scheme.At the second stage, an implicit Euler scheme is used for the treatment of the pressure head term, T j  , of Eq.( 24), which results into a generalized Poisson's equation for  by satisfying the transformed continuity equation as well.The transformed dynamic (normal stress) free-surface condition and nonpenetration bottom condition are imposed at this stage.At the third stage, the remaining viscous terms,  , of Eq. 24, are treated by an Euler implicit scheme satisfying the transformed dynamic (tangential stress) free-surface and bottom conditions.According to the hybrid scheme for the spatial discretization, each flow variable f (velocities and pressure) is approximated by the following expression where mn f  is the Chebyshev-Fourier transformation of f , M is the number of Fourier modes, L 2 = M • Δ 2 , is the length of the computational domain in s 2 , and T n is the Chebyshev polynomial of order n, and N is the highest order of the Chebyshev polynomias.The transformations between physical and spectral space are performed by a Fast Fourier Transform algorithm (Press et al., 1992).
Application of Eq. ( 25) for the discretization of Eq. ( 24), leads to the formation of a system of algebraic equations, with the general form    , for each of the transformed flow variables.
The aforementioned system may be divided into M independent subsystems (   for each Fourier mode), which can be solved in parallel, due to the decoupling of the Fourier modes.Each subsystem is solved at each time-step, using an iterative generalized Gauss-Seidel method.The matrix of coefficients   m A is band diagonal, and is decomposed once at the beginning of the computation by using the LU-decomposition method.
In the present work, the propagation, transformation and spilling breaking of incoming secondorder Stokes waves over a constant slope bed is simulated.As shown in the sketch of the computational domain (Fig. 2), a flat bed region of length L I and constant depth d I , which ensures the development of the incoming waves, is followed by the inclined region of the bed.Then, a flat region of length L E and constant depth d E << d I (the formulation allows the outflow depth d E to be small but nonzero) is considered in order to simulate the swash zone of a coastal area, where the waves are completely damped.For this reason, two overlapping zones are placed in the outflow region: a wave absorption zone of length L A ≈ L E , which ensures that waves are not reflected by the outflow boundary (Dimas and Dimakopoulos, 2009), and a velocity attenuation (slowdown) zone of length L D ,.For the numerical solution of Eq. ( 24), a reduced value of Re d is used within the slowdown zone, which corresponds to an increased value of the kinematic viscosity.

RESULTS
The validation of the inviscid version of LWS methodology, coupled with the Euler equations, was performed by Dimakopoulos and Dimas (2011), and one of the main results of their work, was the calibration of the parameter C η , used in the model for the wave SGS stresses.The chosen value, C η = 0.4, resulted from the comparison of their numerical results to corresponding experimental measurements, presented in Ting and Kirby (1994;1996), for the case of spilling wave breakers, propagating normal to the shoreline over a beach of constant slope tanβ = 1/35.
In the present study, the accuracy and efficiency of the viscous version of LWS methodology, coupled with the Navier-Stokes equations, is investigated, performing numerical simulation for the normal to the shoreline propagation, transformation and spilling breaking of incoming second-order Stokes waves over a bed of constant slope tanβ = 1/35.Our numerical results are, also, compared to the experimental measurements conducted by Ting and Kirby (1994), while the suggested value of C η = 0.4 (Dimakopoulos and Dimas 2011) is adopted.It must be noted that all of the numerical results presented in this study are spanwise averaged.
The experimental flow parameters (Ting and Kirby 1994) for the case of spilling breaking are summarized in the following: wave inflow depth, d I = 0.4 m, wave height and period, H I = 0.125 m and Τ = 2 s, respectively, which correspond to wave height and wavelength at deep water, H ο = 0.127 m and λ ο = 6.245 m, respectively.In our case, the wave parameters at deep water are identical to those of Ting and Kirby (1994), but a larger inflow depth d I = 0.7 m is considered, since the Stokes wave theory is utilized for the incoming waves, which are of height H I = 0.118 m.The parameters of the incoming waves are rendered dimensionless by d I , and g, and the resulting values are H I = 0.168, Τ = 7.487 and λ = 6.605.The Irribaren number is ξ ο = tanβ(λ ο /Η ο ) 1/2 = 0.2, which corresponds to a spilling breaker of medium strength, while a value of Re d = 250,000, is considered.In the slowdown zone, a value of Re d divided by 100, is utilized.The total length of the computational domain is L = 60, the flat inflow region has length L I = 15, while the swash zone is of length, L E = 11.05, and depth, d E = 0.03.The wave absorption zone and the slowdown zone have lengths L Α = 11 and L D = 2, respectively.The numerical parameters are: Δ 1 = 0.04, N = 128, M = 32, Δ 2 = 0.02 and Δt = 10 -4 .
In Fig. 3, snapshots of the resolved free-surface elevation,  , at several time instants after 20 wave periods, are presented and compared to the experimental measurements of maximum (wave crest) and minimum (wave trough) values of the free surface elevation (Ting and Kirby 1994).The numerical model predicts accurately the breaking depth d b = 0.28, which corresponds to the position x 1 = 40.2,but underestimates the breaking height, as indicated by the deviation of about 9%, of the breaking freesurface elevation, 0.176  .This prediction is still better than the one in Lin and Liu (1998), Bradford (2000) and Christensen (2006).At the outer surf zone (40 < x 1 < 44), the numerical model underestimates the wave height dissipation, as opposed to the inner surf zone (x 1 > 44), where the prediction of the model for the height dissipation is very good.Generally, the numerical results for the minimum free-surface elevation, are in a very good agreement with the trough envelope of the experimental data.In the outer coastal zone, the numerical results for wave shoaling agree adequately with the corresponding experimental data, however, LWS predicts a monotonic wave height increase during shoaling.In the surf zone, the numerical results predict very well the wave setup.In Fig. 4, three typical snapshots of the spanwise vorticity distribution, ω 2 , in the surf zone, are presented at the time that incipient wave breaking occurs, and at time instants t = 0.5T and 0.8Τ after | incipient breaking.Negative vorticity is generated at the breaking wavefront during incipient wave breaking (wave crest at x 1 = 40.2),which corresponds to clockwise recirculation of the fluid.The strength of the surface roller increases, as the spilling breaker propagates towards the inner surf zone, until it reaches its full development (Fig. 4b), and subsequently, as indicated in Fig. 4c, vorticity is advected and diffused in the roller wake.The vorticity distribution close το the sloping bed, due to the bed shear stress, is also shown in Fig. 4, with its amplitude being considerably larger (up to 10 times) than the one developing in the surface roller.According to the LWS methodology, wave breaking and dissipation in the surf zone are coupled with the generation and combined action of the eddy and wave SGS stresses.The distribution of the wave stress 13   , which is the most significant SGS stress in terms of its magnitude, in the surf zone at two time instants (at the time of incipient breaking and 0.5T later), is presented in Fig. 5.It is indicated that the development of the surface roller (see Fig. 4) is connected to the continuous increase of the magnitude of wave stress 13   at the breaking wavefront, which takes place during a time interval Δt = 0.5Τ after the incipient breaking, corresponding to the region 1 ≥ d/d b ≥ 0.78.For time greater than 0.5T after breaking, its strength attenuates gradually before it vanishes in the inner surf zone (x 1 > 45).Similar behavior is also exhibited by the other SGS stresses on the x 1 -x 3 plane, while stresses in the x 2 direction are an order of magnitude smaller.
As mentioned in the introduction, the development of the undertow current, resulting from wave breaking, is constrained by the fact that the period-averaged cross-shore water flux is zero.Fig. 6 presents a typical period-averaged velocity field in the surf zone where it is indicated that the numerical model is able to capture the occurrence of the undertow current.The period-averaged velocity distribution indicates the presence of an onshore directed current, which is due to the mechanisms of the Eulerian drift and the surface roller, in the upper layer of the water depth, and an offshore directed current, i.e., the undertow current, near the bed, which balances the onshore flux.Very close to the bed, a steady current of weak strength, the so-called wave boundary layer streaming, exists offshore to the breaking region and part of the outer surf zone, and is directed towards the shoreline.In Fig. 7, LWS-predicted profiles of the undertow current at four positions in the surf zone are compared to corresponding experimental measurements presented in figure 5 of Ting and Kirby (1994) for the case of spilling breaker.The period-averaged horizontal velocity, U 1 , is normalized with respect to the breaking depth, d b .Overall, the LWS prediction is deemed adequate, since the order of magnitude as well as the gradient of the numerical profiles, agree well with the experimental ones.A significant deviation related to the depth of the minimum of the velocity is observed between numerical and measured profiles, as is indicated in Figs.(c) and (d), which may be attributed to the large difference between our Re d value and the one of the experiments, which is about 7 times larger.
Finally, in Fig. 8, the evolution of the resolved free-surface elevation and the bed shear stress, τ b , distribution at several time instants during a wave period, are presented.The amplitude of the bed stress variation is found to be substantially increased over the sloping bed, especially in the surf zone, becoming up to six times larger (at the region around x 1 = 42.5)than the corresponding amplitude in the flat region (x 1 < 15).The magnitude of τ b decreases in the inner surf zone, following the wave height attenuation, while the decrease of the wavelength with decreasing water depth is indicated by the boldlined snapshot.The position of maximum bed stress does not coincide with the breaking position, where the maximum free surface elevation is located, presenting a phase difference of about 0.5T.

CONCLUSIONS
A numerical model for the simulation of wave propagation and spilling breaking over a constant slope bed is presented.The model is formed by the coupling of the large-wave simulation (LWS) method with a numerical solver of the three-dimensional Navier-Stokes equations.According to the LWS methodology, the wave and eddy subgrid (SGS) stresses are modeled, by use of eddy-viscosity models, and then applied to the viscous flow solver, in order to capture wave breaking and wave energy dissipation in the surf zone.In general, the LWS model predictions related to the characteristics of incipient breaking (breaking depth -height) are very well.The validation of our results is based on the comparison with corresponding experimental measurements conducted by Ting and Kirby (1994).The development of the surface roller in the breaking wavefront is connected to the increase of the strength of the SGS stresses in the outer surf zone and their successive decrease at shallower depths close to the shore.The period-averaged velocity field in the surf zone predicts very well the qualitative characteristics of the undertow current generated along the cross-shore direction by the wave breaking, while the quantitative comparison to the corresponding experimental data (Ting and Kirby 1994) is good.Finally, it is found that the magnitude of the bed shear stress increases substantially in the surf zone becoming up to six times larger than the corresponding one in the inflow flat region.

Figure 1 .
Figure 1.Decomposition of free-surface elevation for the case of a spilling breaker.

Figure 2 .
Figure 2. Sketch of the computational flow domain.

Figure 3 .
Figure 3. Snapshots of the resolved free-surface elevation during shoaling and in the surf zone, over bed of constant slope (tanβ = 1/35).Symbols correspond to the experimental free-surface envelope presented in Ting and Kirby (1994).

Figure 4 .
Figure 4. Vorticity field in the surf zone, at three time instants during wave period, T. Snapshot (a) corresponds to incipient breaking, (b) at time instant t = 0.5T and (c) at t = 0.8T after incipient breaking.Dashed contour lines correspond to negative vorticity, while solid lines to positive one.

Figure 5 .
Figure 5. Snapshots of the wave SGS stress,   13 , that correspond to (a) and (b) of Fig. 4. Note that the first one includes two wave crests (at x 1 ≈ 40 and 44.5), while contour lines are plotted at equal intervals from 0 to 0.001 with a spacing of 0.0002.

Figure 6 .
Figure 6.Period-averaged and spanwise-averaged velocity distribution in the surf zone.Breaking occurs at x 1 ≈ 40.

Figure 7 .
Figure 7. Normalized period-averaged horizontal velocity profiles at four positions in the surf zone, compared to the corresponding experimental data (symbols) presented in figure 5 (c)-(e) of Ting and Kirby (1994).

Figure 8 .
Figure 8. Snapshots of the resolved free-surface elevation and the bed shear stress distribution.The bold lines correspond to the same time instant, while bed slope starts at x 1 = 15.