Sea level variability and coastal inundation over the northeastern Mediterranean Sea

ABSTRACT Sea level is a key element of global-scale climatic changes with significant coastal impacts. The northeastern Mediterranean Sea (Aegean, Ionian, and Cretan Seas: AICS) consists of extended low-lying coastal zones exposed to coastal inundation due to seawater level increases. The variability of the Sea Level Anomaly (SLA) in the AICS has been evaluated with the use of satellite-derived observations from 1993 to 2021. The general SLA trend over the AICS domain was 3.6 mm/year, which is higher than that of the overall Mediterranean Sea. The coastal inundation variability along the AICS coastline, specifically in 20 characteristic study areas, was investigated based on the long-term satellite-derived sea level observations, high-resolution land elevation data, and numerical simulations, conducted with a coastal inundation (CoastFLOOD) model. The coastal zones were categorized based on the frequency of occurrence and percentage coverage of flooding incidents during the 29-year period, while the interannual trends and the seasonality of inundated littoral and inland areas also showed intense spatial variation. High flooding levels occurred at coastal areas of either extended low-lying areas and/or high sea level elevations. The coastal zones were categorized based on a heuristic Coastal Inundation Hazard Index (CIHI: 1–5 ranking) that was determined by the occurrence frequency of severe floods (temporal component) and their impact area (spatial component). Our study findings provide vital information about the coastal inundation conditions and the respective hazard level of the AICS, useful to support enhanced coastal zone management for the improvement of coastal protection and decisions about coastal hazards’ mitigation.


