Numerical Study over the E ﬀ ects of a Designed Submerged Breakwater on the Coastal Sediment Transport in the Pescara Harbour (Italy)

: In 1997, in front of the Pescara Harbour (Italy), a detached breakwater was constructed. In the successive years, the sediment transport due to the combined action of waves and coastal currents, in the area between the detached breakwater and the entrance of the Pescara Harbour, produced an accumulation of about 40,000 m 3 of sediment per year. In this paper, the causes of the accretion of the bottom elevation in front of the Pescara Harbour entrance and the e ﬀ ects produced by the existing detached breakwater are investigated. The e ﬀ ects on the sediment transport of the introduction of a new submerged breakwater designed to protect the entrance of the harbour from sediment siltation are investigated. In particular, the ability of the designed submerged breakwater, located orthogonally to the longshore current, to intercept the aforementioned solid material and to signiﬁcantly reduce the accretion of the bottom in the area in front of the harbour entrance, was numerically veriﬁed. Numerical simulations were carried out by means of a model of the bottom-change composed of two sub-models: a two-dimensional phase resolving model that is used to calculate the ﬂuid dynamic variables changing inside the wave period and a second sediment transport sub-model to simulate the bottom changes, in which the suspended sediment concentration is calculated by the wave-averaged advection–di ﬀ usion equation. The equations of motion, in which the vector and tensor quantities are expressed in Cartesian components, are written in a generalised curvilinear coordinate system. The fully nonlinear Boussinesq equations are written in an integral form and used to simulate the velocity ﬁelds.


Introduction
In the literature [1,2], the simulation of sea bottom changes produced by a coastal defence structure are usually carried out by simulating the hydrodynamic fields and the concentration of the suspended solid particles. In the framework of coastal sediment transport, 3D models for the simulation of the hydrodynamic fields [3][4][5][6] require considerable computational time, due to the long-time scale of the sea bottom variations. For this reason, in the context of morphodynamic simulations, the hydrodynamic equations [7][8][9][10][11][12][13] and the suspended sediment concentration equation are usually depth averaged [14].
In this context, the two-dimensional phase-resolving models, in which motion equations are not averaged over the wave period, employ the so-called Boussinesq equations, obtained with the depth integration of the motion equations and with the definition of the depth dependence of the variables. Shock-capturing schemes adopted for numerical integration of the fully nonlinear Boussinesq equations (FNBE) [15][16][17][18] allow explicit simulations of the wave propagation from deep water to the shoreline, including the wave breaking.
In the coastal regions, sea bottom modifications are caused by complex hydrodynamic processes; in these processes, the variation of the hydrodynamic quantities over the wave period causes the resuspension, transport and settling of the solid particles. The variability of the hydrodynamic field and the wave-current interactions over the wave period are taken into account in FNBE models.
The undertow is one of the phenomena that play key roles in the solid particles transport. In the undertow phenomenon, the horizontal velocities that are closest to the bottom are off-shore directed in the surf zone, while the horizontal velocities that are closest to the free-surface are on-shore directed. Lynett [19] proposed a correction to the vertical distribution of the horizontal velocity to simulate the undertow in the context of models that solve depth-averaged motion equations.
In the standard forms of FNBE, in the continuity equation, dispersive terms are present. In particular, in several works in the literature [20][21][22], located in the context of hybrid (finite volumefinite difference) numerical schemes, the dispersive term that is present in the continuity equation is discretised by means of a second-order cell-centred finite difference scheme. As demonstrated in [15], in the continuity equation, the truncation errors related to the discretisation of the dispersive term may reduce the accuracy of the numerical scheme and cause oscillations in the numerical solution. The dispersive terms in the continuity equation can be avoided by written FNBE in a proper conservative form.
In this paper, we present a simulation model for the sea bottom modifications, composed by two sub-models: the hydrodynamic modifications are simulated by means of a two-dimensional phase-resolving model that considers the variation of the variables inside the wave period and takes into account the undertow and the sea bottom modifications are simulated by means of a second sub-model that calculates the suspended sediment concentration by means of the wave-averaged advection-diffusion equation.
The motion equations and the suspended sediment concentration equations, in which the vector and tensor quantities are expressed in Cartesian components, are written in a boundary conforming generalised curvilinear coordinate system. In the region extending from deep water up to the beginning of the surf zone, the wave and hydrodynamic fields are calculated by an integral form of the FNBE; in the surf and swash region, the governing equations reduce to the nonlinear shallow water equations, in order to directly simulate the wave breaking. In the proposed model, in which the continuity equation is written in an integral form, dispersive terms are avoided. Therefore, the continuity equation is entirely discretised by a shock capturing finite volume scheme and the errors caused by the low-order discretisation of the dispersive terms are avoided.
The near-bed velocity, boundary layer thickness, friction velocity and bed shear stress are evaluated from the integration of the momentum equation over the wave boundary layer. The turbulence contribution of the wave boundary layer is considered for the evaluation of the eddy viscosity vertical distribution under breaking waves.
Following Gallerano et al. [16], the hydrodynamic numerical scheme employs genuinely 2D WENO (Weighted Essentially Non-Oscillatory scheme) reconstructions. As stated by Gallerano et al. [16], by means of this choice, it is possible to avoid numerical errors that would arise with one-dimensional reconstructions, in the context of generalised curvilinear coordinate systems.
The bottom modifications are calculated whereby the advection-diffusion equation for the suspended sediment concentration. The advective sediment transport terms in the advection-diffusion equation are formulated by using a quasi-three-dimensional (Q3D) approach [23,24]. The depth integration of the product of the horizontal velocity and the concentration vertical distribution allows calculating the advective sediment transport terms and taking into account the sediment transport related to the undertow. To calculate the bottom elevation variation over time, the spatial variation of the bed load transport and the product between the settling velocity and the difference between reference and actual sediment concentration (i.e., the concentration at a given distance from the sea bottom) are used.
The proposed model for the simulation of the sea bottom changes was applied to the abovementioned case study of the Pescara Harbour (Italy). In this paper, the effectiveness of a submerged breakwater designed to protect the entrance of the Pescara Harbour from sediment siltation is investigated by the proposed numerical model.
The present work is structured as follows. In Section 2, the hydrodynamic model and morphodynamic model are presented. In Section 3, we present the case study of Pescara Harbour (Italy). In Section 4, the discussion of the results is presented. The conclusions of the study are drawn in Section 5. At the end of the paper, there are appendices.

