Modes of the BiOS-driven Adriatic Sea thermohaline variability

In this study the impact of the Adriatic-Ionian Bimodal Oscillating System (BiOS) on the interannual to decadal variability of the Adriatic Sea thermohaline circulation is quantified during the 1987–2017 period with the numerical results of the Adriatic Sea and Coast (AdriSC) historical kilometer-scale climate simulation. The time series associated with the first five Empirical Orthogonal Functions (EOFs) computed from the salinity, temperature and current speed monthly detrended anomalies at 1-km resolution are correlated to the BiOS signal. First, it is found that the AdriSC climate model is capable to reproduce the BiOS-driven phases derived from in-situ observations along a long-term monitoring transect in the middle Adriatic. Then, for the entire Adriatic basin, high correlations to the 2-year delayed BiOS signal are obtained for the salinity and current speed first two EOF time series at 100 m depth and the sea-bottom. Finally, the physical interpretation of the EOF spatial patterns reveals that Adriatic bottom temperatures are more influenced by the dense water circulation than the BiOS. These findings confirmed and generalized the known dynamics derived previously from observations, and the AdriSC climate model can thus be used to better understand the past and future BiOS-driven physical processes in the Adriatic Sea.


Introduction
A decade ago, Gačić et al. (2010) proposed the first explanation of the physical mechanisms driving the Adriatic-Ionian Bimodal Oscillating System or BiOS, more than half a century after the existence of quasi-decadal variability of salinity and water masses in the Adriatic was noticed (Buljan 1953;Buljan and Zore-Armanda 1976). This explanation provides a scientific framework which connects the quasidecadal reversals of the Northern Ionian Gyre (NIG) circulation to the salinity variability and the dense water formation in the Adriatic Sea (Fig. 1). Namely, the southern Adriatic Sea salinity is decreased by the advection of the less-saline Modified Atlantic Water during the anticyclonic phase of the NIG and increased by the advection of the highly-saline Levantine/Eastern Mediterranean waters during the cyclonic phase of the NIG. Additionally, the reversals from anticyclonic to cyclonic (or cyclonic to anticyclonic) of the NIG are explained by the decreased (or increased) density of the Adriatic dense water outflow in the northern Ionian Sea, where the isopycnal surfaces are deepened (or raised) and thus the upper water layer is stretched (or squeezed) leading to a change in vorticity and hence a NIG reversal. Later, the BiOS regimes have also been connected to the salinity variability and the dense water formation of the Aegean basin Theocharis et al. 2014;Velaoras et al. 2014;Reale et al. 2017) and wind stress curl variability over the Ionian Sea (Pinardi et al., 2015).
On the one hand, the BiOS phases in the northern Ionian have been widely studied based on in-situ and remote sensing observations (e.g. Malanotte-Rizzoli et al. 1997;Larnicol et al. 2002;Pujol and Larnicol 2005;Borzelli et al. 2009;Bessières et al. 2013;Gačić et al. 2010Gačić et al. , 2014. However, due to the well-documented difficulties that Regional Climate Models (RCMs) have to capture the hurricane-strength bora events driving the dense water formation in the northern Adriatic Dunić et al. 2018;Denamiel et al. 2021a), the BiOS reversal mechanism-as defined by Gačić et al. (2010)-remained partially unproven till recently. Indeed, with a rotating tank experiment, Rubino et al. (2020) and Gačić et al. (2021) demonstrated that the injection of dense water on a sloping bottom could trigger the reversal of their near-surface gyre circulation. Further, the 100-year long realistic numerical simulation of Liu et al. (2021) reproduced the known BiOS phases between 1985 and 2010, despite not properly capturing the highest salinities and densities measured in the northern Ionian. By using an artificial cooling over the entire Adriatic Sea, Liu et al. (2021) also confirmed the results obtained by Rubino et al. (2020), and provided evidences that the NIG reversal from a cyclonic to an anticyclonic phase occurs 1-2 years after a major cooling event driving the Adriatic dense water formation.
On the other hand, many of the processes observed in the Adriatic Sea have been connected and even correlated to the decadal variability of the BiOS. For example, Civitarese et al. (2010) demonstrated the impact of the BiOS on the nutrient content in the Adriatic Sea (i.e. decrease/increase of the nitracline depending on the cyclonic/anticyclonic circulation of the NIG); Lavigne et al. (2018) explained the phytoplankton phenology in terms of the decadal reversals of the Ionian circulation and Mihanović et al. (2015) extracted BiOS-driven phases from long-term measurements of salinity, temperature and dissolved oxygen in the middle Adriatic. The abundance and species within the Adriatic zooplankton communities have also been connected to the BiOS by Batistić et al. (2014). More recently, Peharda et al. (2018) found a consistent anti-correlation between the bivalve chronology in the northern Adriatic Sea and the BiOS regimes while Vilibić et al. (2020) and Ciglenečki et al. (2020) strongly correlated bottom salinity and the surface active substances fraction of the dissolved organic carbon in the northern Adriatic to the delayed BiOS index. However, no numerical study has ever quantified the impact of the BiOS regimes on the decadal variability of the Adriatic Sea as the available long-term ocean runs were not properly simulating the dense water formation, spreading and storage within the Adriatic basin .
In this study, the newly evaluated 31-year long ocean historical results (Pranić et al. 2021) of the kilometer-scale Adriatic Sea and Coast (AdriSC) climate model-which has demonstrated a decent level of skills in reproducing the Adriatic overall thermohaline properties and dynamics, including generation of dense waters-are thus used to connect the BiOS phases to the long-term variability of the Adriatic Sea. The presented work is structured as follows. First, the AdriSC model, the available observations along two long-term monitoring transects in the Adriatic Sea and the methods used hereafter are briefly presented in Sect. 2. Then, in Sect. 3, along the long-term monitoring transects, the AdriSC results are evaluated against the observations and the modelled BiOS-driven phases are compared with those already derived in Mihanović et al. (2015). This section is continued by connecting and interpreting the temporal and spatial patterns of the AdriSC model results to the BiOS regimes. To conclude in Sect. 4, the validity of the presented methods is discussed and the obtained results are further connected to the previous findings derived, in the last decade, from observations or simple model experiments.

