Sea-Level Change along the Emilia-Romagna Coast from Tide Gauge and Satellite Altimetry

Abstract

Coastal flooding and retreat are markedly enhanced by sea-level rise. Thus, it is crucial to determine the sea-level variation at the local scale to support coastal hazard assessment and related management policies. In this work we focus on sea-level change along the Emilia-Romagna coast, a highly urbanized, 130 km-long belt facing the northern Adriatic Sea, by analyzing data from three tide gauges (with data records in the last 25–10 years) and related closest grid points from CMEMS monthly gridded satellite altimetry. The results reveal that the rate of sea-level rise observed by altimetry is coherent along the coast (2.8 ± 0.5 mm/year) for the period 1993–2019 and that a negative acceleration of −0.3 ± 0.1 mm/year is present, in contrast with the global scale. Rates resulting from tide gauge time series analysis diverge from these values mainly as a consequence of a large and heterogeneous rate of subsidence in the region. Over the common timespan, altimetry and tide gauge data show very high correlation, although their comparison suffers from the short overlapping period between the two data sets. Nevertheless, their combined use allows assessment of the recent (last 25 years) sea-level change along the Emilia-Romagna coast and to discuss the role of different interacting processes in the determination of the local sea level.

1. Introduction

Sea level lies among the 54 Essential Climate Variables (ECVs) that represent the Earth system’s key parameters, as defined and periodically assessed by the Global Climate Observing System (GCOS), in support of the United Nations Framework Convention on Climate Change (UNFCCC) and of the Intergovernmental Panel on Climate Change (IPCC). Direct observation of the sea level dates back to the 18th century [1] and this made possible realistic models for the sea-level variation and its interconnection with climate change. The ongoing sea-level change is the product of a continuous, mutable, and complex interaction among several components through the whole Earth system, combining natural and anthropic causes. It is commonly recognized that the increase in anthropogenic greenhouse gasses emission in the atmosphere, plus a small increment of natural solar irradiance, has led to a progressive effective radiative forcing (ERF) growth since the 1970s [2], hence to climate warming. More than 90% of the heat increase in the climate system, derived from the ERF, has been stored by the ocean, with the consequential thermal expansion and sea-level rise. At global scale, there is a general consensus that sea level has been rising at 1.7 ± 0.2 mm/year over the last century [2], and at 3.3 ± 0.5 mm/year over the last 25 years [3]. This is not a steady process since at both temporal scales, a positive acceleration (0.084 ± 0.025 mm/year²) was also confirmed [4,5]. Presently, there is the virtual certainty that sea-level rise and climate change will bring to society and the environment significant consequences, especially in low-elevation coastal zones where more than 620 million people live presently, a number that is expected to double by 2060 [6,7]. This enhances the climate-related concern about the sea-level rise impact at global scale and at coasts.

Satellite radar altimetry (SA) missions started in the early 1990s with the goal of providing, among other data, a fortnightly measure for the global sea surface height (SSH) in the reference system of the Earth’s mass center. Given its reference frame, this is commonly referred to as the absolute sea level (ASL). SSH represents the difference between the satellite altitude above the reference ellipsoid and the height of the satellite above the instantaneous sea surface, measured by transmitting microwave radiations toward the surface which are partly reflected to the satellite [8]. As a consequence of this configuration, SSH measurements are affected by several sources of uncertainties. Such uncertainties are resolved by compensating for the effect of instrumental drift, signal refraction passing through the atmosphere, Glacial Isostatic Adjustment (GIA), tides and atmospheric pressure influence whose correction is commonly known as the Dynamic Atmospheric Correction (DAC). However, in shallow water and at the coast SA data are influenced and contaminated by both the presence of land in the satellite footprint and the spatio-temporal scale reduction of atmospheric and oceanographic processes with respect to the open ocean [8,9,10]. The need for reliable coastal altimetry data has led over recent years to the development of new sophisticated altimetry products, where the improved reprocessing and dedicated corrections aim to extend valid measurement to within < 20 km from the coast [11,12,13].

On the contrary, tide gauges (TGs) provide a measure of the sea level with respect to the pier at which the TG is attached and for this reason it is commonly referred to as relative sea level (RSL). The largest TG database is the Permanent Service for Mean Sea Level (PSMSL), where almost two thousand TG time series are collected, periodically updated, and reduced to a common datum, producing the Revised Local Reference (RLR) dataset [14,15,16]. Developed for different purposes, i.e., the observation of the harbor tide, currently TG data used for long-term sea-level studies suffer from the contamination by the vertical land movement (VLM) of the pier itself. This problem was sorted out in the recent past by the co-location of GNSS (Global Navigation Satellite System) at the TG site. However, this practice only started since the late 1990s [17,18] and for this reason it does not cover the whole time series at most of the TG sites. Despite the long period coverage (since 19th century for some specific cases) TGs data are geographically limited because of their heterogeneous spatial distribution. TGs measurement has increased since the 1950s, primarily throughout the Northern Hemisphere [19]. This limitation generates a bias that may hinder accurate and reliable sea level consideration at global scale.

Despite SA and TG time series currently represent the two main source of information for the sea level, the contrast between their observational ability remains an unresolved point. RSL, measured by the TGs, reflects the local variations in sea level and it is influenced by the VLM effect and the coastal processes, differently from the SA ASL records that provide geocentric sea-level variations. SA and TG data sets are, moreover, characterized by low and high acquisition rate, respectively [11]. Only SA measurement allows observation of the variability on sea level from the open ocean towards the coastal area and it is considered essential for giving perspective at regional and global scale without resolving, however, some key coastal processes, which instead are recorded by TGs. RSL and ASL variation and VLM are linked by means of the Sea Level Equation (SLE) that, in its simplest form, reads [20]:

S = N − U

(1)

where S is the RSL variation, N represents the ASL variation and U is vertical displacement [21,22]. This equation tells that the sea-level change at the pier is the combination of the absolute seal level variation in a geocentric reference frame minus the vertical displacement at the TG coastal location.