Hydrodynamic Model
Let h be the local still water depth and η be the free-surface elevation with respect to the still free surface. Let H = h + η be the total local water depth.
The fundamental idea at the basis of the Boussinesq equations is the elimination of the vertical coordinate from the flow equations and retaining some of the effects of the vertical structure of the flow. The mathematical procedure commonly used in the standard Boussinesq approximation is based on the assumption of irrotational flow and can be summarised by the following steps. A Taylor expansion is made of the horizontal and vertical velocity around a certain elevation. The Taylor expansion is truncated to a finite number of terms. The conservation of mass for an incompressible flow and the zero-curl condition are used to replace vertical partial derivatives of hydrodynamic quantities in the Taylor expansion with horizontal partial derivatives.
In the recent literature, to extend the applicability of the Boussinesq-type models to the surf zone, as well as to simulate the breaking-generated, horizontal, rotational flow, some authors derived new form of the Boussinesq equations by a procedure that allows taking into account the vertical vorticity [20][21][22]. In fact, as demonstrated in [21,22], fully nonlinear Boussinesq equations derived from the assumption of potential flow can be readily converted to equations that allow for the partially rotational flow with vertical vorticity. Wei et al. [20] retained the leading-order vertical vorticity in their equations. Chen et al. [21] and Chen [22] improved the model of Wei et al. [20] to obtain the conservation of vertical vorticity correct to second order.
Let z be the vertical coordinate. Let us employ a Taylor expansion of the velocity about an arbitrary distance from the still water surface; by assuming zero horizontal vorticity, the vertical distribution of the Cartesian based horizontal velocity → U(z) can be written as where → u + is the Cartesian based horizontal velocity at the chosen arbitrary distance σ from the still free surface and  Let ∇ be the two-dimensional differential operator defined as ∇ = (∂/∂x, ∂/∂y) in a Cartesian coordinate system. We define the symbol ⊗ as the outer product in a Cartesian coordinate system.
Let G be the gravity acceleration and We define η c as an arbitrary constant value. In this paper, to obtain a "well-balanced" numerical scheme, the surface gradient term is decomposed as follows In presence of irregular seabed, this decomposition ensures a balance between flux and source terms and allows the numerical scheme to satisfy the C-Property.
Furthermore, the term → V on the right-hand side of Equation (4) is rewritten as follows where → V and → V are defined in Appendix C.
Let us define an auxiliary variable Let x l = x l (ξ 1 , ξ 2 ) (with l = 1, 2) be a transformation from the Cartesian coordinate system (with m = 1, 2). The Jacobian of the transformation is √ g = det(g lm ). Let g LetĤ be the cell average value of H given bŷ By using Equations (5) and (7) and the chain rules [25] that transform the Cartesian based differential operators in the differential operators defined in a generalised curvilinear coordinate system, Equation (4) can be rewritten as The integral form of the continuity Equation (6) reads Equations (10) and (11), in which the vector and tensor quantities are expressed in Cartesian components, are a novel integral form of the fully nonlinear Boussinesq equations written in a generalised curvilinear coordinate system. Equations (10) and (11) are accurate to O µ 2 , ε 3 µ 2 for the dispersive terms and conserve vertical vorticity with a leading-order error of O µ 4 .
In Equations (10) and (11) the conserved variables are H and → M. Therefore, the momentum equation is different from the one presented by Cannata et al. [16] and, unlike Gallerano et al. [16], the continuity equation does not contain dispersive terms. This allows us to entirely discretise the continuity equation by means of a high-order shock capturing finite volume scheme. By doing so, errors in the numerical solution due to the discretisation of the dispersive term in the continuity equation, are avoided. The surface gradient term has been rewritten, in order to discretise it by a finite volume technique and obtain a "well balanced" numerical scheme.
Lynett [19] proposed a correction to the vertical distribution of the horizontal velocity calculated by depth-averaged motion equations ( → u + ), in order to represent the undertow phenomenon. Let → u B (z) be the corrective velocity vector [19]. The horizontal velocity → u (z) is given in Appendix D. Ω(t) is the thickness of the boundary layer, U Ω (t) is the horizontal velocity at the top of the wave boundary layer and u f (t) is the friction velocity. u f c is the current friction velocity and it can be found in Appendix E. By integrating the momentum equation in the boundary layer and from the logarithmic law of the velocity profile, we calculate: Ω(t), U Ω (t), u f (t) and u f c . We indicate by the mark [ ] the wave average over the wave period T of the instantaneous quantities. This calculus procedure is shown in Appendix E.
The eddy viscosity, within the boundary layer, is given by The eddy viscosity, outside the boundary layer, is given by With breaking waves, the turbulence is produced by the current, the wave boundary layer and the wave breaking. Let l be the turbulence length scale, k t = k t (z, t) the turbulent kinetic energy induced by the wave breaking and c d equal to 0.08. The kinetic energy production P k is evaluated according to Deigaard et al. [26]. In the calculation of the eddy viscosity ν t, f (z, t) related to the wave breaking, the turbulent kinetic energy equation by Deigaard et al. [26] is considered The value of k t , obtained with Equation (14), is used to evaluate the eddy viscosity produced by the wave breaking, as follows The total eddy viscosity ν t (z, t) can be defined as the quadratic sum of two terms: the eddy viscosity related to the current and wave breaking and the eddy viscosity produced by the wave boundary layer. Therefore, ν t (z, t) is evaluated as follows By solving the system of equations composed by Equations (10) and (11), we obtain the cell averaged valueĤ and the cell averaged value→ u + . By introducing WENO reconstructions, we obtain the point values H and → u + . By using these point values, we obtain the vertical distribution of the horizontal velocity → u (z). By using u f (t), we obtain the current friction velocity u f c . From Equations (12) and (13), the vertical distribution of the eddy viscosity, ν t,r (z), is calculated inside and outside the boundary layer. The turbulent kinetic energy is calculated by Equation (14). The eddy viscosity produced by the wave breaking, ν t, f (z, t), is calculated by Equation (15). By using ν t,r (z) and ν t, f (z, t) in Equation (16), we obtain the total eddy viscosity ν t (z, t).
The instantaneous values of → u (z), u f , ν t and H are used as input variables of the morphodynamic model that is shown in Section 2.2.

