Numerical simulation of the three-dimensional wave-induced currents on unstructured grid

By coupling the three-dimensional hydrodynamic model with the wave model, numerical simulations of the threedimensional wave-induced current are carried out in this study. The wave model is based on the numerical solution of the modified wave action equation and eikonal equation, which can describe the wave refraction and diffraction. The hydrodynamic model is driven by the wave-induced radiation stresses and affected by the wave turbulence. The numerical implementation of the module has used the finite-volume schemes on unstructured grid, which provides great flexibility for modeling the waves and currents in the complex actual nearshore, and ensures the conservation of energy propagation. The applicability of the proposed model is evaluated in calculating the cases of wave set-up, longshore currents, undertow on a sloping beach, rip currents and meandering longshore currents on a tri-cuspate beach. The results indicate that it is necessary to introduce the depth-dependent radiation stresses into the numerical simulation of wave-induced currents, and comparisons show that the present model makes better prediction on the wave procedure as well as both horizontal and vertical structures in the wave-induced current field.


Introduction
In coastal nearshore zone, the formation of fluid flow caused by the attenuation of the wave height is known as the wave-induced current, i.e., the longshore current which is developed while the breaking wave is oblique to the shoreline, the rip current occurred under the irregular bathymetry, and the cross-shore current known as undertow on the slope. All of them develop in the surf zone, and have their own morphology as the longshore current is roughly parallel to the shoreline (Longuet-Higgins, 1970); the rip current is approximately shore-normal, which originates within the surf zone and broadens outside the breaking region (MacMahan et al., 2006); the undertow is oriented shoreward near the surface, and directed seaward close to the bottom, while the flow pattern is opposite outside the surf zone (Nobuoka et al., 2005). In view of this, the waveinduced currents have significant impact on the coastal sediment transport and pollutants movement.
In order to study the wave-induced currents, many works have been accomplished by several researchers in terms of theoretical (e.g. Longuet-Higgins and Stewart, 1964;Svendsen and Lorenz, 1989), physical (e.g. Okayasu et al., 1988;Ting and Kirby, 1994) and numerical model (e.g. Dong and Anastasiou, 1991;Li et al., 2007). Because of its economic and adaptability, the numerical model has been fully developed and applied to study current in the surf zone, which includes two numerical solutions. Firstly, based on the extended Boussinesq equations (BES), the wave propagation and breaking are incorporated in the numerical model, and the wave-driven currents are calculated by timeaveraged water velocity over the wave period. Wei et al. (1995) and Chen et al. (2000) constructed the famous numerical module of Fully Nonlinear Wave model (FUN-WAVE) to simulate the wave and current in the surf zone. Secondly, according to Longuet-Higgins and Stewart (1964), the radiation stresses are defined as the excess momentum flux, which represents the momentum transferring from the wave to current. Therefore, the wave-induced currents are simulated by coupling the hydrodynamic model with the wave model through the radiation stresses, and the wave is also affected by the current through the velocity and direction. Some numerical solutions of the wave action have been developed and taken to simulate the wave-induced current. For instance, based on different types of mild-slope equations (MSE, Berkhoff, 1972), Yoon et al. (2004), Johnson et al. (2005), Tang et al. (2008), Tong et al. (2010) and Xie (2011) simulated the wave-driven currents on structured grid by the finite difference method; Wu and Liu (1985), Roland et al. (2009) and Tang and Wei (2010) established the unstructured wave-induced current models respectively. But the BES and MSE models require the prohibitively high spatial resolution over the entire computational region, which will lower the large-scale current computational efficiency. Castelle et al. (2006), Warner et al. (2008) and Choi et al. (2009) simulated the wave-induced current by combining different current modules with the Simulating Waves Nearshore model (SWAN), which can simulate the wave propagation in large-scale area, but the coupled processes are only one-way and the wave-current interactions are not considered, and the three-dimensional structure of the wave-induced current due to three-dimensional radiation stresses has not been discussed in the above research.
It seems clear that the development of the wave-induced current model which can simulate the three-dimensional structure of the wave-induced current in the complicated coastal geometry region will be extremely useful. In the paper, the eikonal equation derived from the extended mildslope equation is coupled with the wave action equation to simulate the diffraction effect (Hong, 1996). Furthermore, a three-dimensional hydrostatic primitive-equation model (Chen et al., 2003) is coupled with the wave model bi-directionally to simulate the three-dimensional wave-driven currents. The coupled model is solved numerically over the unstructured triangular grids, which provides great flexibility for modeling the currents in the complex geomorphology of the nearshore zone. Several numerical cases with the wavedriven currents over the surf zone are performed to validate the proposed wave-induced current model.

