Drivers of the decadal variability of the North Ionian Gyre upper layer circulation during 1910–2010: a regional modelling study

A long simulation over the period 1901–2010 with an eddy-permitting ocean circulation model is used to study the variability of the upper layer circulation in the North Ionian Gyre (NIG) in the Eastern Mediterranean Sea (EMed). The model is driven by the atmospheric forcing from the twentieth century reanalysis data set ERA-20C, ensuring a consistent performance of the model over the entire simulation period. The main modes of variability known in the EMed, in particular the decadal reversals of the NIG upper layer circulation observed since the late 1980s are well reproduced. We find that the simulated NIG upper layer circulation prior to the observational period is characterized by long-lasting cyclonic phases with weak variability during years 1910–1940 and 1960–1985, while in the in-between period (1940–1960) quasi-decadal NIG circulation reversals occur with similar characteristics to those observed in the recent decades. Our simulation indicates that the NIG upper layer circulation is rather prone to the cyclonic mode with occasional kicks to the anticyclonic mode. The coherent variability of the NIG upper layer circulation mode and of the Adriatic Deep Water (AdDW) outflow implies that atmospheric forcing triggering strong AdDW formation is required to kick the NIG into an anticyclonic circulation 1–2 years later. A sensitivity experiment mimicking a cold winter event over the Adriatic Sea supports this hypothesis. Our simulation shows that it is the multi-decadal variability of the salinity in the Adriatic Sea that leads to periods where low salinity prevents strong AdDW formation events. This explains the absence of quasi-decadal NIG reversals during 1910–1940 and 1960–1985.