Morphodynamic Model
Let C (z) be the wave-averaged value of the suspended solid concentration. Let the mark [ ] indicate the depth averaged of quantities. Let C be the depth averaged of C (z), defined as where a is a distance from the bottom which delimits the bed load transport region. Moreover, letĈ be the cell average value of C. H is the wave-averaged value of H and Ĥ is the wave-averaged value ofĤ.
In the context of a quasi-three-dimensional methodology, the integral form of the solid particle concentration equation, expressed in a generalised curvilinear coordinate system, reads In Equation (18), u x j (z) is the wave-averaged Cartesian component projected onto the x j -axis of the vector → u (z), ν t is the wave-and depth-averaged total eddy viscosity, De is the rate of the sediment deposition and P is the rate of turbulent sediment pick-up. Let C a (t) and C R (t) be, respectively, the actual concentration and reference concentration. We define C a and C R as the wave-averaged value of the actual and reference concentration, evaluated at height from the bottom a = 2d 50 (d 50 is the sediment mean diameter). C R is evaluated in Appendix F.
De and P are given by In Equations (19) and (20), w sed is the sediment fall velocity.
To solve Equation (18), the solutions of Equations (19) and (20) are needed. Let ν t (z) be the wave-averaged value of the total eddy viscosity ν t (z). The value of C a is taken as the lower boundary condition of the turbulent suspended sediment diffusion equation, defined as follows and as lower extreme of the integral that gives the depth-averaged value, C, defined by Equation (17). An iterative procedure is used to calculate C a , using Equation (17), in which the values of ν t (z) and C are known from the previous time step. The total sediment transport is evaluated as the sum of the suspended load transport and of the bed load transport, which takes into account the near bed transport mechanism. Let p be the sediment porosity, z f the bottom elevation, and ( q b ) x j the Cartesian based component projected onto the x j -axis of the bed load transport vector The equation of the bottom modification over time, expressed in a generalised curvilinear coordinate system, reads A sequence of five steps schematises the way in which the wave and current model is coupled with the model for the simulation of the bottom changes: (1) The instantaneous hydrodynamic quantities are computed by the hydrodynamic model presented in Section 2.1. The simulation time has a duration of about 150 times the mean wave period. (2) The instantaneous hydrodynamic quantities obtained by Step 1 are averaged over the period T * (which is 150 mean wave periods [27]).
(3) The wave-averaged hydrodynamic quantities obtained by Step 2 are used as input for the solid particles' concentration in Equation (18). Equation (18) is numerically integrated and gives the suspended sediment concentrationĈ .
(4) The concentration valuesĈ (that are depth and wave averaged) are used to compute (by WENO reconstructions) C and the values wave-averaged actual concentration C a . (5) The wave-averaged reference concentration C R (computed by Step 1), the wave-averaged actual concentration C a (computed by Step 4) and the bed load transport → q b are used as input for the equation of the bottom modification (Equation (22)). We define the time step of the integration of the bed change equation as "morphological time step". By numerically integrating Equation (22), it is possible to update the bathymetry for the successive morphological time step. The morphological time step is greater than the simulation time interval of the hydrodynamic model. The morphological time step is evaluated, by means of a trial-and-error procedure, by a posteriori verifying the numerical results.