AdriSC climate model
The Adriatic Sea and Coast modelling suite (AdriSC; Denamiel et al. 2019a) has been recently developed with the aim to study atmospheric and oceanic processes in the Adriatic region at various temporal and spatial scales ranging from the impact of climate change on extreme events (Denamiel et al. 2020a(Denamiel et al. , b, 2021a to the operational forecast of extreme sea-levels along the Croatian coasts (Denamiel et al. 2019b;Tojčić et al. 2021). The modelling suite is built around three main modules that can be used either together or separately, depending on the type and scale of studies performed in the Adriatic region. First, the basic module uses a modified version of the Coupled Ocean-Atmosphere-Wave-Sediment-Transport modelling system (COAWST; Warner et al. 2010) and produces atmospheric results  Skamarock et al. 2005) at up to 3-km resolution and ocean results with the Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams 2005) at up to 1-km resolution (hereafter referred as AdriSC ROMS 1-km; Fig. 1). Second, the extreme event module is coupling offline the WRF model at 1.5-km resolution with the unstructured barotropic ADvanced CIRCulation model (ADCIRC; Dietrich et al. 2012) at up to 10 m resolution along the Adriatic coasts. And, third, the stochastic module based on uncertainty quantification methods (Denamiel et al. , 2020c is providing, along the Croatian coast, extreme sea-level hazard assessments during meteotsunami events (i.e. tsunami-like waves driven by atmospheric disturbances).
Within the AdriSC modelling suite, the climate component-based solely on the basic module-is dedicated to the study of long-term kilometer-scale processes occurring in the Adriatic region. The complete description of the AdriSC climate component, including the set-up of the atmospheric and oceanic models, is presented in Denamiel et al. (2020aDenamiel et al. ( , 2021b and Pranić et al. (2021). Additionally, the performances of the AdriSC climate component have already been thoroughly assessed for both the atmosphere (Denamiel et al. 2021b) and the ocean (Pranić et al. 2021). Overall, the evaluation revealed that the AdriSC climate component is better suited to study the long-term processes over the Adriatic region than the classical Regional Climate Models (RCMs) used in the Mediterranean Sea. In this study, the daily results from the evaluation run performed with the AdriSC climate component and spanning over a 31-year long period  are thus analyzed. These results can be accessed and retrieved via the web interface https:// vrtlac. izor. hr/ ords/ adrisc/ inter face_ form (Ivanković et al. 2019;Denamiel et al. 2021b).

Data and methods
Even though the AdriSC ROMS 1-km performances have already been assessed by Pranić et al. (2021), temperature and salinity from in situ temperature and salinity qualitychecked observations collected along the long-term monitoring Palagruža Sill transect (by the Institute of Oceanography and Fisheries, Croatia) and the northern Adriatic transect (by the Ruđer Bošković Institute, Croatia) are used to evaluate the thermohaline properties of the model (Fig. 1). These two datasets present the advantage of having been collected regularly (i.e. nearly every month since 1979 for the northern Adriatic transect and twice a year to seasonally since 1952 for the Palagruža Sill transect) and sampled at predefined stations (i.e. 6 for the northern Adriatic transect and 5 for the Palagruža Sill transect; Fig. 1) for standard depths (i.e. 0, 5, 10, 20, 30, 35 m depth for the northern Adriatic transect and 0, 10,20,30,50,75,100,120,150, 170 m depth for the Palagruža Sill transect). These transects have been established as being the most convenient for maintenance by the monitoring institutions (that are positioned at the transect ends), while covering two important areas for the Adriatic basin-wide dynamics. Namely, the northern Adriatic transect is overlaying the dense water formation site (Bergamasco et al. 1999;Janeković et al. 2014;Vilibić et al. 2019) where the densest Mediterranean waters are being generated (Zore 1963), bringing oxygen to the bottom layers and driving the thermohaline circulation of the Adriatic-Ionian region (Orlić et al. 2006). The Palagruža Sill transect is covering a key region for the exchanges of water masses along the Adriatic Sea, and is a passage where quasi-permanent topographically-driven Southern Adriatic Gyre and Middle Adriatic Gyre are exchanging water masses (Artegiani et al. 1997;Martin et al. 2009). Additionally, dense water outflow from the northern Adriatic can be traced at the deepest and southern parts of the sill (Artegiani et al. 1987), while the acrosssill transport driven by the Western Adriatic Coastal Current is confined to the western coast ). Consequently, the evaluation of the AdriSC ROMS 1-km salinity and temperature daily results along these transects provides a representative picture of the quality of the model during the 1987-2017 period, which is relevant for a study aiming to connect the long-term variability of the Adriatic Sea with the northern Ionian BiOS.
Except for the evaluation along the long-term monitoring transects for which the AdriSC ROMS 1-km model results are directly compared to the observations, the hereafter presented sea-level, temperature, salinity and current speed analyses are all performed with the monthly means of the detrended daily anomalies. Along the two transects, the AdriSC 1-km model results are extracted at regular depths (i.e. every meter for the northern Adriatic transect and every 5 m for the Palagruža Sill transect) as linear segments defined by two stations. The results are then presented as one single section oriented from west to east for the northern Adriatic transect and from south to north for the Palagruža Sill transect. Finally, long-term variabilities are derived from Empirical Orthogonal Functions (EOFs) which are used to compare, in space and time, the most important variability patterns in the Adriatic and northern Ionian seas. Indeed, Gačić et al. (2010Gačić et al. ( , 2011Gačić et al. ( , 2014 have demonstrated that the BiOS-consisting in the decadal switch from cyclonic to anticyclonic of the circulation in the northern Ionian Sea (and vice versa)-is well described with the change of sign of one of the main EOF components derived from sea-level products in the Ionian Sea. In this study, for consistency, the BiOS signal is derived from the sea-surface height field (SSH) of the MEDSEA reanalysis (Simoncelli et al. 2019) distributed by the Copernicus Marine Environment Monitoring Service (CMEMS). MEDSEA reanalysis was also used to force the open boundaries of the AdriSC ocean model in the evaluation run. All the presented spatial EOFs (also known as Principal Component Analysis or Eigen Analysis) are obtained via a covariance matrix and are normalized (i.e. the sum of squares for each EOF pattern equals one). The time series of the amplitudes (also known as principal components or expansion coefficients) associated with each eigenvalue in the EOF are derived via the dot product of the data and the EOF spatial patterns, and the mean is subtracted from the value of each component time series. Consequently, the BiOS signal derived from one of the MEDSEA sea-level EOFs can be correlated with the long-term variability of temperature, salinity and current speed given by the EOF time series derived from the AdriSC 1-km model results. All presented correlations are significant with levels below 5% following the null-hypothesis and the statistical significance of the correlation is conducted by a Monte Carlo simulation using random phase resampling (Ebisuzaki 1997). It should be noted (1) that all correlations and time series comparisons between the BiOS signal and the long-term Adriatic Sea thermohaline variabilities have been made for signals smoothed with a 3-month running mean and (2) that, for visual purposes, EOF time series have been rescaled in order to have similar amplitudes in the comparison plots.

Evaluation along the long-term monitoring Adriatic transects
First, the performances of the AdriSC ROMS 1-km model along the long-term monitoring northern Adriatic and Palagruža Sill transects are investigated with scatter plots (Fig. 2). Overall, they show that the hexagons with the largest number of points are following the reference line for both temperature and salinity, which indicates that the vast majority of the AdriSC ROMS 1-km results corresponds well to the observations in both intensity and timing. However, for the northern Adriatic transect, the model systematically overestimates the salinity observations below 37. This highlights that the AdriSC ROMS 1-km model is not accurately reproducing the intensity, the complex variability, and the extent of the Po River fresh water plume known to strongly influence the thermohaline circulation in the northern Adriatic Sea (Kourafalou 1999; Manzo et al. 2018). It should also be highlighted that, along the northern Adriatic transect, the temperatures above 20 °C (i.e. temperatures in the surface layer modelled during summer) are also systematically underestimated by up to 2.5 °C while they can be overestimated by up to 5 °C between 10 and 20 °C. The AdriSC ROMS 1-km has been demonstrated to present a cold summer temperature bias linked to various factors such as a deficit of solar radiation by the AdriSC atmospheric model during the summer, or the fact that the river temperatures are imposed by taking the ERA-Interim skin temperatures the closest to the river estuaries (Pranić et al. 2021). Additionally, the northern Adriatic transect is located in shallow waters (below 40 m depth) where the optical properties of the water, not well parameterized in the Adriatic Sea, are playing a crucial role in modelling the turbidity and hence vertical mixing. The turbidity is responsible for most of the downward shortwave radiation absorption in the upper layer and hence highly influences all the temperatures along the transect. By contrast, for the Palagruža Sill transect, despite some minor scattering for few occurrences up to ± 0.25 for the salinity and ± 2 °C for the temperature, the AdriSC ROMS 1-km model is in good agreement with the observations and thus provide reliable results concerning the thermohaline circulation.
Second, as the AdriSC ROMS 1-km model is showing some skill in representing the thermohaline properties along the two long-term monitoring Adriatic transects, an EOF analysis is performed in order to find potential long-term variabilities. For this section, the five main EOF components (representing the highest percentages of the signal) of salinity and temperature monthly detrended anomalies are extracted and presented along both transects in Fig. 3 for the first component (hereafter referred as EOF 1) while the remaining components (EOFs 2 to 5) are not shown in this study. Along the northern Adriatic transect, for both temperature and salinity, none of the EOF components captures interannual to decadal variabilities. This predominant lack of long-term variability can be probably explained by the influence of the Po River plume, which affects the full transect, even the eastern stations (Vilibić et al. 2019). For example, dissolved oxygen and organic carbon measured at the easternmost station have been strongly correlated to the long-term variability of the Po River discharge Ciglenečki et al. 2020). The impact of the Po River plume is also seen on the spatial patterns of the salinity EOF 1 (66.2% of the signal) which highlight a difference in the strength of EOF signals between the upper layer up to 10 m depth and the rest of the profile (i.e. presence in the upper layer of a stronger signal for salinity along almost entire northern Adriatic transect). For the temperature, EOF 1 (79.5% of the signal) contains no decadal variability, reflecting the fact that the northern Adriatic transect is located along a shallow shelf (at less than 40 m depth), highly affected by and rapidly responding to local atmospheric processes and tides and the respective heat distribution within the water column.
By contrast, the Palagruža Sill transect located in the middle of the Adriatic Sea in deeper waters less influenced by river discharges, tides and atmospheric conditions, is more likely to capture long-term variabilities. Still, for the AdriSC ROMS 1-km temperature results, none of the EOF components displays any decadal variability even though EOF 1 (representing 58.4% of the signal) presents some interannual variabilities, particularly strong in the northern shallower part of the transect. However, for the AdriSC ROMS 1-km salinity results, the time series of amplitude associated with the EOF 1 (representing 78.5% of the signal) clearly shows some well-defined interannual to decadal oscillations. Additionally, the EOF 1 spatial patterns highlight some differences: the weakest signal can be seen in the deepest part of the transect while the strongest signal is present in the northern shallower part. More importantly, the oscillations obtained with the AdriSC ROMS 1-km salinity results are in good agreement with the BiOS driven phases defined by Mihanović et al. (2015) with a Self-Organizing Maps (SOM) method applied to temperature, salinity and dissolved oxygen observations along the Palagruža Sill transect during the 1952-2010 period. Indeed, the mostly negative amplitudes are obtained for the 1987-1990, 1999-2006 and 2012-2017 periods instead of the 1987-1990 and 1999-2005 periods found in Mihanović et al. (2015). And, the mostly positive amplitudes are occurring for the 1991-1998-2011periods instead of 1991-1996periods described in Mihanović et al. (2015. The fact that certain phases are shifted by a year or two compared to Mihanović et al. (2015) can be largely attributed to their definition of intermittent phases for the 1997-1998 and 2006-2008 periods. This definition is directly linked to their interpretation of the results derived with the SOM method and does not apply to the EOF method presented in this study nor, more generally, to the BiOS indices (Gačić et al. 2010(Gačić et al. , 2014Civitarese et al. 2010;Liu et al. 2021). Further, Mihanović et al. (2015) used an observational dataset stopping in 2010 and the phases obtained with the SOM method may have been slightly shifted if the records extended till 2018 (as seen in the sensitivity experiments provided Fig. 5 of their study).
In brief, while the AdriSC ROMS 1-km model struggles to reproduce the complex variability of the Po River plume influencing the thermohaline properties of the entire long-term monitoring northern Adriatic transect, it performs well along the Palagruža Sill transect where salinity and temperature scatter has been shown to be minimal. Additionally, despite not seeing substantial long-term variability along the northern Adriatic transect due to the dominance of local processes such as Po River plume spreading, the EOF analysis of the AdriSC ROMS 1-km salinity results along the Palagruža Sill transect does reproduce the phases previously extracted from the observations by Mihanović et al. (2015). Further, as both Modified Atlantic Water and Levantine/Eastern Mediterranean waters are mainly characterized by their salt content, the BiOS signal is thus expected to be principally seen in the salinity EOFs. Consequently, both the AdriSC ROMS 1-km results and the EOF methodology presented in this study are reliable enough to be applied over the entire Adriatic Sea.

Interannual to decadal variability of the Adriatic Sea thermohaline circulation
In this section, the AdriSC ROMS 1-km salinity, temperature and current speed monthly detrended anomalies are extracted over the entire Adriatic Sea at the sea-surface, the sea-bottom and at 100 m depth. First, the five main EOF components (representing the highest percentages of the signal) derived from these results are analyzed by correlating the Ionian BiOS signal (hereafter simply referred as BiOS signal) with the time series of EOF amplitude (Table 1 and Fig. 4).
In this study, the temporal variability of the BiOS signal is defined as the second normalized spatial EOF component and associated time series of amplitude derived from the MEDSEA SSH fields in the Ionian Sea during the 1987-2017 period (top panels, Fig. 4). The spatial extent as well as the interannual to decadal variabilities of the obtained signal are in good agreement with indices calculated from altimetry data (Gačić et al. 2010(Gačić et al. , 2014, cruise observations , and model results (Liu et al. 2021). More precisely, the anticyclonic phases (positive sign of the EOF amplitude) are clearly seen for the 1987-1997 and 2006-2009 periods, while the cyclonic phases (negative sign of the EOF amplitude) are present for the 1998-2005 and 2010-2017 periods. However, small discrepancies concerning the timing of the reversals exist between the different BiOS indices. For example, the obtained 2006-2009 anticyclonic phase is 1 year shorter compared than the one derived by Gačić et al. (2010) and Bessières et al. (2013). In fact, uncertainties concerning the precise timing of the reversals are known to be introduced by the presence of mesoscale eddies during the transition periods of the BiOS (Gačić et al. 2014;Liu et al. 2021) and different methods can give slightly shifted phases. Consequently, it can be safely concluded that the BiOS regimes are appropriately reproduced by the second normalized EOF of the MEDSEA SSH fields in the Ionian Sea.
The correlation coefficients obtained between the BiOS signal and the time series of temperature, salinity and current speed EOF amplitudes at the surface, 100 m depth, and the bottom of the Adriatic Sea are presented in Table 1 for the BiOS signal with 0-year to 3-year lags. The re-scaled time series of the BiOS signal on top of the re-scaled EOF amplitudes with the highest correlation coefficients and representing the highest percentages of the signal are also plotted in Fig. 4 (bottom panels). For the temperature, the highest correlations are obtained between the BiOS signal with 1-year lag and EOF 2 amplitudes at the sea-surface (0.30) and at 100 m depth (0.58) as well as EOF 3 amplitude at the bottom (0.57). For the salinity, the highest correlations are obtained between the BiOS signal with 2-year lag and EOF Table 1 Correlation coefficients between the time series of normalized EOF amplitudes of the BiOS signal extracted from MEDSEA (EOF 2 of SSH) with 0-year to 3-year lags and the temperature, salinity, current speed signals at the surface (Surf.), the bottom (Bot.) and 100 m depth (100 m) extracted from the AdriSC ROMS 1-km model (EOFs 1-5) Non-significant correlations following the null hypothesis are marked with the "/" symbol. The highest correlations are highlighted in bold Correlation analyses are also performed between the time series of temperature and salinity EOF amplitudes extracted along the long-term monitoring northern Adriatic and Palagruža Sill transects and the BiOS signal with 0-year to 3-year lags, and are presented in Table 2. As before, the only meaningful correlation is found between the BiOS signal with a 2-year lag and the salinity EOF 1 amplitude time series along the Palagruža Sill transect (i.e. correlation coefficient of 0.62). The fact that the Adriatic Sea haline response is delayed by approximately 2 years compared to the BiOS reversals sheds some light on the previously presented results.
For example, the salinity EOF 1 amplitude sudden and short reversal from negative to positive seen in 2014-2015 along the Palagruža Sill (bottom right panel, Fig. 3) can be associated to the observed premature inversion of the BiOS signal in late 2012 caused by substantial generation of dense waters (Gačić et al. 2014). Other sudden and short reversals from negative to positive or variability of the salinity EOF 1 amplitude along the Palagruža Sill transect could also be the 2-year delayed response to BiOS signal variability driven by extreme bora events, like the ones documented between 1998 and 2017 in Denamiel et al. (2021a). However, it should be noted that these reversals are not present in the obtained BiOS signal (top right panel, Fig. 4). Yet the 2000-2001 and 2013-2014 periods are marked with a pronounced increase of the SSH EOF amplitude, from approximately − 6.0 to − 1.5 and − 4.0 to − 0.25. This is probably linked to the lack of accuracy of the atmospheric forcing used to produce the MEDSEA fields. Indeed, the change of vorticity in the northern Ionian Sea leading to the change of sign of the SSH EOF highly depends on the dense water travelling through the strait of Otranto and resulting from the intensity of the bora-induced cooling over the Adriatic Sea.
At this point it must also be mentioned that correlations should be interpreted with some caution. As seen in Sect. 3.1, the temperature EOFs did not show any decadal variability along the long-term monitoring transects and thus may not be such a strong indicator of the Adriatic Sea response to the BiOS signal, contrarily to the salinity EOFs. Further, as mentioned before, both Modified Atlantic Water and Levantine/Eastern Mediterranean waters-which are the main tracers of the BiOS signal phases in the Adriatic Sea-are mainly characterized by their salt content. That can be assessed also from much lower percentages of explained variance in temperature vs. salinity for different EOFs correlated with the BiOS signal: EOF 2 in temperature (18.3%) vs. EOF 1 in salinity (60.8%) at 100 m and EOF 3 in temperature (6.5%) vs. EOF 1 in salinity (37.1%) at the bottom. Indeed, the wintertime cooling, in particular during severe bora events (Janeković et al. 2014;Denamiel et al. 2021a), may substantially uptake the heat from the deep Adriatic (i.e. up to 900 m) to the atmosphere (e.g. Gačić et al. 2002;Cardin et al. 2020). Consequently, temperature might not be as a conservative tracer in the deep Adriatic as salinity. Coming from correlations with salinity EOFs, the response time of the Adriatic thermohaline circulation to the BiOS phases is thus likely to be of approximately 2 years.
As one of the major drawbacks of any EOF analysis is the potential lack of physical meaning of the obtained components, the spatial patterns of the temperature, salinity and current speed EOF components with the highest correlation to the BiOS signal (including a 2-year lag) are now analyzed with the aim to connect them to the known dynamical properties of the Adriatic Sea (Figs. 5, 6, 7).
At the sea-surface where the influence of the atmospheric forcing and the river fresh water discharges are the highest (Fig. 5) and masking the BiOS signal, the time series of the EOF amplitudes tend to display important interannual variabilities for temperature, salinity and current speed. Additionally, the associated spatial patterns highlight some major differences. First, the spatial patterns of the temperature EOF 2 (representing only 5.5% of the signal) is strongly negative in the south and strongly positive in the north and associated with time series where the delayed BiOS reversals cannot be easily distinguished. Second, the spatial patterns of the salinity EOF 3 (representing only 9.5% of the signal) are homogeneously negative over the entire Adriatic Sea, except along the coast where the influences of the river discharges are the strongest (the Po River plume, other northern Adriatic rivers, Neretva river, Albanian rivers). Such a multipolar structure might resemble river plumes being Comparison of the re-scaled BiOS signal with the time series of rescaled normalized EOF amplitudes-with the highest correlations to the BiOS signal including a 1-year or 2-year lag-derived during the 1987-2017 period from the AdriSC ROMS 1-km temperature, salinity and current speed monthly anomalies at the surface, the bottom and 100 m of depth (3 bottom panels) ◂ more confined towards river mouths and coastlines during high salinity conditions in the Adriatic and vice versa. Third, the spatial patterns of the current speed EOF 4 (representing only 3.6% of the signal) are extremely complex and the associated time series disconnected to the BiOS signal. Still, it seems that the BiOS sporadically influences the size of the South Adriatic Gyre, which can be confined either closer to the coastlines or to the center of the pit as the result of different BiOS regimes, while also affecting the transport rates in the Strait of Otranto.
At 100 m depth (Fig. 6) where, contrarily to the seasurface, the thermohaline circulation is dominant, (1) the obtained EOFs all represent a large or major part of the signal (18.3% for the temperature, 60.8% for the salinity and 9.0% for the current speed), (2) the associated time series all display the 1 year (for temperature) or 2 year (for salinity and current speed) delayed response in reaching maximum correlations to the BiOS signal and (3) the spatial patterns with a mostly negative signal are all nearly homogeneous over the entire Adriatic Sea, in particular for the salinity. Additionally, for all the variables, the weakest signal is found in the middle of the South Adriatic Pit which seems less influenced by the BiOS signal than the rest of the domain.
Indeed, the center of the South Adriatic Pit is the center of the cyclonic gyre, where upwelling and deep convection take place (i.e. vertical processes are bringing more water masses from the deep and surface Adriatic than at the perimeter of the gyre) while being somehow separated from the advected waters. In contrast, the perimeter of the South Adriatic Pit is characterized by a persistent and strong current, bringing waters from the Ionian Sea and keeping about 90% of them in the loop (Gačić et al. 2014). This can be seen in the spatial patterns of the current speed EOF 2 where the South Adriatic Gyre is highly sensitive to the BiOS signal despite the associated time series of amplitudes displaying a BiOS related signal strongly embedded in the interannual variabilities. Further, the spatial EOF of the temperature at 100 m clearly shows a distinct area with positive signal in the western part of the Jabuka Pit and southern Palagruža Sill which are known collectors of the dense waters formed during extreme bora events. It should also be noted that the brief inversions of the BiOS signal, already mentioned for the salinity EOF analysis along the Palagruža Sill transect, can clearly be seen in the time series of both temperature and salinity EOF amplitudes. At 100 m depth, it even seems that the temperature inversions are stronger than the ones of the salinity.
Finally, at the sea-bottom (Fig. 7), the obtained EOFs most correlated to the BiOS signal are spatially homogeneous (negative over the entire Adriatic Sea) and represent a large part of the signal (more than 30%) for the salinity and the current speed. Precisely, salinity EOF 1 spatial amplitude (accounting for 37.1% of the signal) is the largest in the northern Adriatic and along the Western Adriatic Coastal Current. Such a strong effect of the delayed BiOS signal to the northern Adriatic bottom salinity indicates that wintertime cooling and dense water production are preconditioned by the salinity content driven by the BiOS. Interestingly, this even affects the shallow northern Adriatic somehow contrasting with the findings along the northern Adriatic transect indicating no BiOS driven salinity variability (Figs. 3). Yet, the transect analyses included the surface layer strongly shaped by local processes and masking the effects of the basin wide processes like the BiOS. Regarding the bottom temperature (Fig. 7), EOF 3 is correlated to the delayed BiOS phases. But it represents only 6.5% of the temperature signal and, consequently, the local cooling due to extreme bora events probably affects the bottom temperature changes much more than the advection of waters driven by BiOS. Still, the EOF 3 spatial patterns present an interesting distribution. They are strongly negative along the south-and middle-eastern Adriatic coast till the Kvarner Bay and along the south-western Adriatic. But they are positive in areas where the northern Adriatic dense waters are known to be either (1) formed during extreme bora events like the northern Adriatic shelf (Bergamasco et al. 1999) or the Kvarner Bay (Janeković et al. 2014;Vilibić et al. 2018), (2) travelling towards the northern Ionian Sea-i.e. along the Italian coast (Artegiani et al. 1987;Vilibić and Mihanović 2013) or (3) collected within the Adriatic Seai.e. the Jabuka Pit (Marini et al. 2006) and the South Adriatic Pit (Querin et al. 2016). Therefore, this EOF is describing the connection between the BiOS signal and the Adriatic dense water formation and spreading, which is not as intuitive for temperature as for salinity. Precisely, the cyclonic BiOS signal is resulting in an advection of warmer waters, reaching maximum correlations after 2 years, while simultaneously connected with colder-than-usual winters that produced dense waters with lower temperatures. In fact, such a coincidence between these two unrelated processes is not frequent in the Adriatic-for that reason accounts just for a small percentage of EOF solutions-but still recognized as the third most frequent EOF driving the bottom temperature. For the current speed, the EOF 1 spatial patterns (consisting in 30.5% of the signal) are similar to the one described at 100 m depth which highlights that the South Adriatic Gyre is mostly affected by the 2-year delayed BiOS signal below 100 m depth. In summary, except at the sea-surface where the atmospheric-and river-driven processes are dominant, the Adriatic Sea thermohaline circulation is well correlated with the 2-year delayed BiOS phases derived from the sea-surface height variability in the Ionian Sea. Overall, salinity is found to be the best indicator of the response

Discussions and conclusions
The principal novelty of the presented study consists in using, for the very first time, a kilometer-scale climate model to investigate, via an EOF and correlation approach, the impact of the BiOS on the Adriatic Sea thermohaline circulation at different depths over the entire basin. The main findings are twofold. On the one hand, the AdriSC ROMS 1-km model is proven to be, to this date, the only numerical model capable to reproduce the BiOS-driven phases observed along the Palagruža Sill long-term monitoring transect on a climate scale. On the other hand, over the entire Adriatic basin, the BiOS signal is demonstrated to be only weakly correlated to the sea-surface circulation and the two main temperature EOF time series (except at 100 m depth with a 1-or 2-year lag), but better correlated with a 2-year lag to the salinity and current speed two main EOF time series at 100 m depth and the sea-bottom. However, due to the uncertainties associated with the numerical results and the methods used in this work, the validity of this approach is further discussed hereafter.
First, EOF methods have been widely used within the climate community in order to identify the modes of variability and the predictability of the Earth system (Navarra and Simoncini 2010). Some known examples are the extraction of the El Niño Southern Oscillation (ENSO) interdecadal variabilities by Zhang et al. (1997) and Newman et al. (2003) or the study of the Atlantic thermohaline circulation variability by Hawkins and Sutton (2007). However, EOF analyses suffer from major drawbacks such as the orthogonality of the EOF patterns and the uncorrelatedness of the associated time series (e.g. Jolliffe 2002;Hannachi 2007;Monahan et al. 2009). Consequently, the interpretation of EOF results in physical terms is not necessarily straightforward, except when an outstanding mode is present (e.g. ENSO in global sea-surface temperature or North Atlantic Oscillation in mean sea-level pressure). In the presented results of the Adriatic thermohaline circulation response to the BiOS regimes, such a predominant mode exists particularly for the salinity and current speed results at 100 m depth and the sea-bottom. Consequently, the correlation between the EOF time series obtained in the Adriatic Sea and the BiOS signal extracted in the northern Ionian Sea provides relevant information concerning the BiOS mechanisms, with relatively high values obtained for the correlation coefficients (i.e. maximum value of 0.66).
Second, during the studied 1987-2017 period, the BiOS reversals are found to impact the Adriatic thermohaline circulation with approximately 2-year delay. Mihanović et al. (2015) observed similar response time during 1993-2014 period, based on the long-term monitoring Palagruža Sill transect measurements and satellite altimetry data. Additionally, the 2-year delay corresponds the decay time of 26 months of the Adriatic Deep Water estimated by Vilibić and Orlić (2002) with a simple box model. Thus, it may be the natural response time of the southern Adriatic basin to a change in water masses, which was proving to be driven by the BiOS in this study. However, Mihanović et al. (2015) also highlighted that during exceptional conditions such as the Eastern Mediterranean Transient (EMT; Roether et al. 2007) in the 1990s, the delay may be reduced to 1 year, while prolonged to more than 2 years during slow BiOS reversals. It should be noted that EMT conditions were extraordinary during this period as the dense water flow of Aegean origin was almost an order of magnitude larger than the dense water flow of Adriatic origin. Consequently, these conditions resulted in the strongest anticylonic phase of the BiOS ever recorded.
Finally, spatial results of the EOFs obtained for the AdriSC ROMS 1-km salinity, temperature and current speed, clearly display three different types of circulation. At the surface, the circulation is nearly totally disconnected from the BiOS signal. This can be explained by the known influence of the temporal changes of both the atmospheric forcing (Béranger et al. 2010;Janeković et al. 2014) and the river freshwater discharges (Orlić et al. 1992;Lazar et al. 2007) on the Adriatic surface thermohaline circulation. In intermediate and deep layers-here represented by EOFs at 100 m and the bottom layer-the haline circulation displaying strong interannual and decadal oscillations associated with nearly homogenous spatial patterns is known to be driven mostly by the BiOS and the EMT (Roether et al. 2007;Gačić et al. 2010Gačić et al. , 2014. Further, the bottom temperature EOF 3 spatial patterns-interpreted as dense water formation, spreading and storage locationswere connected with a delayed BiOS signal. In the Adriatic basin, dense waters were observed at the exact locations (Vilibić et al. 2004;Wang et al. 2006) and are known to drive the deep Adriatic thermohaline circulation (Orlić et al. 2006;). However, this signal is an order of magnitude less pronounced than the BiOS-induced intermediate and deep salinity signal, nominating the salinity as the major tracer of the incoming Adriatic water masses.
To conclude, the analysis of the AdriSC ROMS 1-km salinity, temperature and current speed EOF modes has confirmed and generalized the known influence of the BiOS on the Adriatic Sea thermohaline circulation derived from in-situ observations. Consequently, further analysis of the model results can now be performed in order to study the BiOS-driven physical processes within the Adriatic Sea as well as the effect of the Adriatic dense water on the NIG reversals. Additionally, the impact of extreme warming (i.e. Representative Concentration Pathway or RCP 8.5 scenario) on the Adriatic-Ionian BiOS regimes can also be derived from the recently completed AdriSC simulation for the 2070-2100 period following the Pseudo-Global Warming (PGW) methodology presented in Denamiel et al. (2020a). Kilometer-scale climate modelling in the Adriatic basin has thus open the door to a better understanding of the complex processes driving the atmospheric and oceanographic baroclinic circulation but also to a better assessment and, eventually, mitigation of the coastal hazards under past, present and future climate warning conditions. Acknowledgements The contributions of the Institute of Oceanography and Fisheries (Croatia) as well as the Ruđer Bošković institute (Croatia) which provided the observations along the two Adriatic long-term monitoring transects are greatly appreciated. The contributions of the European Centre for Middle-range Weather Forecast (ECMWF) which provided the computing and archive facilities used in this research as well as the support of their staff are also greatly appreciated.
Authors contributions IV and CD contributed to the study conception and design. Material preparation was done by IT, PP and CD. Set-up of the model and simulations were performed by CD. Analysis of the results and production of the figures were performed by IV and CD. The first draft of the manuscript was written by CD and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Conflict of interest
The authors declare that they have no conflict of interest.

Data availability
The model results and the programs used to produce this article can be obtained under the Open Science Framework (OSF) FAIR data repository https:// osf. io/ 7em8h/ (https:// doi. org/ 10. 17605/ OSF. IO/ 7EM8H).