The Mediterranean Sea Overturning Circulation

The time-mean zonal and meridional overturning circulations of the entire Mediterranean Sea are studied in both the Eulerian and residual frameworks. The overturning is characterized by cells in the vertical and either zonalor meridional planes with clockwise circulations in the upper water column and counterclockwise circulations in the deep and abyssal regions. The zonal overturning is composed of an upper clockwise cell in the top 600m of the water column related to the classical W ü st cell and two additional deep clockwise cells, one corresponding to the outﬂow of the dense Aegean water during the Eastern Mediterranean Transient (EMT) and the other associated with dense water formation in the Rhodes Gyre. The variability of the zonal overturning before, during, and after the EMT is discussed. The meridional basinwide overturning is com- posed of clockwise, multicentered cells connected with the four northern deep ocean formation areas, located in the Eastern and Western Mediterranean basins. The connection between the W ü st cell and the meridional overturning is visualizedthrough the horizontal velocities vertically integrated across two layers above 600m. The component of the horizontal velocity associated with the overturning is isolated by computing the divergent components of the vertically integrated velocities forced by the inﬂow/outﬂow at the Strait of Gibraltar.


Introduction
The overturning circulation of the global ocean plays a key role in setting the stratification of different basins because it regulates the ocean carbon budgets and provides the mechanism for the supply of oxygen and other tracers from the surface to the deep ocean. The conceptual framework of the global overturning circulation has recently advanced, indicating a middepth clockwise overturning cell that connects the Northern Hemisphere deep water formation areas to the wind-driven Southern Ocean upwelling (Toggweiler and Samuels 1993;Gnanadesikan 1999;Lumpkin and Speer 2007;Marshall and Speer 2012;Talley 2013) and an abyssal counterclockwise overturning cell driven by bottom-enhanced diapycnal mixing that balances Antarctic bottom water formation. In other words, the global meridional Denotes content that is immediately available upon publication as open access. overturning circulation has been depicted as composed of counterrotating meridional cells, and a recent overview is offered by Cessi (2019).
This new picture has emphasized the importance of adiabatic, along-isopycnal motion in the clockwise intermediate-depth overturning, with Ekman transport and eddy-flux of buoyancy in the Antarctic Circumpolar Current playing a fundamental role (Marshall and Radko 2003;Wolfe and Cessi 2010). In contrast, diapycnal mixing is important in the dynamics of the counterclockwise abyssal cell (Nikurashin and Vallis 2011).
Much of the understanding of the global overturning circulation dynamics has been advanced by using the concept of residual circulation (Andrews et al. 1987;Marshall and Radko 2003). Taking the meridional overturning as an example, the velocity is averaged in density layers of thickness h and can be decomposed into its mean and eddy components, that is, where the overbar indicates the zonal and time average, and the prime indicates departures from the zonal and time average. The left-hand side of Eq. (1) is the thickness-weighted or residual velocity in the chosen layer; the first term on the right-hand side (rhs) is the traditional, Eulerian-average velocity; and the second term is the velocity due to large-scale circulation gyres, standing waves, and transient eddies. The advantage of the residual framework is that it includes the contribution of eddy fluxes due to velocity structures with zero time or zonal averages, but nonzero fluxes [the second term on the rhs of Eq. (1)]. Thus, the residual velocity is representative of the transport of scalars, such as temperature, salinity, and density (Plumb and Mahlman 1987;Bachman and Fox-Kemper 2013). From the thicknessweighted velocity in Eq. (1) a transport streamfunction is calculated in density and latitude coordinates. In this paper the residual overturning circulation of the Mediterranean Sea is studied and compared to the Eulerian overturning circulation. The geometry and bathymetry of the Mediterranean Sea is described in Fig. 1. Several ;500-m-deep sills subdivide the basin into areas where different processes dominate. The narrow (;7-km width) and shallow (300-m depth) Gibraltar Strait sill connects the Mediterranean to the Atlantic. The exchange at the Strait is characterized by a baroclinic two-layer flow related to the energy and buoyancy budget of the semienclosed basin (Cessi et al. 2014). This two-layer flow provides low-salt waters to the Mediterranean, balancing the increasing salt time tendency due to evaporation at the air-sea interface, thus helping to maintain a stratification in the basin. The wide and shallow (;500-m depth) Sicily Strait sill divides the Western Mediterranean (WMED) from the Eastern Mediterranean (EMED), allowing only surface and intermediate waters to be exchanged. The EMED has two marginal seas of its own, the Adriatic and the Aegean Sea, where deep waters are formed. Deeper sills (;800-m depth) connect the Adriatic and Aegean Seas to the EMED through the Otranto and four Cretan Sea Straits, among which the Kasos and Antikythira Straits are 1000 and 700 m deep, allowing dense waters to exit the Cretan Sea.
The Mediterranean Sea has been called a miniature ocean for climate studies (Bethoux et al. 1999;Tsimplis et al. 2006) because it is a basin with deep and intermediate mass formation processes generating a vigorous vertical circulation. The WMED overturning circulation is connected to the deep water formation areas of the Gulf of Lion Houpert et al. 2016) and an overall understanding of its structure is still lacking. The semienclosed geometry of the Mediterranean Sea enables a zonal overturning cell connecting the two-layer flow at the Strait of Gibraltar to remote areas of the Eastern Mediterranean basin . The zonal circulation of the Mediterranean Sea was described for the first time by Wüst (1961) analyzing a vertical section of salinity, and by Zavatarelli and Mellor (1995) in terms of the Eulerian zonal transport streamfunction obtained with a coarse-resolution numerical model. The traditional picture of the EMED overturning considered the Adriatic Sea as the only source of deep water formation (Zavatarelli and Mellor 1995), but after 1995 it has become clear that two EMED meridional cells are possible, one with the downwelling branch in the Adriatic Sea and the other in the Aegean Sea. Between 1988 and 1998, a climatic event, termed the Eastern Mediterranean Transient (EMT; Roether et al. 1996;Klein et al. 1999Klein et al. , 2000Manca et al. 2003;Gertman et al. 2006) was observed, characterized by the formation of very dense waters in the Aegean Sea outflowing through the Cretan Straits and filling the abyssal plain of the EMED. The overturning circulation linked to the Adriatic Sea has been the subject of recent studies (Verri et al. 2018;Somot et al. 2006), and the overturning circulation generated from both the Adriatic and Aegean Sea sources has been discussed by Amitai et al. (2017). There is no description of the overturning circulation generated by the intermediate and deep water formation in the Rhodes Gyre area (see Fig. 1), although this formation process has been documented (Pinardi et al. 2015;Hecht and Gertman 2001). In this work, we show for the first time the role of the Rhodes Gyre water mass formation area in the EMED overturning circulation.
The first multiscale descriptions of the meridional overturning circulation that included the Aegean deep water source were given by Robinson et al. (2001), Tsimplis et al. (2006), and Pinardi and Masetti (2000). These early studies outline a Mediterranean Sea overturning composed of three main spatial scales: basin scale, subbasin-scale gyres, mesoscale eddies, and standing waves. However, this amount of work produced a qualitative scientific synthesis and not a quantitative estimate of the transport streamfunction. Although we are aware that other scales might also be important, such as the submesoscales (Pinardi et al. 2016;Pascual et al. 2017) and the baroclinic tides generating mixing (Morozov et al. 2002), there is no dataset currently available to investigate these high-frequency components.
The zonal and meridional overturning circulations in the Mediterranean Sea have received attention in the context of simulations of paleoceanographic scenarios, where atmospheric forcing and sill depths were varied in simplified and coarse-resolution general circulation models (Myers 2002;Alhammoud et al. 2010). These studies showed the importance of the Gibraltar sill depth and the value of the net water balance (today with evaporation largely overcoming precipitation) in determining the overturning circulation. Non-eddy-resolving simulations (Pisacane et al. 2006;Amitai et al. 2017) were also used to study the stability of the EMED overturning cells as a function of the strength of the two different deep water formation sources, that is, the Adriatic and Aegean Seas, showing that multiple equilibria are possible depending on the dominance of one deep water source over the other.
A comprehensive description of the zonal and meridional overturning circulation of the Mediterranean Sea for the present climate is still missing. The goal of this work is to describe the observational evidence of the time-mean meridional and zonal overturning system in the Mediterranean Sea with an eddy-resolving model that assimilates data. We use part of an eddy-resolving reanalysis dataset (Simoncelli et al. 2017), spanning the period from 1987 to 2013, to characterize the basin-scale overturning circulation and compute the residual transport streamfunction to understand the role played by eddies and permanent gyres. The key questions are: 1) What is the structure of the zonal and meridional overturning circulations in the Mediterranean Sea? 2) What is the contribution of eddying/gyre motion to that structure? 3) What are the connections between the zonal and the meridional overturning circulation? 4) What is the variability of the Mediterranean overturning before and after the EMT?
The paper starts with section 2 briefly describing the reanalysis dataset. Section 3 presents the analysis of the zonal overturning circulation compared with Wüst (1961) schematic, documenting the vertical density structure of the basin. Section 4 analyses the meridional overturning circulation in the WMED and EMED. Section 5 explores the connection between the zonal and meridional overturning circulation. Section 6 presents the EMT influence of the Mediterranean zonal overturning circulation. Section 7 offers a discussion and conclusions.