Pescara Harbour Case Study
The case study of the Pescara Harbour, located at the mouth of the homonymous river, represents an example of how coastal structures can modify the nearshore current patterns and consequently the sediment transport and siltation phenomena.
The Pescara Harbour is located in the central Italy and overlooks the Adriatic Sea, as shown in Figure 1.     The preeminent directions of the incoming waves representative of Pescara annual wave climate [28] are shown in Figure 3. In this figure, it is possible to notice that there are two preeminent directions: 345 • N-15 • N (primary sector) and 65 • N-95 • N (secondary sector). The primary sector incoming waves (blue arrow) produce a longshore current coming from the northwest direction (green arrow) and the sediment transport produced by the above current comes from the same direction. The secondary sector incoming waves (orange arrow) produce a longshore current coming from the southeast direction (red arrow) and a sediment transport coming from the same direction. Thus, the potential sediment transport is bimodal and the total resultant is directed from northwest to southeast.
The wave fronts coming from north direction (before the construction of the detached breakwater) caused the navigation to be difficult at the entrance of the Canal Port. To protect the entrance of the Canal Port from the wave motion and ensure good navigability conditions, in 1997, the detached breakwater was built. After the construction of the detached breakwater, a significant increase of the sea bottom level has occurred in the area in the front of the entrance of the Canal Port. The impact of the abovementioned detached breakwater on the bottom modifications in front of the entrance of the Canal Port is negative. The wave motion coming from north direction induces longshore currents coming from northwest direction. The solid material, put into suspension by the wave breaking, is transported by the longshore current into the area in front of the entrance of Canal Port. In this area, the energy of the wave motion is consistently damped by the presence of the detached breakwater. Consequently, the solid material silts near the entrance of the Canal Port, where it produces an increment of the sea bottom level and makes it difficult to navigate. The sediment transport coming from southeast is mostly intercepted by the Touristic Port and by the Eastern Jetty, before reaching the proximity of the entrance of the Canal Port. Consequently, the solid material coming from southeast does not significantly contribute to the silting of the area in front of the Canal Port. The silting of this area is due, in minor part, also to the sediment transport coming from the Pescara River during flood events.
In summary, the area in front of the entrance of the Canal Port, bounded by the detached breakwater, the Eastern Jetty and the Touristic Port, undergoes a process of silting, mainly caused by the sediment transport coming from northwest. The access to the Canal Port, in absence of periodic appropriate dredging operations, is prevented in the northwest harbour gate, due to the formation of extended areas with still water depth less than 1 m. The access to the Canal Port is possible, with difficulty, only by routes approaching the entrance from east-northeast and only by vessels with draft less than 2 m.
In Section 3.1, we present the numerical simulations carried out to validate the hypothesis that the solid material put into suspension by the wave breaking is transported by the longshore current coming from northwest and silts near the entrance of the Canal Port.
Recently, the Italian Ministry of Public Works has considered the possible construction of a submerged breakwater that is able to act as an obstacle to the flow of solid material carried by the longshore current coming from northwest. Such solid material, as demonstrated in the present paper, settles to the northwest of the submerged breakwater. The reasons at the basis of the decision to verify the functionality of the submerged breakwater lie in the fact that (in the present configuration, without the submerged breakwater) the solid material coming from northwest mixes with the one coming from the Pescara River and settles in front of the entrance of the Canal Port, by increasing the sea bottom level. The sediment coming from northwest is generally of good quality, while the one coming from the Pescara River is polluted. The continuous dredging that should be done in front of the entrance of the Canal Port would not provide good quality sediment to be used for possible beach nourishment. According to the Italian legislation, the polluted dredged sediment must be treated, with consequent very high costs. The presence of the submerged breakwater should block the pollution-free sediment coming from northwest, which could be used for the nourishment of the touristic beach at the northwest of the submerged breakwater. Furthermore, as stated above, the submerged breakwater, by blocking the sediment coming from northwest, would limit the increase of sea bottom level at the entrance of the Canal Port, which is mainly due to sediment carried by the above mentioned longshore current. Another design solution, which consists in the construction of a system of long groins that are perpendicular to the shoreline and placed on the northwest side of the Canal Port entrance, has been considered. This solution has been rejected by the Italian Ministry of Public Works after a cost-benefit analysis because it is too expensive and is characterised by a main environmental impact.
In Section 3.2, we present the results of the investigation over the effectiveness of the introduction, in the coastal area of the Pescara Harbour (Italy), of the submerged breakwater. Figure 4a shows the present configuration (Configuration 1), in which there is only the detached breakwater (drawn in red). Figure 4b shows the designed configuration (Configuration 2), in which the submerged breakwater (drawn in yellow) connects the detached breakwater to the shoreline. The dimensions of the submerged breakwater are: crest width equal to 3 m, larger base equal to 23 m and slope of banks equal to 1 : 2. The maximum depth where the breakwater is placed is equal to 5 m and the freeboard is equal to 0.5 m. The submerged breakwater is modelled as a bottom increase.
In both figures, the still water depth measured from mean water level is reconstructed based on bathymetric data which have been obtained after dredging operations. The still water depth in Figure 4 (measured from mean water level) is assumed as initial conditions for the simulations. The primary sector incoming waves are characterised by a significant wave height of 1.5 m, which has an occurrence frequency at least 320 h/year. By the proposed model, the simulations of the wave motion and the sea bottom changes that occur over three years in both Configurations 1 and 2 were carried out. To reproduce the primary sector wave forcing, in the numerical simulations, we generated, as initial conditions, random waves incoming from 0 • N, represented by a spectrum belonging to the JONSWAP type characterised by a significant wave height equal to 1.5 m, which act for 320 h/year for a total of 960 simulated hours. Figure 5 shows an instantaneous wave field related to the Configuration 1 whose initial still water depth is shown in Figure 4a. In Figure 5, it is possible to notice that the waves, between Y = 100 m and Y = 700 m, at first undergo an increase of the wave height and a reduction of the wavelength (shoaling) and then a decrease of the wave height due to the wave breaking. In particular, we note that the waves, intercepted by the west edge of the detached breakwater, undergo a rotation (diffraction) and, proceeding to the coast, a reduction of the wave height due to the breaking. In the same figure, it is possible to see that, between Y = 700 m and Y = 1500 m, the waves are intercepted by the detached breakwater and the Eastern Jetty and undergo a reflection effect due to such structures.  Figure 6 shows a three-dimensional view of the wave field that is shown in Figure 5. In Figure 6, it is possible to notice that the numerical simulation produces wave heights equal to 0.40 m in correspondence of the entrance of the Canal Port: these reduced values of the wave height underline that the detached breakwater improves the navigability conditions at the entrance of the Canal Port.  Figure 7 shows wave-and depth-averaged velocity field produced by the wave field in Figure 5. In Figure 7, it is noted that, in Configuration 1, between Y = 400 m and Y = 1000 m, the simulation produces a longshore current that is parallel to the shoreline, characterised by maximum values of the velocity equal to 0.45 m/s. Furthermore, between Y = 1000 m and Y = 1500 m, the abovementioned coastal current passes through the detached breakwater and the Eastern Jetty, where it reaches maximum values of the velocity equal to 0.6 m/s.  Figure 8a shows the still water depth (measured from mean water level) produced by the wave field in Figure 5, at the end of the first simulated year in Configuration 1. The still water depth initial conditions used in this numerical simulation are the ones shown in Figure 4a. Compariing  Figures 4a and 8a, it is noted a decrease of the still water depth (with consequent sea bottom level increase) in the area between about X = 500 m and Y = 700 m and in the area between X = 450 m and X = 550 m, Y = 950 m and Y = 1050 m. In this area, at the end of the first simulated year, the 3 m still water depth contour line moves toward the north, with respect to the initial conditions. This comparison also shows a still water depth contour line decrease (with consequent sea bottom level increase) in front of the entrance of the Canal Port. Comparing Figures 4a and 8a shows a decrease of the still water depth (with consequent sea bottom level increase) in the area between Y = 700 m and Y = 1000 m, in proximity to the shoreline, where an offshore shift of the 2 m still water depth contour line occurs, from X = 800 m to X = 700 m. Figure 8b presents the still water depth, at the end of the third simulated year in Configuration 1, produced by the wave field in Figure 5. Comparing Figure 8a,b, it is possible to notice a further decrease of the still water depth (with consequent sea bottom level increase) in the area between X = 500 m and between Y = 700 m and Y = 850 m. In this area, at the end of the third simulated year, the 3 m still water depth contour line shifts further toward the north. Furthermore, by the same comparison, it is highlighted that there is a further decrease of the still water depth in front of the entrance of the Canal Port. Figure 8b shows that the area (between X = 675 m and X = 750 m, Y = 925 m and Y = 1100 m) bounded by the 1 m still water depth contour line increases with respect to the one presented in Figure 8a. The same figure also shows that the 2 m still water depth contour line further advances toward the northeast, by reaching Y = 1100 m at the end of the third simulated year. Comparing the results at the end of the first and third simulated years also shows a further decrease of the still water depth (with consequent sea bottom level increase) in the area between Y = 650 m and Y = 800 m and close to the shoreline, where an offshore shift of the 2 m still water depth contour line occurs. The numerical simulations are conducted without taking into account the supply of sediment from Pescara River.  Figures 9-11 show, respectively, the wave field, the wave-and depth-averaged velocity field and the still water depth, measured from mean water level, related to Configuration 2 characterised by the presence of the submerged breakwater.   In Figure 9, an instantaneous wave field, produced by the same initial conditions used for the simulations of the wave field in Figure 5, is shown. Configuration 2 is characterised by the still water depth shown in Figure 4b. In particular, it can be noticed that the waves, intercepted by the west edge of the detached breakwater, undergo a rotation (diffraction), then proceed towards the shoreline and undergo a reduction of the wave height due to the wave breaking. The same figure shows that the waves diffracted by the west edge of the detached breakwater pass over the submerged breakwater and undergo a reduction of the wave height. In Figure 9, it can be noticed that there are no significant modifications of the wave height in correspondence of the submerged breakwater. This is in marked contrast to the real situation. The reduction of the simulated wave heights is due to the lack of refinement of the computational grid near the submerged breakwater. The fundamental aim of the present study was related to the evaluation, from a global point of view, of the effects on the bottom changes due to the submerged breakwater. The detailed analysis on the wave heights in the area close to the submerged breakwater and on the top of the submerged breakwater is not required for the abovementioned global analysis. The detailed analysis requires the use of small spatial and temporal discretisation steps that increase the computational time very much.

Configuration 2
In Figure 10, the wave-and depth-averaged velocity field, produced by the wave field shown in Figure 9, is presented. Figure 10 shows that, in the area between the submerged breakwater and the Eastern Jetty, the maximum value of the velocity is equal to 0.55 m/s. The wave-and depth-averaged velocities on the submerged breakwater are significantly greater than the ones in all the domain and, for this reason, they are not represented in the figure. Figure 11a shows the still water depth produced by the wave field in Figure 9, at the end of the first simulated year in Configuration 2. The still water depth initial conditions used in this numerical simulation are the ones shown in Figure 4b. Comparing Figures 4b and 11a, it can be noticed that, in the area between the submerged breakwater, the detached breakwater and the entrance of the Canal Port, the still water depth is slightly different from the initial conditions. From the same comparison, it can be also noticed that, in the west area of the submerged breakwater, there is an increase of the sea bottom level. In the area between about X = 700 m and X = 1050 m, the 1 m, 2 m and 3 m still water depth contour lines advance in the direction of the detached breakwater, with respect to the initial conditions.
In Figure 11b, the still water depth, at the end of the third simulated year, produced by the wave field in Figure 9, is shown. Comparing Figure 11a,b, it can be noticed that, in the area between the submerged breakwater, the detached breakwater and the entrance of the Canal Port, the still water depth is slightly different from the initial conditions. The only contribution of solid material that could modify such sea bottom level is the one coming from the Pescara River. This solid material was not taken into account in the numerical simulation. Comparing Figure 11a,b, it can be noticed that, in the west area of the submerged breakwater, at the end of the third simulated year, there is an accretion of the sea bottom level. In the abovementioned area and between about X = 650 m and X = 1050 m the Figure 11b shows that the 1 m, 2 m and 3 m still water depth contour lines advance in the direction of the detached breakwater. There is a significant increase of the seas bottom level in the area between X = 800 m and X = 950 m and between Y = 660 m and Y = 750 m, where the still water depth is 1 m.

Discussion
In this section, we present summary considerations on the results of the numerical simulations relative to the sea bottom changes after three simulated years, either in absence or in presence of the submerged breakwater.
As stated above, the wave motion coming from north direction induces longshore currents coming from northwest direction. The solid material, put into suspension by the wave breaking, is transported by the abovementioned longshore current in the area in front of the entrance of Canal Port. In this area, the energy of the wave motion is consistently damped by the presence of the detached breakwater. Consequently, the solid material silts near the entrance of the Canal Port. The siltation increases the sea bottom levels in front of the entrance of the Canal Port where makes it difficult to navigate. Figure 12 shows the still water depth contour lines after three simulated years, for Configuration 1. In this figure, the main accretion areas are shown. In the initial configuration (Figure 4a), in front of the entrance of the Canal Port, the still water depth with respect to the mean sea level is about 3 m. After three simulated years, in the area coloured in blue in Figure 12, the sea bottom levels increase by about 1 m on average, with a maximum increase of about 2 m, and the volume of accumulated sediment is about 25,000 m 3 . The area in front of the beach, at the northwest of the entrance of the Canal Port, indicated by the red colour in Figure 12, shows an increase of the sea bottom levels by about 1 m on average, with a maximum increase of about 2 m, with respect to the initial configuration in Figure 4a. In this area, the volume of accumulated sediment is about 65, 000 m 3 (after three simulated years). Finally, the area coloured in green in Figure 12, located close to the northwest edge of the detached breakwater (after three simulated years), shows an increase of the sea bottom levels by about 1 m, with respect to the initial configuration shown in Figure 4a. In this area, the volume of accumulated sediment is about 20,000 m 3 . After the construction of the detached breakwater, periodic measurement campaigns of the still water depth have been carried out by the Italian Ministry of Environment and Protection of Land and Sea [29] in the area between the detached breakwater and the entrance of the Canal Port. In Figure 13, three areas, indicated by A-C, are shown. This classification comes from the fact that, by the abovementioned measurement campaigns, the Italian Ministry of Environment and Protection of Land and Sea carried out an evaluation of the volume of sediment that have accumulated over time in Zone A. The quantification carried out by the measurement campaigns indicates a volume of accumulated sediment of about 40,000 m 3 /year.
The volumes of accumulated sediment obtained by the numerical simulations (carried out without the sediment contribution coming from the Pescara River) are: 34,200 m 3 /year in Zone A, 24,500 m 3 /year in Zone B and 22,500 m 3 /year in Zone C. The volume of sediment coming from Pescara River is estimated as 10,000 m 3 /year by [29]. By taking into account this contribution, the results obtained by the proposed model are consistent with the measurements carried out by Italian Ministry of the Environment and Protection of Land and Sea [29].    Figure 15 shows the areas in which the volumes of accumulated sediment are calculated by the numerical simulations. In this configuration, the area previously called A is divided into two sub-areas, A1 and A2, because of the presence of the submerged breakwater. The volumes of accumulated sediment obtained by the numerical simulations (carried out without the sediment contribution coming from the Pescara River) are: 7000 m 3 /year in Zone A1, 11,100 m 3 /year in Zone B and 16,100 m 3 /year in Zone C. The designed submerged breakwater represents an obstacle to the solid material transport by the longshore currents. This solid material tents to settle in Zone A1, B and C (corresponding to the area bounded by the designed submerged breakwater in the east, by the detached breakwater in the north and by the shoreline in the south). In absence of the submerged breakwater, the sediment coming from northwest mixes with the sediment coming from the Pescara River and settles in front of the Canal Port. The sediment coming from the northwest is of good quality, while the one coming from the Pescara River is polluted and cannot be used for beach nourishment. The presence of the submerged breakwater blocks the pollution-free sediment that can be used for beach nourishment. Furthermore, as stated above, the submerged breakwater, blocking the solid material coming from northwest, significantly limits the increase of the sea bottom level in front of the entrance of the Canal Port and, consequently, makes it possible to maintain good navigability conditions.

