3D Numerical Simulation of the Interaction between Waves and a T-Head Groin Structure

: The aim of coastal structures for the defense from erosion is to modify the hydrodynamic ﬁelds that would naturally occur with the wave motion, to produce zones of sedimentation of solid material, and to combat the recession of the coastline. T-head groin-shaped structures are among the most adopted in coastal engineering. The assessment of the e ﬀ ectiveness of such structures requires hydrodynamic study of the interaction between wave motion and the structure. Hydrodynamic phenomena induced by the interaction between wave motion and T-head groin structures have three-dimensionality features. The aim of the paper is to propose a new three-dimensional numerical model for the simulation of the hydrodynamic ﬁelds induced by the interaction between wave ﬁelds and coastal structures. The proposed model is designed to represent complex morphologies as well as coastal structures inside the domain. The numerical scheme solves the three-dimensional Navier–Stokes equations in a contravariant formulation, on a time-dependent coordinate system, in which the vertical coordinate varies over time to follow the free-surface elevation. The main innovative element of the paper consists in the proposal of a new numerical scheme that makes it possible to simulate ﬂows around structures with sharp-cornered geometries. The proposed numerical model is validated against a well-known experimental test-case consisting in a wave train approaching a beach (non-parallel with the wave front), with the presence of a T-head groin structure. A detailed comparison between numerical and experimental results is shown.


