**Maria Gabriella Gaeta 1,\*, Davide Bonaldo 2, Achilleas G. Samaras 3, Sandro Carniel <sup>2</sup> and Renata Archetti 1,2**


Received: 4 July 2018; Accepted: 28 September 2018; Published: 2 October 2018

**Abstract:** This work presents the results of the numerical study implemented for the natural area of Lido di Spina, a touristic site along the Italian coast of the North Adriatic Sea, close to the mouth of River Reno. High-resolution simulations of nearshore dynamics are carried out under climate change conditions estimated for the site. The adopted modeling chain is based on the implementation of multiple-nested, open-source numerical models. More specifically, the coupled wave-2D hydrodynamics runs, using the open-source TELEMAC suite, are forced at the offshore boundary by waves resulting from the wave model (SWAN) simulations for the Adriatic Sea, and sea levels computed following a joint probability analysis approach. The system simulates present-day scenarios, as well as conditions reflecting the high IPCC greenhouse concentration trajectory named RCP8.5 under predicted climate changes. Selection of sea storms directed from SE (Sirocco events) and E–NE (Bora events) is performed together with Gumbel analysis, in order to define ordinary and extreme sea conditions. The numerical results are here presented in terms of local parameters such as wave breaking position, alongshore currents intensity and direction and flooded area, aiming to provide insights on how climate changes may impact hydrodynamics at a site scale. Although the wave energy intensity predicted for Sirocco events is expected to increase only slightly, modifications of the wave dynamics, current patterns, and inland flooding induced by climate changes are expected to be significant for extreme conditions, especially during Sirocco winds, with an increase in the maximum alongshore currents and in the inundated area compared to past conditions.

**Keywords:** climate changes; sea-level rise; TELEMAC; natural beach; flooded area

#### **1. Introduction**

Strong anthropic pressures, together with the effects of a changing climate, are contributing to the recent increase in the vulnerability of coastal areas. This is particularly significant in the case of extreme events, which, even in the case of a possible decrease of the storminess, are expected to be superimposed to an increased sea level, with an overall intensification of the flooding hazard in coastal regions (for the Mediterranean Sea see for instance [1,2]).

The development of efficient tools to accurately represent nearshore wave and current-induced dynamics is essential for today's operational and forecasting applications in coastal zones [3,4]. High-resolution wave and hydrodynamics modeling offer an extensive range of capabilities to support

coastal planning and decision assessment, accounting for typical features at a coastal engineering scale, such as nonlinear processes of wave propagation and interactions between offshore and coastal structures and the inclusion of inshore boundary conditions, such as river run-off.

In recent years, research efforts have focused on the development of methodological frameworks based on advanced numerical modeling [5,6], which can be used to study the effects of future climate change scenarios affecting both the intensity and frequency of storm-surge events, wave climate, currents, sea-level rise, and riverine sediment discharge. Since the above-mentioned phenomena may increase the flood risk for coastal areas, the understanding of their dynamics at coastal scale becomes essential for the design of climate-change resilience protection and, in general, spatial planning activities.

Thus, the development of multipurpose measures mitigating erosion and inundation and increasing coastal defense efficiency requires a challenging prediction of sea forcings variation induced by the estimated effects of climate change.

Recently, regional future scenarios accounting for the Intergovernmental Panel on Climate Change (IPCC) sea-level projections at 2100 characterized by Representative Concentration Pathways equal to +8.5 W/m2 (hereafter RCP8.5) were applied along the Italian peninsula at DTM (Digital Terrain model) scale [7,8], in order to define the most vulnerable coastal areas in the country and to assess coastal mitigation and adaptation strategies in response to climate change [9]. This static GIS approach can be applied to indicate flood-prone areas at regional [10] or national [11] scales, which can potentially become inundated under future sea-level conditions [12]. However, the implementation of coupled wave-hydrodynamics models, in addition to this static approach, may be advisable on a local scale to obtain more accurate results in terms of nearshore dynamics connected to vulnerability levels and to account for processes of wave propagation, littoral drift, and coastal flooding.

