Monitoring and Modelling Coastal Vulnerability and Mitigation Proposal for an Archaeological Site (Kaulonia, Southern Italy)

Abstract

This paper presents a Coastal Vulnerability Assessment (CVA) of a microtidal beach located on the Ionian Sea in Calabria region (southern Italy) in order to examine the influence of the different run-up equations on CVA score and propose mitigation measures for the most vulnerable parts of the beach. The coastal area has been severely eroded by extreme wave storms, which have also damaged important archaeological structures located on a nearby cliff. A typical 1 year return period (T_r) storm, associated with the recent criticalities, was chosen to test the different run-up formulas (Holman (1986), Mase (1989) Stockdon et al. (2006) and Poate et al. (2016)) on a number of beach profiles in order to check the sensitivity of the CVA calculation with regard to the different run-up equations. The obtained results provide evidence that different run-up levels often give rise to different CVA scores. Based on vulnerability results, some mitigation measures have been proposed for the beach in front of the archaeological area, based on submerged detached breakwater and an adherent gabion wall for the cliff defence.

1. Introduction

The complex multi-dimensional behaviour of coastal systems to a wide range of stresses and perturbations requires a strategic approach for the management of the coastal environment. Continual flooding, coastal erosion and loss of livelihood of coastal communities appear in fact to be increasing in intensity given accelerated global climate change [1].

The Coastal Vulnerability Index (CVI), introduced by Gornitz [1], was based on seven physical variables, associated with inundation and erosion hazard related in a quantitative manner. This index, suitable for regional studies, was used by a number of researcher both in Mediterranean context, as De Leo et al. [2] and Besio et al. [3] and in ocean coastal sites as in Di Paola et al. [4]. Moreover, with the addition of demographic variables, the CVI was transformed in the Social Vulnerability Index (SVI) [5,6].

A more specific coastal hazard index, which needs local information in the Coastal Vulnerability Assessment (CVA) [4], was used in this paper. We applied this index for examining strategies that can reduce the vulnerability, particularly in severely damaged coastal areas. The mechanism of wave attack on the beach is schematized with two contributions, the first linked with coastal inundation (and relative wave run-up evaluation) and the second associated with the short-term beach erosion. These contributions have already been used to evaluate the coastal vulnerability and risk [7], due to wave storms, sea level rise [8] and subsidence [9].

In CVA, the wave run-up, defined as the time-varying location of the shoreward edge of water on the beach face, is statistically processed, associated with its maximum length X_max, and normalized through the beach width L to give a beach inundation index I_Ru. The sensitivity of the CVA index with respect to the different run-up equations on several profiles of a micro-tidal beach located on the Ionian Sea in the Calabria region (southern Italy), was examined. Moreover, based on vulnerability results, a mitigation strategy for the archaeological area above a nearby cliff, subjected to severe wave attack during several storm events occurring in 2013 and 2014, was conceived calculating the lowering of the CVA with a proper index recently introduced [10].

In the following, we examine the main equations for run-up evaluation, which is a major parameter for assessing coastal vulnerability. In most equations, the run-up is statistically parametrized in terms of exeedance percentage Ru_2% normalized by H₀ deep water significant wave height, and linked with Iribarren number [11]

$${\xi }_0=\frac{tan\beta }{\sqrt{\frac{H_0}{L_0}}},$$

(1)

where $\beta$ is the beach slope and L₀ is the deep water wave length. Low Iribarren numbers (${\xi }_0$ < 1) indicate dissipative and intermediate conditions, while higher values (${\xi }_0$ > 1) indicate reflective conditions.

Holman [12] proposed an empirical formula to obtain Ru_2%, based on Iribarren number $\xi$_f constrained with the foreshore slope $\beta$_f.

$$\frac{R{u}_{2\%}}{H_0}=0.83{\xi }_f+0.2,$$

(2)

in which the first term is related to the infra-gravity component while the second one is linked to the swell component.

Similarly, Stockdon et al. [13] expressed Ru_2% as a function of wave setup and swash, as given in Equation (3):

$$R{u}_{2\%}=1.1\left(0.35tan{\beta }_f{\left(H_0{L}_0\right)}^{0.5}+\frac{{\left[H_0{L}_0\left(0.563tan{\beta }_f^2+0.004\right)\right]}^{0.5}}{2}\right),$$

(3)

where $\beta$_f is the foreshore slope.

Another empirical formula was proposed by Mase [14] who obtained the run-up level Ru_2% as a function of the empirical coefficients a and b, set as a = 1.89 and b = 0.71:

$$\frac{R{u}_{x\%}}{H_0}=a{\xi }^b.$$

(4)

More recently, Poate et al. [15] introduced a new simplified formula that depends on both local significant wave height H_s, the total beach slope $\beta$ and mean wave period T_m

$$R{u}_{2\%}=Ctan{\beta }^{0.5}{T}_m{H}_s,$$

(5)

where C is a constant fixed to 0.49 by Poate et al. [15].

In order to investigate the CVA sensitivity to these wave run-up formulas, Equations (2)–(5) are examined and discussed with regard to their influence on the whole CVA calculation. Then, a mitigation strategy that could reduce the high vulnerability in the area occupied by the archaeological site was proposed, taking into account the structure typology and position.

The analysis started with the wave simulations based on the WaveWatch III (WWIII) wave model, implemented by Campania Center for Marine and Atmospheric Monitoring and Modelling (CCMMMA)—University of Naples “Parthenope” [16]. This model, in use in southern Italy since 2006 [17,18,19], was first verified in offshore conditions by in situ measurements provided by the Italian Sea Wave Measurement Network buoy of Crotone [20] and then implemented for the tested micro-tidal beach. The run-up levels associated with a typical one-year return period storm were used in the standard CVA evaluation, as conceived by Di Paola et al. [21] and Benassai et al. [22]. To this aim, the offshore wave simulations were first validated with field measurements [23,24] and then used as initial conditions for run-up calculations.

The paper is organized as follows: the study area is presented in Section 2, the methodology is reported in Section 3, results are given in Section 4, and the adaptation strategy is reported in Section 5. Finally, the discussion and conclusions are given in Section 6 and Section 7, respectively.

2. Study Area and Wave Climate

The study area coincides with a ca. 3.5 km long beach stretch, located along the southern Ionian coast of Calabria in the municipality of Monasterace Marina (province of Reggio Calabria, Italy), bordered by the Assi and Stilaro torrents (Figure 1a).

This coastal area is part of the sector Punta Stilo—Capo Spartivento [25], which forms a slightly curved, low-gradient, wave-dominated coast associated with micro-tidal ranges. The sandy to mixed sand-gravel beaches are fed by numerous so-called Calabrian fiumare such as the Assi and Stilaro torrents, i.e., ephemeral streams with typical braided riverbeds that drain steep and rather short hilly slopes, with a hydrological regime closely dependent on rainfall distribution [26].

The longitudinal sediment transport is very variable and constrained by the prevailing winds. In more detail, drift is from north to south during the storms coming from the east and southeast, while it changes direction when storms come from the south. However, the prevailing drift is from north to south, with a quantity of mobilized material equal to 70.000 m³/yr [27].

The well-known archaeological area of Kaulonia, a Magna Graecia colony founded in the 8th century BC and inhabited up until the 2nd century BC, lies just landward of the beach. The remnants of Kaulonia are located on a terraced coastal plain (Figure 1b), overhanging the beach and separated from it by a steep ~5 m high sea-cliff (Figure 1b,c).

The coast is open to waves coming from east and south-east (40° N–220° N), while it is protected from the northeastern waves. The main fetches are found between the directions 40° N and 220° N (in the following S2), while the secondary fetches are found between the directions 0° N and 40° N (in the following S1) (Figure 2a,b).

The offshore wave climate was determined through a statistical analysis of the wave parameters recorded by the Crotone buoy (39°01′06″ N; 17°13′18″ E—Italian Sea Measurement Network—http://dati.isprambiente.it) [20]. Wave data refer to the period January 1990–December 2013 and comprise a total of 115.651 wave records, each characterized by H_s, T_m and D_m. Figure 2c provides the entire data set of pairs of H_s, D_m, divided into secondary fetch S1 and main fetch S2, subsequently subdivided into three sub-sectors: 40–110° N (S2a), 110–160° N (S2b) and 160–220° N (S2c). Each sub-sector is characterized by different wave conditions: H_s lower than 5.5 m for S2a and S2c and H_s higher than 5.5 m for S2b (Figure 2). The historical series of H_s was elaborated to obtain a representative sample of extreme wave events, based on of the Peak Over Threshold method [28], as prescribed by Mathiesen et al. [29].

Table 1. Extreme wave heights corresponding to different return periods for each homogeneous sector calculated using the Gumbel (Fisher–Tippett I) cumulative probability distribution function [30,31]. H_Tr is the wave height for different return period T_r and y_r is the relative reduced variable.

Tr (yr)	Sector	yr	HTr (m)	Tr (yr)	Sector	yr	HTr (m)
1	S2a	3.32	5.00	20	S2a	6.33	7.50
	S2b	3.77	5.95		S2b	6.78	8.74
	S2c	3.97	3.37		S2c	6.98	4.58
5	S2a	4.94	6.35	50	S2a	7.25	8.26
	S2b	5.39	7.45		S2b	7.69	9.59
	S2c	5.59	4.02		S2c	7.89	4.95
10	S2a	5.64	6.93	100	S2a	7.94	8.84
	S2b	6.09	8.10		S2b	8.39	10.23
	S2c	6.29	4.30		S2c	8.59	5.23

provides the offshore wave heights obtained for each part of the main sector, according to the relative return period T_r, which are maximum for the external directions (sector S2b), with H_s > 5.5 m.

In order to evaluate wave run-up levels associated with frequent wave conditions, a recent storm was selected (February 2014) from eastern directions with significant wave heights occurring at least once a year (i.e., a return period of one year).

For this storm, the numerical wave simulations have been performed using a WWIII model; the wave condition has been a starter point to model run-up calculations using the empirical Formulas (2)–(5).

Archaeological Aspects

The study area, owing to the presence of the Archaeological Park of the Magna Grecia colony of Kaulonia, has a strong historical and cultural identity. The ancient city of Kaulonia was founded ca. 700 years B.C. and frequented until the 2nd century B.C.. Its ruins are located on a soil covered, terraced alluvial plain that is separated from the beach by a cliff, up to 5 m high. Kaulonia became a subject of great interest and systematic archaeological activities starting with the excavation campaigns from 1911 to 1916 and the discovery of the basement of the Doric temple, which represents the major archaeological structure of Kaulonia.

One of the most important subsequent archaeological discoveries dates to the end of the 1960s and consists of a mosaic of remarkable beauty representing the famous Dragon of Kaulonia. Since the 1980s until today, numerous excavations and surveys proceeded successfully [32]. In particular, the very recent excavation campaigns in 2012–2014 added very important archaeological findings. Among them, the largest and most articulated Hellenistic mosaic ever found in southern Italy, and the bronze tablet reporting the longest Achaean inscription ever discovered in Magna Graecia.

3. Methodology

3.1. Topo-Bathymetric and Sediment Survey

In order to characterize the beach morphology, topographic measurements were carried out in July 2014 with a DGPS (Differential Global Positioning Systems) system (Trimble R6), characterized by very high accuracy (Root Mean Square Error − RMSE = ∓0.3 m), which consisted of: (i) shore-normal beach profile surveys, extending from the dune crest/retaining wall down to 1 m water depth, i.e., approximately 4 m distance offshore; and (ii) shoreline delineation. Obtained topographic data refer to 10 beach profiles, from north to south P1–P10 (for location, see Figure 1a and Figure 3a). The shoreline position was calibrated taking into account the micro-tidal excursion determined by ISPRA [33]. A bathymetric survey was performed after the February 2014 wave storm, in order to obtain data also on the submerged beach. A digital single beam gauge coupled with a GPS receiver was used (differential mode DGPS post processing) with an acquisition frequency of 0.5 Hz, obtaining 15 shore-normal transects in calm sea conditions reaching a maximum depth of 12 m.

Grain size analysis was carried out on samples collected in July 2014 by applying the ASTM international standard D421 and D422 [34,35] and the Wentworth scale. The medium to long-term shoreline evolution was investigated considering two periods: 1958–2014 and 1996–2014 in order to provide a view of the overall shoreline evolution since the 1950s as well as of the more recent shoreline changes that have occurred during approximately the last 20 years.

3.2. Numerical Wave Model

The simulation of the offshore wave characteristics has been done with the working scheme [36] implemented by the Campania Center for Marine and Atmospheric Monitoring and Modelling (CCMMMA, [37]) of the University of Naples “Parthenope” [38,39]. The operational configuration is based on the Weather Research and Forecasting (WRF) numerical model [40,41], which gives the atmospheric forcing (10 m wind field) needed to estimate the offshore waves. The wave simulations were carried out using the National Oceanic and Atmospheric Administration (NOAA)/National Centers for Environmental Prediction (NCEP) third-generation wave model WaveWatch III (WWIII) [42,43].

The configured model spatial domain has an extension in latitude from 6.33° N to 20.88° N, and in longitude from 36.42° E to 46.98° E. The model outputs include, among others, gridded fields of significant wave height (H_s), wave direction (D_m), and mean wave period (T_m). The grid points of the WWIII model near the coastline were used as "virtual buoy" (VB) to calculate the wave coastal transformation, with the ultimate goal of simulating the run-up on the beach.

3.3. CVA Method

The coastal vulnerability assessment (CVA) is based on the methodology firstly proposed by Di Paola et al. [4], who calculated this index according to the following equation:

$$CVA=I_{Ru}+I_R+I_D+E+T,$$

(6)

where I_Ru is an index associated with the wave run-up distance, I_R is the maximum beach erosion potential index for the shoreline, I_D is the backshore coastal protection structures stability index, E is the multi-year beach erosion rate index and T is the tidal range that is negligible in a micro-tidal environment.

The wave run-up height index I_Ru evaluates the beach potential inundation during wave storms through the run-up distance X_max, calculated as the horizontal distance that the wave travels during the run-up process, given by the run-up level divided by the foreshore slope ($\beta$_f):

$$X_{max}=\frac{R{u}_{2\%}}{tan{\beta }_f}.$$

(7)

I_Ru is obtained from X_max, normalized to the beach width (L). The index of potential erosion I_R provides a measure of the potential retreat of the beach. The I_R values depend on the percentage associated with the maximum potential retreat (R_max) normalized to the width of the beach L. R_max is evaluated as the maximum value of the general solution associated with the convolution method of Kriebel and Dean [44]:

$$\frac{R\left(t\right))}{R_{\infty }}=\frac{1}{2}\left\{1-\frac{{\gamma }^2}{1+{\gamma }^2}exp\left(-\frac{2\sigma t}{\gamma }\right)-\frac{1}{1+{\gamma }^2}\left[cos\left(2\sigma t\right)+\gamma sin\left(2\sigma t\right)\right]\right\},$$