Introduction
Sea level is a key element of global scale climatic changes with potentially significant coastal impacts (Cazenave et al. 2014), especially on the low-lying areas of the Mediterranean Sea, which is one of the most vulnerable regions to Sea Level Rise (SLR) worldwide (Adloff et al. 2015;Calafat, Chambers, and Tsimplis 2012;Church et al. 2013;Diffenbaugh and Giorgi 2012;Gualdi et al. 2013).The increase of the sea level is one of the most crucial drivers of natural hazards (floods, encroachment, erosion, etc.) putting intense pressure on the coastal zone (Bosom and Jiménez 2010;Nicholls 2002;Vousdoukas et al. 2017).The impact might become more intense during the 21 st century, depending on the prevailing future climatic scenario (Neumann et al. 2015;Vitousek et al. 2017).The sea level variability controls a) the inundation levels of the coastal zone, depending on the local topographic characteristics (Calafat et al. 2022;Denamiel and Vilibić 2022;Ferrarin et al. 2021), b) the shoreline retreat, coastal erosion, and land loss (Enríquez et al. 2019;Jiménez et al. 2017;Monioudi et al. 2017;Sharaan and Udo 2020;Spencer et al. 2016), and c) the salinity levels in estuaries and freshwater aquifers (Bijlsma et al. 1996;Krvavica and Ružić 2020;Shalem et al. 2015).The efficient assessment of the coastal zone risk under flooding conditions (Satta et al. 2017;Vousdoukas et al. 2016) is strongly associated to the accurate estimation of the sea level temporal and spatial variability and the respective past, current and future trends (Bevacqua et al. 2019;Hallegate et al. 2013;Lionello et al. 2021).
The water level characteristics are associated with three main environmental factors: a) the steric effect due to thermal and salinity changes (Adloff et al. 2015;Carillo et al. 2012), b) the addition of water mass (e.g.glacier ice melting and changes in ground water storage) (Galassi and Spada 2014;Jordà and Gomis 2013), and c) the atmospheric component (e.g.storm surge formation) (Androulidakis et al. 2015;Jordà and Gomis 2013;Jordà et al. 2012;Levitus et al. 2012;Makris et al. 2023;Marcos, Tsimplis, and Shaw 2009;Šepić et al. 2012, 2015).The North Atlantic Oscillation (NAO; Hurrell 1995) impact on the weather conditions varies between the eastern and western parts of the Mediterranean Sea (Lionello 2012) altering the respective influence on the sea level variability that is linked primarily to atmospheric pressure changes and secondarily to local wind field changes (Tsimplis et al. 2013).The extreme positive Sea Level Anomaly (SLA) that occurred in 2010 was mainly associated to the negative phase of the NAO (Haddad, Hassani, and Taibi 2013;Landerer and Volkov 2013) causing about 12 cm of basin rise in the Aegean Sea (Tsimplis et al. 2013) and the broader eastern Mediterranean (Mohamed and Skliris 2022).The meteorologically induced component of the sea level (storm surge), which is associated to the Sea Level Pressure (SLP) and winds, reveals strong spatial variability, with stronger inverse barometer effect in the Aegean Sea compared to the Adriatic and Ionian Seas, where storm surges are mainly determined by the wind state (Androulidakis et al. 2015(Androulidakis et al. , 2023;;Denamiel, Tojčić, and Vilibić 2021;Denamiel and Vilibić 2022;Krestenitis et al. 2011;Makris et al. 2016Makris et al. , 2023)).The Mediterranean cyclogenesis can be categorized (Flaounas et al. 2022) to (i) cyclones formed over the Atlantic ocean and intrude into the basin (61-85% of total tracks; Neu et al. 2013), (ii) to Medicanes that form in western and central Mediterranean (weak frequency but strong magnitude; Cavicchia, von Storch, and Gualdi 2014), (iii) to north African cyclones that form over the deserts, (iv) explosive cyclones or weather bombs that form over northern areas (Kouroutzoglou et al. 2011), (v) Vb cyclones that affect mainly the western basin (Messmer, Gómez-Navarro, and Raible 2015), and (vi) secondary cyclones that form in the periphery of existing cyclones (Ziv et al. (2015).The AICS domain is mainly affected by the first and the second type that may cause intensive coastal flooding due to storm surges (Androulidakis et al. 2023).
In the Mediterranean Sea, the recorded strong positive mean sea level trends, associated to climate change impacts, have been confirmed by long-term tide-gauge measurements (e.g.Adloff et al. 2018;Dieng et al. 2021;Marcos and Tsimplis 2008;Tsimplis et al. 2013) and especially by remote sensing data during the satellite altimetry era after 1993 (e.g.Bonaduce et al. 2016;Fenoglio-Marc 2002;Mohamed et al. 2019;Vigo, Garcia, and Chao 2005).Zanchettin et al. (2021) investigated both historic and future trends of Sea Level Rise (SLR) in Venice (Italy) finding an average shift of about + 2.76 ± 1.75 mm/year, between 1993 and 2019, estimated with tide-gauge data after the removal of the contribution of calculated land subsidence.The Mediterranean Sea is characterized by strong regional variability with large differences between the sub-basins, showing stronger interannual trends in the central region, in the Adriatic and along the Tunisian coasts (Bonaduce et al. 2016;Cid et al. 2016;Meðugorac et al. 2018).Mohamed and Skliris (2022), using satellite-derived SLA gridded data, showed that the total sea level trend between 1993 to 2019 is around 3.1 ± 0.61 mm/ year over the eastern Mediterranean domain.Herein, we expand the results of Mohamed and Skliris (2022) using the same satellite product, extended to 2021 (29years) and focused on the Aegean, Ionian and Cretan Seas (AICS), which are divided to seven sub-regions (Figure 1a), following the categorization introduced by Androulidakis and Krestenitis (2022).The AICS domain includes various urban settlements, environmentally protected areas, extensive touristic infrastructure containing a large number of recreational areas with sandy beaches, and many regions for several other activities (e.g.aquaculture, fisheries, navigational transportation, seaport commerce, athletics, etc.).The ultimate goal is to derive the coastal interannual SLA characteristics that will be used as boundary conditions (input) in the high-resolution coastal inundation assessment of characteristic lowlying areas (Figure 1b), located along the AICS coastline.
The worldwide socioeconomic and environmental impacts of SLR are significant.Desmet et al. (2021) have estimated the possible projected spatial shifts in economic activity until 2200 due to future SLR scenarios, including the reduction of global real Gross Domestic Product (GDP) and significant displacements of population from coastal localities.As a countermeasure, many researchers have introduced several approaches for an Integrated Coastal Zone Management (ICZM).For example, Clemente et al. (2022) have proposed a simplified risk assessment model for coastal built environments in terms of potential damages, direct and tangible, economic losses based on the Copernicus "Coastal Zones" database and exposure estimates by damage scenarios due to extremes of SLR in Italian littorals.The evaluation of estimations for future (projected) SLR is supported by a vast literature.For instance, Jeon et al. (2021) investigated sea level fingerprints and their relation to regional sea level changes, based on measured mass results from the Gravity Recovery And Climate Experiment (GRACE) and Argo float data (measuring thermal expansion), agreeing well with satellite altimetry.The most prominent natural hazard induced by episodic or long-term SLR is coastal flooding (shortterm severe events) and/or coastal inundation (more long-term or perpetual situations) of low-lying littoral floodplains with various significant implications (Alvarado-Aguilar, Jiménez, and Nicholls 2012;Aucelli et al. 2017;Krestenitis et al. 2011;Nicholls and Hoozemans 1996;Refaat and Eldeberky 2016;Reimann et al. 2018;Rizzo et al. 2022;Shaltout, Tonbol, and Omstedt 2015;Snoussi, Ouchani, and Niazi 2008;Vandelli et al. 2022;Wolff et al. 2018). Favaretto et al. (2021) characterized storm events (by sea level variability and mean sea level rise) for the sake of coastal flood hazard assessment in north-central Mediterranean exposed littorals (Adriatic Sea; Venetian lagoon).Except from robustly defining the sea level components, Kulp andStrauss (2016, 2019) discussed the necessity to minimize errors in Digital Elevation Models (DEM) to avoid underestimations of coastal vulnerability due to sea level-induced flooding.Kulp and Strauss (2017) traced the rapid escalation of exposure of populations and infrastructure to coastal flooding in US municipalities due to rising sea levels, assessed by LiDAR elevation data.Hauer et al. (2021) also assessed the exposure of US population to coastal flooding due to SLR, showcasing how the choice of zone results in differential spatiotemporal exposure estimates.
The scope of the study is to evaluate the interannual and spatial variability of the sea level over the AICS region and evaluate the respective impact on the coastal inundation characteristics along the topographically complex Greek coastline for 29 years .The main motivation is to assess the state and severity of coastal hazards, related to sea level increases, over the AICS.Coastal inundation estimations and trend analysis for the last three decades are based on the observed sea level increasing trends associated to changes in climatic patterns.The SLA differences among the seven AICS sub-regions are investigated based on a continuous long-term satellite dataset that covers the study region.The derived information will be used to estimate the inundation variability during the same 29-year period over the Greek coastal zone with the use of numerical simulations with CoastFLOOD model (Androulidakis et al. 2023;Makris et al. 2022Makris et al. , 2023;;Skoulikaris et al. 2021) and an updated high-resolution dataset of land elevation derived from a DEM that covers the entire Greek coastline.We identified twenty (20) low-lying coastal areas (Figure 1b; based on the Copernicus DEM of 25 m resolution) that have been impacted by intense flooding events in the past and expand the analysis with the use of the finest available land elevation data (5 m) at the 20 study areas.

Sea level observations
The altimeter satellite SLA data (sea surface height above the mean sea level), computed with respect to a twentyyear (1993-2012) mean, has been collected from the Copernicus Marine Service (CMEMS; https://data.marine.copernicus.eu/;accessed in June 2023; Landerer and Volkov 2013;Pujol and Larnicol 2005) to investigate the sea level variability over the AICS region and provide the boundary conditions of the coastal flooding simulations (see Section 2.4) covering a 29-year period (1/1/1993-31/ 12/2021).
The DT-2021 European product (SEALEVEL_EUR_PHY_L4_MY_008_068), used in the study, has a spatial resolution of 1/8° (approximately 13 km), and it is provided in a daily step, while the initial coverage of the dataset extends over all European Seas.The Level L4 gridded data are produced by merging observations of Topex/Poseidon, ERS1/2, Jason 1-2-3, Sentinel (3A/B and 6A), HaiYang-2A/B, Saral[−DP]/Altika, Cryosat-2, ENVISAT, and GFO altimetry missions.These satellite-derived SLA fields have been previously used to evaluate the sea level variability in the Mediterranean Sea (e.g.Bonaduce et al. 2016;Landerer and Volkov 2013;Mohamed and Skliris 2022;Mohamed et al. 2019).Pujol et al. (2023) validated the recent version of the gridded L4 SLA product against independent SARAL-DP/Altika along-track measurements.The error is around 4 cm in the Mediterranean Sea and increased in high variability areas (e.g.Gulf Stream; >6 cm).Moreover, further comparisons with tide gauges showed that the improved altimeter product DT-2021, used herein, is characterized by the reduction of variance of the differences between altimetry and field measurements ranges between 0.2% to more than 5% of the tide gauge signal.The average results of different missions are used to calculate the SLA L4 interpolated fields of the SEALEVEL_EUR_PHY_L4_MY_008_068 product; the daily mapping SLA error is also provided in the product for every grid point (see Supplementary material).All standard geophysical and environmental corrections have been applied to this product, including the Dynamic Atmospheric Correction (DAC; Landerer and Volkov 2013) to improve the representation of high-frequency (<20 days) atmospheric forcing considering both pressure and wind effects.More information and quality control of the SLA dataset can be found in the product user manual (Pujol et al. 2023; https://catalogue.marine.copernicus.eu/documents/QUID/CMEMS-SL-QUID-008-032-068.pdf,accessed in June 2023).
In this study, we perform an additional evaluation of the satellite altimetry product with field observations, covering 11 continuous years (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012), derived from in situ measurements by four tide gauges of the Hellenic Navy Hydrographic Service (HNHS; https:// www.hnhs.gr/en/,accessed in June 2023) in the Aegean and Ionian seas (Thessaloniki, Alexandroupolis, Chios, Lefkada; Figure 1a).To compare the satellite altimetry product with tide gauge datasets, the respective altimetry grid point that covers the location of each station was used.Moreover, the tide gauge sea level timeseries were also corrected based on the DAC product provided by AVISO (ftp.aviso.altimetry.fr,accessed in June 2023), following the methodology presented by Mohamed and Skliris (2022).Daily and monthly averages were derived from the original hourly field data and were compared to the respective daily and mean monthly timeseries derived from the satellite gridded fields (Section 3.1).

Digital elevation model data
Two sources of land elevation were used to detect and describe the low-lying coastal areas along the AICS coastline, serving also as the topography input in the inundation estimations and simulations (Section 2.4).The first dataset was derived from the EU-DEM v1.1, which covers the entire European domain with a 25 m spatial resolution, including the AICS region (Figure 1b).The EU-DEM v1.1.is distributed by the Copernicus Land Monitoring Service (CLMS; https:// land.copernicus.eu/,accessed in June 2023).We have used the EU-DEM v1.1 data to select 20 characteristic low-land study areas along the AICS coastal zone (Figure 1b) based also on the flooding events of these areas as reported in scientific literature and mass-media press (see Section 2.6).The second dataset of coastal elevation has higher resolution (5 m) and was derived by the available geospatial data from the DEM of the official Greek service for comprehensive recording of real estate and property's metes-andbounds, i.e. the Hellenic Cadastre (https://www.ktimatologio.gr/en,accessed in June 2023).This highresolution dataset was used to investigate the variability of the inundation levels over the 20 coastal areas that have been affected by intense flooding events in the last 3 decades.The geometric accuracy of the DEM is less than 0.70 m, while its absolute accuracy is less than 1.37 m with a 95% confidence level (Chrisafinos and Kavvadas 2016).This DEM product is highly efficient to perform numerical simulations of flooding over the Greek coastal zone (Xafoulis et al. 2023).Maps of land elevation of the 20 study areas, derived from the high-resolution DEM, are provided (together with respective flood cover maps; Section 4) as supplementary material (Figures S2-S20).

Atmospheric conditions
The prevailing meteorological conditions during the 29-year period were evaluated based on the ERA5 reanalysis dataset produced by the European Centre for Medium-Range Weather Forecasts (ECMWF).The hourly data have a spatial resolution of 0.25° and are freely distributed by Copernicus Climate Change Service (https://cds.climate.copernicus.eu,accessed in June 2023).The SLP fields and the two (zonal and meridional) components of the wind at 10 m height are used to investigate the effects of the atmospheric forcing on the variability of the SLA over the AICS domain.

Numerical simulations of coastal inundation
Herein, we implement modeling of coastal inundation induced by a slowly varying component of the nearshore sea level with CoastFLOOD (Makris et al. 2022(Makris et al. , 2023)).CoastFLOOD is an Aristotle University of Thessaloniki (AUTh) in-house numerical model, applied in littoral inundation hydrodynamic studies focused on local-scale selected areas of the Greek coastal zone (Androulidakis et al. 2023;Makris et al. 2022Makris et al. , 2023;;Skoulikaris et al. 2021).These study coastal regions may include a) urban areas with engineered waterfronts, ports and coastal structures, road, highway and railway networks, b) rural areas with agricultural farmlands, natural or wild vegetated fields, forests, bare lands, stony element lands, pastures and grasslands, and c) coastal water bodies such as estuaries, lagoons, torrents, river deltas, and natural beaches.The model is fed by the SLA dataset (Section 2.1) as hydraulic load boundary conditions (i.e.floodwater level) on the study areas' coastlines represented by the marginal dry cell of the computational domain during Still Water Level (SWL) conditions.Hence, it provides estimations of possible extended flooding of lowland coastal areas, rural plains, farmlands, urban sites, engineered spaces, deltaic regions, natural beaches, and estuaries along littoral stretches due to sea level increase, e.g.storm surges/tides (Androulidakis et al. 2023;Makris et al. 2023).Two (realistic and idealized) types of numerical simulations were developed to derive i) daily realistic estimations of coastal inundation covering the 29-year period, based on the satellite-derived SLA levels at each study area (jobs: 20 areas x 10,592 days), and ii) estimations of the extreme floods, based on four idealized scenarios of SLA (0.5, 1.0, 1.5 and 2.0 m; jobs: 4 scenarios x 20 areas) along the coastline boundary of each study area.The extreme values are based on crude thresholds of extreme storm surge levels or storm tides (surge + highest astronomical tides) as derived by Makris et al (2016) and Galiatsatou et al (2019Galiatsatou et al ( , 2023)).and Total Water Level (TWL) as defined by Galiatsatou et al (2019Galiatsatou et al ( , 2021)).along the Greek coastline.Detailed information (conceptual approach, equations, grid discretization techniques, numerical schemes, etc.) about CoastFLOOD model are thoroughly presented by Makris et al. (2023a), who also discuss the benefits and limitations of specific model applications in coastal lowlands; a concise summary is provided in the Appendix.

Statistical methods
The magnitude of the interannual trends was evaluated based on the Sen's Slope estimation (Sen 1968) that provides the annual change of the variable under investigation.This slope corresponds to the median of all the slopes calculated between each pair of points (x i , x j ) in the series: Sen's Slope is considered efficient to detect the linear relationship as it is not affected by outliers in the data (Ray et al. 2021).The significance, which represents the threshold for which the hypothesis that there is not trend is accepted, is determined based on the computation of the p value .A significance level of α = 1% was used in a two-tailed test, in which the null hypothesis is that there is no monotonic trend is the time series (H 0 : τ = 0; p value >α), and the alternative hypothesis is that there is a significant monotonic trend in the time series (H a : τ≠0; p value <α).Here, the Sen's Slopes and the respective p values were computed with the XLSTAT version 2018.1 for the timeseries and with the MATLAB toolbox "Mann-Kendall Taub with Sen's Method (version 1.15.0.2;Burkey 2023)" to produce horizontal maps.The detection of abrupt changes during specific points in the long-term timeseries was also assessed with XLSTAT toolbox, based on the Pettitt's homogeneity test (Pettitt 1979).This test detects shifts in the average and calculates their significance in a hypothesis test.The null hypothesis is that the data are homogeneous, as against the alternative hypothesis that there is a datum at which there is a change in the data.The empirical significance level (p value ; significance level of 1%) was performed to test the hypothesis of a shift in the mean in the dataset.If a shift exists (p value <0.01), the difference between the two means, before (μ1) and after (μ2) the shift, is also computed.Pettitt's test is a widespread nonparametric method to detect inhomogeneity, the change point in a time series (Bickici Arikan and Kahya 2019).The Pettit's test was mainly applied in this study to detect the changes of flooding levels during the 29-year period.
The correlation between timeseries was tested with the computation of the Pearson product-moment correlation coefficient (Pearson 1903); it is a value that varies between −1 and 1 measuring the strength and direction of the relationship between two timeseries.The statistical significance of the correlations was also derived based on the Mann-Kendall (Kendall 1975;Mann 1945) correlation test (p value ); if p value is lower than a conventional value (e.g.5%; p value <0.05), the correlation coefficient is statistically significant.The Willmott Skill Score (WSS) or Index of Agreement (Willmott, Robeson, and Matsuura 2012) is also used to evaluate the agreement between two timeseries; the higher the WSS (with ≤1 as a limit), the better match is reached between the timeseries.Another metric to define the potential need for bias-correction of satellitederived timeseries, compared to field observations, is the Hit-Rate-of-Percentiles (HRP) index (Makris et al. 2016;Schoetter et al. 2012).It calculates the differences between sorted percentiles of two timeseries, then defined as the sum of all categorical fractions (viz., differences compared to an allowed deviation, usually the standard deviation of observed timeseries).In general, if HRP index is greater than 0.95, then bias correction of derived SLA is not necessary.
The characteristics of SLA over each AICS subregion was also analyzed based on the development of heat maps.A heat map depicts values for a main variable of interest across two axis variables as a grid of colored squares.The axis variables are divided into ranges like a bar chart or histogram, and each cell's color indicates the value of the main variable in the corresponding cell range.It shows the relationships and potential patterns' similarities between two variables (e.g.SLA over each year and over each sub-region) one plotted on each axis.By observing how cell colors change across each axis (e.g. years in x-axis and sub-regions in y-axis), one can observe if there are any patterns in value for one of both variables (e.g.sampling of years and areas with similar patterns of SLA).Clustering of the heat maps is also applied to create dendrograms for both rows and columns.The dendrogram is a tree-structured graph used in heat maps to visualize the result of a hierarchical clustering calculation.Herein, the heat maps were derived with the use of the XLSTAT toolbox to investigate both interannual and seasonal variability of SLA.

Case study areas
The twenty (20) selected case study areas (and the aquatic basins they belong to), located along the AICS coastal zone, are presented in Figure 1b.These are generally established lowland coastal areas, based on the EU-DEM v1.1, that are well-known to be prone to flooding events.We chose them based on a series of recorded coastal inundation events due to storm surge (and/or compound flooding) in AICS littoral regions, which were reported rather either in scientific literature (Androulidakis et al. 2023;Diakakis, Deligiannakis, and Mavroulis 2011;Ferrarin et al. 2023;Ghionis et al. 2015;Kombiadou et al. 2012;Kontogianni et al. 2012;Krestenitis et al. 2011;Makris et al. 2023;Martzikos et al. 2021;Monioudi et al. 2016;Poulos et al. 2022;Skoulikaris et al. 2021;Stavropoulos et al. 2020;Tragaki, Gallousi, and Karymbalis 2018;Tsoukala et al. 2016) or in journal press, mass and social media.A list of such indicative events and the respective weather characteristics is provided in Table S1 (included as Supplementary Material).The extent of each study area was selected to include its most important environmental and socioeconomic features such as the coastal urban infrastructure, seaport facilities, environmentally protected areas (lagoons, wetlands, watersheds, etc.), sandy beaches, touristic facilities, holiday residencies, lowland agricultural areas, natural bare lands in the coastal plain, etc.For example, Area 1 (Thermaikos Gulf) includes Thessaloniki city at its northern coasts, which is the second largest city of Greece, touristic regions and recreational beaches in the east, river deltas and extensive agricultural farmlands along the west coasts.The 20 study areas are representative of the entire AICS coastal zone with respect to topographical, environmental, and social characteristics.

Interannual trends and spatial variability of Sea level Anomaly
Initially, we evaluated the satellite-derived SLA observations against field observations at 4 coastal areas (Section 3.1).The distribution of the SLA over the 7 sub-regions of the AICS (Figure 1a) was examined to detect the interannual trends (Section 3.2) and seasonal cycles (Section 3.3), highlighting the differences between the AICS sub-regions.Finally, the contribution of the two main atmospheric components (SLP and winds) on the SLA variability of each sub-region was also investigated to evaluate the spatial variability of air pressure and wind states over the AICS domain (Section 3.4).

Comparison between satellite and tide-gauge observations
The daily variability of the SLA at four locations of the AICS domain, derived from both satellite and tidegauge observations during an 11-year period, is presented in Figure 2. The Pearson correlation coefficients in most of the stations are higher than 0.80 (Table 1) and they are statistically significant based on the MK correlation test (p values <0.01: statistical significance level 99%).The respective monthly averages show even higher correlation in Thessaloniki and Chios areas (R P >0.90).The agreement between the two datasets is good in terms of percentiles too, showing quite high HRP-index (>0.99), reaching 1.00 in the Ionian Sea (Lefkada; Figure 2d), a region with significantly high storm surge frequency, especially due to severe tropical-like atmospheric cyclones (e.g.Medicanes; Androulidakis et al. 2023;Lagouvardos et al. 2022).These values signify that there is no imperative need for bias-correction of the SLA timeseries, at least for the present scope of feeding a large-scale coastal inundation model for flood hazard mapping.The Root Mean Square Error (RMSE) is between 2 to 5.5 cm given the difficulty to achieve agreement in point-to-point values' comparisons between gridded SLA values and in situ measurements.The ratio between RMSE and maximum SLA is less than 20% supporting the efficiency of the satellite-derived data to measure the high sea level peaks in the coastal zone.The WSS are significantly high, reaching close to unity (>0.90) in most of the cases (Table 1).The maximum SLA values (>0.40 m) were well recorded by both techniques, especially during the extended winter months (December -March) of 2009-2010 and 2010-2011, which were periods with high SLA values in the eastern Mediterranean Sea under negative NAO phase conditions (Bonaduce et al. 2016;Mohamed and Skliris 2022;Tsimplis et al. 2013).Differences between tide gauge measurements and satellite observations can also be related with the vertical land motion of coastal areas affected by geodynamic processes and landscape evolution due to tectonic movements, eustatic factors, etc.For example, in the coasts of central Ionian there is a slow land subsidence rate ≤2 mm/year (Anzidei et al. 2014).This may result to approximately 2 cm per    1).Near-null movements have been detected in the central and northern Aegean and even a slight uplift in Crete (Anzidei et al. 2014;Fenoglio-Marc, Dietz, and Groten 2004).
The SLA error associated to the transmission from the L3 to L4 product used herein is generally larger near the coast (up to 2.5 cm; supplementary Figure S1b) than offshore, in the open sea of the AICS region (less than 1.5 cm).It is noted that the averaged SLA error over the entire region reduced throughout years during the last three decades from 2.3 cm in 1993 to less than 1.6 cm in 2021 (supplementary Figure S1a) probably due to the several correction improvements and ongoing addition of extra satellite tracks during the 29-year period.Thus, the comparison results support the usage of the remote sensing dataset to assess the sea level interannual variability along the coastal zone of the AICS domain, especially during periods of SLA peaks.These serve as forcing conditions of the coastal flooding simulations (see Section 4).

Interannual variability
The interannual evolution of the mean annual SLA is characterized by clear positive trends in all AICS subregions over the 1993-2021 period (Figure 3a).The general trend over the entire area, based on the computation of the Sen's Slope, was 3.6 mm/year (Table 2; black line in Figure 3a), higher than the overall Mediterranean (2.44 mm/year; Bonaduce et al. 2016) and the eastern Mediterranean (3.1 mm/year; Mohamed and Skliris 2022) trends.All trends are statistically significant (p value <0.01;Table 2).A clear north-to-south increasing gradient of mean SLA was computed in the Aegean Sea (Figure 4a) with relatively low values south of Chalkidiki and over the north-eastern islands (<4 cm), while higher levels were found in the southern Aegean (>4.5 cm).The spatial distribution of the interannual SLA trends (Figure 4c) shows statistically significant (Figure 4d) strong Sen's Slopes over the Aegean Sea (>3.5 mm/year) and the northern Cretan Sea.The highest Sen's Slopes were derived for the SW Aegean (4.1 mm/year; Table 2 and Figure 3a); specifically, over its northernmost coastal area (south Attika and north-eastern Peloponnesus; Figure 4c), the interannual increase exceeded the value of 4 mm/year.The same area also revealed the highest mean SLA of the entire 29-year period (Figure 4a).Very strong trends were also computed along the eastern Aegean islands and Minor Asia coasts as well as the southern Dodecanese (Figure 4c).Weaker slopes of the annual mean SLA trends were computed offshore in the Ionian and southern Cretan Seas (Figure 4c), with even negative and statistically insignificant values (p values >0.05; Figure 4d) in the broader area of the Ierapetra Anticyclone (Menna et al. 2022), which is a semipermanent ocean eddy in the southeastern Cretan Sea.This agrees with the findings of Mohamed and Skliris (2022) and Bonaduce et al. (2016).The Ierapetra Anticyclone that is usually formed during the summer months (Ioannou et al., 2017), also revealed the highest standard deviation values (>0.1 m; Figure 4f).The characteristics of the Ierapetra Anticyclone (e.g.size, location, intensity, lifetime) are determined by the negative wind stress curl and the blocking of the strong Etesian summer winds by the Cretan orography (Horton et al. 1994;Mkhinini et al. 2014).Generally, the interannual standard deviation showed different behavior between the several sub-regions (Figure 3d), but without revealing any significant trend (p values >0.01;Table 2).
The highest 99 th percentiles in the Aegean Sea were derived for the same regions over the SW Aegean (>0.22 m; Figure 4b), where the strongest increasing trends of the maximum values were also computed (>3.4 mm/year; Table 2; Figure 3b).Although all trends were statistically significant, the smallest p values (<0.0001;Table 2) were found in the SW Aegean and Cretan Seas.High 99 th percentiles of SLA were also observed over the Cretan Sea, along the central Anatolia coasts, and around the eastern Aegean islands (e.g.Lesvos and Chios; Figure 4b).The northern Aegean sub-regions revealed the highest 99 th percentiles during the "extreme" year of 2010 (>0.25 m; Figure 3b), when the maximum SLA values of the entire study period occurred (Figure 4e).Two additional peaks of annual 99 th percentiles were observed in 2019 and 2021 (Figure 3b); large areas of the Ionian Sea and central and southern Aegean sub-regions revealed their highest maximum values during 2021 (Figure 4e).Another area with very high 99 th percentiles (Figure 4b), strong positive trends (Figure 4c), and large standard deviation values (Figure 4f) is located south of the Peloponnesus, where the Pelops Gyre is formed (Bonaduce et al. 2016;Menna et al. 2022;Pinardi et al. 2006).On the contrary, the Ionian Sea showed both weaker interannual variations (StDev <7 cm; Figure 4f) and trends (<3 mm/year; Table 2) compared to the Aegean and Cretan seas.Besides the mean and maximum values, very strong increasing trends were also computed for the minimum SLA values (10 th percentiles) for all regions (Figure 3c) with steeper Sen's Slopes in the Aegean Sea (>3.5 mm/year; Table 2) and milder ones in the Ionian and Cretan Sea (<3 mm/year; Table 2).This result is associated to the increase of the lowest annual 10 th percentiles that are even positive at several years during the last decade (e.g.2013, 2016, and 2018; Figure 3c).
The heat maps and the respective dendrograms of the mean and the 99 th percentile of SLA for each year (row) and each sub-region (columns) are presented in Figure 3e and f, respectively.The 29 annual means are grouped in four categories (Figure 3e).The highest-SLA group (most red in Figure 3e) includes the years of 2010, 2013, 2014, 2016, 2018, 2019, 2020, and 2021, all during the last 12-year period, confirming the strong inter-decadal increase.The other three categories can be grouped into one cluster, including all the years before 2010 (most yellow in Figure 3e).As far as the highest SLA values (99 th percentiles) are concerned, four main categories are also formed.A similar pattern of lower values is observed in two of them for the years before 2010 (most yellow in Figure 3f).The highest 99 th percentiles of SLA, similarly to the annual averages, were also detected in 2009, 2010, 2012, 2013, 2014, 2015, 2016, 2018, 2019, 2020, and 2021 (Figure 3f).The areas can also be grouped in two categories with similar patterns of mean annual SLA for the Ionian sub- regions (1 st cluster) and Aegean-Cretan sub-regions (2 nd cluster).

Seasonal variability
The seasonal interannual variability for all AICS subregions is characterized by similar increasing trends but reveal remarkable differences among the seasons (Figure 5).The highest SLA levels were detected primarily in autumn (Figure 5d) and secondarily in winter (Figure 5a), while the lowest seasonal means occurred in spring (Figure 5b) in agreement with previous findings, derived from tide-gauges (Tsimplis and Shaw 2010) that showed maximum values in most Mediterranean stations in winter and autumn and smaller levels in spring and summer.The main contribution to the sea level cycle in the Mediterranean Basin is the steric component (Gomis et al. 2008).The maximum seasonal values are represented by the seasonal 99 th percentiles that confirm the smaller variance of the autumn values around high mean levels (Mean = 14 cm, StDev = 4 cm and narrow distribution between the 1 st and 3 rd quartiles; Figure 5f).The autumn mean SLA from all years showed unified variability characteristics (1 st group; Figure 5e), while a second group of seasons with similar variability includes winter and summer.Spring also constitutes a separate group that has more similarities with winter and summer.The peak of 2010 is apparent in the winter's timeseries with higher values in the northern Aegean (Figure 5a).
Despite its high mean and maximum values, the interannual Sen's Slopes of autumn were the weakest among all seasons, especially for the Aegean Sea (Table 3).The strongest trends at all regions were detected in the summer, revealing the highest Sen's Slope among all areas and seasons for the SW Aegean (4.7 mm/year) in agreement with the strong annual trends for the same area (Table 2).The Ionian Sea revealed the weakest trends for winter and spring (2.2-2.7 mm/year; Table 3), while both of its northern and southern parts followed the stronger general trend of the AICS domain in the summer and autumn (2.7-3.3 mm/year) due to the high occurrence frequency of the atmospheric low pressure systems that usually form over the central Mediterranean (between the northern African coast and the Ionian; Cavicchia, von Storch, and Gualdi 2014;Makris et al. 2023) in autumn, contributing on the storm surge formation along the western coasts of Greece (Androulidakis et al. 2015(Androulidakis et al. , 2023)).The spring SLAs were the lowest among all seasons (Figure 5b) showing also very small maximum values (99 th percentiles; Figure 5f).Until 2010, the mean spring SLAs were mainly negative, while the highest spring mean levels occurred in 2018 (7 cm; Figure 5b).

Atmospheric contribution
The barotropic effects of wind and atmospheric pressure on the SLA variability over the AICS region are separately analyzed and presented in Figure 6.The mean annual wind speed, SLP and SLA were computed and the Pearson correlation coefficient together with the MK statistical test were computed for all domain's cells covering the entire 29-year period.The highest SLP and SLA correlation coefficients, derived from annual means of the 1993-2021 period, were detected for the Aegean Sea (Figure 6a); the Pearson coefficients were higher than 0.5 and statistically significant (p value <0.01) over the central (Cyclades Archipelago) and eastern Aegean (along the Minor Asia coasts) Sea.On the contrary, the atmospheric pressure has weaker impact on SLA variability over the Ionian Sea (R P <0.3), where the highest correlation coefficients between wind speed and SLA were derived, especially around the central Ionian coasts (R P >0.5; Figure 6b).The wind-SLA coefficients were smaller in the Aegean Sea in agreement with previous studies (i.e.Androulidakis et al. 2015;Krestenitis et al. 2011;Marcos and Tsimplis 2007), who showed that high surges are more associated to the inverse barometer effect and favorable winds in the Aegean and Adriatic (including Ionian) Seas, respectively.The northern Aegean coasts were characterized by lower correlation coefficients between SLA Table 2. Sen's slopes (mm/year) and the respective MK statistical significance thresholds (p value ) for the mean, the 99 th percentiles (99 Perc), the 10 th percentiles (10 Perc), and the standard deviation (StDev), derived from annual mean SLA values, averaged over all regions (Figure 1) for the 1993-2021 period.and SLP in comparison with the rest of the Aegean Sea (Figure 6a), while the wind effect was relatively stronger (Figure 6b).The effect of the wind direction (wind component) is examined in Figures 6c (westerlies), 6d (easterlies), 6e (southerlies) and 6f (northerlies).The positive and statistically significant correlation coefficients (>0.5) between the eastward wind component (positive) and SLA (Figure 6c) represent the increase of SLA under westerlies which are mainly detected over areas of the S Ionian, in SE Aegean, along Minor Asia coasts, and in NE Aegean coastal areas.The negative coefficients between the westward wind component (negative sign) and the high SLA (positive sign) highlight the effect of easterly winds on the SLA increases (Figure 6d); the northern Cretan Sea, the SW and SE Aegean and the broader Thermaikos Gulf (NW Aegean) revealed high negative coefficients (<-0.5).The increased SLA levels of the N Ionian, areas of northern Aegean and SE Aegean were mainly affected by southerly winds (positive correlation coefficients: >0.5; Figure 6e).On the contrary, the northerly winds (negative wind meridional component; Figure 6f) had an opposite effect on the coastal areas of Thermaikos, northern Aegean, and especially along Minor Asia (and eastern Aegean islands) enhancing the reduction of SLA due to the offshore removal of waters (positive correlations: >0.5).The strong Etesian northerly winds (i.e.Meltemia; Poupkou et al. 2011) that prevail over the Aegean Sea during summer months contribute to the coastal upwelling over the eastern Aegean, reducing the sea surface level and temperature (Androulidakis et al. 2017;Mamoutos et al. 2017) and controlling the biochemical processes of the broader area (Chaniotaki et al. 2021).

Coastal inundation at selected areas along the AICS coastline
The numerical simulations of coastal flooding were performed with the CoastFLOOD model (Section 2.4 and Appendix; Makris et al. 2023) at the 20 study areas (Figure 1b), covering the long-term period from January 1 1993 to December 31 2021 (29 years).The simulated results were used to provide daily estimations of Flooded Area (FA) based on the SLA derived along the respective shoreline (boundary condition) of each coastal area (realistic simulations).Moreover, apart from the realistic conditions, we have simulated the inundation at all study areas during several reference values of extreme cases for sea level elevation (idealized simulations).These are characteristic threshold values of storm-induced SLA on the Greek coastal zone (Galiatsatou et al. 2019(Galiatsatou et al. , 2021(Galiatsatou et al. , 2023)), including e.g. the combined influence of storm surges, astronomical tides, mean SLR, and wave setup or conditionally mild estimations of the TWL.To this end, for intercomparison between different scenarios, the idealized numerical experiments involved values of: (i) SLA = 0.5 m, equivalent to a rough upper limit of the 98% storm surge quantile level corresponding to a 50-years return period within a stationary Extreme Value Analysis (EVA) context in a yearly time interval under climate change conditions (SRES scenarios) during 1951-2100 (Makris et al. 2016); this also coincides with the upper threshold of timedependent (nonstationary EVA) estimates of the most likely event of nearshore storm surge in the Aegean Sea, based on fitted parametric trends of bivariate data analytics (Galiatsatou et al. 2019) or simpler EVA approaches along the entire Mediterranean coastline under more recent (RCP) climate change scenarios (Galiatsatou et al. 2023).(ii) SLA = 1.0 m, equivalent to a rough-rounded and common (for Greece) upper limit of the annual return levels of extreme storm surges (0.65-0.75 m), based on the upper 95% confidence bounds under climate change conditions during 1951-2100 (Makris et al. 2016), combined with typical values of the (geographically) mean highest (astronomical) tidal range (≤0.35 m; HNHS, 2015).(iii) SLA = 1.5 m, corresponding to the above estimate (ii) enhanced by adding a rough upper limit estimate of the wave set-up (0.5 m) based on time-dependent (nonstationary EVA) estimates of the most likely event of extreme nearshore waves incident to the Greek coasts (Galiatsatou et al. 2019).(iv) SLA = 2.0 m, corresponding to a conservative lower threshold (to avoid overexaggerated estimations) of the time-dependent estimates of the 100-years return level of extreme TWL (including tides, mean SLR, and wave run-up) for typical Greek mild cross-shore profiles (Galiatsatou et al. 2023).
Please note that the above cases (iii) and (iv) refer to enhanced SLA values by the wave-induced components on the coast (wave set-up and wave run-up, respectively), which are outside the scope of this paper and are presented only for the completeness of coastal inundation simulations.They are not included as upper reference in the coastal inundation hazard assessment calculations below to avoid underestimations of the storm surge-induced hazard.
To allow the comparison between coastal areas with different sizes and to produce a more unbiased  index for the estimation of the coastal flooding for each area, the daily Flood Coverage Percentages (FCP, %) was computed as: where FA i,j is the flooded area for each study region j and simulation day, i, based on the respective satellite-derived actual SLA i,j ; FA ext, j is a considered maximum of the theoretically exposed area to flooding derived for each study region, j, yet based on a common extreme case scenario for the driver of coastal inundation (i.e.SLA).Herein, the FA i,j is derived for all days of the 1993-2021 period (thus producing ixj = 10,632×20 values for FCP i,j ), while the FA ext, j is derived for each area from the CoastFLOOD simulations based on the inundation scenario forced with SLA = 1 m (thus producing j = 20 values for FA ext, j ).
The upper threshold FA ext , used as a reference, corresponds to a diverse (different for each area) hazard (response) magnitude, i.e. flooded area FA i,j , but based on a common value (characteristic to the Greek littoral zone) for the environmental driver of coastal inundation, i.e.SLA ext , in all areas.This was done to secure a common approach for the calculation of unbiased hazard ranking in the 20 areas of study.This way we can compare study areas against each other in terms of flood hazard (please note not risk; this would require a proper quantification of exposure and vulnerability), albeit their sizes differ a lot further adding to quite diverse percentages of lowland parts in the 20 study cases.

Interannual and spatial variability of coastal inundation
The FCP levels are related to two main factors: the sea level high peaks (e.g.99 th percentile of SLA) and the topographic characteristics of each coastal area (e.g.distribution of low-lying areas).Figure 7 presents the distribution of both factors for each study area, dividing the land elevation percentages (lower than 1 m) at 6 categories with a 20 cm step, derived from the highresolution DEM dataset (5 m; Section 2.2).The maximum SLA spatial variation (red line in Figure 7) is related to the majority of the study areas' FCP values (similar variation between FCP and SLA); however large differences (e.g. Area 3) occurred due to the different land terrain of the 20 cases.The FCP levels for all study areas, averaged over the 29-year period (black line in Figure 7), show strong spatial variability with very low percentages (<4%) at Areas 3 (Alexandroupolis), 8 (Agiokampos), and 15-20 (central and northern Ionian Sea).
Respectively, the 99 th percentile of SLA also varies between 0.18 and 0.225 m among the 20 study regions confirming the clear spatial variability of maximum sea levels.It is shown that although relatively high SLA values were recorded for the central and northern Ionian areas (Areas 15-20; >0.21 m 99 th percentile), the low FCP levels were mainly associated with the very small percentages of low elevation cells of the respective areas.The majority of the land cells (>40%) are between 0.8-1.0m, while very high percentages were also computed for categories between 0.4-0.8m, confining the flooding extension at these coastal areas.
Similarly, Area 3 (Alexandroupolis) was also characterized by high SLA levels but with quite few lowlying land cells (only <10% of all cells was less than 20 cm and more than 45% was over 0.8 m), reducing the inundation impact over the inland coastal zone (FCP <4%).Area 4 (Chios) also revealed very small FCPs (<5%) due to the large area characterized by high land elevations (<30% is lower than 0.4 m).Although, in Area 6 (Rethymno), the maximum SLA ranged at similarly high levels as Areas 3 and 4 (99 th Percentile SLA >0.22 m), the FCP was three time larger (>10%) due to the more extended lowing-area; the area with elevation over than 0.8 m was less than 23%, while almost 40% was lower than 0.4 m.The very high FCPs computed for the southern Ionian areas (Areas 11 and 12) were related to both relatively high SLA levels (>0.20 m) and the extended low-lying land (e.g.>40% of Area 11 lay even below the mean sea level: <0 m), allowing extended flooding, especially during strong storm surge events (>16% mean FCP).
The inundation extends under four hypothetical sea level elevation scenarios (0.5-2.0 m) and the maximum realistic record, derived during the study period over the Thermaikos Gulf's coastal zone (Area 1), are presented by overlaid discretely colored areas in Figure 8.It is shown that under the maximum realistic sea level conditions for Area 1, derived from the satellite SLA observations (SLA = 0.292 m; Figure 8a), the flooded waters covered large areas located at the western coast of the gulf, especially over the lowland agricultural areas, and around the two, environmentally protected, main river deltas (Androulidakis et al. 2021).Note that Area 1 is characterized by a large percentage of land elevation below zero (35%; Figure 7), which is mainly located along its western coast (Figure 8).Under higher sea levels (e.g.SLA = 0.5 m; Figure 8a), the simulated flooded area is expanded along the entire western coast (overall FA = 56 km 2 ).The entire western coastal zone up to 4 km inland is even more exposed to extended flooding under extreme conditions of SLA between 1-2 m (FA = 120-180 km 2 ; Figure 8b).Under these extreme conditions, potential inundation may also occur over the urban area of the Thessaloniki port in the north and around the airport of Thessaloniki (Macedonia Airport) at the eastern coast.The daily evolution of the FCP for Area 1 (days with zero FA excluded) during the entire study period is presented in Figure 9a, showing a clear positive trend from approximately 15% in 1993 to 17% at the end of the 29-year period.A peak close to 25%, associated to the SLA = 0.292 m (Figure 8a), is observed at the end of 2009 (Figure 9a).The same year is also the datum of the upward shift (Figure 9b), dividing the 29-period in two separate periods based on the Pettitt's homogeneity test.The FCP averaged over two periods, before (μ1) and after (μ2) the datum of the upward shift was also computed and shown in Figure 9.The difference between the two FCP means (before and after 2009; μ2-μ1) for Area 1 is around 1.1% for the mean annual levels (Figure 9b).The Sen's Slope is almost double for the maximum annual FCP levels (Sen's Slope = 0.11%/year; Figure 9c) showing that the increasing trend is stronger for the more severe flooding conditions than the mean annual levels, confirming the increase of the maximum inundation extension throughout the years, especially after 2009.The μ2-μ1 for the two periods (before and after 2009) is also double (2.4%) for the maximum values compared to the mean levels.The results of Area 1 agree with the general Sen's Slope (Table 4), derived from all 20 study areas' slopes (0.07%/ year for the mean levels and 0.15%/year for the maximum levels).In all study cases, the trends are positive and statistically significant (p value <1%) with stronger slopes for Area 6 (Rethymno) considering both mean (0.25%/year) and maximum (0.40%/year) annual FCP levels.These follow the respective SLA trends, which showed that the Cretan Sea revealed the strongest interannual trend of the SLA 99 th percentiles (3.5 mm/year; Table 2).In all cases, the upward shift took place after 2009, ranging around 1.39% for the mean and 2.86% for the maximum levels (Table 4).The empirical significance  level (p value ; significance level of 1%) was performed for all study areas confirming that the hypothesis of a shift in the mean in the dataset is valid.As a matter of fact, inundation maps for all study areas, similar to Area 1 (Figure 8) are provided as supplementary material (Figures S2-S21).Before 2009, the daily FCP values were usually lower than the thresholds derived from the summation of the median with once (median+σ) or twice (median +2σ) the standard deviation (σ) representing two kinds of extreme levels for the overall FCP values during the 1993-2021 period (Figure 9a).These two FCP thresholds were used to detect occurrence frequency of the severe flooding events during all years of the study period.In particular, the occurrence frequency of the FCP over the median+σ for all study areas is presented in Figure 10.The FCP exceeds the extreme threshold at several dates during the most recent 13 years (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021), when the highest SLA occurred.The interannual variability of both mean SLA and percentages of FCP for each case ranges at different levels between the 20 study areas, confirming the strong spatial variability of both sea level and flooding levels shown in Figure 7.In most of the areas (except Area 7), the interannual trend of the frequency is positive and statistically significant showing very high frequencies during several years after 2009, in agreement with the trends of the mean annual SLA; the variability between the two variables (SLA and FCP) is very similar in all cases.As shown in the regional SLA fields (Section 3), the stronger SLA trends were detected for the southern regions of the AICS (Figure 3) and especially in the Cretan Sea (Table 2).Thus, high coastal Sen's Slopes for Area 5 (Ierapetra; 4.1 cm/decade) and Area 6 (Rethymno; 4.1 cm/decade) resulted to stronger and more frequent flooding events (>median+σ) throughout the years (10%/decade; Figure 10).Very strong trends of severe flooding were also computed for several areas of the Ionian Sea (e.g.Areas 15 and 16, ~10%/decade) that are more exposed to storm surge events associated to the impacts of deep atmospheric depression systems (Makris et al. 2023), e.g.Medicanes (Nastos, Papadimou, and Matsangouras 2018) that form over the central Mediterranean and propagate toward the Ionian Sea (Androulidakis et al. 2023;Lagouvardos et al. 2022).Although the mean FCP levels were relatively low for these two areas mainly due to their topographic characteristics (Figure 7), the extreme flooding events (>median+σ) were more frequent in the last decade resulting to strong interannual trends from 1993 until 2021.Very strong flooding trends were also computed along the northern Aegean coastal zone, particularly in Thermaikos Gulf (Area 1; 9%/decade) and especially along the north-western coastal zone (Areas 8 and 9; >11%/decade), which are exposed to major coastal erosion processes (Kombiadou and Krestenitis 2012;Tsoukala et al. 2015) that may downgrade the natural protection of the coastal zone against sea-originated inundation.

Seasonal variability of coastal inundation
The FCP daily variability presented in Figure 9a for Area 1 is characterized by a clear seasonal signal with higher values in autumn (September-November) and winter (December-February) and lower levels during spring (March-May) and summer (June-August).The seasonal distribution of FCP for all study areas is presented in Figure 11a.In all cases, the higher FCP percentages were computed for autumn, reaching very high levels (>20%) for Areas 11 (Katakolo) and 12 (Laganas) located in the south-central Ionian Sea.The higher autumn FCP values are related to the occurrence of storm surges, related to low-pressure systems (e.g.Medicanes) which are mainly formed in the autumn and early winter months (September to January; Nastos, Papadimou, and Matsangouras 2018).The autumn mean SLA from all years showed unified variability characteristics (1 st group; Figure 5e) determining the high FCP autumn levels.There are a few cases during the spring months, and even fewer during summer, which implies that the "warm" season of the year is not a period of Medicane-induced storm surges (Zhang et al. 2021).In agreement with the Medicanes' seasonality, high FCP levels were also detected in winter with increased values in Katakolo and Laganas (~17%; Figure 11a).Higher FCP values, comparatively to winter levels, were also computed in the summer period for the northern Aegean areas (e.g.Thermaikos, Nestos, Alexandroupolis).The mean seasonal cycle of sea level peaks between October and November everywhere in the Mediterranean Sea, except in the Aegean Sea where it peaks in August.
The interannual variability of the seasonal FCP values, spatially averaged over the selected 20 study areas, confirms the higher flooding level during autumn, which is characterized by the highest Sen's Slope (1.3%/decade) among all seasons.Although the lowest FCP levels were detected for spring in relation with the lowest spring SLA values (Figure 5f), the respective trend slope was also high (1.3%/decade).Winter and summer trend slopes were also positive, but weaker, ranging around 1%/decade and 1.1%/decade, respectively.Almost zero spring values were only detected during the first decade (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003) indicating the total absence of flooding events (e.g. 1993, 1995, 1997, and 2003).This result agrees with the low mean spring SLA levels observed during the same years (between −0.05 to −0.1 m; Figure 5b).The SLA increase during spring after 2003 also enhanced the spring flooding conditions in agreement with Cavicchia et al. (2014) who showed that beside autumn and winter a low number of Medicanes may also form in the spring months.The interannual variability of the seasonal SLP 10 th percentiles over the entire AICS domain shows that the lowest atmospheric pressure levels for all seasons occurred after 2003 (Figure 12), and especially after 2010 for spring (Figure 12b).These low SLPs are associated with the formation of low-pressure systems that control the storm surge variability over the region.12c).This increase is probably related to frequency increases (or intensification) of more ordinary low-pressure systems than to Medicanes which are mainly formed during autumn and winter.Androulidakis et al. (2015), based on climatic predictions, showed that there is a possible 10% increase of the future (21 st century) storm surges during the summer months.Although cyclogenesis in the Mediterranean is more intense from November to March, there are types of cyclones, like Sahara cyclones, whose frequency is larger in summer (Lionello et al. 2006).The FCP peak of 2010 (Figure 11b) occurred during autumn and winter (>12%), when the highest SLA levels occurred (>0.15 m; Figure 5).

Coastal inundation hazard index
To identify a generic categorization of coastal impacts from the flood inundation modeling output, a casespecific and heuristic Coastal Inundation Hazard Index (CIHI) was defined as the product of the Frequency of Occurrence (FO, %) of the days with extreme FCP i,j multiplied with the temporal mean of the FCP i,ext (from the daily values, i, of the 1993-2021 period), and calculated for each study area (j) as shown in Equation 3. The average (FCP i,ext ) j is expressed by the mean of the FCP i,j values above the sum of the median and standard deviation (>median+σ, %) derived from the entire study period for each study area.
The CIHI values were normalized over their min -max range for the 20 studied areas in Greece, to produce individual hazard index results for each of them directly comparable to each other (i.e.ranked values are case-specific), by following a Likert scale classification ranging from very low to very high corresponding to a 1-5 ranking system.The lower and higher CIHI values were chosen as the lower and upper limit of the Likert scale, equally dividing the intermediate ranks of the scale.The ranking of the 20 study areas along the AICS coastal zone, based on the CIHI index, is presented in Table 5.Although this is a case-specified index derived from the 20 study areas, it can be considered representative for the entire Greek coastal zone, since the flooding conditions, the topography characteristics and the sea level of the 20 study regions are characteristic of the entire Greek coastline.The lowest CIHI (1) was derived for several study areas that both FO and mean FCP i,ext were low (e.g.FO = 0.6% and mean FCP = 9.2% for Pineios).Five areas (Areas 1,5,6,11,12) with CIHI larger than 3 were detected, characterized by a higher hazard of coastal inundation, as derived for the 1993-2021 period.Thermaikos Gulf (Area 1) and especially its western coastal area (Figure 8) revealed high FO (10.4%) of extended flooding (FCP i,ext ), while the mean FCP i,ext was around 19% resulting to CIHI = 4.More hazardous flooding conditions (CIHI = 5) were detected for the two areas in Cretan Sea (Ierapetra and Rethymno) showing FO > 14% and mean FCP i,ext >16%.The highest 99 th percentiles of SLA (Figure 4b) and the strongest interannual trend of the 99 th percentiles (3.5 cm/ decade; Table 2) among all AICS sub-regions were computed for the Cretan Sea.Similarly, high CIHI values were computed for two areas of the South Ionian (Katakolo and Laganas), mainly associated with the extended flooding area (mean FCP i,ext >26%) and less with the temporal occurrence frequency of extreme floods (FO was around 10%); the larger FCP values among all study areas were computed for Laganas (mean FCP i,ext = 32%).Very low CIHI levels were derived for almost half of the study areas ( 9), located at all sub-regions of the AICS domain.Nestos and Alexandroupolis (Areas 2 and 3) in the northern coast of the Aegean were ranked with CIHI = 1 mainly due to the low occurrence frequency of extreme floods (2.8%); very weak interannual trends were also computed for these two areas (1.2%/decade and 3.9%/ decade; Figure 10).Low (CIHI = 2) and very low (CIHI = 1) hazardousness was also estimated for several areas of the central northern Ionian Sea (Manolada, Patra, Argostoli, Livadi, Vassiliki, Palairos, Preveza, Igoumenitsa) although the occurrence frequency for most of them was relatively high (>10%); the low values of CIHI in these cases are mainly related to the modeled flood extents, i.e. smaller inundated areas (<10%).Although low CIHI values indicate a weak hazard level of the overall selected ensemble of coastal areas, flooding impacts at specific low-lying parts of the study's coastal zones can be quite high due to a large occurrence frequency of extreme events.Although very low overall CIHIs were derived for Argostoli and Livadi (CIHI = 1), the flooding hazard for the low-lying areas of these cases remains relatively high due to the high frequency occurrence levels (>16%).Moreover, there are areas with very high FCPs (e.g.Nestos) that resulted to low CIHI indices due to the very low frequency of occurrence (<3%).
Although it is not very common, this type of areas may experience extended flooding as derived by the computed FCP i,j .Therefore, although the CIHI index may provide very useful information about the hazard threat of a coastal area, it is necessary to evaluate separately the levels of both flooding coverage (spatial component) and frequency (temporal component).

Discussion
The variability of the SLA over the AICS has been evaluated with the use of satellite-derived observations during a 29-year period (1993 to 2021).The SLAs along the coastal zone of the AICS domain were used to force numerical simulations of coastal inundation at 20 selected low-lying areas to estimate their hazard levels due to increase in sea level elevation.Targeted comparisons of the satellite-derived SLA estimations against coastal tide-gauges showed that especially the high sea level peaks were well captured by the satellites, efficiently describing the extreme conditions, responsible for the formation of littoral inundation events over the coastal zone of the AICS domain.Nevertheless, the gridded coastal values of SLA over the entire Mediterranean basin and its nearshore marine areas, derived from interpolations of satellite altimetry over certain tracks, are known to have contained discrepancies and underestimations of actual sea levels in the past induced by the smoothing effects of gridding.However, several researchers have recently pointed out the ongoing quality enhancement of such products (e.g.Cipollini et al. 2017).
It is furthermore noted that the sea level dataset does not include the influence of the high-frequency sea-surface dynamics controlled by wave-induced Table 5. Means and ranking (1-5) of frequency of occurrence (FO, %), mean Flood coverage percentage (mean-FCP i,ext ) j (%) and coastal inundation hazard index (CIHI*) derived for the 20 coastal areas and the 1993-2021 period.The FO and FCP ext levels were derived from the daily values with FCP>median+σ.Ranking: 1 (very low), 2 (low), 3 (medium), 4 (high), 5 (very high).CIHI over 3 (high and very high) are highlighted with red color.processes that also contribute to coastal flooding (e.g. by wave runup and overtopping).The SLA observations contain though the most important components that control the sea level variability related to longand mid-term coastal inundation, such as the sea level rise (due to climatic changes, e.g.global warming effect) and especially the contribution of storminess.The latter controls the formation of storm surges and the respective coastal flooding during severe lowpressure meteorological systems.Storm surge is the main controlling factor of inundation variability over the coastal zone, especially for floods that extend further inland than the proximal area to the coastline.Waves can contribute to the total water level by means of wave setup through the radiation stress in the surf zone (Wolf 2009), and mainly affect the coastal area by onshore wave runup and overtopping of protection structures.These can cause damages on the built and natural environment along the coast, but mainly in a zone rather close to the shoreline, at least for the Greek coastal zone that is not affected by high and long swells.
The numerical simulations of coastal inundation were performed with CoastFLOOD model (Makris et al. 2023) at 20 study areas of the AICS domain, covering the same 29-year period.The simulated results were used to provide daily estimations of flooded area based on the satellite-derived SLA along the respective shoreline (boundary condition) of its coastal area.CoastFLOOD is a reduced complexity hydraulic flood model, based on the Manning-type mass balance approach for coastal inundation by a wet/dry storage cell concept on fine-resolution raster grids.Therefore, it is susceptible to inconsistencies of the implemented DEM that may induce inaccuracies to flooded area derivation (Karamouz and Fereshtehpour 2019;West, Horswell, and Quinn 2018).We consider that the use of a quite high resolution raster grid (dx = 5 m) may include the most prominent topographic features in the coastal zone, such as the man-made infrastructure, urban settings, ports, coastal structures, roads, and the natural land formations in the floodplain, such as farmlands, beach slopes, dunes, barriers, geodetic anomalies, and other rural land elevation peculiarities (Kahl et al. 2022;Murdukhayeva et al. 2013).Their terrain roughness is also considered, thus including detailed representations of flood flow bottom friction, by taking advantage of freely available official EU-scale land cover data (see Appendix).Nonetheless, CoastFLOOD does not account for the effects of porous bed percolation or floodwater mass flow in sewage and drainage systems, land level conduits, bridge culverts, wells, shafts, etc. in the urban environment.The coastal inundation levels are related to three main factors: a) the sea level high peaks, and b) the topographic characteristics of each coastal area, and c) the flood-prone land's terrain roughness.
We have categorized the study areas based on the computation of a heuristic Coastal Inundation Hazard Index (CIHI) that is determined by the occurrence frequency of extreme floods (temporal component) and their coverage area (spatial component).The CIHI ranking of the coastal areas range from 1 to 5 showcasing the hazard level (very low to very high) for each area, and it is case-specific (i.e.implementing study-related min-max thresholds).High flooding levels occurred at coastal regions of either extended low-lying areas and/ or high sea level rise.Although, high storm surges may occur in the Ionian Sea, the topography characteristics of several coastal areas in the central and northern parts (e.g.Argostoli, Livadi, Vassiliki, Palairos, Preveza, Igoumenitsa) controlled the flooding extents, resulting to relatively low CIHI (<3).Two areas of the southern Ionian (Katakolo and Laganas) were characterized by extended inundation incidents during the 1993-2021 period due to both high sea levels and extended low-lying areas in the coastal zone (CIHI = 5).The Sen's Slopes of the inundation levels computed for all 20 study areas were positive (increasing trends) and statistically significant (0.07%/year for the mean levels and 0.15%/year for the maximum levels).In all cases, a clear upward shift took place after 2009, ranging around 1.39% for the mean and 2.86% for the maximum levels.The highest SLA for all areas occurred in the years of 2010, 2013, 2014, 2016, 2018, 2019, 2020, and 2021, all during the last 12-year period, confirming the strong inter-decadal increase that enhanced the flooding levels after 2009.Moreover, in most of the areas, the interannual trend of the frequency of extreme floods (>median+standard deviation) is positive and statistically significant, revealing high frequencies during several yearly cycles after 2009 in agreement with the trend of the mean annual SLA.Very strong trends of extreme flooding were computed for the Cretan Sea, where the strongest trends of maximum SLA occurred.More hazardous flooding conditions (CIHI = 5) were detected for the two areas of the Cretan Sea (Ierapetra and Rethymno) showing both high frequency (>14%) and large inundation areas (>16%), agreeing with findings in relevant literature (Martzikos et al. 2018;Makropoulos et al., 2015).Strong trends of extreme flooding were also derived for several areas of the Ionian Sea (e.g.Argostoli and Livadi, ~10%/decade) that are more exposed to storm surge events associated to the impacts of severe lowpressure systems (e.g.Medicanes; Nastos, Papadimou, and Matsangouras 2018).Romera et al. (2017) investigated the future formation of Medicanes up to 2100, and estimated a future reduction in their number but an increase in the intensity of the strongest cyclones that may lead to more severe storm surges and respective coastal floods during the 21 st century.The area between Sicily and Greece extending southward to the coast of Libya is the second preferred area where Medicanes usually form (Cavicchia, von Storch, and Gualdi 2014).These extreme meteorological phenomena are triggered by the simultaneous formation of cold air aloft and the presence of very warm surface waters.The increased sea surface temperature (SST) trends during the last decades over the Mediterranean Sea (Androulidakis and Krestenitis 2022;Darmaraki et al. 2019) formed more favorable conditions (warmer water pools), strengthening their intensity and the accompanying storm surges.
To extend the findings of our research, the next step is to develop climatic projections of coastal inundation under different climate-change scenarios (Makris et al. 2023) and to further test if the current increasing trends in the AICS domain are also probable in the 21 st century (Krestenitis et al. 2021).Moreover, for the complete assessment of integrated seawater flood probability at coastal sites of the AICS domain that are at risk due to littoral inundation (Masina, Lamberti, and Archetti 2015), the next step is to include satellite observations of waves.This will allow for a better estimation of the combined effect of waves and water levels on coastal flooding, especially in shoreline proximal zones along low-lying littoral areas.