Introduction
Since the 1980s variability has been observed in the upperlayer circulation in the North Ionian Gyre (NIG), which quasi-periodically changed from anticyclonic to cyclonic and vice versa on a decadal time scale (Bergamasco and Malanotte-Rizzoli 2010;Gačić et al. 2010). The first reversal of the basin-wide circulation, from cyclonic to anticyclonic, was discovered in 1986-1987 by in situ CTD data (Malanotte-Rizzoli et al. 1997). Since then, the occurrences of two anticyclonic phases have been identified by altimetry data during 1987and 2006(Bessières et al. 2013Gačić et al. 2010). In addition, two cyclonic phases have been documented: the first during 1998-2005 (Borzelli et al. 2009;Larnicol et al. 2002;Pujol and Larnicol, 2005) and the second starting in 2011 (Gačić et al. 2014). This second cyclonic phase was terminated by an unexpected 'premature' inversion in 2012 due to the extremely harsh winter (Gačić et al. 2014). This inversion lasted only until the beginning of 2013 after which a cyclonic circulation was restored. The 'premature' inversion illustrates potential irregularity in the timing of phase reversals.
Understanding the variability of the NIG circulation is of importance, since reversals between the two modes, affecting the spreading of the Atlantic Water (AW) eastward and salty Levantine waters westward, have a direct influence on the physical and biogeochemical properties of the adjacent regions Gačić et al. 2011Gačić et al. , 2013Mihanović et al. 2015) and potential impacts on the broad marine ecosystem functioning (Lavigne et al. 2018;Peharda 1 3 et al. 2018). Moreover, similar shifts of circulation regimes have also been observed in other regions of the world ocean (e.g. Montoya et al. 2011). Thus, related studies in the Mediterranean Sea, where many processes are similar to those occurring in the global ocean but at much smaller spatial and temporal scales (Bethoux et al. 1999), will give valuable insights which are transferable from the miniature ocean to a larger-scale context (Rubino et al. 2020).
The driving mechanisms behind the quasi-regular reversals has sparked an intense discussion in the Mediterranean scientific community. Such reversal was once attributed to changes in the wind stress curl over the northern Ionian Sea (Demirov and Pinardi 2002;Korres et al. 2000;Pinardi et al. 1997). However, later findings suggested that the local wind stress is not the primary forcing of the surface circulation by arguing that the geostrophic vorticity rate of change is of the opposite sign of the vorticity source contributed by the wind stress (Borzelli et al. 2009;Gačić et al. 2010). Gačić et al. (2010) explained the observed cyclical reversal since the late 1980s in terms of a feedback mechanism called the Adriatic-Ionian Bimodal Oscillating System (BiOS) that links the Ionian and the Adriatic Seas according to the following mechanisms: (1) During the anticyclonic phase of the NIG circulation, a large portion of the AW is advected into the northern Ionian interior and further deflected into the southern Adriatic Sea, leading to a decrease of the salinity. As a result, the newly formed Adriatic Deep Water (AdDW) has a lower density while it spreads into the Ionian Sea, leading to a deepening of the isopycnal surfaces along the AdDW pathway. The consequent stretching of the upperlayer column leads to a reversal of the circulation from anticyclonic to cyclonic conditions. (2) During the cyclonic phase, the high-salinity Levantine waters propagate into the northern Ionian interior and subsequently into the southern Adriatic Sea, making the water saltier and denser. The increased density of the AdDW outflow induces a rising of the isopycnal surfaces at the flanks of the Ionian Sea and thus a squeezing of the upper-layer column. The vorticity is consequently decreased, leading to a more anticyclonic circulation pattern. The BiOS feedback mechanism was further extended to a coupled Adriatic-Ionian-Aegean system, emphasizing the role of the internal thermohaline pumping mechanism in regulating the salinity distribution and the dense water formation processes in the Adriatic and Aegean basins Reale et al. 2016Reale et al. , 2017Theocharis et al. 2014;Velaoras et al. 2014).
The above mechanisms suggest that the reversal of the NIG circulation is directly linked to the intensity of the AdDW formation events. We have to take into consideration that deep water formation does not only depend on the preconditioning salinity field but also on the buoyancy loss induced mainly by winter air-sea heat flux (Beuvier et al. 2010;Lascaratos 2004, 2008). Therefore, the winter meteorological conditions can modify the efficiency of the AdDW formation and in turn perturb the intrinsic variability of the NIG circulation that is sustained by internal processes as proposed by the theories mentioned above. The impact of such perturbation is evidenced by the NIG reversal in 2012, when a large amount of very dense AdDW was produced consecutive to strong heat loss (Gačić et al. 2014). Another possible source of perturbation arises from the freshwater flux over the Adriatic Sea, which is a key factor influencing the haline properties of this region (Skliris et al. 2007). Considering the potential impact of these perturbations, one might speculate that the NIG circulation in the last century might show a different temporal variability from that documented in the recent past. However, the variability of the NIG circulation prior to the comprehensive observations could not be ascertained, mainly due to the sparseness of observations in the Eastern Mediterranean Sea (EMed) (Rixen et al. 2005). Consequently, the role of the external forcing (in this paper referred to atmospheric forcing) in modulating the NIG variability in that period has not yet been assessed. To fill this gap of knowledge, it has been suggested that long-term model simulations with realistic external forcing would be a good approach (e.g. Roether et al. 2007). A few modelling efforts have tackled the NIG variability back to 1960 Theocharis et al. 2014). However, comparable simulations covering earlier periods are not available.
Accordingly, in this paper, we aim at: (1) identifying the variability of the NIG circulation over the period from 1910 to 2010, with a focus on the period before the 1980s, during which the observations are quite sparse, and (2) assessing the roles of the atmospheric forcing and variations of large-scale atmospheric circulation (e.g. the North Atlantic Oscillation, NAO) in modulating the variability of the NIG circulation. We perform and analyze simulations with an eddy-permitting regional ocean circulation model, driven with consistent forcing over the period from 1901 to 2010. To evaluate the model performance, we assess the simulated variability against that identified from observations. Then, we analyze the simulated variability for earlier periods. Based on an additional sensitivity experiment we test the hypothesis that NIG reversals can be triggered by strong AdDW formation events and we deduce the major dynamical features involved in this phenomenon.