The reanalysis dataset
A reanalysis is a three-dimensional retrospective estimate of the ocean dynamical variables obtained through a data assimilative numerical experiment (Masina and Storto 2017). In this work, we used the Copernicus Marine Environment Monitoring Service reanalysis of the Mediterranean area (Simoncelli et al. 2017), which provided daily mean fields for over 27 years, from 1987 to 2013.
The reanalysis modeling component uses the Nucleus for European Modelling of the Ocean (NEMO), version 3.2, ocean general circulation model (Madec 2008) that solves the primitive equations on the sphere using the Jackett and McDougall (1995) equation of state. The horizontal grid has a resolution of 1/168 3 1/168, and the vertical grid has 72 unevenly spaced levels from the surface to a maximum depth of 5000 m. Since the model uses vertical partial cells, the thickness of the bottom layer is allowed to vary as a function of the geographical location to yield a better representation of the bathymetry. The model domain covers the entire Mediterranean Sea and a portion of the Atlantic Ocean in order to resolve salinity and heat exchanges through the Gibraltar Strait (Oddo et al. 2009). Specifically, the model domain spans from 18.1258W to 36.258E in longitude and from 30.18758 to 45.93758N in latitude. The Mediterranean part of the domain is illustrated in Fig. 1. The model is forced by interactive momentum, heat, and freshwater fluxes forced by the ERA-Interim atmospheric fields (Dee et al. 2011) and climatological precipitation from Xie and Arkin (1997). The data assimilation system uses a three-dimensional variational method developed by Dobricic and Pinardi (2008). The assimilated data consist of in situ observations from conductivity-temperature-depth (CTD), expendable bathythermograph (XBT), mechanical bathythermograph (MBT), bottle, and Argo data as well as remotely sensed data from satellite sea surface temperature and along-track altimetry (TOPEX/Poseidon, ERS-1 and -2, Envisat, Jason-1 and -2). The in situ data were contained in the MEDAR/MEDATLAS dataset (Maillard et al. 2005) and other datasets collected by the authors along the years. These data are now part of the SeaDataNet European archive (https://www.seadatanet.org/). Gridded maps of sea surface temperatures derived from satellite data are not directly assimilated into the model, but were used to correct iteratively the heat flux at the air-sea interface. The background error multivariate correlation matrix was estimated from a historical model simulation and varied seasonally in 13 regions of the Mediterranean Sea each having different characteristics (Dobricic et al. 2005).
The reanalysis was initialized in January 1985 using a January temperature and salinity climatology obtained from the historical dataset of in situ data collected from year 1900 to the beginning of the reanalysis. The years 1985-86 are considered as a spinup period for the system and are not used in the following analysis.
3. The Mediterranean overturning circulation from salinity and density mapping a. Vertical salinity mapping In a landmark study, Wüst (1961)  The 1987-2013 salinity mean field is computed at approximately the same section as Wüst's and the comparison is provided in Figs. 2 and 3 (bottom panels). The similarities between the old description and the new estimate are striking, except the new estimate gives a more homogeneous deep water column and a weaker LIW tongue in winter. Salinity values are consistent between the two analyses even if the LIW and deep salinity values are 0.1 psu higher in the new estimate than in Wüst's dataset both in the WMED and the EMED, probably due to the positive salinity climate trend documented by many authors of 0.001 psu yr 21 (e.g., Vargas-Yáñez et al. 2017).

b. Vertical density mapping
One approach to determine the pathways of the overturning circulation is to examine maps of potential density and the associated velocity field (Talley 2013). We introduce our analysis by considering a vertical transect in the Mediterranean zonal direction and vertical meridional sections across the four deep water formation areas. In their study of a previous Mediterranean Sea reanalysis dataset, Pinardi et al. (2015) confirmed the existence of four well-defined water mass formation areas: the Gulf of Lion Gyre in the WMED , the southern Adriatic, the Cretan Sea, and the Rhodes Gyre in the EMED (Velaoras et al. 2014). Since the densest waters in the Mediterranean are found as shallow as 500 m below the surface, the potential density referenced to the surface, s 0 , is a useful variable to use in this analysis.
A zonal transect of s 0 similar to Wüst's longitudinal section of Figs. 2 and 3 is shown in Fig. 4. West of longitude 68W, fresh Atlantic Waters (s 0 , 27.3 kg m 23 ) enter the Mediterranean Sea, while east of 208E intermediate isopycnals rise to the surface. The densest waters (s 0 . 29.2 kg m 23 ) are stored in the abyssal plains of the EMED. Below the depth of 500 m, isopycnals of s 0 ; 29.15 and s 0 ; 29.225 kg m 23 are flat, with a larger stratification in the deep EMED water column. The circulation that can be deduced from this figure corresponds to that obtained with the salinity sections of Figs. 2 and 3, except that the subsurface return flow is now along the 28.1-29.1 kg m 23 isopycnals, crossing the basin east-west at depths between 150 and 400 m.
Meridional transects connecting the deep water formation areas and the southern boundary of the basin are shown in Fig The overall density analysis of meridional sections shows the emergence of two distinct meridional areas: 1) above 200 m, the northern regions are dominated by the outcropping of isopycnals and the southern regions are dominated by AW flow, and 2) below 300 m, the deep isopycnals downwell from the open ocean formation areas toward the coasts and upwell on the southern shores, with the exception of the Cretan Passage section. From these density sections it is not possible to draw a complete picture of the deep circulation of the basin. Therefore the smoothed time-mean vertical velocity was calculated at 1000 m, that is, an interface deep enough to be affected only by the vertical motion of the deep water masses (a qualitatively similar picture is found at 500 m).
The model vertical velocity field is displayed in Fig. 6. Along the four sections shown in Fig. 5, the mean vertical velocity shows downwelling along the northern shores and upwelling in the southern coasts, strongly enhanced near the boundary. This is not surprising, since enhanced mixing and friction is found near rough topography, associated with vertical motion. Thus, large downwelling motions are likely to be found near areas of intense mixing, such as deep water formation regions. It is interesting to note that there are subregional differences, especially in the Alboran Sea, Ionian Sea, and the Cretan Passage.
Recently, Waldman et al. (2018) has also looked at the areas of downwelling that were most likely associated with the overturning circulation of the Mediterranean Sea. They report measurements of downwelling motion along the northern escarpment of the Gulf of Lion, which is in agreement with our model results. However, they cannot estimate all areas with observations. Our work extends their finding that overturning circulation downwelling is localized in the northern boundary current areas, that is, the Liguro-Provençal current for the Gulf of Lion and the Asia Minor current in the northern Levantine basin, [for the nomenclature, see Pinardi et al. (2015)]. Deep upwelling is instead found along most of the southern Mediterranean coastlines, which might be connected to abyssal overturning or to the return flow of the upper water column overturning.
To better understand the basin-scale overturning, we will now compute the streamfunctions in the vertical plane and relate them back to this downwelling-upwelling system along the boundaries of the domain.