Introduction
Coastal defense structures are aimed at modifying the free-surface elevation and velocity fields produced during wave events, in a way to damp the wave energy and to prevent the nearshore wave-induced currents transporting the suspended solid material offshore. T-head groin-shaped coastal structures are one of the most used types of structures that are able to induce the sedimentation of solid material near the shore, by physically blocking the long-shore currents that transport the sediment (see [1,2] for details of these structures). Therefore, the study of the hydrodynamics around T-head groin-shaped structures is essential to properly design them, with the purpose of obtaining the desired response in terms of hydrodynamic fields and the morphological modifications of the coastline.
As recalled by [3][4][5], coastal currents and hydrodynamic phenomena related to the interaction of waves and coastlines have a three-dimensional structure that can be locally important. One of the most important aforementioned three-dimensional phenomena is the undertow current [6], i.e., a circulation in the vertical plane consisting in offshore-directed velocities near the bed and onshore-directed velocities near the free surface. The importance of these near-bed current velocities relies on their substantial impact on the morphology modification. Furthermore, the presence of coastal structures can induce vortices and velocity fields with a significant vertical component, which can cause local phenomena of erosion and scour around the structures.
The study of the interaction between waves and coastal structures can be carried out by means of numerical models. Models that are based on the nonlinear shallow water equations (NSWEs) [7,8] or Boussinesq equations [9][10][11] represent the velocity fields exclusively in terms of horizontal components that can be depth-averaged (NSWE-type models) or defined at a given depth (Boussinesq-type models). In particular, in Boussinesq-type models, the vertical distribution of the horizontal velocity field is taken into account by means of a Taylor series expansion of the velocity. In order to allow NSWE-type or Boussinesq-type models to adequately represent undertow currents, correction procedures of the horizontal velocity fields are necessary [12]; in any case, the vertical components of the velocity are not represented and the corrections of the horizontal velocities are effective for wave fronts parallel to the coastline and simple coastal morphologies. Furthermore, with NSWE-type or Boussinesq-type models, it is not possible to represent the fully three-dimensional structures of the hydrodynamic fields induced by the interaction between the wave and structure.
In order to overcome the limitations imposed by models based on depth-averaged equations, a number of recent models have been devised to solve the three-dimensional Navier-Stokes (NS) equations [13][14][15]. In the solution of NS equations in the context of free-surface flows, the tracking of the free-surface is a challenging issue. A number of models [16][17][18][19] solved this issue by adopting the so-called volume of fluid (VOF) technique [20]. Several numerical models available in the literature employ the VOF technique, as it is a robust and stable one; OpenFOAM® [18] and IHFOAM® [21] (a model derived from OpenFOAM®specifically designed for coastal and ocean hydrodynamics) are among the most popular numerical models in which the VOF technique is implemented. With this technique, a correct assignment of pressure and kinematic boundary conditions at the free surface is difficult, because the vertical fluxes can arbitrarily cross the computational cell. Some recent models are based on a particular coordinate transformation, called sigma coordinate transformation [22][23][24], by means of which the Cartesian vertical coordinate is transformed into a moving vertical coordinate in order to follow the free surface. Hence, at every instant, the free surface is always located at the upper boundary of the computational grid and pressure boundary conditions can be correctly assigned at it. Cannata et al. [25] solved an integral form of NS equations with a shock-capturing scheme over a boundary-conforming moving grid; this allows the models to properly treat the free-surface boundary, as well as to precisely take into account lateral boundaries in complex geometries.
All the above-described numerical models fall under the general class of models in which the free-surface elevation is represented by means of a function of horizontal coordinates and time, so that at a given position in the horizontal plane and in time, a single value of the free-surface elevation is allowed. This approach does not make it possible to directly simulate a wave front that curls over, forming a tube. This limitation is currently overcome by smoothed particle hydrodynamics (SPH) meshless-type models [26,27], in which the position of the water particles that belong to the free surface is tracked one by one, by means of a Lagrangian approach. The main drawback of SPH models is that they are computationally very expensive, so their use is limited to cases in which the spatial scales are lower than the typical scales of a coastal engineering problem.
Three-dimensional models with shock-capturing schemes that employ a vertical coordinate transformation to follow the free surface allow the wave front to be directly simulated up to the point at which the wave front inclination overcomes the vertical direction inclination. After that, the successive wave breaking and particles' turbulent mixing phase is represented with a moving vertical front. This vertical front is treated as a discontinuity in the water depth and in the product between it and the flow velocity and no criterion is needed to set the breaking location and features [24,25]. With these shock-capturing schemes, non-plunging breakers can be directly simulated while plunging breakers are represented by means of a moving discontinuity in the solution.
Cannata et al. [28] recently proposed a three-dimensional non-hydrostatic model that solves an integral form of NS equations in contravariant formulation, on a time-dependent curvilinear coordinate system. In models by [28,29], the time dependence of the curvilinear coordinate system is related to the free surface movement, following the above-mentioned approach. The model proposed by [28], thanks to the use of boundary conforming grids, is suitable for complex geometries that are present in coastal regions, but it is not indicated to represent domains in which structures with a sharp-cornered geometry are present. In fact, in the model, the possibility to introduce internal boundary conditions is not considered. One of the most used strategies in order to represent internal boundary conditions in computational domains, alternative to the use of boundary conforming grids, is the immersed boundary (IB) method [30,31], by means of which a structure can be treated as a virtual body inside the computational domain. In particular, Roman et al. [32] proposed an implementation of the IB method into the context of curvilinear grids. With the use of a boundary conforming grid, the IB method is not necessary to impose internal boundaries, since the computational grid can be chosen so that the coordinate lines follow the boundaries of the structures.
In this work, we propose a novel numerical scheme for the solution of the integral contravariant form of the NS equations presented in [28]. In the proposed numerical scheme, a two-steps predictor-corrector method is employed: In the predictor step, the governing equations are solved with the hydrostatic pressure assumption; in the corrector step, a correction of the hydrodynamic field is made, by means of the solution of a Poisson-like equation. Differently from the numerical scheme presented by Cannata et al. [28], the proposed numerical scheme allows the representation of three-dimensional velocity fields around sharp-cornered coastal structures. In particular, we implement a proper strategy for the imposition of boundary conditions in the numerical solution of the Poisson-like equation. By doing so, the proposed three-dimensional non-hydrostatic numerical model is able to reproduce internal objects inside the computational domain. Furthermore, we adopt a novel procedure to assign the dynamic pressure boundary conditions at the free-surface. The proposed model is used to numerically simulate the wave field and the wave-induced three-dimensional velocity fields produced by the interaction between random wave trains and a T-head groin-shaped sharp-cornered structure. The numerical results are compared against experimental measurements by [33].
The paper is structured as follows: In Section 2, we describe the model equations; in Section 3, we present the numerical model; in Section 4, the numerical results obtained by the proposed model are compared with the experimental measurements by [33]; and in Section 5, we present the conclusion of the study.