At global scale, the drivers of sea-level change are primarily related to variations in steric contribution [23,24] and ocean mass component [25,26,27]. By contrast, at regional (e.g., the Mediterranean Sea) and local scale the sea-level change could be amplified or mitigated by several other drivers [28,29] such as VLMs, ocean dynamics (which become significant in semi- enclosed basins [30,31]) and the GIA (i.e., the visco-elastic crustal response to the inhomogeneous loading redistribution of ice melted water throughout continents and ocean [32,33,34]). The sea-level variability observed in the Mediterranean Sea is the consequence of several factors among which the circulation forced by heat and buoyancy flux at the surface, wind stress, and water exchange through the Strait of Gibraltar, varying periodically the water mass amount within the entire domain [30,35,36]. Regional sea-level variations are, in fact, influenced also by the balance between evaporation processes and precipitation–river runoff in the basin [37,38]. Wind stress can also lead to large Mediterranean sea level, generating periodic non-steric fluctuations due to mass transportation through Gibraltar from Atlantic [39,40]. A critical role of atmospheric forcing on sea-level change since the 1960s, not constant in time, has been documented for the Mediterranean Sea [41]. Indeed, several studies found a direct influence of inter-annual, decadal, and multidecadal climate fluctuations on sea level, corresponding to several tenths of mm/year [42,43]. This non-stationary variability is correlated with natural modes of the coupled ocean–atmosphere system such as El Niño Southern Oscillation (ENSO), North Atlantic Oscillation (NAO) and the Atlantic Multidecadal Oscillation (AMO) [44,45].

The Emilia-Romagna (E-R) coast is a highly urbanized, 130 km-long belt of low-elevated area (1.45 m maximum average elevation), located at the eastern margin of the Po Plain in Northern Italy (Figure 1). Coastal anthropization started in the 1950s [46], in the following decades the E-R coast became one of the most attractive tourist destinations in Italy and an important source for the regional economy. The combined effect of subsidence, sediment supply deficit due to the marked decrease in fluvial transport, removal of the dune ridge and increased anthropic pressure, cause the E-R coastal tract to be extremely sensitive to coastal flooding due to storm surges and sea-level rise [47]. In particular, a considerable portion of the E-R coastal area is currently below sea level and this would exacerbate the impact of future sea-level rise. As shown by Perini et al. [48] the combined effect of the IPCC worst scenario [2] and of the current local subsidence would cause, at 2100 with respect to 2012, a further extension (346 km²) of the areas below mean sea level. Also, for this reason, the determination of the relative sea-level (RSL, i.e., the local sea level relative to a benchmark on land) variation at local scale for the E-R coast is paramount for the management of the territory and of its fragile environment.

In this paper, we analyze the sea-level variation for the E-R coastal area by selecting three sites along the E-R coast: Porto Garibaldi, Marina di Ravenna, and Rimini (Figure 1). These are the only sites where a TG in operation currently exists. For each site we also identify the closest SA altimetry cell to put the RSL observed at TGs in the framework of the ongoing ASL for the region. TG times series cover different time spans, and this complicates the analysis. The main objectives are: (i) to verify the consistency of different instrumental sea-level records (TG/SA) at each site, by comparing in situ data with the latest CMEMS (Copernicus Marine Environment Monitoring Service) altimetry products (see Section 3); (ii) to observe and analyze differences among obtained sea-level data at the three sites, i.e., at the coastal scale, with a focus on the local variability, heavily influenced here by VLM; (iii) to assess the recent (last 25 years) sea-level change and if it can be assumed as a reference for the E-R region, also highlighting the variability of sea-level rates which strongly depends on the records length and (iv) to discuss this figure in the framework of similar studies at the basin-scale (Adriatic/Mediterranean), and of climate change-related processes. The determination of current sea-level rise at regional scale remains a complex task especially when only short time series are available. However, since this remains a crucial point in support of coastal hazard evaluation and related management/adaptation policies, efforts should be made to understand sea-level variability and to describe its intricacy from the coastal and basin perspective. The paper is organized as follows: Section 2 introduces the study area, Section 3 and Section 4 are dedicated to the description of data and analysis, respectively. Results are presented in Section 5 and discussed in Section 6. Finally, in Section 7, conclusions are drawn.

2. Study Area

The coastal belt of E-R is located in the eastern part of the Po Plain, south of the Po Delta (Figure 1). In geological terms, the Po Plain represents a foredeep basin suffering from natural subsidence. VLMs are mainly due to the consolidation of Quaternary alluvial deposits, with estimated rates along the coastal area of about −1 mm/year and −2.5 mm/year to the south and north of Ravenna municipality respectively, and twice as much at the Po river Delta [50]. Long-term VLMs due to tectonic component have been taken into account by Ferranti et al. [51] and Antonioli et al. [52,53] with rates in the order of −0.95 mm/year since the Last Interglacial (last 125 kyrs) for the northern E-R coast, while −0.22 ± 0.05 mm/year represents the glacio-hydro-isostatic contribution estimated by different GIA models [48]. Since the second half of the 20th century, the subsidence rate in the area has been markedly enhanced by the withdrawal of groundwater and gas extraction [54], with rates locally exceeding the natural subsidence by one order of magnitude. This caused land settlement locally exceeding 100 cm in a few decades [55]. The following gradual decrease in subsidence observed for almost the whole E-R coast is commonly attributed to the groundwater withdrawal cessation, occurred since the 1970s e.g., [54,56], and it has become more evident in the last two decades [55,57]. Current estimates of VLM based on the regional high-precision levelling network, regularly monitored since 1983 by the Regional Agency for Prevention, Environment and Energy of the E-R region (ARPAE), indicate average subsidence rates in the order of 5 mm/year for the last 15 years, decreasing to 3–4 mm/year in the period 2011–2016 [55]; however, peak values in the order of about 20 mm/year are still locally observed in the Ravenna coastal tract. The Coastal Geodetic Network is now the base for the definition of the orthometric heights along the E-R coast [58], integrated with InSAR techniques and GNSS measurements. Recently, Montuori et al. [59] integrating GNSS and multitemporal DInSAR techniques for the monitoring of land subsidence processes along the coastal plain, obtained average rates >3 mm/year, increasing to over 5 mm/year for Ravenna and some tracts along the SE coast. Rates of land subsidence for the Marina di Ravenna dock (from 1970 to present) based on geometric levelling, InSAR techniques and GNSS measurements, are compared by Cerenzia et al. [60].

From a morphological point of view, the E-R coast is represented by fine to medium-sand low-gradient dissipative beaches, 60% of which is the seat of hard protection structures [55]. The E-R coast is in general affected by low-energy waves, with Hs commonly below 1.25 m and microtidal regime [61,62]. The main storms come from the Bora (northeast) and Scirocco (southeast) winds, resulting in storm surge levels >1 m on a 100-year return period [48]. This is a huge concern for such low-elevation coastal zones, exacerbated by forthcoming climate change effects.

