Revisiting Vertical Land Motion and Sea Level Trends in the Northeastern Adriatic Sea Using Satellite Altimetry and Tide Gauge Data

Abstract

We propose a revisited approach to estimating sea level change trends based on the integration of two measuring systems: satellite altimetry and tide gauge (TG) time series of absolute and relative sea level height. Quantitative information on vertical crustal motion trends at six TG stations of the Adriatic Sea are derived by solving a constrained linear inverse problem. The results are verified against Global Positioning System (GPS) estimates at some locations. Constraints on the linear problem are represented by estimates of relative vertical land motion between TG couples. The solution of the linear inverse problem is valid as long as the same rates of absolute sea level rise are observed at the TG stations used to constrain the system. This requirement limits the applicability of the method with variable absolute sea level trends. The novelty of this study is that we tried to overcome such limitations, subtracting the absolute sea level change estimates observed by the altimeter from all relevant time series, but retaining the original short-term variability and associated errors. The vertical land motion (VLM) solution is compared to GPS estimates at three of the six TGs. The results show that there is reasonable agreement between the VLM rates derived from altimetry and TGs, and from GPS, considering the different periods used for the processing of VLM estimates from GPS. The solution found for the VLM rates is optimal in the least square sense, and no longer depends on the altimetric absolute sea level trend at the TGs. Values for the six TGs’ location in the Adriatic Sea during the period 1993–2018 vary from −1.41 ± 0.47 mm y⁻¹ (National Research Council offshore oceanographic tower in Venice) to 0.93 ± 0.37 mm y⁻¹ (Rovinj), while GPS solutions range from −1.59 ± 0.65 (Venice) to 0.10 ± 0.64 (Split) mm y⁻¹. The absolute sea level rise, calculated as the sum of relative sea level change rate at the TGs and the VLM values estimated in this study, has a mean of 2.43 mm y⁻¹ in the period 1974–2018 across the six TGs, a mean standard error of 0.80 mm y⁻¹, and a sample dispersion of 0.18 mm y⁻¹.

1. Introduction

Sea level is a crucial index of climate change [1]. Its variation reflects critical processes involving our environment of both natural and anthropogenic origin. Examples include subsidence, eustatism, thermal expansion of seawater, and melting of land-based ice sheets and glaciers. Other phenomena (e.g., oceanic circulation changes due to modified atmospheric circulation at synoptic scale and other coastal processes such as storm surges) that have a high impact on coastal settlements may also contribute to worsen the rising sea level scenario in the mid to long term [2,3].

The traditional technique for measuring sea level is to use a tide gauge (TG), which observes sea level height changes relative to the solid land surface where it is grounded. As a relative measurement, it accounts for both sea level change and land vertical displacement. Satellite radar altimetry (ALT) provides an alternative way to measure sea level changes from space [4], supplying the absolute part of the sea level measurement, and dependent only on the height of the sea surface in the geocentric reference frame, i.e., the absolute sea-level change, called sea surface height (SSH). To close the circle, Global Navigation Satellite System (GNSS) devices such as Global Positioning System (GPS) supply measurement of the land surface absolute geocentric height above the Earth’s center of mass. Indicated by $\dot{s}$, $\dot{g}$, the relative and absolute sea-level change rates measured by TGs and satellite altimetry, respectively, and $\dot{u}$, the rate of geocentric vertical land motion (VLM), the relation $\dot{u}=\dot{g}-\dot{s}$ holds true.

The sea level anomaly (SLA), the variable part of SSH, is usually derived by subtracting the mean sea level (MSL), obtained by averaging the local SSH, from the SSH itself. The satellite-borne altimeters measure the absolute sea level change up to a few kilometers from the coast [5]. As satellite altimeters have a limited ability to measure sea level height in proximity to the coast, where TGs are usually settled, the VLM affecting the TG measurements is more easily estimated with GPS. GPS co-located with TG allows separation of land- and sea level vertical changes. When GPS instruments are permanently installed at a site and operated regularly, continuous, real-time sampling of that site’s motion is obtained (CGPS).

The three independent measuring techniques can be used in synergy to make an integrated study of sea-level variability in order to identify significant trends over a long time. Tide gauges alone cannot distinguish between rising of the sea surface and subsidence of the land. Radar altimetry measurements, which are independent of land, can replace GNSS when not available at tide gauge places. This is under the assumption that other trend sources (e.g., corrections in altimetry) are negligible. This can be achieved only by using datasets with processing optimized for climate studies.

Many studies have focused on sea level rise and land subsidence throughout the world. Cazenave et al. [6] computed the simple difference between relative sea-level changes measured at tide gauges and nearby absolute sea-level changes from TOPEX/Poseidon satellite radar altimetry and found a striking correlation with vertical crustal movements time series obtained with Doppler Orbitography and Radiopositioning Integrated by Satellites (DORIS) measurements, co-located with TGs. Since then, other authors have used altimetry-TG sea-level difference techniques in several regions of the globe: Fennoscandia (northern Europe) [7,8], southern Europe [9], the Mediterranean Sea [10,11], the Adriatic Sea [12], and globally [13]. Our study focuses on the Adriatic Sea, an ideal laboratory for understanding sea-level rise, as long time series of coastal sea level from tide gauges are available, as well as altimetry and coastal altimetry datasets, and for which a comprehensive review of all aspects of sea level research has recently been compiled by Vilibić et al. [14]. Many civil settlement and commercial activities bound to the shipping, trading, tourism, chemical, and naval industries are spread along the low-lying coastal plains, often threatened by storm surges on the northeastern side, near Venice [15]. Thus, the study has a twofold benefit, as it focuses on a strategic region with regard to natural hazards, commercial and civil relevance, tourism, and vulnerability to climate change, and at the same time the region has been monitored for a long time, and plenty of literature is available. The positions of the six TGs used in this study, along with a schematic map of the Adriatic Sea, are shown in Figure 1.

Sea level variability in the study area has been investigated at different time scales using tide gauges, satellite altimetry, and GPS data. Most of these studies start from the same data but use different customized processing and time periods, thus making it difficult to compare the results.

Some studies using TGs have been performed in the Adriatic Sea, including in Venice, Trieste, and several Croatian sites (Rovinj, Zadar, Bakar, Split, Dubrovnik) [16,17], but a regional and integrated assessment of sea level change is still limited. Broader analyses focused on sea level trends of southern Europe and the Mediterranean Sea using TGs and altimetry do include the Adriatic Sea region [9,10,18]. Such studies use approaches like the one adopted here, but owing to different time periods and aims, it is difficult to perform meaningful comparisons. The research to date has demonstrated that trends are site-dependent and the largest differences are arguably introduced by local causes.

Cazenave et al. [19], using the last operational product generated within the Copernicus Climate Service (C3S), provides a global estimate of 3.3 ± 0.5 mm y⁻¹, which means a global mean sea level rise (GMSLR) of about 8 cm during the altimetry era. This altimetric product adopts the climate-quality processing methodology of the European Space Agency (ESA) Sea Level Climate Change Initiative (SLCCI) project and extends the GMSLR time series until February 2019 (26 years).

However, MSL does not rise (or fall) uniformly, and regionally the geocentric sea level rates can differ considerably from the GMSLR [20]. In the Mediterranean Sea, there is marked spatial variability of MSL rates (Figure 3), with the entire area of the Adriatic Sea exhibiting positive sea level rates above the GMSLR, with the highest values observed in the northern part. Despite different time spans and methods of comparison, the difference in sea level rates around Venice available in the bibliography is not significant, with values of 5.9 ± 1.3 mm y⁻¹ during 1993–2008 using data from the Radar Altimeter Database System (RADS) [12]; 4.18 ± 0.92 mm y⁻¹ using data from Archiving, Validation and Interpretation of Satellite Oceanographic data (AVISO); 3.40 ± 0.99 mm y⁻¹ using the SLCCI dataset during 1993–2013 [18]; and 4.25 ± 1.25 mm y⁻¹ during 1993–2015 with SLCCI data [21]. The results around Trieste showed a sea level rate of 5.9 ± 1.3 mm y⁻¹ during 1993–2008 [12] using the along-track RADS product, similar to Venice in the same period, supporting a uniform sea level rate in the Northern Adriatic Sea at that time. Later, lower values (4.01 ± 1.07 mm y⁻¹) using the AVISO product (1993–2014) were shown, and even much lower (2.76 ± 1.40 mm y⁻¹) using the SLCCI V1 product (1993–2013) [18]. A reassessment made using the AVISO product (1993–2012) showed 3.19 mm y⁻¹ [22]. The SLCCI V2 product during 1993–2015 showed sea level rates largely lower than previous estimates, with a value of 1.14 ± 1.35 mm y⁻¹ at the nearest gridded ground point in Trieste [23]. This rate is significant, although with larger uncertainty. The cause of the lowest value of sea level rate is still unclear. In this case, there is evidence of a larger difference between RADS, AVISO, and SLCCI products, in addition to the differences in terms of time series length and temporal span.