Large-scale wave model
To simulate the wave propagation from deep to shallow water, the action balance equation is adopted and reads as: N = E/σ where ; C θ and C σ are the propagation velocities of N in θ and σ spaces, respectively; C g are the velocities of wave group; S tot is the source and sink of N, which can represent the process of the wave generation, dissipation, nonlinear wave-wave interaction, etc. Energy dissipation mode in the wave breaking is derived from Zheng et al. (2008), the other source and sink of energy in this paper refers to the solution in the SWAN model.
The eikonal equation can be obtained from the mild-slope equation to simulate the wave diffraction, which reads as: where k is the wave number and satisfies the dispersion relation ; is the wave number vector affected by the wave amplitude and current velocity; U is the flow velocity; and the parameter . For the sake of considering the wave diffraction effect in the action balance equation, the wave number k is replaced by the wave number vector K, and C, C g , C θ , C σ are modified as: (3) Based on the numerical solution of Eqs.
(1)-(4), a largescale wave model with the refraction and diffraction effect is established. The effect of simulating wave diffraction over the large-scale region has been validated in literature (Wang and Zhang, 2015) by use of the above model, which will not be discussed in this paper.

Modified 3-D hydrodynamic model
The 3-D hydrodynamic model derived from FVCOM has been modified to simulate the current fields in the surf zone, including the wave-induced turbulent effect and the wave-averaged forcing due to the wave breaking, and written as: is the vertical sigma coordinate; i and j represent the x and y directions; w is the velocity component normal to the sigma surface; D is the water depth; ρ 0 is the water density in barotropic mode; ρ is the water density in baroclinic mode; K m and A m are the vertical turbulence and horizontal diffusion coefficient generated by the current and wave; and F i is the radiation stress component.
Many expressions (e.g. Zheng et al., 2000;Xia et al., 2004;Zhang, 2004;Mellor, 2005) of the depth-dependent radiation stresses have been developed; the equations derived by Zhang (2004) are adopted: where E is the wave energy density; δ ij is the Kronecker symbol. The turbulence and diffusion coefficient generated by the wave-current interaction are the sum of results which are caused by current and wave separately, and have been solved by Xia et al. (2004) as: The wave-induced horizontal turbulence coefficient is calculated by derived by Larson and Kraus (1991), is the maximum velocity at the bottom of wave; and the current-induced horizontal turbulence is defined as Eq. (10), is the area of the control volume.
The wave-induced vertical eddy viscosity is defined as Eq. (11), and the current-induced vertical eddy viscosity is calculated through the numerical solution of the turbulent kinetic energy and macro-scale equations, which has been described by Mellor and Yamada (1982). . (11) The calculation of the wave and current model are coupled bi-directionally in each step, the radiation stresses and the wave-induced turbulence coefficients are provided to the 3D current model, and the current field and surface elevation are provided to the wave module. Additionally, the boundary conditions of the current module will be modified to account for the wave-current interaction. If only consider the factor of wave, the wave-induced current field will be obtained by the coupled model.

Discretization procedure and solution method
The 3-D hydrodynamic model (FVCOM) is solved by the finite volume method in the integral form of Eqs. (5) and (6) over non-overlapping, unstructured triangular grids, and the discrete formats of the control equations of FV-COM have been constructed by Chen et al. (2003). To facilitate the coupled procedure of the wave and current model, the vector variables: such as k, K, u, v, c and C g are placed at the centroid, and the scalar variables: η, D, A, E, a and C are placed at the node, as shown in Fig. 1.
The control equations are numerically solved by the Euler upwind format in time, and using Green formula converts the surface integration to the line integration in space, written as: where A is the area of the control body, S is the edge of the control body and m is the total edge number, is the value of in the edge e, n e and is the normal vector and length of the edge.
According to different factors, the spectral action balance Eq. (3) has been split into four equations for numerical simulation, as follows: The discrete way of action density in frequency is solved by the FCT (Flux-Corrected Transport) method (Boris and Book, 1997), and in the wave direction is calculated through the second-order accurate implicit Crank-Nicolson scheme, and in space is solved numerically by the unstructured grid finite volume approach in the control volume TCE (Tracer Control Element, Fig. 1). The influence of source and sink is numerically dispersed by the second-order and semi-implicit difference scheme. The eikonal equation is also solved by using the unstructured-grid finitevolume approach of the cell center format in the control volume MCE (Momentum Control Element, Fig. 1). The concrete discrete formats of the above equations are given in literature (Wang and Zhang, 2014).
The radiation stress items are solved by the unstructured grid finite volume approach of the cell center format and using Green formula converts the surface integration to the line integration in the control body MCE, given as:

