A Simpliﬁed Method for an Evaluation of the E ﬀ ect of Submerged Breakwaters on Wave Damping: The Case Study of Calabaia Beach

: Erosion processes threaten the economy, the environment and the ecosystem of coastal areas. In addition, human action can signiﬁcantly a ﬀ ect the characteristics of the soil and the landscape of the shoreline. In this context, pursuing environmental sustainability is of paramount importance in solving environmental degradation of coastal areas worldwide, with particular reference to the design of complex engineering structures. Among all the measures conceived to protect the shoreline, environmentally friendly interventions should be supported by the stakeholders and tested by means of mathematical models, in order to evaluate their e ﬀ ectiveness in coastal protection through the evaluation of wave damping and bedload. This study focuses on protected nourishments, as strategic interventions aimed to counteract coastal erosion without a ﬀ ecting the environment. Here, we develop a simpliﬁed method to provide a preliminary assessment of the e ﬃ ciency of submerged breakwaters in reducing wave energy at a relatively low computational cost, if compared to the standard 2D or full 3D mathematical models. The methodology is applied at Calabaia Beach, located in the southern Tyrrhenian Sea (Italy), in the area of the Marine Experimental Station of Capo Tirone. The results show that the simpliﬁed method is proven to be an essential tool in assisting researchers and institutions to address the e ﬀ ects of submerged breakwaters on nourishment protection.


Introduction
The erosion process threatens urban settlements, the environment and the ecosystem located in the Mediterranean Sea [1,2]. Coastal areas involve strong relationships between the needs of human communities and the environment [3]. The increasing demand of coastal use should aim to conserve the resilience of the environment, allowing the ecosystem to absorb the unavoidable human impacts, keep functioning, and provide goods and services to the population [4][5][6]. However, coastal areas regularly deal with complex environmental and ecosystem challenges, which are worsened by mean sea level rises driven by climate change, subjecting these fragile lands to large economic and environmental losses and damage produced by flooding and erosion processes [7][8][9][10][11][12][13][14]. Increasing urban pressure will further exacerbate these damages, requiring an integrated approach to conserve the land and the environment [15,16], evolving from sea defence interventions only aimed at protecting human and goods, to coastal protection measures focused on the overall needs of the land [17].

Materials and Methods
The design of submerged breakwaters as a support for protected nourishments first requires the evaluation of the effect of the barrier on the reduction in the wave height from offshore to the shore by means of the wave transmission coefficient Kt [37].
where H out and H in stand for the wave height upstream and downstream of the barrier, respectively [38]. Although this simplified approach does not describe the variability of the wave climate and the possible effect of the anthropic structures, it is consistent with the aim of this study.

