Bottom changes in coastal areas with complex shorelines

Amodelforthesea-bottomchangesimulationsincoastalareaswithcomplexshorelinesisproposed.Indeepandintermediatewaterdepths,thehydrodynamicquantitiesarecalculatedbynumeri-callyintegratingthecontravariantBoussinesqequations,devoidofChristoffelsymbols.Inthesurfzone,thepropagationofthebreakingwavesissimulatedbythenonlinearshallowwaterequations.Themomentumequationissolvedinsidetheturbulentboundarylayerinordertocalculateintra-wavehydrodynamicquantities.Anintegralformulationforthecontravariantsuspendedsedimentadvection-diffusionequationisproposedandusedforthesea-bottomdynamicsimulations.TheproposedmodelisappliedtotherealcasestudyofPescaraharbor(inItaly).


Introduction
In literature, sea-bottom evolution simulations are generally carried out by adopting two different models: one to calculate the fluid dynamic fields and another for the sediment transport and the dynamics of sea-bottom changes. In this introduction, some considerations are developed, firstly concerning the hydrodynamic models, and secondly concerning the sediment transport models. The depth-averaged equations of motion permit representation of the velocity fields in bidimensional form. In order to reduce the computational effort, the abovementioned two-dimensional models can be seen as an interesting compromise between classical methodologies and recent methodologies based on three-dimensional approaches (Chen, Shi, Hsu, & Kirby, 2014;Keshtpoor, Puleo, Shi, & Ma, 2015;Ma, Chou, & Shi, 2014).
The wave oscillatory features are not explicitly taken into account by the models based on depth-and waveaveraged motion equations (hereinafter called 2DWA). In these models, the effects produced by the waves on the current circulation are taken into account via the use of radiation stresses.
Explicit representation of the periodicity of the wave motion can be obtained using bidimensional models that are not averaged over the wave period (hereinafter called 2DPR). These models use the Boussinesq equations, obtained by defining the depth dependence of the variables, and by depth integrating the equations of motion (Ghadimi, Jabbari, & Reisinezhad, 2012;Karambas & Koutitas, 1992;Ouahsine, Sergent, & Hadji, 2008).
CONTACT F. Gallerano francesco.gallerano@uniroma1.it In these models, there is no decoupling between wave simulations and current simulations. Shock-capturing schemes for numerical integration of the hybrid Boussinesq/nonlinear shallow-water equation models allow explicit simulations of wave breaking (Pu, Shao, Huang, & Hussain, 2013;Roeber & Cheung, 2012;Shi, Kirby, Harris, Geiman, & Grilli, 2012). In fact, these models solve the Boussinesq equations where both the nonlinear and dispersive effects are relevant, and adopt conservative numerical schemes to solve the nonlinear shallow water equations in which nonlinearity prevails over dispersion. The only criterion introduced in the hybrid models is based on the passing of a threshold value in the ratio between the wave height and the water depth, in order to establish when and in which computational cells the solution switches from one set of equations to the other. This ratio is generally assumed to equal 0.8, as suggested by Tonelli and Petti (2010).
The interactions between wave and current and undertow are taken into account by 2DPR models. Moreover, Wenneker, van Dongeren, Lescinski, Roelvink, and Borsboom (2011) emphasized that the wave period variability of the fluid dynamic variables is taken into account by 2DPR models. Deigaard, Justensen, and Fredsøe (1991) underlined that bottom changes in the coastal region are produced by complex hydrodynamic processes: among these, undertow plays a fundamental role in the transport of solid particles in the offshore direction. Fredsøe (1984) highlighted that the hydrodynamic quantities that vary in the wave period drive the resuspension of the solid particles, and their transport and settling. Briganti and Dodd (2009) and Briganti, Dodd, Pokrajak, and O'Donoghue (2011) demonstrated that the swash zone fluid dynamics quantities, especially the duration of swash events, the flow velocity and the maximum uprush, affect morphological bed variations. Such quantities generally come into the net sediment transport rate formulas that represent the swash zone sediment contribution to the surf zone.
Many authors (Kalinske, 1947;Karambas & Koutitas, 2002;Menéndez, Laciana, & García, 2008;Murillo & García-Navarro, 2010) have used the Exner equation for bottom change simulations: in this equation, the bottom time variations are expressed as a function of the long-shore and cross-shore variations of the bed load and suspended sediment load. The Exner equation admits the hypothesis that only the hydrodynamic local conditions drive the suspended sediment load.
In order to dynamically simulate the sediment load, the solid particle concentration equation must be used (these models are indicated hereinafter as 2DHC). When a large portion of the water column is involved via distribution of the suspended sediment load, there is no state of equilibrium for the sediment particle (Nicholson et al., 1997). In the equation adopted in 2DHC models, the spatial variation of the product between the wave-and depth-averaged concentration and wave-and depth-averaged horizontal velocity gives the advective transport. As a consequence, in 2DHC models the transport of solid particles in the offshore direction produced by the undertow is not taken into consideration. In order to simulate the transport of solid particles due to the undertow, vertical profiles of the horizontal velocity and concentration must be assumed not to be uniform. Extant studies (Drønen & Deigaard, 2007;Kim & An, 2011) have simulated the bed evolution dynamics using a quasi-three-dimensional approach (Q3D) in which the advective terms are related to the vertical profile of the horizontal velocity and the vertical profile of the suspended sediment concentration. Other authors (Burger, Kumar, & Ruiz-Baier, 2015;Incelli, Briganti, & Dodd, 2015;Khosronejad, Kang, Borazjani, & Sotiropoulos, 2011;Li & Duffy, 2011) have calculated the morphological changes by means of fully coupled models, which simultaneously solve the nonlinear shallow water and Exner equations. Such full coupling allows one to easily implement different closure laws or to solve an advection equation to represent the suspended sediment transport. Some fully coupled models use a time step for the bed change equation, which is a multiple of the time discretization step of the hydrodynamic model (Lesser, Roelvink, van Kester, & Stelling, 2004;Roelvink, 2006). As shown by Lesser et al. (2004) and Ranasinghe et al. (2011), in these models the time step used for the bed change equation is rarely greater than 20 times the time discretization step of the hydrodynamic model. In this work, we adopt a weak-coupling approach, in which the bed changes are calculated once the hydrodynamic model equations are solved and the hydrodynamic quantities are time averaged. In accordance with the approach of Rakha, Deigaard, and Brøker (1997), the proposed model does not simulate the hydrodynamic phenomena in realworld time; the hydrodynamic phenomena are calculated over a simulation step equal to a multiple of the period of the waves. This approach provides the opportunity to obtain acceptable results, similar to those obtained using fully coupled models, and entails reduced computational efforts.
We present a bottom-change simulation model composed of two submodels: a two-dimensional phaseresolving model by which the fluid dynamic variables changing inside the wave period, the undertow and the swash zone dynamics are calculated; and a second submodel to simulate the bottom changes. in which the suspended sediment concentration is calculated by the wave-averaged advection-diffusion equation with a Q3D methodology. The fluid motion equation and the concentration equation are expressed in a new contravariant formulation. Many authors (Cioffi & Gallerano, 2006;Hu & Shu, 1999;Mandal & Rao, 2011;Sørensen, Schäffer, & Sørensen, 2004;Titarev & Drikakis, 2011) have used unstructured computational grids in order to simulate the fluid dynamic fields on morphologically articulated domains. An alternative strategy is performed by writing the equations of motion in contravariant formulation (Cannata, Lasaponara, & Gallerano, 2015;Gallerano, Cannata, & Lasaponara, 2016a;Luo & Bewley, 2004;Rossmanith, Bale, & LeVeque, 2004). In studies by Shi and Sun (1995), Shi, Dalrymple, Kirby, Chen, and Kennedy (2001), Shi, Kirby, and Hanes (2007), Shi, Kong, andDing (1998), andShi, Svendsen, Kirby, andSmith (2003), the fluid dynamic simulations are carried out by numerically integrating the contravariant Boussinesq equations. In a generalized curvilinear coordinate system, the contravariant components are vectorial components relative to a set of basis vectors that are locally tangent to the coordinate lines. Usually, the Christoffel symbols are present in the contravariant formulation of the equations of motion. These appear as additional terms that arise due to the variability in space of the base vectors and that do not make it possible to express the convective terms in a conservative form. As underlined by Yang, Habchi, and Przekwas (1994), the convergence to weak solutions of conservation laws is not ensured for the numerical methods based on the nonconservative form of the convective terms. In this paper, the contravariant equations of motion are presented in a new formulation characterized by the absence of the Christoffel symbols. This makes it possible to directly simulate the discontinuities related to the breaking waves. The abovementioned momentum equations retain the second-order vertical vorticity term. The propagation of breaking waves is represented by the numerical solution of the nonlinear shallow-water equations: a weighted essentially nonoscillatory (WENO) reconstruction technique is involved in a shock-capturing scheme with high-order. In this paper, we present two main numerical contributions: (1) the proposal of using the Boussinesq equations in a new contravariant formulation, in which dispersive terms are not present in the continuity equation; and (2) the proposal of using an upwind WENO scheme with a high order of accuracy, for reconstructing the conserved variables on the cell faces, that is genuinely two-dimensional. Unlike the commonly used WENO schemes based on a sequence of consecutive one-dimensional reconstructions (Gallerano, Cannata, & Tamburrino, 2012), our proposed scheme is based on two-dimensional interpolating polynomials. Thus, we avoid numerical errors due to the one-dimensional reconstructions, which in generalized curvilinear systems could reach the importance of the terms that arise from Taylor expansions in the Boussinesq equations.
A new methodology for representing the backwash and uprush phenomena is adopted to consider the sediment transport. Calculations of the fluid dynamic variables changing inside the wave period are performed by solving the equations of motion in the boundary layer, according to Deigaard, Fredsøe, and Hedegaard (1986), Fredsøe, Andersen, and Silberg (1985), Rakha (1998) and Rakha et al. (1997). The suspended sediment concentration equation is expressed in contravariant formulation. The Q3D approach is used in order to formulate the advective terms of the suspended sediment concentration equation. The swash zone solid particle production is calculated by following the approach of Larson and Wamsley (2007).
The paper is organized as follows. In Section 2, the hydrodynamic model and the morphodynamic model are described. The numerical scheme is shown in Section 3. In Section 4, validation tests and their application to an engineering case study of the proposed model are presented. The conclusions are drawn in the final section.