Governing Equations
The main feature of the proposed three-dimensional numerical model relies on the fact that the physical domain occupied by the free-surface flow is described by moving curvilinear coordinates defined by a time-dependent coordinate transformation. The physical domain is divided into grid cells bounded by curvilinear coordinate surfaces that move with a velocity different from the fluid velocity; the free-surface coincides with a moving coordinate surface. By means of a time-dependent coordinate transformation, the irregular and time-varying physical domain is transformed into a regular and fixed computational domain, in which the upper boundary represents the free surface. In the proposed integral form of the governing equations, the calculation cells, which have a time-varying geometry, are used as control volumes. Consequently, control volumes, whose boundary surfaces move with a velocity that is different from the fluid velocity, must be defined in the time-dependent curvilinear coordinate system.
Let us consider a coordinate transformation x i = x i ξ 1 , ξ 2 , ξ 3 , τ , from a system of Cartesian coordinates → x = x 1 , x 2 , x 3 , t to a system of moving curvilinear coordinates ξ 1 , ξ 2 , ξ 3 , τ , with vectors. The metric tensor is g lm = → g (l) · → g (m) , with (l, m = 1, 2, 3). The Jacobian of the transformation is given by: The velocity vector of the moving curvilinear coordinates, expressed in the Cartesian coordinate system, is: where ∂ → x /∂t represents the time derivative of the position vector that represents a point that is defined by fixed curvilinear coordinates. As demonstrated by [28], the contravariant component of the velocity vector → v is given by: In the curvilinear coordinate system, the integral contravariant expression of the continuity equation for a control volume that coincides with a material fluid volume ∆V(τ) with a surface ∆A(τ), at time t, which moves with a velocity → v different from the fluid velocity → u , is: where ρ is the density of the fluid and n m is the covariant outward normal vector component.
Let → λ be a constant and parallel vector field. This vector field is represented in the Cartesian coordinate system by constant and uniform components and in the curvilinear coordinate system by constant and space-varying (not uniform) covariant components, λ l . Indeed, since the base vectors of the curvilinear coordinate system vary from point to point, the values of vector components λ l , which are relative to those base vectors, must vary from point to point to represent the same physical direction at every point. From a general point of view, to express the momentum conservation law in an integral form, the rate of change of the momentum of a material volume and the total net force must be projected in a physical direction. As stated by [28], let us define a parallel vector field → λ, which is normal to the coordinate line on which the ξ l coordinate is constant at point P 0 ∈ ∆V: For the sake of brevity, we indicate The integral form of the momentum conservation equation projected in the physical direction defined by the vector field λ k , reads: where T km are the contravariant components of the stress tensor and f k is the external body force per unit mass vector. The stress tensor can be expressed in the form T km = −pg km + 2µS km , where p is the pressure, µ is the dynamic viscosity, and S km are the contravariant components of the strain rate tensor (2µS km is the deviatoric part of the viscous stress tensor). Let us consider ∆V(τ) as a volume element defined by surface elements bounded by curves lying on the coordinate lines. We define the volume element in the physical space as: and the volume element in the transformed space as: It is possible to see that the volume element defined in Equation (7) is time dependent while the one defined in Equation (8) is not. Similarly, we define the surface element in the physical space as ∆A(τ), which is time dependent, and the surface element in the transformed space as ∆A * , which is constant over time. By doing so, Equation (6) can be rewritten as [28]: where ∆A * α+ and ∆A * α− indicate the contour surfaces of the volume ∆V * on which ξ α is constant and which are located at the larger and the smaller value of ξ α , respectively. Here, the indices α, β, and γ are cyclic. By adopting the same control volume, ∆V(τ), the integral contravariant continuity of Equation (4) can be rewritten as: Let x 3 = 0 at the horizontal reference plane defined by the still free surface. Let the water depth be H x 1 , where h x 1 , x 2 is the still water depth and η x 1 , x 2 , t is the free-surface elevation. Let the bottom elevation with respect to the horizontal reference plane be Let us indicate by G the acceleration due to gravity and let us split the pressure p into a hydrostatic part, ρG η − x 3 , and a dynamic one, q. In order to represent the bottom and surface geometry and to correctly assign the pressure and kinematic conditions at the bottom and at the free surface, we introduce the following transformation from Cartesian to curvilinear coordinates, in which coordinates vary in time in order to follow the free-surface movements [28]: Under the transformation (Equation (11)), the components of vector → v are: This coordinate transformation basically maps the time-varying coordinates of the physical domain into a fixed coordinate system ξ 1 , ξ 2 , ξ 3 , where ξ 3 spans from 0 to 1. Let where → k indicates the vertical unit vector and × indicates the vector product. As stated by [28], the Jacobian of the transformation can be written in the form √ g = H √ g 0 . This makes it possible to write the integral contravariant momentum equation expressed by Equation (9) in a three-dimensional conservative form, in which the conserved variables are given by the cell-averaged product between the water depth H and the three contravariant components of the velocity u l with l = 1, 3. To this end, let us define the following cell-averaged values in the transformed space: where ∆A * ξ 1 ξ 2 = ∆ξ 1 ∆ξ 2 is the horizontal surface element in the transformed space. By recalling that ∆V * is not time dependent, by dividing by ∆V * , and by recalling Equation (13), Equation (9), expressed in the time-dependent coordinate system defined in Equation (11), becomes: In order to take into account turbulence, the deviatoric part of the viscous stress tensor is replaced by the turbulent stress tensor, R kα = 2(ν + ν t )S kα , in which ν t is the eddy viscosity and it is estimated by a Smagorinsky sub-grid model.
Let us integrate the continuity Equation (10) over a vertical water column, between the bottom and the free surface, which is bounded by coordinate surfaces. In the transformed space defined by Equation (11), recalling Equation (14), the integral of Equation (10) over the vertical water column reads: in which ∆ξ α+ and ∆ξ α− indicate the contour lines of the surface element ∆A * on which ξ α is constant and which are located at the larger and the smaller value of ξ α , respectively. Equation (16) represents the governing equation that predicts the free-surface motion. Equations (15) and (16) represent the expression of the three-dimensional motion equations as a function of the Hu l and H variables in the coordinate system ξ 1 , ξ 2 , ξ 3 . The numerical integration of the abovementioned equations allows the fully dispersive wave-propagation simulation.