Model description
The model used in this study is the primitive equation ocean general circulation model MPIOM (Max-Planck-Institute ocean model, Marsland et al. 2003). We set up a regional version for the Mediterranean with a model domain covering the entire Mediterranean Sea, the Black Sea and a small Atlantic box (Fig. 1). Horizontally, the model uses a curvilinear grid with a resolution of approximately 9 km, allowing some eddy variability to occur in the simulation. Vertically, 30 unevenly spaced layers are set up with the thickness increasing towards the bottom, with nine layers in the top 100 m and partial cells at the bottom. At the western boundary of the Atlantic Ocean box, a sponge zone of ~ 70 km is set, over which the modelled temperature and salinity are relaxed towards climatological monthly fields from the World Ocean Atlas Boyer et al. 1998) with a relaxation time scale of 31 days. Considering that the Mediterranean Sea is an evaporative basin, we set up a restoring of the sea level towards zero at the sponge zone in the Atlantic box following Mikolajewicz (2011), to prevent continuous sea level sinking in the interior of the model domain. This ensures that the climatological net transport through the Strait of Gibraltar equals the climatological net evaporation in the Mediterranean and Black Seas. The model has a free surface and there is no relaxation on the surface or in the interior of the model domain.
Heat fluxes at the surface are calculated using standard bulk-formulas from prescribed atmospheric forcing and modelled sea surface temperature (SST). The net shortwave radiation is calculated from a constant albedo while the upward longwave radiation is calculated from the modelled SST. The penetration of shortwave radiation into the ocean is calculated according to Jerlov optical water types Ib in the whole model domain (Jerlov 1976). The prescribed atmospheric forcing are 3-hourly 2-m air temperature, dew point temperature, downward shortwave and longwave radiative fluxes, precipitation, wind stress and speed from the twentieth century reanalysis dataset ERA-20C (ECMWF's first atmospheric reanalysis of the twentieth century: 1901-2010, Poli et al. 2016). The forcing data are interpolated onto the model grid, using ocean points only. The choice of ERA-20C data ensures a consistent atmospheric forcing and thus a consistent model performance over the entire simulation period. This is important for this study as we will use the simulation for the years 1985 onwards to validate the model with observations and will draw conclusions from the model results for earlier periods without extensive observation coverage.
Freshwater fluxes are calculated from prescribed precipitation, river runoff and evaporation, with the latter determined from the interactively calculated latent heat flux. 286 rivers surrounding the model domain are taken into consideration. Monthly runoff of the Po River from January 1901 to December 2010 is derived from Zanchettin et al. (2008). For other 13 rivers, including the major ones such as the Nile and Rhone, the values are derived from the UNESCO RivDis database (Vörösmarty et al. 1998). Since measurements of those rivers cover only part of the simulated period (information about those rivers can be found in Appendix A, Table S1), climatic monthly values are calculated from the first 10 years of observations and are applied to each year before the observations. The same treatment is applied to the period after the available observations until December 2010 by using the last 10 years' data. In this way, the limited river runoff data is extended to the entire simulation period. For the other rivers not available from RivDIS we use climatic data from GlobalNEWS2 (Mayorga et al. 2010). The fraction of the total runoff based on observations peaks in 1925-1980 and decreases sharply afterwards, which indicates a lack of observations for the recent decades in the database (Appendix Fig. S1). This method allows reproducing at best the varying river run-offs over the desired period but it underestimates their variability at both the beginning and the end of our simulation.
The model is initialized from MEDATLAS II climatology (MEDAR Group 2002) and run for 747 years (9 cycles) with realistic atmospheric forcing from 1901 to 1983, during which the Mediterranean climate does not show apparent anthropogenic climate change like increasing in air temperature. Constant monthly river runoffs corresponding to the beginning of the twentieth century are used in the spin-up phase. By doing this, the model is able to reproduce the interannual variability that is induced by the atmospheric forcing. Compared to starting directly from climatology, this type of spin-up ensures that all the water masses found in the model are internally generated. It also minimizes the artificial drift of the model, but at the cost of some deviations from the observed climatology. After the spin-up the transient simulation runs from 1901 to 2010 with time varying river runoffs. The first years are considered to be affected by artificial drift and are not included in our analysis. Therefore, the following analysis is restricted to the years 1910-2010.

Evaluation of simulated NIG circulation since the late 1980s
To assess the model's capability in generating the variability of the EMed we evaluate the simulated NIG variability over the last 3 decades against observations. A good model performance would enhance our confidence in the results for periods without sufficient observations. The two phases of the NIG circulation are reproduced in our simulation. We choose two typical years, when the cyclonic and anticyclonic modes are well developed, for a detailed analysis of the corresponding circulation patterns. Velocities are averaged over the upper 124 m (upper 10 model layers) to represent the upper layer circulation. Salinity fields are averaged over the same depth in order to distinguish the AW and Levantine water pathways, as low salinity water masses (< 38.5 psu) potentially have their origin in the Atlantic, while high salinity water masses (> 38.5 psu) are potentially originated from the Levantine basin (Fig. 2).
In 1991, the NIG circulation is characterized by an anticyclonic mode (Fig. 2a), which is well documented by many observational studies (e.g. Malanotte-Rizzoli et al. 1997). The AW is deflected northeastwards from the Strait of Sicily and is clearly visible in the northern Ionian as an anticyclonic meander. The AW encounters surface waters from the Levantine Basin, which flow cyclonically northwards along the Ionian eastern flank, and prevents the salty Levantine waters from spreading into the Adriatic Sea. Due to the dilution by the fresh AW, the northern Ionian Sea and the southern Adriatic Sea are occupied by low salinity waters.
The simulated water mass distribution in this year is in good agreement with that observed in a hydrographic survey carried out in 1991 (Malanotte-Rizzoli et al. 1999). In 2000, a well-developed cyclonic circulation pattern is visible (Fig. 2b). The broad cyclonic flow leads to the intrusion of salty waters from the Levantine Basin into the entire northern Ionian Sea and the Adriatic Sea. The northern branch of the AW flow is suppressed or even stopped, while the eastern branch is reinforced. As a consequence, the Mid Ionian Jet is formed and further propagates into the Levantine Basin. A similar circulation pattern occurring over 1997-2006 was described by Pinardi et al. (2015) based on a reanalysis data set for the Mediterranean Sea.
To investigate the temporal variability of the NIG circulation, we define an index based on the sign of the zonal component of the current averaged over the upper 124 m in the northern portion of the Ionian Sea (black box in Fig. 2). The time evolution of this index shows clearly interannual to decadal variability over the period of 1985-2010 (Fig. 3). The simulated variability is in good agreement with indices calculated from altimetry data (e.g. Figure 3 in Gačić et al. 2014) and cruise observations ( Table 2 in  instance, the geostrophic current shows the appearance of a second anticyclone in 2006, but it appears 2 years earlier in our simulation. This is mainly due to the presence of mesoscale eddies, which are quite active during the NIG phase transition period and thus introduce big uncertainties in determining the time of reversals (Gačić et al. 2014).
The simulated annual mean salinity and density of the water masses in the intermediate layer of the southern Adriatic Sea are compared with that derived from the CTD data of 39 cruises carried out between 1985 and 2008 (red and blue lines and diamonds in Fig. 3). The CTD data is averaged in an area of 1° lat × 1° lon centered on 41.75°N, 17.75°E in the Southern Adriatic Pit for the layer between 200 m and 800 m (Cardin et al. 2011;Gačić et al. 2010). The model results are averaged in the same region (black box in Fig. 1) over 206-727 m depth (model layers 14-20). Both the observations and the simulation show that the salinity and density co-vary with the NIG upper layer circulations. The decreasing of the salinity coincides with the anticyclonic circulation and vice versa. The maxima of salinity and density occur at the time when the increase in the NIG circulation index is strongest or 1 year later. The simulation agrees well with the observations before 1990 and after 2002. Deviations of the simulation from the data are found in the period 1990-2002, during which the simulated salinity is lower than the observations, leading to an underestimation of the water mass density. An unexpected increasing of the salinity is found around 1995, which is not documented by the CTD data. The underestimation of salinity is probably due to the overestimation of the inflow of fresh AW water into the Adriatic during the anticyclonic phase in the early 1990s. Around 1994, the anticyclonic circulation is weakened, suggesting more saline Levantine waters flowing into the Adriatic and thus a mild increase of the salinity. The variation of the simulated density resembles that revealed by the CTD data, even though the increase in 1990s is underestimated to some extent.
The model also reproduces the general distribution of climatological SST, SSS and circulation of the Mediterranean Sea well (see Appendix B) and the evolution of the Eastern Mediterranean Transient (EMT) event (Roether et al. 1996) (see Appendix C). Therefore, we conclude that the model provides satisfactory performance in reproducing the climatology and the variability in the EMed over the time of observations and is suitable for investigating earlier periods without extensive observation coverage.
In the period 1940-1960, the variability of the NIG circulation resembles that of the well assessed period from 1985 onwards. In these two periods, cyclonic circulation modes are substituted by pronounced anticyclonic modes on a decadal time scale. The other two periods, over 1910-1940 and 1960-1985, during which the NIG seems to be quite calm, are characterized by phases of weak decadal variability in the circulation mode. The NIG circulation persists in a relatively stable cyclonic mode over more than 20 years. In these phases, there are tendencies towards a reversal of the NIG circulation, but an anticyclonic mode is barely reached.
Obviously, the variability of the NIG circulation of the recent decades (1985 onwards) does not fully represent all characteristics of the long-term variability. The modeled multi-decade periods with weak decadal variability can not be solely explained by internal feedback mechanisms which were proposed in previous studies (e.g. the BiOS). We infer that the external forcing, referring to the atmospheric forcing, plays an important role in shaping the oceanic conditions and thus modulating the variability of the NIG circulation.
A striking feature is the similarity between the variability of the NIG circulation (Fig. 4a) and of the intensity of the AdDW outflow (Fig. 4b) which is quantified by the meridional overturning stream function across an east-west transect south of the Strait of Otranto (at 39.84°N, see Fig. 1). The transect has a model depth of 1165 m, deeper than the Strait of Otranto which has a depth of ~ 700 m. We refer Where v represents the monthly mean meridional component of the horizontal velocity field at location x, y and depth z , − H is the bottom of the sea, E is the lower bound of the integration in the east while W is the upper integration bound in the west.
The stream function reveals both the amount and by the location of its maximum indirectly also the density of the AdDW outflow (Fig. 4b). In some years (e.g. 1987, 2004 etc.), large amounts of AdDW leave the southern Adriatic Sea through the Strait of Otranto being dense enough to continue as deep flows in the Ionian Sea (see maxima in Ψ at 806 m in Fig. 4a). In the following we will refer to such events as strong AdDW formation events. In other years (e.g. 1924, 1962), the outflow is confined to the upper layers and does not reach the near bottom layers. Furthermore, there are years (e.g. the 1910s, 1965-1980) with no substantial outflow, indicating a weak or even absent deep water formation in the Adriatic Sea.
Overall, the overturning stream function shows quite substantial variability, similar to the variability of the NIG circulation (Fig. 4). Apparently, there is a better correlation between the variability of the NIG index and of the stream function at 806 m depth than at 191 m. In 1963, a relatively large amount of water flows out from the Adriatic Sea at the depth of 191 m but fails to trigger a complete NIG reversal, indicating that both the amount and density of the AdDW outflow affect the NIG circulation. Anticyclonic circulations are established after strong AdDW formation events with a lag of 1-2 years, while cyclonic patterns are linked to relatively calm periods, during which the outflow of the AdDW is weak or even disappears. The time lag between the strongest AdDW outflow and the maximum of the anticyclone is discussed in Sect. 4.
We therefore hypothesize that the occurrence of the anticyclonic mode is directly triggered by strong AdDW formation events, which do not happen every year but are determined by the combined effect of extreme winter heat loss due to the external climate forcing and the salinity increase in the Adriatic Sea. The hypothesis is further tested in Sect. 4 with a sensitivity experiment.

Long-lasting cyclonic phase in history
A main feature, revealed by our long-term simulation, is the existence of two 'calm periods' characterized by a longlasting cyclonic phase and the lack of strong AdDW formation events in 1910-1940and 1960-1985). Prior to this study, this feature has not been reported. Based on the BiOS theory (Gačić et al. 2010), the circulation reversals can be sustained by internal processes. This implies that the NIG circulation reversal and the AdDW formation are interacting without depending on the external forcing. The noticed absence of strong AdDW formation is not able to be explained by the BiOS theory and is difficult to validate due to the lack of observations, but our simulation offers the possibility to investigate the mechanisms leading to these 'calm periods'. We focus on the period 1960-1985, during which the formation of AdDW is rather weak.
In general, the occurrence of strong AdDW formation is a result of the combination of excessive winter convection and preconditioning characterized by a salinity increase in the intermediate layer Lascaratos 2004, 2008). The winter (NDJF) air-sea heat exchange over the Adriatic Sea is a major contributor to the buoyance loss (Josey 2003;Mantziafou and Lascaratos 2008). The surface net heat flux anomalies of our study based on ERA20C re-analysis data and the calculation by Josey (2003) with NCEP/NCAR data show a general agreement (Fig. 5). For the 'calm period' of 1960-1985 we find two parts with quite different fluxes: in the first part, over the 1960s, relatively strong heat loss occurs in consecutive winters. In particular, the heat loss in 1963 is significant, implying an extreme harsh winter that can induce sufficient buoyancy loss at the sea surface to trigger a strong convection. The second period, covering the 1970s and the first half of 1980s, is characterized by positive surface heat flux anomalies, indicating an increase of the stability of the water column, which hampers convection. Therefore, we conclude that the absence of strong AdDW formation over 1970-1985 is ascribable to mild winters but the strong winter heat loss over the 1960s disputes the absence of strong AdDW formations.
The other factor influencing AdDW formation is the preconditioning of the water column, which we deduce from the time evolution of the annual mean salinity averaged over 206 m − 727 m in the Southern Adriatic Sea (Fig. 6). After 1955, a remarkable salinity decrease occurs, leading to a salinity minimum in the late 1960s. The salinity decrease contributes to a decrease of density, which in turn impedes AdDW formation. Such a salinity decrease cannot be simply explained in terms of the BiOS mechanism, by which the salinity fluctuations in the southern Adriatic intermediate layer are attributed to the intrusions of different water masses that are regulated by the NIG circulations. Instead, the change of salt content (SC) of the entire Adriatic Sea is caused by a combined effect of dilution by the net freshwater flux (precipitation + river − evaporation, P + R − E) over the Adriatic Sea and horizontal advection of salinity through the Strait of Otranto (Fig. 7a). In the 1960s, the negative anomaly of the salt content change in the Adriatic Sea is dominated by the strong input of freshwater while the transport of saline Levantine water is less important. Since 1970, the salinity increases in the southern Adriatic Sea (Fig. 6) is then mostly ascribed to the intensive salinity import from the northern Ionian Sea, which outweighs the effect of anomalous freshwater input (Fig. 7a). Strong increase of salinity in the Adriatic starts around 1985, coinciding with the beginning of the period with strong decadal BiOS oscillations.
The horizontal advection of salinity varies well with the NIG circulation patterns, but with an opposite sign (Fig. 7a). Such a relationship has been well explained by the BiOS theory (Gačić et al. 2010). Variations of the net freshwater flux into the Adriatic Sea, mostly attributed to precipitation changes, are out of phase with the winter NAO index (Hurrell 2003) (Fig. 7b). During negative NAO phases, the large-scale atmospheric circulation brings in anomalously high precipitation to the Adriatic Sea, hampering the AdDW formation and thus impeding the reversal of the NIG circulation. The period over the 1960s is quite unique compared to the entire study period as it is characterized by a strong negative NAO signal. Please note that our prescribed river runoff only captures part of the real variability. Due to the lack of observations, approximately 60% of the river runoff into the Adriatic Sea is given by climatological values (see Appendix A). As a result, the variability of the freshwater flux might be underestimated. Nevertheless, we identify a direct link between the NIG pattern and the large-scale atmospheric circulation.

