A Numerical Model of the Wave-Induced Currents in the Turbulent Coastal Zone

A numerical model is developed, validated and applied to the turbulent coastal currents. The currents are driven by the sea surface slope and the radiation stresses of water waves. They are resisted by friction due to turbulent eddies and sea bottom. The k-ε model is used tomodel the turbulent stresses. Five simultaneous nonlinear partial differential equations govern the depth-averaged dynamics in the surf zone. An implicit finite-difference scheme is used to obtain an accurate numerical solution of the resulting initial-boundary value problem. It is tested against the case of straight coast with uniform bottom slope and a protective jetty. To investigate the actual wave-induced currents, the model is applied to simulate the currents for three real case studies. Results show that the model could be used to compute currents caused by the constructing coastal protection measures and could predict the locations of accretion and scouring.


Introduction
As surface waves progress onshore, through the nearshore zone, they transform due to the combined effect of ambient currents and bottom variations.As they approach the surfzone-between the shoreline and the breaker line-the energy of the broken waves is transformed into radiation stresses, which drive turbulent currents.The effects of these currents on exciting bottom sediments and transporting them offshore and longshore were studied by Martinez and Harbaugh [1].
Fahmy [2] showed that if the horizontal dimensions of the flow domain are much larger than the vertical dimension, the depth-averaged and wave period-averaged mass and momentum equations can be used to represent the flow.As a result of large horizontal scales in the surf zone area, the flow can be considered turbulent even for very small velocity values.There are two approaches of mathematical models of wave-induced current: the first is based on the radiation stress theory, while the second is based on the Boussinesq equation.Recently, many researchers dealt with this topic.Significant publications include the one of Mengguo and Chongren [3].They used a simplified model for the resisting force caused by eddy viscosity.Kim and An [4] presented a numerical model that takes the wave-current interaction into consideration and quasi-3D equations were used in circulation analyses.

Model Theoretical Formulation
2.1.Model Assumptions.The depth-average shallow water model equations (SWE) are invoked for coastal regions where hydrostatic approximation is valid [5].These equations are based on the following assumptions.(i) The flow is incompressible.(ii) The vertical acceleration and Coriolis acceleration are neglected.We further assume that the two driving forces of wind stress and horizontal density gradient are neglected.

Model Governing Equations.
Consider the turbulent flow in the surf zone between the shoreline and the breaker line (Figure 1).The currents are affected by four dynamical forces.These are categorized as follows: (b) two resisting forces of bottom friction   and friction due to turbulent eddies ∇ ⋅ .The SWE take the form where V is the depth-average speed,  is the turbulent eddies frictional stress tensor,  is the water depth,  is the water desist,  is the water surface elevation, ∇ = (/, /),  is the time,  is acceleration of gravity, ,  are the offshore and longshore coordinates, respectively, and  is the radiation stress tensor defined by Dean and Dalrymble [6].

Model of Driving Forces.
For a surface gravity wave incident with an angle , energy E, and group velocity-celerity ration , The gradients of the radiation stress components provide driving forces in the nearshore zone.For waves propagating obliquely into the surf zone, breaking will result in a reduction in wave energy and an associated decrease in   , which is manifested as a wave thrust .The components of  are Similarly, where  = Chezy coefficient, given by Handreson [7], where  is the turbulent kinetic energy,  is the time rate of dissipation of , and  V is the production of  due to vertical velocity gradient V is the production of  due to vertical velocity gradient and  1 ,  2 ,   ,   ,   ,   , and   are constants given by Rodi [9]; see Table 1.
is the turbulent energy production due to mean motion * is the shear velocity These are five simultaneous nonlinear partial differential equations ( 1)- (8).They are continuity equation and momentum transport equations, with the turbulence model equations being the governing equations that are used to calculate the velocity vector in the surf zone.The two driving forces that appear in the momentum equations are the gradients of the wave's radiation stress and the mean slope of sea level.To feed the model with the radiation stress components, the authors developed a wave model to study wave propagation from deep water to the shore line; see Fassieh et al. [10].

