The Sources of Sea-Level Changes in the Mediterranean Sea Since 1960

Past sea-level changes in the Mediterranean Sea are highly non-uniform and can deviate significantly from both the global average sea-level rise and changes in the nearby Atlantic. Understanding the causes of this spatial non-uniformity is crucial to the success of coastal adaptation strategies. This, however, remains a challenge owing to the lack of long sea-level records in the Mediterranean. Previous studies have addressed this challenge by reconstructing past sea levels through objective analysis of sea-level observations. Such reconstructions have enabled significant progress toward quantifying sea-level changes, however, they have difficulty capturing long-term changes and provide little insight into the causes of the changes. Here, we combine data from tide gauges and altimetry with sea-level fingerprints of contemporary land-mass changes using spatial Bayesian methods to estimate the sources of sea-level changes in the Mediterranean Sea since 1960. We find that, between 1960 and 1989, sea level in the Mediterranean fell at an average rate of −0.3 ± 0.5 mm yr −1 , due to an increase in atmospheric pressure over the basin and opposing sterodynamic and land-mass contributions. After 1989, Mediterranean sea level started accelerating rapidly, driven by both sterodynamic changes and land-ice loss, reaching an average rate


Introduction
Global mean sea level has been rising more than two times faster in the past three decades (∼3.3 mm yr −1 ) than throughout most of the twentieth century (Dangendorf et al., 2017;Frederikse et al., 2020;Hay et al., 2015). This acceleration has been largely driven by increased rates of both land-ice melting and ocean thermal expansion (Dangendorf et al., 2019;Frederikse et al., 2020). Both of these processes give rise to highly non-uniform patterns of sea-level change, implying that the local rate of sea-level change can deviate significantly from the global average rate. In the case of land-ice melting, the associated patterns of (static-equilibrium) sea-level change are unique to each melt source (Clark & Lingle, 1977) and are known as sea-level fingerprints. This spatial non-uniformity in the sources of global sea-level rise highlights the need for accurate knowledge of how sea level changes at a local and regional scale.
Abstract Past sea-level changes in the Mediterranean Sea are highly non-uniform and can deviate significantly from both the global average sea-level rise and changes in the nearby Atlantic. Understanding the causes of this spatial non-uniformity is crucial to the success of coastal adaptation strategies. This, however, remains a challenge owing to the lack of long sea-level records in the Mediterranean. Previous studies have addressed this challenge by reconstructing past sea levels through objective analysis of sea-level observations. Such reconstructions have enabled significant progress toward quantifying sea-level changes, however, they have difficulty capturing long-term changes and provide little insight into the causes of the changes. Here, we combine data from tide gauges and altimetry with sea-level fingerprints of contemporary land-mass changes using spatial Bayesian methods to estimate the sources of sea-level changes in the Mediterranean Sea since 1960. We find that, between 1960 and 1989, sea level in the Mediterranean fell at an average rate of −0.3 ± 0.5 mm yr −1 , due to an increase in atmospheric pressure over the basin and opposing sterodynamic and land-mass contributions. After 1989, Mediterranean sea level started accelerating rapidly, driven by both sterodynamic changes and land-ice loss, reaching an average rate of 3.6 ± 0.3 mm yr −1 in the period 2000-2018. The rate of sea-level rise shows considerable spatial variation in the Mediterranean Sea, primarily reflecting changes in the large-scale circulation of the basin. Since 2000, sea level has been rising faster in the Adriatic, Aegean, and Levantine Seas than anywhere else in the Mediterranean Sea.
Plain Language Summary The rate of sea-level rise varies considerably across the Mediterranean basin. Characterizing this geographic variability is essential for the design of coastal adaptation plans. This, however, is no easy feat because of the sparseness of the observational sea-level record. In the Mediterranean Sea, long tide gauge records are restricted to a few locations along the European coastlines. Satellite altimetry provides sea-level observations with better spatial coverage but only since 1992. Past studies have addressed the difficulties related to data sparseness by combining data from tide gauges and altimetry to produce sea-level reconstructions. However, existing sea-level reconstructions have difficulty capturing long-term changes and do not distinguish between the causes of the changes. Here, we analyze sea-level observations in combination with fingerprints of land-ice melting using probabilistic methods to show that the average rate of sea-level rise in the Mediterranean Sea increased from −0.3 mm yr −1 in the period 1960-1989 to 3.6 mm yr −1 in 2000-2018, as a result of increased ocean warming and land-ice melting. Furthermore, we find that, since 2000, sea-level in the Adriatic, Aegean, and Levantine Seas has been rising faster than anywhere else in the Mediterranean Sea, due to changes in the circulation of the basin.

