River runoff influences on the Central Mediterranean overturning circulation

The role of riverine freshwater inflow on the Central Mediterranean Overturning Circulation (CMOC) was studied using a high-resolution ocean model with a complete distribution of rivers in the Adriatic and Ionian catchment areas. The impact of river runoff on the Adriatic and Ionian Sea basins was assessed by a twin experiment, with and without runoff, from 1999 to 2012. This study tries to show the connection between the Adriatic as a marginal sea containing the downwelling branch of the anti-estuarine CMOC and the large runoff occurring there. It is found that the multiannual CMOC is a persistent anti-estuarine structure with secondary estuarine cells that strengthen in years of large realistic river runoff. The CMOC is demonstrated to be controlled by wind forcing at least as much as by buoyancy fluxes. It is found that river runoff affects the CMOC strength, enhancing the amplitude of the secondary estuarine cells and reducing the intensity of the dominant anti-estuarine cell. A large river runoff can produce a positive buoyancy flux without switching off the antiestuarine CMOC cell, but a particularly low heat flux and wind work with normal river runoff can reverse it. Overall by comparing experiments with, without and with unrealistically augmented runoff we demonstrate that rivers affect the CMOC strength but they can never represent its dominant forcing mechanism and the potential role of river runoff has to be considered jointly with wind work and heat flux, as they largely contribute to the energy budget of the basin. Looking at the downwelling branch of the CMOC in the Adriatic basin, rivers are demonstrated to locally reduce the volume of Adriatic dense water formed in the Southern Adriatic Sea as a result of increased water stratification. The spreading of the Adriatic dense water into the Ionian abyss is affected as well: dense waters overflowing the Otranto Strait are less dense in a realistic runoff regime, with respect to no runoff experiment, and confined to a narrower band against the Italian shelf with less lateral spreading toward the Ionian Sea center.


Introduction
The annual mean freshwater budget in the Mediterranean Sea, composed of evaporation minus precipitation and river runoff, has been found to be positive, corresponding to a net surface loss of 0.54 ± 0.15 m/year as the average of several long term annual estimates proposed by Sanchez-Gomez et al. (2011). Moreover, the net heat budget of the basin is well known to be negative according to most estimates (Bethoux 1979;Pettenuzzo et al. 2010;Sanchez-Gomez et al. 2011;Schroeder et al. 2012). The combination of the positive freshwater budget and the negative heat budget of the Mediterranean Sea gives a net buoyancy loss, thus sustaining a vigorous mean kinetic energy in the basin (Cessi et al. 2014).
This forcing induces also a vigorous meridional antiestuarine overturning circulation. Pickard and Emery (1982) differentiate between the semienclosed seas with horizontally and vertically separated inflow and outflow regimes at the strait. Semienclosed basins with strait flows that Abstract The role of riverine freshwater inflow on the Central Mediterranean Overturning Circulation (CMOC) was studied using a high-resolution ocean model with a complete distribution of rivers in the Adriatic and Ionian catchment areas. The impact of river runoff on the Adriatic and Ionian Sea basins was assessed by a twin experiment, with and without runoff, from 1999 to 2012. This study tries to show the connection between the Adriatic as a marginal sea containing the downwelling branch of the anti-estuarine CMOC and the large runoff occurring there. It is found that the multiannual CMOC is a persistent anti-estuarine structure with secondary estuarine cells that strengthen in years of large realistic river runoff. The CMOC is demonstrated to be controlled by wind forcing at least as much as by buoyancy fluxes. It is found that river runoff affects the CMOC strength, enhancing the amplitude of the secondary estuarine cells and reducing the intensity of the dominant anti-estuarine cell. A large river runoff can produce a positive buoyancy flux without switching off the antiestuarine CMOC cell, but a particularly low heat flux and wind work with normal river runoff can reverse it. Overall by comparing experiments with, without and with unrealistically augmented runoff we demonstrate that rivers affect the CMOC strength but they can never represent its are vertically separated can have two vertical circulation modes, estuarine and antiestuarine.
Traditionally, the estuarine and anti-estuarine circulation in semienclosed seas has been classified on the basis of the surface water flux and the conservation of salt (Knudsen 1900;Sverdrup 1947;Pickard and Emery 1982). More precisely, estuarine and anti-estuarine circulations can then be represented by meridional overturning cells which form into the semi-enclosed basin and which are connected to the external regions by the Strait inflow/outflow system.
The general characteristics of the Mediterranean overturning circulation are schematized in Fig. 1  ). This circulation is characterized by interannual as well as multi-decadal time scales and it is composed of three major conveyor belts: the Zonal Overturning Circulation (ZOC) in the Southern Mediterranean propelled by the Gibraltar stream flow and Levantine Intermediate Water (LIW) formation processes, the Western Mediterranean Meridional Overturning Circulation (WMOC) originating in the Gulf of Lion, the Central Mediterranean MOC originating in the Adriatic Sea (CMOC). These overturning cells are triggered by buoyancy losses and water mass sinking occurring in the regions offshore the Gulf of Lions, in the Southern Adriatic and in the Northern Levantine Basin, enhanced by the presence of large scale permanent cyclonic gyres driven by wind stress curl . The Aegean Sea is marked as another dense water site in Fig. 1: Roether et al. (1996) and Gertman et al. (2006) demonstrated that dense waters could form there during intense winters and large LIW inflow from the Levantine basin. This Aegean Meridional Overturning Circulation (AeMOC) characterizes the Eastern Mediterranean Transient (EMT) which started at the end of the eighties and ended around the middle of the nineties. Gacic et al. (2011) pointed out that the EMT could be a repeating phenomenon, and part of its cause is connected with inversions of the Ionian circulation.
In this paper we would like to exhamine for the first time the question of what is the influence of river runoff on the CMOC. The question is relevant because more than 30% of the Mediterranean runoff is concentrated in the Adriatic Sea and the hypothesis that large freshwater runoff in the Adriatic could shutdown the CMOC is realistic at least for the past. As far as we know this is the first study on the river influence on the CMOC, similar to Rahmstorf's study (Rahmstorf 1995) on the freshwater role in the Northern Atlantic overturning circulation. Previously, Skliris et al. (2007) and Vervatis et al. (2013) investigated the impact of reduced discharge of Ebro and Nile rivers on the dense water formation in the Eastern and Western Mediterranean sub-basins.
More generally, the role of freshwater inputs (due to both rivers and precipitations) on the dynamics of a marginal sea, such as the Adriatic basin, has been widely investigated in the literature. Spall (2012) demonstrates that, in marginal sea areas, an increase in surface freshwater gain can lead to a shutdown of dense water formation and sinking, and the marginal sea MOC switches from antiestuarine to estuarine mode. Recently, Cessi et al. (2014) established that the estuarine/anti-estuarine character of a semi-enclosed sea with a two-layer flow at the strait is determined by both wind and buoyancy forcings. The wind forcing is normally a source of mechanical energy for the circulation, while the buoyancy forcing could be either an energy source or a sink depending on the sign of the net buoyancy flux at the surface. For estuarine basins, such as the Baltic and Black Sea, the positive buoyancy flux (dominated by precipitation and runoff exceeding evaporation), is a net energy sink for the circulation, thus producing a less Fig. 1 The conveyor belts of the Mediterranean Sea. The red and yellow dashed streamlines in the zonal direction stand for the zonal overturning circulation in the surface-intermediate layers that is forced by the Gibraltar stream flow and Levantine Intermediate Water (LIW) formation processes. The red spirals indicate the preferential sites for strong heat losses during wintertime and dense water formation processes. Two anti-cyclonic meridional overturning circulation patterns can be distinguished (white spirals): the Western Mediterranean MOC originating in the Gulf of Lion, and the Central Mediterranean MOC originating in the Adriatic Sea. Reproduced from Pinardi et al. (2006) vigorous meridional circulation than in the anti-estuarine basins.
The Adriatic Sea is a "dilution basin" with a negative annual freshwater budget, equal to the difference between evaporation and precipitation and runoff, of about −1 m year −1 (Artegiani et al. 1997a, b) but a negative annual mean heat budget. The buoyancy flux could eventually then be positive, thus determining a net sink of energy and possibly a net estuarine character of the circulation. Pinardi et al. (2006) show that, due to river runoff, the Adriatic Sea could be characterized by zero net buoyancy flux, thus producing a basin where the circulation is mainly powered by the wind stress. However, the energetics proposed by Cessi et al. (2014) cannot be applied satisfactorily to the Adriatic Sea, since the flow at Otranto is not just a two-layer flow. Thus, a comprehensive analysis of surface buoyancy and the CMOC is needed to fully establish the estuarine/anti-estuarine character of the Adriatic Sea circulation.
In order to answer the effects of runoff on the CMOC and the Adriatic Sea deep water formation mechanisms, a high resolution general circulation model was set up, forced by realistic fluxes of water, heat and momentum. The impact of river runoff on the circulation is assessed by using a mechanistic approach: a full forcing and dynamics experiment is compared to a zero-runoff case. The full forcing experiment is considered to be the present state of the Adriatic Sea and it is validated with existing data sets.
The paper is organized as follows. Section 2 outlines the experimental design, and details the parameterization of the rivers as surface boundary conditions and the validation of the model performance. Section 3 describes the experiments highlighting the role of the river runoff forcing. A summary and conclusions are presented in the last section.