The proposed model
The proposed model embraces the two-phase flow one-way coupling assumption, according to which the fluid flow affects the velocity of the solid particles, but not vice versa. Under this assumption, first, the hydrodynamic fields are calculated by a two-dimensional phase-resolving model (hydrodynamic model), and second, the sediment transport and bed morphological changes are calculated by another model (morphodynamic model).

Hydrodynamic model
By indicating with η the elevation of the free surface and h the depth of the local undisturbed water, the total local water depth is defined as H = h + η. Let σ be a distance from the surface of the undisturbed water. The horizontal velocity at z = σ is indicated by u. The horizontal velocity U(z), as proposed by Chen (2006), Nwogu (1993) and Wei, Kirby, Grilli, and Subramanya (1995) is given by where (2) where ∇ = ((∂/∂x), (∂/∂y)) in a Cartesian coordinate system. We indicate with v the depth-averaged value of v(z) correct up to order O(μ 2 )) and O(ε 2 μ 2 ) (with ε = a 0 /h 0 and μ = h 0 /L 0 nondimensional coefficients, in which a 0 is the wave amplitude, L 0 the wave length and h 0 the depth), given by The Boussinesq equations in a two-dimensional Cartesian system are where G is the constant of gravity, V is due to the vertical vorticity approximation, and W and T are the dispersive terms. The bottom resistance term is represented by R. W, T and V are given in Appendix.
The transformation from the curvilinear coordinate system ξ to the Cartesian coordinate system, x, is given by x k = x k (ξ 1 , ξ 2 ) (henceforth the apices indicate components and not powers, with k = 1, 2). The contravariant base vector and the covariant base vector are, respectively, b (k) = ∂ξ k /∂ x and b (k) = ∂ x/∂ξ k . The metric ten- It is known that a generic vector, s, which, with respect to a basis of Cartesian vectors i, j is expressed as s = s x i + s y j, can be expressed with respect to a system of covariant basis vectors, b (1) , b (2) , as s = s 1 b (1) + s 2 b (2) , in which s 1 and s 2 are the contravariant components of vector s given by the expressions s 1 = s · b (1) and s 2 = s · b (2) . By writing the differential operators, which appear in Equations (4) and (5), in the curvilinear coordinate system as and performing the scalar product between Equation (5) and the contravariant base vector b (k) , after simple passages we obtain the following Boussinesq equations, expressed in contravariant formulation, in the curvilinear coordinate system: In which a repeated index in a lower and upper position is held to be summed over its two values, 1 and 2; u k ,v k , V k , W k , T k , R k are the contravariant vectors, which represent, in the curvilinear system of coordinates, u, v, V, W, T and R. The covariant derivative is indicated by an index preceded by a comma. Application of the conservative form of the convective terms of the Boussinesq-type equations is necessary for the definition of a shock-capturing numerical procedure. Gallerano, Cannata, and Villani (2014) wrote the Boussinesq equations in terms of conservative variables H and Hu l . The use of these variables introduce, in the equation of continuity, a new term. A finite difference approximation accurate to the second order was used in order to discretize this new term. In this paper, we present an alternative strategy, according to which H (total local depth) and M k , given by are the conserved variables. With this choice, the continuity equation (Equation (4)) reads The continuity equation (Equation (9)) is written without any source term; consequently, a high-order shock-capturing finite-volume scheme can be used to solve this equation. The contravariant formulation of the momentum equation, in which M k is the conserved variable, is where D k is the '∧' refers to the exponent of the power; η c is a constant (according to Xing & Shu, 2006) and the term W k has been split into where W k and W k are given in Appendix.
In the momentum equation (Equation (10)), covariant derivatives are present. These derivatives produce the Christoffel symbols (given by Equation (A.3)). The Christoffel symbols come into in the momentum equation as source terms; consequently, it is not possible to write the convective terms in conservative form. In order to overcome this drawback, we propose the Boussinesq equations in contravariant formulation without the Christoffel symbols. We equate the net force to the material volume momentum rate of change in a direction defined by a parallel vector field (Gallerano & Cannata, 2011). Let b (k) (ξ 1 0 , ξ 2 0 ) be the contravariant base vector at point P 0 (ξ 1 0 , ξ 2 0 ) ∈ A. In Figure 1 a graphic sketch of the generic surface element A in the curvilinear coordinate system is shown. In this figure, a pair of covariant base vectors b (1) , b (2) and a pair of contravariant base vectors b (1) , b (2) defined at the generic point P(ξ 1 , ξ 2 ) are shown; We denote b (k) (ξ 1 0 , ξ 2 0 ) by b (k) , and b (l) (ξ 1 , ξ 2 ) by b (l) . The integration of Equation (10) over the surface element of area A and the projection of this equation along λ l gives ∂D l ∂t dA Figure 1. Graphical sketch of a computational cell in the curvilinear system of coordinates ξ 1 , ξ 2 . Gray colored area: surface element A bounded by four segments lying on coordinate lines of the curvilinear coordinate system. g (1) , g (2) : covariant base vectors defined at the generic point P(ξ 1 , ξ 2 ). g (1) , g (2) : contravariant base vectors defined at the generic point The continuity equation becomes where n r is the covariant vector, which is normal for the contour line L of A. The accuracy of Equations (14) and (15) is of the order of O(μ 2 ) and O(ε 3 μ 2 ); furthermore, Equation (14) retains the second-order vertical vorticity. The use of H and M l , as conserved variables, produces a momentum equation that is different from the equation in Gallerano et al. (2014). This difference is related to the different form of the convective terms. Furthermore, the continuity equation (15) can be solved by a numerical scheme based on a shock-capturing methodology because there are no dispersive terms. Equations (14) and (15) are numerically integrated by a high-resolution hybrid finite volume-finite difference scheme, as shown in Section 3. The shock-capturing method proposed in this work makes it possible to intrinsically model the wave breaking; therefore, no additional terms are needed to account for the energy dissipation due to the wave breaking in the surf zone.