The E-R coast faces the northern Adriatic basin (Figure 1), a shallow sub-basin (a few tens of m of depth) with a low topographic gradient. The Adriatic is a narrow epicontinental shelf (800 by 200 km), communicating with the Ionian Sea through the Otranto Strait (Figure 1). The sea-level variability in the Adriatic Sea significantly differs from the other Mediterranean Sea sub-basins because of its setting [63,64,65]. It is dominated by cyclonic circulation driven by wind and thermohaline currents, mainly represented by three gyres centered in the southern, middle, and northern sector of the basin; northerly and southerly currents flow along the eastern (Eastern Adriatic Current, EAC) and western (Western Adriatic Current, WAC) coast, respectively (Figure 1) [64,65,66,67]. Previous analyses of sea-level rise in the altimetry era found an average rate of 3.2 ± 0.3 mm/year over the 1993–2008 for the Adriatic Sea, not uniform on the whole period, but positive (9 ± 0.5 mm/year) between 1993 and 2000 and negative (−2.5 ± 0.5 mm/year) between 2001 and 2008 [68]. Similar positive and negative rates before and after 2001 result for the (basin-average) steric component of sea level, which is correlated with Adriatic Sea level. For the northern sector of the Adriatic, rates of 4.25 ± 1.25 mm/year for the period 1993–2015 [69] and of 3.00 ± 0.53 mm/year for the period 1993–2018 [70] have been recently estimated.

3. Data

In situ sea-level signals are achieved from the three TGs currently operational along the E-R coast (Figure 1). These are installed at Porto Garibaldi (TG_PG), Marina di Ravenna, also known as Porto Corsini (TG_RA), and Rimini (TG_RN). In particular:

TG_PG has been operational since July 2009 and it is co-located with a permanent GNSS (Figure 1) station (GARI) that provides VLM information for this site. Different sources are freely available for geocentric surface velocity data at GARI [71]. In this work we considered the GARI time series produced with the INGV (Istituto Nazionale di Geofisica e Vulcanologia) solution by employing the GAMIT/GLOBK software and expressed in the IGS14 (International GPS Service) reference frame [72]. Only TG_PG and a portion of TG_RA are RLR data retrieved from the PSMSL.

At Marina di Ravenna, different TGs almost continuously collected relative sea-level data since 1873 [73,74], with a main gap from 1922 to 1933. Presently, data are managed by the Italian National Tide Gauge Network (Rete Mareografica Nazionale, RMN) which refers TG data to a local benchmark, located at the lighthouse, 150 m apart. TG_RA has a large gap between 2016 and 2019 and, for this reason, the time series considered in this work terminates in December 2015. A permanent GNSS station has operated since 1996 close to TG_RA [73] and is currently managed by the Department of Physics and Astronomy of the University of Bologna (no open data).

TG_RN has been in operation since July 2012, it is managed by Hera Group corporation, and it is not equipped with a GNSS antenna. Despite the shortness of the data and the lack of a co-located GNSS antenna, TG_RN is valuable since it provides the only datum for the southern portion of the E-R coast.

Both TG_PG and TG_RA operate with radar systems technology (the former coupled with a float system sensor, following the Intergovernmental Oceanographic Commission (IOC) directives [75]) while TG_RN operates with a differential pressure transducer. Since only float systems directly provide a SSH measurement, seawater density and gravitational acceleration are considered to convert pressure measured by TG_RN into sea level, while radar systems are supplied with dedicated hardware and software that directly convert the measures [76]. Regarding the SA, in this study we use the CMEMS daily gridded dataset (SEALEVEL_MED_PHY_L4_REP_OBSERVATIONS_008_051), having a spatial resolution of 0.125 degree x 0.125 degree [77]. This product is processed by the DUACS (Data Unification and Altimeter Combinations System) multi-mission altimeter data processing system, merging data from all the satellites available at a given time [77,78,79,80,81]. We selected the three closest pixels (SA_PG, SA_RA, SA_RN) to each of the TG (Figure 1), located at about 5, 13, and 9 km from Porto Garibaldi, Marina di Ravenna, and Rimini TGs, respectively. Each time series was downsampled from daily to monthly by computing the average for the whole samples of each month.

4. Analysis

To minimize the contrast between SA and TG time series, some procedures need to be considered, depending on the purpose of the study. DAC correction in SA (see Section 1) accounts for the low-frequency response of the sea surface to the atmospheric loading, known as the Inverse Barometer (IB) effect [82], and for the high-frequency (HF) response to wind and pressure forcing effect [83]. Since CMEMS products data are corrected by the DAC, the IB and HF effects must be subtracted from TG data for comparison purposes. However, the HF effect is negligible on periods longer than 10 days, hence no correction is needed when monthly averaged data are used [23]. The correction for the IB effect allows the subtraction of the atmospheric loading contribution from the observed RSL at TG. The IB effect correction considered in this work is based on the Wunsch and Stammer [82] formula:

IB(θ,λ,t) = −9.948 (Δrdry(θ,λ,t) – 1013.3)

(2)

where Δrdry(θ,λ,t) and −1013.3 are the local and global (standard) atmospheric pressure (in millibar) respectively, and −9.948 a constant that accounts for the water density and gravitational acceleration. The monthly mean time series of atmospheric loading at each TG was recovered from the ERA5 reanalysis [84] developed by the European Centre for Medium-Range Weather Forecasts (ECMWF). Furthermore, the application of the GIA correction on TG time series is considered of great importance when reconstructing the ASL, since no worldwide TG record can be assumed to be completely unaffected from the GIA effect [85]. However, since we are interested in studying the local RSL change with respect to the ASL change, and since we are considering short-lived time series, we can neglect to apply this correction.

Assessing the rate of sea-level variations is a complicated task exacerbated when time series are short and possibly affected by multi-annual and decadal oscillations, difficult to determine on the short term. However, we stick on standard trend assessment by means of ordinary least squares (using NumPy’s Python library [86]) to ease the comparison with previous results even though we are conscious that more sophisticated analysis could provide different results and that this approach neglects the autocorrelation of the time series. This weakness is compensated by the characterization of the periodic signals using two separate approaches, first we compute the Lomb–Scargle periodogram [87,88] (LSP) for each time series to provide a complete description of the frequency content by using the Lomb package [89] implemented in R [90]. Then, we apply the Empirical Mode Decomposition (EMD) method [91,92], a standard method for splitting non-linear and non-stationary time series into a limited sequence of empirically orthogonal “intrinsic mode functions” (IMFs) that describe cyclic modes not necessarily characterized by a constant amplitude nor phase. Differently from traditional spectral methods, the EMD is not requiring assumptions on the functional expression of the regression model. This allows us to decompose each time series in sequence of IMFs plus a residual. The latter, normally monotonic, helps to define how a linear trend is a good model for the non-periodic component of each time series. The significance of the resulting trends is discerned by means of standard test [93] that checks the probability of the null hypothesis to be true. This test is reported beside each trend with its significance code: *** p < 0.001; ** p < 0.01; * p < 0.05; ° p < 0.1; ‘ ‘ (blank) p < 1.