Subsidence variability along the Adriatic Sea coast is dominated by short timescale processes, mainly of anthropogenic origin. In Venice, VLM rates varied in the range 1–4 mm y⁻¹ at the lagoon basin scale in the period 2000–2015, with the notable exception of some centimeters per year detected at the three lagoon inlets [24]. The most probable cause of such a high rate of displacement is the construction of structures housing MOSE, the system of barriers intended to insulate Venice from floods by isolating it from the sea during high waters. In Marina di Ravenna, 110 km south of Venice along the western coast of the Adriatic Sea, the withdrawal of ground water and other fluids from underground has greatly enhanced the VLM of natural origin, inducing a combined rate of subsidence up to 30–40 mm y⁻¹ over several decades of the 20th century [10]. Likewise in Venice, in the same period, the subsidence of both origins has been estimated to be around 150 mm in total [25]. Trieste is considered to have not been subjected to sensitive subsidence in the last century [26]. However, a recent study reported the existence of a GPS station named TRIM very close to the tide gauge [18] where, over a period of almost 8 years (2007–2014), an uplift of 0.94 ± 0.10 mm y⁻¹ was observed. In the same study, an analysis of the 18-year GPS time series of Marina di Ravenna (PORT) revealed a subsidence of −5.67 ± 0.10 mm y⁻¹ (1996–2014). Positive relative VLM in the central and southern east coast of the Adriatic Sea with respect to the northeastern coast was noticed in the second half of the last century [27] using only TG observations. The latter were in Dubrovnik, Split, and Bakar, and the methodology used was a comparison of sea level records of two tide gauges, obtaining the movement of the second with respect to the first, assuming that both experienced the same absolute sea-level change. Unfortunately, no GPS measurements near Bakar are available, and the sea level (SL) record at the Permanent Service for Mean Sea Level (PSMSL) ends in 2013, so we did not include it in our study.

The recent availability of outstanding altimetry datasets with specific coastal processing (SLCCI and C3S) made it possible to revisit climate-scale sea-level variations in the Adriatic Sea, and they were used here. This study is based on previous work by Kuo et al. [7], who used a new method of combining long-term TG records with satellite altimetry measurements of sea-level change rates. With respect to the traditional approach, in which altimetry and TG observation time series are differentiated at each TG to obtain the local vertical land motion (VLM) solution independently from other locations [6], Kuo et al.’s methodology uses a linear system in matrix form populated from several TG and altimetry time series and solved simultaneously. The linear system is a realization of an inverse linear model where absolute vertical land motion rates are constrained with relative sea-level change rates obtained by differentiating some couples of TG sea-level time series. Constraints were implemented in the linear system through the Lagrange multipliers (LMs) technique [28]. This approach, applied to the Adriatic Sea, allows us to obtain a twofold objective: First, the best estimates of the absolute VLM rate at each TG can be obtained in an ordinary least square (OLS) sense, with accuracies and errors determined from the overall consistencies of all TG and altimetry time series considered. Second, the uncertainty in the VLM results is much lower than in previous studies [29], due to the accurate constraints of absolute vertical motion derived from long-term TG–TG relative sea level rate estimates. The constraints are valid under the assumption that the absolute sea level change rate is the same in both TGs paired in each constraint. This limitation was not overcome by Wöppelmann and Marcos (W&M) [9], who nevertheless further refined the methodology formulated in [7], and still remains today.

Our study stems from previous works [7,9] and focuses on a small basin of the Adriatic Sea that the authors are familiar with, and tries to overcome the limitation posed by the structure of the constraints, generalizing the usefulness of the methodology to regions with unrestricted absolute sea level rates.

The rest of the paper is organized as follows. Section 2 describes the TG and satellite altimetry data used in this study. The methodology outlined in [7] and refined in [9] is explained briefly at the beginning of Section 3, and in the rest of the section we formulate our idea of going beyond the intrinsic limitations of the TG–TG constraints. The obtained results are given in Section 4 and discussed in Section 5, where the conclusions of this study are also summarized.

2. Data

2.1. Altimetry Sea Level Datasets

Absolute sea level change rates were calculated from gridded datasets of altimetry SLA provided by two processing chains. The first relies on the SLCCI project funded by ESA [30]. The CCI initiative is ESA’s response to the need for systematic monitoring of the global climate system and the essential climate variables (ECV) as expressed by the United Nations Framework Convention on Climate Change (UNFCCC), directed at contrasting the impact of anthropogenic activities on the global climate. ESA’s CCI Programme aims to exploit the full potential of long-term consolidated global terrestrial observation archives to realize state-of-the-art datasets of ECVs in collaboration with the international scientific community.

The Copernicus Climate Change Service (C3S) provides the second SLA altimetry processing chain. Copernicus is the European Union’s Earth Observation Programme and offers information services based on satellite Earth observation (EO) and in situ data [31]. The C3S actualizes the Copernicus effort in the climate research field, combining observations of the climate system with the latest science. As an example, Figure 2 shows the Trieste record of the two altimetric datasets at the chosen grid point. Data were deseasoned by fitting 12- and 6-month harmonic curves, then the dynamic atmospheric correction (DAC) was re-added to make them comparable with TG observations. The general pattern of the two altimetry records is similar, but differences are evident. Both also have a similar linear correlation coefficient with the Trieste TG record: 0.89 for SLCCI and 0.90 for C3S.

2.1.1. SLCCI Altimetry SLA

The global sea level essential climate variable (ECV) v2.0 is a merged multi-satellite product covering the period January 1993 to December 2015. It consists of time series of gridded monthly SLA means covering the whole globe, and other variables. The SLA grids are calculated after merging all altimetry mission measurements together into monthly grids with a spatial resolution of 1/4 of a degree [32]. The SLCCI dataset is an improved set of reprocessed satellite-based sea level products intended to provide a reference for climate studies. It merges data from all altimetry missions: TOPEX/Poseidon, Jason-1, Jason-2, ERS-1, ERS-2, GeoSat Follow-On (GFO), Envisat, SARAL/AltiKa, and CryoSat-2. The SLCCI project has developed consistent altimeter corrections in order to produce a homogeneous and stable global sea level climate data record (CDR), designed to be a reference for climate-related sea level studies. Information on the SLCCI ECV SLA product can be found in [20,33]. Unfortunately, the latest SLCCI ECV gridded product data stop with 2015. However, the operational production of the climate-oriented global sea level product has now been taken over by the C3S [34]. The main difference is that all available satellites have been included in the SLCCI product, whereas a stable number of two altimeters is used for the C3S product; this contributes to increased stability of the sea level record, especially on a regional scale, where some MSL jumps can be introduced in an “all-satellite merged” product when a new mission is added to (or removed from) the constellation [34].

2.1.2. C3S Altimetry SLA

