Iteration Coupling Simulation of Random Waves and Wave-Induced Currents

Abstract

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-induced wave 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.

1. 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 as phase-resolving models and phase-averaged models. The well-known Boussinesq models belong to the category of “phase-resolving models” [1–3]. These models solve the detailed wave characteristics but are more CPU demanding, based on the mass and momentum conservation equations. They are often used to study the problems with small-scale domains. On the other hand, phase-averaged models solve the wave energy balance equation or wave action equation in a frequency domain. The models of this type are widely applied to large study areas from a practical sense. Examples are CMS-Wave [4, 5], SWAN model [6], STWAVE model [7], and WABED model [8, 9]. WABED is a practice-oriented random spectral wave model which is capable of simulating wave shoaling, refraction, diffraction, reflection, depth/current-limited breaking, energy dissipation, and wave-current interaction. It is originally developed by Mase [8] and then improved by numerous researchers with respect to the model capabilities of wave diffraction, breaking, and wave-structure interaction [9–12].

The models for wave-induced currents are mostly based on the concept of radiation stress [13]. Svendsen et al. [14] established a surf zone circulation model SHORECIRC that is quasi-three-dimensional using the long wave theory for the radiation stress formulation in the vertical direction. Warner et al. [15] developed a fully three-dimensional ocean circulation model ROMS which includes a three-dimensional formulation of the vertically varying wave excess momentum flux. A comprehensive comparison between SHORECIRC and ROMS was presented in Haas and Warner [16], who pointed out that the models were very similar for the depth-integrated flows and qualitatively consistent for the vertically varying components. Kumar et al. [17] further improved ROMS for a better model performance in surf zone and rip-current applications. Other relevant studies include Wang and Shen [18], Bennis et al. [19], and Xie [20]. However, while these studies mainly focused on the wave impacts on current (e.g., radiation stress, surface roller, and turbulence mixing), the reverse effects of currents on the deformation of waves are not considered in enough detail, thus only partially representing the wave-current interaction. As highlighted by Yoon et al. [21] and Choi et al. [22], the effects of wave-induced current are very important and should not be neglected for accurate simulation of wave transformation over an irregular bathymetry. Yoon et al. [21] investigated the effects of breaking-induced currents on refraction-diffraction of irregular waves over a submerged shoal, using an iterative feedback scheme between wave and current models. They found that the one-way approach is insufficient to obtain the correct distribution of wave height behind the shoal. However, the shallow water equation was used which cannot deal with the vertical variations of radiation stress and current velocity. Choi et al. [22] carried out a numerical study of wave transformation due to wave-induced currents, based on a spectral wave model and a quasi-three-dimensional circulation model. Although they successfully simulated the breaking wave-induced strong jet-like current and the current-induced wave defocusing along the centerline of a shoal, the surface roller effects were not taken into account in their model. Zheng and Tang [10] and Zhang et al. [23] have shown that the inclusion of roller contribution leads to a significantly improved description of the long-shore and cross-shore wave-induced currents.

In this study, the wave model WABED and the quasi-three-dimensional wave-induced current model WINCM are coupled and used to investigate the effects of two-way wave-current interaction on irregular wave transformation over an elliptic shoal. The surface roller effects are represented in the formulation of surface stress, and the roller characteristics are solved by a roller evolution model. A two-way coupling algorithm is developed to deal with both the generation of wave-induced current and the current-induced wave transformation. A model test is conducted on the basis of the experiment of Vincent and Briggs [24], and results are discussed with respect to the effects of two-way wave-current interaction on the wave height distribution near the shoal.

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

2.1. Wave Model