5. Results

5.1. Sea-Level Assessment

Satellite altimetry time series (SA_PG, SA_RA, SA_RN) are shown in Figure 2a. The visual comparison among the three SA time series highlights excellent coherence that is confirmed by the correlation coefficient [94] r = 0.98. To better characterize the frequency content of these time series we compute the Lomb–Scargle periodogram (LSP; Figure 2b) for each of them [87,88]. The time series are dominated by the annual and semi-annual oscillation, both related to the seasonality, with peak-to-peak amplitude varying in the order of 150–300 mm (Figure 2a). The resulting LSPs (Figure 2b) also confirm the coherence between the three sites and demonstrate the presence of multi-annual oscillation whose most relevant periods are about 4 and 12 years (see Section 6.1 for detail). LSP significance analysis indicates a small probability (<10%) for the peak at about 12 years to be artefact even though the large sample spacing prevents an accurate definition of the period itself. Conversely, the random peak probability is larger than 10% for the one at about 4 years.

Standard linear regression provides a coherent rate of 2.8 ± 0.5 mm/year over the timespan January 1993–May 2019 (details in

Table 1. Linear trend rate and acceleration estimated in this work from SA data. Symbol *** indicates that the model is strongly significant with p < 0.001 and symbol * refers to p < 0.05.

	SA Rate (mm/year)	SA Acceleration (mm/year2)
Porto Garibaldi (January 1993–May 2019)	2.8 ± 0.5 ***	−0.3 ± 0.1 *
Marina di Ravenna (January 1993–May 2019)	2.8 ± 0.5 ***	−0.3 ± 0.1 *
Rimini (January 1993–May 2019)	2.9 ± 0.5 ***	−0.3 ± 0.1 *

). These values are smaller but comparable to the global average (in the range of 3.1–3.4 mm/year [95]) over the same timespan. The rate assessed in our analysis is smaller than what observed for the northern sector of the Adriatic for the period 1993–2015 (4.25 ± 1.25 mm/year) by Vignudelli et al. [69] even though the two error bars partly overlap, and consistent with what observed for the period 1993–2018 (3.00 ± 0.53 mm/year) by Mohamed et al. [70]. We in fact observe that when fitting the same time series with a second-order model (quadratic, Figure 2a), we get a slight but statistically significant negative acceleration of −0.3 ± 0.1 mm/year². This local negative acceleration, in contrast with that at the global scale for which a slight positive acceleration (0.084 ± 0.025 mm/year²) has been observed [4], is discussed in Section 6.3.

As summarized in Section 3, the three TG time series do not cover the same timespan and are shorter than SA (Figure 3 and Figure A1). This limits the comparison between different sites and, for the shorter one, the reliability of the RSL rate assessment. The longest time series is TG_RA (January 1993–December 2015), showing a RSL rate of 5.8 ± 0.8 mm/year (

Table 2. Linear trend rate (for raw data and IB-corrected data) for the tide gauges considered in this work. The significance of the resulting trend varies according to the symbol beside each rate (significance code: *** p < 0.001; * p < 0.05; ° p < 0.1; ‘ ’ p < 1).

	TG Rate (mm/year)	TG (IB-Corrected) Rate (mm/year)
Porto Garibaldi (July 2009–December 2019)	1.4 ± 2.8	2.8 ± 2.2
Marina di Ravenna (January 1993–December 2015)	5.8 ± 0.8 ***	5.5 ± 0.8 ***
Rimini (July 2012–September 2020)	−7.7 ± 3.6 *	−5.1 ± 3.0 °

). This value is consistent with previous studies e.g., [60,68] and significantly higher than the ASL rate observed offshore (January 1993–May 2019, Table 1), as a consequence of the strong vertical land motion expected in this coastal sector. The rate from TG_PG has a large uncertainty that leaves the null hypothesis (absence of trend) likely, while from TG_RN it is largely negative. Despite the large disparity in the resulting RSL rate is also a consequence of the different time spans, we note that the correlation coefficient, restricted to the overlapping periods, is 0.80 between Porto Garibaldi and Marina di Ravenna, 0.94 between Porto Garibaldi and Rimini, and 0.77 between Rimini and Marina di Ravenna.

Since one of the goals of this study is the comparison between RSL and ASL to better understand the different governing mechanisms and observational approaches, we apply the IB correction to the raw TG data to compare them with SA time series (corrected for DAC). This correction gives a valid reprocessing of the sea-level variability without the meteorological effect influence [96]. We can observe that the IB correction leads to a reduction, in amplitude, for most of the extrema (Figure 3 and Figure A1) and this reflects in smaller variance for the residual of the OLS regression and, indeed, in a smaller error associated with the trend (Table 2). As expected, the IB-corrected TG time series result in different RSL rates and smaller uncertainties. We also note that the quadratic fit for TG_RA returns a positive acceleration (0.7 ± 0.3 mm/year²), not present in the ASL data (Figure 3).

The SA time series have been shortened to match each of the three TG time series (Figure 3). Also in this case the correlation coefficient is high (0.89) for Porto Garibaldi e Rimini and moderate for Marina di Ravenna (0.66). We notice that the rates for the shortened SA time series (

Table 3. Linear trend rate comparison between SA and TG for the common timespan covered. The significance of the resulting trend varies according to the symbol beside each rate (significance code: *** p < 0.001; * p < 0.05; ‘ ’ p < 1).

	SA Rate (mm/year)	TG (IB-Corrected) Rate (mm/year)
Porto Garibaldi (July 2009–May 2019)	−3.1 ± 2.0	1.2 ± 2.3
Marina di Ravenna (January 1993–December 2015)	3.5 ± 0.6 ***	5.5 ± 0.8 ***
Rimini (July 2012–May 2019)	−3.3 ± 3.2	−7.6 ± 3.7 *

) differ from what observed for the complete time series (January 1993–May 2019). In particular, in the case of SA_PG and SA_RN, for which we consider only the last part of the time series, rates turn negative. In details, for SA_PG (timespan July 2009–May 2019) the trend is negative, while the corresponding rate of RSL is positive (Table 3). This result for SA linear trends is not unexpected since it is consistent with the observed negative acceleration (Table 1) that suggests a decrease in ASL rate of about −5 mm/year over the removed 16 years (1993–2009). For the same reason, at SA_RA the trend for the period January 1993–December 2015 is larger than the one computed for the whole SA timespan, and lower than the RSL rate (Table 3). Finally, for the case of Rimini, SA_RN trend is negative (Table 3) but the large error bar demonstrates that the null hypothesis (no trend) could not be rejected with 0.1 < p < 1. In other words, the large error bar demonstrates the low significance of the linear model and that the rate could also be null or even positive. As in the case of altimetry, the analysis of shorter-term time series also modifies the TG trend values (compare Table 2 and Table 3), except for TG_RA whose time series is fully covered by the altimetry and therefore does not need to be cut. Furthermore, an important contribution to the local sea-level variation at different sites is represented by the vertical movements (VLM) component included in the TG data. This will be taken into account and discussed in Section 6.2.