Conclusions
In the current study, the Sea Level Anomaly (SLA) variability over the Aegean, Ionian, and Cretan Seas (AICS) and the associated coastal inundation along the coastline were investigated based on long-term sea level observations, high-resolution land elevation data, and numerical simulations of coastal flooding.
The general SLA trend over the entire AICS domain was 3.6 mm/year, which is higher than the overall Mediterranean (2.44 mm/year; Bonaduce et al. 2016) and the eastern Mediterranean (3.1 mm/year; Mohamed and Skliris 2022) trends.The highest Sen's Slopes were derived for the Aegean Sea with strongest trends for the SW Aegean (4.1 mm/year).Moreover, the Cretan Sea revealed the strongest interannual trend of the maximum SLA values (99 th percentiles; 3.5 mm/ year).The highest SLA levels were detected primarily in autumn and secondarily in winter, while the lowest seasonal means occurred in spring in agreement with Tsimplis and Shaw (2010).The inundation levels and the respective positive trends were higher in autumn associated with the low-pressure systems (storm surge formation) that mainly develop in the autumn and early winter months (September to January).Although the lowest inundation levels were detected during spring season related to the lowest spring SLA values (caused by the coincidence of neap tides with extended calm weather periods of high barometric systems), the respective trend slope of spring floods was also high.This is in agreement with Cavicchia et al. (2014), who showed that besides autumn and winter a significant number of severe storms may also form during spring.The Ionian Sea revealed weak SLA trends for winter and spring, while both of its northern and southern parts followed the stronger general trend of the AICS domain in the summer and autumn seasons.This is due to the high occurrence frequency of the atmospheric lowpressure systems that usually form over the central Mediterranean (between the northern African coast and the Ionian; Cavicchia, von Storch, and Gualdi 2014) during early-to mid-autumn, contributing to the storm surge formation along the western coasts of Greece (Androulidakis et al. 2015(Androulidakis et al. , 2023)).
The meteorological residual of the sea level, controlled by Sea Level Pressure (SLP) and winds, shows a strong spatial variability over the AICS domain, with a stronger SLP impact over the Aegean Sea and more significant wind effects in the Ionian Sea.Moreover, the sea level variability is also controlled by the direction of the prevailing wind, whose effect varied along the topographically complicated AICS coastal zone.The occurrence of storm surges in the Ionian Sea is strongly related to the wind-induced accumulation of waters along the Ionian coasts, especially during severe low-pressure systems that usually form over the central Mediterranean.Androulidakis et al. (2023) showed that the accompanying winds of an extreme meteorological system in September 2020 (Ianos Medicane) increased the water level even before the landfall of the cyclone's core to the central and southern Ionian coastal zone.Similar conditions have mainly affected the large low-lying areas of the southern Ionian Sea, such as Katakolo and Laganas, characterized by both large inundation extents and high occurrence frequency of floods.On the contrary, the inverse barometer effect was more profound in the Aegean Sea due to the morphological complexity of the archipelago that reduced the wind impact on storm surge formation in comparison to the direct impact by the SLP (Krestenitis et al. 2011;Makris et al. 2023).Strong flooding trends were also computed along the northern Aegean coastal zone, particularly in Thermaikos Gulf (9%/decade) and especially along the north-western coastal zone (>11%/decade).These are exposed to coastal erosion processes (Kombiadou and Krestenitis 2012;Tsoukala et al. 2015) that might have downgraded the natural protection of the coastal zone against seawater inundation.
The 20 study areas represent different parts of the AICS coastal zone that are characterized by different morphological features and sea level conditions.The study areas were categorized based on the Coastal Inundation Hazard Index (CIHI), based on the occurrence frequency and flood coverage of inundation events.The interannual trends and the seasonality of the flooding also showed intense variation between the study areas.Our approach and study findings provide useful information about the hazard level of sea level-induced coastal flooding in the AICS domain that can hopefully support more robust coastal zone management.We set the full assessment of coastal flood risk as a future goal, by further quantifying the exposure, vulnerability, and resilience of the study areas.