The Mediterranean Sea overturning system
To understand the mean vertical circulation of the Mediterranean Sea, we compute the time-mean zonal and meridional Eulerian and residual streamfunctions associated with the velocity field estimate provided by the reanalysis.
The Eulerian meridional (zonal) streamfunction is calculated by integrating the meridional (zonal) velocity y (u) first in the vertical direction, then in the zonal (meridional) direction, and finally averaging over the reanalysis period. The resulting streamfunctions are where y B 1 and y B 2 are the meridional boundaries, x B 1 and x B 2 are the zonal boundaries, T 5 t 1 2 t 0 is the temporal averaging interval, and H is the bottom bathymetry. With these definitions, positive (negative) values represent clockwise (counterclockwise) circulations in the vertical plane.
To include the transport associated with the contributions from eddies, gyres, and standing waves, the residual streamfunctions c * are defined from the velocities integrated in density (s 0 ) space, rather than depth, then integrated in one horizontal direction and then averaged in time, as follows: c * mer (y,s) 5 where H is the Heaviside function. This is the streamfunction corresponding to the thickness weighted averaged velocity defined in Eq.
(1) and, formally, in Young (2012). With these definitions, c * is the transport occurring below the isopycnals, which in turn is a function of all three spatial dimensions plus time. Within the residual framework, the natural vertical coordinate iss, but we remapped c * onto a depth-like coordinate z in order to ease the comparison with the Eulerian representation, whose natural vertical coordinate is depth. The definition of z is This mapping implies a distribution ofs which is given bys where x 5 x or x 5 y for the zonal and meridional streamfunctions, respectively (see the appendix). The quantitys is contoured together with c * to show the location of potential density surfaces used in the mapped coordinates.
a. Zonal overturning circulation Figure 7 shows the Eulerian mean c zon and residual c * zon averaged over the entire period of the reanalysis. The Eulerian zonal overturning (top panel of Fig. 7) is characterized by three structures: the first is a shallow clockwise cell, above the 29.1 kg m 23 potential density surface, corresponding to AW flowing eastward and LIW flowing westward. We called this shallow overturning circulation the Wüst cell. Below 300-700 m, two other cells dominate, one in the WMED and the other in the EMED. In the WMED, a counterclockwise cell occupies the region below 700 m, while in the EMED a multiple centers clockwise cell extends to the bottom. The WMED deep zonal counterclockwise cell fills the region east of 58W, occupying the deep regions of the Alboran Sea, the Algerian basin, and the Tyrrhenian Sea, the latter located at about 108-158E.
The EMED deep clockwise cell maximum is centered at the longitude of the Aegean deep water outflow, that is, 228E, and it is connected with the surface intensified Wüst cell. Several other subsurface clockwise cells exist in the Levantine, disconnected from the surface. The clockwise cell centered at approximately 288E is notably associated with the dense water formation processes occurring in the Rhodes Gyre area.
The residual streamfunction (bottom panel of Fig. 7) shows the different dynamical balances dominating the deep WMED and EMED overturning cells due to the transport by permanent gyres, standing waves and transient eddies [the second term in Eq. (1)]. The WMED deep counterclockwise residual overturning is stronger than its Eulerian counterpart indicating that the transport by stationary gyres and eddies reinforces the mean abyssal vertical circulation. The residual velocities (tangent to the streamfunction by definition) cross the isopycnals indicating that diapycnal mixing is important for the maintenance of this cell. In the Tyrrhenian Sea, the residual overturning circulation mostly flows along isopycnals indicating a more adiabatic balance. But it is in the EMED where the major differences between Eulerian and residual zonal overturning appear. The residual circulation is weaker at depth than the Eulerian counterpart, and a detectable residual counterclockwise circulation emerges. The large and deep Eulerian clockwise cell is broken down in two parts in the residual representation: the deep overturning at 228E is along isopycnals, while the circulation around the secondary maxima located approximately in the Rhodes Gyre area of 288E is across isopycnals. The deep to abyssal zonal overturning circulation of the EMED is strongly influenced by the transport due to the eddy/permanent gyres component. In our analysis, the Aegean deep water source is apparent in the zonal overturning cell of the Mediterranean Sea, because our dataset straddles the period of the EMT event. This feature will be further discussed in section 6. Values of the zonal overturning circulation counterclockwise and clockwise cells are on the order of 1 Sv (1 Sv [ 10 6 m 3 s 21 ) in both the WMED and the EMED. The surface Wüst cell is stronger in the Eulerian mean than in the residual form, supporting the conclusion that the eddy field in the first 300 m weakens the zonal transport of the Wüst cell.
In synthesis, the major outcomes of the zonal overturning circulation analysis are 1) the revisited Wüst cell occupies the first 500 m of the water column and is weakened by the gyre/eddy components of the transport field; 2) a clockwise deep cell is present in the region of the Aegean outflow, dominated by along isopycnal motion, and a weaker, less deep clockwise cell is present in the Rhodes Gyre area; and 3) the permanent gyre/eddy transport field strengthens the deep overturning counterclockwise cells of the zonal overturning.