The Mathematical Model
In recent decades, several models have been developed to compute the morphodynamic evolution of coastal environments, characterized by several computational schemes, scales, levels of detail and computational costs [39][40][41][42][43][44]. The spatial scale of these models ranges from meters to kilometres. The temporal scale typically ranges from hours to decades. Most of the existing large-scale models are limited to an assessment of the evolution of the shoreline, while the smaller space and time scale models are mainly used to reproduce specific forcing events, involving explicitly reductionist methodologies where the conservation of momentum forms the explicit means for the evolution of the system. [45][46][47].
Here, we use MIKE 21-3 Coupled Model FM, based on an unstructured grid which consists of an integrated system capable of reproducing the coastal processes and dynamics of the shoreline at all scales. MIKE evaluates the effectiveness of the sea defence interventions, such as the optimisation of the beach nourishments and costal protection structures, and the impact of the buildings located in the active beach.
The Hydrodynamic (HD) module of MIKE solves the 2D Navier-Stokes equations of incompressible fluids, under the hypothesis of hydrostatic pressure [48]. The numerical solution of shallow water equations is achieved through the approximate Riemann solver [49], which calculates the convective flows at the interface of the element of the grid. Second-order precision is achieved through the use of a linear gradient reconstruction technique. The average gradient is computed by means of the Jawahar and Kamath [50] approach.
The Spectral Wave (SW) module describes the wave phenomena by means of the wave action conservation equations [39,40], which can be written by using the Cartesian system as: where N (θ, σ) is the density of the wave, which is function of θ, wave direction and wave frequency σ. S is the source term accounting for different processes: where S in is the wind action; S nl is the waves non-linear interaction; S ds is the white capping effect; S bot is the bottom dissipative action; S surf is the dissipation due to wave breaking phenomena.
Among all these parameters, S surf is essential in representing the wave breaking, which occurs when the waves propagate in shallow water and can no longer be supported by the depth. This process is described following the formulation of Battjes and Janssen [51].
If the shoreline is characterized by shallow waters surrounding the shoreline, the maximum wave height can be calculated as: where d is the water depth and γ a coefficient which ranges between 0.5 and one, depending on the slope of the bottom and the wave characteristics [52]. The formulation of γ can be chosen according to Nelson's formulations [53,54], or by means of the more recent Ruessink approach [55], valid both for steep seabeds without bars and flat bottoms. The Ruessink approach calculates γ at each cell of the domain as: where k is the local wave number and d the local depth. The bedload is solved by means of the Sand TransPort Quasi-3-Dimensional module (STPQ3D) [56], which computes the transport of the sediments along two horizontal directions by processing the input data provided by the hydrodynamic model by time-averaging the results.
Waves affect the motion of sediments in shallow waters due to wave breaking, with particular reference to the transport of sediments transversal to the coast, which is mostly involved in the erosion process [57,58]. Morphological changes due to erosion and deposition are accounted in terms of bed elevation, expressed with respect to the centre of each triangular element of the grid as δz δt . This parameter can be computed by means of the continuity equation of the sediment proposed by Exner [59], where n is the porosity: The bedload is computed as the divergence of the sediment flow with respect to the boundary of the elements and it is equal to the sum of all the flows that cross the boundaries, identifying whether the element is receiving or losing sediment.
where S in is the sediment flux normalized to the face of the element; ds i is the lenght of the face of the element; m is the index of the element; The continuity equation of sediments affects the bed elevation of each element, based on the total sediment flux crossing the boundaries of each element.
The coupling of the models is based on the establishment of an overall time step of the system, which is necessary to match the instant in which the models exchange the information. Each model has its own time step, which is a multiple of the overall time step. The exchange of information occurs when the time of each model match with an overall time step (e.g., the HD module acquires a radiation stress field from the SW module).
Since the information sharing between each module of MIKE is very complex, we refer to the manuals for a more detailed description [48,56,60]. A summary of the main steps is provided in Figure 1.
which is necessary to match the instant in which the models exchange the information. Each model has its own time step, which is a multiple of the overall time step. The exchange of information occurs when the time of each model match with an overall time step (e.g., the HD module acquires a radiation stress field from the SW module).
Since the information sharing between each module of MIKE is very complex, we refer to the manuals for a more detailed description [48,56,60]. A summary of the main steps is provided in Figure  1.