As noted on the SLCCI product web page [35], C3S, through the Climate Data Store, has taken over production of the altimetry SLA global and regional datasets, with updated temporal and spatial coverage. For the aims of this study, we downloaded the gridded daily SLA means over the Mediterranean Sea at 0.125° resolution, covering 26 years (01/1993–12/2018). C3S provides this state-of-the-art, climate-oriented SLA dataset. Up-to-date altimeter standards are used to estimate the SLA with a mapping algorithm specifically dedicated to the Mediterranean Sea. The sea level product is a time series of gridded sea surface height and derived variables obtained by merging two satellite altimetry measurements. It is generated by the Data Unification and Altimeter Combination System (DUACS) processing system of the French Space Agency (Centre National d’Etudes Spatiales, CNES) and Collect Localisation Satellites (CLS, a subsidiary of CNES), and includes data from several altimetry missions. The product mainly focuses on retrieving the long-term variability of the ocean, which is only obtained using a stable altimeter constellation and homogeneous corrections and standards in time. The later constraints are addressed by using a two-satellite constellation during the whole altimeter period [36].

Monthly means were obtained from daily means, averaging available daily values in a month at each node of the grid. This CDR is designed to be an up-to-date extension of the SLCCI SLA dataset to the present.

2.1.3. Dynamic Atmospheric Correction and TOPEX-A Drift

The original two altimetry SLA datasets were formed by removing the sea level variability due to the atmospheric wind and pressure forcing. Such a contribution is calculated by modelling the effects of wind drag and atmospheric loading at the sea surface and is usually referred to as dynamic atmospheric correction (DAC). As the objective of the study is to compare the SL time series at tide gauge, which includes the effect of the atmospheric wind and pressure, and the nearby altimetry SLA, from which the DAC was removed, we have two choices: subtract the DAC from the TG time series, or re-add the DAC to the altimetry datasets. We opted for the latter, as it is more consistent to rebuild the dataset that was initially deconstructed. CNES, through its reference portal for the altimetry AVISO+, supplies a consistent DAC correction for the C3S altimetry SLA product based on the MOG2D-G finite element model [37]. The C3S DAC MOG2D-G model is forced by the European Centre for Medium-Range Weather Forecast (ECMWF) Re-Analysis (ERA-Interim) model reanalysis data until April 2002 (TOPEX/Poseidon/ERS-1/ERS-2), and by ECMWF model analysis after that date [34], with the same temporal and spatial coverage as the C3S SLA. Similarly, SLCCI supplies consistent DAC for its SLA dataset, always using the MOG2D-G model, but forced by the ECMWF ERA-Interim re-analysis dataset. Thus, the DAC used to de-correct the C3S SLA (C3S DAC) was assembled with the SLCCI DAC from 01/1993 to 04/2002 and the AVISO+DAC after 04/2002.

The TOPEX-POSEIDON mission suffered from a drift of the TOPEX-A instrument in 1993–1998. The current drift has been estimated through a sea level budget closure approach and by comparison with tide gauges. Reprocessing of the TOPEX-A measurements is currently being finalized by a joint initiative of CNES and the NASA Jet Propulsion Laboratory (JPL), but the final step of having the along-track data of TOPEX-POSEIDON corrected for this drift is still far away; it is believed that this drift is uniform at regional scale and the global estimate can be applied regionally [34]. The status of the TOPEX-A instrumental drift has been described in a dedicated publication [38]. Because of the uncertainty of the correction to be made for the TOPEX-A drift in the SLA data, neither the C3S nor the SLCCI datasets were corrected for it in this study.

2.2. Tide Gauge Sea Level Dataset

Six tide gauges are considered in this study: Venice (TG name: VENEZIA), Venice off-shore (TG name VEPTF), Trieste, Rovinj, Split, and Dubrovnik (TG names: TRIESTE, ROVINJ, SPLIT, DUBROVNIK). VENEZIA and TRIESTE are among the TGs with the longest sea level records in the world. The SL time series of the historic TG in the city center of Venice dates back to 1872. The VENEZIA TG is operated by the Venice Tide Forecast and Early Warning Center (Centro Previsioni e Segnalazioni Maree, CPSM) of Venice Municipality. CPSM owns a network of TGs composed of several stations inside the Venice Lagoon, three at the lagoon’s inlets and one off-shore, installed on the Acqua Alta off-shore platform (45°18′51.29″ N, 12°30′29.69″ E) of the National Research Council of Italy (Consiglio Nazionale delle Ricerche, CNR). The ongoing offshore record, VEPTF, started in 1974. It is particularly used in routine forecasts and early warning of high waters in central Venice, since tides usually occur at the platform at levels similar to those of Venice, but about 50 min earlier, thus offering the possibility of alerting citizens of significant differences with expected levels. A significant part of the VENEZIA SL record, and all of the Acqua Alta platform, was made available by CPSM as hourly time series. The VENEZIA SL record has been reconstructed by merging records from CPSM and PSMSL. We will indicate SL records reconstructed from different sources with an asterisk after the TG name. The principal characteristics of the TG records are reported in

Table 1. Tide gauges considered in Adriatic Sea. Column 1: location name; column 2: TG name; columns 3–4: latitude and longitude; column 5: percentage of available data; column 6: total timespan; column 7: length of record in years. Asterisks in tide gauge (TG) names indicate that the sea level (SL) record has been reconstructed from different sources.

Location	TG Name	Lat (° N)	Long (° E)	Data (%)	Time Span (Years)	Record Length (Years)
Venice	VENEZIA *	45°25′51.45″	12°20′13.39″	97	1872–2018	147
Venice off-shore	VEPTF	45°18′51.29″	12°30′29.69″	100	1974–2018	45
Trieste	TRIESTE *	45°38’50.00″	13°45’33.90″	89	1875–2018	145
Rovinj	ROVINJ	45°05’01.18″	13°37’44.86″	99	1955–2018	64
Split	SPLIT *	43°30’23.88″	16°26’18.44″	100	1952–2018	67
Dubrovnik	DUBROVNIK	42°39’28.40″	18°03’38.84″	99	1956–2018	63

.

As an example, Figure 3 shows four records from different stations inside the Venice city center, covering periods from 1872 to 2019. Two are available from the Permanent Service for Mean Sea Level (PSMSL) [39,40], VENEZIA (PUNTA DELLA SALUTE) and VENEZIA (S. STEFANO). A third record (VENEZIA2) was kindly provided by CPSM; it is formed by different datasets from the same historic tide gauge of VENEZIA PUNTA DELLA SALUTE (1923–2020), initially operated by the National Public Works Office (Genio Civile), whose data were distributed to the Local Administration Authority, and whose ownership passed to the Hydrographic Bureau of the Magistrato alle Acque (Magistrate for the Waters, the public works office of northeastern Italy) [41], and more recently to the Istituto Superiore per la Protezione e la Ricerca Ambientale (Italian Institute for Environmental Protection and Research (ISPRA)), which also inherited the scientific and technical functions of the Hydrographic Bureau of the Magistrato alle Acque. The fourth (VENEZIA1) is available online [42] from ISPRA. Together, the four records cover a longer period than they would separately, and their quality has been checked and assured by the respective institutions. For this study, the four records were inter-calibrated by eliminating the relative biases calculated in overlapping periods, finally forming a complete record from the partial ones, with priority given to the most affordable of each record set. The same technique was adopted for the records of Split (SPLIT RT MARJANA and SPLIT–GRADSKA LUKA) [39] and Trieste, whose merged record was obtained from the PSMSL record of TRIESTE and the TRIESTE3 record made available by Fabio Raicich, CNR-ISMAR (CNR—Institute of Marine Sciences) [43].

The time series of the various SL records available for VENEZIA*, TRIESTE*, and SPLIT* TG stations are shown in Figure 3, Figure 4 and Figure 5.

As the VEPTF record has gaps (e.g., 1993 is almost completely missing), we used the SL record registered in another TG of the CPSM, namely Venezia Diga Sud Lido, to integrate the missing periods of the VEPTF record. This procedure is described in [21]. The TGs near Venice cited to date are well maintained and constantly monitored. They report the SL height with respect to the local Venice datum, the Zero Mareografico di Punta della Salute (ZMPS); Punta della Salute is the official zero-height reference against which sea level heights are historically measured in the Venice Lagoon. It corresponds to the mean sea level of 1884–1909 (conventionally referred to 1897). To compare with satellite altimetry records, already organized as monthly time series, TG records available only as hourly time series were reduced to daily averages by applying the Doodson X0 filter [44], following the Intergovernmental Oceanographic Commission guidelines [45]. Then, the monthly means were calculated from the daily means. TG records from the PSMSL were already available as monthly mean time series. The mutual Pearson’s linear correlation coefficient for the six SL records is always higher than 0.87, as shown in Figure 6.