Boundaries of the Study Area.
The surf zone is the flow region that follows the breaking zone where water waves break and transform their energy into turbulent energy whose momentum flux derives coastal currents.The surf zone is bounded by the shore line and three open boundaries, an offshore boundary and two lateral boundaries (see Figure 1).These open boundaries are artificial boundaries through which water flows and are needed for reasons of computational economy.The off shore boundary is placed before waves breaking by about one wave length.The lateral boundaries are perpendicular to the shore line (too far from the study area to affect calculations).
2.6.Boundary Conditions.At the shore line boundary, the normal velocity component vanishes: At the offshore boundary, Dirichlet boundary conditions are imposed via the wave model.The water level is taken to be equal to the local mean sea level.The steady (depth and period-averaged) velocity components due to mass flux defined from the wave model are At the two lateral boundaries, the velocity gradient in the normal direction (y-direction) is taken to be zero to ensure their perfect nonreflective character, (15)

Initial Conditions.
The water surface in the study area is assumed initially to be horizontal (zero surface displacement with no setup or setdown).Since the bottom topography is known, the local water depth is given initially at each point.The initial velocity distribution is chosen such that the mass conservation equation is satisfied everywhere, that is, at each grid point; see Fahmy [11].

Numerical Solution.
The discrete form of any variable  in the (, , ) space will be  ,, , where in which Δ, Δ are the grid spacing in  and  directions, respectively,  and j are the grid point labels in x and y directions respectively, Δ is the time step, and  is the time level.
2.9.The Computational Grid.The transport equations are discretized over a staggered uniform rectangular mesh; see Anderson et al. [12].The staggered mesh is superimposed on the area under study as shown in Figure 2. In the staggered grid technique, the domain is divided into control volumes.The scalar variables H, k, , ],   ,   , and   are defined at the center of the cell, so that the transport equations of the scalar variables  and  and the continuity equation are discretized at the point (i, j), which is the center of the control volume.The velocity component u is defined at the two xconstant boundaries of the control volume, so that the xdirection momentum equation is discretized at point ( + 1/2, ).The velocity component v is defined at the two yconstant boundaries of the control volume, so that the ydirection momentum equation is discretized at point (,  + 1/2).The boundaries of the area under study lie at the controlvolume boundaries to form a stair-step outer boundary.An implicit scheme is used to discretize the five governing equations ( 1)- (8).Centered space difference is used for spatial derivatives, and backward difference is used for time derivatives.Upwind difference is used to discretize the convective terms.The various terms of the governing equations are discretized at the time level ( + 1).The produced system of algebraic equations has the following form: The  coefficients contain viscosity, convective velocities from both directions, and water depth.The  terms are the source terms for each variable.Starting with the initial velocity and water depth, the linear algebraic equations are solved by iteration to get values at the next time step using Gauss Seidel method.The calculated values from the previous time step are taken to be initial guesses for the time next step.