The objective of the present paper is the description of an operational strategy for the development of a multiple-nesting system based on the implementation of open-source numerical models at a high resolution. The natural beach of Lido di Spina, in the vicinity of River Reno mouth, along the eastern coast of Northern Italy (Section 2), is selected to simulate nearshore dynamics induced by the projected effects of climate-change scenarios.

In the specific area of study presented here, we recall the efforts carried out in a 3D context [13–15] and in integrated 2D as well [16], the latter representing a good trade-off between computational costs and results. The proposed study aims to investigate on a site scale how climate changes may impact local parameters such as wave breaking position, alongshore currents intensity and direction and flooded area.

Section 3 describes materials and methods adopted in the present study: In Section 3.1, the framework of the multi-model approach is introduced by specifying the modeling chain of the implemented open-source models, namely, SWAN and TELEMAC. The results of two thirty-year numerical sea state simulations of past (1971–2000) and future (2071–2100) scenarios forced by the climate model COSMO-CLM, in a severe emission future scenario (RCP8.5 in [17]), are presented in Section 3.2, under ordinary and extreme sea events. The implementation of the high-resolution coastal model at the natural beach of Lido di Spina is described in Section 3.3. The model results for both simulated scenarios are discussed in terms of wave nearshore dynamics, current patterns, and flooded areas in Section 4. Conclusions derived from this study are presented in Section 5.

#### **2. Description of the Study Area**

Figure 1 shows the area selected as the representative study case: The natural beach of Lido di Spina and adjacent low-lying zones, located along the coast of Emilia Romagna Region, in the Northern Adriatic Sea. The zone was declared a National (Italian) Nature reserve in 1981, as part of the Po Delta Regional Park, and the Sites of European Community Importance under Directive 92/43/EEC. The site was an extremely dynamic area, under erosion in the last decade and exposed to several climate change related effects, e.g., higher occurrence of flooding and beach retreat [18,19].

The modeled site includes, from North to South, the coastal area between the south-side of Porto Garibaldi and the Reno River mouth (Figure 1), and comprises the mouth of Logonovo channel (zone *a*) at Lido degli Estensi, where sediments transported northwards from the Reno outflow are intercepted by the jetty of Porto Garibaldi, the touristic shoreline of Lido di Spina (zone *b*), and the old outlet plain of Gobbino channel (zone *c*) with low-lying plains.

As is common in all the Northern Adriatic Sea, the wave climate in the study site is characterized by storms mainly generated by northeasterly and southeasterly winds, named Bora and Sirocco, respectively [20,21], the latter inducing the highest surge levels [9].

The investigated area is subject both to a high rate of coastal erosion, with a beach retreat of about 3 km over the last 70 years and a total of about 75 ha of land lost, and to high, human-induced land subsidence due to fluid extraction and off-shore gas platforms activity [22,23]. Moreover, the coastal area presents highly urbanized touristic resorts, known for the wide and sandy open beaches and the natural areas of the Po Delta Park.

The recent local management strategies to protect the area from flooding and erosion have been based on frequent beach nourishments, using depositional materials from the north of the study area, and on the construction of defense structures, such as stone revetments (especially in the area northern the mouth of the River Reno, see [18]). Sea banks along the coastal area strongly preserved it from inundation, although these defense structures may turn out to be rather ineffective in case of poor maintenance and frequently occurring storm events have rendered these defense structures ineffective [24]. In addition, important morphological variations at the mouth of the River Reno, as those related to the deconstruction of the final part of its right riverside, partly compromised the inland safety.

**Figure 1.** Wider view of the study area of Lido di Spina (images from Google Earth [25]; privately processed). The plot also shows the location of areas of interest: Logonovo channel at Lido degli Estensi (zone *a*), Lido di Spina (zone *b*), and Gobbino channel (zone *c*).

#### **3. Materials and Methods**

#### *3.1. Adopted Methodology and Numerical Modelling Chain*

The proposed multiple-nesting methodology is based on a progressive down-scaling, from global climate to coastal models, as described in the following.