Model validation
To verify the validity of the model on the numerical simulation of the wave-driven current fields, the coupled model is firstly used to simulate the longshore currents and tested against the Snell's Law; then the model is applied to simulate the undertow on a sloping beach to demonstrate the ability of the model to produce three-dimensional structure of the wave-induced current; last the rip current and the meandering longshore current under irregular bathymetry are simulated to produce the horizontal velocity distribution of the wave-induced current.
3.1 Longshore current Zou et al. (2003) measured the wave-induced long-shore currents with different slopes in laboratory. The parameters for two experimental cases are listed in Table 1. The computation meshes are plotted in Fig. 2a.
In Table 1, h 0 is the still water depth before the slope, H 0 is the incident wave height, T is the wave period, γ is the parameter of wave breaking, and τ is the bottom friction coefficient under the wave and current action. The wave propagates along the x direction, and the numerical calculations will stop until the current field reach an approximately stable state. The numerical wave direction distributions of Case 1 are shown in Fig. 2b, and the comparisons of the wave height and direction between the numerical and analytical results are plotted in Figs. 3a and 3b, respectively.
The numerical results of wave-induced longshore current distributions for two cases are shown in Figs. 4a and 4b, and the comparisons of the wave height, water elevation and current velocity between the computational and experimental results in two cases are shown in Figs. 5 and 6. From the above results, it is obvious that the coupled module can accurately simulate the wave, wave set-up and longshore current due to the oblique incident wave on the slope beach. The deviation exists between the numerical and experimental velocity results, which may be caused by the surface roller of wave breaking which is not considered in the numerical model, the surface roller (Zheng and Tang, 2009) will lead to the momentum be delivered to the shore, and  cause the position of the peak of the wave-induced current offset to the surf zone. Then, the surface roller will lead to the influence of the longshore current on the wave direction in the numerical result which is wider than that of the experimental result. The effect of the surface roller will be added into the model in future research.