Anticyclonic NIG circulation provoked by a strong AdDW formation event
In order to test the hypothesis (proposed in Sect. 3.2) that NIG reversals are triggered by strong AdDW formation events, we carry out an additional sensitivity simulation. We artificially set the 2-m air temperature and the dew point temperature over the entire Adriatic Sea to zero degree C in January, February, November and December of 1971 and 1972 to mimic two years of extremely harsh winter conditions. We choose these two years for the sensitivity experiment because they are centered in a period  which is rather calm without substantial deep water formation in the Adriatic Sea, and with a persistent cyclonic circulation in the NIG. In the following, the sensitivity simulation is referred to CW ('Cold Winter') while the standard simulation is referred to CTL (Control run). The forcing for the remainder of the simulation is identical to the one used in CTL. The imposed cooling in CW leads to a relatively strong AdDW formation event, which is characterized by a massively dense outflow in 1972 that can be detected down to 1100 m deep in the northern Ionian (Fig. 8b). Consequently, this strong AdDW outflow introduces a change of the NIG circulation from cyclonic to anticyclonic. The anticyclone reaches its maximum in 1973 and starts to weaken in 1974 (Fig. 9). After a short re-adjustment (1974)(1975), the dynamical system in CW approaches the CTL (Figs. 8 and 9). The provoked NIG reversal in CW is also confirmed by the spatial circulation patterns (Figs. 10 and S8 in Appendix E).
Our experimental results resemble the findings from an idealized rotating tank experiment, which shows that the switch of polarity of the near-surface circulation in the NIG can be induced by a mere injection of dense water on a sloping bottom (Rubino et al. 2020).
The comparison of the two simulations demonstrates the mechanisms triggering the NIG circulation reversal by strong AdDW outflows. Comparison of the profiles in a T-S diagram for the first years of CW and CTL illustrates that water masses denser than 28.9 kg m −3 in the northern Ionian Sea (see white box Fig. 1) are replaced by new water masses having lower salinities, indicating the intrusion of newly formed Adriatic dense waters since 1971 (Fig. 11). In the first year of the perturbation (1971), waters in the upper layers are saltier in CW than in CTL. In this year, the anticyclonic circulation has not been established and the AW has not reached the northern Ionian Sea (see Fig.S8 in Appendix E). Instead, an enhanced amount of Levantine-originated water is advected in the upper layers into the southern Adriatic Sea to balance the excessive outflow of AdDW. Since 1972, fresher water masses occupy the upper layers of the water column, indicating that the intrusion of the AW to the northern Ionian Sea is reinforced due to the NIG circulation changes. In the deep layers, the newly formed dense water in CW is mixed with resident water, as the salinity increases along with depth and the properties of the deepest water masses in CW are approaching the properties of the original pre-event deep waters.
Hereafter, the analysis focuses on September, which is in general the month with the most pronounced signal of the anticyclonic pattern (Gačić et al. 2014). In CTL the density at the depth of 671 m does not show a significant horizontal gradient in the northern Ionian Sea (Fig. 12a, c). The isopycnal surfaces are rather flat or show a weak doming in the center (Fig. 13a, c), leading to a cyclonic circulation in the upper layers. In contrast, the density at 671 m depth in CW is subject to substantial differences in both space and time, indicating the spreading of the newly formed AdDW in the northern Ionian Sea. In 1972 in CW, the newly formed dense AdDW flows out of the Adriatic Sea through the western part of the Strait of Otranto and spreads southwestwards along the northwestern flank of the Ionian Sea, whereas there is no AdDW flowing along the northeastern flank    (Fig. 12b). Please note, that the dense water being present at the northeastern flank of the Ionian Sea in this year is not of Adriatic origin and the intrusion of such water mass into the northeastern Ionian Sea is rather small (see upper panels in Fig. S9 in Appendix E). In consequence of the southwestward spreading of the dense AdDW, the isopycnal surface (e.g. 29.06 kg m −3 ) is uplifted at the northwestern flank of the Ionian Sea but shows little changes at the northeastern flank (Fig. 13b). Accordingly, AW is protruding northward into the upper layers of the northern Ionian Sea along the northwestern flank while the anticyclonic flow has not yet been well developed along the northeastern flank (Fig. 10b). In 1973, the dense water found along the northwestern flank of the Ionian Sea starts to spread eastwards through the central Ionian Sea (Fig. 12d, also see lower panels in Fig. S9 in Appendix E), leading to an uplifting of the isopycnal surfaces at the northeastern flank. Eventually a depression of the density distribution in the central part of the transect evolves (e.g. 29 kg m −3 in Fig. 13d). The density gradient in turn facilitates the development of the anticyclonic NIG circulation especially along the northeastern flank (Fig. 10d). The anticyclonic circulation reaches its maximum in 1973, lagging one year behind the strongest AdDW outflow. In CTL, we see a quite similar lag between the strongest AdDW outflow and the timing of the maximum of the anticyclonic NIG circulation. This lag effect has not been identified or discussed in previous studies.

