Assessment of the Coastal Vulnerability to the Ongoing Sea Level Rise for the Exquisite Rhodes Island (SE Aegean Sea, Greece)

Abstract

The foreseeable acceleration of global sea level rise could potentially pose a major threat to the natural charm and functional integrity of the world-renowned tourist coastal attractions of Rhodes Island, as a result of the anticipated increasing frequency of flooding and erosion events. Hence, this study aims to determine the most vulnerable segments (in terms of physical impact) of the Rhodes coastline through the widely accepted coastal vulnerability index (CVI), applying a combination of well-known, broadly used approaches and methods. The frequency distribution of the current CVI along the island’s coastline suggests a rather worrying high to very high vulnerability of 40%. In addition, a CVI projection to the end of the 21st century (based on the Intergovernmental Panel on Climate Change predictive scenarios) indicates an enhancement of the total vulnerability by 48%, mainly focused on the majority of the western coastline. Hence, a considerable number of popular coastal destinations in the island shall remain under unignorable threat and, therefore, coastal managers and decision-makers need to hatch an integrated plan to minimize economic and natural losses, private property damage and tourism infrastructure deterioration from flooding and erosion episodes, which will most likely be intensified in the future.

1. Introduction

The coastal areas are low-lying dynamic landscapes that are affected by a variety of natural processes, such as sea level change, wave and longshore and/or cross-shore current activity, tidal fluctuations, vertical land movements and sediment transport dynamics. Hence, the drastic change or intensification of any of the previous processes due to anthropogenic or natural causes may initiate major coastal degradation, including land loss, severe damages of infrastructure, seawater pollution and decline in marine biodiversity [1,2,3,4,5,6,7,8,9,10].

Annually, the European coastal regions, and especially the Mediterranean islands, attract hundreds of millions of tourists, who contribute substantially to local economic growth [2,4,11]. However, in several cases, mass tourism initiates environmental deterioration when the carrying capacity of a given area is overweighed [12,13]. In addition, the supersaturation of other human activities in coastal regions, such as extensive urbanization, industrial, agricultural and aquaculture/fishery farming, and maritime transport, may also adversely affect the environmental state [10].