5.2. Periodic Signals

The large uncertainty in estimated rates, and in their variation over different time spans, suggests analyzing the periodic components of the time series and to define their relevance in terms of amplitude and frequency to better constrain the characteristics of the sea-level change. With an approach similar to that used in Section 5.1 for SA, we analyze the TG time series by means of LSP and EMD.

Results are shown in Figure 4. LSP analysis (Figure 4a) for the whole period range confirms the prominence of the oscillation at 1 year (especially for longer time series), while at shorter periods (below 1 year) different peaks emerge between 0.2 and 0.5 years (i.e., between 2 and 6 months). For the case of Porto Garibaldi and Rimini we notice that the IB correction enhances the signal for the annual oscillation and that a clear peak emerges at 0.85 years (10 months). This partly contrasts with what observed for the LSPs computed for the SA time series in which the semi-annual oscillation emerges more clearly. The ascending termination of the periodograms at the longest period (right side) for the case of TG_RN and TG_RA could reflect the presence of periodic signals at periods longer than the timespan covered by this (very short) time series. We also notice that IB correction reduces the power amplitude at longer periods.

To better discern the periodic signals from the long-term trend, we also applied the EMD analysis to each IB-corrected time series. From Figure 4b, we note that at high frequency (IMF1–2), large non-sinusoidal oscillations exist, and these show an overall coherence between the three sites. This suggests a regional scale origin. IMF3, whose periodic signals range between 1.5 and 3 years, is not coherent during the overlapping timespan, pointing to a site-related effect. IMF4 exists only for Marina di Ravenna and Rimini but these time series overlap for a very short timespan and this prevents any comparison. IMF5 at even longer periods emerges only for Ravenna and it demonstrates an almost sinusoidal oscillation with a period of about 10.7 years. Residuals confirm what observed from linear regression (Table 3) and introduce some more complexity: Marina di Ravenna shows short-lived positive acceleration around year 2000, Porto Garibaldi is almost linear with two minor inflections, while Rimini shows a parabolic shape contrasting with a linear model. For Rimini, the shortness of this time series limits the distinction between a branch of a parabola or of a sinusoid with period of about twice the length of the time series.

6. Discussion

The above data analyses provide a complex picture for the sea-level variation along the E-R coast. Actually, different data and processing, and the large literature available for the region suggest several perspectives for discussing and interpreting such complexity. In this section, we propose a comprehensive interpretation of our results, also in the framework of previous studies. We focus on factors affecting sea-level variability, on the role of the VLM in the E-R coast, and on the difficulties of providing a figure for the current or recent rate of sea level for the area.

6.1. Sea-Level Time Series and Rates for the E-R Coast

To verify the consistency of different instrumental sea-level records at the coastal scale, TG and satellite altimetry time series were analyzed at three sites along the E-R coast. Altimetry time series at the three sites turn out to be very consistent, showing basically the same rate of sea-level rise (Figure 2a; Table 1) and variability in the last 25 years, also denoted by a perfect correlation among them (r = 0.98). On the other hand, the results provided by the three TGs are more variable in terms of trend and variability both among them and with respect to the altimetry. This can be the consequence of some critical factors. First, the TG time series cover different time spans (Figure 3), and this can cause discrepancies in the observed rates (Table 2 and Table 3). Indeed, the timespan for which the three TG time series overlap is very short, thus preventing a reliable comparison. Secondly, rate assessments based on short records are particularly sensitive to oscillations especially at their extremes, thus they are not suitable for deriving climatological-related information in the long term [2].

It is a common practice to reproduce the sea-level variation over time by a linear model. However, this approach sometimes forces the observed phenomenon into a model that is too simple, neglecting periodic and aperiodic variations of different time scale and size that can bias the resulting trend. The addition or subtraction of a few samples at the time series extrema can cause considerable variations of the trend rate as recognizable by comparing values listed in Table 1 and Table 2 with those of Table 3. For the case of the Adriatic Sea, the rate of sea-level rise strongly varies with time because of the chosen record length [68,85,97,98], and by the contamination of seasonal and inter-annual variability [99]. To constrain this, we combined EMD analysis and LSP periodograms. LSP, among the expected annual and semi-annual oscillations, demonstrates two oscillations at about 4 and about 12 years for SA data, while the interpretation of TG data is more complex given the short timespan for TG_PG and TG_RN. TG_RA shows a more complex periodogram in this band, suggesting that several events contribute to hinder the SA peaks; however, EMD analysis returns a period signal at about 10.7 years. We can speculate that this oscillation and the SA peak at about 12 years are comparable with the 10 years-oscillation at Mediterranean scale observed by Bonaduce et al. [100].

The oscillations that influence the Adriatic Sea are mainly due to cycles that act at annual (stronger energy signal), semi-annual (6 months) and inter-annual (5 years and higher) scales [101], as also demonstrated in our data (Figure 2a and Figure 4a). The whole Mediterranean basin is also influenced by oscillations at lower frequencies driven by natural modes, which produce intermittent anomalies in the long-term rate. As analyzed by Galassi and Spada [45] the NAO and AMO power are mainly concentrated at a period of about 8 and 9 years respectively; they also pointed out that the AMO strongest signal is found at multidecadal periodicities, i.e., between 50 and 70 years. Moreover, one of the most energetic modes of variability in the Adriatic basin is linked to the periodic occurrence (approximately 20 years) of opposite phases of the NAO and AMO [45] which produce a large scale zonal atmospheric gradient [102] and a sea-level non-steric fluctuation directly linked to mass transport through the Strait of Gibraltar driven by wind [39,103]. The last occurrence of this phenomenon was determined as the cause of the large positive anomaly recorded in sea-level signal during 2010–2011 [40,70] by all Mediterranean TG stations [100]. Indeed, the right site upgoing termination of the LSPs (Figure 2b and Figure 4a) could be an indicator of a longer period oscillation that could not be quantified here because of the length of the time series.

One further cause of the variability in SL rate and corresponding uncertainties is the different data processing and different corrections applied. An example of the above cited complexity is the Marina di Ravenna site for which several previous studies exist, based on the same TG time series but using different methodologies and adopted corrections, as shown in Figure 5 where the results of some of the most representative studies conducted on TG_RA are summarized (see also

Table A1. Linear sea-level trend estimates (and related corrections applied) from some of the most representative studies conducted on the Marina di Ravenna TG time series. Time frames considered by various authors are indicated in brackets.

	Rate (mm/year)	Type of Data Considered
Antonioli et al., 2017	8.3 ± 0.3	monthly mean (2000–2013)
Bruni et al., 2019	1.25 ± 0.16	monthly mean (1873–2016), detrended for non-linear VLM
Cerenzia et al., 2016	8.5 ± 0.2	annual mean (1897–2014)
	7.7 ± 0.3	monthly mean (1970–2013)
	9.8 ± 1.3	annual mean (1990–2013)
	(2.2 ± 1.3)	(detrended for linear VLM)
Fenoglio-Marc et al., 2012	6.5 ± 1.5	monthly mean (1993–2008), IB correction, deseasoning
Masina and Lamberti, 2013	8.6 ± 0.5	annual mean (1896–2011)
	3.1 ± 0.7	annual mean (1896–1950)
	10.3 ± 8.0	annual mean (2000–2011)
Tsimplis et al., 2012	3.1 ± 1.3	annual mean (1897–1921)
	1.5 ± 1.3	annual mean (1905–1921), normalized for Trieste
	8.4 ± 1.1	annual mean (1937–1972)
	(7.1 ± 0.8)	(normalized for Trieste)
	3.9 ± 0.6	daily mean
Zerbini et al., 2017	1.22 ± 0.32	annual mean (1896–2012), IB correction, detrended for non-linear VLM
	1.05 ± 0.54	annual mean (1934–2012), IB correction, detrended for non-linear VLM

). There is a relation between the time series length and the error size of the trend resulting from OLS regression, i.e., the error scales inversely with the number of samples of the analyzed time series [104]. Rate for longest datasets generally agree, but only when similar data corrections are applied: Zerbini et al. [73] and Bruni et al. [74] removed from data the non-linear VLM resulting from GNSS data, while Cerenzia et al. [60] accounted for the linear VLM (from various geodetic data), and Tsimplis et al. [41] normalized the data with respect to the TG time series of Trieste (considered to be a tectonic stable site). Rates computed from long time series without VLM correction are significantly higher [41,60,105] but still relatively grouped. By contrast, rates from short time series (i.e., less than 20 years, see [60,68,105,106]) denote extreme variability and inconsistency between them regardless of the corrections applied, e.g., IB correction and deseasoning, whose influence is relatively low here with respect to the crustal movements corrections. This highlights how the natural variability and the contribution of VLMs strongly influences the RSL signal at the TGs, preventing the estimate of a univocal sea-level trend valid at the short-, middle-, and long-term scales. Given the short time spans and the described complexities, we thus cannot establish if and how the observed trends could be representative of the long-term rate of RSL rise. Also, the time series analyzed in this work (Figure 3; Table 2) suggest that a not constant rate of subsidence might play a crucial role in the observed RSL change at this time scale.

Regarding the comparison between the two different techniques, the analyzed SA and TG records have been downsized to cover only the overlapping timespan at each site and this causes differences in the estimated rates at the three sites (Table 2 and Table 3). The results of this comparison are thus expected to improve in the next years/decades when the available time series will be longer for the analyzed sites. An 18-year- [107] or 15-year-long [108] recording is needed to obtain a good correlation. It has been shown, in fact, that over intervals shorter than 10 years the inter-annual variability affects more the rate of global mean sea level from TG than from altimetry [68]. However, despite these considerations, we obtain very high correlation between SA and TG sea-level time series at the three studied coastal sites. In particular, the three downsized SA time series show a strong correlation with the IB-corrected TG time series (Figure 3), especially at Porto Garibaldi and Rimini (r = 0.98). This confirms the good quality of the chosen SA product that, as suggested by Vignudelli et al. [69] for the northern Adriatic Sea, is less affected by degradation at the coast than previous SA data. As a matter of fact, different linear rates result for TG and SA at each of three study sites (Table 3). Since SA provides a coherent rate along the E-R coastal area (Table 1), an important contribution to the spatial variability of RSL at different sites should derive from the VLM component included in the TG data, and this is the focus of the next section.

6.2. The Role of Vertical Movements in Sea-Level Determination

In the E-R coastal area (and on a great portion of the Northern Adriatic coast, Italian side) subsidence is widespread and occurs at high rates [109], while vertical deformations caused by other phenomena as GIA are more than one order of magnitude smaller than subsidence that here dominates the VLM affecting RSL change [60,68,69,105,106,109]. Presently, different integrated techniques, such as high-resolution levelling, e.g., GNSS and InSAR, provide data on VLM with increasing accuracy [57,58]. However, the small size and short-lived variation of vertical deformation processes causes this information to be not always comparable/homogeneous. Furthermore, the common practice of co-locating a GNSS station at the TG site is relatively recent and many TGs do not have associated VLM data with the necessary precision.

As mentioned in Section 3, TG_PG has co-located GNSS station (coded name GARI) since Jul. 2009. GARI station provides well constrained VLM data well reproduced by a linear trend with rate equal to −2.8 ± 0.1 mm/year. At Marina di Ravenna, the VLM rate estimated for the GNSS station located close to TG_RA, estimated through the Gipsy software over the timespan July 1996–December 2015, is −5.6 ± 0.2 mm/year [110]. These VLM rates, when compared with rates of RSL at each site well explain a relevant portion of the sea-level change in the considered time frames (Table 2) while the ASL contribution seems low or negligible. At Rimini, the rate of RSL in the time frame Jul. 2012–Sep. 2020 (the only one available at TG_RN) is negative (−5.1 ± 3.0 mm/year; Table 2), in contrast with what observed for TG_RA and TG_PG. Here, a GNSS station (ITRN) operated by the Coastal Geodetic Network [57] (located at 2.6 km from the tide gauge, at 44°2’54”N 12°34’55.4”E) shows a VLM rate almost null (0.00 ± 0.8 mm/year) for the period 2011–2016 (for further details see [111] and references therein).