b. Meridional overturning circulation
We now discuss the meridional streamfunction, comparing the Eulerian mean c mer and residual c * mer previously defined in Eqs. (3) and (5). We first show the WMED and EMED and then the whole basin-scale meridional overturning system. The streamfunction for the WMED considers a zonal averaging up to a section in the Sicily Strait and the one for the EMED an averaging carried out from the same section up to the Levantine basin coastlines.

CIRCULATION
The Eulerian WMED overturning circulation is shown in Fig. 8 The residual WMED meridional overturning changes this picture, showing only two well-separated counterrotating cells, both stronger than their Eulerian counterparts. For the clockwise cells, the Eulerian streamfunction indicates a volume transport approximately of 0.36 Sv, while the residual counterpart is 0.88 Sv. For the counterclockwise cells, the minimum of 20.22 Sv in the Eulerian framework is increased to 20.7 Sv in the residual streamfunction. The boundary between the two cells is at approximately 398N, that is, the latitude marking the division between the permanent cyclonic gyre of the Gulf of Lion and the eddy-dominated anticyclonic area of the Algerian current described by Pinardi et al. (2015). The clockwise northern cell is associated with the deep water formation areas in the northern basin, normally centered around 418-428N. The downwelling branch of the clockwise cell occurs near along the northern boundary of the domain, where vertical velocities are negative, as illustrated in Fig. 6. The western Mediterranean deep water formation and spreading phenomena is an eddy dominated process, documented in several papers (Madec et al. 1996;Demirov and Pinardi 2007;Send and Testor 2017). Thus, the residual overturning is stronger than the Eulerian mean. The residual southern counterclockwise cell extends from approximately 300 m to the bottom and it is stronger than the northern clockwise cell. The downwelling branch of this abyssal cell is located along the northern African escarpment, below the upwelling branch illustrated in Fig. 6. The residual overturning, for both counterclockwise and clockwise cells, crosses isopycnals, indicating that diapycnal mixing dominates the WMED meridional overturning circulation.