Undertow
Indicating by u k B (z) the corrective contravariant velocity vector (Lynett, 2006), the horizontal velocity u k (z) reads as follows where u k α is the horizontal velocity contravariant vector computed by Equations (14) and (15).

Fluid dynamic quantities
We indicate by U(z, t) the Cartesian horizontal component of the fluid velocity; (t) indicates the thickness of the boundary layer, U (t) the horizontal velocity at the top of the wave boundary layer and u f (t) the friction velocity. The momentum equation integration inside the boundary layer gives where the lower limit in the integral, k/30, represents the characteristic length scale experimentally evaluated by Nikuradse (1933), in which k is the bed roughness that, according to Fredsøe et al. (1985) and Engelund and Fredsøe (1976), is assumed equal to 2.5d 50 , with d 50 the sediment mean diameter. By adopting the approach proposed by Fredsøe (1984) and also used by other authors (Briganti et al., 2011;Fredsøe & Deigaard, 1993;Williams, Briganti, & Pullen, 2014), the vertical profile of the horizontal velocity in the wave boundary layer is assumed to be logarithmic With the introduction of Equation (18) into Equation (17), the variation in time of the friction velocity u f (t) results from the following equation: in which U is the horizontal velocity at top of the boundary layer.
Denoting by f (t) the following quantity, and using the logarithmic velocity profile given by Equation (18), the following expression for the thickness of the boundary layer is obtained: Solving the system composed by Equations (19), (20) and (21), the values f (t), u f (t) and (t) are given. The average of the instantaneous quantities over the wave period T is indicated by the mark []. The turbulence inside the boundary layer is produced by the interaction between wave and current. Let be u fc the current friction velocity given by Within the boundary layer, the eddy viscosity is