Integrated data on vertical deformations (topographic levelling, GNSS and InSAR) of selected benchmarks of the E-R Coastal Geodetic Network and data previously acquired in the framework of the Regional Network for subsidence monitoring [112], provide gradually decreasing subsidence rates for all the coastal sites under study in the last 25 years [55,57] (

Table 4. Ranges for the rate of subsidence estimated over different time frames for each of the study sites from selected benchmarks of the E-R Coastal Geodetic Network [111].

	1992–2000 (mm/year)	2002–2006 (mm/year)	2006–2011 (mm/year)	2011–2016 (mm/year)
Porto Garibaldi	7.5–10	7.5–10	2.5–5	2.5–5
Marina di Ravenna	7.5–10	7.5–10	5–7.5	2.5–5
Rimini	5–7.5	5–7.5	5–7.5	0–2.5

). These data are gridded at 100 m × 100 m and reproduced as isokinetics of VLM at the regional scale; associated error is estimated as ±2 mm/year [57]. For the area that contains the TG data, rates summarized in Table 4 can be considered representative of the subsidence. Actually, these data agree with the rate observed by GNSS at Porto Garibaldi and Marina di Ravenna. For Rimini, subsidence data suggest that this phenomenon is gradually decreasing but not inverted, as the negative rate of RSL rise might suggest.

As a parallel approach, we use the SLE (Equation (1)) to compare the observed sea level or ground variation with the one obtained by the combination of the over two data. It must be pointed out, however, that the accuracy of this method is frequently larger than 1 mm/year whatever its application [113,114,115]. Despite the coarseness of this method (20 years of data should be required to yield a high-precision estimate, [115]) several authors have applied it to regional studies, since understanding VLM for sea-level research is of farthest significance e.g., [21,116,117,118,119,120,121,122,123,124].

For the case of Porto Garibaldi, the three variables that form Equation (1) are available from the different observables considered for this work over the common period July 2009–May 2019. This allows us to evaluate the different estimates obtained through the two approaches, i.e., by comparing the rate obtained from the linear regression of the direct observation (Table 3) with those computed by using Equation (1) (hereafter computed). Visual inspection of Figure 6 suggests a strong correlation between the pairs (SA, U+S) and (TG, N-U), this confirmed by r = 0.90. By contrast, the correlation between (N-S) and the GNSS is low (r = 0.37). This follows the different accuracy of GNSS data (below 0.1 mm/year) with respect to TG data contained in N-S. The trend values related to the computed time series (

Table 5. U, N and S resulting from the SLE (Equation (1)) where N is the ASL variation, U is vertical displacement and S is the RSL variation (values in mm/year). NA = not available. The significance of the resulting trend varies according to the symbol beside each rate (significance code: *** p < 0.001).

	U	N	S
Porto Garibaldi (July 2009–May 2019)	−4.3 ± 3.1 ***	−1.7 ± 2.3	−0.2 ± 2.1
Marina di Ravenna (January 1993–December 2015)	−2 ± 1	NA	NA
Rimini (July 2012–May 2019)	4.3 ± 4.9	NA	NA

) slightly differ from those of the observed one (Table 3). We note that only for U the null hypothesis (no trend) could be rejected. Figure 7 shows that the observed and computed trends match within their error bars, despite the differences between their central values, suggesting a good coherence between the different observations. In particular, the trend observed for N (from U+S) is lower in value than the ASL derived from SA_PG (−3.1 ± 2.0 mm/year) even though their error bars overlap (Figure 7). The trend for S is almost null and its error bar shows good overlap with the rate of RSL given from TG_PG (1.2 ± 2.3 mm/year; Figure 3, Figure 6, and Figure 7). U is affected by larger uncertainties with respect to the GNSS time series, with uncertainties lower in value (−2.9 ± 0.1 mm/year in the time frame Jul. 2009–May 2019) and completely contained in the U error bar (Figure 7).

For the case of Marina di Ravenna and Rimini, only S and N can be directly determined, since GNSS instruments are not co-located at these TGs. Then, Equation (1) can be applied to indirectly estimate U. Figure 8 shows the time series (SA-TG) for Marina di Ravenna and Rimini (with their related errors) whose rate (Table 5), according to Equation (1), should represent (U = N − S). These values are quite unexpected since subsidence at Marina di Ravenna is known to be larger (see Table 4); at Rimini this contrasts with the observed subsidence (Table 4). However, if U for Marina di Ravenna is computed on the reduced timespan for which GNSS data exist (July 1996–December 2015) the rate of VLM is −4.7 ± 0.6 mm/year, smaller than that provided by the GNSS located near TG_RA (−5.6 ± 0.2 mm/year). The two results do not match within their error bar for only 0.1 mm/year. As mentioned above, this GNSS is in the vicinity of TG_RA and, the larger the distance between GNSS and TG, the smaller the representativeness of the real VLM at the TG. It has been observed that, when large heterogeneities exist, a few hundred meters of distance between TG and GNSS stations can cause decoupled movements [18]. This seems to be the case for the Marina di Ravenna harbor where Cerenzia et al. [60] found different rates of subsidence at very short distances, likely due to the weight effect of the structures themselves (e.g., the embankment with respect to the dock and the lighthouse). Fenoglio-Marc et al. [68] reports similar results for the Marina di Ravenna site, where GNSS observations provide a higher VLM rate with respect to the one computed from SA and TG difference, attributing it to anthropogenic causes, to the distance between the TG and the closest SA grid point and to the corrections applied on the RSL time series.

At Rimini, the context appears completely different with respect to Marina di Ravenna and Porto Garibaldi. Despite the high correlation between TG_RN with SA_RN, at Rimini the RSL shows a negative trend (Figure 3), although the small timespan covered by TG_RN time series gives at this data very low statistical significance (Table 2 and Table 3). The time series (SA-TG) shows a positive trend for U (Table 5) indicating a “proxy-derived” local effect of land uplift (Figure 8), in contrast with the VLM rates (0.00 ± 0.8 for the timespan 2011–2016) provided by the ITRN station. GNSS values agree with recent subsidence monitoring activities using InSAR (Table 4). The error associated with InSAR measure is 2 mm/year, suggesting that for the period 2011–2016, some portion of the Rimini area that includes the GNSS and TG stations may undergo opposite VLMs (uplift) with respect to other sectors of the E-R coastal belt. At Rimini, the GNSS station is located 2.6 km apart from the TG and this makes the data comparison even less effective than for the case of Marina di Ravenna [18]. Measured subsidence at Rimini has markedly decreased in the last two decades [57] and it cannot be excluded that VLM could reverse in the next decades; unfortunately, to date there is not enough data to support this hypothesis and more years of data acquisition are required.

6.3. The Sea-Level Trend Acceleration Issue