Submerged Breakwaters
The effects of submerged breakwaters on wave damping is typically evaluated through the formulation of Goda [61,62], which computes the transmission coefficient Kt (1) as: The values of K t, max and K t, min are set according to the boundary conditions; the parameter f (i.e., the free board) represents the height of the breakwater crest minus the surface level; H i the value of the incoming wave height, i.e., the input wave height. If the breakwater is submerged, the wave damping effect is negligible when the wave height is equal to or smaller than the depth of the crest, computed with respect to the mean sea level. For parameter α, Goda proposed a value equal to 2.2, while for parameter β, which depends on the inclination of the breakwater, a value ranging between 0.15 and 0.8 [61,62]. Goda's formula allows to evaluate the performance of the breakwater at different water levels, and it is not affected by the angle of attack of the waves with respect to the breakwater.
2D and 3D modeling of long and narrow structures such as the submerged barriers significantly increases the computational effort to avoid model instability. In this study, we propose a simplified method based on the implementation of Goda's formula directly on the domain, avoiding the 2D modeling of the barrier but keeping the high accuracy of the results. Similar approaches are widely implemented into mathematical models available in the literature by developing a set of specific links to reproduce the presence of long and narrow structures, such as sills and levees [63][64][65][66].
Here, we link the nodes of the elements located upstream and downstream of the submerged barrier by implementing Goda's formula. We calibrated parameter β of Goda's formula, by comparing the results computed by means of the full 2D model, using different depths of its ridge (z) and different wave climates (H m0 ).
In the simulations performed by means of the full 2D model, we reproduced the rigid submerged breakwater by means of 2D elements ( Figure 2). The hydrodynamic of the barrier is solved by the HD and SW modules, without any change in bottom elevation, by disabling the use of the sedimentological model.
implemented into mathematical models available in the literature by developing a set of specific links to reproduce the presence of long and narrow structures, such as sills and levees [63][64][65][66].
Here, we link the nodes of the elements located upstream and downstream of the submerged barrier by implementing Goda's formula. We calibrated parameter of Goda's formula, by comparing the results computed by means of the full 2D model, using different depths of its ridge (z) and different wave climates (Hm0).
In the simulations performed by means of the full 2D model, we reproduced the rigid submerged breakwater by means of 2D elements ( Figure 2). The hydrodynamic of the barrier is solved by the HD and SW modules, without any change in bottom elevation, by disabling the use of the sedimentological model. The simplified approach we develop allows us to evaluate the effect of the submerged breakwater by reproducing the structure by linking the elements of the grid located before and after the barrier, which works as a weir on whose edges we apply Goda's equation (Figure 3). The simplified approach we develop allows us to evaluate the effect of the submerged breakwater by reproducing the structure by linking the elements of the grid located before and after the barrier, which works as a weir on whose edges we apply Goda's equation (Figure 3). The comparison between the simulations performed by means of full 2D modeling and the simplified scheme allows us to identify the optimum value of of Goda's formula that better reproduces the effects of the breakwaters on the wave damping. The results of the simplified methodology will be synthetized in two different formulations, which aim to relate the wave damping only to the ridge depth of the barrier and to the wave climate, without performing any simulation by means of the models.

The study area: Calabaia Beach
Calabaia Beach is located 1 km south of Belvedere Marittimo, into the Lao region, which is extended 25 km from Capo Scalea to Capo Bonifati. The coast is characterized by sandy and pebbly beaches bounded by narrow dune belts and sedimentary rocks with deposits of conglomerates of pebbles, sand and clay beds [67] (Figure 4). The comparison between the simulations performed by means of full 2D modeling and the simplified scheme allows us to identify the optimum value of β of Goda's formula that better reproduces the effects of the breakwaters on the wave damping. The results of the simplified methodology will be synthetized in two different formulations, which aim to relate the wave damping only to the ridge depth of the barrier and to the wave climate, without performing any simulation by means of the models.

The Study Area: Calabaia Beach
Calabaia Beach is located 1 km south of Belvedere Marittimo, into the Lao region, which is extended 25 km from Capo Scalea to Capo Bonifati. The coast is characterized by sandy and pebbly beaches bounded by narrow dune belts and sedimentary rocks with deposits of conglomerates of pebbles, sand and clay beds [67] (Figure 4).

The study area: Calabaia Beach
Calabaia Beach is located 1 km south of Belvedere Marittimo, into the Lao region, which is extended 25 km from Capo Scalea to Capo Bonifati. The coast is characterized by sandy and pebbly beaches bounded by narrow dune belts and sedimentary rocks with deposits of conglomerates of pebbles, sand and clay beds [67] (Figure 4). A granulometric characterization of the site has been carried out by means of several surveys along the shoreline performed in the years 2003, 2006 and 2008, using the ASTM 200 procedure [68].
For each year of monitoring, 24 surveys were carried out, six on each section located at different depths, allowing us to evaluate the composition of the soil ( Figure 5). The sample analysis shows that the materials that constitute the bottom have similar features, with an average diameter d50 = 0.22 mm, a specific weight of the soil γs/γw = 2.65, and a porosity  = 0.4. The results of the surveys, together  [68].
For each year of monitoring, 24 surveys were carried out, six on each section located at different depths, allowing us to evaluate the composition of the soil ( Figure 5). The sample analysis shows that the materials that constitute the bottom have similar features, with an average diameter d 50   During the year 2002, a protected nourishment intervention was built at Calabaia Beach. The submerged breakwater is located parallel to the shoreline, 700-m long, 250-m seaward and with a ridge 2.5 m below the mean sea level ( Figure 6). Two perpendicular semi-submerged groynes complete the structure, delimiting a submerged closed cell aimed at supporting the nourishment. During the year 2002, a protected nourishment intervention was built at Calabaia Beach. The submerged breakwater is located parallel to the shoreline, 700-m long, 250-m seaward and with a ridge 2.5 m below the mean sea level ( Figure 6). Two perpendicular semi-submerged groynes complete the structure, delimiting a submerged closed cell aimed at supporting the nourishment. During the year 2002, a protected nourishment intervention was built at Calabaia Beach. The submerged breakwater is located parallel to the shoreline, 700-m long, 250-m seaward and with a ridge 2.5 m below the mean sea level ( Figure 6). Two perpendicular semi-submerged groynes complete the structure, delimiting a submerged closed cell aimed at supporting the nourishment.

