Numerical investigation of the three-dimensional velocity fields induced by wave-structure interaction

Submerged shore-parallel breakwaters for coastal defence are a good compromise between the need to mitigate the effects of waves on the coast and the ambition to ensure the preservation of the landscape and water quality. In this work we simulate, in a fully three-dimensional form, the hydrodynamic effects induced by submerged breakwaters on incident wave trains with different wave height. The proposed three-dimensional non-hydrostatic finite-volume model is based on an integral form of the Navier-Stokes equations in σ-coordinates and is able to simulate the shocks in the numerical solution related to the wave breaking. The obtained numerical results show that the hydrodynamic phenomena produced by wavestructure interaction have features of three-dimensionality (undertow), that are locally important, and emphasize the need to use a non-hydrostatic fully-three-dimensional approach.


Introduction
Submerged breakwaters for coastal defence are a good compromise between the need to mitigate the effects of waves on the coast and the ambition to ensure the preservation of the landscape and water quality. The main function of this structure is to protect the shoreward area of the breakwater from wave actions by way of attenuating the incoming waves. Submerged breakwaters lower the wave energy in the protected area by provoking the breaking of the incoming waves over the structure.
However, if not properly designed, such structures can force circulation patterns that enhance shoreline erosion rather than shoreline accretion.
According to the simple response-function model proposed by Ranasinghe et al. [1], erosion at the shoreline is expected to occur when the resultant current field contains divergent alongshore currents in the entire protected zone and return off-shore again at the gaps transporting sediments out of the protected area (circulation of an erosive nature). Accretion at the shoreline is expected when convergent alongshore when an additional converging flow closer to the shoreline, promoting sediment deposition, is generated (circulation of an accretive nature).
Coastal currents and, more generally, hydrodynamic phenomena produced by wave-structure interaction have features of three-dimensionality that are locally important [2]. The most important of the above threedimensional phenomena and the cause of offshore sediment transport is the undertow [3] which consists of a circulation in the vertical plane in which the near-bed current velocities are off-shore directed in the surf zone [4].
In literature, the current circulations are generally simulated by two-dimensional Boussinesq models [5,6,7,8] obtained by depth-averaging a simplified form of the three-dimensional Navier-Stokes equations. This approach, based on the depth-averaged motion equations, assumes a simplified distribution of the hydrodynamic quantities along the vertical direction and proves to be valid only in the cases in which a fully three-dimensional representation of the fluid flow is not needed.
In this work, we propose a three-dimensional numerical model in which the non-hydrostatic Navier-Stokes equations are expressed in integral form on a coordinate system in which the vertical coordinate is varying in time [9,10]. The boundary conditions for pressure are placed on the upper face of each computational cell. The solution is advanced in time by using a three-stage Strong Stability Preserving Runge-Kutta (SSPRK) fractional step numerical method, and at each stage a pressure correction formulation is applied in order to get a fluid velocity field which is divergencefree.
A shock-capturing technique based on high-order WENO reconstructions is employed in order to discretize the fluid motion equations. At every cell interface, the numerical flux is computed by solving an approximate HLL Riemann problem.
The proposed model is used to simulate the circulation patterns induced by normally incident waves on a beach with submerged longshore bars and rip channels [11,12].
The results obtained show that features of threedimensionality in the fluid flow induced by wave-structure interaction are locally important and entails the need of a fully three-dimensional model.

