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

Abstract

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 to model 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.

1. 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.

2. 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.

2.2. 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:(a)two driving forces of water wave radiation stress ∇⋅𝑆 and surface slope ∇𝜂,(b)two resisting forces of bottom friction 𝜏_𝑏 and friction due to turbulent eddies ∇⋅𝜎. The SWE take the form 𝜕𝜂𝜕𝑡+∇⋅𝐻𝑣𝜌=0,𝜕𝑣𝜕𝑡+𝑣⋅∇𝑣=−𝑔𝜌∇𝜂+∇⋅𝑆1𝐻−𝜌𝐻𝜏_𝑏+∇⋅𝜎, (1)where 𝑣 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].

Figure 1 

Boundaries and axes of the surf zone.

2.3. Model of Driving Forces

For a surface gravity wave incident with an angle 𝜃, energy 𝐸, and group velocity-celerity ration 𝑛, 𝑆_𝑥𝑥𝐸+12=𝑛1+cos²𝜃,𝑆_𝑦𝑦𝐸+12=𝑛1+sin²𝜃,𝑆_𝑥𝑦𝐸=𝑛sin𝜃cos𝜃. (2) 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 𝐹_𝑦=1𝐻−𝜕𝑆_𝑥𝑦−𝜕𝑥𝜕𝑆_𝑦𝑦.𝜕𝑦 (3) Similarly, 𝐹_𝑥=1𝐻−𝜕𝑆_𝑦𝑥−𝜕𝑦𝜕𝑆_𝑥𝑥.𝜕𝑥 (4)

2.4. Model of Resisting Forces

2.4.1. The Bottom Stress

Consider the following equation: 𝜏_𝑏=𝑔𝐻𝑟²𝑣𝑣, (5) where 𝑟= Chezy coefficient, given by Handreson [7],𝑟=28.6𝑅_𝑛^1/8,𝑅_𝑛<10⁵,32𝑔^1/2log₁₀𝑅_𝑛8𝑔2.5,𝑅_𝑛>10⁵, (6)𝑅_𝑛=Reynoldsnumber=𝑣𝐻/𝜇 and 𝜇 is the eddy viscosity𝜎=𝜇𝐻∇𝑣. (7)

2.4.2. The Turbulent Stress: The 𝑘-𝜀 Model

To evaluate the turbulent stresses, the depth averaged eddy viscosity 𝜇 is calculated from the transport equation of the depth average 𝑘-𝜀 model, Rastogi and Rodi [8], 𝜕k𝜕𝑡+𝑢𝜕𝑘𝜕𝑥+𝑣𝜕𝑘=𝜕𝜕𝑦𝜇𝜕𝑥𝜎_𝑘𝜕𝑘+𝜕𝜕𝑥𝜇𝜕𝑦𝜎_𝑘𝜕𝑘𝜕𝑦+𝐺+𝑃_𝑘𝑣−𝜀,𝜕𝜀𝜕𝑡+𝑢𝜕𝜀𝜕𝑥+𝑣𝜕𝜀=𝜕𝜕𝑦𝜇𝜕𝑥𝜎_𝜀𝜕𝜀+𝜕𝜕𝑥𝜇𝜕𝑦𝜎_𝜀𝜕𝜀𝜕𝑦+𝐶₁𝜀𝑘𝐺+𝑃_𝜀𝑣−𝐶₂𝜀²𝑘,𝜇=𝐶_𝜇𝑘²𝜀, (8) where 𝑘 is the turbulent kinetic energy, 𝜀 is the time rate of dissipation of 𝑘, and 𝑃_𝑘𝑣 is the production of 𝑘 due to vertical velocity gradient 𝑃_𝑘𝑣=𝐶_𝑘𝑢^∗3𝐻, (9)𝑃_𝜀𝑣 is the production of 𝜀 due to vertical velocity gradient 𝑃_𝜀𝑣=𝐶_𝜀𝑢^∗4𝐻², (10) and 𝐶₁, 𝐶₂, 𝐶_𝜇, 𝜎_𝑘, 𝜎_𝜀, 𝐶_𝑘, and 𝐶_𝜀 are constants given by Rodi [9]; see Table 1.

Table 1 

𝐺 is the turbulent energy production due to mean motion 𝐺=𝜇2𝜕𝑢𝜕𝑥²+2𝜕𝑣𝜕𝑦²+𝜕𝑣+𝜕𝑥𝜕𝑢𝜕𝑦², (11)𝑢^∗is the shear velocity 𝑢^∗=𝑔𝑟²𝑢²+𝑣². (12)

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].

