Nonlinear waves and nearshore currents over variable bathymetry in curve-shaped coastal areas

Abstract

Irregular coastlines and variable bathymetry produce nonlinear effects on wave propagation which play a significant role on the formation of nearshore currents. To protect the coastline from the erosional action of nearshore currents, it is usual to adopt coastal defence works such as submerged breakwaters. If properly designed, they give rise to circulation patterns capable to induce sedimentation of suspended material at the nearshore region. To numerically simulate the hydrodynamic effects of submerged breakwaters in irregular coastal areas, we use a numerical model which is based on an integral contravariant formulation of the three-dimensional Navier–Stokes equations in a time-dependent coordinate system. These equations are numerically solved by a non-hydrostatic shock-capturing numerical scheme which is able to simulate the wave propagation from deep water to the shoreline, including the surf zone and swash zone.

1 Introduction

The coastal regions can be characterised by a coastline of irregular shape and high variability in the water depth, where transformation of nonlinear waves takes place due to processes such as refraction, shoaling and wave breaking. Nonlinear interaction between wave motion and variable bathymetries produce gradients in mean water level which drive nearshore currents such as long-shore currents and rip currents. Coastal circulation patterns characterised by offshore-directed currents can occur in curve-shaped coastal areas, due to natural seabed variability or water depth variations produced by submerged breakwaters. The presence of coastal defence structures, such as submerged breakwaters, can modify coastal sediment transport since it can reduce the wave energy and can also modify the coastal circulation patterns. The presence of submerged breakwaters results in a local reduction of the water depth which induces the waves to break farther from the shoreline than in the absence of the barriers and dissipate part of the incident wave energy. Such a phenomenon can have a local positive effect on the coastal protection from erosion since it induces a reduction in the wave capability to resuspend sediment from the bottom. Furthermore, the presence of the submerged breakwaters induces nearshore circulation patterns that in some cases can favour the sedimentation of the sediment carried by the currents, with consequent advancement of the coastline, while in other cases can increase the offshore-directed sediment transport, with consequent retreat of the coastline. A system of submerged breakwaters separated by gaps can induce two different current circulation patterns in the lee of the barriers: an erosive or an accretive circulation pattern. The erosive one is characterised by a two-cell circulation pattern. In this circulation pattern, the wave breaking over the structure produces a wave setup in the lee of the barrier which is greater than the wave setup in the gaps; consequently, the inflow of water over the barrier leaves the system at the gaps by means of diverging currents in the lee of the structure, generating two circulation cells. Such a circulation pattern induces shoreline erosion since the sediment resuspended by breaking waves near the shoreline can be transported out of the protected area by the gaps. The accretive circulation pattern is characterised by a four-cell circulation system. Like the preceding case, the onshore flow over the barrier diverges towards the gaps, generating two circulation cells. Unlike the preceding case, just near the shoreline, the wave setup in the lee of the barrier is lower than the wave setup in the gaps. Consequently, close to the shoreline, the alongshore flow direction may reverse and be directed towards the centreline of the barrier, generating two additional circulation cells rotating in the opposite direction with respect to the previous ones. In this case, the sediment resuspended by breaking waves near the shoreline tends to settle in the protected area. The above makes clear that, to numerically verify the effect of a system of submerged breakwaters, it is necessary to simulate the nonlinear interactions between wave motion and variable bathymetry. This requires numerical models able to represent the evolution of the waveform, wave breaking, energy dissipation and wave-induced nearshore currents. In the literature, many authors used a numerical model to simulate the wave breaking (Roelvink et al. 2009; Yamazaki et al. 2009; Zijlema et al. 2011; Stelling and Duinmeijer 2003). Other authors used numerical models based on the depth-averaged motion equations, obtained by assuming a simplified distribution of the hydrodynamic quantities along the vertical direction. Among these, the most advanced (Roeber et al. 2010; Tonelli and Petti 2012; Shi et al. 2012; Gallerano et al. 2016a; Cannata et al. 2018a, b) make use of the extended Boussinesq equations to simulate the propagation of the waves in deep water, where dispersive effects are dominant, and make use of the nonlinear shallow water equations to simulate shoaling and wave breaking in the surf zone, where the nonlinear effects are dominant. One of the main drawbacks of the above-mentioned models is the fact that they require the introduction of an “a priori” criterion to established where switching from one system of equations to the other. Furthermore, these models are not able to adequately represent fully three-dimensional velocity fields. The simulation of three-dimensional velocity fields and free surface elevation can be carried out by numerically integrating the Navier–Stokes equations in fully three-dimensional form. The main drawback of the numerical models based on this approach is the significant increase of the computational time with respect to the models based on the depth-averaged motion equations. In fact, in general, in the three-dimensional non-hydrostatic numerical models the computational time increases more than linearly with the number of layers used for the vertical discretization of the physical domain: most of the computational time is spent to calculate the dynamic component of the pressure by iterative numerical methods. Some of the most recent three-dimensional non-hydrostatic numerical models make use of coordinate transformation in the vertical direction (σ-coordinate transformation), by which the vertical Cartesian coordinate, z, is expressed as a function of a time-dependent curvilinear coordinate, σ, that follows the free surface movements. The adoption of this strategy allows to obtain numerical models with excellent dispersive properties even using few layers in the vertical direction (Lin and Li 2002; Young and Wu 2010). A recent evolution of the three-dimensional models that adopt a time-varying coordinate makes use of shock-capturing numerical schemes that are able to explicitly simulate wave breaking. In Bradford (2011), Ma et al. (2012), Cannata et al. (2017), and Gallerano et al. (2017), the Navier–Stokes equations are written in a coordinate system in which the horizontal coordinates coincide with the Cartesian coordinates and the vertical one is a time-dependent curvilinear coordinate. These models are able to simulate in fully three-dimensional form the wave propagation from deep water, where dispersive effects are dominant, to surf and swash zones, where the nonlinearity of the motion equations allows to represent the shoaling and the wave breaking, and do not require to switch off any term of the equations. The main drawback of these models is that their application is limited to quite simple geometries that can be represented by Cartesian rectangular grids in the horizontal plane. The simulation of three-dimensional flow fields and free surface elevation in coastal regions characterised by complex morphology requires numerical models that make use of unstructured grids (Sørensen et al. 2004) or boundary conforming grids and generalised curvilinear coordinates. In the context of curvilinear coordinates, the most general approach is the one based on motion equations expressed in contravariant formulation in time-dependent generalised curvilinear coordinates that are able to reproduce curve-shaped coastal areas. The existing three-dimensional models that adopt this approach are based on the numerical integration of the Navier–Stokes equations expressed in differential form (Ogawa and Ishiguro 1987; Luo and Bewley 2004). The differential formulation proposed by these authors does not allow the adoption of shock-capturing numerical schemes and thus cannot be used for the simulation of breaking waves. In this paper, to evaluate the hydrodynamic effects produced by a system of submerged breakwaters in curve-shaped coastal areas, we use a three-dimensional numerical model recently proposed by Cannata et al. (2019) which is based on an integral, contravariant form of the Navier–Stokes equations in a system of time-dependent generalised curvilinear coordinates. The adopted integral, contravariant formulation allows discretization of the irregular physical domain that represents a curve-shaped coastal area by time-varying curvilinear cells which adapt to the free surface evolution. The adoption of a time-dependent coordinate transformation allows to numerically integrate the motion equations in a transformed space in which the computational cells are regular and fixed. The finite-volume shock-capturing numerical scheme adopted in this paper makes it possible to simulate the propagation of the discontinuities in the solution associated with wave breaking. The paper is structured as follows. In Sect. 2 is presented, in synthetic form, the procedure used to obtain the integral and contravariant formulation of the equations of motion in a system of curvilinear coordinates that vary over time; in Sect. 3, the numerical model has been validated with experimental tests and has been applied to simulate numerically the effects of submerged barriers on the propagation of nonlinear waves on variable bathymetries in curve-shaped coastal area; in Sect. 4 are shown the conclusions.