CIRCULATION
We now consider a spatial domain bounded in the zonal direction by the Sicily Strait and the Adriatic, Balcanic, Turkish, and Levantine basin coasts. To have an overall expression of the EMED basin-scale vertical circulation, both deep water formation marginal sea areas are included, that is, the Adriatic and the Aegean Sea. The Eulerian and residual EMED overturning circulations are shown in Fig. 9.
With reference to the Eulerian framework (top panel of Fig. 9), the upper 250 m of the water column are occupied by a counterclockwise surface cell, which we interpret to be associated with the Ekman transport and its shallow return flow. This surface counterclockwise cell is interrupted at 36.58N due to the physical gap of the Sicily Strait at this latitude. Below 250 m a clockwise multicentered overturning cell is present. In general we can distinguish between two clockwise centers, one extending southward of 36.58N and the other north of it. North of 36.58N, the overturning streamfunction is connected to the Adriatic and Aegean deep water formation areas with the maximum centered at about 378N, corresponding to the Aegean deep water formation area (Cretan Sea). The northward maximum extension of this clockwise cell includes the Adriatic Sea deep water formation regions with transports about 0.1 Sv, in agreement with the recent estimates of Verri et al. (2018). It is clear that this northern clockwise center is weaker than the Aegean one, probably due to the specific time period chosen. In fact, from 1989 to 1998 the EMT deep water formation event is known to coincide with smaller southern Adriatic deep water formation rates (Amitai et al. 2017).
South of 36.58N, several clockwise overturning centers exist, with maxima between 250-and 1000-m depth, some of them extending to almost 3000 m. These cells are associated with the Aegean Sea deep water formation area, its outflow during the EMT and the Levantine Deep Water formation processes in the Rhodes Gyre. Velaoras et al. (2014) report that dense water outflow from the Aegean continued between 2000 and 2010 thus contributing to give a large overturning in this area. These clockwise cells add up to approximately 0.5-0.6 Sv at the latitudes of 348-368N.
A deep counterclockwise Eulerian cell occupies the depth below 1000 m in the meridional region of the Ionian Sea (north of 358N). It is bounded by the s 0 ; 29.2 kg m 23 isopycnal surface, below which the stratification is weak. This deep cell has been described before (Zavatarelli and Mellor 1995;Pisacane et al. 2006;Verri et al. 2018) and it is stronger than the Eulerian WMED deep counterclockwise cell, probably due to the deeper bathymetry of the eastern basin.
In the residual framework (bottom panel of Fig. 9) the vertical circulation appears quite different, as it is the case for the WMED. The clockwise part of the streamfunction is subdivided into three major cells: the weakest is connected to the Adriatic Sea water formation areas, the second to the Aegean Sea, and the third, south of 35.58N, is connected to the Aegean dense water outflow and the Rhodes Gyre dense water formation area. The southern residual clockwise cells extend over most of the water column and the deep counterclockwise cell has practically disappeared. While the clockwise overturning cell north of 36.58N is essentially the same in the Eulerian and residual framework, the southern residual circulation is dipolar and stronger than in the Eulerian description, adding up to ;0.7 Sv in each center. The second clockwise center, located between 35.58 and 36.58N, corresponds to the areas where the Aegean Sea dense waters formation occurs, and to the Rhodes Gyre formation area around 36.58N. The southernmost positive maximum, located between 338 and 358N, corresponds geographically to the so-called Cretan Passage, where the dense Aegean Seawater outflows and to the southern Rhodes Gyre. All the clockwise residual circulations are cross-isopycnals indicating that diapycnal mixing processes are important over the whole basin. It is noteworthy that the EMED southern clockwise cell has never been depicted before.