Conclusions
The real case study of Pescara Harbour (Italy) was investigated by means of a numerical bottom-change simulation model, composed by a hydrodynamic sub-model and a morphodynamic sub-model. In particular, the effects of the introduction of a designed submerged breakwater in the existing context were studied in terms of morphodynamic alterations. The used hydrodynamic model is based on an integral form of the fully nonlinear Boussinesq equations, while the proposed morphodynamic model is based on a quasi-three-dimensional approach. It was demonstrated (by numerical simulations) that the detached breakwater, built to reduce the navigation risks at the entrance of the Canal Port, damps the energy of the wave motion and causes the siltation of the solid materials, transported by the longshore current coming from northwest, in front of the harbour entrance. From the numerical investigation carried out, it can be seen that the designed submerged barrier substantially modifies the hydrodynamic field and the morphodynamic evolution of the zone near Pescara Harbour. In particular, the ability of the designed submerged breakwater, located orthogonally to the longshore current, to intercept the aforementioned solid material and to significantly reduce the accretion of the bottom in the area in front of the entrance of the Canal Port entrance, was verified.
Author Contributions: Conceptualisation and methodology: F.G.; Validation, investigation, writing and editing F.P. and B.I.; and Supervision F.G. All authors have read and agreed to the published version of the manuscript.
Funding: This research received external funding.