Model Verification and Applications
Given the basin lateral dimensions and bottom topography, a Cartesian mesh is superimposed to define the computational domain.The water depth is assigned at all the grid nodes of the mesh covering the domain.The wave model [10] simulates the manner by which the incident wave characteristics in deep water are transformed through wave refraction, shoaling, diffraction, and breaking as the wave propagates over the water surface of a basin-varying bottom topography.The wave model is used to compute wave height and wave orthogonal direction at all the grid nodes.These results are used to compute the components of the radiation stress tensor and the particle velocity components at each grid point.This surf zone model uses the radiation stress components as source terms that represent wave forces, for the depth-averaged period-averaged momentum equations.The numerical solution of the coupled system of the five nonlinear partial differential equations of momentum, mass, turbulent energy, and eddy viscosity produces the surf zone velocity field.
3.1.Nearshore Currents around a Protective Jetty.Consider a square basin of one dimensional uniform slope depth, whose formula is as follows: The basin dimensions are 6 m × 6 m.A square computational mesh whose grid size is Δ = 0.1 m, Δ = 0.1 m is imposed on the basin.This grid size gives total of (60 × 60) grid points at which the bottom depth array, (  ,   ), is defined.A 1.5 m long jetty is situated perpendicular to the shoreline, as shown in Figure 3. Regular surface gravity waves are obliquely    incident at 30 ∘ to the shoreline, with height 1.8 cm and wave period 1.2 seconds.Upon breaking, part of their energy is dissipated and the reformed part is transformed into an induced nearshore current that moves around the jetty.An infinitely long straight shoreline is imposed for the case study of the boundary condition utilizing a pure longshore current.
The same case study was studied by Nishimura et al. [13] and Watanabe and Maruyama [14].They also built a physical model and took measurements for the longshore currents around the jetty.Figure 3 shows the computed longshore current obtained by Nishimura et al. [13].The observed current velocity vectors (measured in the physical model) are plotted in Figure 4.The sets of parallel arrows in the figure denote longshore components of local velocities measured by using small current meters.Figure 5 shows the current field computed by this model.The velocity vectors show that two circulation zones are formed: one zone upstream of the jetty (the illuminated region) and another, much larger, circulation zone downstream the jetty (the shadow region).A fair agreement is found between the predictions of this model and the physical model measurements.However, the present model results are more clear and accurate in simulating the current dynamics at the sharp edges of the jetty and the narrow corners near the coast.

Application 1: Current System near the Nile Delta Coast.
The coastal area extending from Gamasa to Ras EL-Bar (about 50 km long) is covered by a mesh of grid size 500 m × 500 m.The basin bathymetry is entered through digitization of the available contour map for the Egyptian Nile Delta coast.Linear interpolation is used to determine water depth at the grid points lying between contour lines.The bathymetry of the study area is shown in Figure 6.
The random surface gravity waves propagating toward the Nile Delta coast have an average period of  = 8 seconds and significant wave height that ranges from 1 m. to 3 m.coming from northwest direction during summer and winter, while coming from north-northwest during fall season.
Figure 7 shows the velocity fields in the study area during winter season.The highest magnitude velocity vectors are observed around the head of Ras El-Bar.

Application 2:
The Flow around Damietta Harbor.The flow around Damietta harbor is studied using a grid of 50.0 m. × 50.0 m.The deep water wave height is 1 meter and coming from the northwest direction.The xx-component of the radiation stress tensor is shown in Figure 8.The velocity vectors around the harbor are shown in Figure 9.The velocity vectors are shown to bend around the harbor walls with decreasing velocity in the front of the main wall due to blocking which may cause accretion at the shore line near the main wall and increasing velocity vectors around the minor wall which may cause a scouring at the shore near the minor wall.

Application 3: The Flow around Detached Break Waters.
This case study is for a set of three detached break waters

Conclusions
In the present study, a nonlinear model of the wave-induced unsteady currents in the turbulent coastal zone was presented.The main conclusions obtained from this study are as follows.
(1) The present model could be applied to accurately compute wave-induced unsteady currents within and outside of the surf zone.(2) The present model could be used to compute currents caused by the constructing coastal protection measures like detached break waters parallel to the shore and protective jetties.
(3) The highest magnitude velocities are observed around the head of Ras El-Bar which may cause a scouring at this location.
(4) The wave-induced current around Damietta harbor shows a decreasing velocity in the front of the main wall which may cause accretion at the shoreline near the main wall, and increasing velocity around the minor wall which may cause scouring at the shore near the minor wall.

Figure 1 :
Figure 1: Boundaries and axes of the surf zone.

Figure 4 :
Figure 4: Observed longshore current around a jetty in a physical model, Nishimura et al. [13].

Figure 6 :
Figure 6: Bathymetry of the study area depths in meters.

Figure 7 :
Figure 7: Velocity field distributions for winter season.