2.5. 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: 𝑣⋅𝑛=0. (13)

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 𝑣𝑦𝜂=Outputofwavemodel,𝑦=Outputofwavemodel. (14) At the two lateral boundaries, the velocity gradient in the normal direction (𝑦-direction) is taken to be zero to ensure their perfect nonreflective character, 𝜕𝑣𝜕𝑦=0,𝜕𝑢𝜕𝑦=0. (15)

2.7. 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].

2.8. Numerical Solution

The discrete form of any variable 𝜙 in the (𝑥,𝑦,𝑡) space will be 𝜙_{𝑖,𝑗,𝑛}, where 𝜙_{𝑖,𝑗,𝑛}=𝜙𝑖Δ𝑥,𝑗Δ𝑦,𝑛Δ𝑡, (16) in which Δ𝑥,Δ𝑦 are the grid spacing in 𝑥  and  𝑦 directions, respectively, 𝑖 and 𝑗 are the grid point labels in 𝑥 and 𝑦 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 𝐻, 𝑘, ε, 𝜈, 𝑆_𝑥𝑥, 𝑆_𝑦𝑦, 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 (𝑖,𝑗), which is the center of the control volume. The velocity component u is defined at the two 𝑥-constant boundaries of the control volume, so that the 𝑥-direction momentum equation is discretized at point (𝑖+1/2,𝑗). The velocity component 𝑣 is defined at the two 𝑦-constant boundaries of the control volume, so that the 𝑦-direction momentum equation is discretized at point (𝑖,𝑗+1/2). The boundaries of the area under study lie at the control-volume boundaries to form a stair-step outer boundary.

Figure 2 

Staggered grid.

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: 𝑢_{𝑖+(1/2),𝑗}=𝐴_1𝑢𝑢_{𝑖+(1/2),𝑗+1}+𝐴_2𝑢𝑢_{𝑖+(1/2),𝑗−1}+𝐴_3𝑢𝑢_{𝑖+1(1/2),𝑗}+𝐴_4𝑢𝑢_{𝑖−(1/2),𝑗}+𝑆_𝑢,𝑣_{𝑖,𝑗+(1/2)}=𝐴_1𝑣𝑣_{𝑖,𝑗+1(1/2)}+𝐴_2𝑣𝑣_{𝑖,𝑗−(1/2)}+𝐴_3𝑣𝑣_{𝑖+1,𝑗+(1/2)}+𝐴_4𝑣𝑣_{𝑖−1,𝑗+(1/2)}+𝑆_𝑣,𝑘_𝑖,𝑗=𝐴_1𝑘𝑘_𝑖,𝑗+1+𝐴_2𝑘𝑘_{𝑖,𝑗−1}+𝐴_3𝑘𝑘_𝑖+1,𝑗+𝐴_4𝑘𝑘_{𝑖−1,𝑗}+𝑆_𝑘,𝜀_𝑖,𝑗=𝐴_1𝜀𝜀_𝑖,𝑗+1+𝐴_2𝜀𝜀_{𝑖,𝑗−1}+𝐴_3𝜀𝜀_𝑖+1,𝑗+𝐴_4𝜀𝜀_{𝑖−1,𝑗}+𝑆_𝜀,𝜂_{𝑖,𝑗,𝑛+1}=𝜂_{𝑖,𝑗,𝑛}+𝑆_𝜂. (17) 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.

3. 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: 𝐻=0.02𝑥. (18) The basin dimensions are 6 m × 6 m. A square computational mesh whose grid size is Δ𝑥=0.1 m, Δ𝑦=0.1m 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.

Figure 3 

Computed longshore current around a jetty obtained by Nishimura et al. [13].

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.

Figure 4 

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

Figure 5 

Computed nearshore currents around a protective jetty.

3.2. 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.

Figure 6 

Bathymetry of the study area depths in meters.

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.

Figure 7 

Velocity field distributions for winter season.

3.3. 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.

Figure 8 

𝑆𝑥𝑥 distribution around Damietta harbor.

Figure 9 

Velocity vectors around Damietta harbor.

3.4. Application 3: The Flow around Detached Break Waters

This case study is for a set of three detached break waters parallel to the shore at a distance of 200 meters from the shore line. The break waters were 200-meter long and 20.0-meter width, with 100 meters spacing between them. A rectangular mesh of 20.0 meters × 5.0 meters is used to simulate the flow in the study area. The radiation stress component 𝑆_𝑥𝑥 distribution is shown in Figure 10, while the flow pattern is shown in Figure 11.

Figure 10 

𝑆𝑥𝑥  distribution around detached break waters.

Figure 11 

Velocity vectors around detached break waters.

4. 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.

Acknowledgment

This work was performed as a part of a research project sponsored by the National Authority for Remote Sensing and Space Sciences (NARSS). The project title is “Assessment of the impact of constructed coastal protection measures on the stability of the Nile Delta shoreline.”