Governing integral three-dimensional -coordinate equations
The integral form of the momentum equations over a control volume ∆V(t) which varies in time is given by where ∆ ( ) is the surface of the control volume, ( = 1,3) is the outward unit normal vector to the surface ∆ ( ), u l ( l =1, 3) and v m ( m =1,3) are respectively the fluid velocity and the velocity of the surface of the control volume, both defined in a Cartesian reference system of coordinates x l ( l =1,3), ρ is the density of the fluid, T lm is the stress tensor and f l ( l =1,3) represents the external body forces per unit mass vector in which δ 13 is the Kroneker symbol and p is the total pressure defined by the sum of the hydrostatic and the dynamic component where G is the constant of gravity, q is the dynamic pressure, η is the free surface elevation, the comma with an index in subscript denotes the derivative as [ ] ,l = ∂[ ]/ ∂x l and (x 1 , x 2 , x 3 , t) is a Cartesian coordinate system. The first integral on the right-hand side of equation (1) can be rewritten as By applying Green's theorem, the right-hand side of equation (4) becomes By introducing equation (5) in equation (1) In order to simulate the fully dispersive wave processes, equation (6) can be transformed in the following way.
Let ( 1 , 2 , ) = ℎ( 1 , 2 , ) + ( 1 , 2 , ) where ℎ is the depth of still water. Let ( 1 , 2 , 3 , ) be a system of curvilinear coordinates which varies in time so as to follow the time variation of the free-surface elevation, the transformation from the Cartesian coordinates ( 1 , 2 , 3 , ) to the curvilinear coordinates is The following relation is valid This coordinate transformation basically maps the varying vertical coordinates in the physical domain to a uniform transformed space where 3 spans from 0 to 1.
By following the procedure proposed by [8] we define the transformation matrix = ⁄ and its inverse ̅ = ⁄ ( , =1,3) . The metric tensor and its inverse are defined by = and = ̅ ̅ , respectively. The Jacobian of the transformation is defined by � = det ( ). It is not difficult to verify that in the above-mentioned transformation � = .
In the curvilinear coordinate system the expression of the momentum equation reds where ���� is the average of over the volume element ∆ * = ∆ 1 ∆ 2 ∆ 3 ; ∆ * + and ∆ * − indicate the contour surfaces of the volume element * on which is constant and which are located at the larger and at the smaller value of respectively. Here the indexes , and γ are cyclic.
In order to ensure conservation of mass over the water column, we define a time-varying control volume, ΔV � = ΔA xy * H where ΔA xy * = ∆ξ 1 ∆ξ 2 , and derive the following integral form of the continuity equation in which + and − indicate the contour lines of the surface element ∆ * on which is constant and which are located at the larger and at the smaller value of respectively. Equation (10) represents the governing equation for the surface movements.
Equations (9) and (10) represent the expressions of the three-dimensional motion equations as a function of the ���� and � variables in the time dependent coordinate system ( 1 , 2 , 3 , ). The numerical integration of the mentioned equations (9) and (10) allows the fully dispersive wave processes simulation.
The turbulent kinematic viscosity in the stress tensor is estimated by a Smagorinsky sub grid model.

