Iteration Coupling Simulation of Random Waves and Wave-Induced Currents

A two-way coupling algorithm for wave-current interaction is developed and implemented into a nearshore circulation model to investigate the effects of fully wave-current interaction on irregular wave transformation over an elliptic shoal. The wave field is simulated by a spectral wave model WABED, and the wave-induced current is solved by a quasi-three-dimensional model WINCM. The surface roller effects are represented in the formulation of surface stress, and the roller characteristics are solved by a roller evolution model. The proposed two-way coupling algorithm can describe both the generation of wave-induced current and the current-inducedwave transformation, which is more physically reasonable than the one-way approaches. The model test with a laboratory experiment shows that wave-induced currents have an important influence on the wave transformation, for example, the wave energy defocusing due to the strong jet-like current along the centerline of the shoal. It is revealed that the accuracy of simulated wave field can be significantly improved by taking into account the two-way wave-current interaction.


Introduction
As waves propagate towards the shore, their characteristics such as heights, lengths, and directions are significantly changed due to shoaling, diffraction, refraction, and breaking.In addition, the wave field would be considerably modified by nearshore currents, which are driven by the gradients of wave excess momentum flux.The numerical simulations of wave, current, and wave-current interaction are important in coastal engineering.
A large number of wave models have been developed to simulate wave propagation in nearshore regions.Wave models are commonly classified into two categories referred to results are discussed with respect to the effects of two-way wave-current interaction on the wave height distribution near the shoal.

Model Description
The multidirectional random wave transformation model WABED is used for wave simulation 9 .The quasi-3D model WINCM is used to simulate wave-induced current 10 .In this section, a brief description of these two models will be given and more details can be found in the original papers.

Wave Model
The governing wave action balance equation is where N is the wave action density, defined as wave energy density divided by the angular frequency σ .x, y are the horizontal coordinates, and θ is the wave direction measured counterclockwise from the x-axis.The first term in the right-hand side of 2.1 represents wave diffraction, where κ is the diffraction intensity coefficient 2.5 8 .C and C g are the wave celerity and group velocity, respectively.ε b is a parameterization of wave breaking function for energy dissipation.C x , C y , and C θ are wave velocity components according to x, y, and θ coordinates, respectively, expressed as where U and V are the current velocity components in x-and y-directions, respectively.k is the wave number.The dispersion relationship between the relative angular frequency σ , the absolute angular frequency ω , the wave number vector k , the current velocity vector U , and the water depth h is

2.3
ε b is the energy dissipation rate due to wave breaking and is parameterized by where ρ is the water density, g is the gravitational acceleration, H rms is the root mean square wave height, σ is the period-averaged relative angular frequency, and D is the total breaking energy dissipation coefficient given by Battjes and Janssen 25 as where α 2 is a one-order constant, f σ/2π is the period-averaged wave frequency, and Q b is the proportionality coefficient of wave heights exceeding the breaking index height H b 0.78 h , estimated according to the Rayleigh distribution

Quasi-3D Wave-Induced Current Model
The governing equations are derived from the Navier-Stokes equations and read as where Λ is a constant with a value around 0.03 to 0.05, u max is the maximum bottom orbital velocity in the x direction, and Γ is a nondimensional parameter with a value around 0.005 to 0.01.

Boundary Conditions
The boundary conditions on the sea bottom are given as where C f is the bottom friction coefficient with a value of order 0.01, and u b and v b are the wave-induced bottom orbital velocities in x-and y-directions, respectively.The free-surface boundary conditions are described by the near-surface shear stress including roller effects 28 , expressed as where τ sx and τ sy are the near-surface shear stresses in x-and y-directions, respectively.T is the wave period, and A is the cross-sectional area of surface rollers.The roller area is solved with the roller energy balance equation proposed by Dally and Brown 29 .The shoreline is treated as a fixed boundary of the computational domain.The seaward boundary is also regarded as fixed by locating it in the deep water where there is no current.Therefore, current velocity components normal to these boundaries are taken as zero.Both the alongshore current velocity and the water level are uniform alongshore at the lateral boundaries.

Numerical Scheme
A forward-marching first-order upwind finite-difference method and the Gauss-Sizel method are used to discrete and solve the wave action balance equation of WABED in a staggered rectangular grid system.The wave velocity quantities are stored at boundaries of the wave action density cell, as shown in Figure 1.Superscripts i, j, k and variables δ x , δ y , δ θ are the grid node number and the spatial steps at x, y, and θ coordinates, respectively.Subscript n represents the nth component of the wave frequency spectrum.Starting from the seaward boundary, wave characteristics are calculated through column to column in the direction of wave propagation.A quadratic upstream interpolation for convective kinematics QUICK is also employed to suppress numerical diffusion in computation.
For the WINCM model, the hybrid method combining the fractional step finite difference method in the horizontal plane with a Galerkin finite element method in the vertical direction is used to discretize the governing equations.As shown in Figure 2, the variables are located on each element according to the staggered method.The computational domain in the horizontal plane is divided into rectangular grid cells.The horizontal velocity components are corrected midway along the cell faces, whereas the vertical components of current velocity and the mean water level are located at the center of each grid cell.The mesh in the horizontal plane is fixed during computation, whereas in the vertical direction, the thickness l k of each layer is automatically updated with time according to the wave setup and setdown.

