Modelling present and future climate in the Mediterranean Sea: a focus on sea-level change

We present results of three simulations of the Mediterranean Sea climate: a hindcast, a historical run, and a RCP8.5 scenario simulation reaching the year 2100. The simulations are performed with MED16, a new, tide-including implementation of the MITgcm model, which covers the Mediterranean—Black Sea system with a resolution of 1/16°, further increased at the Gibraltar and Turkish Straits. Validation of the hindcast simulation against observations and numerical reanalyses has given excellent results, proving that the model is also capable of reproducing near-shore sea level variations. Moreover, the spatial structure of the elevation field compares well with altimetric observations, especially in the Western basin, due to the use of improved sea level information at the Atlantic lateral boundary and to the adequate treatment of the complex, hydraulically driven dynamics across the Gibraltar Strait. Under the RCP8.5 future scenario, the temperature is projected to generally increase while the surface salinity decreases in the portion of the Mediterranean affected by the penetration of the Atlantic stream, and increases elsewhere. The warming of sea waters results in the partial inhibition of deep-water formation. The scenario simulation allows for a detailed characterization of the regional patterns of future sea level, arising from ocean dynamics, and indicates a relative sinking of the Mediterranean with respect to the Atlantic more pronounced than the current one.


Introduction
Future sea level rise (SLR) constitutes a threat for the coastal environment and economies (Hinkel et al. 2014), which is liable to be further exacerbated by the superposition of waves, atmospheric surge, and tides. Climate scientists are therefore becoming more and more committed to improve projections of SLR as to both resolution and accuracy, and to reduce current uncertainties, at the same time increasing our understanding of such a challenging scientific problem and enabling more accurate risk assessments. The latter, in particular, suffer from the difficulty in accounting for the cascading effects from sea-level rise to actual coastal impacts, across a range of spatial and temporal scales that demand complex and differentiated approaches (Thiéblemont et al. 2019). Likewise, the lack of comprehensive projections of extreme sea levels (ESL) that include mean sea level (MSL), tides, waves, and storm surges is lamented by Vousdoukas et al. (2017), who present a first attempt to provide an impact-oriented, regional-scale evaluation of ESL for the European shoreline. Their effort, albeit still inadequate for the design of specific measures, can nevertheless help highlight vulnerabilities and emerging research issues. As a matter of fact, the authors denounce the limitations inherent to their approach, mainly arising from the inadequacy of the regional atmospheric projections available at the time, as well as those of the projections derived from coupled Global Circulation Models (GCMs). In addition, they highlight that explicitly solving the nonlinear interactions between waves, storm surge, and tides should be a key element of any future effort to reliably model ESL, rather than continuing to resort to linear combinations of the ESL components resulting from independent simulations.
Besides being highly vulnerable to SLR due to the potentially disruptive impacts on its coastal economies (Jeftic et al. 1992), the Mediterranean basin is especially challenging from the scientific perspective, due to the inherent diversity of its geological history and the peculiar and complex features of its marine circulation and environment.
In addition to the local dynamics and geological processes, interacting over a broad spectrum of scales, sea level in the basin is also constrained by the water mass exchange across the Strait of Gibraltar, which, in fact, regulates the hydraulic jump between the Mediterranean and the Atlantic Ocean and influences how sea level rise in the Atlantic is transmitted into the Mediterranean. In its turn, the connection with the Black Sea, through the Dardanelles, Sea of Marmara and the Bosphorus, couples the basin hydrology to the land-based hydrological cycle of a vast portion of continental Europe. Nevertheless, despite the recognized inability of coarse resolution GCMs to accurately represent the highly non-linear, small-scale processes in marginal seas such as the Mediterranean (Slangen et al. 2017;Marcos and Tsimplis 2008), the global sterodynamic (following the terminology proposed by Gregory et al. 2019) sea-level projections for the Atlantic area near Gibraltar, are still often used to estimate the basin's internal sea level (Thiéblemont et al. 2019).
The question is still open whether the pronounced global trend observed in the last 25 years through satellite altimetry is the evidence of longer-term tendencies, or reflects more recent changes in the atmosphere-ocean coupled system.
Recent work by Dangendorf et al. (2019) provides evidence that the global trend accelerated at the end of the 1960s, and indicates that the resulting acceleration in the global sea level rise is linked to modifications of the southern hemispheric westerlies, leading to warming and circulation changes in the southern world ocean. Other factors have also been shown to play a key role, such as the change of terrestrial water storage, and the melting of ice sheets and glaciers, which has significantly increased in the last decades [see, e.g., Shugar et al. (2020), which focuses on the growth of glacial lakes, and King et al. (2020), where the mass loss from the Greenland Ice Sheet is analyzed]. The work by Frederikse et al. (2020) provides a thorough review of the various contributions and their relative weight.
As to future scenarios, since 1995 the coordinated worldwide CMIP effort (Coupled Model Intercomparison Project) has provided constantly updated projections of future sea level rise using GCMs and Earth-System Models (ESMs, which explicitly model chemical and biological interactions across sub-systems, thus closing the overall carbon cycle), to which the contributions from continental glaciers and geologic processes are added offline. The Phase 5 models (CMIP5), together with global glacier and ice-sheet models driven by CMIP5 GCM projections and complemented by estimates of land water storage, were used in the Fifth Intergovernmental Panel on Climate Change (IPCC) Assessment Report (AR5, IPCC 2013) to project a global sea level rise between 26 and 97 cm (overall likely range) at the end of this century, depending on the climate scenario, with pretty large inter-model spread. Such estimates were revised, mainly as regards the Antarctic contribution, in the 2019 Special Report on the Ocean and Cryosphere in a Changing Climate (SROCC), which assigns medium confidence to a likely range between 29 and 110 cm (IPCC 2019). For the area of the North Atlantic in proximity of the Strait of Gibraltar, projections can range from about 40 cm to over 100 cm for the most pessimistic RCP8.5 scenario, with an ensemble mean of about 80 cm, and approximately from 30 to 80 cm for the intermediate RCP4.5 scenario, with an ensemble mean of 60 cm (Vousdoukas et al. 2017). Slangen et al. 2014 derived regional sea-level projections and estimated their associated uncertainties up to the end of the twenty-first century, by combining model-and observationbased regional contributions of land ice, groundwater depletion and glacial isostatic adjustment, and CMIP5 projections of the changing ocean circulation, increased heat uptake and atmospheric pressure patterns. They also accounted for gravitational effects due to mass redistribution, ending up with an estimate of about 70 cm for the Mediterranean Basin. However, while concluding that regional variations in sea-level can significantly differ from the global mean (up to 30% above and 50% below) and that the land ice contribution dominates the overall uncertainty, they also stressed the need for dedicated regional downscaling of each global climate scenario.
In general, the inherent uncertainty in SLR projections has induced researchers to adopt a variety of alternative estimates in their impact assessments. Thiéblemont et al. (2019) analyze the effects of a median estimate of about 80 cm and two different high-end scenarios, resulting from two alternative extreme estimates of the sea-level equivalent of melting glaciers and of other individual contributions (i.e. the upper limit of the likely range and the "worst model" projections). Antonioli et al. (2017), while reviewing possible alternative ranges for the projected high-end SLR at 2100, use the 530-970 mm interval reported in the IPCC AR5, and Rahmstorf's semi-empirical estimate of about 1.400 mm (Rahmstorf 2007). Improved projections from the next generation of global models have been recently published by the AR6 WG1 and are available at https:// seale vel. nasa. gov/ ipcc-ar6sea-level-proje ction-tool. However, as remarked above, the spatial resolution of current global models is not sufficient to provide realistic estimates of local sea level rise in areas, such as the Mediterranean, where crucial processes are far from being explicitly resolved, and hardly allow a reliable parameterization (Sannino et al. 2009a, b). Adloff et al. (2018) recently discussed the complexity inherent to projecting sea level in the Mediterranean, and analyzed the performances of four different regional hindcast simulations of the basin circulation (period 1980-2012), driven by realistic atmospheric forcing. The incorporation of improved sea level information at the Atlantic lateral boundary was found to significantly enhance the reliability of results, confirming that the correct representation of the interactions between the two basins is an important requirement for a successful numerical simulation. However, an adequate treatment of the complex, hydraulically driven dynamics across the Strait of Gibraltar was still missing, as well as a more refined treatment of the exchange between the Mediterranean Sea and the Black Sea through the Dardanelles Strait, at the model eastern boundary. In addition, model performances were mostly evaluated by computing basin averages, while altimeter data series clearly indicate that sea level rise has not been homogeneous in the basin over the last decades. By analyzing altimeter data over a 25 years period , Mohamed et al. (2019) showed that the observed increase in the average sea level anomaly (SLA) has been quite different in different subbasins, ranging from a minimum of 1.95 mm/year, in the Ionian Sea, to a maximum of 3.73 mm/year, in the Aegean Sea. At a broader scale, the SLA positive trend appears to be significantly stronger in the Levantine basin than in the Western Mediterranean Sea, as well as characterized by distinctive features, with a fairly regular linear increase in the Western Mediterranean Sea, possibly due to the Atlantic contribution, and a more complex behavior in the Levantine basin. The latter, in particular, is liable to be influenced by the evolution of the Eastern Mediterranean Transient (EMT), the dramatic event occurred in the eastern Mediterranean at the beginning of the 1990s (Roether et al. 1996;Klein et al. 1999;Theocharis et al. 2002). The difficulty in separating climate-change-induced variations from the local circulation variability is therefore evident, as well as the role played by small-scale features, either attributable to internal variability or determined by atmospheric forcing, in generating longlasting differences across the Mediterranean sub-basins, amplifying or mitigating the effects of global sea level rise.
Prompted by these considerations, we developed a regional ocean model for the long-term simulation of the Mediterranean Sea circulation (hereinafter MED16) which we used to obtain high-resolution projections of the Mediterranean sea level. The model represents the climate version of the high-resolution, tide-including ocean model described in Palma et al. (2020). The two models share the same computational domain, which includes the Black Sea, thus allowing to consistently compute water exchanges at the Dardanelles Strait and to avoid any ad hoc assumption at the Mediterranean eastern boundary. The climate version necessarily uses a coarser horizontal grid with respect to its operational counterpart, yet resolution is significantly increased in critical regions such as the Gibraltar, Dardanelles, and Bosphorus Straits. The explicit, high-resolution representation of inter-basin water exchanges at the boundaries constitutes a unique, distinguishing feature of the present implementation.
In the following, we analyze a hindcast simulation of the Mediterranean circulation and the corresponding historical and RCP8.5 scenario, all driven by downscaled atmospheric fields obtained with the SMHI-RCA4 regional model (Strandberg et al. 2014). The overall basin-scale model performance is assessed through comparison with observations and reanalysis data, with special focus on the model ability to reliably represent the local sea level height. In particular, comparison with data from tide gauges allows us to evaluate the model performance in coastal regions, where the impacts of sea level rise really affect human communities and economies, and need to be specifically assessed.
The paper is organized as follows. The main characteristics of the model and of its implementation are described in Sect. 2, whereas Sects. 3-5 are devoted to the detailed analysis of the simulations performed. In particular, Sect. 3 presents a validation of the hindcast and historical simulations in terms of transports at the main straits, hydrology, and circulation, obtained through comparison with observations and numerical reanalyses. The corresponding sea-level reconstructions are discussed in Sect. 4, and validated using satellite altimeter data and coastal tidal gauge observations. Section 5 then describes the results of the scenario simulation, highlighting future changes in the basin hydrology and circulation, and discussing the projected sea-level. Finally, conclusions are drawn in Sect. 6.