2 Governing equations

Let us consider a Cartesian reference where $\left( {x^{1} ,x^{2} ,x^{3} } \right)$ are the spatial coordinates and $t$ is the time. Let us consider a time-varying control volume $\Delta V\left( t \right)$ that coincides with a material fluid volume, i.e. a time-varying volume which moves with the fluid and always encloses the same fluid particles. Let also $\Delta V_{1} \left( t \right)$ be a time-varying control volume bounded by a surface, of area $\Delta A_{1} \left( t \right)$, every point of which moves with a velocity $\vec{v}$ that is different from the fluid velocity.

Let $H\left( {x^{1} ,x^{2} ,t} \right) = h\left( {x^{1} ,x^{2} } \right) + \eta \left( {x^{1} ,x^{2} ,t} \right)$ be the total water depth, where h is the undisturbed water depth and $\eta$ is the free surface elevation with respect to the undisturbed level. We indicate the acceleration due to gravity by $G$ and we split the pressure $p$ into a hydrostatic part, $\rho G\left( {\eta - x^{3} } \right),$ and a dynamic one, $q$. To accurately represent the complex geometry of a curve-shaped coastal area and to follow the wave-induced free surface evolution, we consider the following time-dependent transformation from the Cartesian system of coordinates, $\left( {x^{1} ,x^{2} ,x^{3} } \right)$, to the moving curvilinear system of coordinates, $\left( {\xi^{1} ,\xi^{2} ,\xi^{3} } \right),$

$$\xi^{1} = \xi^{1} \left( {x^{1} ,x^{2} ,x^{3} } \right);\;\xi^{2} = \xi^{2} \left( {x^{1} ,x^{2} ,x^{3} } \right);\;\xi^{3} = \frac{{x^{3} + h\left( {x^{1} ,x^{2} } \right)}}{{H\left( {x^{1} ,x^{2} ,t} \right)}};\;\tau = t,$$

(1)

where the horizontal curvilinear coordinates $\xi^{1}$ and $\xi^{2}$ conform to the horizontal boundaries of the physical domain and the vertical coordinate $\xi^{3}$ varies over time to adapt to the free surface movements. This coordinate transformation basically maps the irregular, varying domain in the physical space to a regular, fixed domain in the transformed space, where $\xi^{3}$ spans from 0 to 1. Let $\vec{g}_{\left( l \right)} = \partial \vec{x}/\partial \xi^{l}$ be the covariant base vectors and $\vec{g}^{\left( l \right)} = \partial \xi^{l} /\partial \vec{x}$ the contravariant base vectors ($l = 1,3$). The covariant and contravariant metric coefficients are defined, respectively, by $g_{lm} = \vec{g}_{\left( l \right)} \cdot \vec{g}_{\left( m \right)}$ and $g^{lm} = \vec{g}^{\left( l \right)} \cdot \vec{g}^{\left( m \right)}$ ($l,m = 1,3$), where the dot “$\cdot$” indicates the scalar product between vectors. The Jacobian of the transformation is given by $\sqrt g = \sqrt {\left| {g_{lm} } \right|}$, where $|\, |$ denotes the determinant of the covariant metric coefficients $g_{lm}$. The transformation relationships between the components of the generic vector $\vec{b}$ in the Cartesian coordinate system and its contravariant and covariant components, $b^{l}$ and $b_{l}$ in the curvilinear coordinate system are

$$b^{l} = \vec{g}^{\left( l \right)} \cdot \vec{b};\;\vec{b} = b^{l} \vec{g}_{\left( l \right)} ;\;b_{l} = \vec{g}_{\left( l \right)} \cdot \vec{b};\;\vec{b} = b_{l} \vec{g}^{\left( l \right)} ,$$

(2)

in which (and hereinafter) the summation convention, where repeated indices are automatically summed over, is employed.

The mass continuity equation and the momentum equation written in integral and contravariant forms, respectively, read as

$$\frac{d}{{{\text{d}}\tau }}\mathop \int \nolimits_{{\Delta V_{1} \left( \tau \right)}}^{{}} \rho {\text{d}}V + \mathop \int \nolimits_{{\Delta A_{1} \left( \tau \right)}}^{{}} \rho \left( {u^{m} - v^{m} } \right)n_{m} {\text{d}}A = 0,$$

(3)

$$\frac{d}{{{\text{d}}\tau }}\mathop \int \nolimits_{{\Delta V_{1} \left( \tau \right)}}^{{}} \overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}_{\left( k \right)} \rho u^{k} {\text{d}}V + \mathop \int \nolimits_{{\Delta A_{1} \left( \tau \right)}}^{{}} \overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}_{\left( k \right)} \rho u^{k} \left( {u^{m} - v^{m} } \right)n_{m} {\text{d}}A = \mathop \int \nolimits_{{\Delta V_{1} \left( \tau \right)}}^{{}} \overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}_{\left( k \right)} \rho f^{k} {\text{d}}V + \mathop \int \nolimits_{{\Delta A_{1} \left( \tau \right)}}^{{}} \overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}_{\left( k \right)} T^{km} n_{m} {\text{d}}A,$$

(4)

where $\rho$ is the fluid density, $u^{m}$ ($m = 1,3$) are the contravariant components of the fluid velocity vector, $v^{m}$ are the contravariant components of the velocity vector with which the points belonging to the surface area $\Delta A_{1} \left( \tau \right)$ move, $n_{m}$ are the covariant components of the outward unit vector normal to the surface of area $\Delta A_{1} \left( \tau \right)$, $f^{k}$ are the contravariant components of the vector of the external body forces per unit mass, $T^{km}$ are the contravariant components of the stress tensor and $\overrightarrow {{\widetilde{g}}}^{\left( l \right)} = \vec{g}^{\left( l \right)} \left( {\xi_{0}^{1} ,\xi_{0}^{2} ,\xi_{0}^{3} } \right)$, being $P_{0} \left( {\xi_{0}^{1} ,\xi_{0}^{2} ,\xi_{0}^{3} } \right)$ a point belonging to the control volume. Further details can be found in Cannata et al. (2019).

Let $\sqrt {g_{0} } = \vec{k} \cdot \left( {\vec{g}_{\left( 1 \right)} {\bigwedge }\vec{g}_{\left( 2 \right)} } \right)$, where $\vec{k}$ indicates the vertical unit vector and ${\bigwedge }$ indicates the vector product. It is not difficult to verify that in the specific case of the above-mentioned transformation, the Jacobian of the transformation can be written in the form $\sqrt g = H\sqrt {g_{0} }$. This makes it possible to write an original three-dimensional integral contravariant conservative form of the momentum equation, in which the conserved variables are given by the cell averaged product between the water depth $H$ and the three contravariant components of the punctual velocity. To this end, let us define the following cell averaged values in the transformed space

$$\begin{aligned} \bar{H} & = \frac{1}{{\Delta A_{0}^{3} \sqrt {g_{0} } }}\mathop \int \nolimits_{{\Delta A_{o}^{3} }}^{{}} H\sqrt {g_{0} } {\text{d}}\xi^{1} {\text{d}}\xi^{2} , \\ \overline{{{\text{Hu}}^{l} }} & = \frac{1}{{\Delta V_{0} \sqrt {g_{0} } }}\mathop \int \nolimits_{{\Delta V_{0} }}^{{}} \overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}_{\left( k \right)} u^{k} H\sqrt {g_{0} } {\text{d}}\xi^{1} {\text{d}}\xi^{2} {\text{d}}\xi^{3} , \\ \end{aligned}$$

(5)

where a restrictive condition has been introduced in the control volume $\Delta V_{1} \left( \tau \right)$ that defines $\Delta V_{1} \left( \tau \right)$ as the volume of a physical space that is bounded by surfaces lying on the curvilinear coordinate surfaces. In the curvilinear coordinate system, $\Delta V_{1} \left( \tau \right) = \mathop \int \nolimits_{{\Delta V_{0} }}^{{}} \sqrt g {\text{d}}\xi^{1} {\text{d}}\xi^{2} {\text{d}}\xi^{3}$, where $\Delta V_{0} = \Delta \xi^{1} \Delta \xi^{2} \Delta \xi^{3}$ indicates the corresponding volume in the transformed space. Analogously, in the curvilinear coordinate system, the area of a surface of the physical space that lies on the coordinate surface in which $\xi^{\alpha }$ is constant is $\Delta A_{1}^{\alpha } \left( \tau \right) = \mathop \int \nolimits_{{\Delta A_{0}^{\alpha } }}^{{}} \left| {\vec{g}_{\left( \beta \right)} {\bigwedge }\vec{g}_{\left( \gamma \right)} } \right|{\text{d}}\xi^{\beta } {\text{d}}\xi^{\gamma }$, where $\Delta A_{0}^{\alpha } = \Delta \xi^{\beta } \Delta \xi^{\gamma }$ indicates the corresponding area in the transformed space. It must be noted that the volume $\Delta V_{1} \left( \tau \right)$ and the surfaces $\Delta A_{1}^{\alpha } \left( \tau \right)$ are functions of time, because they are expressed as functions of the base vectors, $\vec{g}_{\left( l \right)}$, and the Jacobian of the transformation, $\sqrt g$, whose values change over time as the curvilinear coordinates follow the displacements of the free surface. Conversely, the volume $\Delta V_{0}$ and the areas $\Delta A_{0}^{\alpha }$ are not time dependent.

Since $H$ and $\sqrt {g_{0} }$ do not depend on $\xi^{3}$, the cell averaged value of the water depth $\overline{H}$ represents a two-dimensional quantity given by the averaged value of the water depth $H$ on the base area of a water column, $\Delta A^{3} = \Delta A_{0}^{3} \sqrt {g_{0} }$. Since the velocity components $u^{l}$ are dependent on the spatial coordinates $\xi^{1}$, $\xi^{2}$, $\xi^{3}$, the cell averaged value $\overline{{Hu^{l} }}$ represents a three-dimensional quantity that varies along the three spatial coordinates.

By adopting the volume $\Delta V\left( \tau \right)$ defined above as control volume and using the definition of the cell averaged values given in Eq. (5), in the transformed space, Eq. (4) reads

$$\begin{aligned} \frac{{\partial \overline{{{\text{Hu}}^{l} }} }}{\partial \tau } & = - \frac{1}{{\Delta V_{0} \sqrt {g_{0} } }}\mathop \sum \limits_{\alpha = 1}^{3} \left\{ {\mathop \int \nolimits_{{\Delta A_{o}^{\alpha + } }}^{{}} \left[ {\overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}_{\left( k \right)} {\text{Hu}}^{k} \left( {u^{\alpha } - v^{\alpha } } \right) + \overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}^{\left( \alpha \right)} {\text{GH}}^{2} } \right]\sqrt {g_{0} } {\text{d}}\xi^{\beta } {\text{d}}\xi^{\gamma } } \right. \\ & \quad - \left. {\mathop \int \nolimits_{{\Delta A_{o}^{\alpha - } }}^{{}} \left[ {\overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}_{\left( k \right)} {\text{Hu}}^{k} \left( {u^{\alpha } - v^{\alpha } } \right) + \overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}^{\left( \alpha \right)} {\text{GH}}^{2} } \right]\sqrt {g_{0} } {\text{d}}\xi^{\beta } {\text{d}}\xi^{\gamma } } \right\} \\ & \quad + \frac{1}{{\Delta V_{0} \sqrt {g_{0} } }}\mathop \sum \limits_{\alpha = 1}^{3} \left\{ {\mathop \int \nolimits_{{\Delta A_{o}^{\alpha + } }}^{{}} \overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}^{\left( \alpha \right)} {\text{Gh}}H\sqrt {g_{0} } {\text{d}}\xi^{\beta } {\text{d}}\xi^{\gamma } - \left. {\mathop \int \nolimits_{{\Delta A_{o}^{\alpha - } }}^{{}} \overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}^{\left( \alpha \right)} {\text{Gh}}H\sqrt {g_{0} } {\text{d}}\xi^{\beta } {\text{d}}\xi^{\gamma } } \right\}} \right. \\ & \quad + \frac{1}{{\Delta V_{0} \sqrt {g_{0} } }}\mathop \sum \limits_{\alpha = 1}^{3} \left\{ {\mathop \int \nolimits_{{\Delta A_{o}^{\alpha + } }}^{{}} \overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}_{\left( k \right)} \frac{{T^{k\alpha } }}{\rho }H\sqrt {g_{0} } {\text{d}}\xi^{\beta } {\text{d}}\xi^{\gamma } - \left. {\mathop \int \nolimits_{{\Delta A_{o}^{\alpha - } }}^{{}} \overrightarrow {{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}_{\left( k \right)} \frac{{T^{k\alpha } }}{\rho }H\sqrt {g_{0} } {\text{d}}\xi^{\beta } {\text{d}}\xi^{\gamma } } \right\}} \right. \\ & \quad - \frac{1}{{\Delta V_{0} \sqrt {g_{0} } }}\mathop \int \nolimits_{{\Delta V_{0} }}^{{}} \overline{{\widetilde{g}}}^{\left( l \right)} \cdot \vec{g}^{\left( m \right)} \frac{\partial q}{{\partial \xi^{m} }}H\sqrt {g_{0} } {\text{d}}\xi^{1} {\text{d}}\xi^{2} {\text{d}}\xi^{3} , \\ \end{aligned}$$

(6)

where $\Delta A_{o}^{\alpha + }$ and $\Delta A_{o}^{\alpha - }$ indicate the contour surfaces of the volume $\Delta V_{0}$ on which $\xi^{\alpha }$ is constant and which are located at the larger and the smaller values of $\xi^{\alpha }$, respectively. Here the indices α, β and γ are cyclic. $T^{k\alpha }$ is now the stress tensor in which the pressure is omitted, the gradient of the hydrostatic pressure is split into two parts using $\eta = H - h$ and the last integral on the right-hand side of Eq. (6) is related to the gradient of the dynamic pressure, $q$. Let us integrate the continuity Eq. (3) over a vertical water column (between the bottom and the free surface) which is bounded by coordinate surfaces. By applying the impenetrability condition, $u^{3} = 0$, to the bottom (where $v^{3} = 0$, because the bottom is fixed), and by applying the kinematic condition, $v^{3} = u^{3} ,$ to the free surface (that moves with the same velocity as the fluid), in the transformed space, we obtain

$$\frac{{\partial \overline{H} }}{\partial \tau } = \frac{1}{{\Delta A_{o}^{3} \sqrt {g_{0} } }}\mathop \sum \nolimits_{\alpha = 1}^{2} \left[ {\mathop \int \nolimits_{0}^{1} \mathop \int \nolimits_{{\Delta \xi_{o}^{\alpha + } }} u^{\alpha } H\sqrt {g_{0} } {\text{d}}\xi^{\beta } {\text{d}}\xi^{3} - \mathop \int \nolimits_{0}^{1} \mathop \int \nolimits_{{\Delta \xi_{o}^{\alpha - } }} u^{\alpha } H\sqrt {g_{0} } {\text{d}}\xi^{\beta } {\text{d}}\xi^{3} } \right],$$

(7)

where $\Delta \xi_{o}^{\alpha + }$ and $\Delta \xi_{o}^{\alpha - }$ (with $\alpha = 1,2$) indicate the contour line of the surface element $\Delta A_{o}^{3}$ on which $\xi^{\alpha }$ is constant and which are located at the larger and at the smaller values of $\xi^{\alpha } ,$ respectively.

Equation (7) represents the contravariant form of the governing equation for the free surface elevation, in which the unknown variable is the cell averaged water depth $\overline{H}$. Consequently, the unknown variables of integral Eqs. (6) and (7) are the cell averaged values of the water depth $\overline{H}$ and the product between the water depth and the three-dimensional flow velocity contravariant components $\overline{{{\text{Hu}}^{l} }}$.

It must be noted that in Eq. (6), the momentum balance equation is expressed in a contravariant integral form in which Christoffel symbols are absent. This result has been achieved by adopting the procedure proposed by Gallerano and Cannata (2011) for the depth-averaged motion equations and extended to the three-dimensional motion equations by Cannata et al. (2019). According to this procedure, the integral form of the vectorial momentum balance equation is projected onto the directions of the contravariant base vectors defined at the centre of the curvilinear control volume. As demonstrated in previous works (Gallerano and Cannata 2011; Gallerano et al. 2012, 2016b), the absence of the Christoffel symbols in the governing equations significantly reduces the sensitivity of the numerical solution to the irregularity of the employed curvilinear grid.

Equations (6) and (7) are discretized on a curvilinear non-orthogonal grid in which a vertical staggering proposed by Stelling and Zijlema (2003) is adopted. According to this grid arrangement, the flow velocity components are defined at the centre of the grid cells, while the dynamic pressure is defined on the upper surface of each grid cell. In this way, the pressure boundary conditions are directly set to zero at the free surface and the dynamic component of the pressure is calculated by a predictor–corrector method: in the first step (predictor step), an approximate velocity field is calculated by eliminating the dynamic (non-hydrostatic) pressure term from Eq. (6) and integrating the resulting equation by a shock-capturing scheme in which an HLL-Riemann solver is adopted. In the second step, to take into account the dynamic pressure, a corrector velocity field is obtained by solving a Poisson-like equation by an iterative method that makes use of an alternating zebra line Gauss–Seidel method and a multigrid technique. The sum of the predictor and corrector velocity fields gives a final, divergence-free, velocity field. This corrected non-hydrostatic velocity flow field is used to update the water depth (and the free surface elevation) by means of Eq. (7). More details about the numerical procedure can be found in Cannata et al. (2019).

3 Applications

3.1 Wave train propagating at a varying depth

To check the ability of the proposed model to simulate the shoaling and breaking wave processes, the experimental test performed by Stive (1980) is here numerically reproduced. In this test, we simulate the propagation of waves whose initial height and period are, respectively, 0.156 m and 1.79 s. In the experimental test, the wave flume was 55 m long and 1 m wide; monochromatic waves were generated on a flat bottom where the water depth was 0.85 m. From the deep water region, a 6-m-long plane slope reduced the water depth to 0.7 m; after a 10-m-long flat bottom a second plane of constant slope 1:40 was placed (Fig. 1). We impose closed boundary conditions (zero normal velocity and zero gradient for tangential velocity and free surface elevation) at the lateral boundaries. No-slip bottom boundary conditions are used. The computational grid adopted for the numerical simulations consists of 104,448 grid cells in the plane, while in the vertical direction 6 layers are used. The time step is 0.002 s. The turbulence stress tensor is estimated by the Smagorinsky sub-grid model in which the Smagorinsky coefficient is set to $0.1$. The numerical simulation is carried on for 250 wave periods. After the first 50 wave periods, during which a quasi-stationary condition is achieved, the time average of the hydrodynamic quantities is performed over 200 wave periods. Aiming to demonstrate the independency of the results from the grid distortion, the numerical simulations are conducted on two different grids: a Cartesian grid and distorted curvilinear one. This distorted curvilinear grid is obtained by deforming a regular Cartesian grid by following the procedure proposed by Visbal and Gaitonde (2002). Figure 2 shows a plan view of the distorted grid and a detailed view of the computational domain. The root means square error, $\sigma_{{v_{y} }}$, obtained by the simulation carried out on the distorted grid exceeds by 1% the one obtained on the Cartesian grid online. Equally good results are obtained also for $\sigma_{wh}$ and $\sigma_{mwl}$ that represent, respectively, the root mean square error of the difference between the wave height and the mean water level carried out in the distorted grid and the ones carried out in the Cartesian grid (see Table 1).

Fig. 1

Topography of the bottom

Fig. 2

Plane view and detailed plane view of the highly distorted curvilinear computational grid

Table 1 Numerical root mean square error of the velocity component ($\sigma_{{v_{y} }}$), of average wave height ($\sigma_{wh}$) and of still water level ($\sigma_{mwl}$)

Figure 3 shows a 3D view of the wave propagation on the sloping beach obtained by the proposed model on the distorted computational grid shown in Fig. 2. In Fig. 4, a plane view of the above-mentioned wave field is shown. By observing these figures, it is possible to deduce that the wave fronts are maintained normal to the x-axis along the channel even in the presence of distorted computational cells.

Fig. 3

Instantaneous wave field obtained with the highly distorted curvilinear computational grid

Fig. 4

Detailed plane view of the instantaneous wave field obtained with the highly distorted curvilinear computational grid

Figure 5 shows the wave height and mean water level obtained by the numerical model on the distorted grid and the corresponding experimental data.

Fig. 5

Comparison between the experimental data (circles) and the numerical results (solid line) obtained with the highly distorted curvilinear computational grid in terms of time-averaged wave height and mean water level

3.2 Monochromatic waves over submerged breakwaters

Submerged breakwaters are designed to reduce the incident wave energy to protect the coasts against erosion and port channels from sand deposition. Performance of a submerged breakwater particularly depends on its geometrical and structural characteristics, the location of the breakwater (the water depth at the toe and its distance from the coastline), wave–structure interaction, sediment transport and local morphodynamics. Such structures absorb some of the wave energy by causing the waves to break prematurely. The remaining energy is partly reflected and partly transmitted shoreward. Furthermore, submerged breakwaters with rip channels are designed in such a way to induce, immediately downstream the breakwaters, gradients of the average water surface elevation such as to produce accretion (Ranasinghe et al. (2010)). The design characteristics of a breakwater structure are important in determining how the structure will impact the shoreline as evidenced by Ranasinghe (2006–2010).

From Fig. 6, it can be seen that close to the submerged breakwaters, due to the breaking, there is an increase in the free surface elevation compared to the still water level (setup $\eta_{1}$) and a reduction in the energy of the wave motion.

Fig. 6

Schematic cross section and profile of the setup

Immediately behind the breakwaters, the wave breaking stops because of the increase in the depth, then the wave (which is now characterised by a new height $H_{2}$) restarts to shoal until it breaks near the coastline. At the coastline, the total setup with respect to the still water level is given by the sum of the two successive setups which the wave has been subjected to ($\eta_{\text{tot}} = \eta_{1} + \eta_{2}$).

At the rip channels placed between the barriers, the incident waves propagate without being directly affected by the presence of the barriers and break as they approach the coastline. Let $\eta_{H}$ be the setup, evaluated near the coastline, behind the rip channel. The modifications produced on the incident wave motion by the presence of the barriers induce circulations which can be summarised in two different circulation patterns. Given that the behaviour immediately behind the barriers does not vary between these two configurations, always determining an increase $\eta_{1}$ > $\eta_{H}$ and, therefore, an average water surface elevation which is greater at the back of the barriers than at the rip channel, a primary rip current is generated towards the offshore region. On the contrary, in the area protected by the breakwaters, two different types of induced circulation patterns can be identified. The first type is characterised (as shown in Fig. 7a) by the fact that, near the shoreline, the average water surface elevation is lower at the rip channel than behind the breakwaters ($\eta_{H} < \eta_{\text{tot}} \left( {\frac{{\eta_{\text{tot}} }}{{\eta_{H} }} < 1} \right)$), so as to induce a secondary circulation with direction equal to that of the primary circulation. This type of circulation pattern is known as two-cell circulation. The two-cell circulation produces a dragging of the sediments towards the offshore region, thus favouring the erosion of the seabed at the coastline with possible regression of the latter.

Fig. 7

Erosive circulation pattern (a) and accretive circulation pattern (b)

The second type is characterised (as shown in Fig. 7b) by the fact that, near the shoreline, the average water surface elevation is greater at the rip channel than behind the breakwaters ($\eta_{H} > \eta_{\text{tot}} \left( {\frac{{\eta_{\text{tot}} }}{{\eta_{H} }} < 1} \right)$), so as to induce a secondary circulation with opposite direction to that of the primary circulation and thus favouring the accumulation of sediments in the protected area and the advancement of the coastline. This type of circulation pattern is known as four-cell circulation.

In order to validate the proposed model, the laboratory test conducted by Haller et al. (1997, 2002) is numerically reproduced. The experimental arrangement consists of a basin which is 17.0 m long and 18.0 m wide (Fig. 8). The beach has an initial steep slope (1:5) followed by a lower slope (1:30). There are three submerged breakwaters: the central bar is 7.32 m long, while the two side bars are 3.66 m long (about half of the central one) and spaced each other by gaps which are 1.82 m wide. We impose closed boundary conditions (zero normal velocity and zero gradient for tangential velocity and free surface elevation) at the lateral boundaries. No-slip bottom boundary conditions are used. The computational grid adopted for the numerical simulations consists of 88,320 grid cells in the plane, while in the vertical direction 6 layers are used. The time step is 0.0025 s. The turbulence stress tensor is estimated by the Smagorinsky sub-grid model in which the Smagorinsky coefficient is set to $0.1$. The numerical simulation is carried on for 365 wave periods. After the first 65 wave periods, during which a quasi-stationary condition is achieved, the time average of the hydrodynamic quantities is performed over 300 wave periods.

Fig. 8

Plane view (a) and cross section of the basin (b)

The numerical model is tested by numerically reproducing tests B and D. The reproduced wave train is monochromatic, regular and normally incident to the shoreline. Table 2 shows the initial wave height, the wave period, the angle of incidence, the average water depth at the bar crest and the coastal position.

Table 2 Experimental conditions: the initial wave height, the wave period, the angle of incidence, the still water depth and the coastal position

In Fig. 9, the time-averaged flow velocity field (in which only one out of every four vectors is represented) obtained near the bottom in accord with the experimental measurements, by the proposed three-dimensional numerical model, is shown. The comparison between the proposed three-dimensional numerical model and the experimental results for two different wave conditions is shown in Fig. 9. The numerical simulation related to test B returns a four-cell circulation pattern (Fig. 9a). On the contrary, test D shows a two-cell circulation pattern (Fig. 9b). Moreover, Fig. 9a, b shows an excellent correspondence between the numerical results and the experimental data.

Fig. 9

Time-averaged velocity fields for Test B (a) and Test D (b). Comparison between the numerical results and the experimental measurements

3.3 Nonlinear waves in a curve-shaped coastal area

In this subsection, we verify the ability of the proposed model to numerically reproduce wave propagation, wave breaking and induced nearshore circulation due to the variable bathymetry in a curve-shaped coastal area. To this end, we reproduce a laboratory experiment carried out by Hamm (1992). The tank used by Hamm (1992) measured 30 m by 30 m; the sea bed consisted of a plane beach sloping at 1 in 30, with a rip channel excavated in the centre which produces a curve-shaped coastline. The bathymetry is given by