Model Set-Up
The domain of the model is an unstructured mesh, with triangles of different sizes, which are larger offshore and smaller near the barrier and the shoreline, in order to represent the area of concern in more detail (Figure 9).

Model Set-Up
The domain of the model is an unstructured mesh, with triangles of different sizes, which are larger offshore and smaller near the barrier and the shoreline, in order to represent the area of concern in more detail (Figure 9).

Model Set-Up
The domain of the model is an unstructured mesh, with triangles of different sizes, which are larger offshore and smaller near the barrier and the shoreline, in order to represent the area of concern in more detail (Figure 9). Wave forcing has been imposed on the west side of the domain. The north and the south side of the grid are characterized by lateral boundary conditions [70][71][72], stated in the set-up manual as good approximations when the boundary line is almost straight and when the depth contours are almost perpendicular to the line [60]. The east side of the mesh, on the other hand, has the condition of land imposed in the MIKE mesh generator. For each run, we set a condition for gradual offshore wave formation, obtaining a linear growth in the significant wave height up to the value we simulated, in order to avoid wave energy peaks that would generate instability in the model. When the targeted wave height has been achieved, the simulation begins, with a duration of 24 h and a print step of 600 s, for a total of 144 pieces of output data. For each run, we consider the forcing as a constant wave (in time and space), setting the wave model with a quasi-stationary formulation, which is suitable for shorelines not longer than 10 km [60].

Results
We simulate the efficiency of the submerged barrier for five different offshore waves (Table 1), in order to cover the significant return periods illustrated in Figure 8. A first set of simulations was carried out by representing the submerged breakwater by means of full 2D modeling (Figure 2). The effects of the breakwater were further evaluated by means of a second set of simulations performed by applying the simplified approach through Goda's formula instead of the 2D modeling of the barrier. Using the same boundary conditions, five values of β have been tested with the aim being to achieve the same results computed through full 2D modeling ( Figure 10).
formation, obtaining a linear growth in the significant wave height up to the value we simulated, in order to avoid wave energy peaks that would generate instability in the model. When the targeted wave height has been achieved, the simulation begins, with a duration of 24 hours and a print step of 600 seconds, for a total of 144 pieces of output data. For each run, we consider the forcing as a constant wave (in time and space), setting the wave model with a quasi-stationary formulation, which is suitable for shorelines not longer than 10 km [60].