outside the boundary layer the eddy viscosity is
Under breaking waves, the turbulence is given by the contributions produced by current, wave boundary layer and wave breaking. The turbulent kinetic energy equation (Deigaard et al., 1986) comes into in the calculation of eddy viscosity ν t,f (z, t) related to the breaking of the wave: where l is the turbulence length scale, k t = k t (z, t) is the kinetic energy of the wave breaking-induced turbulence and c d = 0.08. The kinetic energy production is where CSM(t) is a coefficient (Kennedy, Chen, Kirby, & Dalrymple, 2000), β p is the portion of the period of the wave in which the kinetic energy of turbulence is produced, δ = −z is the distance from the free surface,H w is the height of the wave and E loss is the kinetic energy dissipation.
The integration of Equation (25) gives the instantaneous value of k t which is used in order to calculate the eddy viscosity produced by the wave breaking as Consequently, the total eddy viscosity, ν t (z, t), is the quadratic sum of the eddy viscosity due to the current and wave breaking and the eddy viscosity produced by the wave boundary layer: The solution of the system of equations composed by Equation (14) (momentum balance equation) and Equation (15) (continuity equation) gives u k α and the water depth H. From these values, from Equation (16) we obtain the modified velocities u k (z, t), as proposed by Lynett (2006). The solution of the system of Equations (19), (20) and (21) gives the instantaneous values of the friction velocity u f , the thickness of the boundary layer wave and the nondimensional quantity f . Once u f is known, from Equation (22) we obtain the friction velocity due to the current u fc , from which it is possible to calculate, from Equation (23), the vertical eddy viscosity distribution ν t,r (z, t)within the wave boundary layer and, from Equation (24), the vertical eddy viscosity distribution outside the wave boundary layer. From Equation (26), we calculate the turbulent kinetic energy production, PRODU, under breaking waves. Once PRODU is known, Equation (25) gives the instantaneous values of the turbulent kinetic energy, k t , under breaking waves, by means of which we obtain the contribution to the eddy viscosity vertical distribution due to the wave breaking, ν t,f (z, t) (from Equation (27)). From Equation (28), the total eddy viscosity ν t (z, t) is calculated. The instantaneous values of u k , u f , ν t , H and are used to calculate the input variables of the morphodynamic model, as shown in the following subsection.

Suspended load
The differential form of the solid particle concentration equation, in a Cartesian formulation, in which a Q3D methodology is used, reads in whichC(z) is the solid particle concentration averaged on the wave period,H is the water depth,ũ x (z) andũ y (z) are the wave averaged vertical distributions of the components of the horizontal velocity vector u(z) defined in a Cartesian system of reference. The distance, a, from the bottom defines the region in which the bed load transport develops; ν t is the depth and wave-averaged eddy viscosity; D is the rate of the sediment deposition and P is the rate of turbulent sediment pick-up. It must be noted that the Q3D approach is used in Equation (29): on the lefthand side of this equation, the second and the third terms are calculated by the vertical distributions of velocity and concentrations. By writing in contravariant form (in the curvilinear system of coordinates) the velocity vector ũ that appears in Equation (29), and by integrating this equation over a generic surface element A (whose contour line is L), we obtain the following contravariant concentration equation: in which n is the outward normal vector to the contour line L andũ r is the horizontal contravariant velocity vector averaged over the wave period. D and P are expressed by in which w sed is the sediment fall velocity,C a andC R are, respectively, the actual and reference concentrations, which are evaluated at height a = 2d 50 . A threshold value of the sediment particle motion comes into the calculation ofC R . In order to integrate Equation (30), the calculation ofC a andC R is needed and shown hereinafter. The value ofC a is taken as the lower boundary condition of the turbulent suspended sediment diffusion equation, and as the lower extreme of the integral that gives the depth-averaged value ofC(z).
Thus,C a is calculated by an iterative procedure using Equation (34), where the values ofC and ν t (z) are known (from the previous time step). The value ofC R is obtained by wave-averaging its instantaneous values C R (t), which are calculated as (Rakha et al., 1997) in which C m is the maximum volumetric concentration that can be reached, θ cr is the parameter of stability of Shield and | θ| = | θ(t)| is the parameter of mobility of Shield, where θ(t) is the bed shear stress induced by current and wave. It must be noted that the instantaneous value of the reference concentration (C R (t)) results in a nonnull value only if the bed shear stress goes over the θ cr threshold value. The turbulence induced by the breakers generally occurs in the upper part of the water column, so it should not influence the bottom concentration significantly (Nielsen, 1992). Furthermore, the turbulence generated in the water surface by the wave breaking results in more sediment being carried in suspension away from the bed, while the sediment concentration near the bed is still determined by the wave boundary layer (Deigaard et al., 1986).