Two-Way Coupling Algorithm
In the two-way coupling algorithm of wave-current interaction, every explicit time step outer time step consists of a few inner steps of two-way iteration Figure 3 .At first, WABED is used to calculate the wave field without current input.By use of the updated radiation stress, wave height, wave direction, and wave orbital velocity, WINCM is then run to obtain the current field.The depth-integrated mean current velocity and mean water level will be fed back to WABED as the input for the next inner step of iteration.This two-way iteration lasts until a steady state is achieved for both wave and current fields, that is, the maximum relative difference was less than 0.001 between any two consecutive iterations.This two-way coupling algorithm can take into account both the generation of wave-induced current and the current-induced wave transformation, which is indeed more physically reasonable than the one-way methods.With this approach, a more accurate solution of coupled wave and current fields can be obtained and it is thus possible to investigate the effects of fully wave-current interaction.6.1 m in x-direction and a dimension of 7.92 m in y-direction.The wave maker was located at x 0 m.The incident wave was generated using TMA frequency spectrum and the normal distribution function for the directional spectrum.The significant incident wave height is 0.19 m, and the peak period is 1.3 s.Nine wave gauges were arranged along a line of x 12.2 m.The computational domain was 25 m × 30 m with an uniform grid spacing of 0.1 m in both horizontal directions.The water depth was divided into 3 layers along the vertical direction.The TMA-type unidirectional wave spectrum consisting of 18 frequency bins and 36 direction bins was used for the incident wave boundary.The model parameters were set as Λ 0.03, Γ 0.005, and Cf 0.01.

Model Tests
Figure 5 shows the horizontal distributions of calculated wave heights without twoway iteration, after the first and the third two-way iterations.At first, the wave field is calculated by WABED without the current effects.Then, the wave-induced current field starts to be fed back to WABED.After the third two-way iteration, the solutions of wave and current fields tend to converge; thus the effects of two-way wave-current interaction could be distinguished.
As shown in Figure 5 a , when the current effects are not considered, waves are transformed by refraction, diffraction, shoaling, and breaking when propagating over a submerged shoal.Significant wave shoaling can be observed, and the maximum wave height appears above the top of the shoal.The wave energy entering the central part of the shoal is dissipated by wave breaking, and the focusing of wave energy is developed behind the shoal.After waves pass the shoal, wave heights decrease due to the combined effects of breaking, refraction, and diffraction.
However, this pattern of wave height distribution is considerably modified and a significant redistribution of wave energy occurs when wave-induced currents are taken into account, as shown in Figure 5  currents in the direction of wave propagation, and these breaking-induced currents, in turn, change the wave characteristics.Wave heights gradually decrease over the shoal, and wave energy is effectively defocused behind the shoal due to the Doppler shift.This currentinduced refraction results in a narrow zone along the centerline of the shoal where wave heights are relatively smaller than both sides.On the other hand, the strong opposing currents lead to the increasing wave heights near both lateral boundaries.With the proceeding of the two-way coupling, the wave field is further scattered by some local flow circulations but the overall pattern maintains Figure 5 c .This phenomenon is consistent with the finding of Yoon et al.21 and evidently can only be simulated with the two-way coupling algorithm.
Figure 6 shows the horizontal distribution of the calculated current field after the third two-way iteration.The model reasonably simulates the generation of a strong jetlike current along the centerline of the shoal caused by the breaking and energy transfer of waves, which is also observed in the laboratory.On the other hand, due to the impermeable boundaries, two symmetric flow circulations form and expand behind the shoal to maintain the mass conservation over the whole basin.Figure 7 shows the horizontal distribution of the calculated mean water level after the third two-way iteration.Remarkable setdown and setup can be observed above and behind the shoal, respectively, affected by both wave transformation and wave-induced currents.
Figure 8 presents comparison of the measured and the calculated wave heights without two-way iteration, after the first and the third two-way iterations.If the effects of wave-induced currents are excluded, the calculated wave heights have a peak value in the center because of the energy concentration behind the shoal, showing a reversed wave height distribution to the observation.If the two-way wave-current interaction is considered, the calculated wave heights in the center are reduced due to the strong following currents and increased at both sides due to the converging waves refracted by currents.Two local peaks appear laterally from the center.This is in a much better agreement with the observation.It is clear that the two-way wave-current interaction plays an indispensable role in accurately modeling the wave transformation over the shoal.

Concluding Remarks
Based on a spectral wave model WABED and a quasi-three-dimensional wave-induced current model WINCM, the effects of wave-current interaction on wave transformation over an elliptic shoal are investigated with a two-way coupling algorithm.The present two-way coupling algorithm considers the current-induced wave field changes, which is the most important difference from the one-way approaches.Results show that the observed wave height distribution near the shoal, for example, the defocusing of waves due to the strong jet-like current along the centerline of the shoal, can only be explained and reproduced when the two-way wave-current interaction is taken into account.It is revealed that the twoway wave-current interaction is important and should not be ignored in numerical models for accurate simulation of wave transformation with an irregular bathymetry.The present two-way coupling algorithm can significantly improve the model predictability in the wavecurrent coexisting field.

Figure 4 :
Figure 4: Experimental setup of Vincent and Briggs 24 .

Figure 5 :
Figure 5: The horizontal distributions of the calculated wave heights a without the two-way iteration, b after the first two-way iteration, and c after the third two-way iteration.

Figure 6 :
Figure 6: The horizontal distribution of the calculated current fields after the third two-way iteration.

Figure 7 :
Figure 7: The horizontal distribution of the calculated mean water level after the third two-way iteration.
Without the two-way iterationAfter the first two-way iteration After the third two-way iteration

Figure 8 :
Figure 8: Comparison of the observed and the calculated wave heights along line x 12.2 m.
, and W are the mean current velocity, u w , v w , and w w are the periodic velocity, v h and v v are the components of turbulent eddy viscosity in horizontal and vertical directions, and η is the mean wave level.The radiation stress terms are calculated according toZheng 13 .The horizontal and vertical turbulent eddy viscosities are calculated according to Larson and Kraus 26 and Tsuchiya et al. 27 , respectively, expressed as