3) THE MERIDIONAL BASIN-SCALE OVERTURNING SYSTEM
In Fig. 10 we analyze the Eulerian and residual meridional overturning circulation for the whole Mediterranean Sea. The Eulerian basin-scale streamfunction (top panel of Fig. 10) shows the superposition of the WMED and EMED clockwise cells occupying the depths between ;100-250 and 1000 m with cell boundaries reaching 2000 m. The maximum value of 1 Sv is reached in the regions south of 388N. The WMED clockwise cells (north of 388N) are weaker than their EMED counterparts, while the WMED abyssal counterclockwise circulation is the largest.
The residual basin-scale meridional overturning brings a balance between the southern and northern clockwise cells. The WMED clockwise maximum is 0.9 Sv, larger than in the Eulerian framework and comparable to the southern one, 0.76 Sv, corresponding to the Cretan Passage and Levantine basin overturning. The basin-scale residual counterclockwise deep to abyssal circulation has an absolute value of the transport around 0.5 Sv, about half of the clockwise counterparts.

Connection between the zonal and the meridional overturning
There is evidence that the LIW, occupying the layer between 150 and 500 m, impacts the deep water formation processes (Wu and Haines 1996) in the WMED. Published results showed that LIW increases the salinity of the waters exposed to the winter cooling influencing the deep water mass formation in the Aegean, Adriatic, and Gulf of Lion areas (Klein et al. 2000;Theocharis et al. 2002;Schroeder et al. 2016). The Western Mediterranean Transition (WMT) has been recently documented that connects the changes of the WMED deep water formation rates to the saltier LIW arriving from the EMED, after the EMT (Zunino et al. 2012).
In addition, the AW, occupying the upper-layer (0-150 m), also influences the deep water properties and formation rates. Malanotte-Rizzoli et al. (1997) found that, during the EMT the diversion of AW toward the Adriatic Sea was an important factor determining the deep water formation in the Rhodes Gyre and the subsequent formation of high-density waters in the Aegean Sea. The AW diversion was due to a specific gyre in the Northern Ionian Sea, so-called Northern Ionian Gyre (NIG) that lasted a decade, from 1987 to 1996, documented for the first time by Pinardi and Navarra (1993) and recently discussed in great details by Reale et al. (2016). After the 1987-96 period the NIG reversed and AW reached again the Levantine while the EMT ended. Since the flow of AW and LIW composes the Wüst cell, there could be important interactions between the zonal and meridional overturning dynamics of the basin. In this section we start to explore the connection between the Wüst cell and the meridional cells by looking at the horizontal flow field in the upper and lower branches of the cell and separating the divergent and rotational components of the flow field. The divergent component is connected to the inflow-outflow system at Gibraltar and it contributes to the vertical downwellingupwelling components of the meridional overturning of the basin.
Considering the inflow-outflow system at the Gibraltar Strait (5.58W, figure not shown), we found that the zero crossing of the zonal velocity is located around 150 m and the outflow extends slightly lower than 600 m. The time-mean value of the inflow is 0.92 Sv, while the outflow is 0.88 Sv for a net transport of 0.04 Sv, which is in close agreement with the literature (Baschek et al. 2001). We analyze the horizontal circulation in two vertical layers, the first extending from the surface to z 1 5 150 m and the second from z 1 to z 2 5 600 m. Specifically, we define the transports and the vertically averaged velocities: where (u, y) is the horizontal velocity field; Dx and Dy are the model grid cells in the longitudinal and latitudinal directions, respectively; k 5 1, 2 are the vertical layer labels (increasing downward) with z 0 equal to the sea surface, taken approximately at z 0 5 0; and Dz k is the layer thickness. The transport field and the vertically averaged velocities are shown in Fig. 11. The upper-layer flow is characterized by intensified current segments that, starting from the Gibraltar Strait, can be traced up to the Levantine basin (top panel of Fig. 11). We can recognize some of the wellknown structures of the surface circulation of the Mediterranean Sea, as described in details by Pinardi et al. (2015). The westward lower-layer flow of the Wüst cell is composed of jet segments that bring LIW and other intermediate waters toward the Strait of Gibraltar (bottom panel of Fig. 11). To note is the intensified branch of the lower-layer flow along the Libyan coasts.
To understand the effects of this two-layer flow field on the overturning circulation we need to disentangle the divergent and rotational component of the flow field.
The divergent component is associated with the vertical velocity component of the overturning cells (Fig. 6).
We subdivide the two-layer transport vector in the two components: where (U k,f , V k,f ) is the divergent component of the transport, while (U k,x , V k,x ) is the divergenceless or rotational component. The divergent and rotational components of the transport are related to the velocity potential f k and streamfunction x k by and respectively. The symbols d i and d j indicate the (dimensionless) difference between two adjacent points in the longitudinal and latitudinal direction respectively. The equation for the velocity potential is defined as while for the rotational component, The solution of Eq. (14) is subject to nonhomogeneous boundary conditions, given by the normal velocity at the Gibraltar Strait, that is, at all the points corresponding to the Gibraltar entrance, while the normal velocity to the solid boundary is zero everywhere else: all of the Gibraltar inflow-outflow is attributed to the divergent component of the velocity. The solutions of Eq. (15) are subject to homogenous boundary conditions everywhere, including the Gibraltar Strait. This implies that all the x k streamlines are closed.
The potential function f k and the velocity vecto, (1/Dz k )[(U k,f /Dx), (V k,f /Dy)] are shown in the top panels of Fig. 12. The potential field is large-scale, except at Gibraltar and the Sicily Strait where it intensifies. Because the velocity associated with the velocity potential singles out the component with horizontal divergence, it describes the overturning circulation, both zonal and meridional. However, the total velocity is the sum of the potential and rotational components, and the latter steers the flow into large-scale gyres subdividing the basin into several opposite circulation components (bottom panels of Fig. 12). In the upper layer (bottomleft panel of Fig. 12), the flow field is characterized by cyclonic large-scale gyres in the WMED and the Levantine Sea, and one anticyclonic in the Ionian Sea. The rotational velocity component of the NIG is weak and broken into opposite sign eddies or subbasin-scale gyres probably due to the time averaging in the period of NIG reversal -96 NIG was anticyclonic, the 1997-2013NIG is cyclonic, as shown in Pinardi et al. (2015]. In the lower-layer flow field (bottom-right panel of Fig. 12), the gyres are at smaller scales and multicentered. It is clear that the flow shown in Fig. 11 is equally composed by both the potential and rotational components of the velocity field.