Data and methods
In order to produce the downscaled Mediterranean sealevel field under present climate  and the RCP8.5 future scenario (2006-2100), the MED16 model was forced using the correspondent regional downscaling experiments performed with the SMHI-RCA4 atmospheric model (Strandberg et al. 2014), alternatively constrained by (a) the ERA-Interim reanalysis data (Dee et al. 2011), for the hindcast simulation, (b) the HadGEM2-ES global model (Collins et al. 2011) for the historical simulation, and (c) the HadGEM2-ES RCP8.5 for the future scenario. The current section describes in detail the atmospheric data used (Sect. 2.1) and the MED16 model characteristics and specific setup (Sect. 2.2).

Atmospheric forcing
The atmospheric forcing is derived from the dynamically downscaled regional atmospheric fields produced by the Rossby Centre regional atmospheric model RCA4 (Strandberg et al. 2014), at 0.11° resolution (i.e. approx. 12.5 km grid spacing), over the EURO-CORDEX domain (Giorgi et al. 2009), whose boundary conditions are either provided by the ERA-Interim reanalysis, or by the atmospheric component of the CMIP5 global model HadGEM2-ES. HadGEM2-ES is especially suitable as a global driver for downscaling studies over Europe due to its outperforming other ESMs along the lateral boundaries of the Euro-COR-DEX domain (Brands et al. 2013). No bias-correction is applied to the atmospheric forcing fields. The hindcast experiment covers the period 1981-2010, the historical experiment covers the period 1981-2005, while the RCP-8.5 scenario spans years from 2006 to 2100. Up to 2005, observed greenhouse gas concentrations have been prescribed, which are then substituted by those indicated by Meinshausen et al. (2011) for the future business-as-usual scenario RCP8.5. To reproduce the surface heat fluxes, shortwave radiation from the atmospheric models is imposed, whereas the wind stress and the other heat flux components are computed via bulk formulas considering the sea surface temperature, the winds at 10 m height, the dry air temperature and the air pressure at 2 m, and the relative humidity as inputs. In particular the long-wave radiation is computed according to the formula proposed by Bignami et al. (1995), while the latent heat flux, the sensible heat flux and wind stress are computed according to the Large and Yeager (2004) bulk formula. Cloud cover is taken from the atmospheric model. The net freshwater flux is computed as precipitation (taken from the atmospheric model) minus evaporation (computed from the latent heat). To reduce salinity drift in the hindcast simulation the sea surface salinity (SSS) is restored toward the MEDHYMAP hydrological monthly climatology (Jordà et al. 2017a) at a timescale of two days over the topmost model layer of 2 m thickness. At the surface, the model is forced by 6-hourly wind, and 3-hourly surface pressure, heat and freshwater fluxes computed via bulk formulae. The fresh water flux is prescribed in conjunction with the non-linear free surface numerical scheme.