of 13
The Mediterranean Sea is a perfect illustration of why regional assessments of sea-level rise are crucial to coastal adaptation planning. For one thing, the Mediterranean is among the world's most vulnerable regions to climate change, with 86% of the region's World Heritage sites already at risk from coastal flooding and erosion (Reimann et al., 2018). For another, past sea-level changes in the Mediterranean Sea have been found to exhibit considerable spatial variation as well as to deviate significantly from both the nearby Atlantic and the global average (Calafat & Gomis, 2009;Marcos & Tsimplis, 2008;Orlić et al., 2019;Tsimplis & Baker, 2000). In particular, there is observational evidence that, between 1960 and 2000, sea level in the Mediterranean remained flat or slightly fell (depending on the study), in contrast to the Eastern North Atlantic where sea level rose at a rate of more than 1 mm yr −1 (Marcos & Tsimplis, 2008;Tsimplis et al., 2005;Tsimplis & Baker, 2000). The small rates of sea-level change in the Mediterranean over 1960-2000 were partly due to a persistent increase in atmospheric pressure over the basin related to the North Atlantic Oscillation (NAO) (Gomis et al., 2008;Tsimplis et al., 2005;Tsimplis & Josey, 2001). Since the early 1990s, the Mediterranean has seen a significant increase in the rate of sea-level rise, leading to an average rate of 2.7 mm yr −1 for the period 1993-2017 (Mohamed et al., 2019), which is on a par with the global average rate. The relative contributions from sterodynamic changes (i.e., changes in ocean density and circulation) and land-ice melting to this recent increase in the rate of Mediterranean sea-level rise remain unclear. Characterizing these spatiotemporal sea-level variations in the Mediterranean Sea and understanding their causes is a priority, in order to enable more effective coastal adaptation. However, achieving this goal is challenging because of the sparseness of the observational sea-level record in both time and space.
Tide gauge records are our main source of information on long-term sea level changes, but long records are restricted to a few coastal locations. In the Mediterranean Sea, the tide gauge network is heavily biased toward the European shorelines, with only one relatively long record (∼35 years) available on the African coasts. Satellite altimetry provides observations with better spatial coverage but only since 1992. Past studies have attempted to address the difficulties related to data sparseness by statistically combining tide gauge records with altimetry data to produce sea-level reconstructions. This is typically done by computing empirical orthogonal functions (EOFs) from the altimetry data (or sometimes from modeled data) and then fitting a subset of those to the tide gauge data to estimate the temporal amplitude associated with each retained EOF (Calafat et al., 2014;Calafat & Gomis, 2009;Calafat & Jordà, 2011;Meyssignac et al., 2011). This approach yields sea-level fields with the same spatial coverage as the altimetry data and spanning the same period as the tide gauge records. These fields provide valuable insights into the spatial structure of long-term sea-level changes and can also be used to compute spatially averaged changes. However, EOF-based reconstructions have difficulty capturing the long-term trends (Calafat et al., 2014), do not adequately account for differences between relative sea level (RSL, as observed by tide gauges) and geocentric sea level (GSL, as observed by altimetry), and do not have the ability to separate the causes of sea-level changes.
Here, we combine data from tide gauges and altimetry with sea-level fingerprints of contemporary land-mass changes using a Bayesian hierarchical model (BHM) (see Banerjee et al., 2004 for a general description of hierarchical models) to reconstruct sea-level changes in the Mediterranean Sea since 1960 and quantify the source contributions to such changes. Our approach has three key advantages over EOF-based reconstructions: (a) it explicitly distinguishes between RSL and GSL; (b) it separates long-term sea-level changes from short-term variability as well as large-scale signals from small-scale changes; and (c) it quantifies the individual contributions from sterodynamic changes, land-mass changes, and glacial isostatic adjustment (GIA, Tamisea, 2011) to long-term sea-level changes. The investigation of the new Bayesian solutions allows us to answer fundamental questions about sea level in the Mediterranean Sea, such as how much sea level in the basin has risen since 1960, whether the rate of change has been constant, and the origin of spatial variation in the sea-level changes.