Model configuration and experiments set-up
The numerical model used is the Nucleus for European Modelling of the Ocean, NEMO (Madec 2008), that is a three-dimensional finite difference numerical model adopting the Boussinesq and hydrostatic approximations and using the linear explicit free surface formulation. The area covered by the model grid is the Central Mediterranean Sea from 12.2 to 21.0°E and 30.2 to 45.8°N with a horizontal resolution of about 2.2 km (2.5 km in the meridional direction and 1.7-2.2 km in the zonal direction). Figure 2 shows the bathymetry of the domain, river mouth grid points and the three sub-regions into which the Adriatic Sea is conventionally subdivided on the basis of its bottom morphology: the Northern (NAd), the Middle (MAd) and the Southern Adriatic (SAd). The connection with the Ionian Sea occurs at the Otranto Strait where the sill is 800 m deep, located at approximately 40°N.
This model is the first implementation of the NEMO code over the Central Mediterranean Sea area, with a representation of almost all the rivers flowing into the Central Mediterranean Sea. In this work the regional model is forced by a 1/16° resolution daily analyses from the operational Mediterranean forecasting system, MFS (Tonani et al. 2008;Pinardi and Coppini 2010) at the lateral open boundaries. The numerical model configuration is explained in details in "Appendix 1".
A set of experiments were performed with present day river runoff, without river runoff and with augmented runoff (Table 1), spanning the period from 1 January 1999 to 31 December 2012. The time series of the kinetic energy integrated over the basin volume (not shown) shows that the spin up period consists of the first few months of 1999, thus this year has been eliminated from the analysis. Such short spin up is due to the fact that the model was initialized by close-to-present-day fields from MFS model.
It is worth to stress that the conceptual paradigm of this study is to explain the "theoretical role" of river inflow since the case without or with increased runoff corresponds to extreme conditions with respect to the present day. EXP1 corresponds as closely as possible to reality and this is our reference against which we analyze the effects only of river runoff changes, leaving all the other forcings identical because we know that the estuarine water balance of the Adriatic is due to river runoff.

River runoff datasets
River runoff into the Mediterranean Sea is mainly concentrated in the Central Mediterranean sub-basin, with rivers flowing into the Adriatic Sea providing almost 1/3 of the total runoff (Struglia et al. 2004;Ludwig et al. 2009). Overall the Mediterranean Sea counts on a great number of very small rivers (Milliman 2001), owing to strong topographic relief favouring the formation of small watersheds. The construction of dams in the Nile and Ebro have reduced their runoff in the Mediterranean (Nile from 2700 to 150 m 3 s −1 , Ebro from 1500 to 400 m 3 s −1 ) so that currently the only major runoff sources located out of the Central Mediterranean Sea are the Rhone and the Ebro rivers flowing into the Western sub-basin.
The freshwater discharge into the Central Mediterranean Sea is almost totally concentrated along the Adriatic coastlines: the Po river is the main freshwater source of the Many of the records were downloaded by Mediterranean Hydrological Cycle Observing System, Medhycos, data server (Medhycos 2001), other from the Global River Discharge database RivDis (Vörösmarty et al. 1996), or from the Hydro database at the French Ministry of Environment (Hydro 2006). Struglia et al. (2004) provides an extensive representation of historical river release into the Mediterranean basin, timeseries of 67 rivers are derived from the Global Runoff Data Center (GRDC) hydrological database and the Medhycos database.
Most of the datasets we adopted consist of observations taken at the hydrometric stations nearest to river mouths, and few of them are estimated values (Pasaric 2004;Malacic and Petelin 2009;Raicich 1996). Time series for the various rivers cover different periods. However, the time series of the major rivers, accounting for the most of the Central Mediterranean discharge, overlap at least for 20 years. The databases we considered have been provided by the GRDC, the regional agency Autorita' di Bacino Basilicata, the hydrological division of the National Research Council CNR IRPI, the regional agencies for the environmental protection ARPAs, the Albanian Hydrometeorological Institute.
Overall this study considers 67 Adriatic and Ionian rivers in total, 52 flowing into the Adriatic Sea and 15 into the Ionian Sea. Table 2 lists the adopted climatological datasets for river runoff, the time range for computing the monthly climatologies and the annual mean discharges as useful reference values.
All river mouths are "point sources" except for 2: Marecchia to Tronto rivers (Tronto excluded) in the Marche region and Vibrata to Fortore rivers (Fortore excluded) in the Abruzzo and Molise regions which are "diffused sources", i.e. the runoff was split across several grid points. These diffused sources and the rivers in Puglia are the only rivers of the model based on Raicich's (1996) climatologies (see Fig. 2 for locating river outlets).
The Po river runoff, referred to the Pontelagoscuro station, usptream of the delta, is unequally subdivided between the nine grid points representing the nine branches of the delta (Po di Goro, Po di Gnocca, Po di Tolle, Po di Bastimento, Po di Scirocco, Po di Bonifazi, Po di Dritta, Po di Tramontana, Po di Maistra) according to percentages in Provini et al. (1992).
Daily time series of total discharge in the model domain during the simulation period (1 January 1999-31 December 2012), are shown in the top panel of Fig. 3, with overlapped Po daily discharge as well as the contribution of all the other rivers (i.e. monthly climatologies interpolated on daily basis). Middle panel highlights the added value of working with daily observations of Po river instead of considering monthly climatologies. Finally the bottom panel shows how the Adriatic river release, Po river excluded, is represented by Raicich's climatologies (1996), Ludwig's climatologies (Ludwig et al. 2009) and the new climatologies we chose in this study on the basis of more updated and reliable datasets. To note that our dataset includes 52 rivers flowing into the Adriatic Sea, while Raicich (1996) considers only 45 and Ludwig et al. (2009) only 32 Adriatic rivers.
Adding all the contributions, the annual average runoff in the Central Mediterranean Sea is equal to 4.72 × 10 3 m 3 s −1 , 29.7% coming from the Po river only, 94.6% is the contribution of all the Adriatic rivers and 5.  Mean annual discharge augmented of 50% with respect to Exp1 EXP 4 Mean annual discharge augmented of 100% with respect to Exp1 (5.19 × 10 3 m 3 s −1 ), 2009 (5.37 × 10 3 m 3 s −1 ), and 2010 (5.26 × 10 3 m 3 s −1 ). The daily and annual maximum values are due to Po river only, as the other rivers consist of monthly climatologies. As detailed above, the Po river represents the largest contribution to the total discharge of the Central Mediterranean Sea, however we should consider we are missing the daily and interrannual variability coming from the other rivers. The interannual variability described above is the result of only Po river runoff interannual differences because daily data sets for other rivers are not available for such an extended time period.