$$z\left( {x,y} \right) = \left\{ {\begin{array}{*{20}l} {\begin{array}{*{20}l} {0.5 } \hfill \\ {} \hfill \\ \end{array} } & {x \le 7} \\ {0.1 - \frac{18 - x}{30}\left[ {1 + 3{ \exp }\left( { - \frac{18 - x}{30}} \right){ \cos }^{10} \left( {\frac{{\pi \left( {15 - y} \right)}}{30}} \right)} \right]} & {7 < x < 25} \\ {0.1 + \frac{18 - x}{30}} & {x \le 25} \\ \end{array} } \right..$$

Since the basin is symmetric with respect to the y-axis, in the numerical simulation, it is sufficient to reproduce only one half of the basin. Consequently, at the line of symmetry (the centreline of the rip channel), we impose a reflective boundary condition. At the opposite lateral boundary, we impose closed boundary conditions (zero normal velocity and zero gradient for tangential velocity and free surface elevation). No-slip bottom boundary conditions are used. The computational grid adopted for the numerical simulations consists of 102,400 grid cells in the plane, while in the vertical direction 6 layers are used. The time step is 0.001 s. The turbulence stress tensor is estimated by the Smagorinsky sub-grid model in which the Smagorinsky coefficient is set to $0.1$. The numerical simulation is carried on for 250 wave periods. After the first 50 wave periods, during which a quasi-stationary condition is achieved, the time average of the hydrodynamic quantities is performed over 200 wave periods. To this end, we use a curvilinear boundary conforming grid which in the horizontal directions reproduces the curve-shaped coastline. Figure 10a shows a plan view of the curvilinear computational grid and bottom variation, in which only one out of every five coordinate lines is shown. Figure 10b, c shows the two vertical sections of the computational domain, one located in the rip channel ($y_{\text{A}} =$ 14.9625 m) and one in the plane beach ($y_{\text{B}} =$ 1.9875 m). We numerically reproduce a regular wave train with period $T =$ 1.25 s and height $H =$ 0.07 m.

Fig. 10

a Bathymetry (only one out of every five coordinate lines is shown). b, c Bottom profiles in sections A and B

Figure 11 shows the comparison between the wave heights obtained by the numerical simulations and the experimental data by Hamm (1992) in Section A and Section B, as presented by Sørensen et al. (1998). Figure 11 shows a good correspondence between experimental measurements and numerical results in both sections.

Fig. 11

Comparison between the computed (solid line) and measured cross-shore variation of the wave height (crosses represent significant wave heights $H_{1/3}$ and circles represent variance-based wave heights $H_{\sigma } /\sqrt 2$ obtained by Hamm 1992)

In Fig. 12, the time-averaged flow velocity field (in which only one out of every four vectors is represented), obtained near the bottom, in accord with the experimental measurements, by the proposed three-dimensional numerical model, is shown. As it can be seen in Fig. 12, the differences in the wave elevation between the plane beach and rip channel drive a current that moves parallel to shore (long-shore current) up to rip channel position where the current turns offshore (rip current). From this figure it is easy to deduce that this circulation pattern represents an erosive condition.

Fig. 12

Plane view detail of the time-averaged velocity field. Only one out of every four vectors is shown

3.4 The effect of submerged breakwaters on nonlinear waves in a curve-shaped coastal area

In this section, we verify if the presence of submerged breakwaters in the curve-shaped coastal area, described in the previous section, is able to modify the wave-induced circulation pattern and the hydrodynamic conditions from erosive to accretive. We chose two submerged breakwaters, separated by a gap, similar to the ones used in the experimental test shown in Sect. 3.2: the barrier length is $3.6\;{\text{m}}$, the distance between the barriers is $1.8\;{\text{m}}$ and the average water depth at the bar crest is $2.67\;{\text{cm}}$. We use the same curvilinear computational grid of the previous test. Figure 13a shows the new bathymetry in the presence of the two submerged breakwaters. Figure 13b, c shows two significant cross sections: one in the gap between the barriers ($y_{A} = 10.5\;{\text{m}}$) and another over the barrier ($y_{B} = 7.5\;{\text{m}}$). For this test, we adopt the same boundary conditions, Smagorinsky coefficient and input wave conditions of the previous test.

Fig. 13

a Bathymetry (only one out of every five coordinate lines is shown). b, c Bottom profiles in sections C and D

A three-dimensional instantaneous wave field is shown in Fig. 14 where the nearshore currents are fully developed. In Fig. 14, it is possible to notice that the breaking waves are induced by the presence of submerged breakwaters because they cause a reduction in the water depth with respect to the one on the plane beach. The wave height which occurs in the gap between the barriers is increased by the presence of an offshore-directed rip current.

Fig. 14

Three-dimensional view of an instantaneous wave field at the time when the breaking-induced circulation is fully developed

The wave height evolution along the two sections, C and D, is shown in Fig. 15a. In Fig. 15b, it is possible to notice the time-averaged cross-shore velocity component, near the bottom, along section C in the gap between the barriers.

Fig. 15

a Wave height comparison between section C (solid line) and section D (dashed line); b Time-averaged cross-shore velocity component along the section in the gap between the barriers (section C). The positive values represent offshore-directed velocities

A plane view of the time-averaged velocity field, near the bottom in accord with the experimental measurements, is shown in Fig. 16. It is possible to notice, from this figure, the presence of a rip current in correspondence with the gap between the barriers (primary circulation) and the presence of convergent currents at the shoreline in the lee of the structure which produce a circulation cell (secondary circulation) with direction opposite to the primary one. This circulation pattern represents hydrodynamic conditions favourable to the accumulation of suspended sediments (accretive conditions). From these considerations, it is possible to deduce that, for these incident wave conditions, the introduction of the submerged breakwaters in this configuration is able to modify the hydrodynamic conditions from erosive to accretive.

Fig. 16

Plane view detail of the time-averaged velocity field. Only one out of every four vectors is shown

Figure 17 shows that the time-averaged cross-shore velocity component, under breaking waves, is offshore directed near the free surface and onshore directed near the bottom (undertow). From this figure, it is possible to notice that the three-dimensional circulation that occurs in the presence of variable bathymetry in a curve-shaped coastal area is well represented by the adopted three-dimensional numerical model.

Fig. 17

Vertical time-averaged cross-shore velocity profile, respectively, at $x = 18.0\,{\text{m}}, \;y = 7.5\,{\text{m,}}$ and $x = 20.0\,{\text{m}},\; y = 7.5\,{\text{m}}$

4 Conclusions

In this paper, a three-dimensional non-hydrostatic shock-capturing numerical model for the simulation of nonlinear wave propagation and nearshore currents over variable bathymetry has been presented. Wave transformation processes, from deep water to the shoreline, including wave breaking and swash zone dynamics have been simulated by integrating the integral contravariant formulation of the three-dimensional Navier–Stokes equations on a time-dependent curvilinear coordinate system. The model has been validated by numerically reproducing experimental tests of breaking waves over submerged breakwater. The good agreement between experimental and numerical results demonstrates that the proposed model is able to reproduce circulation patterns that represent erosive or accretive conditions. The proposed model has been applied to simulate the hydrodynamic effects produced by a system of submerged breakwaters on wave propagation and wave-induced nearshore currents in a curve-shaped coastal area. The numerical results show that, for the simulated incident wave conditions, the introduction of a system of submerged breakwaters is able to modify the hydrodynamic conditions from erosive to accretive.