Observational Data and Sea-Level Fingerprints
Monthly mean values of sea level from tide gauge records are obtained from the Revised Local Reference data archive of the Permanent Service for Mean Sea Level (PSMSL, https://www.psmsl.org/) (Holgate et al., 2013). Monthly values flagged as suspicious by the PSMSL are rejected. The inverse barometer (IB) effect is removed from the monthly time series based on sea-level pressure data from the NCEP/NCAR reanalysis (Kalnay et al., 1996). Annual and semi-annual cycles are also removed through ordinary least squares. The adjusted monthly data are then converted into annual time series for our analysis, noting that only years with at least 10 3 of 13 valid monthly values are considered. The final tide gauge data set comprises 2292 annual values from a total of 107 tide gauge records spanning the period 1960-2018. The location of the tide gauges is shown in Figure 1.
The satellite altimetry data are obtained from the multi-mission gridded sea surface height product provided by the Copernicus Marine Environment Monitoring Service (available at http://marine.copernicus.eu/). The data are provided as weekly fields on a 1/4° × 1/4° near global grid covering the period from January 1993 to December 2018. These weekly fields are averaged into annual fields for our analysis. The data are made available with all standard corrections applied, including corrections for tropospheric (wet and dry) and ionospheric path delays, sea state bias, tides (solid earth, ocean, loading, and pole), and atmospheric effects as quantified by the Mog2D barotropic model (Carrère & Lyard, 2003) for periods <20 days and the IB effect for longer periods.
The sea-level fingerprints associated with contemporary land-mass changes as well as the GIA estimates are the same as derived by Frederikse et al. (2020). The sea-level fingerprints represent the gravitational, rotational, and deformation (GRD) effects on sea level due to contemporary mass changes of the Greenland and Antarctic ice sheets, glaciers, and terrestrial water storage. The fingerprints are provided as an ensemble (5000 members) of annual sea-level fields for the period 1960-2018. Because GRD effects in our BHM are described in terms of rates rather than levels (to facilitate separation from the short-term variability), we convert the fingerprints from sea levels to yearly rates of change. This is done by first smoothing the sea-level time series at each grid point using an integrated-random-walk-plus-noise model that we fit through a Kalman smoother and then by taking first differences of the smoothed time series. This procedure removes any short-term variability that might exist in the original time series (see Figure S1 in Supporting Information S1 for an example). The GIA data are provided as an ensemble (5000 members) of time-invariant patterns of sea-level change rates. The ensembles allow us to derive probability distributions for the land-mass and GIA contributions, which we use as prior distributions in the BHM as described later. We note that both the sea-level fingerprints and the GIA patterns are different for RSL and GSL, and we take account of this when building the BHM.

BHM of Sea Level
Our goal here is to reconstruct sea-level changes in the Mediterranean Sea for the period 1960-2018 at any arbitrary location and quantify the individual source contributions to such changes. We achieve this by analyzing observations from tide gauges and satellite altimetry using a spatiotemporal probabilistic model. We adopt a Bayesian hierarchical approach, meaning that we use probability distributions to describe uncertainty in all model quantities (i.e., observations, source contributions, and unknown parameters). In designing the model, it is important to recognize that while tide gauges measure sea level relative to the land on which they reside (i.e., RSL), satellite altimetry measures sea level relative to the Earth's center of mass (i.e., GSL). This means that tide gauge measurements reflect changes in both sea surface height (SSH) and the underlying solid Earth (i.e., vertical land movements, VLMs), while altimeters only sense changes in SSH. Our model rigorously accommodates these differences between tide gauges and altimetry.
Another point worth noting is that the sea-level observations are subject to significant uncertainty arising from measurement error, uncertain altimetry corrections, etc. To account for this uncertainty, we assume that the observations are noisy measurements of a latent process which represents the true changes in sea level. The true latent sea-level process is then modeled as the sum of five contributions (or processes), namely sterodynamic changes ( ), contemporary GRD effects ( ), GIA ( ), small-scale changes ( ), and short-term variability ( ). The process captures long-term trends due to small-scale processes, including VLM from tectonics, sediment compaction, and groundwater. The model is specified as a BHM with three levels: (a) a probability model that describes the distribution of the sea-level observations conditional on the true latent process (data model); (b) a probability model that describes the dynamics of the five source contributions conditional on a set of parameters (process model); and (c) a prior distribution that describes the uncertainty in the model parameters and encodes our prior knowledge about the data and the processes (parameter model). Inferences in our model are made by sampling from the posterior distribution of the processes and parameters given the observations, which is proportional to the product of the three probability models that form the hierarchy. We note that the BHM assimilates sea-level data corrected for the IB effect, and thus it does not estimate the contribution from changes in atmospheric pressure. Such contribution is quantified separately based on data from the NCEP/NCAR reanalysis and assuming an IB effect. Next, we describe the three layers of the BHM.
, denote the set of measurement sites (tide gauges for = and altimetry for = ). Further, denote , ( ) as the sea level observation at time step t and location , for ∈ S and ∈ { } . Then, the data level of the BHM can be written as where , ( ) denote the true latent sea level (RSL or GSL), ( ) represent data offsets which are modeled as i.i.d. normal random variables ( ) ∼ 0, 2 , and the , ( ) are assumed to be serially and spatially uncorrelated observation errors , ( ) ∼ 0, 2 . Equation 1 is evaluated for = 1, . . . , when = but only for = 34, . . . , when = since the altimetry record starts in 1993 (this is also the case for any of the time-varying processes described below whenever these are subscripted by ). The altimetry data are subsampled at every five grid points in order to reduce the computational cost of the BHM.
Next, we describe the process level of the BHM. Except for the sea-level variability, all processes are described in terms of rates of sea-level change, not sea levels. This facilitates the separation of the variability from the long-term processes since it constrains the processes associated with long-term changes to be smooth. We assume that there is no difference between RSL and GSL in what refers to sterodynamic changes and short-term variability, whereas the contributions from GRD effects, GIA, and small-scale changes may be different for RSL and GSL. We recognize that sterodynamic changes can induce GRD effects through an influence on ocean bottom pressure, but such effects are generally small (Richter et al., 2013) and thus we ignore them.
The instantaneous rate of change due to sterodynamic changes ( ) is modeled as the sum of a spatiotemporal process that captures the effect of sterodynamic changes inside the Mediterranean ( ) plus a spatially uniform time-varying process that captures the effect of changes outside the basin ( ). Note that both processes ( and ) can affect the basin-average rate of change. The use of a spatially uniform process ( ) is motivated by the physical expectation that, on long time scales, any significant difference in sea level between the Mediterranean and the nearby Atlantic will tend to be rapidly compensated via changes in the Atlantic inflow through the Strait of Gibraltar (Calafat, Chambers, & Tsimplis, 2012). This compensation takes the form of a spatially uniform barotropic signal. The process is assumed to evolve according to a spatiotemporal random walk: where ( ) is a zero-mean Gaussian process (GP) (Rasmussen & Williams, 2006) ( ) ∼ GP (0, c ( , ′ ; , )) , with c(⋅, ⋅) being a covariance function and and denoting, respectively, the standard deviation and length scale that define the covariance function. The initial state of ( 0 ) is modeled as 0( ) ∼ GP 0, c , ′ ; 0 , 0 . For the covariance function of the GPs we assume a Matérn kernel with smoothing parameter = 5∕2 .
The spatially uniform process, , is assumed to follow a univariate random walk: where ∼ 0, 2 . The initial value, 0 , is modeled in the parameter layer of the BHM by placing a prior distribution on it (in practice, Equation 3 is evaluated using backward recursion, and thus 0 ∶= = ).
With this, the sterodynamic contribution to the rate of sea-level change is ( ) = ( ) + .
The rate of sea-level change due to contemporary GRD effects is modeled as: CALAFAT ET AL.
10.1029/2022JC019061 5 of 13 where ∈ { } , , ( ) is the ensemble-mean fingerprint of the rate of sea-level change associated with GRD effects from Frederikse et al. (2020) and ( ) is an error term that accounts for uncertainty in the fingerprint and is modeled as a spatiotemporal random walk. The initial value of is modeled as ,0( ) ∼ GP 0, Σ ,0 . The covariance matrices Σ and Σ 0 are extracted from the ensemble of RSL and GSL fingerprints (i.e., they capture variation across the ensemble members after removal of the ensemble mean).
The contribution from GIA is modeled as a stationary Gaussian process: where ∈ { } , and ( ) and Σ are, respectively, the mean and covariance matrix of the ensemble of GIA rates from Frederikse et al. (2020). The rate of sea-level change due to small-scale processes is modeled as i.i.d. normal random variables (i.e., time-invariant trends): where ∈ { } . Note that Equation 8 assumes the small-scale trends to be spatially uncorrelated. We tested another version of the BHM in which was allowed to vary in time as VLMs are not always well described by a linear trend, but such version of the model provided a worse fit to the data according to convergence diagnostics (e.g., larger R-hat values) and the solutions were very similar to those presented in this paper.
The sea-level rates due to sterodynamic changes, GRD effects, GIA, and small-scale processes can be integrated and combined to estimate the total long-term sea-level changes: where the initial value of is set to zero, 0( ) = 0 .
Finally, the short-term sea-level variability, , is modeled in terms of EOFs: , ∼ 0, 2 where is the temporal amplitude associated with the nth leading EOF. The EOFs are determined a priori from the full-resolution altimetry data after removal of a quadratic polynomial fit, whereas the temporal amplitudes are unknown and need to be estimated from the observations as part of the BHM. The first two EOFs alone capture 50.6% of the variance in the detrended altimetry data, while the first five EOFs explain 69.1% of the variance. Beyond five EOFs, the fraction of variance explained increases slowly with each subsequent EOF, indicating that such higher order EOFs may not be statistically significant. Indeed, past studies have found that there is no upside to including more than five EOFs in Mediterranean sea-level reconstructions (e.g., Calafat & Jordà, 2011). Based on this, we decide to use only the five leading EOFs. The spatial structure of the first three EOFs along with their temporal amplitudes are shown in Figure S2 (in Supporting Information S1).
We have tested the sensitivity of our results to the number of retained EOFs by comparing Bayesian solutions based on 2, 5, and 10 EOFs. We find that estimates of RSL rates are consistent across the three models within one posterior standard deviation (see Table S1 in Supporting Information S1). However, we also note that when using 10 EOFs convergence of the Markov chains (i.e., posterior samples) is worse than in the other two cases, indicating that including more than a few EOFs degrades the performance of the BHM. It is also important to recognize that using EOFs is not the only way to describe the variability. In particular, during the design phase of the BHM, we considered an alternative approach in which the variability was modeled as a spatial first-order autoregressive (AR1) process (similar to Piecuch et al., 2017). We found that when using EOFs the BHM captures a larger fraction of the variability and provides a more adequate fit to the data (according to convergence diagnostics) than when using an AR1 process.
With all the contributions now defined, the true latent sea-level containing all contributions can be written as: CALAFAT ET AL.
10.1029/2022JC019061 6 of 13 The BHM is completed by specifying the parameter model. Most of the unknown parameters are given prior distributions (i.e., they are modeled), while a few parameters are set to a fixed value based on prior knowledge or domain expertise. The prior distributions ascribed to the modeled parameters as well as the posterior mean and credible interval for such parameters are summarized in Table 1. The fixed parameters include , , , and . The standard deviation is set equal to 1 m, which implies a diffuse prior distribution on the data offsets ( ) (the time series from both tide gauges and altimetry have been processed in such a way that offsets are <1 m). Both and are set equal to 0.3 mm yr −2 , which encapsulates our prior belief that the sterodynamic rate is unlikely to change by significantly more than 0.3 mm yr −1 between consecutive years. Finally, the length scale is set equal to 5°, which is necessary to ensure that reflects large-scale changes, otherwise the model will favor smaller length scales owing to the relatively high resolution of the altimetry data for the period 1993-2018 (small length scales are not resolvable prior to 1993 due to the sparseness of the tide gauge data).
To make the problem computationally tractable, we approximate all GPs using nearest neighbor Gaussian processes with 20 neighbors, following the formulation given by Datta et al. (2016). The BHM is fitted using the No-U-Turn sampler (a variant of Hamiltonian Monte Carlo with adaptive optimization) as implemented by the Stan probabilistic programming language (Carpenter et al., 2017). We run the sampler with four chains of 1250 iterations each (warm-up = 500) for a total of 3000 post-warm-up draws.
The BHM provides estimates of sea-level changes, the yearly rate of sea-level change, and the contribution from each source for the period 1960-2018, both at tide gauge sites and on a 1/4° × 1/4° grid. Because our estimates take the form of samples from the joint posterior distribution of the processes at all locations and times, we are able to compute spatial and temporal averages of any quantities of interest with quantifiable and meaningful uncertainty. In the following, our estimates are summarized by posterior means and standard deviations (SDs).
As mentioned above, the contribution from atmospheric pressure changes is computed outside the BHM assuming an IB effect. The rates of change associated with such contribution are computed by fitting an integrated-random-walk-plus-noise model separately at each grid point via Stan.

Model Evaluation
We evaluate the BHM in the light of two questions: (a) does the sampler accurately characterize the posterior distribution?; and (b) does the BHM provide an adequate fit to the data? The first question is addressed through a number of Markov chain Monte Carlo (MCMC) diagnostics that aim to diagnose problems with the sampler and assess convergence and mixing. Such diagnostics include the potential scale reduction statistic (̂ ), the effective sample size per iteration ( ef f ∕ ), divergent transitions, and tree-depth saturation. We find that ̂ < 1.05 for all parameters in our model (̂ should be close to 1 at convergence), suggesting that the Markov chains have converged to the equilibrium distribution and are providing a good approximation to the posterior distribution.
We estimate ef f ∕ to be >0.1 for all scalar parameters, indicating low autocorrelation and good mixing (The Note. The potential scale reduction statistic (R-hat) and the effective sample size per iteration ( ef f ∕ ) are also shown. In general, R-hat should be close to 1 at convergence, whereas ef f ∕ 0.003 indicates low autocorrelation. lower the autocorrelation the smaller the MCMC standard error). There were no divergent transitions during sampling and none of the iterations saturated the maximum tree depth, indicating that the sampler is able to fully explore the posterior distribution. The values of ̂ and ef f ∕ for the scalar parameters of the BHM are summarized in Table 1.
To assess the adequacy and skill of the BHM, we begin by examining the model residuals at tide gauge sites (i.e., − − ). The measurement model (Equation 1) assumes that the residuals are normally distributed, and so evaluating the validity of this assumption is an important self-consistency check. To this end, we test the null hypothesis that the residuals come from a normal distribution using the one-sample Kolmogorov-Smirnov test. The test returns a p-value = 0.4, thus failing to reject the null hypothesis and, in turn, indicating that the residuals conform to the normality assumption. We also note that the time-mean of the residuals across the tide gauges is distributed symmetrically around zero, with a mean value of 0.002 mm and a standard deviation of 0.08 mm, indicating no systematic deviations between our estimates and the observations. The excellent skill of the BHM in capturing the general structure of the data can be appreciated more clearly by plotting estimates of RSL against tide gauge records at a selection of four stations (Figure 2). The long-term changes are very well captured by the BHM at all sites, even though such changes can differ significantly across sites due to spatial variations in sterodynamic sea level, GRD effects, and GIA as well as small-scale VLM. All of the records show considerable inter-annual to decadal variability, with peak-to-peak fluctuations as large as 10 cm, and this variability is generally well reproduced by the BHM, albeit with a tendency to underestimate its magnitude. The degree of underestimation varies across locations and largely depends on the extent to which the local sea-level variability can be explained by the five leading EOFs over the Mediterranean basin. This is not an issue since most of the unexplained variability is captured by the measurement error term ( ), thus ensuring consistency and inflating the uncertainty accordingly.
For comparison, we have also plotted estimates of sea-level changes from an EOF reconstruction at the four selected tide gauge sites (Figure 2). The EOF reconstruction is based on the same tide gauge data as the BHM and follows the formulation described by Calafat et al. (2014). At sites where VLM is small, the two reconstructions give similar time series of sea level, however at sites with strong VLM signals, such as Khalkis South (Figure 2c), the long-term changes are very different. This demonstrates an important point about EOF reconstructions. Because such reconstructions combine altimetry data with tide gauge observations without adequately accommodating the differences between GSL and RSL, it is unclear whether and to what extent the reconstructed long-term changes contain contributions from VLM. This is an important issue that has been addressed in our BHM.
We can form a more complete picture of the skill of the BHM in capturing the local sea-level variability by plotting the correlation between our estimates of detrended GSL changes and altimetry observations at each grid point (Figure 3a). In interpreting the correlation map, it is important to recall that the altimetry data are subsampled at every five grid points before being assimilated into the BHM while estimates are made on a 1/4° × 1/4° grid (the full altimetry grid). This means that the correlations with altimetry data largely reflect model skill at unobserved locations. The average correlation over the whole Mediterranean basin is 0.83, indicating a very good agreement. However, there is significant spatial variation in the correlations, with values close to one in the Adriatic, Aegean, and Ionian Seas, and relatively low values (∼0.6) in regions of high mesoscale activity (Fusco et al., 2003) such as the Alboran Sea, along the Algerian current, the Bonifacio Gyre, and along the Mid-Mediterranean Jet south of Crete. Mesoscale processes are characterized by small length scales (10-100 km), and so their effects on sea level are not well represented by the leading EOFs, nor can they be estimated with confidence at unobserved locations, which explains the lower correlations in the corresponding regions. Long-term sea-level changes are likely better captured by the BHM than the variability given that generally they are associated with larger length scales.
It is interesting to produce a similar correlation map but for the EOF reconstruction (Figure 3b). The regions of high and low correlation are the same as for the BHM, but, overall, the correlations are much lower, with an average value of 0.60. This is not surprising because the BHM assimilates both tide gauge observations and altimetry data, whereas the EOF reconstruction relies only on the tide gauge data to estimate the temporal amplitude of the altimetry-derived EOFs. In this regard, it is important to remark that the high correlations with the altimetry data in the case of the BHM (Figure 3a) are unlikely to be reflective of a similarly good performance over the pre-altimetry period. Looking now at the time series of sea level at an arbitrary point along the Algerian current (Figure 3c), we see that, unlike the BHM, the EOF reconstruction is, indeed, unable to reconstruct the sea-level variability in regions of high mesoscale activity. This can be interpreted as to mean that the existing tide gauge network does not provide enough information to capture the variability in such regions.

of 13
As an additional check of the BHM, we compare basin-average sea-level changes since 1960 as estimated by the BHM with those derived from the EOF reconstruction (Figure 3d). The two reconstructions agree very well on both the magnitude and the phase of the inter-annual variability, with a correlation of 0.92 for detrended time series. The long-term changes are also very similar in the two reconstructions, with an initial period of slow change from 1960 to 1990, followed by a period of rapid sea-level rise after 1990. Finally, we note that both reconstructions closely match the observed sea-level changes from altimetry for the period 1993-2018, giving us additional confidence in the skill of the BHM.

Results
We now investigate the BHM solutions with focus on long-term RSL changes. Hereafter, where we refer to total sea-level rates this includes the IB effect and excludes the contribution from small-scale processes ( ). We begin We show the posterior mean for both the total RSL changes (blue line) and the long-term RSL changes due to sterodynamic changes, gravitational, rotational, and deformation (GRD) effects, glacial isostatic adjustment (GIA), and small-scale processes (black line). The gray shading represents ±1 standard deviation (SD) associated with the estimate of long-term changes. 9 of 13 by analyzing time series of the basin-average rate of sea-level change for the period 1960-2018 (Figure 4a). The time series show that the total rate was relatively small and mostly negative before 1990, at which point it started to increase rapidly reaching a maximum value of 4.6 ± 0.6 mm yr −1 (posterior mean ±1SD) in 2000 that has largely persisted to this day. The average RSL rate for the periods 1960-1989 and 2000-2018 was, respectively, −0.3 ± 0.5 mm yr −1 and 3.6 ± 0.3 mm yr −1 (Figure 4b). Our approach also provides estimates of the contributions from individual sources to the total RSL rate. For the period before 1990, we find that sterodynamic changes and GRD effects made small contributions of opposed sign that largely canceled one another out, with average rates of −0.2 ± 0.5 mm yr −1 and 0.3 ± 0.2 mm yr −1 (Figure 4b), respectively. The contribution from GIA, which is constant over the whole period, is relatively small with a positive value of 0.2 ± 0.1 mm yr −1 . The IB effect made the largest contribution over this period with an average negative rate of −0.6 ± 0.1 mm yr −1 . For the period 2000-2018, sterodynamic changes and GRD effects both made positive and comparable contributions to the total rate of change, with values of 1.7 ± 0.3 mm yr −1 and 1.5 ± 0.2 mm yr −1 (Figure 4b), respectively. The contribution from the IB effect was relatively small with a value of 0.2 ± 0.2 mm yr −1 .
We now turn our attention to the local rate of RSL change, focusing on the period 2000-2018 as this is when the highest rates occur. The spatial structure of RSL rates is characterized by a single-sign pattern of positive rates that displays considerable geographical variation (Figure 5a). The highest rates (>4 mm yr −1 ) are found in the Tyrrhenian, Aegean, and Levantine Seas, whereas the lowest rates (<2 mm yr −1 ) occur mainly along the Cretan ) with those computed using an EOF reconstruction (green line). We show the posterior mean for both the total relative sea level (RSL) changes (blue line) and the long-term RSL changes due to sterodynamic changes, gravitational, rotational, and deformation (GRD) effects, glacial isostatic adjustment (GIA), and small-scale processes (black line). The gray shading represents ±1 standard deviation (SD) associated with the estimate of long-term changes. The basin-average sea-level changes from altimetry are also shown (red line).
Passage. Hence, since 2000, RSL has been rising everywhere in the Mediterranean at a rate that varies significantly across space.
Focusing now on the individual contributions, we find that the spatial structure of the sterodynamic rates closely resembles that of the total RSL rates, with values that range from −0.1 to 3.1 mm yr −1 (Figure 5b). The GRD rates are positive everywhere (Figure 5c), with values ranging from 1 mm yr −1 in the Northern Adriatic Sea and the Eastern Levantine Sea to >1.6 mm yr −1 in the Alboran Sea and south of the Ionian Sea. The relatively large gradient in the GRD rates in the North Western Mediterranean and the Adriatic Sea is associated with melting of the Greenland ice sheet. Indeed, the zero-line in the Greenland sea-level fingerprint is situated over the North Sea, with a consequent sharp gradient in the Northern Mediterranean. The sharp gradient in the eastern sector of the Levantine Sea reflects land uplift due to groundwater depletion in the Middle East (Rodell et al., 2018;Wada et al., 2010). The GIA contribution is relatively small and displays a pattern of meridional variation wherein rates increase from small negative values (−0.2 mm yr −1 ) along most of the African coast and in the easternmost part of the Levantine Sea to positive values of up to 0.5 mm yr −1 in the Western Mediterranean (Figure 5d). The contribution from the IB effect is small and positive over most of the basin, with values ranging from −0.1 mm yr −1 in the Eastern Levantine Sea to just over 0.4 mm yr −1 in the Central Mediterranean and the Aegean Sea (Figure 5e).
To provide a more quantitative description of the spatial structure of the rates, we have computed spatially averaged rates of RSL and those of the individual contributions over the five regions depicted in Figure 1 for the period 2000-2018 (Table 2). The Adriatic, Aegean, and Levantine Seas show the highest RSL rates, with spatially averaged values of 3.9 ± 0.5 mm yr −1 , 4.1 ± 0.4 mm yr −1 , and 3.9 ± 0.4 mm yr −1 , respectively. The lowest RSL rates are found in the central Mediterranean (3.5 ± 0.4 mm yr −1 ). The difference in RSL rates between these regions is largely explained by spatial variations in the sterodynamic rates (Table 2), which likely reflect a steric effect associated with changes in the thermohaline circulation of the Mediterranean Sea (as we discuss in the next section). Indeed, the spatially averaged sterodynamic rate over the Levantine Sea is almost twice that over the central Mediterranean (2.5 ± 0.3 mm yr −1 compared to 1.4 ± 0.3 mm yr −1 ). The GRD rates are similar and positive across the five regions (∼1.5 ± 0.2 mm yr −1 ), whereas the GIA rates are relatively small (up to 0.3 mm yr −1 ). The contribution from the IB effect is similar in magnitude to that from GIA, with a maximum value of 0.4 ± 0.2 mm yr −1 in the Aegean Sea and a minimum value of 0.1 ± 0.2 mm yr −1 in the Levantine Sea.

Conclusions and Discussion
Our analysis of data from tide gauges and altimetry has shown that spatially averaged sea level over the Mediterranean Sea fell at a rate of −0.3 mm yr −1 between 1960 and 1989 (Figure 4), in agreement with findings from past studies (Marcos & Tsimplis, 2008;Tsimplis & Baker, 2000;Tsimplis et al., 2005). This contrasts with a rate of global mean sea-level rise of 0.9 mm yr −1 in the same period (based on data from Frederikse et al., 2020). The fall in sea level during 1960-1989 was largely due to a positive anomaly in atmospheric pressure over the Mediterranean pushing water out of the basin, which previous studies have linked to the NAO (Gomis et al., 2008;

Subregion
Total ( Figure 1 12 of 13 Tsimplis et al., 2005;Tsimplis & Josey, 2001). However, even after accounting for the effect of this increase in atmospheric pressure, the rate of sea-level change in the Mediterranean (0.4 ± 0.5 mm yr −1 ) over that period was still less than half the global average rate (0.9 mm yr −1 ). The relatively low rate of IB-corrected sea-level rise in the Mediterranean in the period 1960-1989 has been noticed by previous studies (e.g., Marcos & Tsimplis, 2008;Tsimplis & Baker, 2000), but the reasons for this remained unclear. Here we have shown that the low rate was due to the small and opposing contributions from sterodynamic changes and GRD effects. We have also found that, since 2000, Mediterranean sea level has been rising at an increased rate of 3.6 ± 0.3 mm yr −1 , which is in line with the rate of global mean sea level rise for that period. This increased rate of rise has been driven by similar contributions from sterodynamic changes and GRD effects.
We have also shown that sea level changes in the Mediterranean Sea are highly non-uniform, with local rates that can deviate from the basin-mean rate by as much as 2 mm yr −1 for the period 2000-2018. This spatial non-uniformity is largely controlled by changes in sterodynamic sea level and is particularly remarkable in the Eastern Mediterranean where the sterodynamic contribution is characterized by rapid sea-level rise in the Levantine and Aegean Seas and significantly slower rise in the Ionian Sea and along the Cretan passage. This trend pattern emerges more clearly when plotting deviations from the basin-mean sterodynamic rate ( Figure 6). The pattern of deviations shows a striking contrast between negative deviations in the central Mediterranean and positive ones in all the adjacent seas (Adriatic, Aegean, and Levantine). This pattern is highly consistent with the changes in the thermohaline circulation and sea level that followed the end of the Eastern Mediterranean Transient (EMT) in 1995 (Calafat, Jordà, et al., 2012;Roether et al., 1996). The EMT started in 1987 and was characterized by a temporary shift in the primary sites of dense water formation in the Eastern Mediterranean from the Adriatic Sea to the Aegean Sea. This shift ended in 1995 when the Aegean stopped exporting significant amounts of dense water and the Adriatic Sea became again the primary deep-water source (Theocharis et al., 2002). These post-1995 changes in water density and circulation had a strong signature in sea level marked by a sea-level fall in the Ionian Sea and a rise in the Aegean and Levantine Seas (Calafat, Jordà, et al., 2012). The pattern of sterodynamic trends for the period 2000-2018 ( Figure 6) likely reflects some of the post-EMT changes, although confirming this would require a detailed hydrographic analysis, which is beyond the scope of this paper.
Our new approach provides coastal planners with a valuable tool to obtain robust and up-to-date estimates of local sea-level changes, while distinguishing long-term changes from short-term variability and informing about the causes of the changes. This will allow planners to make more confident and timely decisions regarding coastal adaptation for sea-level rise in the Mediterranean Sea. Additionally, since our new estimates are only based on observations, they can be used to evaluate simulated sea-level from climate models, on which sea-level projections are based, thus contributing to efforts to reduce model uncertainty and, in turn, further benefiting coastal planning decisions.

Data Availability Statement
The Bayesian solutions as well as the input data needed to run the BHM and replicate the results shown in this paper are available at Zenodo via https://doi.org/10.5281/zenodo.7046702 . The code that implements the BHM is also available at Zenodo via https://doi.org/10.5281/zenodo.6797648 (Calafat, 2022).