Summary and conclusion
Our long-term hindcast simulation is the first effort to allow investigating the variability of the NIG upper layer circulation over 101 years , with the main focus on earlier periods without good observation coverage. The model reproduces the observed major variability in the EMed thermohaline circulation, e.g. the quasi-decadal reversals of the NIG upper layer circulation since the late 1980s and the evolution of the EMT event in the early 1990s. Consistent realistic atmospheric forcing over the entire simulation period ensures a consistent performance of the model. Hence, our simulation is suitable to explore the hitherto unknown behavior of the NIG circulation before the observational period.
The simulated NIG circulation over 1910-1985 does not oscillate between two opposite modes periodically as proposed by the BiOS mechanism (Gačić et al. 2010). Only in the period of 1940-1960, the variability resembles that observed since the late 1980s, while in the other periods, over 1910-1940 and 1960-1985, the circulation stays in the cyclonic phase with rather weak variability. The BiOS feedback mechanism, implying that internal processes are able to sustain the decadal reversals of the NIG circulation independently from external forcing, is inadequate to explain the type of variability of the two calm periods.
Our results suggest that sufficiently strong AdDW formation events are essential for the reversal of the NIG circulation to the anticyclonic mode. For periods when the Adriatic is not salty enough, strong AdDW formation events are practically not possible and thus the switch to an anticyclonic circulation in the NIG cannot be provoked. During these periods the circulation stays in a weak cyclonic mode because the internal feedback processes are not strong enough to reverse the circulation. Under such low-salinity Adriatic conditions, extremely harsh winters with drastic heat loss are required to provoke a NIG circulation reversal. Our results show that there is a time-lag of 1-2 years between the strongest AdDW outflow and the timing of the maximum of the anticyclonic NIG circulation. The multi-decadal variability in preconditioning the Adriatic water can be explained by the combined effects of net freshwater flux and horizontal advection of salinity. The variations of the net freshwater flux are influenced by multi-decadal changes in the large-scale atmospheric forcing, e.g. the NAO, while the variations of the horizontal advection of salinity are modulated by the NIG circulation patterns.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.