As anticipated in the introduction, at global scale the sea level has been rising at a rate of 3.3 ± 0.5 mm/year [3], with a positive acceleration of 0.084 ± 0.025 mm/year² [4] over the altimetric era (1993–present). Global maps show that the sea-level rise is not spatially uniform, with a diffused pattern of null rate and even showing a sea-level fall, see [125]. Then, it is not surprising that this heterogeneity also reflects in different local sea-level acceleration as for the case of the small spot target of this study in which we observe, for the same period, a negative acceleration of –0.3 mm/year². Acceleration could be the consequence of global short-lived phenomena, as for the case of the Pinatubo eruption in 1991 [126] or of climatological events. Moreover, given the length of the SA time series, accelerations could be also a consequence of the contribution of multidecadal oscillations for which only a portion of it is sampled. Because of this consideration, the observation of a full cycle (at least) is required to properly model the oscillation itself, determine its origin, and its plausible relation with climate change. Unfortunately, since the satellite altimeter era only started in 1993, and long-term TG observations are contaminated by an unknown non-linear VLM (measured only since co-located GNSS acquisition started, i.e., in the late 1990s), it is currently complicate to properly attribute to the ongoing observed regional sea-level variations to climatological origin. In our case, we observe a localized negative acceleration for the E-R coast, whose interpretation should be put in the broader context of the Mediterranean Sea in which a marked spatial variability of sea-level trend is observed, for instance, by Bonaduce et al. [100]. The ASL negative acceleration observed at the E-R coast in recent decades (Figure 2a; Table 1) has been previously noted by other authors. After a positive acceleration of sea-level rise observed in the Mediterranean Sea during the 1990s [127] since the early 2000s, in fact, the sea level stopped rising [128,129,130]; this seems to match well with the quadratic fit shown for SA in Figure 2a. As previously explained in Section 6.1, ocean–atmosphere oscillations are responsible for inter-annual and interdecadal cyclic sea-level variations in the Mediterranean; the main engine of these phenomena has been associated, at a global scale, with natural modes (e.g., the NAO and the AMO) inducing temporary reductions or magnification of sea-level trend amplitude at the basin to regional scale [103]. Besides this inter-annual variability, changes in sea-level acceleration on multidecadal time scales have been documented by several authors for the whole 20th century which is generally characterized by a greater rate in the first half, followed by an overall negative acceleration [43,131,132,133]; however, the processes causing this variability, and their link to climate change, are not completely understood yet.

We showed in Section 5.1 that the SA time series, shortened for matching with the TG time series (Figure 3), provide rates differing from what observed for the complete SA time series (January 1993–May 2019, Figure 2a). In particular, in the case of SA_PG and SA_RN, rates turn negative (Table 3; although for Rimini the hypothesis of no trend could not be rejected), consistently with the observed negative acceleration (Table 1). For the Marina di Ravenna site, the overlapping timespan (1993–2015) is long enough to obtain, from the quadratic regression, a significant acceleration in the time series (Figure 3). Namely we get 0.7 ± 0.3 mm/year² for TG_RA and −0.2 ± 0.2 mm/year² for SA_RA. The opposite signs suggest that two different phenomena govern the two accelerations and a good candidate could be VLM that, as mentioned above, is governed here by subsidence. From Equation (1), the VLM acceleration is the difference between those of the two terms mentioned above, and it results equal to −0.9 ± 0.4 mm/year². Actually, subsidence is slowing down in the region and locally in the surrounding of TG_RA as resulting from Table 1, and this should reflect in a positive acceleration in the ongoing negative VLM. We can then argue that what we are observing at TG_RA is something different from the effects of the local subsidence or, more likely, other phenomena are governing the acceleration at TG_RA.

At the regional scale, the decrease of subsidence measured by the E-R Coastal Geodetic Network over recent decades (see Table 4) is a favorable indication in terms of reducing future flooding vulnerability for this area. However, as already pointed out by Cerenzia et al. [60] and Fenoglio-Marc et al. [68], and discussed at Section 6.2, different rates of VLM are available, depending on the precise site of measurement and not necessarily representative of the VLM at the TG location. We thus remark the need for local and repeated VLM measurements at the TG site to correctly interpret RSL series and for detecting possible ASL positive or negative accelerations.

7. Conclusions

In this paper, we have studied the sea-level variation along the E-R coast by focusing on monthly sea-level observations from the three available TGs located at Porto Garibaldi, Marina di Ravenna, and Rimini. Data are compared with CMEMS SA data from the three grid points closer to each of the TGs. The combined use of altimetry and tide gauge data (supplemented by GNSS acquisition) is considered, in fact, a promising approach, in terms of precision and cost-effective implementation, to better define various physical processes in the coastal domains and to anticipate the impacts of future rise in sea levels [10,11,22,69,134,135,136,137].

Our results show that the ASL time series (from SA) are coherent along the coastal tract, providing a rate of 2.8 ± 0.5 mm/year for the timespan Jan 1993–May 2019, smaller but comparable to the global average (3.3 ± 0.5 mm/year). Moreover, we note in the E-R coast a negative acceleration of −0.3 ± 0,1 mm/year² that contrasts with the positive global acceleration. The comparison with the IB-corrected TG time series suffers from the short overlapping period between the two different datasets (Jul. 2009–May 2019 at Porto Garibaldi; Jan 1993–Dec 2015 at Marina di Ravenna; Jul. 2012–May 2019 at Rimini). Despite their shortness, the SA and TG time series show high correlation coefficient (0.89) at Porto Garibaldi and Rimini, and a moderate one (0.66) at Marina di Ravenna. The lower correspondence achieved at Marina di Ravenna, as well as the difference between ASL and RSL rates for this site (Table 3), are mostly associated with the local rate of VLM. Actually, VLM can hamper the detection of signals associated with sea-level change and related spatial variations [22]. For all the three analyzed sites, this remains an important contribution to quantify. It is remarked the need for co-located monitoring of the ground motions at the TG sites, since the determination of VLM is a key element in understanding how sea level has changed over the past century and how future sea levels may impact coastal areas [22].

Finally, as already observed for other coastal areas [11,13,69], the good match between altimetry and in situ observation confirms that high-standard satellite product can provide consistent data even at short distance from the coast. Linking satellite altimetry measurements with tide gauge, by bridging the open-ocean measurements with those in proximity of the coast [12] is part of a new challenge in monitoring sea-level rise along coastal tracts. In this view, the maintenance and implementation of observation stations for continuous monitoring of sea level in the coastal zone is a fundamental tool for coastal management and defense strategies.