The EMT influence on the Mediterranean zonal overturning
Now that the general overturning structure is defined over the whole 1987-2013 reanalysis period, we can analyze some aspects of the interannual variability of the zonal residual overturning circulation shown in the bottom panel of Fig. 7. One of the deep clockwise cells in the EMED is centered at the longitude of the Aegean deep water outflow, that is, 228E, and it is connected to the EMT event. Roether et al. (2014) defined the characteristics of the EMT: the Aegean Sea, which up to 1987 had only formed intermediate-depth water masses, began to discharge unusually dense waters in the EMED. The large Aegean water outflow started in 1992 and lasted until 1995 (Roether et al. 2014), so that we consider the interval between 1991 and 1996 to be the EMT peak period, with 1987-90 the onset period and 1997-2013 the decay phase. The Aegean outflow during the peak phase was 10 times that in the two neighboring periods.
The average residual zonal streamfunction for these three periods is shown in Fig. 13. In the onset phase the deep clockwise cell at the Aegean Sea longitude of 228E is weak (;0.9 Sv) and relatively shallow, only down to 750 m. During the EMT peak phase, the clockwise cell moved deeper, down to 2750 m and more, and its amplitude almost tripled, up to ;2.4 Sv. During the decay phase the anticlockwise cell detached from the Wüst cell, stabilized between 500 and 1750 m and weakened back to ;1 Sv.
We then conclude that the EMT is the source of the deep anticlockwise Wüst cell in the EMED and that the interannual variability of deep water mass formation events in the Aegean is an important forcing for the clockwise cells of the zonal overturning circulation.

