Relative Sea-Level Rise Scenario for 2100 along the Coast of South Eastern Sicily (Italy) by InSAR Data, Satellite Images and High-Resolution Topography

Abstract

The global sea-level rise (SLR) projections for the next few decades are the basis for developing flooding maps that depict the expected hazard scenarios. However, the spatially variable land subsidence has generally not been considered in the current projections. In this study, we use geodetic data from global navigation satellite system (GNSS), synthetic aperture radar interferometric measurements (InSAR) and sea-level data from tidal stations to show the combined effects of land subsidence and SLR along the coast between Catania and Marzamemi, in south-eastern Sicily (southern Italy). This is one of the most active tectonic areas of the Mediterranean basin, which drives accelerated SLR, continuous coastal retreat and increasing effects of flooding and storms surges. We focus on six selected areas, which show valuable coastal infrastructures and natural reserves where the expected SLR in the next few years could be a potential cause of significant land flooding and morphological changes of the coastal strip. Through a multidisciplinary study, the multi-temporal flooding scenarios until 2100, have been estimated. Results are based on the spatially variable rates of vertical land movements (VLM), the topographic features of the area provided by airborne Light Detection And Ranging (LiDAR) data and the Intergovernmental Panel on Climate Change (IPCC) projections of SLR in the Representative Concentration Pathways RCP 2.6 and RCP 8.5 emission scenarios. In addition, from the analysis of the time series of optical satellite images, a coastal retreat up to 70 m has been observed at the Ciane river mouth (Siracusa) in the time span 2001–2019. Our results show a diffuse land subsidence locally exceeding 10 ± 2.5 mm/year in some areas, due to compacting artificial landfill, salt marshes and Holocene soft deposits. Given ongoing land subsidence, a high end of RSLR in the RCP 8.5 at 0.52 ± 0.05 m and 1.52 ± 0.13 m is expected for 2050 AD and 2100 AD, respectively, with an exposed area of about 9.7 km² that will be vulnerable to inundation in the next 80 years.

1. Introduction

Sea-level rise (SLR) is one of the major consequences of global warming, driving the melting of ice sheets and the thermal expansion of the oceans, increasing the vulnerability of coastal areas to flooding. Over the past two centuries, the global sea-level (GSL) has risen at faster rates than in the last two millennia [1,2], with values from 1.7 mm/year in the late 20th century up to 3.2 mm/year over the last decades [3,4,5,6]. The latest report of Intergovernmental Panel on Climate Change (IPCC, Special Report on the Ocean and Cryosphere in a Changing Climate SROCC [7]) shows that in 2100, GSL projection has an upper limit of about 1.1 m higher than present sea-level, while about 400 million of people will be highly exposed to coastal hazard. SLR can be exacerbated by land subsidence from natural (i.e., tectonic and volcanic activity and sediment compaction) and anthropogenic causes (i.e., fluid withdrawal, dams building), accelerating the submersion process of low-elevated continental [8,9,10] and insular coasts [11,12,13]. Coastal flooding is likely to be the biggest socioeconomic impact of SLR in the 21st century [14] exceeding the damage caused by earthquakes and volcanic eruptions.

Recently, Toimil et al. [15] delineated the requirements for a scientific approach for future projections of shoreline changes, which should consider changes in the mean sea-level, storm surges, coastal erosion and other additional factors, with the aim of reducing uncertainty in coastal change estimates. Indeed, the rate of relative SLR (RSLR) can vary significantly due to other processes such as isostatic adjustments, ocean currents, seismic and volcanic activity, besides the already cited natural and anthropogenic land subsidence. These phenomena are already affecting the global coasts and large coastal cities [16], including the Mediterranean [17] (www.savemedcoasts.eu and www.savemedcoasts2.eu, accessed on 4 February 2021 [18]). SLR, in conjunction with extreme storm surge events, caused dramatic flooding, such as in 1953 for the coast of Netherlands or more recently in New Jersey (2012, Hurricane Sandy), Louisiana (2015, Hurricane Katrina), Florida (2016, Hurricane Matthew), Louisiana and Texas (2017, Hurricane Harvey) and in the Mediterranean Sea in 2019 and 2020 [19,20,21].

Because relative sea-level change along the coast is given by the sum of eustatic, thermo-steric, isostatic, sediment compaction and tectonic factors [22,23], the projection of the future coastline positions in specific zones require to consider in the analysis: (i) the local rate and trend of sea-level change and (ii) the rate of vertical land movements (VLM), including tectonic and isostatic contributions and anthropogenic influence. The first are provided by the latest IPCC reports [7,24] and can be downscaled for specific regions, while the second are provided by geodetic data, such as sparsely GNSS stations and InSAR observations [25]. The latter can much improve the spatial coverage and resolution to characterize and detail land subsidence, whereas GNSS provide precise measurements in well-defined, ground-based, reference frame, to which InSAR measurements can be referred. Assuming that VLM rates estimated from geodetic data will remain constant for the next few decades in a specific coastal zone, i.e., in absence of any relevant tectonic event, relative sea-level rise projections and flooding scenarios for the coming years (2100) can be reasonably realized. Scenarios can be mapped on high resolution 3D topography of the investigated zones, supporting the analysis of the flooding extension and the assessment of the coastal impacts. Digital Elevation Models (DEM) extracted by LiDAR data are the most suitable tools to realize realistic scenarios.

Here, we report on accurate observations of VLM along the coast of southern Sicily between the city of Catania and the village of Marzamemi, from a combination of GNSS, InSAR and tide gauge observations at an unprecedented resolution and accuracy (Figure 1 and Figure 2). This allowed us to provide the first multitemporal RSLR projections up to 2100 for this area, providing high resolution inundation maps based on the IPCC projections downscaled for the Mediterranean Sea. Projections of revised SLR for the central Mediterranean provided by the IPCC 2019 (SROCC report, Oppenheimer et al. [7,26]) predict a rise of 0.15 ± 0.04 to 0.20 ± 0.05 m by 2050 and 0.33 ± 0.07 to 0.72 ± 0.13 m by 2100, relative to 2016, based on different RCPs of greenhouse-gas emission scenarios.

IPCC projections have already been used to identify areas of potential inundation in coastal regions of the Mediterranean [11,17,18,30,31,32] and especially in USA for planning flood resilience strategies for large coastal urban areas such as in San Francisco, the Los Angeles Basin and New York [16,33,34].

In this study, we focus on six selected areas characterized by low elevated topography, relevant coastal infrastructures and natural areas, where the expected SLR could be a potential cause of large inland flooding and morphological changes of the coastal stretches (Figure 2).

To reconstruct the retreating coasts in the recent past, optical images from satellites over the last two decades have been analysed in conjunction with tide gauge data and current rates of VLM estimated in this study. An integration between RSLR and effective horizontal shoreline migration allowed to map the future submersion surfaces on high-resolution Digital Elevation Models, thus to predict the submersion scenarios up to 2050 and 2100.

2. South-Eastern Sicily: Geodynamic and Geological Setting

South-eastern Sicily, between the Catania Plain in the north and Capo Passero in the south, is shaped by compressional, extensional and strike slip movements related to the relative motion of the Africa and Eurasia tectonic plates, where geodynamic processes associated with subduction of the Ionian lithosphere beneath the Calabrian arc cause a complex volcano-tectonics [35,36] (Figure 1). The geology of this region shows thick Mesozoic to Quaternary carbonate sedimentary sequences and volcanic layers forming the emerged foreland of the Sicilian–Maghrebian thrust belt [37]. The most important tectonic domains of this area are the Hyblean Plateau and the Malta escarpment. The latter is a Mesozoic boundary separating the continental domain from the oceanic crust of the Ionian basin, reactivated during the Quaternary [38,39,40].

Since the Early-Middle Pleistocene, active faulting has contributed to extensional deformation along the coastal sector of south-eastern Sicily, where NNW–SSE trending normal faults control the Ionian shoreline [40,41,42]. These structures are mostly located off-shore and their Quaternary activity is probably associated with the recent reactivation of the Malta Escarpment system [43]. This area is marked by a high level of crustal seismicity that released earthquakes of a magnitude of about 7, such as the destructive events occurred in 1169, 1542 and 1693 [44,45]. The seismogenic sources of these historical events is still debated but they are likely located in the Malta Escarpment, between Catania and Siracusa [36,40,42,43,44,46,47,48,49]. The last major earthquake (ML 5.4) occurred on December 13, 1990 in the Augusta off-shore [50,51,52].

Due to geodynamic processes, such as the retreating of the subducted Ionian slab [53,54,55] or asthenospheric flow beneath a decoupled crust [56,57], and off-shore normal faulting [43,58], during the Late Quaternary, this area has been affected by regional uplifting. Uplift progressively decreases north of southern Calabria and south of north-eastern Sicily, as shown by flights of marine terraces developed along the coasts [53,59,60,61]. In the northern sector of south-eastern Sicily, the long-term uplift has been estimated at rates of 0.2–0.7 mm/year from Middle-Upper Quaternary marine terraces and paleo-shorelines [40,58,62,63], gradually decreasing to zero toward the stable areas of the south-eastern corner of Sicily [59,64].

Toward the south-eastern corner of Sicily, the observed uplift also decreases during the Holocene. In fact, archaeological and borehole evidence show vertical land stability or weak uplift during the late Holocene [65,66,67,68].

On the other hand, levelling surveys performed since 1970 [61,69] and more recent GNSS data [29,70,71] and InSAR analysis [72], highlighted a diffuse and interseismic low rate land subsidence between Catania and Siracusa, which is in contrast with the long term geological uplift of this region.

The studied coastal area has been affected by several marine extreme events in historical times. Effects of several tsunamis have been reconstructed from the analyses of boulder accumulation [73], high-energy deposits [74] and cores performed inside lagoons [75,76]. The analysed area also experienced the effects of several storms and Medicanes occurred over the last decades, mainly represented by boulders dislocation along the coastal area [19]. This is particularly important considering that in a SLR scenario the effects of the extreme marine events will probably impact on the coastal landscapes currently emerged [31,32,68,77,78,79,80].

3. Material and Methods

To realize high resolution maps of expected coastal flooding by 2100, we have analyzed the following data sets: (i) geodetic data from GNSS networks and InSAR observations to estimate the current rates of VLM for the last two decades; (ii) time series of sea-level data collected in the time span 1992–2020 at the tide gauge station of Catania; (iii) optical satellite images to evaluate the intervening coastal retreat occurred in the time span 2001–2019; (iv) sea-level projection released by the IPCC (Report SROCC [7]) corrected for the Mediterranean Sea to obtain RSLR projections for 2050 and 2100 in the RCP 2.6 and RCP 8.5 climatic scenarios. Finally, we analyzed LiDAR data to extract high resolution DEMs to obtain detailed maps of flooding scenarios.

3.1. GNSS Data

We used the available raw observations provided by GNSS-RING network [81], managed by INGV [82], integrated by other active GNSS stations belonging to regional and national networks [83], in order to estimate the VLM rates at a set of stations located within 10 km from the coastlines of south-eastern Sicily (Figure 2). In our analysis, we considered all the data available from 1995 to 2020, but, specifically, the sites shown in Figure 2 span the 2006–2020 time-interval. The analysis was performed by the GAMIT/GLOBK software, following the three step procedures described in Serpelloni et al. [29].

The raw GPS observables have been analysed using the GAMIT/GLOBK software (V.10.70) [84,85] adopting IGS standards [86]. In this step, the satellites orbit parameters are tightly constrained to the IGS final values. The GAMIT software estimates station positions, atmospheric delays, satellite orbits and Earth orientation parameters from ionosphere-free linear combination GPS phase observables using double differencing techniques in order to eliminate phase biases caused by drifts in the satellite and receiver clock oscillators. GPS pseudo-range observables are also used to constrain clock timing offsets and improving automated editing of phase data, helping in the resolution of integer phase ambiguities. GPS phase data are weighted adopting an elevation-angle-dependent error model [85]. The IGS absolute antenna phase centre model is used for both satellite and ground-based antennas, allowing us to improve the accuracy of vertical site positions by mitigating reference frame scale and atmospheric mapping function errors. While the first-order ionospheric delay is eliminated by the ionosphere-free linear combination, the second-order ionospheric corrections are applied using IONEX files from the Centre for Orbit Determination in Europe (CODE). The tropospheric delay is modelled as piecewise linear model and estimated using the Vienna Mapping Function 1 (VMF1). We use the Global Pressure and Temperature 2 (GPT2) model to provide a priori hydrostatic delays. The pole tide was also corrected in GAMIT adopting IERS standards. The Earth Orientation Parameters (EOP) are tightly constrained to priori values obtained from IERS Bulletin B. Ocean loading is corrected using the FES2004 model. The International Earth Rotation Service (IERS) 2003 model for diurnal and semi-diurnal solid Earth tides was set. Because of the large number of stations included in our Euro-Mediterranean GNSS data processing (~3000), the GAMIT analysis is performed for several sub-networks, each made by <50 stations, with each sub-network sharing a set of high-quality IGS stations that are later used as tie-stations in the combination step.

The daily sub-nets, loosely constrained, solutions are combined using the GLOBK software, which adopts a Kalman filter estimation algorithm, simultaneously realizing a global reference frame by applying generalized constraints. Specifically, we define the reference frame by minimizing the velocities of the IGS core stations (http://igscb.jpl.nasa.gov, accessed on 6 January 2021 [87]), while estimating a seven-parameter transformation with respect to the GPS realization of the ITRF2014 frame, i.e., the IGS14 reference frame.

In the third step, we analyse the position time series in order to estimate the 3-components’ (east, north and vertical) of linear velocities and uncertainties. Changes in stations positions are modelled using the following functional model:

$$x\left(t\right)=x_0+bt+\alpha \sin \left(\omega t+\phi \right)+\sum _{j=1}^n∆x_{j}H\left(t-t_j\right)$$

(1)

where x is the position of a point, t is the time, x₀ is the initial position bias, b is the secular rate, α and φ are the amplitude and phase of the annual and semi-annual seasonal signals, respectively, and H is the Heaviside step function defining coordinate jumps (Δx) at a given time t_j. Only stations with a minimum time-span of 4.5 years are retained in this and subsequent analyses, in order to avoid biases due to unreliable estimated seasonal signals and underestimated velocity uncertainties due to absorbed correlated noise content in estimated trends of short time series. We estimate the stations rates uncertainties assuming a white–flicker noise model, using the approach implemented in the CHEETAH software [88].

3.2. InSAR Technique and Data

Since the 1990s, the Synthetic Aperture Radar Interferometry (InSAR) technique allowed the detection and measurements of crustal movements from space with unprecedented accuracy and spatial coverage [89,90]. The InSAR is a non-invasive technique [91], suitable to monitor large areas of the Earth’s surface, measuring the projection of the deformation vector onto the Line of Sight (LoS) direction.

In the last 20 years, the increased availability of new SAR instruments and satellite constellations, has stimulated a steady improvement of processing algorithms. Several multi-temporal InSAR techniques have been proposed, exploiting the redundancy offered by hundreds of image pairs, i.e., hundreds of interferograms, to retrieve mean ground velocity and time-series of relative ground displacements. The existing algorithms fall into two broad categories, based on the Permanent Scatterer (PS) [92] and the Small Baseline (SB) [93] approaches, although more recently algorithms exploiting the basic principles of both methodologies have also been proposed [94].

In this study, we adopted the SBAS methodology applied along both the ascending and descending orbit of the Sentinel-1A satellites (C-Band SAR sensor wavelength=5.6 cm) operated by the European Space Agency (ESA) in the TOPSAR acquisition mode (VV polarization) and free-distributed (https://sentinel.esa.int/web/sentinel/missions/sentinel-1, accessed on 6 January 2021 [95]). Both the ascending and descending SAR datasets have been processed using the Geohazard Exploitation Platform (GEP, https://geohazards-tep.eu/, accessed on 6 January 2021 [96]) and the P-SBAS service, implemented therein. An initial multi-looking operation was applied to the Single Look Complex images equal to 20 and 5 looks along the range and azimuth direction, respectively, resulting in a ground pixel resolution of 90 meters. The main processing steps of the SBAS approach consist of differential interferograms generation from the formed SAR image pairs with a small orbital separation (spatial baseline) to reduce the spatial decorrelation and topographic effects. The Shuttle Radar Topography Mission (SRTM) [97] with 1 arc-second (~30 m pixel size) DEM from NASA and precise orbits from the European Space Agency (ESA) were used for the co-registration and for topographic phase removal from the interferometric phase on each of the compute interferograms. Then, a filtering operation is performed to improve the obtained interferogram’s quality and make easier the following phase unwrapping step according to a fixed coherence threshold (equal to 0.7 for both the ascending and descending processing). Then, an inversion step is performed to put together the ground displacement time series in the time interval covered by the considered SAR acquisitions. At this stage, the atmospheric artifacts are also estimated and removed through a double filtering in the spatial and temporal domain (high and low pass filtering equal to 1200 meters and 365 days, respectively) [93,98]. Finally, the geocoding step is applied by means of the adopted DEM, obtaining a 90m ground pixel in the WGS84 datum to reference the displacements estimated for both the orbits.

In detail, we considered 60 images along the descending orbit (track 124) covering the time interval from 6 July 2017 to 8 March 2020 with a revisiting time of 12 days, enough to detect the movement in the investigated area, and resulting in 163 pairs. The same procedure was considered for the ascending case (track 44) selecting 103 images relative to the 3 October 2016–28 March 2020 temporal span and forming 290 interferograms. The retrieved InSAR ground velocity maps along the satellites LoS are showed in Figure 3 and Figure 4 for the ascending and descending orbits, respectively.

Thanks to the availability of both the ascending and descending InSAR outcomes, we retrieved the surface movement along the vertical and eastward directions. This can only be done for the pixels common to both the ascending and descending LoS maps.

Because of the SAR geometry acquisition, the east–west component of ground motion is proportional to the difference between descending and ascending LoS maps, and it is also a function of the local incidence angles of SAR; the vertical is proportional to the sum of the two datasets, and depends on the local incidence angles, too. Unfortunately, the north–south component of the ground movement has a little effect on the LoS measurements (less than 10% of the north–south movement is captured by SAR); therefore, we cannot estimate it. The analytical expression that we applied to calculate vertical and east-west motions is reported in Dalla Via [99].

Generally, ground movements along the coastal stripes are characterized by subsidence phenomena due to the coastal erosion and soil compaction causing as a consequence the retreat of the coastline [80,100,101]. Moreover, the tectonic regime of the area of interest foresees vertical movements, e.g., rise of marine terraces [58,60,62,65]. Therefore, the horizontal component of the retrieved LoS velocity is negligible in the areas investigated in this study (see also S1—Supplementary Material).

InSAR products were validated and calibrated through a comparison with the GNSS data superimposed and projected on the InSAR LoS directions, for both the orbits (see S1—Additional Material for further information). Such decomposition was made pursuant to the InSAR and GNSS data comparison and validation; hence, the reported Up component already include any post-processing performed operation to calibrate InSAR products (see S1—Additional Material for further details).

3.3. Sea-Level Data

The analysis of sea-level data collected at the tidal networks located in the Mediterranean basin highlighted a mean rate of sea-level rise of about 1.8 mm/year for the last two-three centuries [70], although local VLMs can influence, even severely, the mean trend [102]. The closest tide gauge station in the investigated area of south eastern Sicily is located in the harbor of Catania (ISPRA tidal network, www.mareografico.it, accessed on 1 January 2021 [103]). The analysis of the time series of monthly mean sea-level data collected in the time span 1992–2020, shows a trend at 4.66 ± 0.01 mm/year (Figure 5), which is about 2.9 mm/year higher than the average rate for the Mediterranean Sea. It is worth noting that InSAR analysis detected a local land subsidence at about 2 mm/year in the proximity of the coast of Catania (Figure 3 and Figure 4) that affects sea-level recordings. In addition, RSLR for this location could be exacerbated by additional subsidence due local instabilities of the tidal station location being placed at the top of a large pier in the Catania harbor (Station coordinates Lat 37°29′44.42′′N; Long 15°05′42.11′′E).

In order to project the expected sea-levels at 2050 and 2100 A.D. for the coast of south eastern Sicily, the novel sea-level projections based on the IPCC special report “The Ocean and Cryosphere in a Changing Climate (SROCC)” [7], released in 2019, were adopted. SROCC estimates consist of sea-level median values and standard errors on a worldwide grid obtained by summing up the contributions of geophysical sources driving long-term sea-level changes. The considered sources of sea-level variations are: the thermosteric/dynamic contribution (from the same 21 CMIP5 coupled atmosphere-ocean general circulation models AOGCMs used for AR5 IPCC projections), the surface mass balance and dynamic ice sheet from Greenland and Antarctica, the glacier and land water storage, the glacial isostatic adjustment (GIA) and the inverse barometer effect. As for the AR5 IPCC, SROCC projections are provided for three different RCPs 2.6, 4.5 and 8. 5: in the RCP 2.6 scenario the emissions of greenhouse-gas emission peak around 2020, in RCP 8.5 they continue to grow throughout the century [104]. SROCC and AR5 IPCC sea-level projections only differ in the estimates of the contributions from Antarctica: new ice-sheet modelling results, not available at the time of AR5, have been incorporated into the SROCC assessment. The difference between global mean sea-levels in SROCC and AR5 IPCCS is negligible under RCP 2.6 and RCP 4.5 scenarios while it is about 10 cm in 2100 under RCP 8.5.

The predicted sea-level rise in the south-eastern margin of Sicily has been estimated by combining the rate of sea-level rise at the grid point closest to Lat 37°2′54.32′′N; Lon 15°17′50.51′′E Gr, as provided by the SROCC projections, with the rate of vertical land motion inferred from InSAR data. In details, the measured rate of vertical land motion, which includes both GIA and tectonic components of both natural and anthropogenic origin, replaces the GIA contribution used in the original SROCC projections. Errors coming from the InSAR data are included in the overall uncertainty estimation.

3.4. Optical Satellite Images

We used the time series of optical satellite images Landsat 7–8 [105] acquired in the time span 2001–2019 and multispectral WorldView 2 acquired in 2011 to evaluate the horizontal migration of the sandy shorelines that underwent to significant retreat. The effective shoreline migration was obtained following the method described in Scardino et al. [80]. From all the selected images, the shorelines were digitalized in GIS environment for each year reported in the satellite images. In the coastal areas affected by several meters of coastal retreat, an estimation of the rate of total shoreline retreat was performed.

To this aim, we applied a linear regression of the shoreline position digitalized with respect to years of the satellite images. This approach involves the coastal retreat as a result of all contribution that cause the sediment movements and shoreline migration. Firstly, to assess the coastal retreat only due to SLR contribution, we assumed the investigated coasts to be in steady-state approach, without significant sediment movements and the shoreline migration caused only by land subsidence estimated by the geodetic measurements. Secondly, the eustatic sea-level component was added to the local land subsidence analysis to obtain horizontal displacements in function of the coastal slope. The latter has been extracted from DEM of LiDAR data. The most relevant shoreline changes have been observed on the coastal plain of Catania, the sandy coast of Siracusa and the coast of Vendicari, near Marzamemi (Figure 6).

To project the shoreline migration at 2050 and 2100 epochs, we incorporated the sea-level projections in the model reported in Scardino et al. [80]. Starting from the present-day coastline extracted from high-resolution satellite image of ESRI ArcGIS service (year 2019), the coastline in 2050 and 2100 was mapped.

3.5. Digital Elevation Model (DEM)

The present-day topography was obtained by the analysis of airborne LiDAR (Light Detection and Ranging) observation data collected by the Regione Sicilia in 2008–2009 [106] and georeferenced in EPSG 3004 Monte Mario Italy 2. The DEM was extracted with 2 × 2 m cell width and a mean vertical resolution of about 20 cm. The rates of VLM estimated by InSAR data have been represented as contour lines on the Digital Elevation Models (DEMs), by means of ArcGIS^® and Global Mapper software (www.globalmapper.com (accessed on 1 January 2021) [107]. 3D high-resolution maps of the investigated area were produced, on which the position of the present-day coastline and its potential position in 2100 as a result of relative sea-level rise are shown by contour lines and colour-shaded options, defining the upper limits of the exposed zones at potential sea flooding by 2100 for different climatic scenarios. DEMs were analysed to obtain the elevation of the coastal areas and the slope of the mobile coastal systems. To model the changes of the local topography at 2050 and 2100, elevations on DEM cells were corrected with the displacements of digital model VLM in 2050 and 2100. 3D high-resolution maps of the investigated area were produced for the topography in 2050 and 2100, on which the position of the present-day coastline and its potential future position in response to the relative sea-level rise, defining the upper limits of the exposed zones at potential sea flooding.

4. Results and Discussion

Since the last decade, new studies primarily based on IPCC projections, started to provide future sea-level rise assessments along the global coasts, generally without considering the contribution of VLM [1,108]. Lambeck et al. [30], Antonioli et al. [31]; Marsico et al. [79] and Antonioli et al. [32], constrained the spatially variable VLMs based on long-term geological data (mainly the elevation of the marine terraces of the MIS 5.5) to obtain RSLR scenarios for 2100 along low-lying coasts of the Italian peninsula. Recently, instrumental data from GNSS networks, InSAR observations and tide gauges have been considered in these projections for the realization of inundation hazard maps in several areas such as the coasts of California [34] and specific coasts of Italy [11,17,32,71,77,80,109,110], assuming that VLM rates measured by geodetic instruments will remain constant in the coming decades in a specific coastal zone.

Our results underline the importance of a multidisciplinary approach for the RSLR assessment to predict the future submersion of investigated coasts. The multi-temporal maps of the expected submersion for 2050 and 2100 have been obtained through an integrated analysis of SLR projections, VLMs and the detailed representation of the shoreline migrations on high resolution DEMs. The conceptual approach is based on the worst conditions for a mobile coastal system, in which the coastal resilience is not sufficient to contrast the sea-level rise. Such conditions are particularly evident in the coastal areas characterized by low sediment supply and retreatment of the river mouths (for example the Simeto river mouth, Figure 7). Moreover, strong subsidence could act locally, preventing the dynamical response of the coastal system. A negative contribution that affects the Southeastern Sicily coasts is due to the anthropogenic factors [111], which determined the decrease of sediment supply and fluid withdrawal. Obviously, in the attempt to propose a future scenario, it is necessary to make some hypotheses and, in our case, it has been assumed that the ground deformation is linear over time. This velocity rate is considered to be constant and equal to that obtained from 4 years of monitoring with SAR data. This assumption is considered valid as the GNSS data also confirm that over a much longer temporal period (more than 20 years). Furthermore, we showed how the vertical displacement retrieved by the two remote sensing techniques, e.g., InSAR and GNSS, are comparable as well confirming our starting assumption. Finally, considering the tectonic regime of the area, we can state that most of the observed vertical velocities recorded at the GNSS stations show tectonic stability, with the exception in the north of Catania, near the Etna volcano. On the other hand, the low rate of land subsidence recorded south of Catania, fit the GIA values for this region, which have been estimated between −0.5 and −1 mm/year [29].

Results highlight the contribution of sea-level rise to the shoreline retreat at a level lower than 35% (Figure 8). Below are reported the estimated shoreline changes for a set of coastal zones:

Catania coastal plain. A variable shoreline change was estimated between 6.34 and 3.23 m/year (with a linear regression of 4.78 ± 2 m/year). In the proximity of the Simeto, Gurnazza and San Leonardo river mouths a fast shoreline retreat, up to 10 m/year, was observed [111]. The southern part of the coastal plain influences large subsidence, reaching a rate up to 8 ± 2.46 mm/year in the Lentini area. The mobile coastal system showed phase of decrease of sediment supply, particularly over the last decades. The reduction of suspended load from Simeto, Gurnazza and San Leonardo determined a negative sedimentary balance, which is exacerbated by the intense beach-dune erosion.

The sandy coast of Siracusa. This coastal zone shows a rate of shoreline changes between 4.01 and 3.9 m/year (linear regression rate of 4 ± 1.9 m/year). The main evidence of shoreline changes occurred south of the Ciane river mouth. This area, which is characterized by salt-marshes, has been affected by a coastline retreat of about 70 m, with local breaching of the beach-dune system (Figure 6). The salt-marsh is subsiding at about 2 ± 2.46 mm/year, while the northern sector of this coast, displays higher subsidence rates mainly located in the Siracusa harbor, with velocities of about 5 ± 2 mm/year. The narrow beaches that separate the salt-marsh from the sea are not nourished by inland sediments due to dams that highly reduce the supply of sediments on the coast.

The coast of Vendicari. This area shows a rate of shoreline retreat between 0.95 and 0.4 m/year, with a mean rate of 0.68 ± 0.4 m/year. The InSAR analysis showed significant rates of local land subsidence, at about −4 ± 2.46 mm/year. Foredune system showed a given stability over the last few decades, while the dynamic of the coasts seems to be correlated to the longshore drift toward south.

Figure 9 and

Table 1. Relative sea-levels projections for RCP 2.6 and RCP 8.5 scenarios for 2050 and 2100 with respect to 2016, for different VLM rates observed in this study in south-eastern Sicily. Uncertainties define the 90% confidence limits. Different cases are considered for: (i) null vertical land motion (only the contribution of GIA already included in the SROCC projections is considered); (ii) mean VLM at −0.7 ± 0.35 mm/year derived by GNSS data and (iii) observed values of VLM estimated by InSAR analysis along the coastal zones reported in Figure 5.

VLM (mm/year)	RCP 2.6 2050 (m)	RCP 2.6 2100 (m)	RCP 8.5 2050 (m)	RCP 8.5 2100 (m)
Null VLM	0.15 ± 0.03	0.33 ± 0.06	0.20 ± 0.04	0.72 ± 0.13
GNSS: −0.7 ± 0.35	0.16 ± 0.03	0.36 ± 0.06	0.21 ± 0.04	0.75 ± 0.13
InSAR: −2.0 ± 2.5	0.20 ±0.04	0.47 ± 0.07	0.26 ±0.05	0.86 ± 0.13
InSAR: −3.4 ± 2.5	0.25 ± 0.04	0.59 ± 0.07	0.30 ± 0.05	0.97 ± 0.13
InSAR: −5.2 ± 2.5	0.31 ± 0.04	0.74 ± 0.07	0.36 ± 0.05	1.12 ± 0.13
InSAR: −7.6 ± 2.5	0.39 ± 0.04	0.94 ± 0.07	0.44 ± 0.05	1.32 ± 0.13
InSAR: −8.0 ± 2.5	0.40 ± 0.04	0.97 ± 0.07	0.46 ± 0.05	1.35 ± 0.13
InSAR: −8.7 ± 2.5	0.42 ± 0.04	1.03 ± 0.07	0.48 ± 0.05	1.41 ± 0.13
InSAR: −10.0 ± 2.5	0.47 ± 0.04	1.13 ± 0.07	0.52 ± 0.05	1.52 ± 0.13

show reference RSLR projections for RCP 2.6 and RCP 8.5 for the investigated coastal zones most prone to SLR in 2050 and 2100.

Because the GNSS stations are often too sparse and far from the investigated areas (Figure 2), we preferred to use the VLM rates obtained from the InSAR analysis to assess the local VLM trend and to calculate the expected RSLR values. It is worth noting that the VLM velocities obtained from the InSAR data coincides with the GNSS velocities at the individual stations with an uncertainty of ±2.5 mm. The expected flooding scenarios reported in Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15 are based on VLM values shown in Figure 5. Maps are referred to the RCP 2.6 (low emission) and RCP 8.5 (high emission) climatic scenarios reported in the SROCC Report (Oppenheimer et al., [26]), for which SLR is accelerated by land subsidence as estimated by geodetic analysis. The multi-temporal coastal positions projected on the high-resolution topography extracted from LiDAR data, allow to detail the potential flooding extension and the related impacts on the investigated coasts (

Table 2. Expected RSLR (m) and extension (km²) of the exposed flooded areas for 2050 and 2100 A.D. for RCP 2.6 and RCP 8.5 emission scenarios.

Area	RSLR RCP 2.6 (m)	RSLR RCP 8.5 (m)	RCP 2.6 2050 (km2)	RCP 2.6 2100 (km2)	RCP 8.5 2050 (km2)	RCP 8.5 2100 (km2)
Catania plain	2050: 0.15 ± 0.03	2050: 0.20 ± 0.04	1.21	2.634	1.762	4.092
	2100: 0.33 ± 0.06	2100: 0.72 ± 0.13
Augusta	2050: 0.31 ± 0.04	2050: 0.36 ± 0.05	0.275	0.668	0.433	1.659
	2100: 0.74 ± 0.07	2100: 1.12 ± 0.13
Siracusa	2050: 0.31 ± 0.04	2050: 0.36 ± 0.05	0.393	0.951	0.544	2.378
	2100: 0.74 ± 0.07	2100: 1.12 ± 0.13
Asinaro river mouth	2050: 0.47 ± 0.04	2050: 0.52 ± 0.05	0.015	0.021	0.017	0.043
	2100: 1.13 ± 0.07	2100: 1.52 ± 0.13
Vendicari	2050: 0.47 ± 0.04	2050: 0.52 ± 0.05	0.675	1.047	0.837	1.528
	2100: 1.13 ± 0.07	2100: 1.52 ± 0.13
Marzamemi	2050: 0.15 ± 0.04	2050: 0.2 ± 0.05	0.007	0.013	0.010	0.045
	2100: 0.33 ± 0.07	2100: 0.72 ± 0.13
Total extension			2.58	5.33	3.60	9.75

). For the RCP 8.5 emission scenario, the maximum extension of expected flooded area for 2100 A.D. in the investigated zone, is about 9.7 km².

The potential extension of the flooded area depends on the topographic features and the expected RSLR. The largest flooding is expected in the coastal plain of Catania and the Lentini area, where the complex topography and the San Leonardo channel can drive a marine ingression up to 5 km inland, both for the RCP 2.6 and 8.5 scenarios (Figure 10).

In the Augusta bay, the expected extension of the potential flooded area is relevant in the northern sector, reaching up to 500 m and 1 km inland for the RCP 2.6 and RCP 8.5 scenario, respectively (Figure 11). The most important flooding is expected for the coast of Siracusa, where the flooded area will possibly extend for more than 2 km inland in both RCP scenarios (Figure 12), involving the local salt-marsh, part of the city and the railways. The RSLR scenario for the pocket beach of Piccio, located in the proximity of Asinaro river mouth (Figure 13, Figure 14 and Figure 15), is less critical both for the RCP 2.6 and RCP 8.5 scenarios, but still with a partial submersion of the backshore.

The natural reserve of Vendicari is exposed to a potential flooding extending up to 1 km inland, involving the lagoons of Pantano Roveto and Pantano Piccolo (Figure 12). The expected RSLR in Marzamemi bay, may cause a flooding extension of about 300 m inland (Figure 13), determining the potential loss of the foreshore and backshore and leading to the abandonment of the historical fisherman village.

Concerning the role of the dune system in counteracting the SLR, field evidence and optical satellite images show that over the last two decades, the beach-dune systems have been subjected to a continuous retreat in consequence of multiple factors (

Table 3. Shoreline changes for different slopes (column C) for three sandy coasts located at Catania, Siracusa and Vendicari (Column A). The difference between the total observed shoreline retreat (due to the sea-level rise, the increasing intensity of extreme marine events, VLM and deficit on the sediment supply in a coastal system—column D) and the retreat due to RSLR only (column E, geometric horizontal retreat), defines the residual shoreline change (D–E). The coastal erosion triggered by RSLR is the predominant effect in the shoreline change process for the analysed gentle sloping sandy coasts which are retreating at up to 4.78 m/year.

A	B	C	D	E	F
Sandy coasts	VLM InSAR, (mm/year)	Coastal slope (degree)	Total shoreline retreat (m/year)	RSLR shoreline retreat (m/year)	Residual shoreline change (D–E) (m/year)
Catania	−8 ± 2.5	2.5	−4.78	−0.27	−4.51
		6.7		−0.09	−4.69
		9.45		−0.07	−4.71
Siracusa	−2 ± 2.5	0.7	−4.01	−0.46	−3.55
		5.38		−0.06	−3.95
		8.3		−0.04	−3.97
Vendicari	−4 ± 2.5	2.5	−0.68	−0.17	−0.51
		3.3		−0.13	−0.55
		10		−0.08	−0.6

), mainly, the sea-level rise, the increasing intensity of extreme marine events, VLM and deficit on the sediment supply in a coastal system. The high rates of land subsidence highlighted by InSAR, suggest that future scenarios could be driven by overtopping processes, rollover or dune breaching [112,113], probably similar in timing and architecture of resulting deposits, to retrogradation processes occurred during the marine transgression after the last glacial maximum (LGM) in different low-gradient settings [114,115,116]. In particular, the time series of satellite images of the sandy coasts at Siracusa between 2009 and 2020 (Figure 8, Table 3), show the occurrence of dune breaching near the salt-marshes over the last decades. Considering the expected accelerated SLR trend in the RCP 8.5 emission scenario, the beach-dune system will be more vulnerable in the next years with respect to the past decades due the combined effects of SLR, VLM and negative mass balance on the coasts.

If sea-level will continue to rise at accelerated rates, it is reasonable to expect a corresponding accelerated erosion that will strike severely the coastal system. The increasing coastal hazard will facilitate the marine flooding with subsequent impacts of the coastal zone, including saltwater contamination of surface and underground waters, wetland losses, and increased flooding which may lead to relevant socioeconomic impacts.

An additional factor, acting in combination with RSLR, is the occurrence of extreme marine events connected to storm and tropical-like system (Medicanes). Several models showed that the future storm events will be more intense than today [117,118,119,120], and the combined effect of wave flooding and sea-level rise could increase the coastal vulnerability, with loss of surfaces greater than those shown in Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15. Moreover, the intensity of tropical-like system is increasing in the Mediterranean [19,121,122,123], leading to extreme marine events that in conjunction with RSLR could determine several meters of inland flooding.

5. Conclusions

The methodological approach here presented to account for SLR and coastal land subsidence is transferable to other coastal zones and can be used to inform land planners and decisions makers that should take into account similar scenarios for a cognizant coastal management. The ongoing global climate change is affecting the inundation risks both through accelerating ice sheet melting (i.e., increasing the rate of eustatic SLR) and through more intense droughts, leading to unprecedented groundwater overexploiting and associated localized coastal land subsidence. Because climate change and SLR are posing unprecedented threats to the coastal environment, urbanization and population, multidisciplinary studies that include Earth Observations in combination with ground data and SLR projections will provide crucial information on the evolution of the coastal zones and their future expected shape.

Our results have shown that the area underwent to coastal retreat up to 70 m in a few years only while diffuse land subsidence is locally exceeding 10 ± 2.5 mm/year, due to compaction of artificial landfill, salt marshes and Holocene soft deposits. Given the ongoing land subsidence, the estimated high end of RSLR at 0.52 ± 0.05 m and 1.52 ± 0.13 m expected for 2050 AD and 2100 AD in the RCP 8.5 scenario, will lead up to about 10 km² of marine flooding. The inundation scenarios we have presented should be considered for a cognizant management of the coastal zone in response to the ongoing climate changes.

Major consequences of inundation risk for the coastal areas of south-eastern Sicily include saltwater contamination of surface and underground waters, accelerated coastal erosion, wetland and saltmarshes losses and, in general, increased flooding. These effects are all related with the ongoing global warming and may lead to unprecedented socioeconomic impacts. Our analysis can improve flood resilience plans for an area of the Mediterranean with significant populations and infrastructures exposed to flood risk that can be exacerbated by storm surges, earthquakes, tsunamis and volcano-tectonic activity.