The regional ocean model MED16
The numerical ocean model MED16 used to simulate the downscaled sea level is the updated version of the Mediterranean ocean model developed by Sannino et al. (2015). It is based on the MITgcm kernel developed by Marshall et al. (1997a, b), and solves the Navier-Stokes equations under the Boussinesq approximation for an incompressible fluid. In the present study, the hydrostatic version of the model has been implemented, using a finite-volume spatial discretization on a staggered Arakawa-C grid, partial step topography, rescaled vertical height (z*) coordinate , and an implicit non-linear free surface formulation ). The source code and documentation are available at the following web site: https:// github. com/ MITgcm/ MITgcm (last access: 20 April 2021).
Model bathymetry for both the Mediterranean and the Black Sea was derived from the European Marine Observation and Data Network (EMODnet) 2016 dataset (https:// www. emodn et-bathy metry. eu), while the high-resolution digitalized chart of Sanz et al. (1991) provided data for the SoG, and the high-resolution bathymetry for the Bosphorus and Dardanelles Straits was made available by Gökaşan et al. (2005Gökaşan et al. ( , 2007, with the permission of the Turkish Navy, Navigation, Hydrography and Oceanography Office. The datasets were combined through a bilinear interpolation on the computational grid, followed by manual check for isolated grid points, islands, and narrow passages (see Sannino et al. 2015Sannino et al. , 2017. Both equilibrium tide (i.e. the generating potential is incorporated in the momentum equations as a body force) and tide propagating from the Atlantic Ocean through the Atlantic open boundary, are explicitly applied as in Naranjo et al. (2014), Sannino et al. (2015), and Palma et al. (2020). The equilibrium tide is composed of four tidal components: M2, O1, S2, K1, the semidiurnal and diurnal principal lunar tides, the principal semidiurnal solar tide, and diurnal lunisolar declination tide, respectively. The tidal values used to prescribe the Atlantic tide are derived from the OTIS global tide inverse model (Egbert and Erofeeva 2002). Differently from the model configuration of Naranjo et al. (2014), and Sannino et al. (2015) the updated version used in this study includes the Black Sea that is interactively connected to the Mediterranean Sea through the Turkish Straits (Dardanelles and Bosphorus). The model domain therefore covers the whole Mediterranean-Black Sea system and a small part of the Atlantic Ocean west of the Strait of Gibraltar (SoG hereafter), whose western limit corresponds to the only open-boundary condition prescribed. The horizontal computational grid has a uniform resolution of 1/16° (about 7 km), which is significantly increased in correspondence of the three straits where higher resolution is needed to reliably represent the local dynamics, i.e., the SoG (maximum resolution of about 200 m), and the Dardanelles and Bosphorus Straits, where a smoothly varying refinement in the latitudinal and longitudinal directions allows to reach a maximum resolution of 1/200° (about 555 m) (Fig. 1). The vertical domain is discretized using 100 z-levels, with grid spacing increasing from 2 m near the surface to 62 m at a depth of 1500 m; below 1500 m a uniform grid spacing of 62 m is used down to the sea floor. The time-step used is one minute. To calculate the vertical eddy viscosity and diffusivity coefficients, we used the 1.5 order Turbulent Kinetic Energy (TKE) closure scheme by Gaspar et al. (1990), adapted from the atmospheric case developed by Bougeault and Lacarrere (1989). The background vertical eddy viscosity and diffusivity were set to 1.5 10 -6 m 2 /s and 10 -7 m 2 /s, respectively. The maximal value of diffusivity allowed to be generated by the turbulence scheme is 100 m 2 /s, in order to let it handle unstable vertical stratification and avoid the use of an enhanced vertical mixing parameterization. A spatialdependent horizontal viscosity is obtained from the turbulence closure scheme by Leith (1968). Differently from the well-known scheme by Smagorinsky (1963), Leith's scheme focuses on resolving the direct enstrophy cascade towards the smaller scales, which is characteristic of 2D turbulence (Fox-Kemper and Menemenlis 2008). A constant horizontal diffusivity coefficient (2 m 2 s −1 ) is applied with a laplacian operator for the tracers (T,S). The no-slip conditions were used at the lateral and bottom solid boundaries, along with a quadratic bottom drag. The latter is calculated as a function of the velocity close to the bottom, with a dimensionless coefficient of 0.0025. The selected tracer advection scheme is a third-order direct space-time flux limited scheme.
As the low spatial resolution of the global model did not allow an accurate representation of the dynamics within the Strait of Gibraltar, the Mediterranean basin was completely detached from the Atlantic Ocean in the global simulation, by artificially closing the Strait. The Mediterranean was therefore virtually represented as a lake in the global projection, thus impairing its use to initialize the regional ocean simulation. Initial conditions (namely temperature and salinity) were thus derived from MEDHYMAP climatological data. To reduce spurious drift, MED16 was spun-up 35 years, using the hindcast forcing for the 1980-1986 period in a five-time loop. The span-up model was then used as initial condition for both historical and hindcast simulations.
At the Atlantic open boundary, consistently with the surface atmospheric forcing, the hindcast simulation was forced with ORAS4 global reanalysis (Balmaseda et al. 2013), while for historical and RCP8.5 scenario simulations monthly data from the HadGEM2-ES global model were used. In particular, the lateral boundary conditions consist of monthly mean temperature, salinity and sea surface height (SSH) which are interpolated onto the MED16 grid. Temperature and salinity are applied through a 3D relaxation with a relaxation time varying linearly from 2 h at the western limit of the domain to 30 days over the first 30 grid points. The prescribed SSH on the Atlantic box includes different contributions depending on the forcing. For the hindcast simulation, the prescribed ORAS4 SSH includes sea level contributions from ice sheet mass loss, glaciers ice melt, changes in land water storage, as well as global thermal expansion. As ORAS4 underestimates the regional SSH seasonal cycle in the near Atlantic region compared with the multi-satellite products provided by the Climate Change Initiative (CCI-SLA doi: https:// doi. org/ 10. 5270/ esa-sea_ level_ cci-MSLA-1993_ 2013-v_ 1.1-201412), we have applied a 12-month correction based on the difference between CCI-SLA and ORAS4 as described in Adloff et al. (2018). For both historical and scenario simulations, the prescribed SSH includes the dynamic sea level (DSL) and the global mean thermal expansion. DSL is the height of the ocean with respect to the time-invariant geoid that is determined by the dynamical balance associated with ocean density and circulation. DSL includes the regional variability of dynamic topography changes due to water mass advection, thermohaline circulation and to the wind-driven circulation (Gill and Niller 1973) and are part of the standard output in the CMIP5 simulations (variable zos). zos was initially dedrifted by removing the linearly fitted HadGEM2-ES control run (which is forced by non-evolving pre-industrial conditions) from each grid point individually. This step removes any artificial signals associated with ongoing spin-up deep ocean and/or limitations in the representation of energy conservation in the model domain, as discussed by Sen Gupta et al. (2013). As the ocean component of HadGEM2-ES is based on the Boussinesq approximation and conserves volume rather than mass (Greatbatch 1994), the value of zos was further corrected (at each time step and grid point) by subtracting its time dependent global mean. This correction guarantees that the resulting sea-level patterns only reflect fluctuations due to the joint effects of ocean density and circulation (Gregory et al. 2019) and thus, it is comparable to those resulting from non-Boussinesq models (Losch et al. 2004;Griffies and Greatbatch 2012).
The corresponding global thermal expansion time series (variable zostoga) was also corrected for control drift by removing a quadratic fit to the control run's thermal expansion time series before being added to the detrended and zero global mean zos field. The resulting sea level variable was then used as SSH lateral boundary condition for the MED16 historical and scenario simulations. As in Naranjo et al. (2014), Sannino et al. (2015), and the higher resolution Mediterranean version of Palma et al. (2020), SSH and the additional 4 tidal constituents were directly prescribed at the western open boundary.
The remaining sea level components for the historical and scenario simulations, namely the surface mass balance and dynamic ice sheet contributions from Greenland and Antarctica, the glacier and land water storage contributions, and the Glacial Isostatic Adjustment (GIA), were extracted from the Integrated Climate Data Center (ICDC) at Hamburg University (https:// icdc. cen. uni-hambu rg. de/ en/ ar5-slr. html), interpolated on the MED16 grid, and added offline.
The ICDC dataset is based on the ensemble mean of CMIP5 climate models that were used in the regional sea-level projections of the AR5 (Church et al. 2013a) of the IPCC. ICDC provides data on a global grid at 1° × 1° spatial resolution for different RCP scenarios. In this study, we used central estimation data of the RCP8.5 ensemble mean.
River discharge at the river outlets is prescribed by transporting the daily total runoff computed by the forcing atmospheric regional model, along the river network of the WBMplus model (Vörösmarty et al. 1998;Wisser et al. 2008Wisser et al. , 2010, by means of a Muskingum-Cunge type scheme that solves the Saint-Venant flow equations (Ponce 1994). Only for rivers discharging in the Black Sea, the so computed discharge was bias-corrected to reproduce monthly climatological values reported in the literature (Kara et al. 2008).

Model validation
Before analysing the model sea level reconstructions (see Sect. 4), here we perform a first validation of the hindcast and historical runs. In particular, we compare the simulated water transports at the straits and the Mediterranean hydrology and circulation with the observations, and with the results of a recent reanalysis of Mediterranean water properties by Escudier et al (2020). The latter is a product from the Copernicus CMEMS database (http:// marine. coper nicus. eu) that covers the period 1987-2019, and results from the assimilation of temperature and salinity vertical profiles and satellite Sea Level Anomaly along track data into a numerical system, consisting of the Nucleus for European Modelling of the Ocean (NEMO), and of a variational data assimilation scheme (OceanVAR).
As our focus is on the Mediterranean Sea, we will mainly show results for this basin, and report the projected water exchanges at the Dardanelles as a boundary condition, determined by the Black Sea effectively acting as a reservoir. Nevertheless, it is worth noting that the model correctly reproduces the permanent basin-scale, cyclonic boundary current that characterizes the Black Sea circulation.

Transports at the main straits: Gibraltar, Dardanelles, Bosphorus
The explicit treatment of water inflow and outflow at the open boundary of the Mediterranean Sea, coupled with explicit tidal forcing, allows to realistically constrain the mass, energy and momentum exchanges with the Atlantic, as well as to finely modulate the temperature and salinity properties of either flow, in both time and space (see Sannino et al. 2015, for a detailed description).   (2002)]. The net transport for the historical experiment is in excellent agreement with the hindcast value, although exhibiting reduced variability, yet it results from higher inflow and outflow.
The net transport across both Turkish Straits (TSS) is on average negative, consistently with the observed net inflow of fresher waters into the Aegean Sea from the Black Sea, where precipitation and river discharge exceed evaporation. The magnitude of the inflow appears to be lower than previous observational and numerical estimates, which range from 0.006 to 0.01 Sv (Ünlüata et al. 1990;Simonov and Altman 1991;Peneva et al. 2001;Kara et al. 2008;Sannino et al. 2017). However, fluxes through the TSS are difficult to accurately measure, due to the large variability in the observed currents and the consequent low signal-to-noise ratio, which might need longer-term monitoring to be adequately increased (Sannino et al. 2017). On the other hand, the simulated net transport is expected to significantly benefit by further refining the grid in correspondence of both Straits, as demonstrated by Sannino et al. (2017), where a very-high-resolution simulation was performed to accurately reproduce transports with the prescription of climatological barotropic forcing at the two open boundaries. Increased resolution together with a specific treatment of local viscosity and diffusivity is therefore liable to enhance water exchanges across the straits, and to further improve the already good model performance.

Basin hydrology and circulation
We first compare model climatologies of the winter (DJF) and summer (JJA) sea surface temperature (SST) and SSS with reanalysis data; mean fields are computed over the common period 1987-2005. In Fig. 2, panels (a, b) are for the Sea Surface Temperature Copernicus reanalysis panels (c, d) and (e, f) show the bias with respect to reanalysis for the hindcast, and the historical simulation (i.e., the difference between the seasonal climatological means), respectively, and panels (g, h) illustrate the RMSE for the hindcast simulation only. RMSE is not shown for the historical simulation as temporal disalignments might artificially inflate the error estimate. Model simulations are in overall good agreement with the reanalysis, even though MED16 systematically produces lower temperatures in winter (more so in the historical simulation). In summer, the SST from the hindcast is in excellent agreement with the reference field, whereas the historical simulation displays slightly warmer temperatures in the western basin, in the central part of the eastern basin and in the Adriatic Sea, and lower ones in the Ionian, and in the easternmost part of the Levantine. For the hindcast experiment, the RMSE pattern is generally very similar to that of the bias, showing that discrepancies with respect to the reference reanalysis mainly correspond to areas of enhanced variability, and confirming the systematic winter underestimation of SST.
In Fig. 3, the SSS model climatologies (panels c-d, hindcast; panels e, f, historical), are compared with the Copernicus reanalysis (panels a-b). In the western basin, the hindcast salinity appears to slightly overestimate the observations in both seasons, although less notably so in summer.
In the eastern basin, the salinity field is generally well reproduced in winter, while in summer some local negative biases are present (e.g., in the Gulf of Gabes, in the Adriatic Sea, and in the north Aegean). The historical experiment, on the contrary, exhibits larger salinity values in both seasons and basins, yet reducing the positive bias around the Balearic Islands, coherently with a limited intrusion of the Atlantic Water (AW) into the Tyrrhenian Sea, the Sicily channel and the Gulf of Gabes, and with a more pronounced northward spreading of the Atlantic inflow after exiting the Alboran Sea. Specific local differences between the numerical experiments and the reanalysis can be attributed to the different treatment of rivers and/or to the different magnitudes of the runoff provided by the driving atmospheric simulations (e.g. along the eastern Adriatic coast and in the North Adriatic Sea), while the different behaviour of the two MED16 runs in the North Aegean is probably related to the representation of the local circulation, as no significant difference was detected in the representation of the net inflow from the Dardanelles between the hindcast and the historical experiment. This discrepancy will be the object of future research. The RMSE pattern (panels g, h) confirms such findings.
Vertical means of the modelled and reanalysed temperature and salinity were computed over the 1987-2005 time window, for three different vertical layers, and averaged over the whole Mediterranean basin (MED), and the western and eastern sub-basins (WMED and EMED). Table 2  The current experiments are in good agreement with the available observations and the bias is coherent with that of similar experiments. In view of the discrepancies among regional models in describing heat and salt redistribution within the Mediterranean basin (Llasses et al. 2018), and of the differences in the spatial patterns highlighted above, the agreement between the hindcast and historical simulations and the reference data is indeed acceptable. A semi-quantitative discussion about this point is presented in the Appendix.
The modelled temperatures are always slightly lower than the observations in the two upper layers and slightly higher in the deep layer. For both simulations and all geographical areas, maximum discrepancies are anyway found in the surface layer, with largest differences in the Eastern basin. When the whole basin is considered, modelled salinity also generally presents negative anomalies, with maximum negative differences in the surface and intermediate layers of the hindcast simulation, while slightly positive anomalies can be observed in the surface and bottom layers for the historical simulation. When separately considering the two sub-basins, the hindcast experiment exhibits more marked negative biases in the intermediate layer of the Western basin and in the surface layer of the Eastern basin, while the historical simulation projects higher salinity in the surface layer of the Eastern basin with respect to the reanalysis.
In order to evaluate the ability of the hindcast experiment to reproduce the time evolution of water temperature, Fig. 4 compares the modelled time series  of the annual temperature anomalies to those of the reanalysis, in the surface and intermediate layers, for the whole basin and the western (WMED) and eastern (EMED) sub-basins. Anomalies are computed relative to the respective means over the whole period. Comparison with the historical simulation is here omitted as time coherence with the reanalysis could not be expected. The observed interannual variability is generally well captured by the hindcast simulation in both layers, The EMED time series, in particular, captures the positive trend in both the upper and intermediate layers.
Moreover, the simulation well represents the strong cooling during the years 1992 and 1993, that is the signature of the Eastern Mediterranean Transient (EMT). Table 3 reports the values of the RMSE and of the correlation coefficient (r) of the temperature annual time series with respect to the reference reanalysis. The RMSE is systematically larger in the surface layer than in the intermediate layer, while the relative weight of contributions from each sub-basin to the overall basin mean is reversed in the two layers, with larger (smaller) errors in the surface (intermediate) layer of the EMED with respect to the WMED. In view of the similar conclusions derived from the analysis of the tabulated mean biases, and in consideration of the high values of the correlation coefficient reported in Table 3, the RMSE indeed appears to be bias-dominated, consistently with the above discussion of the spatial patterns presented in Fig. 2. To this regard, the very low value of r in the WMED intermediate layer is not significant, as here the two temperature curves are virtually constant, only exhibiting the small-amplitude, either uncorrelated or phase shifted variability depicted in the figure, while the overall correlation at the basin scale is determined by the EMED trend. The presented RMSE values can be qualitatively compared to their analogues reported in the Copernicus Reanalysis Quality Information Document, which are of the same order of magnitude, yet systematically larger, as arising from  (a, b), difference between hindcast and reanalysis (c, d), difference between historical simulation and reanalysis (e, f). Rmse between hindcast and reanalysis (g, h). Winter (summer) climatology on the left (right) panels ◂ bi-weekly statistics that are intrinsically noisier, despite the reference observation being only "quasi-independent" from the reanalysis.
Analogous considerations apply to the time series of the mean annual salinity anomalies (Fig. 5), and to the correspondent RMSE and correlation coefficients (Table 3). For the whole Mediterranean Sea, the RMSE is now of comparable magnitude across the two layers, yet somewhat larger in the intermediate layer. The two salinity curves are generally well correlated, except for the intermediate layer in the WMED, where they are, again, quasi-constant over the period yet affected by uncorrelated low-frequency variability, causing r to be very small and ineffective at the basin scale, where the systematic EMED trend emerges. In the upper layer of the WMED, the model captures both the sudden drop in salinity around 1990 and the maximum in 2005. The first event is due to an excess of freshwater input and is therefore confined to the upper layers, whereas the latter is related to a well-known event of intense dense water formation in the Gulf of Lion, which is however too localized to be captured in the spatial mean of the intermediate layer.
A good representation of the interannual variability is also achieved in the EMED surface layer, with a local maximum during the 90s and a positive trend in the second half of the simulation.
The only notable discrepancy between the modelled salinity and the reanalysis occurs in the years from 2001 to 2004, during which the reanalysis displays a relative minimum, particularly strong in the WMED, which is not reproduced by the model. However, a discrepancy with respect to observations for this period is in fact reported in the Copernicus Reanalysis Quality Information Document, while the corresponding data from the MEDHYMAP dataset are in good agreement with MED16 also for years 2000-2005. Figure  The main Mediterranean surface current systems, either driven by the local wind forcing, or by the inflow of AW at the Gibraltar Strait, are present both in our simulations and in the reanalysis, with similar shapes. This also applies to most of the robust cyclonic and anticyclonic circulation structures, even though there are a few features, such as the quasi-permanent cyclonic circulation to the southeast of Sardinia, and some smaller vortices in the eastern basin, that are not well defined. Overall, considering that our results derive from almost "free" simulations spanning more than two decades, whereas the reanalysis relies on the assimilation of observed data, the representation of the average Mediterranean circulation provided by our models can be considered quite satisfactory.
Three hindcast sub-regional surface circulation maps are displayed in Fig. 7. Panel (a) of the figure highlights the two permanent cyclonic circulations that characterize the northern portion of the western Mediterranean, that is, the wide gyre formed by the Liguro-Provençal current, occupying the western Ligurian Sea, and the Bonifacio cyclone, in the northern part of the Tyrrhenian Sea, to the east of Corsica. Interestingly, between the two cyclonic circulations, there is the weak signature of an anticyclonic circulation, in the area of the Corsica Channel, between Cap Corse (the long peninsula in the northern tip of Corsica) and the Italian coast. This structure, which was denoted as Ligurian anticyclone in Ciuffardi et al. (2016), appears to be a robust feature of the local dynamics during the summer months (see also recent works by Poulain et al. 2020, and by Iacono and Napolitano 2020).
As shown in Fig. 7b, the model correctly reproduces the wide cyclonic structure, with two poles, that occupies the southern, deeper portion of the Adriatic Sea [see, e.g., recent work by Palma et al. (2020), and references therein].
Finally, we note that the surface dynamics of the Marmara Sea (Fig. 7c) are dominated by a central, elongated gyre, more attached to the northern coast, in agreement with what was found in a recent high-resolution numerical investigation of the Turkish Straits System by Sannino et al. (2017).
The time evolution of the Mixed Layer Depth (MLD- Fig. 8) highlights how deep-water-formation processes in the Gulf of Lion and the Adriatic basin are well represented both in the hindcast and in the historical experiments.

Hindcast and historical sea-level analysis
The Mediterranean Sea level average and variability simulated in the hindcast experiment are compared to their observed counterparts, derived from altimetric satellite sea level anomaly (SLA) measurements. The model ability to reproduce nearshore sea level is assessed by means of tide gauge data.  (a, b), difference between hindcast and reanalysis (c, d), difference between historical simulation and reanalysis (e, f). Rmse between hindcast and reanalysis (g, h). Winter (summer) climatology on the left (right) panels ◂

Comparison with satellite observations
The altimeter data used are daily maps of gridded L4 SLA and derived variables reprocessed (Copernicus portal, http:// marine. coper nicus. eu, last access April 2021), obtained from the merging of track data from different satellites. These maps, especially developed for the Mediterranean Sea, have a nominal resolution of 1/8° × 1/8°, which is twice the resolution of the products that cover the world ocean. They have been available since 1993, overlapping the last 18 years of the simulation. A different altimeter dataset was used in Adloff et al. (2018) (see Ablain et al. 2015, for details), which consists of monthly SLA maps with a spatial resolution of 1/4°.
The average seasonal cycle reproduced by the hindcast is compared to that of the observations in Fig. 9, where the solid lines with black dots refer to model results and the dashed lines with diamonds to altimeter data (the averages are over the period 1993-2010). Panel (a) shows that the model cycle over the whole basin is fairly close to the observations, with some overestimation in winter and some underestimation at the end of summer and in part of the autumn.
Overall, the hindcast model cycle is close to that obtained from the NEMO-MED12 simulation discussed in Adloff et al. (2018), which was found to provide the best SL reconstruction among the four considered in that work.
The two other panels in the figure show the same comparison for the western and eastern sub-basins, respectively. The cycle for the western Mediterranean is even closer to the observations than the basin cycle, whereas in the eastern basin, the model underestimation during SON is also larger. This may be due to the fact that a fraction of the seasonally varying steric contribution to SL, which is in principle absent in most state-of-the art circulation models, is here indirectly introduced through the Atlantic boundary conditions. Its effects in fact appear to be confined to the Western basin, while additional basin-wide sea-level components induced by local heat fluxes are not explicitly represented in the regional model. The correspondent curves for the historical experiment (dot-dashed line) are also reported in Fig. 9, and show a behaviour similar to that of the hindcast run, although the overall maximum is anticipated to the midsummer and the positive phase of the cycle appears to be longer.
The interannual variability of the sea-level anomaly of the hindcast simulation is shown in Fig. 10, where the same notations of the previous figure are adopted. The plotted values are yearly anomalies with respect to the mean over the period covered by the altimeter (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010), averaged over four different regions: the Atlantic box, the whole Mediterranean basin, the WMED, and the Levantine sub-basin. The first panel is relative to the Atlantic box, i.e., the small domain to the west of the Gibraltar Strait where boundary conditions are prescribed, making use of the available altimetric data (see Sect. 2 for details). From 1993 on, the model time series appears well correlated with the observations, showing that the boundary conditions are correctly applied in our model. In fact, the cross-correlation coefficients of Table 4 show that the model time series are very well correlated with the corresponding altimeter-derived time series in all the basins in consideration.
The time series of the average Mediterranean SL anomaly (top right panel in Fig. 10) reproduces the growing basinscale trend displayed by the altimetric observations (see Bonaduce et al. 2016, and the more recent investigation by Mohamed et al. 2019, where trends of the SL in the Mediterranean were discussed, based on 25 years of altimeter data). It also captures short-time variations, such as the sharp increase seen in 2010, which results from a large wintertime Table 2 For each geographical area in the leftmost columnthe whole Mediterranean (MED) and the Western (WMED) and Eastern (EMED) sub-basins-the first row reports the average temperature (°C) and salinity (Psu) at various depths, from the Copernicus Reanalysis dataset; the second and third rows (labelled "hindcast" and "historical") respectively report the associated differences between the hindcast and the historical simulations and the reanalyses. All spatial averages are taken over the period [1987][1988][1989][1990][1991][1992][1993][1994][1995][1996][1997][1998][1999][2000][2001][2002][2003][2004][2005] Landerer and Volkov (2013), this was due to a basin-wide barotropic fluctuation induced by the variation of the wind stress in the Gibraltar area and in the neighboring Atlantic sector, associated with inflow of Atlantic water through the Gibraltar Strait. A similar event occurred during the following winter of 2011. Comparison with Figure 4 of Adloff et al. (2018) shows that the model time series of Mediterranean anomalies also captures the main interannual variability observed before 1993; the reference for this period is provided by two reconstructions of the Mediterranean SL by Jordà (2011), andMeyssignac et al. (2011), based on data from coastal tide gauges. The series is also in good agreement with that obtained from the NEMO-MED12 simulation of Adloff et al. (2018). This is due to the fact that, despite the many differences, the two models apply the same Atlantic boundary conditions, and use as local atmospheric forcing two different downscaling of the same dataset (ERA-INTERIM). In some cases, the model series displays stronger SL variations; this may be ascribed to model differences, but could also partly be due to the contribution of mesoscale features that are better resolved by the model, thanks to the higher horizontal resolution (1/16°). This holds a fortiori, when compared with the altimeter data, whose nominal resolution is 1/8° × 1/8°. The lower row of Fig. 10 shows the SL anomalies for the sub-basin that is closest to the Atlantic (western Mediterranean) and for that more far away from it (Levantine basin, to the east of 22° of longitude). The signals are quite different in the two sub-basins; the series in the western Mediterranean resembles that in the Atlantic box, with a fairly regular increase of the SL over the simulation period, whereas that in the Levantine shows a decrease of the SL at the end of the 80s and beginning of the 90s, followed by a sharp increase, leading to a different state over the last 10 years of the simulation. The latter behavior is of course related to the Eastern Mediterranean Transient (EMT), the major event that has affected the circulation and water mass properties of the Mediterranean in the last decades. The hindcast is capable of reproducing the rapid SL rise as a result of strong anomalous mass modification in the deep layers of the EMED during the end of the 80s and the beginning of the 90s; the effects of the EMT on sea level has been discussed for the first time by Tsimplis and Josey (2001). The model correlates very well with the altimeter observations after the transient, which display a definitely nonlinear trend, first evidenced in Amitai et al. (2010).
We can conclude that the hindcast simulation not only captures the basin scale variability of the SL over the hindcast period, but also some expected, significant differences between the main Mediterranean sub-basins.
In Fig. 11, the modelled (both hindcast and historical) SSH fields, averaged over the period 1993-2011, are compared with the corresponding climatological average of the absolute dynamic topography (ADT) from altimeter data (AVISO, http:// www. aviso. ocean obs. com, last access May 10th 2021). In all cases, the respective regional climatological means have been subtracted to ease the comparison of patterns. Being a long-term average, the field derived from the AVISO data is expected to be very close to the Mean Dynamic Topography (MDT), that is, to the mean signal, representative of the average circulation, which produces the ADT when added to the zero-mean remotely measured sea level anomaly (SLA). The reference observational field shown in Fig. 11 is indeed very close to that of Rio et al. (2014), which was constructed using observations and model results, and validated using independent observations. The spatial pattern of the modelled SSH field generally appears to be in good agreement with the observations, and is also similar to that obtained by the NEMO-MED12 model, the best performing among those considered in Adloff et al. (2018). As to the hindcast experiment, the MED16 model outperforms NEMO-MED12 in the western portion of the basin, indicating that the stretched grid at Gibraltar indeed allows for a better representation of the complex dynamics inside the Strait. Nevertheless, the signatures of the Rhodes gyre and of the western Cretan gyre are weaker than expected. On the other hand, neither the MED16 simulation nor the four hindcasts analyzed in Adloff et al. (2018) resolve the Ierapetra gyre, possibly due to the still insufficient spatial resolution or to the quality of the downscaled wind forcing in the area, the Ierapetra gyre being a wind driven feature.
The SSH values obtained for the Black Sea (not shown here) exhibit larger spatial variability with respect to those of the Mediterranean Sea, approximately ranging from − 6 to 18 cm, and compare well with the MDT reconstructions by Knysh et al. (2002) and by Kubryakov and Stanichny (2011). In particular, Kubryakov and Stanichny (2011) find localized, deeper minima in the core of the dominant cyclonic cell, while MED16 produces significantly higher SSH values than either reconstruction in the westernmost portion of the same cell.

Comparison with tide gauges
Hindcast sea surface heights along the coast have been compared with tide gauge data retrieved from the Permanent Service for Mean Sea Level (PSMSL) (Holgate et al. 2013).
In particular, we used time series of sea level monthly means from the Revised Local Reference dataset (RLR), which are independently calibrated for each station in the PSMSL according to the tide gauge history provided by each local supplying authority. However, no correction is made by PSMSL for vertical land movements (VLM), although tide gauges without known benchmarks are omitted from this data set by PSMSL prior to release (Hamlington et al. 2016). These data therefore measure the relative movement of the sea surface with respect to the land, disregarding the multiple processes over a variety of time scales that can affect MSL estimates, some of natural origin, such as the glacial isostatic adjustment (GIA) and plate tectonics, and some related to human activities (e.g. ground-water pumping). In the following, RLR data have been used without any further correction for VLM, selecting all the available Mediterranean stations which covered at least 50% of the simulated period. In order to avoid relative biases due to different zero-level assumptions, monthly data from the model grid point nearest to the tide gauge have been compared with the monthly mean tide gauge data after subtraction of the respective mean values, over the common time interval, from both modelled and tide gauge time series. The 34 locations for which the comparison was carried out are shown in Fig. 12.
The cross-correlation coefficients obtained are detailed in Table 5, where the geographic location of the stations  (17,23,25,26,28,29,31,34) located in the Aegean Sea or nearby, and one in the Black Sea. This is quite encouraging, considering that the model values are not at the stations' sites, which are very close to the coast, and that the model freely evolves for three decades. It shows that, in part of the basin, the model is capable of providing a fairly accurate description of the seal level variability on long time scales, not only over large spatial scales, but also nearshore, where the consequences of sea level change are stronger. The fact that the model performs less well in the Aegean Sea is likely partly due to difficulties in reproducing small scale details of the near coast dynamics, induced by the very complex bathymetry of the area. Sea surface height time series for six stations (red circles in Fig. 12) where VLM can be considered to be negligible are plotted in Fig. 13. The main interannual variability is reproduced fairly well at all stations, and especially at Malaga, Trieste, Preveza and Alexandria.

Future basin hydrology and circulation
For the RCP8.5 scenario, Fig. 14 shows the differences in the future surface temperature climatologies with respect to the reference period. Averages are computed over years 2046-2065 (left) and 2081-2100 (right), for winter (top) and summer (bottom). The sustained increase in temperature is evident in both seasons. As for the correspondent salinity fields (Fig. 15), the salinification of the Eastern basin is apparent, especially in the north, as well as the northward displacement of the Atlantic stream in the Western basin, yet accompanied by a more efficient penetration of the AW  The temperature-induced additional buoyancy is evident in the mixed layer depth which exhibits significantly lower relative maxima at dense water formation sites. The deep convection in the Gulf of Lion is almost neutralised, whereas in the Adriatic Sea is characterised by enhanced variability (Fig. 16).
Basin-wide averages of T and S are summarized in Table 6 Atlantic Ocean (https:// icdc. cen. uni-hambu rg. de/ en/ ar5-slr. html, Slangen et al. 2014). Prior to spatial averaging, all the variables have been interpolated onto the MED16 grid via a bilinear interpolation, with nearest-neighbor interpolation near the coasts, results from the MED16 regional projection are also shown. The zero level is set at the 1986-2005 mean that Slangen et al. (2014) consider as the common zero SSH anomaly value.
The most relevant contribution is that from the sterodynamic term (dark red; solid line = GCM, dashed line = MED16, shaded area = GCM 90% spread), which accounts for changes in the circulation and density distribution of sea water and for the inverse barometer (IB) correction associated with atmospheric pressure patterns (Gregory et al. 2019). Positive changes (here listed in decreasing order) also derive from distributed glacier melting (aside from Greenland and Antarctica ice sheets-green), Antarctic dynamic ice loss (pink), Greenland surface mass balance (light blue), Greenland dynamic ice loss (blue) and glacial isostatic adjustment (red), ground water storage (yellow), while the Antarctic surface mass balance (dark green) determines a reduction in sea level. At the end of the XXI century, these components are projected to respectively account for 51%, 23%, 13%, 10%, 6%, 2.5%, 2%, − 7% of the total, with an overall contribution of 16% from Greenland and 6% from  Antarctica. Recent works (Spada, 2017;Santamaría-Gómez et al. 2017) have dwelled on the necessity to consistently include the GIA effect in updated estimates of SLC in the Mediterranean, based on the alternative datasets provided by Lambeck et al. (1998) and Peltier et al. (2004), whose limits and reliability have been the subject of a lively scientific discussion (Peltier 2002;Lambeck et al. 2002). Here, the GIA contribution is the one computed by Church et al. (2013b) by averaging results from the ICE-5G model (Peltier et al. 2004) and the ANU ice sheet model (Lambeck et al. 1998 and subsequent improvements), using the updated SELEN code for the sea level equation (Spada andStocchi 2006, 2007). In the Mediterranean basin, however, despite its comparative importance for the assessment of relative sea level rise for present times, the GIA contribution appears to progressively lose relevance by the end of the XXI century, unless ongoing research revises the long-term ranking of individual contributions. Although only representing one of the numerous possible realizations that can be obtained by different GCM-RCM modelling chains, the MSL trend estimated via the MED16 scenario simulation is very similar to that of the GCM ensemble mean from Slangen et al. (2014), although it slightly accelerates from 2040 onwards, veering towards the high-end global scenarios. With respect to the coarse-resolution global average, the MED16 downscaling better characterizes the regional patterns of SLC arising from ocean dynamics, and it consequently retains enhanced variability.
The time evolution of the difference in SSH between the Mediterranean basin interior and the neighbouring Atlantic Ocean (average over the region just outside the Gibraltar Strait, from 9 to 6 W) is reported in Fig. 18, for both the MED16 model (left panel) and the CMIP5 multi-model mean (right panel). The red curves correspond to the sterodynamic term (i.e., to the dark red curves of Fig. 17), the yellow curve in the left panel represents the historical reference experiment, while the AVISO observations are plotted in black. Beside the good agreement between the AVISO data and the historical experiment, it is worth noting that (Fig. 19) the regional simulation is capable of capturing the hydraulic jump across the Gibraltar Strait (Sannino et al. 2009a(Sannino et al. , b, 2015 that is not accounted for in the global time series, which only describes the evolution of anomalies with respect to the reference period 1986-2005. Both experiments exhibit the relative "sinking" of the Mediterranean with respect to the Atlantic, although the global model follows a trend that is half that of the regional experiment, and unsurprisingly characterized by lower amplitude and frequency variability. Some coherence is yet detectable between the respective typical time-scales of variability, possibly when driven by large-scale patterns, as well as a comparable signal-to-noise ratio. However, the global projection seems to be reaching a new stable state towards the end of the XXI century that is not discernible in the regional experiment, a longer run being needed to confirm such a hypothesis. So far, the change in the relative level between the two basins can be estimated to be within centimetres, over an expected sea level increase on the order of tens of centimetres by the end of the century. Figure 20 illustrates the spatial distribution of the progressive increase in sea level under RCP8.5 with respect to the historical period 1985-2005, which qualitatively retains the overall pattern of the dynamic topography under current climate, where the largest increments apparently correspond to the Atlantic signal, mirroring the AW progress Fig. 12 Location of the tide gauge stations whose data are used to construct Table 5. Red circles indicate stations used to extract the time series displayed in Fig. 13 into the basin interior. The main notable difference appears to be the sustained amplified signal around the Balearic Islands, presumably itself a consequence of the northward displacement of the Atlantic inflow that numerical models associate with the newly detected Western Mid-Mediterranean Current (WMMC). This current corresponds to the northward meandering stream arising from the bifurcation of the Atlantic inflow after exiting the Alboran Sea and overcoming the Almera-Oran front (Arnone et al. 1990;Pinardi et al. 2015). The WMMC branch of the Atlantic inflow has also been detected in the hindcast and historical experiments (Fig. 6), yet it has never been well documented  in observational datasets. Pinardi et al. (2015) tentatively interpret it as a residual current resulting from multi-year averaging. Figure 21 reports the differences in SLC between MED16 and the GCM ensemble mean for the two targeted time horizons. The overall pattern is maintained in time, with localized yet roughly constant negative biases (i.e., CGM estimates are higher) in the Levantine basin, and positive biases in the Ionian Sea and the Western basin, generally not exceeding 10 ÷ 12 cm, and an increasingly more pessimistic regional projection for the Balearic region, going from a deviation of 15 cm÷17 at mid-century, to one of 17 ÷ 22 cm in 2100.

Conclusions
We have here analysed the results of three climatic simulations of the Mediterranean Sea circulation and sea level, performed with MED16, an updated version of the tide-including, climatic ocean model presented in Sannino et al. (2015). These consist of a hindcast simulation, covering three recent decades (1981-2010), a historical run , and a The hindcast simulation has been used to assess the ability of the new model to reproduce the past evolution of the basin circulation, with excellent results. This simulation has a value in itself, because it is the first pluri-decadal, eddyresolving simulation of the Mediterranean Sea circulation that includes the main tidal effects, together with an accurate treatment of the complex dynamics occurring in the Gibraltar Strait area. We have therefore preliminarily analysed the simulated transports at the Gibraltar, Dardanelles and Bosphorus straits, and the basin hydrology and circulation.
We found that the model correctly represents the transport evolution at Gibraltar, the average circulation in the surface and intermediate layers, and the main variability of the Mediterranean hydrologic structure during the simulation period. The simulated SLA was found in good agreement with satellite observations, exhibiting a positive trend in the Western Mediterranean Sea, and a more complex behavior in the Levantine basin, where sea level variability has been strongly influenced by the EMT and its evolution.
The spatial structure of the hindcast average SSH field is well reproduced over the entire basin, especially in the western Mediterranean where results are much closer to that of the MDT, due to the fact that the use of improved sea level information at the Atlantic lateral boundary significantly enhances the reliability of results. A further improvement has been obtained through the adequate treatment of the complex, hydraulically driven dynamics across the Gibraltar Strait. The evolution of the hindcast SSH near coast compares very well with tide gauge observation along most of the basin's coasts. This shows that MED16, if adequately forced, can provide useful information on the sea level variability in coastal areas, which are the more exposed to SLR, even though the current horizontal resolution of the model does not allow to exactly resolve the locations of the tide  The historical simulation, which sets the reference baseline for the future scenario projections, also correctly reproduces the main mean features of the Mediterranean circulation and hydrology, a result that is all the more satisfactory when considering that no specific constraint or relaxation is prescribed.
Under the RCP8.5 future scenario, the temperature is projected to generally increase. The spatial pattern of variations in surface salinity is affected by the penetration of the Atlantic stream into the basin interior, with the westernmost regions exhibiting appreciable decreases. The warming of sea waters results in the inhibition or enhanced convection variability in the main deep-water formation sites.
As regards future sea level, the overall DSL trend estimated by MED16 is comparable to the ensemble mean computed from the Global Circulation Models analysed in Slangen et al. (2014). Nevertheless, our regional dynamical downscaling provides a better characterization of regional patterns arising from ocean dynamics, and of their variability, and gives a more pronounced relative sinking of the Mediterranean with respect to the Atlantic. The regional differences in SLR between the presented downscaling and the GCM mean projection are generally constant in time, yet an acceleration emerges in the regional projection for the Balearic area.
Indeed, state-of-the-art numerical estimates show that future sea level rise in the Mediterranean basin will crucially depend on the evolution of the sterodynamic  This also provides an indirect check of the general quality of the atmospheric forcing from the SMHI downscaling. This suggests that further improvement of the forcing may be needed to obtain an even more realistic description of the circulation structure. Fig. 19 Evolution of salinity in a section of the Strait of Gibraltar at intervals of 1 h component, as determined by the propagation of the Atlantic signal through the Strait of Gibraltar and by the local circulation. MED16 effectively provides an accurate and reliable representation of the basin's DSL, by accounting for the numerous and complex processes that concur to determine it, including explicit tidal forcing and an accurate resolution of the Gibraltar Strait. Nevertheless, in view of the significant numerical resources demanded, the assessment is recommended of the minimum requirements models must comply with, in terms of resolution and explicit process representation.  Slangen et al. (2014). Mean over the period 2046-2065 (upper panel), and over the period 2080-2099 (lower panel) observational products-in this case, RIXEN from Rixen et al. (2005) and EN 4 from Good et al. (2013)-which for the period of interest reach up to ~ 0.1 °C, compared to an associated uncertainty that can be twice as large.
The distance between the reanalysis and their partly constraining reference observations over the period 1987-2021 is reported in Table 8, as derived from the Quality Information Document for the Copernicus Reanalysis (https:// catal ogue. marine. coper nicus. eu/ docum ents/ QUID/ CMEMS-MED-QUID-006-004. pdf). Results for the deepest layer suffer from scarcity of data and more pronounced model error. For the sake of completeness, we recall here that MEDHYMAP data are derived via a variational interpolation of observation, as the RIXEN data are, while the EN 4 dataset is the result of Optimal Interpolation. The Copernicus Reanalysis is computed by assimilating observations in the operative NEMO model runs via a 3DVar scheme, and therefore represents their dynamical interpolation.