Conclusions
This paper analyzes the Mediterranean Sea largescale overturning circulation in its zonal and meridional components. The basin-scale vertical circulation was studied using the Eulerian and residual circulation frameworks, with the latter including the contribution from permanent gyres, eddies and standing waves. This analysis was made possible by the availability of an eddy-resolving reanalysis of the whole basin (Simoncelli et al. 2017).
Starting with the zonal overturning circulation, we documented the Wüst cell, well reproduced by the reanalysis in the salinity and density sections, thus giving a general validation of the reanalysis data. The zonal Eulerian streamfunction shows a new, deep clockwise cell in the region of the Aegean deep water outflow, where the Wüst cell goes down to 3000-m depth. In the residual framework this cell flows largely along isopycnals indicating an almost adiabatic buoyancy balance of deep water outflow from the Aegean. In the Levantine Sea, two other clockwise cells are present, one at the location of the EMT dense water overflow and the other in the Rhodes Gyre area. These clockwise zonal overturning cells are new and we argue that one of them is connected to the EMT dominating the deep water formation processes in the Levantine Sea for almost a decade, from 1989 to 1998. The meridional overturning circulation of the whole Mediterranean has been analyzed for the first time. In summary, the meridional overturning circulation can be characterized as follows: 1) the WMED meridional overturning circulation is dominated by two clockwise cells, that is, a clockwise cell connected to the deep formation area in the Gulf of Lion and a counterclockwise cell located in the abyssal regions of the southern basin, and 2) the EMED is dominated by three clockwise cells associated with the three deep and intermediate water formation areas and a weak counterclockwise cell.
In the Eulerian framework, the clockwise overturning is multicentered occupying the depths between 250 and 2500 m. The EMED clockwise vertical circulation is dominated by the Aegean, Cretan Passage, and Rhodes Gyre overturning, with a weak Adriatic Sea contribution. The WMED clockwise overturning reaches a maximum of 0.36 Sv while the EMED meridional clockwise overturning reaches 0.5-0.6 Sv. Conversely, the counterclockwise abyssal overturning is stronger in the WMED than EMED.
In the residual framework, the clockwise meridional cells are characterized by large cross-isopycnal flows, leading to the conclusion that diabatic mixing is an important component of the basin buoyancy balance. This differs from the middepth Atlantic meridional overturning circulation balance, which shows a largely adiabatic residual flow (Marshall and Speer 2012). The limited meridional extent of the Mediterranean Sea might be the cause of the major role of diapycnal mixing.
The Wüst cell has been connected for the first time to the meridional overturning by studying the divergent and rotational components of the two-layer flow field composing the clockwise vertical flow field. The divergent part of the transport that produces the vertical velocities composing the downwelling and upwelling branches of the meridional circulation is shown to be at large scales and intensified in the Alboran Sea and Sicily Strait. Such a divergent flow field is steered by the rotational transport in the basin, giving rise to jet segments and local divergences. In accordance with the results of Waldman et al. (2018) downwelling areas are found along the northern shelf of the basin where intense boundary currents develop, that is, the Asia Minor current and the Liguro-Provençal current.
Last, we discuss the EMT influence on the residual zonal clockwise cell: it is found that the strength of the deep clockwise zonal cell in the area of the Aegean Sea outflow changes in phase with the different phases of the EMT, reaching 2.4-Sv maxima during the peak phase of the climatic event.
In conclusion, our analysis of the Mediterranean Sea global overturning circulation has revealed several new features, not discussed before to the best of our knowledge. Specifically, we highlight the multiple meridional clockwise cells of the EMED, the abyssal counterclockwise vertical circulation of the WMED and the structure of the Wüst cell that is composed of the classical shallow cell in the upper ;600 m of the water column and two deep clockwise zonal cells, one out of the Aegean and the other corresponding to the Rhodes Gyre. This study gives rise to many new questions that should be answered in the future, for example, what is the role played by the deep vertical counterclockwise cell in the WMED? What is the seasonal and interannual variability of the divergent part of the flow field? How is such large buoyancy mixing balance produced?
Finally, future investigations of the overturning circulation of the Mediterranean Sea should analyze how this active vertical circulation is connected to the basinscale air-seawater, heat and momentum exchanges and how did this change in the past. This would finally elucidate the potential for oxygen depletion in the Mediterranean Sea as well as in other semienclosed seas of the World Ocean.
water, h ck 5 0 a cell totally inside the topography, and 0 , h ck # 1 a partial cell near the bottom boundary.
The value of the potential density carried by the model represents the average of s over the kth grid cell: where dz f k 5 (Dz f k 2 Dz f k11 )h ck . 0: Fields are remapped from height into isopycnal coordinates using the piecewise parabolic method (PPM; Colella and Woodward 1984;Carpenter et al. 1990).

a. Building the interpolator
To specify uniquely a parabola describing the variable distributionŝ k (z) over the kth volume cell z f k11 # z # z f k , three quantities are necessary: