The Role of Sea State to the Morphological Changes of Prasonisi Tombolo, Rhodes Island, Greece

Abstract

Wave-induced morphodynamic processes that cause formation, preservation, and destruction of the Prasonisi tombolo in Rhodes Island are investigated, based on satellite image analysis and numerical modeling. A new method is developed for extracting wave events that consist of successive wave data of similar characteristics. The wave events refer either to wind seas or swell seas. This process combined with the satellite image analysis is then utilized for the derivation of the most representative wave scenarios that affect tombolo and salient formation. In particular, the main factors that play a significant role in tombolo and salient evolution are the offshore wave conditions, the location and width of the surf zone, the maximum value of the wave breaking index in the study area, and the initial bottom bathymetry before the study area is exposed to a new sea state. In general, the proposed method provides a realistic insight into tombolo morphodynamics and can be used to provide a cost-effective approach and a wave data-reduction technique for coastal engineering studies.

1. Introduction

A tombolo, also known as a “sandy isthmus”, is a dynamic landform of great geomorphological and socioeconomic interest that is altered by wind, waves, and coastal currents in a continual process of accretion and erosion. A tombolo is defined as a pair of opposing beaches or a sandbar that joins an island with a mainland or another island as a result of longshore drift or the migration of an offshore bar toward the coast [1]. It is a constructive feature occurring above a shoal—a natural submerged (rocky or sedimentary) ridge—which is well protected from large waves due to the presence of islands.

The formation of a tombolo is caused by the wave diffraction at the ends of a natural or artificial obstacle originally detached from the mainland, the subsequent convergence of opposing flows behind the obstacle, and the final deposition of sediment at the area between the obstacle and the mainland. The last stage of this process is the emergence of a sedimentary bar connecting the obstacle and the mainland.

A salient is a dynamic landform developed in response to patterns of diffraction and energy reduction in the lee of an island or a coastal structure [2] and as a result the beach forms a sinusoidal shape in this region. In the case of a salient, the island or coastal structure and the mainland are disconnected; this is the main difference between the states of salient and tombolo. A salient can also be formed by the eroding effects of high-energy waves attacking a tombolo. The transition from tombolo to salient usually has a shorter duration than tombolo formation as a consequence of diffraction and low-energy conditions.

The processes of a tombolo or a salient formation, though detectable, are not so clearly understandable and predictable [3]. This is mainly due to the complex interaction between the underwater and the terrestrial morphology of the tombolo [4], the shape and dimensions of the island and its distance from the mainland [5,6], the availability and texture of the sediments [7], and the sea state [8,9].

Although some tombolos, formed earlier than or during Holocene, remain stable for a long time, in terms of their evolution, such as those located on the coast of Western Australia [4], the Perachora Peninsula, the Eastern Gulf of Corinth, Greece [10], Tyre, Lebanon [11], and Büyüksivri Tepe of Ephesus, Turkey [12], the majority of them are active (ephemeral) features formed, modified, and destructed within a short time period of less than a year to a few centuries [13].

Additionally, based on the existing literature [14,15,16,17,18,19], the main morphometric features that can be measured in tombolos are the tied island’s length (L), the tombolo’s length from the initial shoreline to the island (D), and the smallest tombolo width (b), whereas for salients, b is the width at the inflexion point and D is the distance from the tip of the salient to the initial shoreline. Furthermore, there is a demarcation line between tombolos and salients, referring to the ratio L/D. However, there is a significant scatter in the results provided by many researchers [14,15,16,17,18,19]. An approximate value of this ratio at 1.0 is in accordance with the relevant international literature, i.e., L/D ratios of usually less than 1.0 lead to the formation of a salient, while ratios larger than 1.0 lead to tombolo formation. Referring now to our case study, the Prasonisi island has a ratio of L/D equal to 1.7, which agrees well with the fact that this is a tombolo.

Referring to the touristic, economic, and environmental value of tombolos, the broader area of Prasonisi is of great importance for the local community. The tourist perception for the Mediterranean region is concentrated into three significant factors: the sun, sea, and sand (3S) market. Especially, a high percentage of international tourists visiting Mediterranean countries (Spain, France, Italy, Greece) are interested in beach-related tourist activities [20]; a characteristic also observed at the international level that beaches are a major factor in the tourist decision market. Tombolos consist of sandy sediments [21] and produce two beaches in most cases, placing them between the top priorities for tourism [22]. This has also been justified by different research approaches concerning the suitability of tourism exploitation of tombolos, which showed significant results [23]. Moreover, they produce beautiful scenery, being one of the top considerations of tourists in order to visit a region. Tombolos usually belong to protected areas and this fact acts in many cases positively in the decision-making procedure for tourists [22].

The Prasonisi tombolo is the largest geomorphologically active tombolo in Greece, linking the Prasonisi Isle with the Rhodes Island. It is subject to wave activity from two completely opposing directions, which causes strong short-term morphological changes over the years.

The present study investigates the morphodynamic processes that cause formation, preservation, and destruction of the Prasonisi tombolo. Some remarks derived from the above investigation will be generalized for all tombolos. In particular, the main goals of this paper are: (a) the determination of the sea state characteristics that affect the tombolo’s morphological changes (i.e., tombolo formation, salient formation, tombolo or salient preservation) using satellite image analysis and numerical modeling; (b) the development of a methodology for extracting the most representative wave scenarios that affect tombolo and salient formation; and (c) the comparison of the representative wave scenarios of tombolo and salient formation.

2. The Study Area

The Prasonisi tombolo is located in the southeastern Aegean Sea, Greece, connecting the Prasonisi islet with the southern landmass of Rhodes Island (Figure 1). The six aerial photos [24] (of years 1960, 1971, 1975, 1980, 1990, and 2005) displayed in Figure 2 show the changeable shape of the Prasonisi tombolo’s shoreline.

The surface area of the Prasonisi Islet is about 1.5 km², having a maximum length and width of 2 km and 1.3 km, respectively. It consists of Plio–Pleistocene limestones and its coast is abrupt and rocky. The southernmost coast of Rhodes Island consists of Pleistocene calcareous and Pliocene clastic (conglomerates) formations drained by small ephemeral streams (Figure 1).

There is not any significant and permanent sediment source feeding the wider area of the tombolo. The available solid material has been derived probably from the long-term erosion of the coastal cliffs and the runoff from the southernmost watersheds of Rhodes Island.