Conflicts of Interest:
The authors declare no conflict of interest.
The right-hand side of Equation (A1) is the second-order term in power series expansion of the

Appendix B
The expressions of → V and → T are respectively: = v x 1 , v x 2 , 0 be the vectors defined in a three-dimensional Cartesian coordinate system; let the symbol ∧ be the vector product in a three-dimensional Cartesian coordinate system; and let ∇ 3 be the three-dimensional differential operator given by ∇ 3 = (∂/∂x, ∂/∂y, ∂/∂z ) in a Cartesian coordinate system. We define By using Equation (9) (17) and (19) in Section 2.

Appendix E
Let k be the bed roughness; let K be the von Karman constant; and let k/30 be the characteristic length scale.
By integrating the momentum equation in the boundary layer and from the logarithmic law of the velocity profile, we have The thickness of the boundary layer can be obtained by Therefore, the values of u f (t) and Ω(t) are obtained by solving the system composed by Equations (A10) and (A11). Inside the boundary layer, turbulence is produced by the interaction between waves and current. Let u f c be the current friction velocity, which is given by

Appendix F
Let C m be the maximum volumetric concentration that can be reached, θ cr , and let → θ = → θ(t) be the parameter of the stability and mobility of Shield. Let → θ(t) be the bed shear stresses induced by wave and current.
C R is the wave-averaged value of the reference concentration and it is evaluated by averaging over the wave period the instantaneous values of C R (t), which is calculated according to [27]:

Appendix G
Let β be the dynamic friction coefficient and let ρ s /ρ w be the ratio between the sediment density and water density.
The wave-averaged value of → q b is given by