The governing wave action balance equation is 𝜕𝐶_𝑥𝑁+𝜕𝐶𝜕𝑥^𝑦𝑁+𝜕𝐶𝜕𝑦^𝜃𝑁=𝜅𝜕𝜃2𝜎𝐶𝐶_𝑔cos²𝜃𝑁_𝑦𝑦−12𝐶𝐶_𝑔cos²𝜃𝑁_𝑦𝑦−𝜀_𝑏𝑁, (2.1) where 𝑁 is the wave action density, defined as wave energy density divided by the angular frequency (𝜎). (𝑥, 𝑦) are the horizontal coordinates, and 𝜃 is the wave direction measured counterclockwise from the 𝑥-axis. The first term in the right-hand side of (2.1) represents wave diffraction, where 𝜅 is the diffraction intensity coefficient (= 2.5) [8]. 𝐶 and 𝐶_𝑔 are the wave celerity and group velocity, respectively. 𝜀_𝑏 is a parameterization of wave breaking function for energy dissipation. 𝐶_𝑥, 𝐶_𝑦, and 𝐶_𝜃 are wave velocity components according to 𝑥, 𝑦, and 𝜃 coordinates, respectively, expressed as 𝐶_𝑥=𝐶_𝑔𝐶cos𝜃+𝑈,_𝑦=𝐶_𝑔𝐶sin𝜃+𝑉,_𝜃=𝜎sinh2𝑘ℎsin𝜃𝜕ℎ𝜕𝑥−cos𝜕ℎ𝜕𝑦+cos𝜃sin𝜃𝜕𝑈𝜕𝑥−cos²𝜃𝜕𝑈𝜕𝑦+sin²𝜃𝜕𝑉𝜕𝑥−sin𝜃cos𝜃𝜕𝑉,𝜕𝑦 (2.2) where 𝑈 and 𝑉 are the current velocity components in 𝑥- and 𝑦-directions, respectively. 𝑘 is the wave number. The dispersion relationship between the relative angular frequency (𝜎), the absolute angular frequency (𝜔), the wave number vector (𝑘), the current velocity vector (𝑈), and the water depth ℎ is 𝜎²=𝑔^→𝑘tanh^→𝑘ℎ,𝜎=𝜔−^→𝑘⋅𝑈^→. (2.3)

𝜀_𝑏 is the energy dissipation rate due to wave breaking and is parameterized by 𝜀_𝑏=⟨𝐷⟩𝜌𝑔𝐻²_rms/8𝜎, (2.4) where 𝜌 is the water density, 𝑔 is the gravitational acceleration, 𝐻_rms is the root mean square wave height, 𝜎 is the period-averaged relative angular frequency, and ⟨𝐷⟩ is the total breaking energy dissipation coefficient given by Battjes and Janssen [25] as ⟨𝐷⟩=𝛼₂𝜌𝑔4𝑄_𝑏𝑓𝐻²_𝑏, (2.5) where 𝛼₂ is a one-order constant, 𝑓=(𝜎/2𝜋) is the period-averaged wave frequency, and 𝑄_𝑏 is the proportionality coefficient of wave heights exceeding the breaking index height (𝐻_𝑏=0.78h), estimated according to the Rayleigh distribution 1−𝑄_𝑏ln𝑄_𝑏𝐻=−^rms𝐻_𝑏². (2.6)

2.2. Quasi-3D Wave-Induced Current Model

The governing equations are derived from the Navier-Stokes equations and read as 𝜕𝑈+𝜕𝑥𝜕𝑉+𝜕𝑦𝜕𝑊𝜕𝑧=0,𝜕𝑈𝜕𝑡+𝑈𝜕𝑈𝜕𝑥+𝑉𝜕𝑈𝜕𝑦+𝑊𝜕𝑈𝜕𝜕𝑧=−𝑔𝜂−𝜕𝜕𝑥𝑢²_𝑤−𝑤²_𝑤−𝜕𝜕𝑥𝑢_𝑤𝑣_𝑤−𝜕𝜕𝑦𝑢_𝑤𝑤_𝑤+𝜕𝜕𝑧𝜈𝜕𝑥^ℎ𝜕𝑈+𝜕𝜕𝑥𝜈𝜕𝑦^ℎ𝜕𝑈+𝜕𝜕𝑦𝜈𝜕𝑧^𝑣𝜕𝑈,𝜕𝑧𝜕𝑉𝜕𝑡+𝑈𝜕𝑉𝜕𝑥+𝑉𝜕𝑉𝜕𝑦+𝑊𝜕𝑉𝜕𝜕𝑧=−𝑔𝜂−𝜕𝜕𝑦𝑣²_𝑤−𝑤²_𝑤−𝜕𝜕𝑦𝑢_𝑤𝑣_𝑤−𝜕𝜕𝑥𝑣_𝑤𝑤_𝑤+𝜕𝜕𝑧𝜈𝜕𝑥^ℎ𝜕𝑈+𝜕𝜕𝑥𝜈𝜕𝑦^ℎ𝜕𝑈+𝜕𝜕𝑦𝜈𝜕𝑧^𝑣𝜕𝑈,𝜕𝑧 (2.7) where 𝑡 is time, 𝑈, 𝑉, and 𝑊 are the mean current velocity, 𝑢_𝑤, 𝑣_𝑤, and 𝑤_𝑤 are the periodic velocity, 𝑣_ℎ and 𝑣_𝑣 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 to Zheng [13].