Results
We simulate the efficiency of the submerged barrier for five different offshore waves (Table 1), in order to cover the significant return periods illustrated in Figure 8. A first set of simulations was carried out by representing the submerged breakwater by means of full 2D modeling (Figure 2). The effects of the breakwater were further evaluated by means of a second set of simulations performed by applying the simplified approach through Goda's formula instead of the 2D modeling of the barrier. Using the same boundary conditions, five values of β have been tested with the aim being to achieve the same results computed through full 2D modeling ( Figure 10).  We compared the results in a control section 1200-m long, perpendicular to the shoreline and located at the mid-point of the length of the barrier, in order to avoid boundary effects ( Figure 10, white line). Figure 11 compares the damping of the wave height produced by the submerged breakwater along the control section, computed by means of full 2D modelling and the simplified 2D approach, with entering waves with H m0 values of 3 m, 5 m and 7 m (Table 1) as boundary conditions (Figure 11).
The results show that simplified 2D modeling is capable of reproducing wave height damping. With the aim of identifying the value of β that better reproduces the effect of the barrier, we compared the wave height reduction due to the breakwater for all the input waves reported in Table 1 (Figure 11d). We compared the height decay efficiency E h (9), setting the height of the offshore significant wave as the average of the first 25 m of the control section (H m0 in ), and the height of the outcoming significant wave (H m0 out ) as the average of the last 25 m of the control section ( Figure 10, white line).
The β coefficient that better represents the effects of the breakwater is 0.30. Table 2 compares the root mean square error (RMS) between the results of the two different sets of simulations, confirming the optimum value for parameter β as equal to 0. 30 We compared the results in a control section 1200-m long, perpendicular to the shoreline and located at the mid-point of the length of the barrier, in order to avoid boundary effects ( Figure 10, white line). Figure 11 compares the damping of the wave height produced by the submerged breakwater along the control section, computed by means of full 2D modelling and the simplified 2D approach, with entering waves with Hm0 values of 3 m, 5 m and 7 m (Table 1) as boundary conditions ( Figure  11). The results show that simplified 2D modeling is capable of reproducing wave height damping. With the aim of identifying the value of β that better reproduces the effect of the barrier, we compared the wave height reduction due to the breakwater for all the input waves reported in Table 1 ( Figure  11d). We compared the height decay efficiency Eh (9), setting the height of the offshore significant wave as the average of the first 25 m of the control section ( ̅ 0 ), and the height of the outcoming significant wave ( ̅ 0 ) as the average of the last 25 m of the control section ( Figure 10, white line).
The β coefficient that better represents the effects of the breakwater is 0.30. Table 2 compares the root mean square error (RMS) between the results of the two different sets of simulations, confirming the optimum value for parameter  as equal to 0.30. Table 2. RMS error of the wave height damping computed through simplified 2D modeling compared to the results obtained by means of the full 2D model.  We further investigated the damping efficiency using different ridge depths for the breakwater z, ranging between −1 and −2.5 m, in order to provide a more general result [61]. Figure 12 illustrates the logarithmic relation between the height and power damping efficiency and the wave height normalized to the ridge depth of the breakwater H m0 z . We selected the range Hm0 z , varying from 0.5 (i.e., the breakwater begins to affect the waves) to 3.0, which is the limit constrained by the depth of the breakwater. We further investigated the damping efficiency using different ridge depths for the breakwater z, ranging between -1 and -2.5 m, in order to provide a more general result [61]. Figure 12 illustrates the logarithmic relation between the height and power damping efficiency and the wave height normalized to the ridge depth of the breakwater 0 . We selected the range 0 , varying from 0.5 (i.e., the breakwater begins to affect the waves) to 3.0, which is the limit constrained by the depth of the breakwater. The regressive parameters m and q can be further related to the ridge depth z, showing a linear relationship (Figure 12c and Figure 12d) and allowing us to relate the damping efficiency only to the wave climate and the depth of the ridge of the breakwater, without performing any simulation (10, wave height damping; 11, wave power damping): The regressive parameters m and q can be further related to the ridge depth z, showing a linear relationship (Figure 12c,d) and allowing us to relate the damping efficiency only to the wave climate and the depth of the ridge of the breakwater, without performing any simulation (10, wave height damping; 11, wave power damping): m h = −0.599 ·z + 0.294 Our findings show that the simplified formulation we propose can successfully support the stakeholders in providing an initial assessment of the effectiveness of the submerged breakwater, with particular reference to its design and the further maintenance measures. Although both Equations (10) and (11) have been computed as site specific depending on the surrounding bathymetry, the prevailing environmental conditions and the configuration of the breakwater, the results can have a general validity since, in the area of concern, there are no other structures that can significantly affect the outcomes. The ranges of wave damping we simulated fall within the range of variability provided in the literature [73,74], both in reference to the wave height and to the wave power, showing that the breakwater efficiency is a function of the relative submergence of the crest.

Application at Calabaia Beach
We finally applied simplified 2D modeling to assess the effectiveness of the submerged breakwater located at Calabaia Beach through the evaluation of the erosion process affecting the shoreline. Using the same boundary conditions illustrated in Section 2.4, we used the model to reproduce the erosion process affecting the nourishment at Calabaia Beach ( Figure 13). Since the area has been subject to a large volume of nourishment, we simulated only one erodible layer, which was characterized by the homogeneous granulometry of the material.  (10) and (11) have been computed as site specific depending on the surrounding bathymetry, the prevailing environmental conditions and the configuration of the breakwater, the results can have a general validity since, in the area of concern, there are no other structures that can significantly affect the outcomes. The ranges of wave damping we simulated fall within the range of variability provided in the literature [73,74], both in reference to the wave height and to the wave power, showing that the breakwater efficiency is a function of the relative submergence of the crest.