Table 1 .
Satellite data skill metrics of reproducing the in situ observations of SLA (m; Dynamic atmospheric correction (DAC) product removed) in Alexandroupoli, Thessaloniki, Chios, Lefkada derived from mean daily and monthly values from 2002 until 2012; quantitative comparisons are based on the statistical measures: Pearson correlation R p , root mean square error (RMSE), Willmott skill score (WSS), HRP-index.The asterisk (*) indicates that the hypothesis that the correlation is statistically significant is true (99% MK test of statistically significant trend: p value <0.01).

Figure 3 .
Figure 3. Annual variability (dashed lines) and trends (continuous lines) of: (a) the mean Sea level Anomaly (SLA; m); (b) the 99 th percentile; (c) the 10 th percentile, and (d) the standard deviation (StDev; m), averaged over 8 regions (all regions: black; N Ionian: red; S Ionian: orange; NW Aegean: dark blue; NE Aegean: light blue; SW Aegean: purple; SE Aegean: light green; Cretan Sea: gray (Figure1a), for the period 1993-2021.Heat maps with respective dendrograms of the (e) mean and (f) 99 th percentile SLA derived from the annual and spatial SLA averages for each year and region.

Figure 4 .
Figure 4. Horizontal distribution of (a) mean Sea level Anomaly (SLA; m), (b) 99 th percentile, (c) Sen's Slope (m/year), (d) Mann-Kendall (MK) trend test (p value ), (e) year of the maximum SLA, and (f) standard deviation (StDev) derived from the daily and mean annual SLA over the AICS domain during the 1993-2021 period.

Figure 5 .
Figure 5. Seasonal variability (dashed lines) and annual trends (continuous lines) of the mean Sea level Anomaly (SLA; m) over (a) winter, (b) spring, (c) summer, and (d) autumn, averaged over 8 regions (all regions: black; N Ionian: red; S Ionian: orange; NW Aegean: dark blue; NE Aegean: light blue; SW Aegean: purple; SE Aegean: light green; Cretan Sea: gray (Figure 1), for the period 1993-2021.(e) heat map with respective dendrograms derived from the annual and spatial SLA seasonal averages for each year and region.(f) seasonal box plots of the 99 th percentiles of SLA of all regions, showing the 1 st (lower line) and 3 rd (upper line) quartiles, minimum (low dot), maximum (upper dot), mean (cross) and median (middle line), and standard deviation (StDev) value for each season.

Figure 6 .
Figure 6.Distribution of Pearson correlation coefficients between Sea level Anomaly (SLA) and (a) Sea level pressure (SLP), (b) wind speed, (c) positive zonal wind component (westerlies), (d) negative zonal wind component (easterlies), (e) positive meridional wind component (southerlies), and (f) negative meridional wind component (northerlies), derived from the mean annual values of the 1993-2021 period.The solid and dashed contours indicate the 99% and 95% p values isolines of the correlation statistical test, respectively.

Figure 8 .
Figure 8. Maps of flood coverage (inundation areas) under the (a) maximum Sea level Anomaly (SLA) derived from the 1993-2021 period (SLA max = 0.292 m) and a hypothetical scenario (SLA = 0.5 m)over area 1 (Thermaikos Gulf).(b) Flood covarages under more extreme scenarios (SLA = 1, 1.5 and 2 m) are also presented.Inundation areas induced by lower SLA are overlaid to the ones driven by higher SLA values.The land terrain showing the area's elevations is presented with shades of gray.Similar maps for all study areas (2-20) are provided as supplementary material (Figures S2-S20).

Figure 9 .
Figure 9. (a) daily evolution of Flood coverage percentage (FCP:0-1) for area 1 (Thermaikos Gulf) during 1993-2021.The green dashed and red lines indicate the median+σ and median + 2σ levels.The linear regression is also marked with blue line.The zero FCP values (no flooding) are excluded from computations.Evolution of annual (b) mean and (c) maximum FCP values.The Sen's Slope (%/year) and the linear regression (black line) are also shown.Pettitt's test for homogeneity was also applied to the two annual timeseries showing the shifts in the mean levels (μ1, μ2).

Figure 10 .
Figure 10.Annual evolution of mean Sea level Anomaly (SLA; blue lines) and annual frequency of flooding indices above median +σ level (red lines) for 20 coastal areas during 1993-2021.The linear regressions and the respective Sen's slopes for each case are also shown.
De la Vara et al. (2021) found the highest Medicane activity in the winter months, with an intermediate activity also in spring and autumn.Although De la Vara et al. (2021) also argued that there is no Medicane activity in the summer months, the interannual variability of the summer minimum SLPs during the 1993-2021 period was characterized by a clear (p value <5%) change after 2003, with generally low-pressure levels during the following period based on the Pettitt's homogeneity test (Figure

Figure 11 .
Figure 11.(a) seasonal Flood coverage percentage (FCP, %) for all study areas during 1993-2021 and (b) evolution of mean FCP averaged over all areas and seasons.The linear regressions and the respective Sen's slopes (%/decade) for each season are shown.All trends are statistically significant based on the MK trend test (p values <0.01).

Figure 12 .
Figure 12.Pettitt's tests for homogeneity for all (a) winters, (b) springs, (c) summers, (d) autumns based on the 10 th percentiles of Sea level pressure (SLP) over the AICS domain (all regions), showing the potential shifts in the mean levels(μ, μ1, μ2).The p value (alpha = 5%) of the homogeneity test (if there is a date at which there is a change in data) for each season is also shown.

Table 3 .
Sen's slopes (mm/year)and the respective MK statistical significance thresholds (p value ) for the mean SLA derived from the interannual SLA values for each season, averaged over all regions (Figure1) for the 1993-2021 period.

Table 4 .
Sen's slopes (%/year) and difference of mu2-mu1%) derived from the Pettitt's homogeneity test of the mean and maximum annual FCP (%/year) for the 20 study areas during the 1993-2021 period.All trends are statistically significant based on the MK trend test (p values <0.01).
Note that the CIHI ranked output is case-specific to the 20-area analysis along the Greek coastal zone and not valid for "global" comparisons application.