Bed load
In a Cartesian system of coordinates, the vector q b of the bed load transport is given by (Engelund & Fredsøe, 1976) as: in which β is the coefficient of the dynamic friction and ρ s /ρ w is the ratio between the sediment and water density.
The numerical solution of the concentration equation admits an uprush boundary (see Figure 2) condition given by the swash zone transport of sediment (Larson & Wamsley, 2007;Nam, Larson, Hanson, & Hoan, 2009). Along the uprush boundary, we define a system of Cartesian coordinates in which the line of maximum slope (offshore directed) identifies the x * axis. We indicate by q net the wave-averaged net sediment transport components, as given by Larson and Wamsley (2007) in which m is the moving grain friction angle, K x * and K y * are coefficients, dh/dx and β e are, respectively, the foreshore and the foreshore equilibrium slopes, and u 0 and v 0 are the scalar components in the local Cartesian reference system of vector u 0 , which represents the speed of the wet and dry front.

Bed morphological changes
The total transport is given by the sum of the bed load transport, which takes into account the near bed transport mechanism, and the suspended load transport. The equation of the bed change expressed in a Cartesian reference form is in which p is the porosity of the sediment and z f is the elevation of the bed. The contravariant components of the vector q b areq k b (k = 1, 2), which are obtained by dotting the vector q b with the contravariant base vector b (k) . The final contravariant formulation of Equation (39) in a curvilinear coordinate system is Let us indicate by q S the net sediment transport averaged over T, which represents the entire swash cycle duration, and let t 0 be the swash cycle duration at point x * inside the uprush and backwash zone. Let u 0S be the speed vector (whose components are u 0S , v 0S ) of the wet and dry front calculated by solving a Riemann problem. The components of q S are The bed changes in the swash zone are given by where q k S is the contravariant component of q S . From the friction velocity u f , using Equation (35) we calculate the instantaneous value of the reference concentration C R (t) and, from this, its averaged value over the wave periodC R . Equation (32) gives the sediment pickup rate, P, which appears in Equation (30). From the eddy viscosity ν t (z), by solving the system of Equations (33) and (34), we obtain the wave-averaged concentration of suspended sedimentC(z) andC a . Starting from the water depthH, the velocityũ r (z), the total eddy viscosity ν t (z) and the source terms P, we calculate the wave and depth averaged concentration of suspended sedimentC and the sediment deposition rate, D, by solving Equation (30). From the instantaneous value of the horizontal velocity u k α , using Equation (35) we calculate the instantaneous value of the bed load transport vector, q b , and, from this, its wave-averaged value q b . The values of the terms P, D and q b are inserted into Equation (40) in order to update the bed elevation in the region from the deep-water zone to the surf zone (zones A-B-C-D in Figure 2).
In the swash zone, from the instantaneous values of wave propagation velocity vector u 0S , we calculate the net sediment transport vector q S , which is averaged over several swash cycles. This averaged value is introduced in Equation (41), in order to update the bed elevation in the swash zone (zone F in Figure 2).
The procedure by which the two models are coupled is as follows. A morphological step is equivalent to fulfillment of the following five-step sequence.
(1) The vertical distribution of the horizontal velocities, the free-surface elevation, the friction velocity, the wave boundary layer, the eddy viscosity and the reference concentration are calculated by the hydrodynamic model.
(2) The instantaneous hydrodynamic quantities are averaged over T * , which is 100-150 times the period of the wave (Rakha et al., 1997).

(3) The equation for the concentration of solid particles
is integrated by using the wave averaged hydrodynamic quantities calculated in step (2). (4)C a is computed starting from the values ofC, calculated in step (3).
The bed change equation uses the values of the reference concentration,C R , actual concentration,C a , and bed load transport, q b , and gives the updating of the sea bottom for the new morphological time step. We define as 'morphological time step' the time step of the integration of the bed change equation. This time step is greater than the Boussinesq model time step integration. The morphological time step is obtained via a posteriori verification of the numerical results. The morphological time step is chosen based on a trial-and-error procedure. The abovementioned procedure is explained by the flow chart shown in Figure 3. The value of the morphological time step is the result of a compromise. Values of the morphological time step that are too small entail excessive computational time; instead, values of the morphological time step that are too big cause changes in the bathymetry that are too sudden and do not establish a realistic bed evolution dynamics.

Numerical scheme
Equations (14) and (15) are numerically integrated by an original scheme with a fourth order of accuracy: the shock-capturing methodology is based on genuinely two-dimensional upwind-WENO reconstructions. Let L(M 1 , M 2 ) be the second side of Equation (15). Let D(H, M 1 , M 2 ) be the sum of the first four terms of the second side of Equation (14) and let D B (H, M 1 , M 2 ) be the sum of the remaining terms on the second side of Equation (14). The integration over [t n , t n+1 ] of Equations (14) and (15) reads i;j averaged over the surface element of area A.
The discretization in time of Equations (44) and (45) is performed by means of a Runge-Kutta procedure belonging to the Strong Stability Preserving family (SSPRK) (Spiteri & Ruuth, 2002), which can be written as in which p = 1, 2, 3 and the values of pq and ϕ pq are found in Spiteri and Ruuth (2002). At every step of the Runge-Kutta method, the velocity component,ũ k , is obtained by numerically solving the following equation: The values of the velocity component at the cell centroids,ũ k i;j , are introduced in the second-order finitedifference discretization of D B (H, M 1 , M 2 ). The remaining terms on the second side of Equation (14) and the term on the second side of Equation (15) are discretized by a fourth-order genuinely two-dimensional upwind-WENO scheme. The procedure at the basis of the numerical scheme can be summarized as follows: • A WENO procedure, based on genuinely twodimensional polynomials, reconstructs the point values of the variables at every intercell face of the computational cells.  (47) and (48).