Application at Calabaia Beach
We finally applied simplified 2D modeling to assess the effectiveness of the submerged breakwater located at Calabaia Beach through the evaluation of the erosion process affecting the shoreline. Using the same boundary conditions illustrated in Section 2.4, we used the model to reproduce the erosion process affecting the nourishment at Calabaia Beach ( Figure 13). Since the area has been subject to a large volume of nourishment, we simulated only one erodible layer, which was characterized by the homogeneous granulometry of the material. We simulated input waves with a return period of 200 and 50 years, respectively, and two ordinary events with monthly and daily frequencies (Table 3): The analysis shows a significant loss of material from the nourishment, with a bedload from the We simulated input waves with a return period of 200 and 50 years, respectively, and two ordinary events with monthly and daily frequencies ( Table 3): The analysis shows a significant loss of material from the nourishment, with a bedload from the protected cell to an offshore direction ( Figure 14).  In order to evaluate the magnitude of the erosion, we performed a quantitative analysis of the material eroded based on the change in the elevation of each element of the grid located within the nourishment zone, multiplied by the value of its area (12).
The aggregated data account for the erosion process within the event ( Table 4). The results indicate that, with the submerged breakwater, 6 months of ordinary wave action produces a volume loss of 700 m 3 , comparable to a single event with a return time of 200 years, which produces a loss of 629 m 3 . Submerged barriers are generally not suitable to counteract the erosive process produced by ordinary events, requiring complementary measures. Depending on the In order to evaluate the magnitude of the erosion, we performed a quantitative analysis of the material eroded based on the change in the elevation of each element of the grid located within the nourishment zone, multiplied by the value of its area (12).
The aggregated data account for the erosion process within the event (Table 4). The results indicate that, with the submerged breakwater, 6 months of ordinary wave action produces a volume loss of 700 m 3 , comparable to a single event with a return time of 200 years, which produces a loss of 629 m 3 . Submerged barriers are generally not suitable to counteract the erosive process produced by ordinary events, requiring complementary measures. Depending on the duration and extent of the event, nourishment is largely eroded without the protection of breakwater, with a loss of material six to 40 times larger than in the protected configuration.

Conclusions
This study focuses on protected nourishment, an environmentally friendly sea defence intervention. Protected nourishments prove to be effective against extreme phenomena, without providing a valid support to face the regular action of waves, which can be very detrimental to the shoreline, subtracting huge amounts of material where the nourishment has been designed. In this work, we develop a simplified methodology to address the effects of submerged breakwaters on the reduction in energy and power of the entering wave, at a significantly smaller computational cost if compared to full 2D modeling. We apply this method to evaluate the shoreline evolution at Calabaia Beach, comparing the effectiveness of the local submerged barrier with the different characteristics of the waves.
The need for such simplified methods arises from the large number of coastal interventions that are planned worldwide and from the need to protect them from the consequences of mean sea level rises and changes in the frequency driven by climate change. The proposed methodology can be used for analysing and predicting the shoreline development over a period of decades, proving to be an essential tool in the modeling and design of different submerged sea defences. Moreover, it can be used to achieve an initial assessment of the long-term efficiency of former submerged breakwaters and the effectiveness of further interventions, e.g., the use of a larger diameter of nourishment sand or the implementation of nature-based solutions, such as the planting of native seagrass meadows (e.g., Posidonia oceanica). However, this method cannot replace the use of the full 2D or 3D mathematical models, since several phenomena, such as local turbulence, flocculation, benthic suspensions and entrainment of fluid mud by shear flow, usually require a two or three-dimensional approach to assess their effect on bedload and erosion. In conclusion, the application of the presented methodology provides the first evaluation of the effectiveness of the former sea defences built at Calabaia Beach, and proves to be promising in assisting engineers and/or environmentalists in designing/evaluating the efficiency of coastal interventions not only in this specific area, but also in worldwide coastal areas where protected nourishments have been designed or already implemented.