Numerical scheme
A combined finite-volume and finite-difference scheme with a Godunov-type method has been applied to discretize equations (9) and (10). By following the strategy described by [2][3][4][5][6][7][8][9][10] a staggered grid framework is introduced, in which the velocities are placed at the cell centres and the pressure is defined at the horizontal cell faces. The state of the system is known at the centre of the calculation cells and it is defined by the cellaveraged values ���� and � . ( ) is the time level of the known variables while ( +1) is the time level of the unknown variables. The solution procedure uses a threestage third-order nonlinear strong stability-preserving (SSP) Runge-Kutta scheme for equations (9) and (10).
A pressure correction formulation is applied to obtain a divergence free velocity field at each time level. By this method, an auxiliary velocity field, obtained by numerically integrating equation (9) devoid of the dynamic pressure term, is corrected by introducing a scalar potential Ѱ , which is calculated by the wellknown Poisson pressure equation Equation (11) is defined at the horizontal cell centre and it is discretized by a second order cell centred finitedifference scheme. By this way, equation (11) can be reduced to an algebraic linear system like AΨ = b , where A is the coefficient matrix (with 15 non-zero diagonal coefficient), Ψ is the scalar potential vector and is the vector of constant terms. This algebraic linear system is solved by using an implicit scheme based on a four-colour Zebra line Gauss-Seidel alternate method [13] and a multigrid V-cycle accelerator as described in [14].
The updating of the flow variables ���� and � is based on the following sequence.
(1) High order WENO reconstructions, from cell averaged values, of the point values of the unknown variables at the centre of the contour faces which define the calculation cells. At the centre of the contour face which is common with two adjacent cells, two-point values of the unknown variables are reconstructed by means of two WENO reconstruction defined on two adjacent cells [15,16].
(2) Advancing in time of the point values of the unknown variables at the centre of the contour faces by solving, by an HLL Riemann solver [17], a local Riemann problem with initial data given by the pair of point-values computed by two WENO reconstructions defined on the two adjacent cells.
(3) Calculation of the spatial integrals on the righthand side of equations (9) and (10).
(4) Solution of the Poisson pressure equation by using a four-color Zebra line Gauss-Seidel alternate method and a multigrid V-cycle.
(5) Correction of the hydrostatic velocity field by using a scalar potential Ψ.
(6) Advancing in time of the total local depth (equation (10) by using the non-hydrostatic velocity field.

Three-dimensional nearshore currents induced by submerged breakwaters
The model is used to simulate the circulation patterns induced by normally incident waves on a beach with a longshore bar and rip channels. The experimental tests performed by [11] and the one performed by [4] on a fixed barred beach with periodically spaced rip channels are here reproduced by the presented three-dimensional non-hydrostatic model.
In the validations, we focus on the circulation patterns enhanced by the presence of submerged breakwaters -that produce the erosion or accretion of the shoreline -and the undertow.

Nearshore current circulation patterns
A deep insight about the different circulation patterns that can arise in presence of submerged breakwaters and normally incident waves can be found in [1] and in [12]. They presented a relationship linking environmental and geometrical properties of the system to the shoreline mode of response, i.e. accretive or erosive, thus identifying two different types of structure-induced circulation current.
Because of the presence of the submerged breakwaters, the incoming waves break at different abscissa along the shoreline: the incoming waves directly approaching the shore through the gaps steepen and finally break due to water depth limitations; whilst the reduction of wave energy due to the wave breaking over the structure causes the transmitted waves to break closer to the shoreline than those at the gaps and with a smaller wave set-up. These set-up variations govern the flow patterns in the protected area behind the structures.
Mass conservation requires that the water flowing onshore over the barrier returns off-shore again trow through the gaps. The resulting diverging current circulation system in the lee of the breakwater is composed of two symmetric circulation cells and drive sediments out of the protected area causing shoreline erosion (Fig. 1). Nevertheless, depending on the direction of the alongshore gradient in the mean water level close to the shoreline, the alongshore flow direction may reverse and be directed towards the centreline of the barrier leading to a converging current circulation system that cause shoreline accretion (Fig. 2).  A plan view and a cross section of the wave basin used to perform the test is illustrated in Figure 3. The basin is characterized by a length equal to 17.2m and width equal to 18.2m. The beach has an initial steep slope (1:5) followed by a milder slope (1:30). There are three submerged breakwaters: the two lateral ones are 3.66m long, and the central one is 7.32m. The distance between the submerged breakwaters is 1.82m. The seaward edges of the bar sections are located at approximately =11.1m with the bar crest at =12m and their shoreward edges at =12.3m.
Thanks to the equations [18] it was possible to obtain the dimensions of the hump considered as composed of two branches of a parabola.
In order to study the influence of rip channels and submerged breakwaters on the nearshore dynamics, the test cases B and D, with incident and monochromatic waves, have been selected ( Table 1). The computational grid resolution is ∆x =0.025m, ∆y =0.05m and the time step is 0.025s.  In Figure 4(a) and 4(b) the resulting plane view of the time-averaged velocity fields are shown respectively for test B and test D.
From a comparison with the above-mentioned circulation patterns (Figs. 1-2), the numerical model results clearly show an accretive mode of response for test B and an erosive mode of response for test D, in agreement with [12]. Figure 5 shows the comparison between the resulting time-averaged velocity field for test B and the experimental results obtained by [19]. The numerical results are in good agreement with the experimental data both in the gap between the bars and in the region area between the submerged breakwaters and the coastline  (red arrows) obtained by [19] A more detailed current comparison has been analysed for two different longshore sections: the first placed in correspondence of the onshore side of the bar, at x =12.2m (Figs. 6a and 6b); the second placed between the bar and the shoreline, at x =13.0m (Figs. 6c and 6d). Figure 6 shows the comparison between the crossshore (U 1 ) and longshore (U 2 ) time-averaged velocity components obtained by the proposed numerical model and the laboratory data relative to test B in [11] and [20], at the two different longshore sections. A very good agreement observed at the two sections for both the cross-shore and longshore currents shows the ability of the numerical model to simulate the current variations. The cross-shore current is offshore directed in the whole gap and reaches its maximum value in the central section. As can be noticed from Figure 7, the proposed numerical model is able to simulate the current variations in fairly good agreement with the experimental measurements. Figure 7 shows a small underestimation of the simulated rip current in comparison with the experimental measurements. As underlined in [4], the rip current tends to become a surface current as it flows offshore and, consequently, in this case the small differences between numerical and experimental results could be likely due to the choice of the measurement points in the water column.  [11] in the channel at x =11.5m (a) , x =11.8m (b) and x =12m (c)

Three-dimensional structure of the circulation currents
As shown in this section, for this test case, the proposed model is able to well reproduce the three-dimensional structure of the circulation currents by adopting only eight layers along the water depth.
The hydrodynamic phenomena produced by the wave-structure interaction have features of threedimensionality that are locally important [4].
The most important of the above three-dimensional phenomena -which is the cause of offshore sediment transport -is the undertow [3], which consists of a circulation in the vertical plane in which the near-bed current velocities are off-shore directed in the surf zone.
Here, the vertical velocity profiles of test B obtained by means of the presented three-dimensional model are validated against the laboratory measurement of [4]. Figure 8 shows the comparison between the vertical distribution of time-averaged cross-shore velocity obtained by the three-dimensional numerical model and the experimental measurements presented in [4] at a location 2m offshore of the bar (x = 9.0m, y =13.6m) and a location inside the channel (x =11.75m, y =13.6m).
The results are normalized by the celerity c = √Gh where G is the constant of gravity and ℎ is the depth of still water. From Figure 8, it can be seen that the vertical velocity exhibit strong depth variations, twisting over depth, with the surface velocity going mainly offshore and the bottom current going shoreward.
In particular, at the offshore location (x = 9.0m, y =13.6m) the largest (offshore directed) current is located at the free surface, whereas the weakest (onshore directed) current is near the bottom.
Inside the channel (x = 11.75m, y =13.6m) the cross-shore current shows higher magnitude and a maximum value at mid-depth.
The good agreement between numerical results and experimental measurements shows the ability of the model to simulate the variations of the velocity field along the vertical direction. the experimental measurements (circles) obtained by [4] at points: x = 9.0m , y = 13.6m (a) and x = 11.75m, y = 13.6m (b)

Conclusions
The proposed model is based on an integral form of the Navier-Stoked equations in a time-dependent coordinate system. The equations of motion are discretized by means of a Shock Capturing finite-volume numerical procedure based on high-order WENO reconstructions. The solution procedure for the equations of motion uses an accurate third-order Runge-Kutta fractional step method and applies a pressure correction in order to obtain a divergence free velocity field. At each cell interface an HLL Riemann solver is used to numerically solve a local Riemann problem. Fully three-dimensional numerical results obtained by the proposed threedimensional model has been compared with experimental measurements by [11,4]. Such comparison demonstrates the capacity of the proposed threedimensional model to simulate the strong variations of the velocity fields along the water depth that occur near the submerged breakwaters. The good agreement with the experimental results shows that the proposed threedimensional model can be a useful tool for the simulation of nearshore hydrodynamics fields and the design of submerged burrier defence systems.