The horizontal and vertical turbulent eddy viscosities are calculated according to Larson and Kraus [26] and Tsuchiya et al. [27], respectively, expressed as 𝑣_ℎ=Λ𝑢_max𝑣𝐻,_𝑣=Γ𝐶𝐻, (2.8) where Λ is a constant with a value around 0.03 to 0.05, 𝑢_max is the maximum bottom orbital velocity in the 𝑥 direction, and Γ is a nondimensional parameter with a value around 0.005 to 0.01.

2.3. Boundary Conditions

The boundary conditions on the sea bottom are given as 𝑣_𝑣𝜕𝑈𝜕𝑧=𝐶_𝑓(𝑈+𝑢_𝑏)(𝑈+𝑢_𝑏)²+(𝑉+𝑣_𝑏)²,𝑣_𝑣𝜕𝑉𝜕𝑧=𝐶_𝑓(𝑉+𝑣_𝑏)(𝑈+𝑢_𝑏)²+(𝑉+𝑣_𝑏)², (2.9) where 𝐶_𝑓 is the bottom friction coefficient with a value of order 0.01, and 𝑢_𝑏 and 𝑣_𝑏 are the wave-induced bottom orbital velocities in 𝑥- and 𝑦-directions, respectively.

The free-surface boundary conditions are described by the near-surface shear stress including roller effects [28], expressed as 𝜏_𝑠𝑥𝑔=−𝜌8𝑑𝐻²+𝑑()𝑑𝑥𝐴𝐶/𝑇𝑑𝑥cos²𝜏𝜃,_𝑠𝑦𝑔=−𝜌8𝑑𝐻²+𝑑()𝑑𝑥𝐴𝐶/𝑇𝑑𝑥cos𝜃sin𝜃, (2.10) where 𝜏_𝑠𝑥 and 𝜏_𝑠𝑦 are the near-surface shear stresses in 𝑥- and 𝑦-directions, respectively. 𝑇 is the wave period, and 𝐴 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.

2.4. 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 (𝑖, 𝑗, 𝑘) and variables (𝛿_𝑥, 𝛿_𝑦, 𝛿_𝜃) are the grid node number and the spatial steps at 𝑥, 𝑦, and 𝜃 coordinates, respectively. Subscript 𝑛 represents the 𝑛th 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.

Figure 1 

Arrangement of variables in WABED.

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 (𝑙_𝑘) of each layer is automatically updated with time according to the wave setup and setdown.

(a)

(b)

Figure 2 

Arrangement of variables in WINCM.

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

Figure 3 

Model procedures.

3. Model Tests

Vincent and Briggs [24] performed an experiment (case B5) of irregular wave propagating over an elliptic shoal in Coastal and Hydraulics Laboratory's wave basin of the US Army Engineer Research and Development Center. The experimental setup is shown in Figure 4. The water depth at the flat bottom was 0.4572 m. The center of the elliptic shoal was located at 𝑥=6.10 m and 𝑦=13.72 m with a water depth of 0.1524 m. The shoal has a dimension of 6.1 m in 𝑥-direction and a dimension of 7.92 m in 𝑦-direction. The wave maker was located at 𝑥=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 𝑥=12.2 m.

Figure 4 

Experimental setup of Vincent and Briggs [24].

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 𝐶𝑓=0.01.

Figure 5 shows the horizontal distributions of calculated wave heights without two-way 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.

(a)  

(b)

(c)

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.

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(b). The wave breaking over the shoal induces strong jet-like 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 current-induced 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 jet-like 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 6 

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

Figure 7 

The horizontal distribution of the calculated mean water level after the third two-way iteration.

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.

Figure 8 

Comparison of the observed and the calculated wave heights along line 𝑥=12.2 m.

4. 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 two-way 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 wave-current coexisting field.

Acknowledgments

The work was supported by the National Natural Science Foundation of China (no. 50979033) and the Special Research Funding of State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering (no. 2009585812).