2.3. Geocentric Surface Vertical Velocities from CGPS

Several sources are currently available online for geocentric surface velocity data and trends from CGPS of selected locations, in particular near TGs: Système d’Observation du Niveau des Eaux Littorales (SONEL)/Université La Rochelle (ULR) [46] and Nevada Geodetic Laboratory (NGL, University of Nevada) [47], among others. Often in SONEL both ULR and NGL solutions are available for the same CGPS station. Along with websites, other specific sources of information are available from local and national public agencies.

Table 2. Geocentric surface vertical velocities at four locations from CGPS stations. Column 1: station name; columns 2–3/4–5/6–7: VLM rate and time span calculated by NGL/SONEL/ISPRA; column 8: pooled means of given rates for the four CGPS stations.

CGPS	NGL (mm y−1)	Span (Year)	SONEL (mm y−1)	Span (Year)	ISPRA (mm y−1)	Span (Year)	Pooled Mean (mm y−1)
PSAL VENEZIA	−1.70 ± 0.80	2014–2020	−	−	−1.46 ± 0.09	2010–2015	−1.59 ± 0.65
TRIE TRIESTE	−0.52 ± 0.45	2003–2020	0.20 ± 0.26	2003–2013	−	−	−0.25 ± 0.52
SPLT SPLIT	0.45 ± 0.68	2004–2012	−0.25 ± 0.34	2004–2012	−	−	0.10 ± 0.64
DUBR DUBROVNIK	−1.99 ± 0.80	2000–2012	–	–	–	–	−1.83 ± 0.70 1
DUB2 DUBROVNIK	−1.94 ± 0.89	2012–2020	–	–	–	–

¹ This value is the pooled mean of DUBR and DUB2.

reports the trend and error values of four of the tide gauge locations used in this study (VENEZIA, TRIESTE, SPLIT and DUBROVNIK), collected from both websites and published reports. For VENEZIA, we include solutions made available online by NGL and those published by the Istituto Superiore per la Protezione e la Ricerca Ambientale (ISPRA), which maintains and conducts periodic high-precision leveling of the geodetic network benchmarks around the CGPS station of VENEZIA PSAL [48]. PSAL is almost co-located with the VENEZIA PUNTA DELLA SALUTE TG, also owned by ISPRA, which results in very good agreement with the sea level data of the TG owned by CPSM. For SPLIT, data from the CGPS station of SPLT were acquired. For more information on this CGPS station, see [49], but here it is worth mentioning that SPLT is, with PSAL in VENEZIA, among the few CGPS co-located with TGs. The TRIE CGPS station, associated with the TRIESTE TG, is around 6.5 km from the TG, over the hills surrounding Trieste, at an altitude of about 280 m a.s.l. Finally, the DUBROVNIK TG is linked to two CGPS stations, DUBR and DUB2, which cover complimentary periods of time from 2000 to 2020. They are 1 km from the TG and 400 m a.s.l. VLM values are sometimes very different from center to center, and often calculated on a limited and variable timespan. Nonetheless, combining trend values from different centers and/or periods of time is the best option to increase the time coverage. This is the case of PSAL and DUBROVNIK, whose longest timespan is formed using vertical velocities calculated by different centers or the same center but using different CGPS and periods. Such measurements have been combined using a pooled mean algorithm, which weights values and errors to give a unique estimate without the need for original time series data. Following [50], the pooled mean and its standard deviations are defined as ${\mu }_{\mathrm{p}}=\frac{{{\displaystyle \sum }}_i{n}_i{\mu }_i}{{{\displaystyle \sum }}_i{n}_i}$ and ${\sigma }_p=\sqrt{\frac{\left({{\displaystyle \sum }}_i\left[\left(n_i-1\right){\sigma }_i^2+n_i{\mu }_i^2\right]-\left[{\mu }_p^2{{\displaystyle \sum }}_i{n}_i\right]\right)}{{{\displaystyle \sum }}_i{n}_i-1}}$ respectively.

3. Synergistic Use of Tide Gauge and Satellite Altimetry Data to Estimate Vertical Displacement Rates

In this section, we look for estimates of the vertical displacement rates at the six tide gauges listed above. To do so, we follow the approach that was first used to estimate vertical crustal motion velocities in Fennoscandia [7], and then in the Great Lakes and Alaska [8]. It relies on the traditional approach described in [6] of deriving crustal vertical velocity rates by subtracting relative tide gauge sea level data from neighboring geocentric satellite altimetry observations. The inverse theory [28] is used to estimate the parameters of the (linear) model that best describes the mutual relationships between relative sea level rates, absolute sea level rates, and vertical land motion at the tide gauges simultaneously. Here, we briefly report on the approach to facilitate the reader’s understanding.

3.1. Linear Inverse Problem with Constraints

The rate of absolute vertical land movement at tide gauge i is given by the difference between the absolute and relative sea level change rates at the same place. Using a notation similar to that adopted by Kuo et al. in [7], we define $g_i$ as the sea level height as measured by satellite altimetry near the $i$th tide gauge, and $s_i$ as the (relative) sea level height observed at the $i$th tide gauge; their difference is indicated by $u_i$. As we are interested in sea level rates, we consider time derivatives of these quantities, indicated by dot

$${\dot{u}}_i={\dot{g}}_i-{\dot{s}}_i$$

(1)

This equation is sufficient to obtain estimates of the VLM rates at each tide gauge, as long as we ensure time consistency for the rates involved in the equation; that is, both ${\dot{g}}_i, {\dot{s}}_i$ refer to the same period of time and have negligible inherent drift and error. An average estimate of the VLM rates among the chosen TG locations can be obtained by using a least square procedure applied to Equation (1) [7]. However, as the timespan of altimetric observations roughly covers two decades, the formal error of VLM rates is heavily dependent on seasonal, interannual, decadal, and interdecadal signals. Long time series of historic tide gauges can supply low uncertainty in estimating relative vertical motion (RVM) between adjacent tide gauges. Kuo et al. [7,8] used an extension of the Gauss–Markov model with stochastic constraints [51] to fully exploit the low uncertainty introduced by RVM estimates. A similar approach was followed in [9] by resorting to the linear inverse problem with constraints (LIP) formulation in the framework of discrete inverse theory [28] to analyze the non-climate contribution to sea level rise in southern Europe [9]. In both approaches, a linear model is used to optimally and simultaneously interpolate VLM estimates of a network of tide gauges and satellite altimetry observations, constraining the model with the RVM estimates between adjacent TGs forcing smaller uncertainties in the variance–covariance matrix of the formal errors.

In such approaches, the rate of relative vertical motion, ${\dot{ru}}_{ij}=\left({\dot{g}}_i-{\dot{s}}_i\right)-\left({\dot{g}}_j-{\dot{s}}_j\right)$, between two nearby tide gauges is considered, as it reduces to ${\dot{ru}}_{ij}= {\dot{s}}_j-{\dot{s}}_i$ when neighboring tide gauges experience the same absolute sea level (in such case, ${\dot{g}}_{ i }-{\dot{g}}_j=0$). Consequently, the RVM definition can be extended to the complete period for which the two neighboring TG observations overlap, regardless of the temporal coverage of satellite altimetry observations. To distinguish the relative sea level rates in longer periods of overlapping TG observations from those related to the shorter periods of overlapping TG–ALT, we use the Greek letter $\mathsf{\zeta }$, and the RVM is thus better defined by

$${\dot{ru}}_{ij}={\dot{\mathsf{\zeta }}}_j-{\dot{\mathsf{\zeta }}}_i$$

(2)

Equation (1) for a number (N) of TGs in matrix form is