2D WENO reconstructions
The WENO procedures commonly adopted in literature for reconstructing the variable point values on the cell faces are not genuinely multidimensional because they are carried out by a sequence of consecutive one-dimensional reconstructions (dimensionby-dimension reconstruction). In the dimension-bydimension approach, a first one-dimensional reconstruction is carried out along the coordinate line ξ 1 in order to pass from cell-averaged values to values that are averaged only over the coordinate line ξ 2 ; a second one-dimensional reconstruction, performed along the ξ 2 coordinate line, is then used to obtain the point values of the conserved variables on the cell face. In generalized curvilinear systems, the above sequence of consecutive one-dimensional reconstructions can produce an amplification of the truncation errors, such that the resulting errors could reach the importance of the terms that arise from Taylor expansions in the Boussinesq equations, thereby compromising the accuracy of the numerical solution. In this paper, in order to avoid this problem, we adopt a WENO technique based on genuinely two-dimensional polynomials for the point values calculation of H and M k .
Starting fromH i;j andM k i;j , (which represent the cellaveraged values), the point values of the conserved variables (H i+p;j+q and M k i+p;j+q with p = −1/2, 1/2 and q = −1/2, 1/2) are reconstructed by using a combination of nine two-dimensional biquadratic interpolating polynomials P i+p;j+q (ξ 1 , ξ 2 ) (with p = −1, 0, 1 and q = −1, 0, 1) given by For each two-dimensional polynomial, the nine coefficients b m (with m = 0, . . . , 8) are calculated by simultaneously using the values of the conserved variables in a two-dimensional stencil composed of nine (3×3) computational cells. In order to obtain a fourth-order reconstruction that is devoid of spurious oscillations, H and M k are approximated by a convex combination R i;j (ξ 1 , ξ 2 ) of the nine two-dimensional polynomials defined on a wider two-dimensional stencil composed of 25 (5×5) computational cells: In such a way, the variable point values on each face of the computational cell I i;j are reconstructed by simultaneously taking into account all the 25 cell-averaged values that surround and include the I i;j cell. In Equation (52), ω k,l i,j (where k and l refer respectively to the ξ 1 and ξ 2 coordinates) are the nonlinear weights of the genuinely two-dimensional WENO reconstruction. Such nonlinear weights are a function of the linear weights and indexes of smoothness (calculated in order to avoid spurious oscillations in the presence of discontinuities), according to Gallerano, Cannata, and Lasaponara (2016b) and Levy, Puppo, and Russo (2002).

CRIEPI large wave flume experiments -Test 1-8
In this section, an experimental test case, extracted from the Central Research Institute of Electric Power Industry (CRIEPI) Test proposed by Kajima, Shimizu, Maruyama, and Saito (1982), is numerically reproduced in order to verify the proposed sediment transport model under erosional regular waves. This experimental test is referred to by the same authors as CRIEPI Test 1-8, and reproduces erosive conditions, over a time interval of 21 h, under a regular wave of a height of 0.85 m and a period of 3.0 s. The initial profile is formed as a uniform slope of 1:20 and is composed of coarse sand with a mean diameter of 0.47 mm. For the hydrodynamic model the time step is 0.005 s; the spatial discretization step is 0.1 m. A morphological time step of 1.0 hours is adopted in order to simulate the bottom change. In Figure 4, the comparison between bed elevation change, with respect to the initial bottom profile (dashed line), for the CRIEPI Test 1-8 obtained by the numerical simulation (black line) and experimental data (red line) is shown. As can be noted from the figure, the proposed sediment transport model is able to consistently predict the magnitude and extension of the bed erosion and bed accretion areas.
In particular, approaching the shoreline, two different areas can be seen from the numerical results: a first area characterized by a bottom accretion around X = −22 m, which shows a 0.35 m height bar, and a second area that is more extended than the first and is characterized by a 0.25 m-deep erosion located around X = −14 m.