The tombolo is a relatively large (c. 750 m in length along its longitudinal axis and c. 850 m in width along its maximum lateral axis) delta-shaped tombolo composed of medium sand. However, its shape is constantly changing, as evidenced by a series of available aerial photographs and satellite imagery since 1960 (Figure 2) and recent photographs (Figure 3). Sometimes, the morphological changes are so severe that the tombolo is cut off from the islet, forming a salient. The salient can be maintained for a period of a few months to a few years, switching afterwards to a tombolo again as a result of coastal morphodynamics. This changeable landform seems to adapt its shape to the dominant wave regime, since the tidal range is very low (a few tens of centimeters) [26]. It is quite possible that during the constructive periods, sediment being transported from the nearshore zone to the beach nourishes and feeds the tombolo, while during the destructive periods, sediment moves in the opposite direction to the sea building or expands submerged longitudinal sandbar.

The mean seafloor slope at the northwestern side of the tombolo is 1%, which is about half of that at the southeastern side (2%). The nearshore zone from both sides of the tombolo consists of medium-to-fine sand, while seagrass meadows are absent at depths shallower than 20 m.

The prevailing winds come from the WNW direction and usually have an intensity of 7–21 knots (for 60% of the winds), while waves usually come from the NW direction, having an average significant wave height of about 1 m. Although the frequency of E, S, and SE winds and waves is very low, these wind and wave conditions are relatively strong [27] (see Figure 4).

Moreover, Prasonisi is currently included in the Natura 2000 network for the protected areas in Greece. Particularly, Prasonisi and the broader area is integrated at the Natura_2000 network as a Special Protection Area (SPA), Site of Community Importance (SCI), and Special Area of Conservation (SAC) with the code GR4210031 and site name “Notio Akro Rodou, Prasonisi, Ygrotopos Livadi Kattavias” since 2005. It is also well known for different tourist activities such as windsurfing, kitesurfing, hotel amenities, beach bars, etc., due to its distinctive morphology and since it is exposed to a variety of wind-wave conditions [26].

3. Materials and Methods

3.1. Used Data

The fieldwork, including the terrestrial and underwater geomorphological mapping of Prasonisi’s tombolo, was conducted during numerous missions lasting from 2014 to 2021. This study was based on geological and neotectonic data provided by E. Mutti et al., [29] and Lekkas et al. [25] and bathymetric data collected from nautical maps of the Hellenic Navy Hydrographic Service.