$$G\dot{u}=d;\dot{u}=\left(\begin{array}{c}{\dot{u}}_1\\ ⋮\\ {\dot{u}}_{\mathrm{N}}\end{array}\right);d=\left(\begin{array}{c}{\dot{g}}_1-{\dot{s}}_1\\ ⋮\\ {\dot{g}}_N-{\dot{s}}_N\end{array}\right);G=I_N,$$

(3)

While the M ≤ N Equation (2) constraints to the linear system of Equation (3) are written as:

$$F\dot{u}=h;h=-F\left(\begin{array}{c}{\dot{\mathsf{\zeta }}}_1\\ ⋮\\ {\dot{\mathsf{\zeta }}}_N\end{array}\right);F=\left(\begin{array}{cccccc}1& -1& 0& 0& \dots & 0\\ 0& 1& -1& 0& \dots & 0\\ & \dots & & & \dots & \\ 0& 0& 0& \dots & 1& -1\end{array}\right).$$

(4)

The constraints can be chosen arbitrarily, but they have to be linearly independent so that the rank of matrix $F$ is ≤ N − 1 and the condition expressed in Equation (2) is true (i.e., ${\dot{g}}_i-{\dot{g}}_j=0$). A unique linear system composed by Equations (3) and (4) with the use of Lagrange multipliers [28] is written by minimizing the form $\mathsf{\Phi }\left(\dot{u}\right)=e^{T}e+2{\lambda }^T\left(F\dot{u}-h\right)$ obtained by the sum of the L₂ norm of the prediction error $e=G\dot{u}-d$ and the inner product of the constraint equation $F\dot{u}-h=0$ by the Lagrange multiplier ${\lambda }^T$. The resulting OLS equation in matrix form is

$$\left(\begin{array}{cc}{G}^{T}G& F^T\\ F& 0\end{array}\right)\left(\begin{array}{c}\dot{u}\\ \lambda \end{array}\right)=\left(\begin{array}{c}{G}^{T}d\\ h\end{array}\right)$$

(5)

where $d, h$ are known and $\dot{u}, \lambda$ are unknown (vertical motion rates and Lagrange multipliers, respectively). Such system is of the form $A\cdot X=Y$, with $A=\left(\begin{array}{cc}{G}^{T}G& F^T\\ F& 0\end{array}\right)$. Following [28], the solution for the model parameters is easily obtained by multiplying by the inverse of A

$$X=\left(\begin{array}{c}\dot{u}\\ \lambda \end{array}\right)=A^{-1}Y;Y=\left(\begin{array}{c}{G}^{T}d\\ h\end{array}\right)$$

(6)

The standard errors of $\dot{u}$ were estimated as the diagonal elements of the variance–covariance matrix of the OLS estimator, i.e., from the diagonal terms of the covariance matrix $A^{-1}\mathsf{\Omega } {\left(A^{-1}\right)}^T$, where Ω is given by

$$\mathsf{\Omega }=\left(\begin{array}{cccc}{\mathsf{\sigma }}_1^2& 0& \dots & 0\\ 0& {\mathsf{\sigma }}_2^2& \dots & 0\\ ⋮& & \ddots & ⋮\\ 0& 0& \dots & {\mathsf{\sigma }}_{N+L}^2\end{array}\right)$$

(7)

With N and L, respectively, as the number of parameters ${\dot{u}}_i$ and constraints ${r\dot{u}}_{ij}$. The previous expression for Ω holds assuming no autocorrelation and heteroscedasticity of the regression residuals. The resultant error estimators of the OLS estimators are generally referred to as heteroscedasticity-consistent standard errors or White–Huber robust standard errors [52]. Finally, the estimated parameters are given by $X_i\pm \delta X_i$, where