The sea level rise (SLR) is one of the most critical pressures exerted on the coasts worldwide (e.g., [11,14,15]). This process is expected to evolve with higher rates in the near future due to: (i) the continuous increase of the average global surface air temperature (by the end of this century, a 3–5 °C temperature rise above pre-industrial levels has been predicted if the current trend in global greenhouse gas emissions continues (https://www.reuters.com/article/us-climate-change-un-idUSKCN1NY186; accessed on 26 July 2021)), (ii) the increase of the upper ocean heat content, that will inevitably cause acceleration of the seawater thermal expansion [16], and (iii) the radical change in the Greenland ice-sheet dynamics along with the complete destabilization of the west Antarctic ice-sheet and partial instability of the east Antarctic glaciers [17]. According to the Fourth and Fifth Assessment Reports of the Intergovernmental Panel on Climate Change (IPCC) [18,19], the mean global SLR by 2100 is projected to reach up to 60–90 cm above the level relative to 2000 (according to the RCP2.6, RCP4.5, RCP6 and RCP8.5 scenarios with a 95% confidence interval). This process constitutes a severe hazard for the coastal sustainability, threatening the quality and value of coastal ecosystem services and, consequently, preventing the smooth development of the tourism industry. In addition, since the level of coastal vulnerability to climate change depends on the interaction of the type, magnitude and rate of climatic variation with the coastal system’s sensitivity and adaptive capacity [18,20], the rate of shoreline retreat is anticipated to be enhanced and, therefore, extensive coastal areas (including sandy, dune and sandy cliff systems) of high tourism value, especially in the Mediterranean Sea, are going to be flooded and will shrink.

Rhodes Island is one of the most stunning islands in the Aegean Sea, representing a top-rated attraction for the global tourism and, consequently, being a source of significant income and employment for the local society as well as a key contributor to the country’s welfare. Hence, the present work assesses the coastal vulnerability of Rhodes to the current and future SLR by 2100, using the well-known multi-hazard coastal vulnerability index (CVI) approach. For the performance of the CVI methodology on local, regional and supra-regional scales, essential physical factors affecting the coastline configuration are considered, such as coastal geomorphology, regional coastal slope, rate of sea level change, significant wave height, historical shoreline displacement rate, tidal range, annual tropical storm probability, annual hurricane probability, hurricane frequency and intensity, annual intensity of extratropical cyclones and storm/hurricane surge (see [21] and references therein). The CVI approach has been applied and/or adapted by a plethora of research efforts for the vulnerability assessment of coastlines around the globe, while its robustness is proven by the fact that more than 30% of the studies targeted on coastal vulnerability have used the relevant methodology [21]. The advantages of the CVI application are twofold: (i) the physical change most likely to occur due to SLR can be sufficiently evaluated, despite the fact that the subjective nature of the expert opinion remains a problem, and (ii) a preliminary planning for managing and protecting coastal resources and infrastructure in the areas of interest can be facilitated.

The CVI determination for Rhodes, for the very first time, can be very helpful because the island is one of the hottest travel destinations in the Mediterranean Sea due to its highly diverse landscape, excellent climate conditions, fast developing state-of-the-art infrastructure and great cultural heritage. It should be emphasized that the island’s population almost doubles during summer, reaching a maximum density of 160 inhabitants per km², while the corresponding density during winter is 83 inhabitants per km² [22]. Hence, the potential significant negative impact of climate change on the physical vulnerability of the island’s coasts is of vital concern for the coastal zone managers and decision-makers in order to implement appropriate and effective policies for the current and future investments. In particular, for the development of optimal and rational investment plans concerning the coastal commercial and tourism areas, the assessment of the relative physical vulnerability along the coastline and the identification of hotspots is certainly crucial.

Even though the magnitude of the coastal vulnerability of Rhodes due to SLR cannot be easily and accurately assessed due to additional interacting factors, such as sediment budget availability, near-shore hydrodynamics and local morphological irregularities, which are not always quantifiable (mainly for future climate conditions) [23], the CVI approach represents a widely used and highly indicative tool for a holistic determination of the physical coastal vulnerability for nowadays and the coming decades as well. As beach inundation and subsequent erosion are going to inevitably modify the tourism carrying capacity [1,24,25], it is absolutely essential for the policymakers, spatial planners, managers and stakeholders to be aware about the current and future vulnerability of a given coast in order to find out the appropriate solutions for the mitigation of any severe effects caused by the global SLR.

Finally, besides the assessment of the current and future CVI distribution along the Rhodes coastline, the present study suggests spatially focused and sustainable management solutions in order for the adverse effects of SLR to be mitigated.

2. Area of Concern

Rhodes is the largest island of the Dodecanese insular complex (located in the southeast Aegean Sea), having an area of 1401 km², a coastline length of 253 km and a lenticular shape that extends along an axis of SW–NE direction (Figure 1a). It lies on the eastern end of the Hellenic Arc, where the rigid oceanic Nubian plate is subducted under the continental Aegean microplate. The occurrence of the Pliny and Strabo trenches in the area, resulting from strike–slip faulting [26,27,28,29], and the development of the 4 km deep Rhodes basin to the east of the homonymous island, suggest intense Plio-Quaternary tectonic deformation. The ‘recent’ tectonic activity associated with Holocene seismic events is indicated by the uplifted tidal notches observed along the Rhodes shoreline [30,31,32]. These notches can be found at various altitudes, ranging from a maximum of +3.8 m near the homonymous city (Rodos City; Figure 1) to lower heights towards the southern parts of the island [30].

Regarding the seasonal wave climate characteristics in the broader region of Rhodes, the mean significant wave height values in winter range from 0.6 to 1.3 m, with the maxima occurring in the northeastern Aegean Sea, central Aegean Sea (north of the Cyclades insular complex) and Crete–Kythira and Crete–Karpathos straits [33,34]. During Summer, the dominant waves derive from the north due to the occurrence of the Etesian winds, which diminish in winter when violent cyclonic storms develop [35]. Focusing strictly on Rhodes, the prevailing winds blow from the west, west–northwest and west–southwest. However, the strongest events come from southeast and east–southeast, with the wind speed exceeding 15 m s⁻¹ [36]. Likewise, the dominant waves come from the west, west–northwest and west–southwest, with a frequency of occurrence of more than 25%. In the northern and eastern parts of the island, the highest waves originate from the south–southeast and south–southwest, while in the southern and western parts, the strongest events propagate from the west–northwest and west [36]. Finally, in the southern Rhodes, when strong winds blow over long fetches during extreme storms, the wave heights can exceed 6.5 m, while in the northern Rhodes, the wave heights can reach up to 4.3 m [36].

The island’s littoral zone consists of a succession of rocky coasts and long and wide beaches [37,38]. Rocky cliffs are mainly situated in the central-eastern and central-western parts of the island, near Lindos and Monolithos–Kritinia areas respectively (see Figure 1), formed due to both local vertical tectonic movements and relatively low weathering rates of the alpine and post-alpine substrates [39,40,41,42] (Figure 1b).

The hydrographic network of Rhodes (Figure 1a) is not well-developed due to the regional tectonic activity. It is mostly directed from NW to SE, while V-shaped valleys exhibiting a W–E orientation have developed throughout the island due to the combined action of erosive and tectonic processes. Various beach morphologies occur at the estuaries of the hydrographic network, mainly including sand, pebbles and limited mud fractions, and governed by torrential and marine depositional processes, while sand dunes and beach rocks are distinct features of several coastal areas [27,39,43,44]. The smooth slope and remarkable length (151 km) of the sandy beaches, which cover 60% of the total island’s coastal zone, together with the outstanding environmental conditions, have made Rhodes one of the most visited tourist hotspots in the Aegean Sea [45].

The highest concentrations of tourists occur in the northwestern, north and northeastern coasts of Rhodes, which accommodate the majority of tourism infrastructure and facilities [3,6,9]. For example, busy ports or superior hotel complexes have been built in Rodos City, the neighboring Faliraki and Afantou bays to the eastern coastal zone and the Ialysos area to the west (Figure 1). The island attracts more than 2,000,000 non-resident visitors per year, that contribute to ~10% of the Greek tourism product, with the gross value added by the local tourism sector being estimated at more than 670 million euros per year [9], excluding, however, the current pandemic period. The rapid urban growth of Rhodes, intensifying the human impact on the natural processes, has put an environmental stress on the coastal zone, contributing considerably to the enhancement of erosion processes with inevitable loss of land [9] and, thus, initiating some degradation of the tourism product.

Finally, a great part of the coastal zone in the southern, southwestern and central-western Rhodes, together with the coast of Lindos and the Masari–Archangelos area (Figure 1), are protected by the European network ‘Natura 2000’ due to their high ecological significance. In addition, the remains of the medieval town in the modern Rodos City (the capital of Rhodes, where 48% of the total population is concentrated [22]) is an UNESCO World Heritage Site, enhancing the island’s tourism value.

3. Methodology

For the CVI determination (current and projected to 2100), the required database (being both quantitative and qualitative) was first established through the compilation and numerical analysis of field data and topographic and geo-environmental datasets. Then, the database was processed through specific commands/applications within the ArcGIS platform using (i) the ArcMap software (v. 10.5) together with Python scripting and the ModelBuilder tool [46,47], and (ii) the digital shoreline analysis system (DSAS) software (v. 5.0). Searching in the literature for a coastal vulnerability ranking scheme to better categorize the physiographic and hydrodynamic conditions in Rhodes, a combination of the widely used approaches of Gornitz [23], Gornitz et al. [48] and Thieler and Hammar-Klose [49] was applied to evaluate the relative ranking of a number of crucial CVI variables, such as: (a) mean absolute SLR rate (mm y⁻¹), (b) mean significant wave height (Hs) (m), (c) mean tidal range (m), (d) lithology and geomorphology of coastal landforms (indicating the relative resistance of a given landform to erosion), (e) coastal slope (%), (f) mean shoreline displacement rate (±m y⁻¹, due to erosion or sediment accumulation) and (g) mean vertical coastal land movement rate (mm y⁻¹).

Table 1 presents the CVI variables and the relative ranking scheme with the classic 1–5 linear scale at an order of increasing vulnerability. All ranking scores were combined to determine the composite CVIs at 20,082 cells of 50 × 50 m along the entire coastline of Rhodes, with the relevant index calculated by the broadly used formulation based on the square root of the product mean of the ranking variables divided by their number:

CVI = [(a × b × c × d × e × f × g)/7]^0.5

(1)

Then, the calculated CVI values for all coastal segments were used to derive the relevant frequency distribution. Based on the five quintiles of the CVI distribution, an equal number of vulnerability classes analogous to those applied for the relative ranking of the CVI variables were prespecified and the Jenks (natural breaks) optimization algorithm (provided by the ArcMap software) was applied to determine the best arrangement of the values into the prespecified number of vulnerability classes. The algorithm minimized the variance within the classes and maximized the variance between them, thus producing an optimal classification. Finally, the results were mapped via the ArcGIS platform and allowed the identification of the most physically vulnerable coastal segments.

Cogswell et al. [50] argue that since all CVI variants are sensitive to the number of variables used, i.e., the greater the number of variables used, the higher the increase of the skewness and mean of the CVI distribution, the adoption of the geometric mean could be more appropriate. However, it should be emphasized that the CVI formulation used in this study (with the inclusion of seven physical variables) offers the advantage of expanding the range of the final value, while, in parallel, dampening the extreme ranges. In addition, it is sufficiently sensitive to slight changes of the variables, since the transition of only one variable to the next lower or upper vulnerability level causes a CVI variation ranging from 10% to 42%. Finally, the applied CVI formulation has the benefit of integrating the selected physical variables in a quantifiable way, providing numerical estimates that albeit cannot be directly related to particular physical effects, though they may highlight those coastal segments where the impact of SLR will be the greatest [51].

Table 1. Relative vulnerability classes for the CVI variables, implying the importance of their individual contribution to the expected physical coastal changes due to sea level rise (SLR). The physical variables shown together with their thresholds follow the Gornitz [23], Gornitz et al. [48] and Thieler and Hammar-Klose [49] approaches. The suitability of the applied vulnerability ranking scheme for the CVI determination of Rhodes Island is supported by the observations and measurements made during the performance of a plethora of assessment field surveys in the island’s coastal zone as well as by the concepts presented in the studies of Pendleton et al. [52,53], Kaliraj and Chandrasekar [54], Karymbalis et al. [55], Dimou et al. [56], Tragaki et al. [57] and Kantamaneni et al. [58]. Key: [] represents a closed interval; ( ] and [ ) represent half-open intervals.

CVI Variable	Relative Ranking Scheme for Each CVI Variable
	Very Low	Low	Moderate	High	Very High
	1	2	3	4	5
(a) Mean absolute SLR (mm y−1)	<1.8	[1.8–2.5)	[2.5–3)	[3–3.4]	>3.4
(b) Mean significant wave height (m)	<0.55	[0.55–0.85)	[0.85–1.05)	[1.05–1.25]	>1.25
(c) Mean tidal range 1 (m)	>6	(4–6]	(2–4]	[1,2]	<1
(d) Lithology and geomorphology (including man-made structures)	Rocky cliff coasts, alpine and post-alpine substrate, beach rocks, concrete retaining walls, port facilities, rip raps	Medium cliffs, indented coasts	Low cliffs, alluvial plains, fluvial–alluvial beaches	Cobble beaches, lagoons	Sand beaches, deltas, sandy beaches under erosion
(e) Coastal slope (%)	>20	(7–20]	(4–7]	[2.5–4]	<2.5
(f) Mean shoreline change (±m y−1)	>+2	(+1 to +2]	(–1 to +1]	[–2 to –1]	<–2
(g) Mean coastal vertical velocity (mm y−1)	>5	(4–5]	(2.5–4]	[1.1–2.5]	<1.1

¹ The peculiar classification of the tidal range is explained in Section 4.1.3.

Table 1. Relative vulnerability classes for the CVI variables, implying the importance of their individual contribution to the expected physical coastal changes due to sea level rise (SLR). The physical variables shown together with their thresholds follow the Gornitz [23], Gornitz et al. [48] and Thieler and Hammar-Klose [49] approaches. The suitability of the applied vulnerability ranking scheme for the CVI determination of Rhodes Island is supported by the observations and measurements made during the performance of a plethora of assessment field surveys in the island’s coastal zone as well as by the concepts presented in the studies of Pendleton et al. [52,53], Kaliraj and Chandrasekar [54], Karymbalis et al. [55], Dimou et al. [56], Tragaki et al. [57] and Kantamaneni et al. [58]. Key: [] represents a closed interval; ( ] and [ ) represent half-open intervals.

CVI Variable	Relative Ranking Scheme for Each CVI Variable
	Very Low	Low	Moderate	High	Very High
	1	2	3	4	5
(a) Mean absolute SLR (mm y−1)	<1.8	[1.8–2.5)	[2.5–3)	[3–3.4]	>3.4
(b) Mean significant wave height (m)	<0.55	[0.55–0.85)	[0.85–1.05)	[1.05–1.25]	>1.25
(c) Mean tidal range 1 (m)	>6	(4–6]	(2–4]	[1,2]	<1
(d) Lithology and geomorphology (including man-made structures)	Rocky cliff coasts, alpine and post-alpine substrate, beach rocks, concrete retaining walls, port facilities, rip raps	Medium cliffs, indented coasts	Low cliffs, alluvial plains, fluvial–alluvial beaches	Cobble beaches, lagoons	Sand beaches, deltas, sandy beaches under erosion
(e) Coastal slope (%)	>20	(7–20]	(4–7]	[2.5–4]	<2.5
(f) Mean shoreline change (±m y−1)	>+2	(+1 to +2]	(–1 to +1]	[–2 to –1]	<–2
(g) Mean coastal vertical velocity (mm y−1)	>5	(4–5]	(2.5–4]	[1.1–2.5]	<1.1

¹ The peculiar classification of the tidal range is explained in Section 4.1.3.

The data for the evaluation of the relative ranking of the CVI variables (present and projected to 2100) were collected as described below:

(a) The rate of the absolute SLR during the entire Holocene was obtained from the predictions of the Lambeck’s glacio-hydro-isostatic model applied to the Aegean and Mediterranean Seas [59,60], while recent SLR variations in the Aegean Sea were obtained from the studies of Cazenave et al. [61], Tsimplis and Rixen [62], Calafat and Gomis [63], Tsimplis et al. [64], Landerer and Volkov [65], Woppelmann and Marcos [66] and Bonaduce et al. [67], all based on the analysis of instrumental and satellite data. In addition, an investigation of the SLR changes around Rhodes was carried out by analyzing datasets obtained from (i) the tide gauge station of Rhodes (ID: 1243 and 2167), downloaded from the permanent service for mean sea level (PSMSL) [68] and corrected for the vertical crustal movements, with a completeness of 72.5% for the 1969–2016 time interval, and (ii) monthly sea level maps (1/8° × 1/8° grid) of absolute dynamic topography, downloaded from the Copernicus marine environment monitoring service (CMEMS) [69] and covering the 1989–2016 time interval. Further, the probabilistic SLR projections at Rhodes coasts throughout the 21st century for the RCP2.6, RCP4.5, RCP8.5 and High-End (i.e., the RCP8.5 together with increased ice-sheet contribution) scenarios were obtained from the studies of Jevrejeva et al. [70] and Jackson and Jevrejeva [71]. These projections, distributed on a 1° × 1° raster grid, are related to the tide gauge measurements of the period 1969–2016 (see Figure 2) and include the contributions of ocean self-attraction and loading, ocean steric and dynamic effects, land-ice inputs from Greenland, Antarctic, glaciers and ice caps, land-water storage and glacial isostatic adjustment.

(b) The present mean Hs distribution along the Rhodes coastal zone was estimated using the TELEMAC-based operational model addressing wave action computation [72]. The previously mentioned software includes an open-source code for a third-generation, phase-averaged, directional, spectral wave model that simulates the evolution in space and time of the energy spectrum of the free sea surface elevation at any water depth, based on a triangular finite element mesh. Time series of spectral wave data, obtained from the CMEMS [73] database (covering the period February 2006–January 2020), were utilized as offshore input data for the simulation. In addition, the changes of the mean Hs by the end of the 21st century were roughly estimated from comparisons between the outcomes of the wave simulations driven by the Euro-Mediterranean Centre for Climate Change, reported in the Deliverable 4.2 of the EU project SOCLIMPACT [74].

(c) The current mean tidal range was obtained from the study of Tsimplis [75], while for the projection to 2100 under the RCP8.5 scenario (averaged over the eastern Mediterranean), the study of Vousdoukas et al. [76] was considered.

(d) The coastal lithology and large-scale geomorphological features were determined from the digitization of the geological map (1:50,000) published by Mutti et al. [37]. However, the resolution of the island’s coastal zone was increased to a cell size of 50 × 50 m, using the results of the extensive fieldwork performed for the aims of the Interreg project COAST (2007–2013) [77]. Additional information about the texture of recent sediments occurring in several beaches of Rhodes was derived from the studies of Pyokari [43], Verikiou-Papaspiridakou et al. [44] and Anagnostou et al. [45].

(e) For the calculation of the current coastal slope, the 0 to +5.5 m topographic interval was analyzed, based on the present digital elevation model (DEM) (grid interval of 5 m) provided by the Hellenic Cadastre [78] and applying the endpoint method proposed by Doran et al. [79]. An averaged slope was computed every 50 m by applying a 100 m wide Hanning window. The +5.5 m elevation was chosen as the upper limit of the island’s coastal zone, because this altitude actually represents the maximum wave run-up during weak storm events (class I; see [36]). The coastal slope for 2100 was estimated from the analysis of the +0.65 to +6.15 m topographic interval of the above-mentioned DEM, since according to the extreme RCP8.5 scenario [70,71], the mean global SLR by the end of this century will reach 0.65 m above the level relative to 2020.

(f) The modern shoreline displacement rate was estimated from the differences between the shoreline positions of 2000 and 2018, identified via the analysis and interpretation of a Landsat-7 level 1T image of 22 October 2000 and a Landsat-8 level 1T image of 1 November 2018. The satellite images were downloaded from the official website of the US Geological Survey [80] at a 30 m resolution and a Matlab algorithm was generated for each shoreline extraction at the pixel level, based on Pardo-Pascual et al.’s [81] methodology. The normal differential water index was firstly applied, and then through spline interpolation, a pixel resolution of 7.5 m was accomplished. The produced final shoreline raster image was cropped to display only Rhodes, and then binarized by applying the best threshold. A morphological dilation filter was applied over the land and the binarized image was subtracted from the dilated one, thus obtaining the shoreline with pixel-scale detail. Eventually, the two shorelines (of 2000 and 2018) were processed with the DSAS software to estimate the mean erosion/accretion rates. For the projection to 2100, the extreme RCP8.5 scenario was adopted and the shoreline displacement rate was determined applying the widely used cross-shore (1D) parametric rule of Bruun [82]. The shoreline retreat (R_∞) due to the predicted SLR was calculated by the following equation:

R_∞ = [L_hc/(B + hc)] × S,

(2)

where S (m) is the relative SLR, L_hc (m) is the width of the active profile (distance from the subaerial berm to the closure depth), hc (m) is the closure depth and, finally, Β (m) is the subaerial berm height. The Bruun rule was chosen over other 2D models, e.g., Dean [83] and Kriebel and Dean [84], because it is (i) more suitable for estimations of the equilibrium shoreline retreat on beaches, where there are not any natural or man-made obstacles or constructions at their landward end (e.g., dune fields and coastal walls), (ii) easy to use as the granulometry of coastal sediments is not taken into consideration in the relevant equation and (iii) widely used, even in recent studies conducted on a global, regional or local scale [1,85,86,87].

(g) Τhe Holocene and historical vertical land movements due to tectonic activity were estimated from a plethora of reliable and detailed investigations carried out in Rhodes over the last five decades [30,31,32,88,89,90,91]. The potential coastal uplift for the next decades was based on the fault models applied by Kontogianni et al. [31] and Howell et al. [32], while for the probabilistic seismic hazard assessment in the area, the study of Woessner et al. [92] was considered.

4. Results

4.1. Current and Future CVI Variables

4.1.1. Mean Absolute Sea Level Rise

The combined effects of meltwater discharges into the ocean from glaciers, ice-sheets and ice cups, and seawater thermal expansion due to the prolonged increase of the upper ocean heat content, are the major reasons for the tangible modern SLR in several places around the globe [18,19]. Between 1993 and 2018, the thermal expansion of the oceans contributed to SLR by 42%, while the melting of temperate glaciers by 21%, Greenland ice-sheet by 15% and Antarctic ice-sheet by 8% [93].

However, the sea level change at regional and local scales can diverge substantially from the global SLR because of the natural variability in the strength of winds and ocean currents. In the southeastern Aegean Sea, the mean rate of the absolute SLR over the last 6 ky, based on Lambeck’s model [59,60], is estimated at 0.8 mm y⁻¹. The previous rate reduces to 0.7, 0.5 and 0.4 mm y⁻¹ over the last 3, 2 and 1 ky, respectively.

In contrast, for the 1969–2010 time interval, data recorded from tide gauge stations in the Aegean Sea (Thessaloniki, Alexandroupolis and Leros stations), corrected for vertical land motion, indicate SLR rates of 1.9–2.2 mm y⁻¹ [65]. Likewise, the study of Woppelmann and Marcos [66], carried out on the variability and trends of sea level dynamics in the Mediterranean Sea, suggests for the central and southern Aegean Sea absolute SLR rates of 2.43 and 3 mm y⁻¹ over the period 1993–2012, based on tide gauge measurements (Syros, Chios, Leros, Souda and Kalamai stations) and satellite altimetry data, respectively. The previous rates, however, decrease by 15% when the data of the 2010–2011 time interval are not considered, since during this particular period large positive sea level anomalies (up to 10 cm above the average climatological values [64]) were observed in the eastern Mediterranean, caused by the massive entrance of Atlantic seawater through the Strait of Gibraltar due to strong and long-term wind forcing [66].

However, an even more recent analysis of the available tide gauge measurements (Syros, Chios and Leros stations) and satellite altimetry data in the central Aegean Sea, spanning 1993–2017 [76], provides absolute SLR estimates of 1.77 ± 0.56 and 3.68 ± 0.33 mm y⁻¹, respectively. In addition, Tsimplis et al. [64], analyzing the measurements of the Rhodes tide gauge spanning 1969–2008, suggest a SLR of 0.103 ± 0.405 mm y⁻¹ plus the glacio-isostatic adjustment effect of 0.415 ± 0.19 mm y⁻¹. For comparison with the previous study, our analysis of the same datasets, enriched, however, with the recordings of the period 2010–2016 (see Section 3), indicates an absolute SLR of 0.72 ± 0.24 mm y⁻¹ (see Figure 2), reduced to 0.55 ± 0.21 mm y⁻¹ when the inverse barometer effect due to the low-frequency response of atmospheric pressure is subtracted from the observed sea level.

The absolute SLR rates derived from the analysis of the tide gauge measurements by the previous two studies differ significantly from the SLR of 3.42 ± 0.37 mm y⁻¹ estimated from the CMEMS satellite altimetry (see Section 3). However, such differences are not unusual and have been identified in many Mediterranean tide gauge stations [94], being attributed to (i) the high isostatic uplift of the southeastern edge of the Aegean microplate (which is an effect that has already been detected nearby the islands of Crete [95] and Leros [65]), and (ii) the high uncertainty of the near-shore estimates from satellite altimetry due to the contamination of the radar altimeter signal by the adjacent land [94,96].

Therefore, based on all the above information, the analyses of the tide gauge measurements by Vousdoukas et al. [76], Tsimplis et al. [64] and the current study were regarded as the most reliable sources for the estimation of the current SLR. Hence, an absolute SLR rate lower than 1.8 mm y⁻¹ was considered as a realistic condition for the Rhodes coasts and, therefore, taking into account the threshold values presented in Table 1, the lowest vulnerability score (1: very low) can be assigned to the entire coastline of the island (Figure 3).

Since the global warming, estimated from the present predictive scenarios, is anticipated to be intensified during the next decades, the current SLR shall be dramatically enhanced, thus potentially increasing the coastal vulnerability level. According to the RCP2.6, RCP4.5, RCP8.5 and High-End scenarios (see Figure 2), the mean SLR in the southeastern Aegean Sea, for the 2000–2100 time interval, has been estimated at 45, 53, 75 and 84 cm, respectively [70,71]. For all the previous sea level projections, the corresponding SLR rates are higher than 3.4 mm y⁻¹, which is the threshold of the highest vulnerability class for the mean SLR variable (see Table 1). Therefore, for 2100, the highest vulnerability score (5: very high) should be assigned to the whole coastline of Rhodes (Figure 3).

4.1.2. Mean Significant Wave Height

The spectral wave data obtained from CMEMS [73], covering the period 2006–2020, indicate mean Hs values offshore of Rhodes of 0.63 m to the north, 0.57–0.69 m to the northeast and east, 0.93–1.01 m to the south and southwest and 0.55–0.78 m to the west and northwest. Analogous mean Hs values have been provided by Zacharioudaki et al. [97] and Jadidoleslam et al. [98], performing wave hindcasting for a 42-year (1960–2001) and a 15-year (1999–2013) period, respectively. However, due to the effects of near-shore propagation and transformation of waves, the coastal areas of the island receive lower wave energy (as calculated by the TELEMAC-based wave model) represented by mean Hs values of 0.53 m to the northwest and north, 0.49–0.55 m to the northeast, east and southeast, and 0.74–0.95 m to the west, southwest and south (Figure 4a). Hence, based on the thresholds displayed in Table 1, the larger part (70%) of the Rhodes coastline (extending from Kritinia to Lachania) can be characterized by very low to low vulnerability (scores 1 and 2, respectively) to the current wave activity, while the remaining coastline to the southwest and south (extending from Kritinia to Prasonisi) appears moderately vulnerable to the presently incoming wave energy (Figure 4b).

Concerning the future changes in the intensity of the wave activity in Rhodes coasts, an overall attenuation has been predicted by various regional models (e.g., [99,100]). Nevertheless, this decrease is not expected to be uniform along the island’s coastal zone. In general, the reduction of the mean Hs is anticipated to be higher (10–21%) to the north and northwest coastline, and lower (3%) to the southern coastal segment (see Figure 4a,c). Hence, more than two thirds of the total Rhodes coastline, extending from Kritinia to Lachania, will be very lowly vulnerable by 2100, while the central-western coastal section from Kritinia to Siana will be less vulnerable (score 2) than the current condition (score 3: moderate) (see Figure 4b,d). The rest of the coastal area (from Prasonisi to Siana) will maintain its present moderate vulnerability level (see Figure 4b,d).

4.1.3. Mean Tidal Range

The tidal amplitude is strongly related to both permanent and episodic flooding events. Gornitz [23] suggests that a broad tidal range is associated with stronger tidal currents, which may cause erosion and sediment transport in such a way that macrotidal coasts (subject to a tidal range of >4 m) become more vulnerable than microtidal ones (subject to a tidal range of <2 m). Similarly, Shaw et al. [101], advancing the concept of Gornitz [23], suggested that an enhancement of the SLR would increase the flooding frequency of highly intertidal environments, such as deltas and estuaries, and this could also result in the inundation of areas being above the highest astronomical tide levels. Therefore, based on the perspectives of the previous studies, coastal areas that are subject to large-amplitude tides should be considered as highly vulnerable, while areas exposed to small tidal ranges should be considered as lowly vulnerable. However, Thieler and Hammar-Klose [49] provided a strikingly different concept, classifying the microtidal coasts as highly vulnerable and the macrotidal coasts as lowly vulnerable. The rationale is that on microtidal coasts, the sea level is always near the high tide level, thus, during the occurrence of a storm surge, the flooding on low-lying sandy/gravelly beaches is more likely than on macrotidal coasts [52]. On the other hand, at the time of ebb tide, the sea level on macrotidal coasts during storm surge events can be significantly lower than the high-tide level, thus increasing the probability of reduced flooding [21]. In addition, the microtidal environments are actually starving for sediments because the erosion processes take place in a narrow littoral zone. The concept of Thieler and Hammar-Klose [49] has been adopted by many researchers, such as Pendleton et al. [53] for the northern Gulf of Mexico, Dwarakish et al. [102] and Parthasarathy and Natesan [103] for the southwestern and southern coasts of India respectively, Hereher [104] and Pantusa et al. [105] for the Mediterranean Sea and Musa et al. [106] and Tano et al. [107] for the coast of western Africa.

For Rhodes, a current mean tidal range of 0.14 m [75] was assigned to all CVI cells, having, however, the disadvantage of missing the variations that could contribute to the identification of the most vulnerable coastal segments. Hence, following the abovem-entioned concept of Thieler and Hammar-Klose [49], the entire Rhodes coastline was assigned to the highest vulnerability score (Table 1, Figure 5). The adoption of the previous concept instead of that suggested by Shaw et al. [101] was based on the idea that the overestimation of the coastal vulnerability is wiser than its underestimation, while it should be taken into account that microtidal coastlines are particularly sensitive to SLR [108], with their coastal ecosystems (e.g., wetlands) exhibiting significantly low resilience to water level changes [109]. Indeed, a major segment of the eastern coastal zone of Rhodes comprises alluvial wetlands (see Figure 6a).

Finally, the mean tidal range by the end of this century is not likely to substantially change, because all probabilistic projections under the RCP2.6, R4.5, RCP8.5 and High-End scenarios suggest tide amplitudes not exceeding 1 m [76]. Therefore, the relative vulnerability score will most probably remain similar to the present one.

4.1.4. Lithology and Geomorphology

Previous studies [43,44,110] and our extensive field surveys spanning 2007–2013 have identified that fluvial–alluvial environments as well as sandy (and gravelly) beaches under erosion constitute more than 60% of the Rhodes coastal zone (Figure 6a). The beaches protected by beach rocks are limited, because the majority of these hard formations are uplifted [43,111] and, thus, they cannot substantially hinder coastal erosion. Rocky coasts of alpine substrate are mainly located in the western part (Kritinia–Monolithos area) and southern tip (Prasonisi) of the island, while post-alpine substrate occurs in the eastern and southern coastal segments (e.g., Koskinou, Archangelos, Lindos, Kattavia), mostly as headlands between extended sandy beaches (Figure 6a). Based on the descriptive thresholds of Table 1, 60% of the island’s western coastal zone, i.e., the Ialysos–Kritinia and Apolakkia–Kattavia areas, consists of loose sediments and may be ranked as very highly vulnerable (see Figure 6a,b). Further, 50% of the eastern coastal zone may be characterized as moderately vulnerable (Figure 6b) due to the presence of fluvial–alluvial deposits (Figure 6a) developed in the lowermost parts of the island’s hydrographic network (Figure 1a,b). Finally, very low vulnerability (37% of the frequency distribution) is assigned to most of the central-eastern (Archangelos–Lindos area) and central-western (Kritinia–Monolithos area) coasts, the larger part of the southern Kattavia–Lachania area and the major section of the northeastern coastal zone located between Rodos City and Afantou Bay (Figure 6b). The minimum vulnerability of the previous coastal segments is due to the occurrence of rocky substrates (alpine and post-alpine) and coarse-grained fluvial–alluvial deposits (pebbles, gravels and boulders), as well as the construction of coastal infrastructure and defenses in Rodos City, such as ports, marinas, rip raps and retaining walls (Figure 6a).

The lithology and major geomorphological features of the Rhodes coastal zone are expected to be only slightly modified until the end of the 21st century, since the current configuration has been evolving through long-lasting processes [110]. Changes in various localities are likely to occur due to erosion/accretion processes affecting a variety of beaches and sometimes causing the exposure of rocky substrates. However, in most cases, these changes will be rather ephemeral and locally restricted [110], thus not impacting the existing coastal geology. Therefore, the score of this CVI variable will rather remain consistent over the next decades.

4.1.5. Coastal Slope

Coasts ranked as highly (score 4) and very highly vulnerable (see Table 1) occur at 12% and 2% respectively (Figure 7a), with the coasts of moderate, low and very low vulnerability representing 20%, 25% and 41% respectively, of the frequency distribution. High to very high vulnerability is mainly suggested for some of the northwestern beaches of Rhodes (e.g., Ialysos, Paradeision, soroni), the eastern beaches of Faliraki, Afantou and Kalathos, and the southern beaches of Lachania and Kattavia, all of them being popular tourist attractions. The eastern part of Rodos City mainly displays very low, low and moderate vulnerability due to the existing infrastructure, such as ports and marinas. Finally, the larger portion of the central-western and central-eastern coastal segments of the island (Fanes–Monolithos and Archangelos–Lindos areas, respectively) is structured by cliffs with inclinations ranging from 7% to 20%, thus ranked as lowly vulnerable, or has gradients higher than 20%, thus assigned to a very low vulnerability level (Table 1, Figure 7a).

The vulnerability scores along the Rhodes coastal zone due to the coastal slope configured by 2100, compared with the current situation, are expected to only marginally change (see Figure 7a,b) according to the resolution of the spatial analysis (50 m; see Section 3). A minimal increase of the very highly (0.1%), highly (0.1%) and lowly (0.2%) vulnerable areas as well as a minor decrease by 0.4% of the moderately vulnerable areas are estimated to occur by the end of the century. The slope of the rocky coasts in the next decades should remain analogous to the current condition. Small variations in the low–lying beaches cannot be excluded, but they will hardly affect the future vulnerability level of the island’s coastline, at least at the spatial scale of the present analysis.

4.1.6. Mean Shoreline Change

The analysis of the evolution of Rhodes shoreline for the period 2000–2018 reveals mean erosion/accretion rates ranging from –1 to +1 m y^–1, thus indicating moderate vulnerability for almost the whole coastline length (Table 1, Figure 8a). However, sporadic beaches of high to very high vulnerability are identified in low-inclination (<4%) coastal segments, such as the southern Apolakkia–Prasonisi area, the northwestern Ialysos–Fanes area and the Rodos City–Kalathos area [112] (Figure 8a). Concerning the beautiful and greatly popular northeastern beaches of Faliraki, Afantou and Archangelos, seasonal shoreline variations are a usual phenomenon over the years [110], with the sediment budget, however, being in balance until now. Hence, during spring and summer, the prevailing wave regime establishes longshore sediment transport in a clockwise direction towards the north, accumulating sediments on the northeastern coastline of the island, while during winter, the domination of easterly and southeasterly winds in the eastern coastal zone triggers intense wave activity that drives the longshore sediment movement in a counterclockwise direction back to the western coastline, removing priceless sediment from the previous northeastern beaches.

On the other hand, some low to very low vulnerability beaches, configured by mean accretion rates of <1 m y⁻¹, are identified in the Monolithos, Siana and Kritinia coasts, Kalathos–Lindos area and Prasonisi islet (Figure 8a). Finally, the shoreline appears practically unchangeable in the Malona, Masari and Gennadion beaches.

It should be emphasized that for the accomplishment of a more representative evaluation of the current vulnerability due to the shoreline variations, the approximate width of the beaches was also taken into account, since the intensity of shoreline retreat or accretion does not have the same practical effect in a narrow or a wide beach.

Based on the outcome from the Bruun model [82], under the extreme RCP8.5 scenario, the SLR rise by 2100 will result in the retreat of the entire shoreline of Rhodes. In particular, maximum retreats, up to 38 m (i.e., 0.5 m y⁻¹), are predicted for the southern part of the island (Lachania–Kattavia area), corresponding to a loss of 54% of the presently occurring area. A smaller retreat of 18 m (i.e., ~0.2 m y⁻¹) is estimated for the northwestern part of the island, between Rodos City and Soroni, which due to the limited width of the existing beaches will result in an area reduction by 80%. In the southeastern (e.g., Gennadion) and southwestern (e.g., Apolakkia) coastal segments of the island, the landward shoreline displacement is estimated at 25 m (0.3 m y⁻¹), resulting in their reduction by 52% and 63%, respectively. The two most famous beaches of the island, i.e., Afantou and Faliraki, will suffer a 20 m shoreline retreat (0.3 m y⁻¹), which will decrease their width by 46%, and this shall negatively affect the local economy. Finally, the attractive Kalathos beach will be reduced by 33% due to a landward shoreline movement of 17 m (0.2 m y⁻¹). All aforementioned shoreline displacements have been estimated without taking into consideration future coastal protection works or a significant alteration (albeit unlikely) of the local wind/wave regime during the next 80 years. Even though the retreat of the entire shoreline of Rhodes by the end of this century seems to be inevitable and remarkable, all estimated displacement rates lie within the moderate vulnerability class (i.e., –1 to +1 m y⁻¹) (Figure 8b), which actually reflects the present-day condition.

4.1.7. Vertical Land Movements

Rhodes is an uplifted landmass lying on the Hellenic forearc, where chaotic tectonic displacements take place, deforming the entire region and causing its unceasing elevation with rates of up to 3.3 mm y⁻¹ over the last 5 ky [61]. However, the uplift of Rhodes is not uniform, demonstrating a variability that allows the division of the island into coastal segments characterized by distinct vertical motion rates [30,32,37]. Among the studies taken into account for the calculation of the mean land uplift along the Rhodes coastal zone, the investigations of Pirazzoli [90], Pirazzoli et al. [30,89], Howell et al. [31] and Pavlopoulos et al. [113] provided the majority of the relevant data. Hence, based on the calculated uplift rates, the Rhodes coastal zone may be divided into the areas illustrated in Figure 9a. The mean ground uplifts range from 3.84 mm y⁻¹ at the north (Rodos City–Koskinou area) to negligible values at the south, while the eastern coastal zone demonstrates a more intensive vertical displacement trend (0.38–0.78 mm y⁻¹) compared to the western coastal zone (0.15–0.17 mm y⁻¹) (see Figure 9a). Hence, according to the thresholds provided in Table 1, the coastal vulnerability of the northern part of Rhodes (Ialysos–Koskinou area) may be ranked as moderate (representing 12% of the frequency distribution), while the rest of the coastline can be considered as very highly vulnerable (Figure 9b). The rationale for the vulnerability ranking of an uplifting coastal zone is supported by the fact that coasts subject to strong land uplifts are usually quite accretive [70].

The projection of the magnitude of the vertical ground motions to 2100 is not possible to evaluate due to the chaotic regional geodynamic setting. The area offshore of Rhodes has a very high seismic potential, which may generate, at any time, strong (M > 7.5) earthquakes, causing differential uplift and tilting of the island’s crustal blocks [31]. However, if paroxysmal seismic phases do not occur in the next decades, then the rate and trend of the vertical displacement of the island’s coastal zone will most probably be similar to the current condition.

4.2. Asssessment of the Current CVI for Rhodes and the Projection to 2100

The variation of the current CVI along the Rhodes coastline is illustrated in Figure 10a. Based on the CVI’s frequency distribution, the highly and very highly vulnerable coasts correspond to 18% and 22% respectively, while the very lowly, lowly and moderately vulnerable coasts correspond to 7%, 32% and 21%, respectively. In particular, the majority of the beaches located in the northwestern coastal section of the island, extending from Ialysos to Kritinia, are assessed as highly to very highly vulnerable (Figure 11a–c) due to the occurrence of low gradients (mostly < 4%), loose (sandy or gravelly) sediments and shoreline retreat rates of 0.3 m y⁻¹. In contrast, southeast of Rodos City up to the Koskinou medieval village, very low to moderate vulnerability can be assessed, mainly due to the presence of coastal infrastructure such as ports or marinas. Further south along the eastern coastline, the Faliraki (Figure 11d) and Afantou bays, including the most visited beaches of the island, are assessed as highly to very highly vulnerable, while the segment stretching from Archangelos to Asklipieion demonstrates a greatly variable vulnerability due to the alternation of rocky cliff coasts and alluvial plains (see Figure 6a). The Asklipieion–Apolakkia area, extending from southeast to southwest, comprising several top-rated tourist destinations, appears highly to very highly vulnerable, except the segment near the crowded during summer Prasonsi peninsula/islet, which mainly shows low vulnerability. Finally, the central-western coastline (Monolithos–Kritinia area) is mostly characterized as lowly vulnerable, due to the existence of steep topographic relief (>20%) and rocky alpine substrate.

Concerning the evolution of the CVI up to the end of this century, we superimposed the CVI distribution of 2100 on the present situation, and the results are demonstrated in Figure 10b. It is evident that the vulnerability level of the major part of the Rhodes western coastline will be severely enhanced (48% of the total evaluated cells indicate higher vulnerability compared to the current condition), while for the majority of the eastern coastline, a vulnerability level analogous to the present one (mainly high to very high) will be established. Only 3% of the total cells composing the entire coastline display a trend for lower vulnerability.

The physical variables contributing to the CVI configuration are not independent, thus, they do not equivalently affect the variance of the unweighted index. Four of them, i.e., rates of SLR and shoreline displacement, and tidal range, as well as coastal vertical velocity (along most of the coastline), are uniformly distributed in the study area (see Figure 3, Figure 5, Figure 8 and Figure 9b) and, therefore, they do not affect the variance in both current and future CVI distributions. Concerning the rest of the variables, the combination of geomorphology and coastal slope essentially influences the current vulnerability of the entire coastline, while the significant wave height plays a secondary role, mainly affecting the variation in the present vulnerability of the southwestern and central-western coasts (compare Figure 4b, Figure 6b and Figure 7a with Figure 10a). Likewise, regarding the CVI of 2100, the increase in the vulnerability of the western coasts should be due to both geomorphology and coastal slope, with the impact of the reduced significant wave height in relation to the current condition (see Figure 4b,d) being minimal.

5. Discussion

Alternative methods used to assess coastal vulnerability to climate change include GIS-based decision support systems as well as techniques based on the application of dynamic computer models [51]. Concerning the former approach, the vulnerability-to-risk assessment is a key component, and an integrated regional risk assessment (RRA) methodology has been developed. For the performance of RRA, a chain of numerical models has to be established simulating circulation and morpho-dynamic processes, ranging from models reproducing atmosphere and ocean dynamics as well as circulation, biogeochemical, fate and transport processes in the coastal environment [51]. However, the building of a multi-model chain requires great experience and hard effort in terms of time and resources to make the approach applicable for a large area. On the other hand, regarding the latter approach, even though the use of dynamic computer models represent important tools for the analysis and mapping of the vulnerability and risks of coastal systems due to climate change, they present complexities arising from issues such as need of high expertise for their running, provision of results not always comprehensible to the decision-makers, provision of results of low spatial resolution, need of intensive testing and hard validation effort, incorporation of a high number of variables, need of very detailed site-specific data, focus mainly on the vulnerability of particular coastal sites (e.g., wetlands) and, finally, high cost [51].

In contrast to the above-mentioned methods, the CVI approach still remains a convenient management tool for shedding light on the potential susceptibility of coasts against a number of climatic and geological hazards. The CVI can efficiently simplify and assess the impact of fundamental physical factors or processes (acting at different scales), and otherwise should demand sophisticated modeling for an evaluation of the induced coastal response. The determination of CVI allows the identification of those coastal hotspots that are potentially jeopardized by climatic hazards over regional scales, which is vital for an effective long-term and large-scale coastal zone planning in order to provide coastal sustainability and successful adaptation strategies [21].

Even though the current absolute SLR in the southeastern Aegean Sea is estimated to occur at a relatively slow rate, i.e., <1.8 mm y⁻¹ (averaged over the last decades), the determination of the current CVI distribution in Rhodes strikingly indicates that 40% of the island’s coastline is highly to very highly vulnerable. We believe that there are two reasons that could explain this finding. The first one should be the formation of a plethora of long, low-gradient fluvial–alluvial/sandy beaches [43,44] in a microtidal setting, while the second one should be related to the sediment supply and coastal sediment transport pathways. Dams and coastal infrastructure (e.g., highways, buildings and walls), intervening with the fluvial network in both the eastern and western watersheds, retain sediment upstream and deprive the beaches of vital coarse-grained material. In addition, despite the fact that the sediment transport trends are complicated and seasonally and locally differentiated, a permanent longshore drift occurs on either side of the island, progressively removing sediment to the north [43,45]. Besides, during extreme storm surge/wave events, cross-shore sediment transport processes shift beach material deeper than the closure depth, which eventually cannot return to the shore in due course during the beach recovery phases [36,112]. Therefore, the majority of the beaches in Rhodes suffer from gradual erosion due to long-term sediment loss, which subsequently configures shoreline retreat conditions. Indeed, field surveys carried out by Anagnostou et al. [45,110] revealed two-decade shoreline landward displacements up to ~1 m y⁻¹ in the Faliraki and Afantou bays, and the Ialysos–Kritinia area, attributed to the erosion processes and not to SLR. The previous studies are in agreement with the results of the present study based on satellite image analysis and explain the shoreline retreat, albeit the very slow SLR that is actually confirmed by the records (since the 1970s) of the tide gauge in Rodos City.

Between the variables that configure the current CVI distribution in Rhodes, the uplift of the island is a crucial factor that counteracts the pressure exerted by the eustatic SLR. Many stratigraphic data derived from marine sedimentary basins and coastal terraces, found nowadays on the Rhodes landmass, have suggested an overall uplift of more than 500 m since the beginning of Quaternary [32]. This long-term ground elevation, albeit not uniform throughout Rhodes, seems to be more effective on the evolution of the island’s coastal zone than the glacio-eustatic cycles. Uplift rates as high as 12 mm y⁻¹ have been estimated for the broader setting of the fragmented forearc of the Hellenic subduction zone [44,88,95,114,115] as a result of the inconsistent slide (rates of <5–35 mm y⁻¹; [116]) of the Nubian (north Africa) plate below the Aegean (south Eurasia) microplate. In general, in the active plate margins, the vertical ground movement is the dominant process for the configuration of the coastal zones [117,118,119], and therefore, it has to be considered as a fundamental CVI variable.

In addition to the uncertainties in the determination of the current CVI for Rhodes, the estimation of the corresponding CVI for the upcoming decades was based on more precarious ‘assumptions’, derived from disparate and doubtful data sources and characteristics that may have significantly affected the reliability of the final scores. Dissimilar information of diverse spatial and temporal resolution, obtained from an amalgam of field surveys, remote sensing analyses and outputs of numerical models, has probably decreased the credibility of the CVI projected to 2100 [120].

The principal uncertainty, which subsequently could have an impact on the rest of the CVI variables projected to the future, is the foreseen rate of the absolute SLR. The IPCC predictions are revised intermittently by the consideration of new, more reliable, higher-quality and better spatially distributed climatic and oceanographic measurements. The new informal estimates imply a higher global SLR of more than 1 m by 2100 and 5 m by 2300, if the global targets on anthropogenic greenhouse gas emissions are not achieved. For the case of Rhodes, all predictions until now suggest a SLR rate higher than 3.4 mm y^–1 by the end of the 21st century. Therefore, the assignment of this CVI variable to the ‘very high’ vulnerability class is rather reasonable and realistic.

The results of wave models applied to the entire Mediterranean basin indicate an initial increase of the mean Hs by the first half of the century, and then a considerable attenuation by 2100 [99]. However, this trend should vary with respect to local peculiarities of the atmospheric circulation and wind fields, as identified in the northwestern Mediterranean Sea [100], Adriatic Sea [121] and Ionian and Aegean seas [122]. According to the EU project SOCLIMPACT [74], around Rhodes, the mean Hs by 2100 is predicted to diminish by 13% in the northern part of the island and by 2% in the southern part. However, Makris et al. [122], applying storm surge and wave models, suggest that an increase of the mean Hs should be expected at the Crete–Karpathos–Rhodes straits (up to 0.2 m locally) by 2100 due to the corresponding increase of the Etesian winds in the central Aegean Sea. The previous contrasting estimates emphasize the high uncertainty of the mean Hs projection to the future. Nevertheless, the present study follows the projection of the EU project SOCLIMPACT [74] because this is in line with the IPCC’s conclusions for the Mediterranean Sea [18,19].

The future shoreline change is another CVI variable with enhanced uncertainty. For the aims of this study, the shoreline response to the future SLR was estimated from the Bruun rule, because the majority of the Rhodes coastal zone comprises relatively low-energy, gently sloping, sandy beaches that are slightly modified by humans [123]. The rate of the predicted beach loss due to the future SLR is not going to exceed the upper threshold of the ‘moderate’ vulnerability class (i.e., −1 m y⁻¹), thus remaining at a level analogous to that of the period 2000–2018. However, the upcoming shoreline retreat could be locally accelerated by sudden and catastrophic events, such as earthquake-induced coastal sediment instabilities (e.g., liquefaction and mass wasting events) and storm-induced coastal erosion, which are presently common in Rhodes [45,124].

Finally, the future vertical movement of the Rhodes coastal areas will mainly be triggered by earthquakes or post-seismic deformation, and will vary from centimeters to several meters, primarily depending on the energy released from the event, the distance of the epicenter form the coast, the type (normal, reverse, strike–slip) and geometry (length, width, depth, orientation, inclination) of the fault and the elastic properties of the medium of the seismic wave propagation [31,125]. Rhodes, located on the active transtensional Hellenic forearc, currently undergoes tectonic uplift caused by (i) infrequent strong earthquakes, (ii) frequent small earthquakes and (iii) continuous gradual uplift due to sedimentary underplating [30,32,126]. The northern, northeastern and eastern parts of the island are more prone to tectonic uplift since they face the active Pliny–Strabo fault zone. Howell et al. [32], applying a variety of tectonically plausible fault models, point out that a single large (M ≥ 7.7) earthquake initiated by reverse faulting can cause in this region a coastal uplift of up to 3 m. However, the previous uplift value cannot be projected to the future, because such devastating earthquakes are rare, and the time of their occurrence is unknown. Therefore, based on the previous explanation, the long-term well-known mean uplift rate, spanning Holocene to today, is reasonably used for the determination of the future CVI distribution throughout the Rhodes coastline. In general, we believe that the tectonic uplift is a crucial factor that can mitigate the coastal inundation in the next decades. In particular, the tectonic uplift of the northern part of the island will most probably reduce the coastal vulnerability to the ongoing SLR by 94%, 74%, 50% and 43% under the RCP2.6, RCP4.5, RCP8.5 and High-End scenarios, respectively.

Considering the above-described uncertainties associated with both current and future CVI determinations, a strict comparison between them is not an ideal approach. Nevertheless, a rough estimation of the CVI evolution up to 2100 can be derived, providing an indication, albeit ambiguous, of the vulnerability condition of the Rhodes coastline in the next decades. Hence, recalling the results already presented in Section 4.2, the vulnerability of the major part of the Rhodes western coastline will be severely enhanced (by 48%), while the majority of the eastern coastline will be under a vulnerability status analogous to the current one. However, it should be emphasized that any changes in the predictions concerning the mean Hs and the gradual uplift of Rhodes coasts, which are actually the factors that may efficiently counterbalance the accelerated SLR, can radically improve or deteriorate the future sustainability of the island’s mesmerizing coastal zone.

The current and future CVI determinations for Rhodes can contribute substantially to an initial planning for the preservation of the high environmental quality of the island’s precious coasts [4,127], where the majority of the local socio-economic activities are concentrated (according to the Hellenic statistical authority; https://www.statistics.gr; accessed on 26 July 2021). The most vulnerable tourist coastal segments to the current SLR appear to be the areas featured by topographic gradients of less than 4% and loose sediments (mostly including sand with granules and pebbles), which are exposed to a mean Hs varying from 0.53 to 0.95 m, such as the northwestern Ialysos–Kritinia area, the southwestern Apolakkia–Kattavia area, the Kalathos, Lachania and Gennadion beaches, and the fascinating Faliraki and Afantou bays (Figure 10a).

In order to protect or at least mitigate the degradation of the most vulnerable tourist coastal hotspots of Rhodes due to the oncoming acceleration of SLR, focusing, in parallel, on the development of multi-seasonal tourism and enhancing the tourism businesses in the southwestern part of the island, which is characterized by less tourism-based economic growth (according to the Hellenic statistical authority in 2012; https://www.statistics.gr; accessed on 26 July 2021), standard coastal engineering solutions could be adopted. These could include the construction of seawalls, revetments, parallel or normal to the shore groins, and detached breakwaters as well as beach nourishment. However, it should be emphasized that the selection of the optimum solution is not a simple procedure, because it incorporates environmental, technical, economic, political, social and aesthetic limitations, and the applied solution should be necessarily considered as part of an integrated management approach [128]. For the optimum choice of near-shore engineering solutions, shoreline management plans and coastal defense practices (among the ones already implemented around the world and particularly in the Mediterranean region), the basic priority should be the harmony with the natural processes. For example, thoroughly studied and quite promising engineering interventions to reduce/minimize the high vulnerability of coastal localities in the Mediterranean region due to SLR are those referred to by: (i) Alexandrakis et al. [2], for the protection of the coastal front of Rethymno City (Cretan Sea) in Crete, including beach nourishment and submerged detached breakwaters as a complementary measure for the dissipation and attenuation of wave energy; (ii) Di Lucio et al. [5], for the control of the flooding episodes (Ionian Sea) in the Kaulonia municipality (Calabria, Italy), recommending the construction of a quite long and wide submerged breakwater; (iii) Torresan et al. [129], for the preservation of the northern Adriatic Sea beaches (Italy), suggesting extensive beach nourishment together with the construction of sea dikes in those sectors of the coastal zone that are easily exposed to severe flooding and erosion.

In addition, an alternative approach to the previously described standard coastal defense practices could be the implementation of measures mostly taking into account crisis management. The planning proposed for the spectacular and highly touristic Galicia coastal zone (Spain) [4], including (i) an effective flood management strategy for the safe tourist transportation from coastal destinations vulnerable to extreme inundation events and (ii) the implementation of adaptive measures such as resilient infrastructure [130] that can lower the damage cost in terms of insurance policy, tourist cancellations and rebuilding of accommodation and recreational facilities, is a representative example. Moreover, setback boundaries could be established along the most vulnerable coasts. In this way, infrastructure development very close to the beaches could be hindered, thus permitting shoreline changes to occur within an area where conflict between coastal development and coastal processes would be avoided [131].

Finally, the protection of the Posidonia oceanica meadows thriving along the major part of the western coastline of Rhodes (i.e., Ialysos–Kattavia area) is of vital importance for the highly to very highly vulnerable western coastal segments, because a potential deterioration of the seagrass habitats may cause an increase in the wave energy and seabed slope, which will further intensify beach erosion.

Overall, in order to improve/validate the assessment of the island’s coastal vulnerability, an investigation should be carried out, including (i) the analysis of extreme flooding episodes according to the requirements of the EU Directive 2007/60/EC for the assessment and management of flood risks [132], (ii) the detailed determination of the vertical ground motions based on the multi-temporal interferometric synthetic aperture radar (MT-InSAR) technique [133] and (iii) the incorporation of socioeconomic variables through a multi-criteria index system, evaluating the relevant risks due to the future SLR [2]. This shall greatly support the sustainability of tourism in coasts threatened by the direct and indirect negative impacts of SLR [134,135,136], which are considered international hotspots and are significant contributors to the local and national economy.

6. Conclusions

CVI applications are robust approaches, playing an essential role in the future coastal management strategy against the increasing threats posed by the current climate trend. However, it is vital to diminish the subjectivity in the computation of the CVI scores. If the vulnerability ranking adapted from databases of the Pacific or Atlantic coasts does not show some kind of analogy with the physiographic and hydrodynamic conditions in the study area, then this may result in misleading vulnerability level estimations. Hence, there is an urgent need to develop region- or, even better, site-specific vulnerability ranking classes.

Even though the present absolute SLR in Rhodes is minimal (<1.8 mm y⁻¹), the current CVI distribution indicates high to very high vulnerability for 40% of the entire coastline length. This cannot be easily ignored (regardless of the CVI determination uncertainties), bearing in mind the strong dependence of the island’s economy on the sustainability of coastal tourism. In fact, in the case of an accelerated SLR (as all present IPCC scenarios predict), the domination of very low coastal relief (<4%), the presence of extensive alluvial–fluvial beaches in the eastern coastal zone and the currently acting erosion processes (landward shoreline retreats up to ~1 m y⁻¹) in the sandy beaches (including granules and pebbles) of the northwestern and southwestern sectors of the island will most probably favor an increase in the vulnerability levels. Indeed, the CVI projection to 2100 indicates a remarkable increase (by 48%) in the vulnerability of the greater part of the western coastline. In contrast, at the same time, the eastern coastline demonstrates a future vulnerability level similar to the current one, without suffering from additional physical stress. However, the maintenance of the high to very high vulnerability for a number of popular eastern beaches cannot be neglected. Even though the vulnerability assessment for Rhodes appears alarming, the magnitude of the actual negative impact cannot be easily quantified due to the complex interaction of many factors, including sediment budget, oceanographic parameters and climatic variables.

For the protection of the most vulnerable coasts of Rhodes, the well-known Integrated Coastal Zone Management strategy [128,137] should be firstly considered as a basis. In addition to that, common and standard widely used hard engineering structures, such as seawalls, revetments, parallel or normal to the shore groins and detached breakwaters, could offer multiple benefits. Even though revetments are normally successful at reducing coastal erosion and holding the beach [130], Cooper and Pilkey [138] argue that these interventions, which are commonly adopted (due to their low maintenance costs) to counteract natural processes, often transfer the problem from one site to another. Alternatively, beach nourishment could be adopted to support the most vulnerable top-rated tourist attractions of Rhodes after, however, the performance of detailed numerical and large-scale physical modeling, while the relocation of inhabitants and tourism infrastructure and facilities in less vulnerable areas (taking into consideration all social implications) could be a wise management option. Finally, as a nature-based solution, the high conservation status of the seagrass meadows flourishing along the western coastline should be a major priority.