(8)

$$R_{\infty }=S\frac{W_b-\frac{d_b}{m_0}}{B+d_b-\frac{S}{2}},$$

(9)

where $\gamma =2\pi {\mathrm{T}}_{\mathrm{S}}/{\mathrm{T}}_{\mathrm{D}}$, i.e., the ratio between the time scale of beach erosion T_S and the storm duration T_D. In Equations (8) and (9), S is the sea level increase due to wave storm; B is the berm height; ${\beta }_f$ is the foreshore slope; d_b is the breaking depth and W_b is the offshore breaking depth distance.

The backshore coastal protection structure index (I_D) takes into account the coastal structures, whose efficiency is linked to their transmission coefficient K_t, according to

Table 2. Scores associated with each variable for CVA (Coastal Vulnerability Assessment) evaluation.

Parameter	Classification
	1	2	3	4
IR(%)	<15	15–30	30–50	>50
IRu(%)	<40	40–60	60–80	>80
E (m/yr)	<−0.5	−0.5–−1.0	−1.0–−2.0	>−2.0
ID	−1	−2	−3	−4
Kt	>0.58	0.41–0.58	0.24–0.41	<0.24

. K_t was calculated through the following equation [45]:

$$K_t=\left(0.031\frac{H_s}{D_{n50}}-0.24\right)\frac{R_c}{D_{n50}}+b,$$

(10)

where D_n50 is the median nominal diameter of rocks for design conditions, R_c is the freeboard, negative for submerged breakwater. The parameter b is given by:

$$b=-5.42s_{op}+0.0323\frac{H_s}{D_{n50}}-0.0017{\left(\frac{W_t}{D_{n50}}\right)}^{1.84}+0.51,$$