Undertow on the slope
In the surf zone, the undertow will occur when the wave propagates onshore on the slope, the flow of the surface water toward the shoreline is driven by the breaking wave, and the mass flux is balanced by the seaward flow at the bottom (Veeramony and Svendsen, 2000;Tajima and Madsen, 2006), and this phenomenon has been observed in the experiment (Ting and Kirby, 1994).
In the experiment, the wave was generated at a constant depth of 0.4 m and propagates on a linear beach with 1:35 slope. The wave height in the depth of 0.4 m was 0.127 m and the period was 2 s. The computational domain was discretized using a grid size of 0.05 m in the horizontal and 30 levels in the vertical, including 28750 cells and 17700 nodes; the model time step was 0.1 s, the bottom friction coefficient of the model was 0.025, the parameter of the wave breaking was 0.76, and the wave-induced current field was stable one hour after the start of calculation. The wave was breaking at 6.4 m away from the shore in the numerical result, which is the same as the experimental result.
The vertical and horizontal profiles of the radiation stresses along the depth and wave direction are plotted in Figs. 7-8, and the vertical distribution of the simulated flow field is shown in Fig. 9. The results show that the three-dimensional nearshore current is caused by the inequality of the vertical and horizontal distributions of radiation stresses. The simulated water set-up and set-down are contrasted to the measured results in Fig. 10. The water set-down derived by the wave breaking is increased from zero meter along the wave propagation direction outside the surf zone, reaches the maximum value at the wave breaking point, and then the surface elevation will begin to rise and change to the water set-up in the surf zone. The velocity variations of numerical   WANG Ping et al. China Ocean Eng., 2017, Vol. 31, No. 5, P. 539-548 and experimental results at eight locations along the depth are shown in Fig. 11. Based on the numerical results, there exist two opposite circulations in the near-shore zone. The water flows shoreward near the surface and seaward near the bottom in the surf zone, and the flow pattern is opposite outside the surf zone.
In order to compare the differences between the 2D and 3D radiation stresses, based on the 2D radiation stresses the wave-induced current in the same slope is simulated. The numerical results at one hour after the start of calculation are plotted in Fig. 12. There is no obvious vertical difference in the wave-induced current field, and the wave-induced current is a reciprocating flow, which is caused by the limit of the computational domain. The undertow on the slope is unstable, thus the vertical structure of the wave-induced currents can only be described by the numerical model with three-dimensional radiation stresses.
Differences of the measured and simulated results in the flow velocity also exist, which mainly due to the radiation stress in the present work are based on the linear wave theory. However, the wave in shallow water is nonlinear, which will affect the vertical distribution of the radiation stress. And the wave breaking mode and surface roller are also the influencing factors of the numerical result.
3.3 Rip current Borthwick and Foote (2002) carried out laboratory studies of the structure of the rip current over a tri-cuspate beach, which was situated on a 1:20 plane beach, and had the planar dimension of 23 m alongshore and 15 m offshore. At the beach, the still water depth was given by , y L =12 m, R= 4 m, and there are three cusps located from and . Fig. 13 illustrates the still water depth contours of the tri-cuspate beach.
In the experiment, two wave conditions were considered: Case A, regular waves with the period T=1.2 s, offshore height H 0 =0.125 m and incident wave angle , which are perpendicular to the shoreline; Case B, regular waves with the period T=1.2 s, offshore height H 0 =0.125 m and incident wave angle , which are oblique to the shoreline. Fig. 14 shows the comparisons of the wave height between the modeled and the observed in different representative sections, which indicates that the model can describe the wave propagation and breaking pro-    cesses accurately.
The vertical velocities at four different points are compared with the observed values, which are illustrated in Fig.  15. Herein the velocity components U and V represent the longshore and onshore directions, respectively. Inside the surf zone, water flows shoreward near the surface and seaward near the bottom, and along the embayment, at the location of x=-3 m and x=-2.5 m, the offshore component has a peak value near the bed boundary layer, which is due to the combined effect of the rip current and undertow. There are few differences between the simulated and observed, which are caused by the wave nonlinearity and reflections from the fixed cuspate beach. The unstable nature of the rip current also made the measurement data extremely scattered. However, the comparisons do prove that the simulated velocity profile captures the major distribution of the experimental currents.
The comparisons of the horizontal current distribution between the simulated and measured in Cases A and B are illustrated in Figs. 16 and 17. In the cases, the regular waves   which are perpendicular to the tri-cuspate beach caused a system of primary nearshore circulation, including the shoreward-directed currents on the cusp horns, and the longshore currents from the cusp ridge to the embayment trough, and the seaward rip currents at the trough. The oblique incident waves produce the steady meandering longshore current at the fixed cuspate beach. Comparisons show that the distribution of the simulated flow field mainly agrees with the observed in the experiment.

Conclusions
Combing the wave model with hydrodynamic model, a three-dimensional model is developed to simulate the waveinduced current on unstructured grid. The cases of the wave set-up, longshore current, undertow on a sloping beach, rip current and meandering longshore current on tri-cuspate beach are evaluated to validate the numerical model. The major characteristics of the model and some conclusions include the following.

546
WANG Ping et al. China Ocean Eng., 2017, Vol. 31, No. 5, P. 539-548 The wave model with the diffraction effect is based on the modified wave action equation and eikonal equation. The hydrodynamics model is coupled with the wave model by considering the impact of the wave radiation stresses, the wave-induced turbulence and the wave-current interaction.
The wave-induced current model is implemented on unstructured triangular grid, which is more suitable to be applied in the estuarine and coastal region characterized with complicated coastline.
Several three-dimensional wave-induced current cases are utilized to validate the established coupled model. The cases cover the variety of wave-induced current phenomena including the wave set-up, longshore currents, undertow on a sloping beach, rip currents and meandering longshore currents due to different incident directions. The numerical results show that the three-dimensional structures of wave-induced currents are mainly due to the horizontal and vertical distribution of the radiation stresses, and the coupled model can effectively perform their major features.
According the above horizontal and vertical flow distributions, the three-dimensional wave-induced currents are the complicated and important hydrodynamic factors in the actual nearshore, and have a significant impact on the coastal sediment transport, coastal deformation and pollutant transport. The coupled module has the potential of modeling and predicting the influence of the wave-induced current on the above phenomenon, but more improvement should be introduced to numerical model in future due to the existence of deviations between the measured and simulated results, and the study of coupling other procedures into the model should be also implemented in future.