River parameterization in the model
Model rivers are parameterized as "surface sources" of water at the estuary border grid points while no temperature information is prescribed. Our assumption of no temperature differences between river inflow and the marine environment is generally valid as river plumes are controlled by the salinity gradients.
The runoff and salinity values are prescribed at river outlets in the vertical velocity and salt flux boundary conditions respectively. Here we follow the natural boundary condition approach (Huang 1993) plus ad-hoc salt values at river outlets, this is the most consistent before rivers will be considered with lateral open boundary conditions. Such boundary condition involves both a water and salt flux condition, as written in Eqs. (20) and (21), in order to conserve the volume integrated salt content in the basin. It can be easly demonstrated that by taking the volume integral of advection/diffusion equation for salinity and by replacing the boundary conditions (20) and (21) we obtain the salt conservation regardless of the constant value of salinity chosen at river mouths. Kourafalou et al. (1996) implemented for the first time the natural boundary condition to the riverine freshwater flux, followed by Skliris et al. (2007), Beuvier et al. (2010), Dell'Aquila et al. (2012 and Vervatis et al. (2013). In addition operational models of the Mediterranean Sea already use the natural boundary condition from several years (Oddo et al. 2009).
Following Simoncelli et al. (2011) we prescribe constant salinity values at all river mouths parametrizing the effects of tidal mixing inside the river estuaries. Values chosen are equal to 15 psu for all rivers, except for the Po river where 17 psu is considered. These constant salinity values are the result of sensitivity tests performed on the basis of salinity profiles measured at river mouths (Simoncelli et al. 2011) and at the center of the basin (Oddo et al. 2005).

Model validation with observations
The performance of EXP1 was evaluated by comparing the simulated fields with available in situ observations. For the in situ data, two Argo profiling floats in the Adriatic and Ionian Seas (Fig. 4) were used to calculate vertical mean profiles Jan99 Jan00 Jan01 Jan02 Jan03 Jan04 Jan05 Jan06 Jan07 Jan08 Jan09 Jan10 Jan11 Jan12 of RMSE and BIAS for temperature and salinity (Fig. 5).
The profiles are in the area of deep water formation and in the northern Ionian exit, so in the most relevant places for our simulations. On the whole, EXP1 is in good agreement with the observed data, comparable RMSE and BIAS values were found for different models of the Adriatic Sea (Guarnieri et al. 2013) and the Ionian Sea (Federico et al. 2016;Kassis et al. 2016). In EXP2 without rivers, water masses in the Adriatic Sea appear saltier and warmer with respect to the observed dataset. EXP2 has then doubled the RMSE and BIAS errors with respect to observations thus showing already a large impact of rivers on the estimate of the circulation in the Adriatic Sea. The mean profiles of RMSE and BIAS do not show significant differences in the Ionian Sea instead. Basin-averaged timeseries of analysed and modelled sea surface temperature are shown in Fig. 6. Analysed SST is obtained by Optimal Interpolation of Pathfinder AVHRR (Advanced Very High Resolution Radiometer) SST data (Pisano et al. 2016). Model results are characterized by a positive BIAS but are capable to reproduce the observed interannual variability. The agreement between model and analyses timeseries is high, with correlations of about 99.6% and RMSE is 0.8 °C for EXP1.

Is the Adriatic Sea an estuarine or anti-estuarine basin?
The literature contains ample evidence that the MOC is primarily driven by wind and tidal stirring (Munk and Wunsch 1998;Paparella and Young 2002;Marshall and Speer 2012). In addition, the relationship between the dense water formation, driven by buoyancy fluxes, and the strength of the overturning circulation has been highlighted in several theoretical, as well as realistic modelling studies (Rahmstorf 1995(Rahmstorf , 1996Pisacane et al. 2006). Similarly, here we focus on the downwelling branch of the Central Mediterranean MOC, which develops inside the Southern Adriatic sub-basin due to the local openocean convection and dense water formation sustained by winter heat losses and a cyclonic gyre driven by wind stress curl. In this section we will describe both the surface forcing and the newest theoretical framework which connects the anti-estuarine and estuarine character of a marginal sea to the vertical overturning circulation.

The surface forcing
In order to assess the Adriatic Sea overturning circulation on the basis of surface forcing, an analysis of both buoyancy and wind stress fluxes was performed. The surface buoyancy flux per unit area, Q B (m 2 s −3 ) is calculated according to Cessi et al. (2014) as follows: where T,S are the coefficients of thermal and haline expansion respectively, 0w is the reference sea surface water density, Q is the net heat flux, C w is the heat capacity of sea water, S 0 is the surface salinity. Finally, (E − P − R∕A). is the freshwater flux with evaporation rate, E, and precipitation rate P, in m s −1 , river discharge R in m 3 s −1 and A representing the grid area of river mouths (m 2 ). Furthermore the following values were assumed: The net heat flux, Q, components are computed according to bulk formulae described in "Appendix 1".
Our results show the the whole Central Mediterranean Sea has a positive freshwater budget, 0.60 m year −1 , while the Adriatic Sea is negative, −0.69 m year −1 , as expected.
The annual time series of Adriatic surface buoyancy flux and wind work is shown in Fig. 7. The realistic  occurred also in 2009 (see Fig. 3) but this year is characterised also by a intense wind work (Fig. 7, bottom panel) which makes river influence on the basin boyancy budget less relevant.
The wind work (m 3 s −3 ) is defined as T w ⋅ u s 0 where u s is the sea surface velocity, 0 is the reference sea surface water density and T w is the wind stress defined in "Appendix 1". The heat, freshwater and wind forcings acting at the air-sea interface may be directly compared by considering the buoyacy flux and the wind work: both quantities represent source/sink terms of energy. In order to compare the surface integrated buoyancy flux with the surface integrated wind work, the latter is integrated over the basin area and divided by the basin volume. Figure 7 shows that the wind work is always positive (10 −8 m 2 s −3 ), implying a net of mechanical energy for the Adriatic Sea that is one order of magnitude larger than the buoyancy flux (10 −9 m 2 s −3 ) in EXP1. Thus the Adriatic Sea results again to be an anti-estuarine basin characterized by a large wind work energy source. Spall (2011Spall ( , 2012 analysed a marginal sea overturning circulation predicting that the freshwater input could stop the anti-estuarine circulation of a theoretical basin. Rahmstorf (1995) speculates that increasing freshwater inflow in the Northern Atlantic may potentially reduce or even shut down the overturning circulation. Cessi et al. (2014) show that both buoyancy forcing and wind stress work are connected to the strength of the circulation and thus, also to the MOC developing between the marginal sea and the open ocean.

Theoretical marginal sea overturning estimates
Following Spall (2012) we assessed the non-dimensional thermal forcing parameter, µ/ε, and freshwater forcing parameter, γ/ε. The former describes the relative balance between heat budget in the interior basin and the lateral eddy fluxes from the opening advecting warm water into the basin. The latter describes the relative balance between freshwater budget in the interior basin and the lateral salt eddy fluxes at the opening again. In the Adriatic Sea lateral eddy fluxes are advecting heat and salt through Otranto thus opposing the net heat loss and water gain (see previous sections): thus it is important to evaluate these nondimensional numbers to know how the lateral open boundary forcings work against the surface buoyancy flux. Small values of these parameters with eventually negative values of the freshwater parameter indicate that the lateral prevail over the surface cooling fliuxes in the basin interior and may trigger the shutdown of deep convection and then the reversal of the MOC.
We computed the two parameters over the whole simulation period and discovered that the thermal parameter is essentially the same in EXP1 and EXP2 experiments, 5 × 10 −5 in EXP2 and 4.9 × 10 −5 in EXP1, while the freshwater parameter is 7 × 10 −4 in EXP2 and −2 × 10 −2 in EXP1 ("Appendix 2" for details on the computation of Spall's coefficients). This means that the Adriatic Sea runoff has the potential to shut down the deep convection and reduce the intensity of the antiestuarine CMOC. Spall (2012) also shows that the non-dimensional ratio ΔS∕ΔT, where ΔT and ΔS are respectively the model-diagnosed temperature and salinity non-dimensional anomalies between the interior basin and the open boundary currents, can be written as a function of the thermal and freshwater forcing parameters (see "Appendix 2" for details on the computation of these values). A ratio ΔS∕ΔT less than 1 means that the general circulation is in "thermal mode", which means that the heat forcing prevails and an antiestuarine MOC develops. A ratio ΔS∕ΔT > 1 indicates the "haline mode" of the marginal sea circulation with the shut down of deep convection and the possibility of an estuarine MOC if the freshwater budget is negative. The collapse of deep convection is demonstrated to be possible also in the thermal mode case if ΔS ΔT > 0.5. We obtain ΔS∕ΔT = 0.28 and 0.12 in EXP1 and EXP2 respectively. Thus the Adriatic Sea is characterized by an anti-estuarine circulation with thermally driven deep water formation processes despite a large runoff budget. Focusing on the year 2002 which showed the largest river runoffs (Fig. 3), ΔS∕ΔT = 0.42 in EXP1 and 0.29 in EXP2, meaning that the Adriatic deep water formation and the antiestuarine circulation characterize both experiments with and without river runoff but EXP1 is close to the collapse of deep convection.

How is the intensity of Central Mediterranean MOC affected by runoff?
In order to better quantify the river influence on the CMOC, an inter-annual analysis of the meridional transport stream function was carried out. The meridional transport stream-function, Ψ, is defined as (Pedlosky 1987): with −H < z < 0 as the depth, x 0 and x 1 the more eastern and more western sea points, v is the time-averaged meridional velocity. The velocity field is tangent to the isopleths of Ψ, and positive Ψ values indicate anti-estuarine cells turning anti-cyclonically or clockwise, while negative values indicate estuarine cells turning cyclonically or anticlockwise.
Top panels of Fig. 8 show the transport stream-function for EXP1 and EXP2, averaged over the whole simulation period. A large anti-estuarine cell down to a 700-800 m depth is detected in both experiments in the Northern Ionian Sea and SAd sub-region, but with different intensities. Interestingly enough many estuarine cells exists in the domain: one at the surface in the NAd, one in the deep layers of the Southern Adriatic Pit around 41-42ºN, SAP, another in the Middle Adriatic Pit around 43ºN, MAP, and the last in the Northern Ionian abyss. The estuarine cells attached to the seabed represent a more stagnant circulation in the Adriatic Pits and Ionian abyss. To note that the estuarine cell which extends to the whole depth of the NAd is the consequence of a strong mixing of the entire water column (i.e. about 35 m depth) due to heat losses and wind stress. In the Ionian abyss the deep estuarine cell attached to the seabed below the CMOC is linked to the SAd dense waters overflowing the Otranto Strait and stratifying below the CMOC cell (Curchitser et al. 2001;Bensi et al. 2013a). Deep estuarine cells exist below the wind driven antiestuarine MOC cells of the mid-to-deep water depths of the global ocean as demonstrated by Vallis (2011) andDe Lavergne et al. (2016). Recently Ferrari et al. (2016) found that a deep overturning cell exists below the antiestuarine Atlantic MOC. This is due to the turbulent boundary layer very close to the sloping bottom of the deep abyssal plains, generating diapycnal upwelling below the Atlantic MOC. In the Mediterranean we can speculate that the same mechanism is present due to the complex topography and the slow abyssal cyclonic motion of deep water masses (Curchitser et al. 2001).
Top panels of Fig. 8 demonstrate that with both a zerorunoff and a realistic runoff case, the anti-estuarine character of the CMOC prevails. Indeed, in EXP1, the secondary estuarine cells of the NAd, MAd, SAd sub-regions and Northern Ionian basin are larger than in EXP2, however the anti-estuarine MOC cell still dominates.
The estuarine component of the MOC may become more evident on a seasonal basis, particularly during summer. Figure 9 focuses on summer 2002 and summer 2009 because these years had the largest spring river discharge (Fig. 3). In summer 2009 (Fig. 9, bottom panels), a welldefined surface estuarine cell characterizes the whole meridional extension of the basin with no differences between EXP1 and EXP2, which means that the wind forcing becomes a dominant contribution to the estuarine secondary cells. During this season the wind work is maximum (Fig. 7, bottom panel) but the buoyancy flux is still negative (Fig. 7, top panel) thus trying to force an  (Fig. 9, top panels), the river influence is more significant due to the weaker wind forcing and the positive buoyancy forcing. In EXP1 the anti-estuarine MOC is weak and restricted to 200-400 m depths, while the secondary estuarine cells are stronger. However even in this case the dominant overturning circulation is still anti-estuarine, in agreement with the Spall's (2012) freshwater and thermal forcing parameters predictions.
Two additional experiments, EXP3 and EXP4, have been performed by enhancing runoff of all rivers of 50 and 100% (Table 1). The bottom panels of Fig. 8 show an abrupt weakening of the CMOC strength in EXP3 and EXP4: the anti-estuarine cell does not disappear but reduces its intensity and extends only in the Northern Ionian sub-basin up to 400 m depth, while the surface estuarine cell is more pronounced and the deep estuarine cell of the Ionian abyss enlarges over the SAd. Overall the multiannual CMOC at middle depths remains an antiestuarine cell in all the experiments. Figure 10 focuses on year 2001 as we found that in 2001 the middle depth antiestuarine cell disappears with both realistic (EXP1) and augmented runoff (EXP3 and EXP4). This is due to the weak wind work (Fig. 7, bottom panel) and the lowest heat budget (not shown). By comparing Figs. 9 and 10 we then conclude that a large river runoff can produce a positive buoyancy flux without switching off the antiestuarine CMOC cell (as shown in Fig. 9 for year 2002) but a low heat flux and wind work with normal river runoff can in fact reverse it (as shown in Fig. 10 for year  2001).
Overall by comparing the experiments with and without increased runoff, we demonstrate that rivers play a relevant role in the CMOC strength but they do not represent the dominant controlling mechanism. Figure 11 provides the time series of the annual MTS for all the experiments, by considering the averaged value over 100-400 m depths in order to focus on the location of the anti-estuarine CMOC cell. These time series corroborate the previous result on river role in the CMOC strength.  As far as we know an analysis of the MOC for the Eastern Mediterranean (EMOC) has been performed only by Pisacane et al. 2006 andAmitai et al. (2016). Pisacane (Pisacane et al. 2006) EMOC cell resembles our CMOC cell intensity and extension with a well defined antiestuarine character but this cell is computed over a different model domain and the interannual average is shown only for the summer time. Amitai et al. (2016) evaluated the MOC in the Adriatic basin: their maximum value of 0.5 Sv is comparable with our results, which show a maximum value of about 0.3 Sv at 300 m depth close to the Otranto Strait (Fig. 8, left panel). However Amitai et al. (2016) differs because the antiestuarine cell is restricted only to the southern part of the SAd Pit while an additional estuarine cell occupies its northern part.

How do rivers influence the formation of dense water in the Adriatic Sea?
The aim here is to establish how rivers impact the dense water formation processes in the Southern Adriatic Sea and thus impact the CMOC. The conditioning factors of opensea convection and dense water formation in the Southern Adriatic are known to be the permanent cyclonic gyre sustained by wind stress curl, the large buoyancy losses, mainly known to be due to large heat losses, and the inflow of LIW through the Otranto Strait. Moreover the SAd deep layers are partially filled by the NAd dense water which partially the MAd Pit and partially flows southward over the Italian western shelf and reaches the Bari Canyon, usually 2 months after its generation. Here the boundary current of NAd dense water mixes turbulently with adjacent warmer and less saline waters and sinks to the Southern Adriatic Pit (Vilibić and Orlić 2001;Bensi et al. 2013b). Figure 12 shows the dense water volume formed in the NAd, MAd and SAd sub-regions computed as the water volume with larger potential density anomaly than the threshold value, 29.2 kg m −3 .
In the NAd, where the major discharge is concentrated, buoyancy gains due to river runoff oppose the preconditioning factors of DW formation i.e. the winter surface cooling, the outbreaks of cold and dry wind like Bora and the autumn/winter NAd cyclonic gyre (Fig. 12, top panel). Rivers may affect DW volume also in the MAd as they reduce the DW advected from the NAd (Fig. 12, middle panel).
Looking at the full dynamics experiment (EXP1) we find that the mean DWF rate of 0.3 Sv is consistent with the literature estimates (Artegiani et al. 1989(Artegiani et al. , 1997aLascaratos 1993;Cushman-Roisin et al. 2002;Curchitser et al. 2001;Manca et al. 2002;Lascaratos 2004, 2008;Pinardi et al. 2015). The interannual formation rates match the recent estimates by Oddo and Guarnieri (2011) and by Gunduz et al. (2013) There is a perfect agreement with DW volumes found by Gunduz et al. (2013) and the SAd DW formation fits the observations by Cardin et al. (2011). Our results also confirm the exceptional dense water formation in winter 2012 as consequence of an extreme wintertime cooling event, which produces DW volumes capable to saturate the volume of NAd, i.e. 2357 km 3 , and MAd, i.e. 4039 km 3 , subregions (Mihanovic ́ et al. 2013;Janekovic ́ et al. 2014). On the other hand our estimated SAd DW volumes are pretty higher with respect to Oddo and Guarnieri (2011). In our full dynamics experiment, EXP1, we found the transport of NAd dense water toward the middle Adriatic had a peak value of 1.6 Sv in 2012 while 0.6 Sv is found by Janekovic ́ et al. (2014), actually defined over a smaller region thus cannot be directly compared. Our SAd DW formation annual rate reaches the peak of 1.52 Sv in 2006, higher then 0.64 Sv estimate by Oddo and Guarnieri (2011) but our estimate is consistent with the known potential peaks of SAd DW rates, following Mantziafou and Lascaratos (2008) and Vilibic and Supic (2005).
We found EXP2, without river forcing, shows 20-30% larger SAd dense water volumes than EXP1. Previous studies suggest that rivers affect SAd dense water volumes because they reduce the lateral advection of NAd dense waters which are known to flow along the western shelf and slide down into the Southern Adriatic Pit near the Bari canyon Orlić 2001-2002;Lascaratos 2004, 2008;Wang et al. 2007). However, we found this mechanism is not sufficient to account for EXP1 and EXP2 SAd dense water volume difference because the latter is larger than the sum of NAd and MAd dense water volume differences between EXP1 and EXP2. This means that rivers affect SAd dense water volumes not only by reducing the lateral advection of NAd dense water but they also oppose the local preconditioning factors of the dense water formation in the SAd, i.e. the "open sea convection" mechanism, as they locally change the vertical stratification characteristics.
To explain the river impact on the open sea convection, Fig. 13 shows seasonal θ-S diagrams in three Adriatic subregion zonal sections for winter 2009. The main difference between EXP1 and EXP2 θ-S diagrams is the degree of stratification of the water column: river runoff determines a well stratified water column in all the three subregions and different water types with characteristic θ and S values can be distinguished. On the other hand the no-river case, EXP2, shows an almost constant salinity value in the three zonal sections.
Thus we conclude that river runoff primarily influences the dense water formation processes by changing the SAd vertical stratification characteristics. However these stratification changes can determine changes in the CMOC characteristics only if the overall buoyancy forcing is positive, as shown for the 2002 case. In the 2009 case, wind stress work is large and buoyancy forcing still negative, thus balancing the large river runoff stratification effects and producing no relevant changes in the CMOC.

How do rivers influence the Otranto water
exchange and the volume of Adriatic dense water that spreads into the Northern Ionian Sea? Figure 14 shows the Otranto inflow-outflow regime for EXP1 and EXP2. EXP2 shows an "horizontally detached" exchange flow with inflowing water on the eastern side and outflowing water on the western one. EXP1 shows a weaker exchange pattern with lower meridional velocity through all depths and strong exchange flux narrowed to the upper levels.
We know that the Modified Levantine Intermediate Water, MLIW, arriving from the Rhodes area of formation occasionally enters the Adriatic Sea mainly in summer and autumn between 150 and 400 m (Malanotte-Rizzoli et al. 1997;Robinson et al. 2001;Manca et al. 2002;Lascaratos 2004, 2008). Moreover it is well known from the literature that the SAd dense water overflows the Otranto Strait mainly along its western part (Curchitser et al. 2001).
Our modeling findings prove for the first time that rivers reduce both the MLIW inflow on the Eastern shelf and the SAd DW outflow on the western side.
The MLIW inflow is not systematic but depends on the cyclonic/anticyclonic character of the Ionian intermediate water circulation: our experiment covers most of the cyclonic decadal phase of the Ionian near-surface circulation, i.e. 1997circulation, i.e. -2007circulation, i.e. (Gacic ́ et al. 2010Bensi et al. 2013a, b), which favors the entrance of MLIW into the Adriatic basin.
Timeseries of the meridional heat and volume transports at the Otranto Strait are provided in Fig. 15. Rivers are shown to reduce the net ingoing heat transport (Fig. 15, upper panel) and to increase the outgoing volume transport (Fig. 15, bottom panel). The annual values of heat and volume transports we found are fully consistent with the annual averaged heat gain, i.e. 2.92 TW, and volume outflow, i.e. −0.003 Sv, computed by Mantziafou and Lascaratos (2004).
Bottom panel of Fig. 15 confirms river role on the volume exchange at the Otranto Strait as already stressed in Fig. 14. Rivers increase the freshwater outflow at the surface, they reduce the inflow of salt and warm water on the eastern flank at the surface, i.e. the Ionian Surface Water (ISW), and middle depths, i.e. the MLIW; finally rivers reduce the outflow of Adriatic DW through the bottom (Cushman-Roisin et al. 2002), which represents only 30% of the total outflow. Thus the overall effect is an increase of the total outflow at the Otranto Strait. Years 2005 and 2006 are the only ones with no differences in terms of volume transport between the two EXPs, because other forcings (e.g. surface heat losses, ISW and MLIW inflow) are prevailing on river freshwater forcing.
The Adriatic dense waters outflow the Otranto Strait and spread into the abyssal Ionian Sea, generally occupying the layer below the salt waters coming from the Cretan Sea and the Levantine basin (Roussenov et al. 2001;Curchitser et al. 2001;Bensi et al. 2013a). Figure 16 shows the seasonal potential density anomaly of the 200 m layer above the seabed in both EXP1 and EXP2 and their differences with a zoom on the Ionian Sea.   . 40°N). Firm lines mean inflow, dashed lines mean outflow, the while line stands for zero meridional velocity one of the largest dense water volume formation (as shown in Fig. 12).
SAd dense water initially spreads into the Northern Ionian Sea following the topography. The Adriatic outflowing dense water is known to move mainly along the isobaths on the Italian shelf at about 600-1000 m in almost geostrophic balance (Budillon et al. 2010), a secondary branch follows a meridional path driven by the local bathymetry (Hainbucher et al. 2006); both branches slowed down and sink due to friction and turbulent mixing with the ambient waters, thereby creating a nearly homogeneous layer below 1200 m (Curchitser et al. 2001). Wu and Haines (1996) found newly formed Adriatic dense waters fill the Ionian abyss in 2-3 years. Figure 16 shows that the bottom boundary current forced by dense water outflow from the Otranto Strait is characterized in EXP1 by less dense waters that intrude offshore with respect to EXP2 as expected. Furthermore the Adriatic dense water outflow in EXP2 is confined to a narrower band against the Italian shelf with less lateral spreading toward the Ionian Seaia. This suggests that river runoff has consequences on the Northern Ionian Sea water mass structure and abyssal mixing processes. It is interesting to connect such outflow differences with the abyssal estuarine cell of Fig. 8: EXP1 shows a stronger deep estuarine cell than EXP2, expecially on seasonal basis. We speculate the slow abyssal circulation of deep water masses in EXP1, as rivers riduce the SAd DW outflow to narrower bands, as well as the complex topography promote the diapycnal upwelling and result in a stronger abyssal estuarine cell.

Conclusions and future developments
In this paper we investigated the influence of river freshwater inflow on the vertical overturning cell of the Central Mediterranean Sea (CMOC). A twin experiment, with and without runoff, was carried out from the beginning of 1999 to the end of 2012. For EXP1 with realistic runoff estimate, the comparison with satellite SST and ARGO profiles shows that the model is capable to reproduce the major characteristics of the observed water masses. We then take this experiment to be the best estimate of reality and we subtract the river runoff from the simulation and we study the differences between the two experiments. Our study is "mechanistic", i.e. we use a full forcing and dynamics simulation with a realistic estimate of river runoff and then we subtract the runoff to understand its impacts.
The anti-estuarine CMOC has its downwelling branch in the SAd deep water formation areas and thus we analyze first the marginal sea overturning circulation characteristics to understand differences due to runoff. We first compute the surface forcing budget and the theoretical marginal sea estuarine/antiestuarine parameters devised by Spall (2012). The long term and surface average of the buoyancy forcing over the Adriatic Sea is negative but during years of large runoff it changes sign, thus supporting the conjecture that the Adriatic Sea could change from anti-estuarine to estuarine vertical circulation. However the computation of Spall (2012) thermal and freshwater nondimensional parameters demonstrates that river runoff cannot reverse the dominant anti-estuarine character of Adriatic circulation or shut down the deep convection in the basin interior because the system is mainly in the thermal mode and wind work is always positive.
Our results for the realistic runoff case show that the multiannual CMOC for the decade 2000-2012 is a permanent anti-estuarine meridional overturning cell, occupying the Southern Adriatic and the Ionian Sea, plus secondary estuarine cells in the NAd and in the MAd, SAd deep layers as well as in the Northern Ionian abyss. The deep estuarine cell in the Ionian abyss is pointed out for the first time in this study and we speculate this is the result of diapycnal upwelling in the turbulent boundary layer close to the sloping bottom of the Ionian abyss.
A key result is that the CMOC is largely driven by wind work and heat fluxes but large and anomalous river runoff can affect its strength, enhancing the amplitude of the secondary estuarine cells and reducing the intensity of the dominant anti-estuarine cell.
Two additional experiments have been performed by enhancing the realistic runoff estimate of all rivers of 50% and 100%. Both experiments show an abrupt reduction of the antiestuarine CMOC cell on multiannual basis: the antiestuarine cell does not disappear but strongly reduces its intensity and extends only in the Northern Ionian subbasin, while the deep estuarine cell of the Ionian abyss enlarges over the SAd.
We found that on 2001 the middle depth antiestuarine cell disappears with both realistic and augmented runoff, due to the weak wind work and the lowest heat budget. This means that rivers play a relevant role in the CMOC strength but they do not represent its dominant forcing mechanism and the potential role of river runoff on the intensity and direction of the CMOC has to be considered jointly with wind work and heat flux, as they largely contribute to the energy budget of the basin.
We focused on the downwelling branch of the Central Mediterranean MOC, which develops in the Adriatic basin Fig. 16 Seasonal potential density anomaly on a 200 m layer above seabed and difference between EXP1 and EXP2 with zoom on the Ionian Sea due to dense water formation processes. Rivers are demonstrated to affect the Adriatic dense water volumes. Previous studies showed that rivers reduce the dense water formation in the Northern sub-region where most discharge is located. Here we show that rivers also directly affect the vertical mixing processes in the Southern Adriatic sub-region by changing the water column stratification in the SAd and thus decreasing the dense water volumes.
Finally we showed that the Adriatic dense waters overflowing the Otranto Strait are less dense in a realistic runoff regime and confined to a narrower band against the Italian shelf with less lateral spreading toward the Ionian Sea center.
Future investigations will point to an advanced implementation of river discharge in the model and the extension to the whole Mediterranean Sea. The next questions could consider the connection between the CMOC and the zonal overtuning cell of the Mediterranean Sea that supports the flow of LIW in the different deep water formation areas of the eastern and western Mediterranean Sea.
variables, but the local values may significantly reduce depending on season and latitude and moving towards the shelf areas. In the Northern Adriatic Sea it reduces up to about 3-5 km in summer and 1 km in winter (Paschini et al. 1993;Masina and Pinardi 1994;Bergamasco et al. 1996). This means our model can explicitly resolve the mesoscale activities in the Adriatic Sea, at least in the Southern sub-basin, on seasonal as well as on interannual basis with the only exception of the Northern Adriatic where the model may result eddy-permitting but not eddy-resolving.
The model bathymetry, covering both the Adriatic and Ionia Sea, is taken from the US Navy 1/60° bathymetric database DBDB1 using bilinear interpolation.
A total of 121 unevenly spaced z-levels with partial steps were adopted in the vertical direction. Partial steps allow a better representation of the bathymetry. The higher resolution in the top layers (23 levels in the top 35 m which is the mean depth of the NAd subregion) leads to an improved simulation of the bottom flow in the NAd and vertical mixing during higher stratification in the summer.
There are two open boundaries on the eastern and western sides of the model domain. Open boundaries data are provided as monthly means and involve the following prognostic variables interpolated on the model grid: zonal velocity (u 3d ), meridional velocity (v 3d ), potential temperature ( ), salinity (S), and the sea surface height ( ).
For both the lateral open boundary conditions, LOBCs, and the initial conditions, ICs, data are taken from daily analysis of the operational Mediterranean forecasting System, MFS (Tonani et al. 2008;Oddo et al. 2009;Pinardi and Coppini 2010) based on the same code, NEMO, and covering the whole Mediterranean basin with 1/16 horizontal resolution.
The numerical schemes adopted for the LOBCs are described below.
Marchesiello's algorithm (Marchesiello et al. 2001) was used for active tracers. It consists of the 2D radiation condition plus a relaxation/nudging term as follows: where φ is the tracer ( or S), nested is the coarser model (MFS) solution for the tracer interpolated on our model grid and provided monthly. The time scale for the nudging term, τ, is constant and equal to 1 day for inward propagation and 15 days for outward propagation.
For outward propagation, i.e. C x > 0 where C x is the component of the phase velocity normal to the boundary, the tangential component is set equal to zero, C y = 0.
For inward propagation, C x < 0, the algorithm prescribes C x = C y = 0 thus reduced to a relaxation condition.
For the horizontal velocity components, u 3d and v 3d , the imposition scheme is used and thus, the incoming and outgoing information is totally determined by the coarser model data, irrespective of the inner solution.
In addition, the horizontal velocity component normal to each boundary is uniformly adjusted according to the "interpolation constraint" procedure (Pinardi et al. 2003) in order to preserve the total volume transport after data interpolation from the coarse to the fine resolution grid.
For the barotropic velocities, u BT and v BT , Flather's scheme (1976) was adopted. The barotropic velocity component normal to the eastern and western boundaries is given by Flather's equation: where nested is the coarser model sea surface height at the boundary interpolated over the finer model grid, is the finer model sea surface height at the boundary and u BT nested is the coarser model normal barotropic velocity over the finer model grid computed as The tangential barotropic velocity is set equal to zero: v BT = 0.
In Flather's formula, values at the boundary follow a "zero gradient boundary condition" which means B = B−1 (subscript B stands for boundary line values). This avoids numerical instabilities.
The bottom boundary condition is applied only on momentum and consists of a quadratic friction.
No slip boundary conditions are adopted along the coastline for tangential velocity.
In order to define the air-sea interaction, the vertical fluxes of momentum, heat and salt and the vertical velocity were parameterized at the sea surface. These parameterizations are the surface boundary conditions (SBCs) of the model. Wind stress and heat flux components are computed by means of "bulk formulae" (Castellari et al. 1998;Maggiore et al. 1998;Oddo and Guarnieri 2011;Madec 2008) using atmospheric data provided by the European Centre for Medium Weather Forecasts (ECMWF). These atmospheric data (2 m air temperature, 2 m dew point temperature, total cloud cover, mean sea level atmospheric pressure, meridional and zonal 10 m wind components) are operational analyses with a 6 h frequency and with 0.5° horizontal resolution up to 2008 and 0.25° thereafter. It should be noted that the increase in spatial resolution of the ECMWF fields could lead to a spurious increase of the wind stress magnitude.
Only the precipitation rate (P) data are extracted from the CMAP (CPC, Climate Prediction Center, Merged Analysis of Precipitation) monthly data set with a horizontal resolution of 2.5° × 2.5°. This coarse precipitation data set is a weakness of our model formulation but precipitation data set from ECMWF was still too poor for the first decade of the XXI century (Haiden et al. 2015) and this would have offset the buoyancy budget in an unphysical manner.
The surface boundary condition for temperature involves a balance between solar short-wave radiation Q s (computed using Reed's formula, 1977), long-wave radiation Q l (computed using Bignami et al. 1995), latent Q e and sensible Q h heat fluxes [by means of bulk formulae proposed by Kondo (1975)].
Reed's formula is: where Q tot is the clear-sky radiation, C is the fractional cloud cover, is the noon sun altitude in degrees, and is the sea surface albedo. The albedo is computed as a function of the sun zenith angle for each grid point from Payne (1972). The Bignami formula is: where is the ocean emissivity, is the Stefan Boltzmann constant, e A is the atmospheric vapor pressure, T S is the sea surface temperature predicted by model, T A is the 2 m-air temperature.
The sensible Q h and latent Q e heat fluxes are parameterized through the Kondo bulk formula: where A = A p, T A , r is moist air density, r is the relative umidity, C p is the specific heat capacity at constant pressure, C H and C E are the turbulent exchange coefficients computed according to Kondo (1975), L e is the latent heat of vaporization, e sat is the vapor pressure, |v| is the wind speed modulus, and p A is the atmospheric pressure.
For the heat flux boundary condition at the surface, we assume: where T r is the Jerlov (1976) transmission coefficient for a "clear" water type and K t is the vertical mixing coefficient for traces.
The wind stress involved in the surface boundary condition for momentum is calculated from the relative winds with the formula: where U rel = u w − u s = (u rel , v rel ) is the relative wind field that is the 10 m wind horizontal velocity u w subtracted from the sea surface horizontal velocity u S , 0a is the density of the moist air and C D T a , T s , u w is the drag coefficient which depends on air temperature, sea surface temperature and wind amplitude according to Hellerman and Rosenstein (1983).
The momentum boundary condition at the surface is: | v rel are the wind stress components and K m is the vertical mixing coefficient for momentum.
The freshwater balance defined as evaporation minus precipitation and runoff (with the latter divided by the cell area of the river mouth), E-P-R/A, is directly involved in the surface boundary conditions for salinity and for vertical velocity. The evaporation rate (E) is calculated by the latent heat flux according to E = Q e L e . The salinity boundary condition at the surface reads: where is the sea surface elevation and S z= is the ocean model surface salinity except prescribed ad-hoc salt values at river mouths.
The surface boundary condition for the vertical velocity is as follows: where w is the vertical velocity.
With regard to the dynamics, the following choices were selected: vector invariant form for momentum advection, bi-Laplacian operator for lateral diffusion and horizontal eddy viscosity coefficient equal −5 × 10 7 m 4 s −1 according to a tuning procedure starting with MFS values, implicit vertical diffusion and TKE turbulence closure scheme (Mellor and Blumberg 2004) to provide the vertical eddy coefficients.
With regard to the active tracers: MUSCL advection scheme, bi-Laplacian operator for lateral diffusion and horizontal eddy diffusivity coefficient equal to −3 × 10 7 m 4 s −1 according to a tuning procedure, implicit and TKE dependent vertical diffusion.
Appendix 2: The computation of Spall parameters for semi-enclosed seas Following Spall's studies (2004, 2010, 2011, 2012 on overturning circulation in marginal seas, we computed two non-dimensional coefficients which represent the relative balance between heat and freshwater budget in the interior of the Adriatic Sea due to air-sea interaction and the lateral eddy fluxes that advect warm and salty water in the basin interior and are detach from the cyclonic boundary current which inflows along the eastern side from the open ocean and encircles the marginal sea. The combinations ∕ and ∕ are called respectively thermal and freshwater forcing parameter and are described below: where A is the area of the Adriatic sea surface (from model domain), is the relaxation constant for the basin sea surface temperature toward the atmospheric temperature (from Spall 2011), f 0 is the Coriolis parameter, T is the thermal expansion coefficient (from Cessi et al. 2014), S is the haline expansion coefficient (from Cessi et al. 2014), H is the depth of the sill (from model domain), P is the perimeter of the basin interior (from model domain), c P is the thermal capacity (from Cessi et al. 2014), L is the L is the width of the sloping topography over which the inflowing boundary current lies (thus computed from model results as the cross-shore width of the inflowing boundary current along the eastern shelf of the Adriatic basin).
The variable = h x −̄x∕̄z = −0.33 represents the topography slope over the mean isopycnal slope in the boundary current, thus both computed along the Southern Adriatic eastern shelf (x-direction stands for zonal direction and z means depth). To note that δ has been computed by considering a zonal transect of potential density anomaly at 40.8ºN (so just north of the Otranto Strait) on annual basis and focusing on the eastern side of the basin. For cyclonic boundary current, δ < 0 and the topography acts to stabilize the boundary current and reduce the amount of lateral eddy flux into the interior. The quantity c = 0.025e 2 = 0.05 is an efficiency coefficient that depends on the bottom slope and regulates the eddy heat flux from the boundary current into the interior (Spall 2004).
The non-dimensional parameter = cP∕L is the ratio of the heat flux toward the basin interior due to lateral eddies compared to that advected into the Adriatic Sea through the inflowing boundary current along the Southern Adriatic eastern shelf. The inflowing boundary current is assumed to be a geostrophic current in thermal wind balance. The value of is very small for stable boundary currents and increases for boundary currents that are sufficiently unstable that they lose all their heat to the interior of the basin before it is carried all the way around the marginal sea.
Moreover the thermal and freshwater forcing parameters required to compute T * , that is the difference between the inflowing temperature and the temperature of the atmosphere over the interior of the marginal sea, as follows : where T 1 is the mean temperature of the inflowing current along the eastern boundary derived from the EXPs, T A is the mean 2 m temperature over the Adriatic basin extracted from ECMWF 25 km dataset.
Finally the surface freshwater flux is defined as follows: All the quantities described above enable to compute the thermal and freshwater forcing parameters in Eqs. (22) and (23), giving: As discussed by Spall (2011), ∕ is a measure of the relative influence of lateral eddy heat fluxes from the boundary current into the basin interior compared to heat loss to the atmosphere. For ∕ ≪ 1, lateral eddy heat flux from the boundary is very strong and leads to a relatively warm basin interior so that T ≈ T 1 , if ∕ > 1 the boundary current is relatively stable and the atmosphere is able to strongly cool the basin interior so that T ≈ T A . Similarly ∕ describes the relative role of surface forcing and lateral eddy fluxes in the salinity balance. Large values of ∕ indicate dominance of atmospheric forcing implying freshwater gains in the basin interior that are not balanced by lateral eddy fluxes of salt from the boudary current, and small values indicate strong lateral eddy fluxes.
In order to evaluate the shutdown of deep convection and the reversal of the overturning circulation, the temperature and salinity anomalies of the basin convective water mass have been computed as normalized differences of T (i.e. ΔT) and S (i.e. ΔS) between basin interior and boundary currents (Spall 2012). A set of two non-dimensional equations has been derived to compute ΔT and ΔS (Spall 2012), these equations include the non-dimensional parameters ∕ and ∕ and describe ΔTand ΔS as function of basin geometry, atmospheric forcing and lateral eddy fluxes.
The simplified formula suggested in Spall (2012) are: where T and S for the basin interior and T 1 andS 1 for the inflowing current have been computed over 50-200 m depth. The ratio ΔS∕ΔT < 1 means the stable circulation state is in the "Thermal mode" with surface heat losses and freshwater gains prevailing the lateral eddy advection of warm and salt water in the basin interior (thermal and freshwater forcing coefficients are significantly high). In this case the density contrast is dominated by the temperature difference and the water in the interior of the marginal sea is more dense than that in the boundary current.
If ΔS∕ΔT > 1, the stable circulation state is in the "Haline mode" with surface heat and freshwater budget of the interios basin favoring the lateral eddy advection of warm and salt water (thermal and freshwater forcing coefficients are low enough and the latter is eventually negative). In this case the density contrast is dominated by the salinity difference and the water in the interior of the marginal sea is less dense than that in the boundary current. Thus the boundary corrent detaches from the eastern shelf and spreads in the interior basin. This means in haline mode the surface boundary current is in the opposite sense, anticyclonic around the coastline, and the deep convection in the basin interior is not longer supported with reversal of the meridional overturning circulation.
According to our findings, both EXPs are in the thermal mode and EXP1 is closer to the threshold limit for shutdown of deep convection than EXP2.
Results collected for EXP1 and EXP2 are summarizes in Table 3 and show that deep convection in the Southern Adriatic, surface cyclonic boundary current and anti-estuarine overturning circulation of the Adriatic basin characterize both experiments but in EXP1, with a realistic parameterization of river runoff, the freshwater forcing coefficient is negative and the ratio ΔS∕ΔT close to 0.5. This corroborates strong river discharge in the Adriatic Sea has the potential to trigger the shutdown of deep convection and the weakening of the anti-estuarine overturning circulation. The model-diagnosed quantities are the thermal forcing parameter ∕ , the freshwater forcing parameter ∕ , the temperature anomaly of the convective water mass ΔT and the salinity anomaly of the convective water mass ΔS Parameter EXP1 EXP2 ∕ 5.0 × 10 −5 4.9 × 10 −5 ∕ −2 × 10 −2 +7 × 10 −4 ΔT 0.35 0.28 ΔS 0.10 0.03 ΔS∕ΔT 0.28 0.12