(11)

where W_t is the width of the top of structure and S_op is the deep-water wave steepness corresponding to peak period T_p:

$$S_{op}=2\pi \frac{H_s}{g{T}_p^2}.$$

(12)

The I_D ranking is reported in Table 2, in which I_D assumes negative values, according to Martinelli et al. [10], reducing the overall vulnerability.

Finally, the multi-year beach erosion rate index E, which takes into account the medium-term (20 years) coastline erosion, was derived from the shoreline variations calculated for period 1996–2015.

The resulting CVA index, evaluated according to Equation (6), was divided into four categories, according to the classification defined by Di Paola et al. [4], coming up to the final ranking: very low (CVA = 3), low (3 < CVA ≤ 6), moderate (6 < CVA ≤ 9) and high (9 < CVA ≤ 12) coastal vulnerability (Table 2).

4. Results

4.1. Beach Characterization

The study beach is characterized by a relatively smooth and regular morphology (Figure 3). The landward limit of the backshore is alternatively represented by a few meters high embryonal dune (P1, P2 and P10, Figure 3), an erosion scarp/seacliff increasing in height from north to south up to max 5–6 m (P3–P4, Figure 3), and a seawall of variable height, up to 5 m (P5–P9), limited seawards by boulder barriers in order to improve its resistance to wave-attack (see for instance P7 and P8 in Figure 3). Both the backshore and foreshore zones are composed of mixed-textured, coarse-grained materials. Internal and external areas of the backshore are composed of coarse/very coarse sands and gravels, (see P1 and P3, Figure 3), highlighting a trend to the coarsening of beach sediments from the landward limit of the backshore towards the summer berm. The foreshore is made of sandy gravels with the exception of the areas around the mouths of the Assi and Stilaro streams that are characterized by gravels and gravelly sands, respectively. Median sediment diameters of foreshore sediments (

Table 3. Grain size and morphological features of the backshore and foreshore zones along profiles P1–P10.

Profile	Backshore Width L (m)	Backshore Slope βb (%)	Foreshore Slope βf (%)	Total Slope $\mathit{\beta }$ (%)	Foreshore Median Sediment Diameter $\mathit{\mu }$ (mm)
P1	36.8	6.1	18.4	7.5	5.18
P2	27.6	9.6	14.2	10.6	6.12
P3	31.8	8.6	17.6	10.2	5.96
P4	28.5	9.7	19.9	11.2	3.45
P5	13.3	10.5	17.1	13.6	2.55
P6	39.4	8.9	20.0	10.1	4.66
P7	17.6	5.7	19.1	7.9	5.93
P8	10.0	2.6	18.2	9.5	4.70
P9	20.3	6.1	19.5	8.6	2.80
P10	29.3	3.1	19.3	5.7	2.62

) range between 2.55 mm and 6.12 mm (P2 and P5 in Table 3). Slope gradients of the backshore $\beta$_b and foreshore $\beta$_f range between 2.6–10.5%, and 14.2–20.0%, respectively.

Based on differences in morphology, grain size and anthropic features, the study beach was sub-divided into the following three stretches.

A northern, 0.6 km beach stretch including profiles P1 and P2, which is 25 to 40 m wide and characterized by a low anthropic impact due to the presence of some tourism structures.

A central 0.75 km beach stretch, with an approximate width of 30–35 m. It includes profiles P3, P4 and the archaeological area and is made up of coarse to very coarse sands in the backshore zone, and of sands and gravels in the foreshore zone. Active dunes are completely lacking and the backshore is limited landwards by the erosion scarp/seacliff with the archaeological area. The cliff is made of sandy-gravelly to silty sediments, rich in archaeological material both sparse and organized in several levels. The beach portion in front of the archaeological area was nearly entirely eroded during the storm event on 2 February 2014 and storm waves attacked the cliff directly, causing its local collapse together with some archaeological structures. Because of the extensive damage, a first gabion wall was built in autumn 2014 immediately south to P4 in order to prevent future cliff collapse and avoid further damage to the archaeological site.

Finally, a southern 2 km long beach stretch that extends approximately between P5 and the Stilaro stream mouth (P10) and is characterized by a a high anthropic load mostly due to dense urbanization. The beach, composed of mixed sand and gravel, is 10 to 40 m wide, therefore locally rather narrow (see especially P8, Figure 3 and Table 3). In the past, some coastal defences were realized along this coastal stretch (P7–P8), whose influence has been taken into account through the I_D parameter (see Section 4.3.2).

4.2. Medium and Long-Term Shoreline Evolution

In order to evaluate the E parameter in Equation (6), we analyzed the long-term (1958–2014) and the medium-term (1996–2014) shoreline evolution using aerial photographs in scale 1:39.000 (IGMI, 1958), with an accuracy of RMSE = ∓3 m, the regional technical map in scale 1:10.000 (CTR, 1996—RMSE = ∓2 m) and the topographic data provided by DGPS (RMSE = ∓0.3 m) and aerial drone surveys carried out in 2014 (RMSE = ∓1 m). The shoreline position was assumed coincident with the high-water line position that represents the landward limit of wave and water level influences in micro-tidal conditions [46]. Given that the tidal elevation at surveying time was unknown, we assumed an uncertainty between 2 m and 2.8 m for the daily water line position, considering both the minimum (P2) and maximum (P6) beach foreshore slope (Table 3) together with the local daily tidal excursion of about ±20 cm [33]. Moreover, the 2014 shoreline was delineated by interpolating the points obtained with the DGPS surveys with a spatial resolution of 2 m.

Shoreline changes were determined by the overlay of digitized shorelines. Rates of shoreline change were calculated using an extension of ArcGis release 10 (ESRI©), the Digital Shoreline Analysis System (DSAS version 4.0, [47]). This software allowed the creation of 339 transects, spaced 10 m apart, along which the positions of the three shorelines dating to 1958, 1996 and 2014, respectively, were delineated.

With regard to the long-term shoreline evolution, most of the shoreline has experienced a more or less marked erosion (Figure 4a), except for the central portion, around transects P6 and P7 (respectively −6.61 and −5.97 m,

Table 4. Total and mean annual shoreline variations along profiles P1–P10 during the periods 1958–2014 and 1996–2014.

Period	Shoreline Variation	P1	P2	P3	P4	P5	P6	P7	P8	P9	P10
1958–2014	Total (m)	−25.43	−8.49	−38.46	−38.17	−24.04	−6.61	−5.97	−43.25	−60.67	−54.86
	Annual mean (m/yr)	−0.43	−0.65	−0.65	−0.65	−0.41	−0.28	−0.10	−0.73	−1.03	−0.93
1996–2014	Total (m)	6.90	−0.89	7.74	−33.26	−9.36	16.46	12.12	23.32	27.19	−9.96
	Annual mean (m/yr)	0.37	−0.05	−0.42	−1.80	−0.51	0.89	0.65	1.26	1.47	−0.53

). Severe erosion, instead, occurred within the southern portion (P8–P10, 40 to 60 m retreat) and within the central-northern portion with values up to nearly 40 m in front of the archaeological site of Kaulonia (P3–P4).

In the medium-term (period 1996–2014, Figure 4b and Table 4), the shoreline underwent prevailing progradation (up to 20 m and more, see P8 and P9), whereas some parts of the beach registered again a more or less strong retreat, especially those located in the central-northern part (up to ~33 m, see P4). These results show that the beach in front of the archaeological site has undergone a strong destabilization just over the last 20 years that is still ongoing, as demonstrated by the storm events of 2013 and 2014.

The positive result in terms of shoreline stability of the southern coastal portion, made evident by profiles P6–P9, is most likely related to the realization in 2008 of some coastal protection structures, which have been taken into account in the I_D index (see Section 4.3.2).

4.3. CVA Evaluation Based on the February 2014 Storm

The WWIII model was run to simulate the wave conditions during the sea storm of February 2014, validated with Crotone buoy wave measurements. Results are reported in Figure 5a,b (simulations) and Figure 5c (validation). The agreement between the two time series is fairly good; however, non-negligible differences can be noted, as reported in

Table 5. Statistical indices, obtained according with Mentaschi et al. [48] and Melby et al. [49], for comparison between simulated and recorded wave of the February 2014 storm (Figure 5c). R = correlation index; BI = normalized bias; SI = normalized scatter index; RMSE = root mean square error; NRMSE = normalized root mean square error; BI_P = BI performance; SI_P = SI performance; NRMSE_P = NRMSE Performance; SPS = Summary Performance Score.

	R	BI	SI	RMSE	NRMSE	BIP	SIP	NRMSEP	SPS
Hs	0.933	0.053	0.135	0.496	0.144	0.947	0.865	0.856	0.889
Tm	0.672	0.220	0.186	1.706	0.283	0.779	0.815	0.717	0.771
Dm	0.826	−0.042	0.105	15.965	0.113	0.958	0.896	0.888	0.914

, which lists the statistical indices introduced by Mentaschi et al. [48] and Melby et al. [49]. It can be noted that the simulated H_s (black line in Figure 5c) is much smoother the one measured by the buoy due to the WWIII incapability to resolve the small scales seen in the buoy records. The correlation coefficient R [48] is 0.933 for H_s, followed by R = 0.826 for D_m and R = 0.672 for T_m, while the Summary Performance Score SPS [49] gives the best agreement for D_m parameter (SPS = 0.914), followed by H_s (SPS = 0.889) and T_m (SPS = 0.771). These results support the reliability of the WWIII simulations, which allow us to use them as a starting point in the run-up model calculation.

4.3.1. Wave Run-Up Levels and Inundation Distances Obtained by Different Equations.

The different behaviour of the run-up Ru_2% formulas, with reference to the case study (February 2014 storm), is reported in Figure 6 as a function of the total beach slope $\beta$ (used in Mase [14] and Poate et al. [15] equations) and the foreshore slope $\beta$_f (in Holman [12] and Stockdon et al. [13] equations), in which the $\beta$_f parameter has been given as a multiple of $\beta$.

The examination of Figure 6 shows that, for examined wave conditions, the Poate et al. [15] equation is an upper limit for Ru_2% for beaches with a maximum slope of 15%. In this range, the Holman [12] and Mase [14] equations give intermediate and similar results, while the Stockdon et al. [13] equation represents the lower limit. For slopes higher than 15%, Holman [12] and Stockdon et al. [13] represent the higher and lower limit, respectively.

In order to provide evidence of the sensitivity of the run-up level Ru_2% to the different equations, this parameter has been reported in Figure 7 for the different transects on the beach. According to the results of Figure 6, the Poate et al. [15] values are the highest in all transects, followed by the Mase [14] and Holman [12] results for transects P2, P3, P4, P5, P6, P8 and P9, and by the Holman [12] and Mase [14] results for transects P1, P7 and P10. The Stockdon et al. [13] results are the lowest, except in profile P10, where they are higher than Mase [14].

The general trend of similarity between the Holman [12] and Mase [14] results is interrupted in transects P2, P3, P4, P5 and P10 in which the ratio $\beta$_f/$\beta$ is different from the one in Figure 7.

The results in terms of floading distance X_max, normalized with the width of the emerged beach, reported in Figure 8, closely follow those of Ru_2%. In particular, for transects P5, P7, P8 and P9, the beach is completely flooded regardless of which formula is considered, while in transects P1, P6 and P10, the maximum X_max is obtained with the Poate et al. [15] equation and the minimum one with the Stockdon et al. [13] equation. In terms of I_Ru, the use of Stockdon et al. [13] or Poate et al. [15] equation increases the classification by one score in transects P1 and P6 (from 2 to 3), and in transects P2, P4 and P10 (from 3 to 4).

4.3.2. Evaluation of I_R, E and I_D

According to Equations (8) and (9), the short-term erosion index I_R shows a maximum in P5, followed by P6, P9 and P10, due to the limited beach width and smaller grain sizes (Figure 9). Finally, according to the results presented in Table 4, maximum values of E are found in P4, corresponding to the archaeological site.

The I_D parameter, linked to the presence and efficiency of coastal defences, assumes −2 score in transects P4, P7 and P8, according to the K_t values calculated considering Equation (10). The coastal structures in P7 and P8, constructed in 2008, demonstrated their efficiency in favoring the beach stability, while the submerged structure in P4 has only recently been realized (see Section 5).

4.4. CVA Sensitivity to the Different Wave Run-Up Formulas