a. Wind fields in the Adriatic Sea are obtained every 6 h at 8 km horizontal resolution from the Regional Climate Model (RCM) COSMO-CLM [26], a climate version of the operational weather forecast model COSMO-LM [27]. The simulation generating this dataset [28] encompasses the Italian peninsula and its marine regions, describing the climate evolution under the IPCC RCP8.5 scenario [29,30] in the period 1971–2100. The model is initialized by ERA-Interim Reanalysis [31], and the same model provided boundary conditions for the period 1971–2005, and by the coupled atmosphere-ocean general circulation model CMCC-CM [32] in the period 2006–2100. Compared to previous implementations [33], the validation of the COSMO-CLM wind fields based on in-situ observations along the Italian coast [34] exhibits a particular skill in capturing the wind directional distribution. In turn, results on climatological projections suggest a tendency towards an overall wind energy decrease (in the range of 0–10% across the whole Adriatic Sea), with a stronger decrease along the Bora jets patterns.


Table 1 shows the numerical characteristics of each of the different modeling chain level of the proposed nesting approach, i.e., the model extension, the mesh-size range, the simulated period, the type of run, and the forcings inducing atmospheric and sea dynamics.


**Table 1.** Numerical characteristics of the proposed nesting approach for each of the different chain levels, i.e., the model extension, the mesh size range, the simulated period, the type of run, and the forcings inducing atmospheric and sea dynamics.

The level (a) of the modeling chain, although described in the table, is not discussed in the present study, since its implementation and results are largely described and examined in [28]. Its numerical outputs, such as wind fields under different scenarios, are used to force the wave model SWAN, developed by [34] for the Adriatic Sea, namely, the chain level (b).

The numerical results of SWAN offshore the study site are (1) statistically elaborated and compared to the transposed observations and under CTR and SCE scenarios, (2) analyzed in order to identify extreme events and estimate Gumbel distributions for each simulated scenario, and (3) used to force the TELEMAC model, implemented at high resolution for the beach of Lido di Spina. Details on the statistics applied on the SWAN outputs and on the sea state conditions simulated by TELEMAC are given in the following sections.

### *3.2. Past and Projected Future Wave Climate Scenarios*

### 3.2.1. Statistics on Wave Forcings

The 10-year-long data acquired by the regional NAUSICAA buoy, offshore the near coastline of Cesenatico, at a water depth of 10 m (data available at http://www.smr.arpa.emr.it/dext3r), are analytically transposed offshore the study site. Since the buoy data are not measured at waters deep enough to neglect the refraction effect on the wave directions [40]—that are supposed to be actually altered by refraction—the recorded waves at NAUSICAA are already slightly rotated towards the perpendicular to isobaths lines, ranging approximately between 55 and 60◦ N.

The wave rose showing the distribution of the significant wave height and the mean wave directions at the study area and as a result of the transposed observations is presented in Figure 2, showing that:


**Figure 2.** Wave rose at the study area of Lido di Spina from the transposed observations.

For the simulations carried out in this work, SWAN outputs are retrieved offshore the beach of Lido di Spina, at the nearest grid nodes to the offshore boundary of the implemented TELEMAC.

The seasonal variability of the monthly average of the significant wave height for Sirocco and Bora events is reported in Figure 3, both for the 30-year-long numerical simulations by the SWAN model results under CTR and SCE scenarios and for the transposed measurements from the 10-year-long NAUSICAA buoy records.

A slightly local increase in wave energy intensity is forecasted for events from Sirocco and especially during winter, while Bora energy tends to decrease on average. By comparing with CTR runs, the available transposed observations show good agreement in case of Sirocco events, while, during Bora events, the numerical results underestimate the monthly average of the significant wave height in winter and fall seasons, despite a slight overestimation during summer events.

**Figure 3.** Monthly averages of the significant wave height Hs,m offshore of the study area of Lido di Spina that resulted from CTR (1971–2000) and SCE (2071–2100) scenarios performed by SWAN implementation [34] and from NAUSICAA—transposed observations (2007–2017) for Sirocco (**left panel**) and Bora (**right panel**) events.

3.2.2. Extreme Event Analysis: Joint Probability of Significant Wave Height, Period, and Sea-Level Rise

Storms are identified following the procedure proposed by [41] and classified based on their generating wind direction (Sirocco or Bora). For each of the two dominant directions, extreme events are analyzed considering the return period (henceforth referred to as TR) associated with their maximum significant wave height and fitting the modeled storms with a return period greater than 1 year by means of Gumbel distributions.

In comparison with the control scenario, the analysis in Figure 4 reveals a significant increase in the extreme events coming from Sirocco for the future projection, which is also expected to be important for return periods greater than 50 years; on the other hand, extreme Bora events are expected to reduce their occurrence at the study site.

**Figure 4.** Statistics of extreme events for CTR (blue) and SCE (red) runs for Sirocco (**left panel**) and Bora (**right panel**) events. Filled dots show extreme values from simulations; tendency curves following Gumbel distribution are plotted in solid lines.

The joint distributions of the significant wave height Hs and the peak wave period Tp along the offshore boundary of TELEMAC model are obtained from the transposed observations and plotted in Figure 5 for Sirocco (left panel) and Bora (right panel) events.

*Water* **2018**, *10*, 1380

The two variables can be directly correlated for extreme events (i.e., Hs > 2 m) with an exponential law (black lines in Figure 5), giving a correlation coefficient of 0.92 and 0.97 for Sirocco and Bora, respectively.

**Figure 5.** Relation between Hs and Tp for Sirocco (**left panel**) and Bora (**right panel**) events for the NAUSICAA-transposed observations at the study site. Tendency curves for the extreme events are also plotted with solid black lines.

The effects of the peak water level, including high tidal level and storm surge, are included in the numerical simulations of the study site, following the accurate analysis performed by [40]. They applied a joint probability analysis by means of a copula-based approach [42] to the registered significant wave heights at NAUSICAA buoy and to the contemporary peak water levels measured at Porto Corsini, located 20 km south of Lido di Spina.

Thanks to the proposed cumulated density functions for the two contemporary variables, peak water levels with return periods TR equal to 1 year and 25 years are extracted for the present study, under the assumption that extreme sea storms are generated and forced by same meteorological events at the two locations, i.e., Cesenatico and Lido di Spina.

The authors of [1,43] suggested the existence of strong correlation in the Northern Adriatic Sea between winds and maximum sea levels induced by meteorological events, and this aspect is particularly evident during Sirocco events. Under future scenario conditions, their validated simulations showed how inter-annual evolution of maximum sea levels is expected to be essentially stable from 2051 to 2100, as well as how unchanged extreme statistics for events characterized by return period of less than 100 years [44]. Following these future projections, no modifications in the extreme statistics connecting peak sea-water levels and significant wave heights are reasonably assumed; therefore, also for SCE runs, the peak water levels contemporary to extreme significant wave heights are also extrapolated by the authors of [40].

Finally, sea-level rise in the SCE runs is taken into account by uniformly increasing the bottom depth by +0.70 m: This is a basin-scale bulk estimate based on the recent projections provided by [8] for the Northern Adriatic Sea, based on the authors of [17,45] adopting RCP8.5 scenarios of climate change and adjusting them to the projected rates of vertical land movements (isostasy and tectonics).

Table 2 summarizes all the conditions used in this work for the TELEMAC simulations, in terms of sets of significant wave height Hs, peak period Tp, mean wave direction Dir, peak water level PWL (comprising higher tidal level and storm surge), and the expected sea-level rise SLR values.


**Table 2.** Characteristics of the simulated scenarios, with the significant wave height Hs, the peak period Tp, the mean wave direction Dir, the peak water level PWL, and the sea-level rise SLR.

Note: (1) The values of Tp were estimated according to the exponential curves in Figure 5. (2) The joint probability analysis presented in [40] was adopted to estimate the peak water levels, as described in the text.

#### *3.3. Set-up of the Coupled Wave-2D Hydrodynamics Numerical Model*

Approaching the coastal region, waves generated offshore are influenced by shoaling, refraction, and loss of energy either due to bottom friction or wave breaking [46]. To simulate all these physical processes, including wave-induced currents, the present study is carried out using the wave and hydrodynamics models of the TELEMAC-MASCARET suite that is distributed under a General Public License and is available at TELEMAC [47].

The suite comprises finite-element-based solvers to simulate shallow water hydrodynamics and wave propagation, and is able to model inshore water levels and wave spectra under different forcings. The different included modules can simulate wind wave propagation, ground water flows, tracer transport, sediment transport, and morphodynamics.

In the proposed approach, the wave and 2D hydrodynamics modules of TELEMAC are implemented in order to propagate offshore waves and currents, and reproduce nearshore dynamics and flood processes.

TOM module is a third-generation spectral wave model and solves a simplified equation for the spectral-angular density of wave action by means of a finite-element type method, in order to describe wave propagation and dynamics in coastal areas [38]. The model takes into account bathymetric wave breaking, bottom friction, non-linear wave-wave interactions, wind wave generation, and white-capping. The authors of [48] showed that it is also able to represent the spread wave fields induced by the presence of sandbanks in the nearshore region.

TEL2D module solves the 2D shallow water equations (also referred to as Saint-Venant equations [39]), derived by integrating the Reynolds-averaged Navier-Stokes equations over the flow depth. Several options for the horizontal diffusion terms (depth-averaged k-ε, or constant eddy viscosity models) and source terms (atmospheric pressure gradients, Coriolis force, etc.) could be chosen in the model setup, while the numerical discretization adopts classical methods for the advection terms, such as characteristics, and distributive schemes. Recently, the implementation of implicit schemes enabling relaxation of the CFL limitation on time steps made TEL2D applicable to the treatment of tidal flats, ensuring positive water depth and mass conservation without extra limitation of the time step [49,50].

Since using the same horizontal discretization with a series of Delaunay triangular unstructured elements, the two described modules can be directly coupled (two-way coupling) to account for the effects of waves on the mean coastal circulation reproducing wave-induced currents leading to littoral drift. The gradients of the radiation stress induced by waves are computed using the theory of [51] as part of the hydrodynamics equations in TEL2D.

The implementation of the numerical models in TELEMAC suite, after a series of preliminary tests on the optimal edge dimension for the mesh by using the freely-available pre-processing tool Blue Kenue™ [52], results in a variable density unstructured mesh, consisting of two density areas, the characteristic triangular element size being set to 25 m close to the shore and on the inland areas, in order to better capture the presence of sandbanks and nearshore seabed profiles, and to 50 m for the rest of the computational domain; the total number of nodes is equal to 40,350.

LIDAR topo-bathymetrical data collected by Emilia Romagna Region in 2012 with a spatial resolution up to 1 m are interpolated on the implemented mesh as presented in Figure 6, allowing a detailed analysis of wave-induced hydrodynamics and floods under different sea-level rise scenarios.

**Figure 6.** Interpolated bathymetry map in UTM-WGS84 coordinates used for the TELEMAC simulations (**left panel**), as well as the position of the two selected transects A-A and B-B (solid lines). Zoom on the implemented variable density unstructured mesh for the study area ((**right panel**) corresponding to the black box in (**left panel**)).

TOM and TEL2D are properly set up for the studied area based on previous experience on coupled wave-2D hydrodynamics runs for the representation of nearshore processes, as presented in [53,54]. The processes included in the wave model simulations are (a) energy dissipation due to wave breaking according to [55], (b) energy dissipation due to bottom friction according to [56], and (c) nonlinear transfer of energy due to triad (three-wave) interactions according to [57]. No movable seabed, no defense breaching, and no past subsidence-induced movements are assumed in the study.

The simulated conditions as described in Table 2 are performed for the duration necessary to reach the equilibrium states for all the spectral and hydrodynamics variables, i.e., 3 h, and the presented results are extracted at the final time step of the computations.

#### **4. Results and Discussion**

Numerical outputs from TELEMAC runs under CTR and SCE scenarios and for Sirocco and Bora events are presented in the following sections, in terms of wave nearshore dynamics, current patterns, and flooded areas.

#### *4.1. Wave and Hydrodynamics*

Model results for SCE runs described in Table 2 are presented in the following figures. For Sirocco and Bora events, the significant wave height distribution and mean wave direction vectors that resulted from TOM are shown in Figures 7 and 8, respectively. The increment of water levels in SCE runs, including sea-level rise induced by climate changes, leads to wave propagation over the initially dry zone of the beach, especially in case of TR = 25 years, revealing an expected increase of sea flooding events in the future.

**Figure 7.** Sirocco events. Distribution of significant wave height and mean wave direction vectors for SCE—S1 (**left panel**) and SCE—S25 (**right panel**) runs, respectively.

**Figure 8.** Bora events. Distribution of significant wave height and mean wave direction vectors for SCE—B1 (**left panel**) and SCE—B25 (**right panel**) runs, respectively.

Since differences between runs may appear less significant than in reality on a scale of representation of the global model domain; a comparison of the results is performed along linear transects from the offshore computational boundary to the shoreline (as represented in the left panel of Figure 6), where the distance is taken from the original shoreline position at 0.0 m of mean sea level (hereafter MSL).

The first transect A-A is located in the zone *a* (in Figure 1) of the beach of Lido di Spina and is characterized by a natural submerged bar at a 160 m offshore distance from the shoreline, at a water depth of around 1.8 m, and by an emerged beach at +1.4 m level with a 2% slope. The transect B-B is instead located in the low-lying area of the old outlet of Bellocchio channel (zone *c* in Figure 1) and presents a submerged natural sand bar at a distance of 220 m from the shoreline, at a water depth of 3 m. The spatial evolution of the significant wave height along these two transects is reported in Figures 9 and 10 for Sirocco and Bora events, respectively.

**Figure 9.** Sirocco events. Cross-shore distribution of the significant wave height Hs as calculated by TOM along the transect A-A (**left panels**) and transect B-B (**right panels**).

**Figure 10.** Bora events. Cross-shore distribution of the significant wave height Hs as calculated by TOM along the transect A-A (**left panels**) and transect B-B (**right panels**).

Wave dissipation occurs as a result of both breaking and bottom friction. TOM model simulates the evolution of the wave height, frequency, and direction as the waves propagate towards the beach. In the presence of sand bars, the cross-shore distribution of significant wave height usually has a typical evolution, influenced by the composite bathymetry. Indeed, wave steadily propagates towards the breaking point, beyond which it starts to decrease in height up to submerged bar crest. Then, it remains constant until water depth reaches again the bar crest depth; then, it slowly decreases up to reach the shoreline.

Simulations of CTR and SCE scenarios produce evident variation in wave breaking position and height, denoting a shift to higher wave heights closer to the shoreline in case of SCE runs (especially in case of Sirocco events) and a significant shift of the breaking line position towards the shoreline.

The wave radiation stresses give rise to the alongshore velocity profiles, which are responsible for inter-annual sediment transport, as reported in Figures 11 and 12, showing Sirocco and Bora (with opposite ordinate axis) events, respectively, where positive values denote currents directed from South to North. The maximum value of alongshore velocity occurs within the breaking zone of each simulation.

**Figure 11.** Sirocco events. Cross-shore distribution of the alongshore velocity as calculated by TEL2D for ordinary S1 and extreme S25 events and of the water depth along transect A-A (**left panels**) and transect B-B (**right panels**).

**Figure 12.** Bora events. Cross-shore distribution of the alongshore velocity as calculated by TEL2D for ordinary B1 and extreme B25 events, and of the water depth along transect A-A (**left panels**) and transect B-B (**right panels**).

In presence of submerged bars, as in the case of the adopted bathymetry, the distribution of alongshore current usually has two peaks, each corresponding to the peaks of wave breaking-points [58,59], with the maximum on the bar crest and the minimum much closer to the shoreline. In the critical part of the surf zone beneath the alongshore-current jet, the intense wave-current stresses occur, while prior to breaking and in the inner surf zone, the current strength is smaller.

The combined effects of the expected sea-level rise and the intensity increase of Sirocco events bring a modest shift in nearshore wave-induced currents, resulting from the respective change in wave dynamics. This is evident for both the analyzed transects in case of Sirocco events but only for the low-lying area of Bellocchio channel in case of Bora events.

### *4.2. Inundated Areas*

Figures 13 and 14 show the comparative effects of changing climate scenarios on the inundated areas that result from Sirocco and Bora runs, respectively: Black lines indicate the initial shoreline at a water depth equal to ±0 m MSL, while blue and red lines represent the upper limits of the inundated area at the end of CTR and SCE simulations.

**Figure 13.** Sirocco events. Inundated areas for ordinary (**left panel**) and extreme (**right panel**) events and for CTR (blue lines) and SCE (red lines) scenarios. Black lines indicate the initial shoreline.

**Figure 14.** Bora events. Inundated areas for ordinary (**left panel**) and extreme (**right panel**) events and for CTR (blue lines) and SCE (red lines) scenarios. Black lines indicate the initial shoreline.

Differences between CTR and SCE scenarios are the most noticeable ones among the presented results on the scale of the representation; they highlight the utility and effectiveness of the implemented numerical approach in the simulation of the respective phenomena, and they are used to make a comparative analysis.

The areas mainly suffering flooding in the scenarios of climate changes are the northern mouth of Logonovo channel (zone *a* in Figure 1) and the southern low-lying plain of the older outlet of Bellocchio channel (zone *c* in Figure 1), which are expected to be the most vulnerable and exposed to flooding hazards induced by sea-level rise, and during ordinary events.

Indeed, looking at the numerical results for the three different zones in which the studied site is divided in Figure 1:


#### *4.3. Discussion of the Effects of Expected Climate Changes*

Changes in deep-water wave climate drive littoral hydrodynamics and morphological variation, such as the patterns of nearshore wave propagation and alongshore currents over local shelf bathymetry.

A quantitative comparison between the different scenarios and events is proposed in Table 3, in which differences between the simulated CTR and SCE scenarios are reported in terms of characteristics of wave-current-induced dynamics with the aim of estimating the effects of expected climate changes at the study site.


**Table 3.** Computed variation induced by climate changes between SCE and CTR runs of wave breaking line position ΔxB, maximum alongshore velocity ΔUls, MAX, and percentage of inundated area ΔΩ.

In particular, the following features are calculated from the numerical results:

• Variation of the wave breaking line position ΔxB (m) calculated as

$$
\Delta \mathbf{x}\_{\rm B} = \mathbf{x}\_{\rm B, SCE} - \mathbf{x}\_{\rm B, CTR} \tag{1}
$$

in which xB,SCE and xB,CTR are the breaking line positions for SCE and CTR runs, respectively. This data would allow estimating wave energy in the surf zone and impacting the shoreline.

• Changes in the maximum alongshore currents calculated ΔUls, MAX (%) as

$$
\Delta \mathbf{U}\_{\text{ls,MAX}} = (\mathbf{U}\_{\text{ls,MAX}} - \mathbf{U}\_{\text{ls,MAX}\,\text{CTR}}) / \mathbf{U}\_{\text{ls,MAX}\,\text{CTR}} \times 100\tag{2}
$$

in which Uls, MAX SCE and Uls, MAX CTR are the maximum alongshore velocity computed in SCE and CTR runs, respectively. This value could give evidence of the estimated variation induced by climate changes in morphological patterns driving the inter-annual area evolution.

• The computed percentage of inundated area variation ΔΩ (-) is calculated as

$$
\Delta\Omega = (\Omega\_{\rm SCE} - \Omega\_{\rm CTR}) / \Omega\_{\rm CTR} \tag{3}
$$

in which ΩSCE and ΩCTR are the estimated inundated areas for SCE and CTR runs, respectively; this information could help in evaluating the flooding hazard increment connected to climate changes.

According to computed results reported in Table 3, modifications of the wave dynamics, currents patterns, and inland water propagation induced by climate changes are expected to be mainly significant for extreme conditions, and especially for Sirocco events, with increases in the maximum alongshore currents and in the flooded area dimensions. Decreases or negligible variations are instead estimated for Bora events in terms of alongshore currents, again revealing how long-term effects of a changing climate along the North Adriatic shoreline may be stronger related to Sirocco events.

Expected effects of increments in wave height and in sea level in case of Sirocco events during SCE runs (as reported in Table 2) may induce no evident variations in the breaking line position, while wave energy impacting the shoreline and inland area may be strongly affected by wave energy flux, which is related to wave height and storm duration. An onshore shift of the breaking line in case of Bora events is instead connected to a balance between wave height decrement and sea-level rise during SCE runs (as reported in Table 2).

An overall increase of the flooded area extension ranging from 3 to 9 times the CTR values is observed in the results of all the simulated scenarios, showing that the greatest threat from future climate changes, especially during ordinary sea storms generated by Sirocco winds, would be the relative rise of sea level, while slight modifications of the littoral currents would drive morphological dynamics in the shoreline stability and evolution.

#### **5. Conclusions**

Wave-current dynamics represent a feature of paramount importance in the northern Adriatic basin, both for triggering local circulation [13,60] influencing connections in terms of amount of dense water produced [59] and driving sediment re-suspension and mobilization [14,61,62]. Modifications induced by climate changes in wave-current dynamics need to be carefully evaluated, since they can be regarded as among the main threats to coastal protection and stability, with important implications for the assessment of littoral management plans and for mitigation efforts [4,63].

In the present study, coupled wave-2D hydrodynamics simulations carried out by means of the open-source TELEMAC suite were implemented in the coastal area of the River Reno mouth, along the eastern coast of Northern Italy. A multiple-nesting approach was developed by adopting COSMO-CML and SWAN results for the Adriatic Sea in order to run past (1971–2000), and future climate change (2071–2100) scenarios were developed at high-resolution coupled wave-2D-hydrodynamics TELEMAC model.

Presented results showed modifications in the area dynamics as a response to expected climate change scenarios and previously integrated carried out studies, limited to sea level dynamics. A slight local growth in the wave energy intensity for Sirocco events, however, inducing an increase in alongshore currents and inundated areas, is observed. On the other hand, during Bora events, negligible changes in wave-induced coastal dynamics are forecasted, while flooded area dimensions again increased under the simulated sea-level rise effects.

The presented numerical results showed that flooding hazards and changes in littoral hydrodynamics at the selected site are nowadays already significant, especially during Sirocco extreme events [24], and are expected to further increase in the future. It is also found that, according to the basin sea level projections elaborated by IPCC [8,17], the highest contribution to the coastal vulnerability of the studied beach is due to the relative rise of sea level, especially when this is combined with extreme sea storms.

Overall, the described methodology is deemed to be of general interest for ocean and coastal modelers involved in the development of procedures for offline multiple-nested high-resolution simulations from global circulation models, the adaptation of coastal models to platforms of operational oceanography, and the implementation of the above in coastal planning and design actions at the light of climate change scenarios.

Future studies of the selected site will focus on the implementation of a validated morphological model to perform erosional risk analysis and the running of long-term scenarios in order to estimate the shoreline evolution and seasonal variability of coastal dynamics in the area.

**Author Contributions:** For the presented research, D.B. performed the climatologic simulations with SWAN, M.G.G. and D.B. analyzed the wave climate and performed statistics, M.G.G. and A.G.S. performed the numerical simulations set-up and analyzed the results from TELEMAC runs, and R.A. and S.C. conceived the research aim and contributed to the discussion on the study outcomes.

**Funding:** This work was supported by the Flagship Project RITMARE—The Italian Research for the Sea, coordinated by the Italian National Research Council and funded by the Italian Ministry of Education, University, and Research within the National Research Program 2011–2013.

**Acknowledgments:** Geological Unit from Emilia Romagna Region is acknowledged to provide DTM data of the study area. D.B. and S.C. acknowledge the support from the UE H2020 program under grant agreement n. 730030 ("CEASELESS" Project) and the INTERREG MED "CO-EVOLVE" Project. The authors thankfully acknowledge Edoardo Bucchignani (CMCC and CIRA, Capua, Italy) for providing the COSMO-CLM wind fields used in the wave model simulations.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