The Numerical Scheme
Equations (15) and (16) are solved by means of a numerical scheme based on the one proposed by [28]. In particular, in order to simulate hydrodynamic fields in the presence of internal objects, in the model presented in this paper, a new numerical procedure for the solution of the Poisson-like equation is proposed. The numerical scheme proposed by [28] consists in a combined finite-volume and finite-difference scheme with a Godunov-type method, and it is briefly described as follows.
A two-stage second-order Runge-Kutta method [34] is used in the advancing of the solution.
Starting from the values of Hu l In Equation (17), the superscripts (0RK), (1RK), and (2RK) indicate values obtained at the end of every Runge-Kutta step, while superscripts (·) (1) and (·) (2) indicate the intermediate value obtained by the predictor-corrector procedure to be described. At each stage of the Runge-Kutta method, the predictor fields are calculated as follows: where D indicates the discretization of the right-hand side of Equation (15), in which the term related to the dynamic pressure gradient is neglected and L indicates the discretization of the right-hand side of Equation (16). D H (s) , u (s) and L H (s) , u (s) are calculated by solving local Riemann problems at cell interfaces. Neglect of the dynamic pressure gradient term leads to a predictor velocity field, In order to obtain a divergence-free velocity field, we add a corrector field u l (s) c to the predictor field u l (s) * : where: where Ψ (s) is a scalar field that is obtained by solving the integral contravariant formulation of a Poisson-like equation: Equation (21) is discretized by means of a finite volume scheme at a control volume, ∆V i,j,k+1/2 , centered at the center of the calculation cells upper faces, in which the pressure is defined: The right-hand side of Equation (22) represents the discretization of the opposite of the divergence of the predictor velocity field. In order to overcome the check-board solution issue [35], Cannata et al. [28] computed the contravariant velocity components defined at the ∆V i,j,k lateral faces from the values obtained by the Riemann solver. By means of this strategy, the pressure-velocity decoupling problems are avoided and the shock-capturing property of the scheme is extended to the corrector step. It must be noted that the terms that appear on the left-hand side of Equation (22) represent the contravariant corrector velocity components: The corrector velocities ∂Ψ (s) ∂ξ m g m1 , and ∂Ψ (s) ∂ξ m g m3 i,j,k are discretized to solve Equation (26), in order to obtain the following set of equations, where the matrix of coefficients has 19 non-null diagonals: where the Q i,j,k+1/2 term on the right-hand side represents the divergence (changed in sign) of the predictor velocity field. See [28] for the 19 a 1 , . . . , a 19 coefficients and for more details about the numerical model.

Internal Boundary Conditions
Cannata et al. [28] adopted a procedure based on ghost cells for the boundary conditions. In particular, on the closed boundaries, a null derivative in the normal direction is assigned for the free surface elevation and a null flux is imposed through the boundary. The lateral boundary condition of the Poisson-like equation is obtained by imposing the Ψ (s) derivative as null in the direction normal to the boundary.
Unfortunately, the procedure proposed by [28] cannot be used for internal boundaries. In fact, in the corner cells of the internal object, it would be necessary to impose two different conditions for a single ghost calculation cell. Let us consider by way of example the case in which in the calculation , an internal object is defined (see Figure 1). The calculation cell [i, j, k] is a ghost cell and the internal boundary is defined on the cells' interfaces defined between [i + 1, j] and [i, j] and between [i, j + 1] and [i, j]. In order to define an internal closed boundary with the procedure proposed by [28], in the calculation cell defined in [i, j, k], a value of Ψ (s) i,j,k has to be computed, in a way to have simultaneously. and: where the , , +1 2 ⁄ term on the right-hand side represents the divergence (changed in sign) of the predictor velocity field. See [28] for the 19 1 , … , 19 coefficients and for more details about the numerical model.

Internal Boundary Conditions
Cannata et al. [28] adopted a procedure based on ghost cells for the boundary conditions. In particular, on the closed boundaries, a null derivative in the normal direction is assigned for the free surface elevation and a null flux is imposed through the boundary. The lateral boundary condition of the Poisson-like equation is obtained by imposing the Ψ ( ) derivative as null in the direction normal to the boundary.
Unfortunately, the procedure proposed by [28] cannot be used for internal boundaries. In fact, in the corner cells of the internal object, it would be necessary to impose two different conditions for a single ghost calculation cell. Let us consider by way of example the case in which in the calculation cells . In order to define an internal closed boundary with the procedure proposed by [28], in the calculation cell defined in [ , , ], a value of (Ψ ( ) ) , , has to be computed, in a way to have simultaneously.
and: These two conditions (Equations (25) and (26)) cannot be satisfied simultaneously, hence a single value of Ψ (s) i,j,k that satisfies both Equations (25) and (26) cannot be generally computed. For this reason, in the present work, an alternative strategy is proposed for the definition of the internal closed boundaries. Firstly, in the calculation cell interfaces in which the internal closed boundary is defined, the corrective velocities, i.e., the components of the gradient of Ψ (s) , which are normal to the interfaces, are explicitly set to zero. Furthermore, the derivatives of Equation (22) are discretized in a way that does not involve the values of Ψ (s) belonging to calculation cells in which the internal object is defined.
In order to better explain the two aforementioned conditions, let us consider by way of example the case shown in Figure 2 Figure 2, regarding the solution of Equation (22), in the calculation cell [i, j, k], the following corrective velocity components will be set to zero: reason, in the present work, an alternative strategy is proposed for the definition of the internal closed boundaries. Firstly, in the calculation cell interfaces in which the internal closed boundary is defined, the corrective velocities, i.e., the components of the gradient of Ψ ( ) , which are normal to the interfaces, are explicitly set to zero. Furthermore, the derivatives of Equation (22) are discretized in a way that does not involve the values of Ψ ( ) belonging to calculation cells in which the internal object is defined. In order to better explain the two aforementioned conditions, let us consider by way of example the case shown in Figure 2 Figure 2, regarding the solution of Equation (22), in the calculation cell [ , , ], the following corrective velocity components will be set to zero: Furthermore, in the example shown in Figure 2, the following derivatives: So, in Equation (22) the following conditions are imposed: Furthermore, in the example shown in Figure 2, the following derivatives: if discretized with a symmetric scheme would involve respectively Ψ (s)  (31) is discretized by means of a composition of two asymmetric derivatives: The derivative Equation (32) is discretized by means of a composition of an asymmetric derivative and a symmetric derivative: By discretizing the other derivatives that are involved in Equation (22) by means of the procedure proposed by [28], the 19 a 1 , . . . , a 19 coefficients of Equation (24) for the specific case of Figure 2 can be computed as follows: Generally speaking, in the discretization of Equation (22), the corrective velocities related to calculation cells' interfaces in which an internal boundary is defined are set to zero. Furthermore, the derivatives that would involve a grid node in which the internal object is defined are discretized asymmetrically in a way as not to involve these grid nodes.

Free-Surface Boundary Conditions
In the proposed model, the free-surface boundary conditions are imposed in order to take into account that, at the free surface, the normal stresses can be different from zero. Let us express with .. A ij the physical contravariant components of the generic tensor A ij . At the free surface, the external and internal normal stresses per unit density are, respectively: .. ..
where p atm is the atmospheric pressure. ..
Recalling that the turbulent stress tensor is related to the strain rate tensor, ..

S mn
, Equation (37) becomes: (38) Equation (38) is used for the calculation of the dynamic pressure at the free-surface, i.e., at the upper face of the top calculation cell.

Results and Discussion
In order to validate the proposed three-dimensional nonhydrostatic numerical model, we reproduced an experimental laboratory test. The test is named T3C1 and it is described in the data report "LSTF Experiments Transport by Waves and Currents & Tombolo Development Behind Headland Structures" [33] and was used as a benchmark test for numerical hydrodynamic and morphodynamic models by several authors [4,36]. The test was originally carried out on a natural beach, with a T-head groin, with a head section parallel to the shore, 4 m long and located 4 m offshore with respect to the shoreline. The T-head groin was constructed by [33] as a rubble-mound structure, with a core built with concrete blocks and the porous external layer built with riprap stones. With the proposed numerical model, it is not possible to represent porous boundaries, so the T-head groin is numerically represented by a non-porous object, with no-slip boundary conditions at its surface. In the experimental test by [33], velocity measurements were taken at approximately two-thirds of the water depth with respect to the free surface, to have a reasonable representation of the depth-averaged velocity. A random wave train was generated, with a 0.26 m significant wave height, 1.5 s period, and an incoming wave angle of 10 • with respect to the shoreline; the deep water level is 0.7 m. The experimental test was aimed to produce waves with a non-zero angle with respect to the shoreline; this causes cross-shore and long-shore currents to be produced. As reported by Gravens and Wang [33], spilling breaking wave conditions were reached in the experimental test.
In order to numerically reproduce the test carried out by [33], we extend the strategy proposed by Gallerano et al. [4] in the context of a two-dimensional Boussinesq for the computational domain, to our three-dimensional model. The plan view of the computational domain is shown in Figure 3. The red line represents the border of the laboratory basin of the experimental test by [33], the bathymetry of 0.7 m represents the line that divides the deep water region (with constant water depth) and the shallow water region, with decreasing water depth. The wave trains are generated at the north boundary in Figure 3. As stated by [4], at the east and west boundaries, the bathymetric lines are perpendicular to the boundary lines: The random wave fronts stay orthogonal to the boundaries over time, so that an open boundary (with null derivative of velocities and free-surface elevation) can be easily imposed.
A simple wet-dry technique is used to simulate the shoreline position (see [37]). The T-head groin structure is sufficiently far from the computational boundaries, so that the experimental test is numerically reproduced with good approximation. The number of cells is 500, 000, with 8 vertical layers. Random wave trains are numerically reproduced by means of the Jonswap spectrum, with a spectral enhancement factor of γ = 3.3. The mean values are calculated in the time interval 60-120 s.
The red line represents the border of the laboratory basin of the experimental test by [33], the bathymetry of 0.7 m represents the line that divides the deep water region (with constant water depth) and the shallow water region, with decreasing water depth. The wave trains are generated at the north boundary in Figure 3. As stated by [4], at the east and west boundaries, the bathymetric lines are perpendicular to the boundary lines: The random wave fronts stay orthogonal to the boundaries over time, so that an open boundary (with null derivative of velocities and free-surface elevation) can be easily imposed. A simple wet-dry technique is used to simulate the shoreline position (see [37]). The T-head groin structure is sufficiently far from the computational boundaries, so that the experimental test is numerically reproduced with good approximation. The number of cells is 500,000, with 8 vertical layers. Random wave trains are numerically reproduced by means of the Jonswap spectrum, with a spectral enhancement factor of γ = 3.3. The mean values are calculated in the time interval 60-120 s.  Figure 4 shows an instantaneous free-surface elevation field at time t = 100 s reproduced with the proposed model. The wave transformation from deep to shallow water can be deducted from Figure 4. The wave height decreases when approaching to the shoreline, thanks to the wave breaking. Furthermore, it can be noted that the effect of the structure over the wave train, in terms of diffraction and reflection, is evident.   Figure 4 shows an instantaneous free-surface elevation field at time t = 100 s reproduced with the proposed model. The wave transformation from deep to shallow water can be deducted from Figure 4. The wave height decreases when approaching to the shoreline, thanks to the wave breaking. Furthermore, it can be noted that the effect of the structure over the wave train, in terms of diffraction and reflection, is evident. bathymetry of 0.7 m represents the line that divides the deep water region (with constant water depth) and the shallow water region, with decreasing water depth. The wave trains are generated at the north boundary in Figure 3. As stated by [4], at the east and west boundaries, the bathymetric lines are perpendicular to the boundary lines: The random wave fronts stay orthogonal to the boundaries over time, so that an open boundary (with null derivative of velocities and free-surface elevation) can be easily imposed. A simple wet-dry technique is used to simulate the shoreline position (see [37]). The T-head groin structure is sufficiently far from the computational boundaries, so that the experimental test is numerically reproduced with good approximation. The number of cells is 500,000, with 8 vertical layers. Random wave trains are numerically reproduced by means of the Jonswap spectrum, with a spectral enhancement factor of γ = 3.3. The mean values are calculated in the time interval 60-120 s.  Figure 4 shows an instantaneous free-surface elevation field at time t = 100 s reproduced with the proposed model. The wave transformation from deep to shallow water can be deducted from Figure 4. The wave height decreases when approaching to the shoreline, thanks to the wave breaking. Furthermore, it can be noted that the effect of the structure over the wave train, in terms of diffraction and reflection, is evident.   Figure 5 shows a comparison between the time-averaged and depth-averaged velocity field obtained by the numerical simulation, and the experimental results obtained by [33]. From Figure 5, it can be noted that in the zone between Y = 18 m and Y = 22 m and between X = 0 m and X = 4 m, the recirculation pattern close to the structure is well rendered by the numerical model. The simulated velocity field is overall in good agreement with the measured field, even if the magnitude of the field is slightly underestimated by the numerical model in the zone between Y = 20 m and Y = 22 m. / = 0.5; and Figure 8 shows the plan view of the time-averaged velocity field near the free-surface, at a relative position from the free surface of / = 0.133. From the comparison between Figures 5  and 6, it can be noted that near the bed the currents are more offshore oriented, while from the comparison between Figures 5 and 8 it can be noted that near the free surface the currents are more onshore oriented. The offshore-oriented near-bed velocity field represented in Figure 6 can be related to an offshore-oriented undertow current. From Figures 6-8, it can be seen that the recirculation pattern located in the zone between = 18 m and = 22 m and between = 0 m and = 4 m preserves its structure over the depth.  In Figures 6-8, three plan views of the time-averaged velocity fields are shown. Figure 6 shows the plan view of the time-averaged velocity field near the bed, at a relative position from the free surface of z/z b = 0.933 (z b is the bed elevation with respect to the still free surface); Figure 7 shows the plan view of the time-averaged velocity field at a relative position from the free surface of z/z b = 0.5; and Figure 8 shows the plan view of the time-averaged velocity field near the free-surface, at a relative position from the free surface of z/z b = 0.133. From the comparison between Figures 5 and 6, it can be noted that near the bed the currents are more offshore oriented, while from the comparison between Figures 5 and 8 it can be noted that near the free surface the currents are more onshore oriented. The offshore-oriented near-bed velocity field represented in Figure 6 can be related to an offshore-oriented undertow current. From Figures 6-8, it can be seen that the recirculation pattern located in the zone between Y = 18 m and Y = 22 m and between X = 0 m and X = 4 m preserves its structure over the depth.     In Figure 9, the significant wave height computed by the proposed numerical model is compared against the experimental results obtained by Gravens and Wang [33]. The results are shown at two different sections, = 22 m (Figure 9a) and = 26 m (Figure 9b). From Figure 9, it can be seen that the decay of the wave height due to the wave breaking is well rendered by the numerical model. From Figure 9a, it can be seen that the numerical model slightly overestimates the significant wave height in the zone between = 10 m and = 14 m while the agreement outside this zone is very good. Moreover, from Figure 9b, a good agreement can be noted between the numerical and experimental results, with a slight overestimation of the significant wave height in the zone between = 8 m and = 14 m. In Figure 9, the significant wave height H S computed by the proposed numerical model is compared against the experimental results obtained by Gravens and Wang [33]. The results are shown at two different sections, Y = 22 m (Figure 9a) and Y = 26 m (Figure 9b). From Figure 9, it can be seen that the decay of the wave height due to the wave breaking is well rendered by the numerical model. From Figure 9a, it can be seen that the numerical model slightly overestimates the significant wave height in the zone between X = 10 m and X = 14 m while the agreement outside this zone is very good. Moreover, from Figure 9b, a good agreement can be noted between the numerical and experimental results, with a slight overestimation of the significant wave height in the zone between X = 8 m and X = 14 m.  the numerical simulation is shown to be in good agreement with the experimental measurements. From Figure 10b, it can be seen that at = 26 m, the long-shore current magnitude is slightly underestimated by the proposed model in the zone between the shoreline and the structure, and slightly overestimated in the offshore zone, close to the structure. Slight differences between the numerical and experimental results, in terms of long-shore currents, are mainly visible close to the structure, in particular for section = 26 m. This can be related to the fact that in the numerical model, the rubble-mound structure is represented as an impermeable object. Figure 11 shows the comparison between the cross-shore currents computed by the proposed model and those experimentally measured by Gravens and Wang [33], at two different sections, = 22 m (Figure 11a) and = 26 m (Figure 11b). Similarly with Figure 10, the results are shown in terms of the depth-averaged velocity (dashed line) and in terms of velocities at different relative positions from the free surface / = 0.133 , / = 0.400 , / = 0.667 , / = 0.933 . The agreement between the simulated depth-averaged velocity and the one measured by [33] is shown to be quite good. From Figure 11b, it has to be noted that at = 26 m, the simulated cross-shore current shows a different behavior from the measured one, in the zone between the shoreline and the structure, and overestimated in the offshore zone; furthermore, the numerical model overestimates the cross-shore current between = 6 m and = 8 m. As for the case of long-shore currents, even  Figure 10 shows the comparison between the long-shore currents computed by the proposed model and those experimentally measured by Gravens and Wang [33], at two different sections, Y = 22 m ( Figure 10a) and Y = 26 m (Figure 10b). The results are shown in terms of the depth-averaged velocity (dashed line) and in terms of velocities at different relative positions from the free surface z/z b = 0.133, z/z b = 0.400, z/z b = 0.667, z/z b = 0.933. The depth-averaged velocity resulting from the numerical simulation is shown to be in good agreement with the experimental measurements. From Figure 10b, it can be seen that at Y = 26 m, the long-shore current magnitude is slightly underestimated by the proposed model in the zone between the shoreline and the structure, and slightly overestimated in the offshore zone, close to the structure. Slight differences between the numerical and experimental results, in terms of long-shore currents, are mainly visible close to the structure, in particular for section Y = 26 m. This can be related to the fact that in the numerical model, the rubble-mound structure is represented as an impermeable object. Figure 11 shows the comparison between the cross-shore currents computed by the proposed model and those experimentally measured by Gravens and Wang [33], at two different sections, Y = 22 m (Figure 11a) and Y = 26 m (Figure 11b). Similarly with Figure 10, the results are shown in terms of the depth-averaged velocity (dashed line) and in terms of velocities at different relative positions from the free surface z/z b = 0.133, z/z b = 0.400, z/z b = 0.667, z/z b = 0.933. The agreement between the simulated depth-averaged velocity and the one measured by [33] is shown to be quite good. From Figure 11b, it has to be noted that at Y = 26 m, the simulated cross-shore current shows a different behavior from the measured one, in the zone between the shoreline and the structure, and overestimated in the offshore zone; furthermore, the numerical model overestimates the cross-shore current between X = 6 m and X = 8 m. As for the case of long-shore currents, even for the cross-shore currents, the deviation between the numerical and experimental results close to the structure (especially for section Y = 26 m) can be related to the difference between the structure and its numerical representation. for the cross-shore currents, the deviation between the numerical and experimental results close to the structure (especially for section = 26 m) can be related to the difference between the structure and its numerical representation. From Figures 10 and 11, it can be deduced that the vertical distribution of the velocity can be of primary importance in the study of this test case. In fact, while for the long-shore currents, the changes at different depths are moderate, in the case of the cross-shore currents, the velocity fields change consistently with the variation of the depth. In particular, the presence of an offshore-oriented current near the bottom and an onshore-oriented current near the free surface is evident from Figure  11.
(a) From Figures 10 and 11, it can be deduced that the vertical distribution of the velocity can be of primary importance in the study of this test case. In fact, while for the long-shore currents, the changes at different depths are moderate, in the case of the cross-shore currents, the velocity fields change consistently with the variation of the depth. In particular, the presence of an offshore-oriented current near the bottom and an onshore-oriented current near the free surface is evident from Figure 11. From Figures 10 and 11, it can be deduced that the vertical distribution of the velocity can be of primary importance in the study of this test case. In fact, while for the long-shore currents, the changes at different depths are moderate, in the case of the cross-shore currents, the velocity fields change consistently with the variation of the depth. In particular, the presence of an offshore-oriented current near the bottom and an onshore-oriented current near the free surface is evident from Figure  11.

Conclusions
In this work, we proposed a novel numerical scheme that aims to simulate coastal flows in complex geometries in the presence of the structure. The numerical model presented is based on the model proposed by Cannata et al. [28] that solves the three-dimensional Navier-Stokes equations in an integral contravariant formulation on a time-dependent curvilinear grid. Differently from the numerical scheme proposed by [28], in the proposed numerical scheme an appropriate boundary conditions treatment allows the model to represent flows around sharp-cornered structures.

Conclusions
In this work, we proposed a novel numerical scheme that aims to simulate coastal flows in complex geometries in the presence of the structure. The numerical model presented is based on the model proposed by Cannata et al. [28] that solves the three-dimensional Navier-Stokes equations in an integral contravariant formulation on a time-dependent curvilinear grid. Differently from the numerical scheme proposed by [28], in the proposed numerical scheme an appropriate boundary conditions treatment allows the model to represent flows around sharp-cornered structures. Furthermore, the proposed scheme adopts an evaluation of the non-zero dynamic pressure at the free surface. The use of boundary-conforming grids combined with the adoption of the proposed numerical scheme makes it possible to represent complex natural morphologies and coastal structures. In order to test the model, a simulation consisting in a random wave train approaching a natural shore with a T-head groin structure was carried out and the numerical results were compared against the experimental measurements. From the comparison between the numerical and experimental results, we verified the ability of the model to reproduce three-dimensional hydrodynamic flows in complex domains with structures.