The calculation of CVA on the examined transects was performed following the classification given in Table 2 and was reported in Figure 10. In agreement with I_Ru results shown in Figure 8, the same value of CVA regardless of the different run-up formulas is reached for the transects P5, P7, P8 and P9. For transects P1, P2, P4, P6 and P10, the use of the Poate et al. [15] equation instead of Stockdon et al. [13] for Ru_2% raises the CVA score by one class, while, for transect P3, the difference is sufficient to change the CVA score by two classes.

5. Adaptation Strategy

Starting from the obtained results on coastal inundation and erosion, we propose an adaptation strategy in order to decrease the vulnerability in front of the archaeological area. The protection is obtained through a longitudinal submerged breakwater (Figure 11), located offshore of transect P4. The breakwater has a submerged berm width of approx. 12 m, a crest height of −0.50 m a.s.l. and an overall longitudinal extension of 460 m.

Moreover, a beach nourishment is proposed between transects P3 and P5 for a total extension of approx. 800 m, with an additional beach width of 15 m and a grain size of 5 mm. In order to promote the stabilization of the beach nourishment and obtain a natural beach recovery, a beach drainage system could be also taken into account [50,51]. The residual energy calculation for the breakwater provides a value of 0.47 for the transmission coefficient K_t, with a consequent I_D index of −2, like the previous submerged structures realized in front of the transects P7 and P8. These structures lower by two units the CVA score in P4, P7 and P8 profiles, based on the cited K_t value (see Table 2). The defense of the cliff on which the archaeological site is located is obtained by means of an adherent structure.

6. Discussion

The behaviour of the different run-up equations show that the Poate et al. [15] and Stockdon et al. [13] equations give the highest and lowest run-up values, respectively, in all transects, while the Mase [14] and Holman [12] run-up values are intermediate.

The sensitivity of CVA to the different run-up equations cannot be noted for the transects characterized by a lower values of beach width and slope (e.g., P5, P7, P8 and P9 in Figure 10) that are always completely inundated. This is the reason why the different performance of the run-up equations can be made most evident for profiles P2, P3 and P4, where the I_Ru trend closely follows the Ru_2% results. The results show that the run-up obtained by Poate et al. [15] and Mase [14] equations are the highest, giving rise to high inundation percentage, followed by the values obtained by Holman [12] and Stockdon et al. [13]. Similar results are found for P1, P6 and P10, with lower I_Ru values due to a greater beach width, while, in the remaining profiles P5, P7, P8 and P9, the influence of the different Ru_2% formulas is masked by lower values of beach width and slope, which lead to inundation of the whole beach regardless of the used formula.

The mitigation proposal for the most vulnerable transect in front of the archaeological area, consisting of submerged breakwater, decreases the CVA score in P4 transect through a negative value of the index I_D, in agreement with Martinelli et al. [10]. The additional defence of the adherent gabion wall facing the cliff, as well as the increase in beach width provided by the beach nourishment in front of transects P3, P4 and P5 (which decrease both indexes I_Ru and I_R) was not taken into account for precautionary reasons.

7. Conclusions

In this paper, the sensitivity of the coastal vulnerability assessment to different wave run-up formulas was examined with reference to a micro-tidal beach in front of an archaeological site. First, we implemented offshore wave simulations based on a WWIII model and validated them by in situ measurements. Then, we configured the run-up calculations with different well known formulas to obtain different run-up levels and examine their influence on CVA calculations. Finally, we proposed a mitigation strategy for the most vulnerable areas, consisting in submerged breakwater for the area facing the archaeological ruins, coupled with an adherent structure for the cliff defence. The following results were obtained:

• The Poate et al. [15] and the Stockdon et al. [13] equations provide a higher and lower limit, respectively, for run-up levels.
• The run-up differences are less evident in CVA evaluation, nevertheless they can lead to different CVA scores, except in the case of limited beach width where the transects are completely inundated, regardless of the used run-up empirical formula.
• The vulnerability mitigation proposal, consisting of a submerged breakwater for the beach protection, decreased the CVA score based on the K_t value.
• The CVA additional decrease due to an adherent gabion wall for the cliff defense and the artificial nourishment placed in front of the Kaulonia site for a longitudinal extension of 800 m was not taken into account for precautionary reasons.