Test T1C1
In this section, the proposed model is validated by reproducing the T1C1 experimental test of Gravens and Wang (2007). These authors carried out the T1C1 test in a laboratory basin in which a breakwater is positioned 4 m away from the shoreline. The length of the breakwater is 4 m and its cross-section can be sketched by a trapezoid whose dimensions are lower base 1.5 m, upper base 0.7 m and high 0.25 m. The experimental test was carried out by generating random waves of significant wave height, H s = 0.26 m; wave period, T = 1.5s; and incoming angle of θ = 6.5°, and with sediment particles with mean diameter d 50 = 150 μm. Figure 5 shows the computational domain (ABCD area) adopted for this test, in which the area of the laboratory basin is schematized by a dashed red rectangle (EFGH area). A grid spacing of x = 0.06 m and y = 0.045 m is used. Simulation of the T1C1 test is performed by generating, along the BC line, random waves characterized by a JONSWAP spectrum with H s = 0.26 m. At the side boundaries (AB and CD lines), a null gradient in the normal direction is imposed for the velocity and the free surface elevation, since the wave fronts are normal to these boundaries during the simulation. Throughout the numerical simulation, in the E-H stretch of coast, a longshore current occurs with velocity values that are practically uniform along the paths. This current velocity field comes very close to that produced in the laboratory during the T1C1 test. Concerning the morphodynamic model, side boundary conditions are defined by a null gradient of concentration and bed load transport.
With reference to the T1C1 test, Figure 6 shows the contour lines of the water depth at the starting conditions and the traces of the two sections where the laboratory measurements are taken. The time step chosen for the two-dimensional phase-resolving model is 0.0058s. The bed change simulation is carried out in 190 morphological steps. The bed elevation variation in every morphological step approximates the one occurring in a real time interval of 60s (morphological time step). Figure 7 shows the simulated wave field for the T1C1 test. In Figure 8, the numerically obtained significant wave height is compared with the experimental results: the black line represents the results obtained with the proposed numerical model, the circles represent the experimental data and the green line represents the initial bottom profile. Figure 8(a) and (b), respectively, refer  to sections Y = 26 m and Y = 22 m, which intersect the extremes of the breakwater. From these figures it can be seen that the breakwater protects the onshore area from the incoming waves, as demonstrated by the  relevant decrease in the wave height onshore of the breakwater. Figure 9 shows the current velocity field, referring to the initial condition of the bottom, obtained by averaging over time the simulated instantaneous velocities. In this figure, the vectors representing the current velocity values experimentally obtained by Gravens and Wang (2007) appear in bold. As shown in the figure, the simulated longshore current is partially intercepted and diverted offshore by the breakwater. The presence of the breakwater also induces a more complex structure of the wave-averaged velocity field in the area between the breakwater and the shoreline. It is in fact known that, in this area, two vortices are formed in proximity of the breakwater: one that rotates anticlockwise with the center placed rear of the breakwater at about Y = 25.5 m and X = 2.5 m, and a second that rotates clockwise, larger than the first, whose center is closer to the coastline (Y = 20.5 m and X = 1.5 m). In the area rear of the breakwater, the longshore current is confined between those two vortex structures. Onshore of the breakwater, the longshore current first undergoes (in the current direction) an increase of velocity in the region closest to the shoreline, and is then diverted, due to the presence of the second vortex, toward the breakwater; from Y = 22 m, the current reverts parallel to the shoreline. As underlined by Nam et al. (2009), the presence of the two vortices, in a position not symmetrical with respect to the center of the breakwater, is a typical phenomenon of current velocity fields produced onshore of the breakwater by waves obliquely incident to the shoreline. In Figure 9, the vectors in bold, representing the experimentally measured current velocities, reveal two zones onshore of the breakwater where the current velocities are close to zero: a first current stagnation area is located close to the breakwater and is approximately centered at Y = 25.5 m and X = 2.5 m; a second current stagnation area is closer to the shoreline and is approximately centered in Y = 20.5 m and X = 1.5 m. The centers of the two vortices produced by numerical simulation described above (in which the current velocity is zero) correspond with the center of the current stagnation areas experimentally measured.
In Figure 10, the longshore component of the simulated velocity field is compared with the experimental results. The numerically simulated current velocity is slightly underestimated in the region between X = 12 m and X = 7 m. From Figure 10(a), it can be observed that onshore of the breakwater (2 m < X < 4 m) the numerically calculated current velocities are slightly negative where the experimental data indicate values almost equal to zero. Close to the shoreline (X < 2 m), in comparison to the laboratory measurements, the simulated longshore current is slightly overrated. Small discrepancies between numerical and experimental results can also be observed in Figure 10(b) for X < 3 m, where the simulated longshore current is slightly overrated in the area closest to the breakwater (X > 2 m) and changes direction near the shoreline (X < 1 m). Figure 11 compares the simulated cross-shore component of the velocity current with the experimental data.  Figure 10. T1C1 test. Comparison between the longshore current calculated by the proposed model (black lines) and the experimental data (red circles). Bottom profile (green line).
(a) (b) Figure 11. T1C1 test. Comparison between the cross-shore current calculated by the proposed model (black lines) and the experimental data (red circles). Bottom profile (green line).
Experimental and numerical results show good agreement. Some differences can be observed onshore of the breakwater where, as mentioned before, the experimental data show a current stagnation area while the numerical simulations produce a small vortex. In order to assess a quantitative evaluation of the reliability of the proposed model, an estimate is performed of the percentage error by which the numerical results (in terms of significant wave height, cross-shore velocity and long-shore velocity) differ from those experimentally obtained in the T1C1 test carried out by Gravens and Wang (2007). For each of the above quantities, the Mean Absolute Percentage Error (MAPE) is calculated as In which N represents the number of measurement locations where the experimental data are provided, and φ i , η i are, respectively, the experimental and numerical values of the generic quantity in the i-th measurement point.
In Table 1, the computed errors in terms of significant wave height, and cross-shore and long-shore velocity are reported. Table 1 shows that the error in terms of significant wave height is 14%, whilst the errors in terms of cross-shore and long-shore velocity are less than 40%. The higher discrepancies in the experimental results observed for the cross-shore and long-shore velocities can be explained by taking into account that the longshore current produced in the experiment is stronger than the one numerically generated. In fact, in the experiment the long-shore current is enforced by means of a pumping and recirculation system, whilst in the simulation it is exclusively generated by the obliquely breaking waves. Contour lines of the water depth after 190 simulated minutes, for the T1C1 test, are shown in Figure 12. Grey lines represent the contour line of the water depth experimentally produced by Gravens and Wang (2007), blue lines represent those obtained by the proposed model and magenta lines represent those obtained with the model proposed by Nam, Larson, Hanson, and Hoan (2011). By comparing Figure 12 (blue lines) and Figure 6, the bed morphological changes produced by the proposed model can be seen. Onshore of the breakwater, between Y = 27 m and Y = 25 m, an advance toward the breakwater of the 0.1 m depth contour line indicates the presence of a bed accretion area: this accumulation is caused by the sediment put into suspension in the area updrift the breakwater near the shoreline that is prone to precipitating by the reduction of the velocity values related to the first vortex. Two troughs can also be seen near the breakwater corners. Figure 12 also shows that the results obtained from Nam et al.'s (2011) model and the experimental data are in quite close agreement, and that the main differences correspond with the advancing shoreline onshore of the breakwater and near the top corners of the same. In these areas, where high sediment deposition (onshore of the breakwater) and intense sediment troughing (near the corners of the breakwater) take place, the depth contour lines obtained by the proposed model better match the experimental data compared to the model of Nam et al. (2011). This better agreement can be justified by the fact that the model of Nam et al. (2011) introduces some simplifications with respect to the proposed model. In fact, this model is composed of a 2DWA model for current velocity simulations (based on the use of radiation stresses), by a sediment transport model which solves the concentration equation that does not adopt a Q3D approach and assumes, for the eddy viscosity, a uniform distribution along the depth.