$$\begin{gathered} X=A^{-1}Y;X_i=\{\begin{array}{c}{\dot{u}}_i \mathrm{if} 0\le i\le N \\ {\mathsf{\lambda }}_i \mathrm{if} N+1\le i\le N+L\end{array} \\ \delta X_i=diag{\left(A^{-1}\mathsf{\Omega } {\left(A^{-1}\right)}^T\right)}_i. \end{gathered}$$

(8)

3.2. A Possible Solution to the Limitation of the LIP Method

Here, we take advantage of the opportunity to extend the analysis with the current temporal availability of satellite altimetry observations included in two state-of-the-art datasets, SLCCI (1993–2015) and C3S (1993–2018), on the narrow basin of the Adriatic Sea. Despite its limited spatial scale when compared to, for example, the Mediterranean Sea or wider ocean basins, there are point-by-point differences in the absolute sea level trends observed by satellite altimetry. Such differences impose a conceptual limit on the suitability of the LIP method, as different TGs do see different absolute sea level rates, thus breaking the requirement that the relative VLM rate experienced by two neighboring TGs is exactly due to the differences of the two relative sea level trend change rates, as stated in Equation (2). In this study, we tried to overcome this limitation by removing from the TG time series and the altimetry time series associated with each TG a linear trend equal to that measured by the altimeter. In other words, we performed a change of variables (COV), which does not alter the statistical properties of the TG and altimetry time series, but eliminates any differences in relative sea level changes due to different absolute sea level changes

$$\begin{gathered} {\dot{g}}_i\to {\dot{g}}_i^{\prime }={\dot{g}}_i-{\dot{g}}_i=0; \\ {\dot{s}}_i\to {\dot{s}}_i^{\prime }={\dot{s}}_i-{\dot{g}}_i; \\ {\dot{\zeta }}_i\to {\dot{\zeta }}_i^{\prime }={\dot{\zeta }}_i–{\dot{g}}_i; \\ {\dot{u}}_i={\dot{g}}_i-{\dot{s}}_i={\dot{g}}_i-{\dot{g}}_i+{\dot{g}}_i-{\dot{s}}_i={\dot{g}}_i^{\prime }-{\dot{s}}_i^{\prime }. \end{gathered}$$

(9)

In this change of variable, two assumptions are made:

(1)

All the time series have a linear trend in every period in which they are considered;

(2)

The absolute sea level rates observed by satellite altimetry in its era can be extended backward in time to cover the timespan of the associated TG relative sea level time series.

The first assumption is easily verified, as we considered only time series with a clear and visible linear trend. The second assumption is needed to permit the third change of variable in Equation (9). Even if it looks restrictive, it is not, as we formed the sea level time series of TGs only from 1974 onward, up to the end of the satellite altimetry dataset considered (2015 for SLCCI, 2018 for C3S).

In this work, all TG sea level trend errors are calculated considering serial correlation and are given with a 95% confidence interval. The altimeter grid point associated with the TG was chosen based on the shortest distance from the TG and the maximum correlation with the TG monthly time series. Another experiment was attempted, in which the altimetry time series associated with the TGs were formed by averaging all the altimeter grid points in a distance range of 10–50 km and farther than 10 km from the coast, to avoid land contamination. No major differences were found in this case, and the results were omitted. For both altimetry and TG time series, annual and semi-annual cyclic variations were subtracted by fitting the monthly time series with a sinusoid made up of appropriate harmonics (annual and semi-annual), whose parameters, phase and amplitude were determined by an OLS procedure [53].

4. Results

4.1. C3S and SLCCI Absolute Sea Level Change Rates 1993–2015

The two altimetric datasets span different periods of time, as SLCCI data is confined to the period 1993–2015, while C3S is continuously updated, and the period covered in this study extends from 1993 to 2018. They differ for the number of contemporaneous satellites used to determine the SLA monthly values: all satellites for SLCCI, two satellites for C3S. Even the regionalization of the products is different, as SLCCI has been released as a unique global dataset at 0.25 × 0.25 km of resolution, while the resolution of the Mediterranean Sea regional dataset of C3S used in this study is 0.125 × 0.125 km. It is therefore natural to ask whether the two datasets can reproduce the same behavior when limited to a common time span. To do so we calculated the sea level change and vertical land motion rates for both datasets in the overlapping period 1993–2015, at the six TG locations used in this study, associating the altimetry grid points and to the TGs with the same approach already described above. The result of the comparison between absolute sea level change rates from altimetry of the two datasets is depicted in Figure 7. Overall, the two datasets present very similar pooled means and standard deviations over the six TG locations: 4.23 ± 1.69 mm y⁻¹ for C3S and 4.14 ± 1.80 mm y⁻¹ for SLCCI. Here the pooled mean is the same as described in [50], while the pooled standard deviation is calculated according to [54]. The pooled standard deviation estimate assumes that a combination of several series of measurements are performed under similar conditions, and an improved estimate of the imprecision of the process is desired. If it can be assumed that all the series are of the same precision although their means may differ, the pooled standard deviations from k series of measurements can be calculated as: ${\sigma }_p=\sqrt{[{\Sigma }_i\left(n_{i-1}\right){\sigma }_i{}^2]/{\Sigma }_i\left(n_i-1\right)}$, where σ_i are the standard deviation of each series of measurements, and n_i the number of samples in each series. There is good agreement on each individual TG value among the two datasets, except for Trieste: in the C3S dataset the TRIESTE SL rate is aligned with those of the other TGs, while in the SLCCI it is lower than the others (see Figure 7). This situation, although possible, does not seem plausible, and we are inclined to believe that the TRIESTE situation is not consistently described by the SLCCI dataset. A possible reason for the different description of the SL change near Trieste by the two datasets is the optimization for the Mediterranean Sea applied to C3S, as well as its finer spatial resolution: both arguments justify the more appropriated trend in Trieste given by C3S with respect to SLCCI.

4.2. Vertical Land Motion Rates Derived with Standard (ALT-TG) and LIP Approaches 1993–2015

During the period of overlap between the two altimetry datasets (1993–2015), the picture emerging from the estimates of the VLM rates at the six TGs along the northern and eastern Adriatic coast by the standard approach ($\dot{g}-\dot{s}$) does not show striking differences, apart from the situation of TRIESTE, already analyzed in the previous section.

Table 3. Results of calculations using C3S altimetry dataset (1993–2015). Column 1: location; columns 2/3: absolute/relative sea level rates in the altimetry era; column 4: relative sea level rates during TG lifetime; column 5: VLM calculated with classical approach (ALT-TG); columns 6–8: VLM derived with LIP and LIP_cov approaches and observed by CGPS. All data are in mm y⁻¹.

Location	$\dot{\mathit{g}}$	$\dot{\mathit{s}}$	$\dot{\mathit{\zeta }}$	$\dot{\mathit{g}}-\dot{\mathit{s}}{}^1$	LIP	LIPcov 2	GPS 3
VENEZIA	4.16 ± 1.71	6.08 ± 2.08	3.29 ± 0.82	−1.93 ± 0.79	−0.94 ± 0.42	−1.01 ± 0.42	−1.59 ± 0.65
VEPTF	4.17 ± 1.73	6.44 ± 2.07	3.83 ± 0.81	−2.27 ± 0.81	−1.46 ± 0.48	−1.52 ± 0.48	-
TRIESTE	4.51 ± 1.84	4.49 ± 1.98	2.40 ± 0.74	0.02 ± 0.71	−0.05 ± 0.37	0.23 ± 0.37	−0.25 ± 0.52
ROVINJ	4.09 ± 1.85	1.91 ± 2.04	1.61 ± 0.73	2.18 ± 1.17	0.73 ± 0.40	0.60 ± 0.40	-
SPLIT	4.44 ± 1.57	4.15 ± 1.87	2.45 ± 0.72	0.29 ± 0.70	−0.10 ± 0.37	0.11 ± 0.37	0.10 ± 0.64
DUBROVNIK	4.01 ± 1.42	4.67 ± 1.70	2.80 ± 0.63	−0.66 ± 0.62	−0.46 ± 0.51	−0.69 ± 0.51	-
Pooled mean	4.23 ± 1.69	4.62 ± 1.96	2.73 ± 0.75	−0.40 ± 0.82	−0.38 ± 0.43	−0.38 ± 0.43	−0.58 ± 0.60

¹$\dot{g}-\dot{s}$ is calculated as the linear trend of the differenced time series $g\left(t\right)-s\left(t\right)$. ² Subscript “cov” indicates results obtained after the change of variables as described in Section 3.2. ³ GPS values reported here and in the following table are derived for time spans given in Table 2.

and

Table 4. Results of calculations using SLCCI altimetry dataset (1993–2015). Same columns as in Table 3. All data are in mm y⁻¹.

Location	$\dot{\mathit{g}}$	$\dot{\mathit{s}}$	$\dot{\mathit{\zeta }}$	$\dot{\mathit{g}}-\dot{\mathit{s}}$	LIP	LIPcov	GPS
VENEZIA	4.47 ± 2.07	6.08 ± 2.08	3.29 ± 0.82	−1.61 ± 0.91	−1.02 ± 0.41	−0.70 ± 0.41	−1.59 ± 0.65
VEPTF	4.47 ± 2.07	6.44 ± 2.07	3.83 ± 0.81	−1.97 ± 0.93	−1.54 ± 0.47	−1.21 ± 0.47	−
TRIESTE	3.42 ± 1.78	4.49 ± 1.98	2.40 ± 0.74	−1.07 ± 0.74	−0.14 ± 0.36	−0.86 ± 0.36	−0.25 ± 0.52
ROVINJ	4.37 ± 1.86	1.91 ± 2.04	1.61 ± 0.73	2.46 ± 1.04	0.65 ± 0.39	0.87 ± 0.39	−
SPLIT	4.15 ± 1.48	4.15 ± 1.87	2.45 ± 0.72	−0.00 ± 0.76	−0.18 ± 0.36	−0.17 ± 0.36	0.10 ± 0.64
DUBROVNIK	3.98 ± 1.45	4.67 ± 1.70	2.80 ± 0.63	−0.69 ± 0.70	−0.55 ± 0.50	−0.71 ± 0.50	−
Pooled mean	4.14 ± 1.80	4.62 ± 1.96	2.73 ± 0.75	−0.48 ± 0.85	−0.46 ± 0.42	−0.46 ± 0.42	−0.58 ± 0.60

report the values of the calculated parameter for C3S and SLCCI, respectively. VENEZIA, VEPTF, ROVINJ, SPLIT, and DUBROVNIK VLM rates are in good agreement across the two datasets, with differences not exceeding 0.3 mm y⁻¹ and a mean standard error of 0.82 mm y⁻¹. In the LIP approach, the differences between the pair of VLM values of each TG are even lower, not exceeding 0.1 mm y⁻¹, with a mean standard error of 0.43 mm y⁻¹, roughly half the values of the estimate by the standard approach. The LIP method is able to handle the different situations depicted for TRIESTE by the two altimetry datasets and provides a very similar solution for both: −0.05 ± 0.37 mm y⁻¹ for C3S and −0.14 ± 0.36 mm y⁻¹ for SLCCI. Finally, using the LIP approach with COV (LIP_cov), the differences between the pair of VLM values of each TG do not exceed 0.3 mm y⁻¹, apart from the TRIESTE case, where the difference between estimates rises to 1 mm y⁻¹ and mean standard error is 0.43 mm y⁻¹ again.

The ability of the LIP approach to manage data errors and find the best solution to the linear problem can be seen in Figure 8, where the VLM estimates for 1993–2015, obtained with the three approaches, are reported by mapping SLCCI against C3S. The LIP scatterplot is indicated by yellow crosses. They are arranged along a straight line parallel to the perfect fit line (in blue), at a very short distance from it, expressing the existence of a constant bias between the estimates obtained from the two datasets. So, the LIP procedure minimized and equally distributed the errors of all estimates. The estimates derived with the standard approach (orange crosses) are instead visibly scattered around the line of perfect fit and span a more extreme range of values. The estimates of the LIP_cov approach are in the middle of these two situations, as they are more scattered around the perfect fit line than those obtained in the LIP approach, but span a shorter range of values than the standard approach. In practice, changing the variable in the LIP_cov approach reinserts the inaccuracies of the absolute sea level change rate estimates in the system, optimizing the errors only for the tide gauge estimates.

Figure 9 shows the plot of the VLM estimates for the six TGs as calculated by the three approaches and two altimetry datasets. In the standard approach, as there is no optimization of errors, we again see a wider range of values reached in both datasets. This is particularly evident for ROVINJ TG, whose $(\dot{g}-\dot{s})$ estimates reach more than 2 mm y⁻¹, while the LIP approach calculates it as less than 1 mm y⁻¹. On the other hand, if the standard approach returns a wider range of values, those calculated by the LIP and the LIP_cov are narrow and concentrated along a common path, except for the TRIESTE location, for which the LIP_cov estimate for the SLCCI dataset is affected by the overestimated rate of absolute sea level change. Figure 9 also shows the GPS estimates for the locations of VENEZIA, TRIESTE, and SPLIT.

We can assess how far the VLM estimates are from the GPS observations for the stations of VENEZIA, TRIESTE, and SPLIT. To do so, we use a simple metric, the root mean square difference (RMSD) of the two time series. The results, reported in

Table 5. Root mean square deviation of GPS observations and VLM estimates (1993–2018). Column 1: dataset name; columns 2–4: difference of VLM estimates from GPS and classical approach, from GPS and LIP, and from GPS and LIP_cov, respectively. All data are in mm y⁻¹.

Dataset	$\mathbf{GPS}-\left(\dot{\mathit{g}}-\dot{\mathit{s}}\right)$	GPS-LIP	GPS-LIPcov
C3S	0.274	0.409	0.435
SLCCI	0.477	0.372	0.642

, are obviously partial for two reasons. First, we can calculate the RMSD only for three of the six locations. Second, the GPS measurements started after the altimetry era began and cover different time periods: 2010–2020 for VENEZIA, 2003–2020 for TRIESTE, and 2004–2012 for SPLIT. Nonetheless, we believe that the GPS estimates can be representative of a wider period because they are caused by geological processes that have longer time scales than sea level variability in the altimetry era. Venice and Trieste belong to a zone of predominant Holocene subsidence, while Split is located in a zone of vertical stability [55].

In Table 5, we see that the RMSD values of the LIP approaches from the GPS estimates are around 0.4 mm y⁻¹ for the C3S dataset, not much higher than that of the standard approach, whose standard error is almost double that of the LIP-derived estimates. For the SLCCI estimates, the RMSD values of the standard and LIP_cov approaches are larger than those of the C3S dataset, for the underestimated rate of absolute sea level change in TRIESTE.

4.3. C3S Rates 1993–2018

When the longest time span of the C3S altimetric datasets is used (1993–2018 in this study), the contribution of the last three years dramatically changes the rates of absolute sea level change at the six TG location in the Adriatic Sea; for all of them the rates are lower by approximately 0.8 mm y⁻¹ than those calculated for 1993–2015, as seen in Figure 10.

However, the relevant modification of absolute sea level change rates, and consequently of relative sea level change rates, observed at the TGs only partially impacts the VLM rates derived with the three methods (see

Table 6. Results of calculations using C3S altimetry dataset (1993–2018). Columns are the same as in Table 3. All data are in mm y⁻¹.

Location	$\dot{\mathit{g}}$	$\dot{\mathit{s}}$	$\dot{\mathit{\zeta }}$	$\dot{\mathit{g}}-\dot{\mathit{s}}$	LIP	LIPcov	GPS
VENEZIA	3.36 ± 1.45	5.15 ± 1.73	3.26 ± 0.73	−1.79 ± 0.65	−0.83 ± 0.39	−0.93 ± 0.39	−1.59 ± 0.65
VEPTF	3.38 ± 1.46	5.50 ± 1.73	3.78 ± 0.73	−2.12 ± 0.67	−1.33 ± 0.47	−1.41 ± 0.47	−
TRIESTE	3.75 ± 1.58	3.56 ± 1.66	2.30 ± 0.67	0.18 ± 0.60	0.13 ± 0.33	0.42 ± 0.33	−0.25 ± 0.52
ROVINJ	3.33 ± 1.58	1.03 ± 1.85	1.36 ± 0.71	2.30 ± 1.06	1.06 ± 0.37	0.93 ± 0.37	−
SPLIT	3.60 ± 1.36	2.92 ± 1.65	2.20 ± 0.66	0.68 ± 0.63	0.23 ± 0.33	0.37 ± 0.33	0.10 ± 0.64
DUBROVNIK	3.34 ± 1.22	3.79 ± 1.48	2.69 ± 0.58	−0.45 ± 0.55	−0.29 ± 0.46	−0.41 ± 0.46	−
Pooled mean	3.46 ± 1.45	3.66 ± 1.69	2.60 ± 0.68	−0.20 ± 0.71	−0.17 ± 0.40	−0.17 ± 0.40	−0.58 ± 0.60

), as the VLM estimates change on average by about 0.2 mm y⁻¹ (minimum 0.1 and maximum 0.4 mm y⁻¹ at Venice and Split, respectively). All results are reported in

Table 7. Comparison of C3S estimated VLMs of 1993–2018 and 1993–2015. Column 1: location; columns 2, 5, 8: VLM estimated with classical, LIP, LIP_cov approach in 1993–2018; columns 3, 6, 9: same for 1993–2015; columns 4, 7, 10: differences between VLM estimates in the two periods. All data are in mm y⁻¹.

	$\dot{\mathit{g}}-\dot{\mathit{s}}$			LIP			LIPcov
Location	2018 1	2015	diff	2018	2015	diff	2018	2015	diff
VENEZIA	−1.79	−1.93	0.14	−0.83	−0.94	0.11	−0.93	−1.01	0.08
VEPTF	−2.12	−2.27	0.15	−1.33	−1.46	0.13	−1.41	−1.52	0.11
TRIESTE	0.18	0.02	0.16	0.13	−0.05	0.18	0.42	0.23	0.19
ROVINJ	2.30	2.18	0.12	1.06	0.73	0.33	0.93	0.60	0.33
SPLIT	0.68	0.29	0.39	0.23	−0.10	0.33	0.37	0.11	0.26
DUBROVNIK	−0.45	−0.66	0.21	−0.29	−0.46	0.17	−0.41	−0.69	0.28

¹ Here 2018 and 2015 indicate whether the C3S altimetry dataset was used in full (1993–2018) or limited (1993–2015) extent. Maximum and minimum value of the differences are obtained for SPLIT (fourth column) and VENEZIA (tenth column) respectively.

.

The different estimates of VLM rates for the C3S altimetry dataset are plotted in Figure 11 for the sake of clarity. Figure 11 also shows the results obtained by W&M in 2012 [9] with the LIP approach and a different altimetry dataset for the period 1992–2010.

Apart from the much lower standard errors of the W&M solution, which we presume were due to the different methodology in obtaining the formal errors of absolute and relative sea level change rates, it appears evident that in the years following 2010, the VLM rates at the five TG locations common to both studies (W&M and the present study) have remained almost steady. The little discrepancies observed between the results of the two studies can be ascribed mainly to the different periods covered by the altimetry datasets (C3S is 44% longer and SLCCI is 28% longer compared to the W&M study). Another important factor is the processing of the altimetry data: in our study, we used more recent gridding than the original W&M dataset, with more accurate corrections and, in the case of C3S, finer resolution and regional tailoring for the Mediterranean Sea. The rates of absolute sea level change at the TGs, calculated as the sum of relative sea level change and VLMs derived in this study for the whole record of the corresponding TGs, are reported in

Table 8. Rates of absolute sea level change from TG records over whole record (1974–2018), corrected for VLM estimated with the LIP_cov approach. All data are in mm y⁻¹.

Location	$\dot{\mathit{\zeta }}-\mathbf{L}\mathbf{I}{\mathbf{P}}_{\mathbf{c}\mathbf{o}\mathbf{v}}$
VENEZIA *	2.33 ± 0.83
VEPTF	2.37 ± 0.86
TRIESTE *	2.71 ± 0.75
ROVINJ	2.29 ± 0.80
SPLIT *	2.57 ± 0.74
DUBROVNIK	2.28 ± 0.74
Pooled mean	2.43 ± 0.80
Sample mean	2.43 ± 0.18

* Asterisks near the TG name indicates that the SL record has been reconstructed from different sources.

.

The uncertainty of the sample mean (last row of Table 8) was obtained by propagating the formal errors, considering the rates as random and independent variables. The absolute sea level change rates vary in a very narrow interval, 2.33–2.71, with a sample mean of 2.43 mm y⁻¹. The standard deviation of the sample is 0.18 mm y⁻¹, which is much lower than the precision of each individual determination of SL change rate at the TGs. W&M, in their study, obtained a slightly lower dispersion for the Adriatic region mean rate: 1.36 ± 0.13 mm y⁻¹ [9]. We agree with W&M, who pointed out that such a low dispersion is unlikely to be determined from estimates of independent random variables and underlined that it is an indication of the high performance of Kuo et al.’s method [7] for determining accurate VLM rates from TG and altimetry differenced time series. On a global scale, the determination of absolute change in SL increases as the upper end of the considered period is more recent. Church et al. [56] calculated the global mean sea level rise from observations of its principal contribution (thermal expansion, continental and Greenland glaciers, land water storage) at a rate of 2.0 ± 0.3 mm y⁻¹ in the period 1971–2010. On the other hand, Hay et al. [57] estimated a global mean sea level of 3.0 ± 0.7 mm y⁻¹, while Church and White [58] reported 2.8 ± 0.5 mm y⁻¹ during 1993–2010. The global estimates supplied by satellite altimetry in the last few years have kept steady around 3.1–3.3 mm y⁻¹ (Cazenave and Le Cozannet [59], Cazenave et al. [19]), indicating an acceleration of the global SL estimate [19]. Similar behavior is observed at a regional scale in the Adriatic Sea:

• 1.36 ± 0.13 mm y⁻¹ (19??–2010), W&M [9] (question marks indicate that the lower end of the period is different for each TG and coincident with the beginning of the observation period);
• 2.27 ± 0.41 mm y⁻¹ (1974–2015), SLCCI, this study (not shown);
• 2.35 ± 0.20 mm y⁻¹ (1974–2015), C3S, this study (not shown);
• 2.43 ± 0.18 mm y⁻¹ (1974–2018), C3S, this study.

Considering the different time periods covered, the results of our study support, though underrate, the estimates given at the global scale, but reinforce the hypothesis that an acceleration in global sea level rate has been occurring in the last two decades. Small differences may exist between local and global representations of the same global process; we consider it reasonable that the SL trend observed in the Adriatic Sea, a semi-enclosed basin with peculiar atmospheric and hydrodynamic conditions [14,60], results in a slight disagreement with the global trend. At the Mediterranean scale, we still observe both a slower regional absolute sea level rate and a tendency toward acceleration:

• 0.7 ± 0.2 mm y⁻¹ (1945–2000) [61];
• 1.60 ± 0.35 mm y⁻¹ (19??–2010), W&M [9];
• 2.44 ± 0.5 mm y⁻¹ (1993–2012) [22];
• 2.87 ± 0.33 mm y⁻¹ (1992–2016) [62].

Some authors attribute the difference between the global and regional absolute sea level change rates of the Mediterranean region to atmospheric and steric processes in the last few decades, specifically the displacement of large-scale atmospheric patterns over the basin, which brought higher atmospheric pressure over the region, progressive water cooling, and increased salt concentration [63].

5. Discussion and Conclusions

We estimated relative and absolute sea level trends as well as vertical land motion rates, with their errors, at six locations in the Adriatic Sea; we used data produced by three measuring systems (tide gauges, radar altimetry, GPS), integrating the information from those systems in order to maximize the knowledge, both qualitatively and quantitatively. We assessed two altimetry products (ESA SLCCI and Copernicus C3S) specifically processed for climate studies. We compared the results with the direct method (subtracting relative from absolute sea level rates) and with linear regression, considering the system as described by the constrained linear inverse problem. Such an approach allowed us to simultaneously solve for the rates of all TGs. We also tested the robustness of the constrained linear inverse problem with respect to the limitation requiring that absolute sea level rates observed by each pair of TGs forming a constraint to the inverse linear problem are the same. We found that the two altimetry products, SLCCI and C3S, supply similar results, but there are differences, mainly in the rates calculated for Trieste, which present a suspicious absolute sea level rate in the SLCCI dataset. Errors in the calculated VLM rates are slightly lower for the C3S dataset, which covers a period 13% longer than SLCCI.

The simultaneous solution for the VLM rates with the constrained linear inverse problem had the best performance as the standard error estimates, in both the standard LIP and LIP_cov linear inverse problem with COV formulation. The errors in the VLM rates are on the order of 0.3–0.4 mm y⁻¹. Overall, we obtained a consistent representation of the absolute and relative sea level change rates on the northern and eastern coast of the Adriatic Sea from altimetry and tide gauges, but with small differences explained by the local VLM, acting at an average rate of −0.38 ± 0.43 mm y⁻¹, but with different rates and signs from place to place, in the range [+0.60, −1.52]. These rates are confirmed by the VLM rates derived from the CGPS stations of Venice, Trieste, and Split. The LIP approach revealed a good ability to handle errors in the data, as it shares errors among all the components, thus minimizing their impact on the final estimates.

We propose a new methodological use of the linear inverse problem with constraints, which suffers from a limitation in its applicability, as it requires that every TG pair used as constraint is affected by the same absolute sea level; such limitation prevents the method from being applicable to wider basins. In our formulation, a simple change in variable in the TG relative SL and altimetry absolute SL time series allows the method to be applicable to every TG pair, whichever local absolute SL is seen by the two TGs. The LIP_cov approach extends the applicability of the LIP to more complex situations in which the absolute sea level rates are not similar among neighboring TGs. However, its sensitivity to inaccuracies in the absolute sea level rates should be assessed in more detail. One could question what benefit the introduction of the COV in the LIP approach brings. The answer is not easy, as the COV was introduced in order to cancel the dependence of the TG relative SL rates at the TGs from the absolute SL influencing them. In any case, the results of the COV in the LIP approach provide a solution that is almost indistinguishable from the corresponding results of LIP without changing variables. Moreover, the COV should grant the validity of the approach irrespective of the choice of TGs paired in the constraints.

Using the VLM rates calculated with this approach, we estimated absolute SL rates at the six TGs by adding VLM rates derived from the C3S dataset to the observed relative SL rate (1974–2018). We found that the absolute SL rates had a sample mean of 2.43 mm y⁻¹ and sample standard deviation of 0.18 mm y⁻¹. The mean absolute SL increase estimated from the SLCCI dataset, valid for the period 1974–2015, is 2.27 ± 0.41 mm y⁻¹, not far from the previous estimate but affected by a wider dispersion.

The very low dispersion of the six estimates around their mean in the C3S case confirms the accuracy of the LIP method in determining VLM rates, and suggests that the change of variable of the LIP allows it to be applicable to wider basins and more heterogeneous collections of TG time series, without requiring a priori splitting of the basin in homogeneous zones of absolute SL change, as in [9].

We stress that the different results of the two altimetry datasets, C3S and SLCCI, are probably due to the coarser resolution and the lack of a regionalization of the second set, as well as its shorter time coverage. Indeed, the SLCCI and C3S products are generated from different processing chains, have different spatial and temporal resolutions, and have different numbers of satellite missions, and C3S relies on a mapping algorithm specifically dedicated to the Mediterranean Sea.

On the other hand, the GPS data used in this study span very different time periods, always shorter than altimetry and TGs time series. Moreover, it is difficult to understand whether CGPS stations effectively reflect the tide gauge movement, as they are often separated from the TG basement and continuous and accurate levelling of their benchmarks is not felt everywhere with the same urgency.

While altimeter data are available globally, tide gauge observations are bound to selected places, but not all of them are equipped with GPS. This study highlights the importance of land movement for sea level change. The synergistic use of altimetry and tide gauges together, in a robust inverse linear approach, provides a surrogate for VLM in the absence of GPS measurements, giving insights for understanding how future sea level scenarios will impact the coast at local scale.