The morphological changes of the tombolo were based on the analysis of free access Landsat satellite images (TM: Thematic Mapper, ETM+: Enhanced Thematic Mapper Plus and OLI: Operational Land Imager) acquired over the period 1995–2021, with an 8-day or 16-day repeat cycle and a spatial resolution of 15 m or 30 m. The Landsat imagery was downloaded in GeoTIFF format from the United States Geological Survey (USGS) Earth Explorer Website (http://earthexplorer.usgs.gov/ (accessed on 15 February 2021)). In addition, Sentinel-2 Level-2A multispectral images were taken between 2017 and 2021 with a 10-day repeat cycle and a spatial resolution of 10 m. The Sentinel-2 imagery was downloaded from the European Space Agency website (https://scihub.copernicus.eu/ (accessed on 15 February 2021)). After the necessary corrections, the shoreline was automatically extracted from the available satellite imagery after applying the Normalized Difference Water Index (NDWI).

The time series of wind and wave data used for the current analysis were obtained from three datasets: the Wind and Wave Atlas of the Hellenic Seas [28], the Copernicus-MEDSEA_HINDCAST_WAV_006_012 [30], and the Copernicus-MEDSEA_ANALYSISFORECAST_WAV_006_017 [31]. The point of wind and wave data sampling is depicted in Figure 1. The time range and time interval of these three datasets are presented in

Table 1. The available datasets used to extract wave scenarios.

s/n	Datasets	From	To	Δt (h)
1	Wave and Wind Atlas of the Hellenic Seas	1 January 1995	31 December 2004	3
2	Copernicus-MEDSEA_HINDCAST_WAV_006_012	2 February 2006	1 January 2020	1
3	Copernicus-MEDSEA_ANALYSISFORECAST_WAV_006_017	2 January 2020	2 June 2021	1

.

3.2. Extraction of Wave Scenarios

A methodology is developed aiming to extract the characteristic wave scenarios from the available wave and wind datasets (Table 1) that cause significant short-term sea bottom morphological changes in the study area.

Particularly, certain subsets of the available wave and wind data are examined, representing specific types of the tombolo’s morphological state for specific periods of time.

Initially, the swell criterion is applied to the datasets, based on the Thompson et al., [32] classification to discriminate swell seas from wind seas. This criterion is implemented because the two wave types have distinct characteristics [33], which in turn produce different effects upon the sea bottom and the beach profile. Wind seas usually have an offshore peak wave steepness higher than 0.025. This steepness is defined as the ratio of offshore spectral significant wave height (Hm0) and wave length Lp, based on peak wave period Tp. The swell steepness diminishes as the wave moves away from the fetch region [34]. Swell seas with steepness values of 0.010–0.025, 0.004–0.010, and <0.004 are classified as young, mature, and old, respectively.

Since sea bottom morphological changes are mainly caused by wave-induced currents and short-crested waves are generated by wind, wind characteristics and their time history have been used here to extract representative wave scenarios for wind seas. The primary aim in the present analysis is the derivation of wind-wave events that are defined here as consecutive wave data of wind seas corresponding to approximately constant wind. It is noted that these wind-wave events could be of different types, such as fetch-limited or duration-limited growing seas or fully developed seas (e.g., [35,36]). It is also noteworthy that the definition of constant wind is commonly used for the estimation and prediction of offshore wave conditions of fetch-limited or duration-limited growing seas or fully developed seas based on wind data.

Moreover, wind intensity and direction in deep waters in combination with fetch length affect the magnitude of wave characteristics of wind seas that are responsible for sediment transport in shallow waters. Thus, the definition of constant wind is utilized [37] to determine the segment of the consecutive data that correspond to each wind-wave event. In this manner, wind duration at i_th (e.g., hourly) wind data is considered equal to the number of successive hours that satisfy the two conditions [36] displayed below:

$$\left|U_{w,i}-{\overline{U}}_w\right|<2.5 {\mathrm{ms}}^{-1}$$

(1)

$$\left|D_{w,i}-{\overline{D}}_w\right|<{15}^o$$

(2)

where${\overline{U}}_w$ and ${\overline{D}}_w$ are the mean values of successive and acceptable (e.g., hourly) wind speed $U_{w,i}$ and direction $D_{w,i}$ data, respectively.

Since wind direction is a circular variable, and thus has a discontinuity at 360°, the mean value ${\overline{D}}_w$ can be estimated via the use of an appropriate formula based on approximations, e.g., [38,39,40,41]. In a similar rationale, the absolute difference angle between $D_{w,i}$ and ${\overline{D}}_w$ is considered the acute angle between mean wind direction and i_th wind direction.

Following the above two conditions of constant wind, each wind-wave event in the present paper is determined by a segment of non-overlapped consecutive wind data that satisfy the statement, noted below:

$$\forall t\in \left[t_1,t_2\right]: \left|U_w\left(t\right)-{\overline{U}}_w\right|<2.5 {\mathrm{ms}}^{-1} and \left|D_w\left(t\right)-{\overline{D}}_w\right|<{15}^o$$

(3)

where $t$ is the time variable and $\left[t_1,t_2\right]$ is the time interval covered by a certain wind-wave event.

Referring now to linear data, such as wave height and wave period, the representative wave characteristics of each wind-wave event are the arithmetic mean values of the corresponding consecutive wave data as follows:

$${\overline{H}}_{mo}=\frac{{{\displaystyle \sum }}_{i=1}^N{H}_i}{N}, H_i=H_{mo}\left(t\right), t\in \left[t_1,t_2\right]$$

(4)

$${\overline{T}}_p=\frac{{{\displaystyle \sum }}_{i=1}^N{T}_i}{N}, T_i=T_p\left(t\right),t\in \left[t_1,t_2\right]$$

(5)

$$D={\displaystyle \sum }_{i=1}^N\Delta t$$

(6)

where ${\overline{H}}_{mo}$, ${\overline{T}}_p$, are the mean values of the successive data of the significant wave heights $H_i$, peak wave periods $T_i$, respectively, included in each wind-wave event, N is the number of successive data, and $t$ is a discrete variable here for time. In addition, $\Delta t$ is the time resolution of the dataset used (see last column of Table 1) and $D$ is the duration covered by a wind-wave event that satisfies conditions described in Equation (3).

As for the circular variable $\overline{MWD}$, this has been decided to be estimated by the following formula (referenced in, e.g., [38,39,40,41]).

$$\overline{MWD}=\left\{\begin{array}{c}\mathrm{atan}\left(\frac{s_a}{c_a}\right), if s_a>0 and c_a>0 \\ \mathrm{atan}\left(\frac{s_a}{c_a}\right)+\pi , if c_a<0 \\ \mathrm{atan}\left(\frac{s_a}{c_a}\right)+2\pi , if s_a<0 and c_a>0 \end{array}\right.$$

(7)

$$\mathrm{where}: s_a={{\displaystyle \sum }}_{i=1}^N\sin (D_i), D_i=MWD\left(t\right), t\in \left[t_1,t_2\right]$$

(8)

$$c_a={{\displaystyle \sum }}_{i=1}^N\cos (D_i), D_i=MWD\left(t\right), t\in \left[t_1,t_2\right]$$

(9)

where $\overline{MWD}$ is the mean value of the successive data of mean wave direction $D_i$ included in each wind-wave event.

It is noted that an assumption is made in the above methodology, i.e., wind characteristics in the point of data collection are the same as those in the wave-generation area, and that wind characteristics in each time step are responsible for wind-wave characteristics in the same time step. The latter means that the time step is adequately large in order that waves generated by the wind in the source area have traveled and reached the check point, but also wave and wind conditions are stable during the whole time interval.

Nevertheless, the condition of constant wind cannot be applied to swell data, due to the fact that swell waves are generated as a consequence of the dispersion of wind waves from distant weather systems and, consequently, their characteristics have been altered as they move out from their source area. It is noted that another assumption is made in the present methodology, for reasons of simplicity; particularly, that all successive values of swell data between two wind-wave events belong to the same swell event.

Moreover, a minimum duration of 3 h is applied to wind-wave and swell events, eliminating in this manner events of short duration that are probably not able to cause significant short-term morphological sea bottom evolution. As for the characteristics of each representative swell scenarios, these are estimated as the mean values of their successive data, similarly to wind-wave events.

It is noted that this methodology differs from wave climate schematization techniques, e.g., [42,43], that are appropriate for long-term morphological simulations (≥1 yr) aiming to assess the potential sediment transport over a long time period. However, a similar approach to the long-term wave climate schematization techniques could also be performed by using the method mentioned above, via grouping of the wind-wave events and swell events.

Moreover, the methodology developed for the extraction of wave events has a main advantage for the investigation of the short-term sea bottom evolution. In particular, wind-wave events consist of wind and wave successive data of similar wind and wave characteristics; thus, from a statistical point of view they can be represented more properly by their mean values compared with wave events extracted from non-successive data or data presenting a significant variability. In this manner, wave events of a variety of wave intensities can be derived and thus the analysis is not only restricted in extreme wave events (storm events) [44,45] as is usually the case in extreme value analysis. Therefore, either coastal erosion (e.g., [46,47,48]) or accretion [49] as well as beach recovery can well be considered and investigated. In addition, the chronology of wave events could also be used when studying the coastal morphodynamics and sea bottom evolution in a sequence of wave events. However, this has not been studied in the present paper, as a preliminary investigation of a tombolo’s most representative wave conditions is conducted.

A more detailed description of the wind-waves events and swell events extraction method from wind and wave data is presented in Figures S1 and S2. Besides, an application example of the methodology is displayed in Table S1 where the dataset is presented, whereas Figure S3 and

Table 2. Time periods covering specific morphological characteristics of the study area.

Start Date	End Date	State of Morphology	Satellite	Symbol
[	)
10 November 1995	24 March 1996	salient formation	Landsat 4-5 TM	1_sf
24 March 1996	28 June 1996	salient preservation	Landsat 4-5 TM	2_sp
28 June 1996	15 August 1996	tombolo formation	Landsat 4-5 TM	3_tf
15 August 1996	9 January 1998	tombolo preservation	Landsat 4-5 TM	4_tp
9 January 1998	30 March 1998	salient formation	Landsat 4-5 TM	5_sf
30 March 1998	13 August 2001	salient preservation	Landsat 4-5 TM	6_sp
13 August 2001	22 August 2001	tombolo formation	Landsat 4-5 TM	7_tf
22 August 2001	9 August 2002	tombolo preservation	Landsat 4-5 TM	8_tp
9 August 2002	25 June 2003	salient formation	Landsat 4-5 TM	9_sf
25 June 2003	18 June 2004	salient preservation	Landsat 4-5 TM	10_sp
18 June 2004	6 September 2004	tombolo formation	Landsat 4-5 TM	11_tf
6 September 2004	11 December 2004	tombolo preservation	Landsat 4-5 TM	12_tp
11 December 2004	15 January 2005	salient formation	Landsat 4-5 TM	13_sf
15 January 2005	7 July 2005	salient preservation	Landsat 4-5 TM	14_sp
7 July 2005	8 August 2005	tombolo formation	Landsat 4-5 TM	15_tf
8 August 2005	4 February 2010	tombolo preservation	Landsat 4-5 TM	16_tp
4 February 2010	11 February 2010	salient formation	Landsat 4-5 TM	17_sf
11 February 2010	12 June 2010	salient preservation	Landsat 4-5 TM	18_sp
12 June 2010	28 June 2010	tombolo formation	Landsat 4-5 TM	19_tf
28 June 2010	19 November 2010	tombolo preservation	Landsat 4-5 TM	20_tp
19 November 2010	28 December 2010	salient formation	Landsat 4-5 TM	21_sf
28 December 2010	12 October 2011	salient preservation	Landsat 4-5 TM	22_sp
12 October 2011	6 November 2011	tombolo formation	Landsat 4-5 TM	23_tf
6 November 2011	?	tombolo preservation	Landsat 4-5 TM	24_tp
?	19 May 2013	salient formation	Landsat 8 OLI/TIRS	25_sf
19 May 2013	30 August 2013	salient preservation	Landsat 8 OLI/TIRS	26_sp
30 August 2013	15 September 2013	tombolo formation	Landsat 8 OLI/TIRS	27_tf
15 September 2013	4 February 2019	tombolo preservation	Sentinel-2 L2A	28_tp
4 February 2019	14 February 2019	salient formation	Sentinel-2 L2A	29_sf
14 February 2019	2 June 2021	salient preservation	Sentinel-2 L2A	30_sp
2 June 2021	13 November 2021	tombolo formation	Sentinel-2 L2A	31_tf
13 November 2021	2 January 2022	tombolo preservation	Sentinel-2 L2A	32_tp
2 January 2022	8 February 2022	salient formation	Sentinel-2 L2A	33_sf
8 February 2022	7 April 2022	salient preservation	Sentinel-2 L2A	34_sp

depict the derived results.

3.3. Numerical Model Setup

3.3.1. Numerical Models, Time, and Sediment Parameters Description

The following TELEMAC-MASCARET modules (http://www.opentelemac.org/) (accessed on 15 February 2021) have been implemented in the present application: TELEMAC-2D (2DH current field prediction), TOMAWAC (wave field prediction), and SISYPHE (sand transport and seabed morphology prediction). The modules are run in parallel and internally coupled in order to simulate waves, currents, and sea bottom evolution, explained in Figure 5. Bathymetry is updated at each time step of coupling between SISYPHE and TELEMAC-2D and the new bathymetry derived feeds the three modules at the next time step.

TOMAWAC is a 3rd generation spectral wave model that simulates the generation and propagation of wind-induced irregular wave fields on triangular finite element meshes [50]. It is a phase-averaged directional spectral wave model that can reproduce the irregular wave shoaling, refraction and depth-limited breaking, the energy dissipation due to bottom friction, and the non-linear triad and quadruple wave–wave interactions, as well as diffraction effects to some degree.

TELEMAC-2D code solves the depth-averaged free surface flow equations as derived first by Barré de Saint-Venant in 1871 (see [51]), derived by integrating the Reynolds-averaged Navier–Stokes equations over the flow depth. The main results at each node of the computational mesh are the depth of water and the depth-averaged horizontal velocity components. The model takes into account the following phenomena: propagation of long waves, including non-linear effects, bottom friction, porosity phenomena, wave-induced currents (by coupling with TOMAWAC module), and coupling with the sediment transport module.

SISYPHE is an open-source sediment transport and bed evolution module of the TELEMAC-MASCARET SYSTEM that can be used to simulate sediment transport processes, i.e., bed load, suspended load, or total load. The prediction of morphological changes of bathymetry at each time step of simulation is derived by consideration of hydrodynamics (conservative laws of mass and momentum), sediment transport (predictors for sediment transport capacity), and bottom evolution (conservative law for sediment mass). The term “bed load” refers to sediment particles in water that are transported along the bed, while “suspended load” is the portion of the sediment that is carried by water flow that settles slowly enough such that it almost never touches the bed. Bed load moves by rolling, sliding, and/or saltating. Suspension is a result of turbulence in the flowing water and consists of particles generally of the fine sand.

The time step of the three models is set equal to 1 s to ensure simulation run stability, enabling the Courant number to be continuously less than unity during the simulation. Thus, the coupling period between TELEMAC2D and the other two models is set to 1. Furthermore, the small and identical time step between the models is expected to provide more realistic results compared with simulations with higher time steps and coupling periods. Nevertheless, it is noted that higher time steps (e.g., 10 s, 100 s, 300 s) and large coupling periods of the models tend to speed up the simulations but may cause stability issues depending on each simulation case and space discretization.

The sediment type is mainly sand, the median diameter dimension of which is considered equal to 0.23 mm, as derived from granulometric analysis applied to sediments of a neighboring sea coastal area. Sand density is considered equal to 2650 kg/m³ and water density equal to 1025 kg/m³.

The setup of the TOMAWAC model has been made as follows: the minimal frequency is considered equal to 0.06 s⁻¹ and the number of frequencies and directions of the directional wave spectrum are set to 24 and 36, respectively. For a growing wind sea, irregular wave modeling has been considered and a JONSWAP spectrum [52] for each sea state with a peak-enhancement factor equal to 3.3, while the directional spread is set equal to 27 deg. [53]. As for swells, in similarity to wind waves, a JONSWAP spectrum has also been used, but with a peak-enhancement factor larger than 3.3 [54] and a directional spread equal to 10 deg. (e.g., [55]). This is because when wind sea evolves into decaying swell, waves tend to become more regular. Furthermore, tide and wind forcing have not been considered in our simulations.

In the TELEMAC-2D model, the friction law on the bottom used in the hydrodynamic calculation is the law of Nikuradse, while in the SISYPHE model the morphological factor is equal to 1. In addition, the bed and suspended load transport formula by Soulsby [56] and van Rijn [57] is chosen among other transport formulas to consider the combined effects of currents and waves.

Additionally, the suspended load, i.e., the portion of sediment carried by the water flow that settles slowly such that it almost never touches the bottom, is maintained by the turbulence in the flow. Suspended sediment transport is considered in the SISYPHE model by solving the 2-D advection-diffusion equation, whereby turbulence is accounted for by the turbulent diffusivity of the sediment often related to the eddy viscosity and Schmidt number. In TELEMAC-2D, the Elder turbulence model is used in our simulations in order that the diffusion coefficient also includes the dispersion due to the heterogeneity of velocities on the vertical direction.

Moreover, as for the sediment deposition rate, which is a function of the settling velocity and the near-bed concentration, evaluated at the interface between the bed load and the suspended load, the SISYPHE model has been calibrated in order that the settling velocity be calculated by the model (as a function of the grain diameter and relative sediment–water density), while a Rouse profile is assumed for the vertical concentration distribution, which is theoretically valid in uniform steady flow conditions [58].

3.3.2. Bathymetry

Bathymetry and coastline information in the study area have been provided by the Hellenic Navy Hydrographic Service (https://www.hnhs.gr/ (accessed on 1 March 2021)). The model domain (Figure 6), including both the coastline and an area of open sea, was represented using a triangular unstructured finite element grid. A high resolution (5–20 m) mesh was applied within the area of the tombolo and also in the surf zone to allow high resolution of the significant bathymetric features, while a coarser grid resolution (up to 250 m) was applied offshore to reduce computational cost. Additionally, a filled contour map of the bathymetry in the intermediate and shallow water area is presented in Figure 7. These data of the study area will be used as the initial bathymetry for all simulated wave scenarios. As we can see in Figure 6 and Figure 7, the narrow area between Prasonisi and Rhodes Island is of the salient type.

4. Results

4.1. Investigation of Long-Term Morphological Changes and Seasonality of the Prasonisi Tombolo and Salient Based on Satellite Images

By using observations from Sentinel and Landsat satellite sensors, we are able to record and qualitatively evaluate the morphological changes of the Prasonisi tombolo (Table 2). It is noted that only qualitative information has been provided by satellite images regarding the state of morphology of the study area in certain time periods; thus, we can search more efficiently in these time periods for the derivation of the tombolo’s most constructive and destructive wave scenarios. Therefore, it is clarified that the bathymetry and the coastline information used in the numerical simulations have not been extracted from the available satellite imagery.

The satellite image analysis reveals the presence of eight cycles in the time period examined from 10 November 1995 to 08 February 2022. In this analysis, tombolo formation is considered the transition from the salient to the tombolo state, whereby the opposite applies for salient formation, as verified by the available satellite images. Tombolo preservation and salient preservation are the states that the landform of tombolo and salient preserved their states, respectively, since their formation.

It is noted that one cycle period is considered the duration of a tombolo formation till the next tombolo formation; thus, each cycle consists of four different morphological states of the study area; the tombolo and salient formation, and the tombolo and salient preservation. The mean cycle period of the eight cycles observed was approximately 43 months, while their minimum and maximum values were 13 months and 94 months, respectively. Moreover, the state of tombolo formation had a mean duration of 2 months, with a minimum value of less than 1 month (9 days) and a maximum value of 5 months. In a similar manner, the duration of salient formation had a mean value of 3 months and a minimum value of less than 1 month (7 days), but a maximum value of 13 months, which was greater than the maximum value of tombolo formation. Tombolo conservation lasted 22 months on average, while its minimum and maximum duration was 2 months and 66 months, respectively. Additionally, salient conservation lasted 13 months on average, while its minimum and maximum duration was 3 months and 41 months, respectively.

The monthly distribution of the four different morphology states of the Prasonisi tombolo shows that the tombolo formation appears mostly from June to November, as expected, since this period is characterized by more fair-weather conditions compared with the remaining months of the year (see Table 2 and Figure 8). The maximum frequency of tombolo formation is observed in July and August. The maximum frequency of salient formation is observed in January, February, and December, when most extreme storms occur. The only month that has zero frequency of salient formation is July. Furthermore, the formation of the tombolo has a lower duration than the salient formation. It is also noteworthy that the tombolo formation does not appear every year but only in 14 months in the time period of about 26 years examined, while the salient formation appears in 23 months in this period. Thus, the tombolo of Prasonisi is not characterized by significant seasonality. Nevertheless, as noted above, the salient formation is observed mainly in winter and the tombolo formation mainly in summer. The relative frequency of tombolo formation is 4.5%, while that of salient formation is 7.3%.

Regarding the other states of morphology, they are observed in all months, though the rate of salient preservation during the year was different from the corresponding trend of tombolo preservation. In particular, the relative frequency of tombolo preservation is 55.4% and the relative frequency of the salient one is 32.8%. Hence, although tombolo formation covered a lower time period than salient formation, tombolo preservation had a higher duration than salient preservation.

At this step, a preliminary statistical analysis of the offshore Hm0 is conducted referring to the two latest cycles, covering the time periods of 30 August 2013–2 June 2021 and of 2 June 2021–8 February 2022, as extracted from the observation of the satellite images (Figure 9). In this preliminary analysis, the mean values of the most crucial wave parameters, Hm0 and MWD, are estimated via Equations (4) and (7), respectively, by considering all wave data in time periods of tombolo and salient formation. In the latest cycle, during tombolo formation, Hm0 has a mean value of 0.99 m ± 0.43 (stdev), a maximum value of 2.54 m, and a prevailing value of 0.71 m. Additionally, the corresponding mean wave direction (MWD) has a prevailing value of 305 deg. (northwest). In the older cycle, during the same state of morphology, Hm0 has a mean value of 1.43 m ± 0.53 (stdev), a maximum value of 2.47 m, and a prevailing value of 0.55 m. Additionally, the corresponding MWD has a prevailing value of 309 deg. (northwest). As seen, as far as the tombolo formation is concerned, the two cycles have similar wave height conditions, and the construction of the tombolo is succeeded from waves coming from the same direction, i.e., from the northwest.

As far as the salient formation is concerned, in the latest cycle, Hm0 has a mean value of 1.63 m ± 0.73 (stdev), a maximum value of 2.68 m, and a prevailing value of 2.36 m. Additionally, the corresponding MWD has a prevailing value of 280 deg. (west). In the older cycle, during the same state of morphology, Hm0 has a mean value of 1.64 m ± 0.77 (stdev), a maximum value of 4.53 m, and a prevailing value of 1.73 m. Additionally, the corresponding MWD has a prevailing value of 140 deg. (southeast). As seen, the two cycles have similar mean Hm0 values, albeit the maximum and mode values are significantly different in the two cycles. Furthermore, the construction of the salient is accomplished due to waves coming from two different directions in the two cycles. Moreover, it is observed that waves coming from the northwest could construct either a tombolo or a salient, depending on their intensity and wave energy.

4.2. Long-Term Wave Climate

The prevailing wave scenarios that correspond to time periods of tombolo formation are characterized by low and (more rarely) moderate wave heights (Figure 10). In periods of tombolo formation, Hm0 has a mean value of 1.05 m, while its maximum value is observed in period 3_tf, being equal to 2.83 m. The period of tombolo formation with the shortest duration (9 days) is the one symbolized by 7_tf and consists of fair-weather waves associated with significant wave heights with a maximum value of 1.81 m and a mean value of 0.56 m. Additionally, those sea states appear to have low variability with respect to mean wave direction as they all come from the NW direction.

In the case of the salient formation, high significant wave heights and extreme wave events are observed, showing higher variability regarding the mean wave direction compared with those of the tombolo formation. Extreme wave events usually come from NW, SE, and SW (Figure 11). In periods of salient formation, Hm0 has a mean value of 1.68 m and a maximum value of 5.67 m. The latter appeared in the 29_sf period, whereas the period of salient formation with the lowest maximum value (4.10 m) of Hm0 is the 9_sf period.

4.3. Extracting Wave Scenarios for the Investigation of the Tombolo and Salient Formation

The implementation of the swell criterion to subsets of data covering the time periods of tombolo formation, as they are defined from satellite image analysis, shows that wind seas and young swells are the more frequent types of sea states (Figure 12). The application of the methodology applied for the extraction of the characteristic wave scenarios and the segments of their successive data reduced the sample size by 85% of the initially available data.

The wave scenarios extracted have a mean Hm0 value equal to 0.85 m, with a standard deviation of 0.28 m, and only the fifth scenario has a representative Hm0 that is relatively high attaining the value of 1.5 m (Figure 13). In all wave scenarios, Tp varies between 4.5 s and 7.8 s, while MWD appears with very low variability since the majority of waves comes from the northwest. Additionally, the maximum value of Hm0 observed in the fifth scenario in the sequence (Figure 13) is equal to 1.85 m.

The application of the aforementioned methodology for the time periods of salient formation reduced the sample size by 90% of the initial dataset and shows that wind seas and young swell are the more frequent types of sea states, similarly to the case of tombolo formation (Figure 14).

The wave scenarios extracted in the period of the salient formation (5_sf) have a mean representative Hm0 value equal to 0.98 m with a standard deviation of 0.57 m, presenting thus higher variability than the wave scenarios in the period of the tombolo formation (7_tf) (Figure 12 and Figure 13). Additionally, the maximum value of the representative Hm0 is 3.68 m and the maximum value of Hm0 is 5.83 m. In all scenarios, Tp varies between 1.6 s and 10.3 s, and the MWD varies between 2 deg and 356 deg from the north, having a prevailing value of 303 deg (southwest) (Figure 15). In comparison with the wind and wave characteristics detected during the tombolo formation, the Hm0, Tp, and MWD are more variable during the salient formation.

Even though a significant data reduction (up to 90% (This is the value of data reduction extracted from the total wind and wave dataset.) can be accomplished via the use of the methodology developed, the most characteristic wave scenarios should be selected among the scenarios extracted. These scenarios are selected regarding the degree of how representative they are, depending on the tombolo and salient formation. The first four scenarios have the same MWD, as they all come from the northwest. The first two scenarios refer to wind-wave events of fair and moderate weather characterized by different values of Hm0 and Tp, aiming at deriving the most effective one regarding tombolo formation with respect to wave intensity. The other two scenarios concern swell seas in order to be compared with the wind-wave scenario 1 of the same offshore Hm0 but different Tp. In this manner, the impact of the wave period on the short-term morphological changes of the tombolo’s formation will be investigated. Moreover, the last two cases are scenarios of more extreme weather than the other scenarios tested, coming from the southeast and the northwest. These MWDs are the most probable ones associated with storm events appearing in periods of salient formation.

It is noted that the characterization of sea states depending on their intensity (i.e., fair, moderate, and extreme wave conditions) has been based on the estimated impact of a large number of different simulations on the sea bottom bathymetry and not on a statistical analysis, albeit the two approaches might be in compliance to a considerable degree with each other. From this numerical modeling analysis of sea bottom evolution in the study area, it has been derived that sea states with Hm0 values lower than 0.5 m can be neglected [49], while Hm0 values ranging between 0.5 m and 1.0 m are considered as fair wave conditions. Additionally, Hm0 values in the range of 1.0 m to 1.8 m are moderate wave conditions, and more intense waves of Hm0 greater than 1.8 m belong to extreme wave conditions in the study area. From a statistical point of view, the 50th percentile of Hm0 is equal to 1.05 m, and the 90th and 95th percentiles of Hm0 correspond to 1.89 m and 2.20 m, respectively.

The characteristic wave scenarios that will be simulated and presented in

Table 3. Wave scenarios to be simulated.

a/a	Hm0 (m)	Tp (s)	MWD of Propagation (deg. from North)	Duration (h)
1	1.00	4.60	120	3
2	1.50	5.55	120	3
3	1.00	6.40	120	3
4	1.00	7.80	120	3
5	3.00	7.70	130	3
6	4.70	9.15	310	3

have been extracted via the aforementioned method from subsets 7_tf and 5_sf, for tombolo and salient formation, respectively. In particular, in the subset of data 7_tf, the prevailing wave event has an Hm0 value close to 1.0 m, and the most intense wave event has an Hm0 value equal to 1.5 m, while all wave events come from the northwest. As for Tp in the dataset of 7_tf, it varies between 4.5 s and 7.8 s. Hence, the first four wave scenarios have been chosen to be simulated based on the wave event variability in this short time period of tombolo formation. Additionally, regarding the salient formation time period 5_sf, the most extreme wave events coming from two different directions have been chosen in order that the most destructive wave conditions of the tombolo be examined. The real time simulation of each wave scenario is set equal to 3 hr, being regarded adequate enough to cause a considerable short-term morphological bottom change. This duration ensures adequately reliable results, but also speeds up the numerical investigation. It is noted that in the present paper, the same duration is applied for all wave scenarios to be simulated in order that they can be compared on equal terms, regarding their effects on short-term sea bottom evolution.

4.4. Numerical Investigation of the Tombolo and Salient Formation

Since the initial bathymetry of the interest area tested is salient, wave scenarios are simulated and tested with respect to whether they result in further erosion of the salient area or whether they cause sand accretion in the salient area. The first four wave scenarios are considered tombolo constructive scenarios, while the last two scenarios (fifth and sixth) are considered salient constructive scenarios. Results are provided regarding the spatial distribution of Hm0 (m), current velocity (m/s), and bottom level change (m) in the sea area of interest for all scenarios simulated (Figure 16, Figure 17 and Figure 18). As observed in the numerical investigation, the first four wave scenarios of fair and moderate wave conditions seemed to have a tombolo constructive action, while the last two ones have an eroding effect upon the salient area. This is obvious, since the fifth and sixth scenarios correspond to extreme wave conditions. Nevertheless, a more insightful description of the impact of each critical parameter inducing the sea bottom change is conducted in the Section 5, illustrating the causal relation of each parameter with the corresponding result.

5. Discussion

From the six representative scenarios examined, four of them (one, two, three, and four) resulting in sand deposition in the area of tombolo are hereinafter characterized as tombolo constructive wave scenarios. In a similar manner, the fifth and the sixth wave scenarios resulting in sand erosion in the area of tombolo are characterized as tombolo destructive, or equivalently salient constructive, wave scenarios. Therefore, the active domain of the study area for tombolo formation (see Figure 16 and Figure 17) and for salient formation (see Figure 18) are defined as follows: the first one is considered the union area of sand accretion of the first four tombolo constructive wave scenarios (

Table 4. Comparison of tombolo constructive wave scenario results of total volume of sand accretion and loss in the active area for tombolo formation.

Wave Scenario a/a	Hm0 (m)	Tp (s)	MWD (deg)	Total Sand Accretion (m3)	Total Sand Loss (m3)
1	1.00	4.60	120	1207	0
2	1.50	5.55	120	25,154	0
3	1.00	6.40	120	4567	0
4	1.00	7.80	120	7191	0

), and the second one is the union area of sand erosion of the fifth and the sixth tombolo destructive scenarios (

Table 5. Comparison of tombolo destructive wave scenario results of total volume of sand erosion in the active area for salient formation.

Wave Scenario a/a	Hm0 (m)	Tp (s)	MWD (deg)	Total Sand Accretion (m3)	Total Sand Loss (m3) (2)
5	3.00	7.70	130	3.00	281,013
6	4.70	9.15	310	4.70	194,724

).

The total volume of sand accretion and sand loss in these two active domains are characteristic measures of the degree that each combination of simulated offshore wave conditions (Hm0, Tp, MWD) contribute to the corresponding tombolo (Table 4) or salient (Table 5) formation. According to this analysis, it is verified that wind-wave events of low intensity and swell events usually result in sand deposition in the tombolo’s active formation area. Specifically, wave events of low Hm0 values (e.g., ~1.0 m in this case study) have a constructive action for tombolo formation, especially when combined with high Tp values. As for moderate waves (1.0 m < Hm0 < 1.8 m), these events are more effective regarding tombolo construction than events of low intensity, as they correspond to higher wave energy conditions capable of initiating sediment transport in the active area.

Offshore wave conditions, i.e., Hm0, Tp, and MWD, are the critical wave parameters that can have a different impact on the tombolo or salient sea bottom evolution. In the general case of tombolo formation, tombolos are formed mainly by wave refraction and diffraction around the former island [59], whereas in the other cases, their formation is usually as a result of a unidirectional longshore drift due to oblique wave incidence on one side of the tombolo. For instance, in the case of Prasonisi, low and moderate values of Hm0 (<1.8 m) coming from a certain directional sector (NW) can result in sand deposition in the salient or tombolo sand bar, while extreme values of Hm0 (>1.8 m) coming from the NW or SE have as a consequence the erosion of the sandbar. As for swell events, when they are associated with low values of Hm0, they could result in sand accretion, especially when they are combined with high values of Tp (e.g., [60,61]).

It is noted that the presence and location of wave breaking plays an important role in coastal morphodynamics. For many decades, the energy principle approach (e.g., [62,63]) is commonly used to investigate the wave-induced longshore sand transport [64]. The quantity of sand transported due to the combined wave and current action is a function of the energy available for transporting the sediments. Thus, wave breaking is a crucial factor for energy dissipation producing the littoral drift. In particular, outside the surf zone, energy dissipation is primarily due to bottom friction, hence the main mode of sediment transport is the bed load, while inside the surf zone is significantly increased mainly due to turbulence, and consequently the primary mode of transport is the suspended load [64]. Moreover, following Shibayama et al., [65], the agitation of the bed materials in the surf zone is strong and the suspended sand concentration in the vicinity of the wave breaking point is extremely high. These suspended sediments, plus some of the sediments on the seabed, are then carried along the shoreline by the longshore current, which has its maximum velocity near the breaker line [66].

Therefore, as in coastal morphodynamics, the location of the surf zone is of major importance for the sea bottom evolution near the area of the tombolo/salient. As a result, following the wave breaking criterion by Thornton and Guza [67], within the saturated (The saturated breaking zone is the zone where essentially all waves are breaking. For gently sloping shorefaces, (typically slope < 0.02), it may be assumed that wave breaking is saturated throughout the surf zone [68].) depth-induced breaking zone for irregular waves, the spectral significant wave height $H_{mo,b}$ may be related to water depth, approximately, as follows:

$$H_{m0,b}=0.6 d$$

(10)

where $d$ is the water depth.

It can be observed that when Tp is increased, and thus sea bottom effects on wave characteristics are met in deeper waters compared with lower Tp values, the width of surf zone is wider, and its location is closer to the area of interest, resulting in further sand deposition in the tombolo area (Figure 19a,c,d). Therefore, wave breaking of fair weather that takes places near the tombolo’s active area contributes to the tombolo’s formation. The condition of fair weather is necessary in order that waves and currents have a combined constructive action, while the appropriate characteristics of the surf zone make this constructive action more enhanced.

Contrariwise, extreme offshore wave conditions coming from the NW and SE (Figure 19e,f) have a destructive action upon the sandbar of the salient/tombolo. This is because the surf zone is located in deeper water depths compared with fair and moderate wave conditions, and the surf zone has a large width. Moreover, the current velocity in the submerged sand bar area of the salient is relatively high, resulting in the latter’s further erosion.

Additionally, as in the present case when a submerged sandbar is located in shallow waters between the mainland and the tied island, waves traveling from offshore towards the sandbar undergo increasing bottom effects (e.g., shoaling, refraction, breaking, diffraction, etc.). However, when they pass the sandbar, the water depth is increased again, which could lead to decreasing bottom effects downstream of the submerged sandbar. The bottom complexity in this zone under certain wave conditions, and mainly when wave breaking takes place and as a consequence the flow becomes extremely turbulent, could cause the presence of the surf zone upstream and downstream of the sandbar. This phenomenon enhances the sea bottom’s morphological changes in this region. This is observed in wave scenarios two, four, and six (see Figure 19b,d,f), respectively. Moreover, it is observed that a breaker index greater than 0.6 results in sand accretion (e.g., Figure 19b) or erosion (Figure 19d,f) to a greater degree than in cases of values lower than 0.6.

Therefore, low non-breaking waves or waves of fair weather whose breaking takes place near the active domain of the study area for tombolo formation contribute to the tombolo’s formation, whereas wave breaking of extreme weather that takes place near the active domain of the study area for salient formation has a destructive action upon the sandbar.

Moreover, a large width of the surf zone is responsible for a sea bottom morphological change to a greater degree than in the case of a smaller width. Particularly, an adequate surf zone width and low offshore conditions are necessary for a considerable sand deposition in a salient/tombolo sandbar (scenarios two and four), and a large width of the surf zone and extreme offshore wave conditions have a more destructive power (as in scenarios five and six).

Additionally, the maximum value of the breaking index in a coastal area is site-specific, since it depends on the bottom slope, the bottom complexity, and the wave conditions (e.g., [69]). However, it is noticed that higher values of the wave breaking index have a more drastic action on the sea bottom. Furthermore, bathymetry is also a significant factor for coastal morphodynamics, since a different coastal bathymetry can alter the results to a certain degree, e.g., the location and width of the surf zone. For example, the different bathymetry in the two sides of the tombolo’s intermediate and shallow water area is the main reason that an extreme wave event coming from the NW, even of lower intensity, can be more destructive (scenario five) than another one of higher intensity coming from the SE (scenario six).

Considering the processes discussed above, the Prasonisi tombolo is not formed mainly by wave refraction and diffraction around the Prasonisi Islet (Figure 20A, but its formation is usually as a result of a unidirectional longshore drift due to oblique wave incidence on one side of the tombolo (see Figure 20B (e.g., [59,70])). The first case refers to the most usual case of tombolo formation induced by the longshore drift from the two sides of the tombolo, resulting in a symmetry of sea bottom bathymetry around the tombolo, while the second one describes tombolo construction as a result of a unidirectional longshore drift from one side of the tombolo. On the contrary, the destructive action of extreme wave conditions could result in salient formation (Figure 20C depicts a salient construction from the eroding effects of extreme waves).

6. Conclusions

The methodology developed for the extraction of wave events has a main property concerning the investigation of short-term sea bottom evolution. In particular, wind-wave events consist of wind and wave successive data of similar wind and wave characteristics; thus, from a statistical point of view they can be represented more properly by their mean values compared with wave events extracted from non-successive data or data presenting a significant variability. In this manner, wave events of a variety of wave intensities can be derived, and thus the analysis is not only restricted in extreme wave events, as is usually the case in extreme value analysis. Therefore, either coastal erosion or accretion as well as beach recovery can well be considered and investigated. Additionally, the chronology of wave events could also be used when studying the coastal morphodynamics and sea bottom evolution in a sequence of wave events. However, this has not been studied in the present paper, as a preliminary investigation of a tombolo’s most representative wave conditions is conducted.

As for the main factors that play a significant role in tombolo formation, these are: the offshore wave conditions, the location and width of the surf zone, the bottom slope, the maximum value of wave breaking index in the study area (a function of the previously mentioned parameters), and the initial bottom bathymetry before the study area is exposed to a new wave event. Additionally, a low ratio of wave breaking height towards breaking depth that takes place near the salient sand bar area contributes to the tombolo formation, while its high ratio in combination with a large surf zone width contributes to the salient formation. Moreover, as expected, it is shown from the present paper that during high-energy conditions, the sand bar of the tombolo/salient, in similarity to the beaches, is attacked severely by the incoming waves, usually resulting in erosion processes, while fair- and moderate-weather waves and swells return the sediments towards the tombolo/salient sandbar.

Referring now to the case of Prasonisi, the area of the Prasonisi tombolo is protected to a great degree by the Prasonisi island, hence opposite flows from diffraction around the obstacle are not the main reason causing tombolo formation, but fair- and moderate-weather waves and swell waves coming from the northwest. Moreover, the study area can be eroded from extreme waves coming from the same direction as waves that cause sand deposition in the study area (northwest) and from the southeast. Therefore, the tombolo seems to have similar morphodynamics to beaches exposed to waves coming from two directional sectors, although the area is more protected from waves coming from the northwest than the southeast.