Long-term bottom changes
In this subsection, the proposed model is applied to the real case study of Pescara harbor (in Italy). In 1997, a detached breakwater was constructed in front of Pescara harbor. In Figure 13, the breakwater is shown and the coastal structures and shoreline are sketched, respectively, by red lines and a dashed black line. In 1997 and 2000, two measurement campaigns were performed by  Figure 14(a) and in Figure 14(b), contour lines of the water depth are shown. In the area between the breakwater and the canal port, the sediment volume (accumulated from 1997 to 2000 and estimated by the two measurement campaigns) is about 40,000 m 3 /year.

MELPS (2008b): in
The proposed model is tested via numerical simulation of the bottom changes in front of Pescara harbor in Italy related to the three years following construction of the detached breakwater. The numerical simulations were performed using, as input data for wave motion forces, the available data on occurrence frequency, period and wave height in the region of Pescara harbor (MEPLS, 2008a). Such data are related to a measurement point located about 500 m northwards of the detached breakwater and are considered, in the abovementioned weather and sea study, as representative of wave climate in the sea region in front of Pescara harbor. Figure 15 shows the polar representation of the direction of the incoming waves representative of Pescara's annual wave climate: two preeminent directions are identified as 345°N-15°N (primary sector) and 65°N-95°N (secondary sector). According to the MEPLS (2008b) study, the sediment transport produced by the secondary-sector incoming waves is almost completely intercepted by the touristic port (Figure 13), so it does not significantly influence the sedimentation phenomena close to Pescara harbor. Furthermore, according to the MEPLS (2008b), the sediment transport produced by the primary-sector incoming waves is bound to deposit onshore of the breakwater, causing silting-up phenomena in the area that stands between the canal port and the breakwater. Such incoming waves are marked by a significant wave height in the range of 1 m and 2 m, which have an occurrence frequency no fewer than 320hours/year. Consequently, in order to achieve results that are representative of Pescara's annual wave climate, the bottom change simulation is performed by generating random waves incoming from 0°N, represented by a spectrum belonging to the JONSWAP type characterized by a significant wave Figure 15. Pescara harbor case study. Directional annual distribution polar histogram of the wave events in the coastal region in front of Pescara harbor.
height equal to 1.5 m, which act for 320 h/year for a total of 960 simulated hours. The equations of motion (Section 2) are discretized on a curvilinear grid whose boundaries have been established by applying a shoaling-refraction model (Cialone & Kraus, 1987) on a wider coastal region with respect to the one used in the hydrodynamic model. The obtained wave rays define the northwest and southeast boundaries of the computational grid, and the deep Figure 16. Pescara harbor case study. Instantaneous wave field two-dimensional representation.
water wave front defines the northeast boundary. The beach is taken into consideration by means of a wet and dry procedure, while over the structure borders reflective conditions are adopted. The time step used for the hydrodynamic model is t = 0.1s. From Figure 16, diffraction phenomena occurring in the proximity of the detached breakwater edges, shoaling phenomena and then wave breaking with a wave height decay can be seen. Furthermore, corresponding to the vertical walls protecting the touristic port, reflection phenomena are evident. Figure 17 shows the fluid dynamic field produced by the proposed model. In the same figure, the three blue circles indicate the points at which a comparison between the simulated and measured velocities is carried Figure 17. Pescara harbor case study. Current velocity field obtained by the proposed model (black vectors) and measurement points (blue circles).  (Table 2). This comparison indicates that the computed longshore current is in good agreement with the measured one. Figure 18 shows the comparison between the contour lines of the water depth produced by the numerical model after three years (blue lines) and the corresponding measured data (black lines). The agreement is good in most of the computational domain. Some differences between numerical results and measurement data are to be seen in the lower part of the computational domain, near the northwest extreme of the detached breakwater. In this part of the computational domain, experimental measurements indicate the permanence of a restricted area characterized by the 6 m-depth contour line. The numerical results indicate a slight overestimation of the silting up in this area and a slight underestimation of the silting up in the part of the computational domain right above the aforementioned area. In order to provide a quantitative evaluation of the reliability of the proposed model for long-term prediction of bed morphology evolution, an estimate is performed of the percentage error by which the numerical results (in terms of water depth) differ Figure 18. Pescara harbor case study. Comparison between the numerically calculated depth contour lines at the end of the third simulated year (blue lines) and the corresponding measured data (black lines). from those obtained by the measurement campaign. The water depth MAPE is calculated by means of Equation (53) and its value is found to be about 10%. Furthermore, the numerically calculated volume of sediment settled in front of Pescara harbor (from 1997 to 2000) is about 37,000 m 3 /year, which is comparable with the measured one (40,000 m 3 /year). These results show that the numerical simulation suitably reproduces the bottom changes in the coastal region of Pescara.
The morphological time step adopted for the Pescara case study was chosen following comparison between the results of different simulations carried out with various values of the morphological time step, ranging between 0.5 and 12 hours. It turns out that when the morphological time step value is less than 4 hours, no significant differences in the bed changes are observed at the end of the 120 simulated hours. For morphological time step values included between 4 and 8 hours, the differences among the numerical results are less than 6%. These differences become significant and rapidly diverge for morphological time step values greater than 10 h.

Conclusions
A model for the bottom change simulations in coastal areas with complex shorelines has been proposed. The contravariant equations of motion have been presented in a new formulation characterized by the absence of the Christoffel symbols, in order to minimize truncation errors that arise in the calculation of the hydrodynamic variables on generalized curvilinear grids. The sediment transport and bed changes have been calculated by a suspended sediment advection-diffusion equation expressed in an integral contravariant form and integrated by a Q3D approach. The proposed model has been validated against experimental tests and by applying it to the real case study of Pescara harbor (in Italy). The close accord between the measured data and the numerical results indicates the capacity of this model to simulate long-term bed changes in sea regions characterized by complex-shaped coastlines.

Disclosure statement
No potential conflict of interest was reported by the author(s). W =