Multi-model analysis of the Adriatic dense water dynamics

10 This study aims to enhance our understanding of the bora-driven dense water dynamics in the Adriatic Sea using different state-of-the-art modelling approaches during the 2014-15 period. Practically, we analyse and compare the results of four different simulations: the latest reanalysis product for the Mediterranean Sea, a recently evaluated fine-resolution atmosphere-ocean Adriatic Sea climate model and a long-time running Adriatic Sea atmosphere-ocean forecast model used in both hindcast and data assimilation (with 4-day cycles) modes. As a first step, we evaluate the resolved physics in each simulation by 15 focusing on the performance of the models. Then, we derive the general conditions in the ocean and the atmosphere during the investigated period. Finally, we analyse in detail the numerical reproduction of the dense water dynamics as seen by the four simulations. The likely prerequisites for proper modelling of the ocean circulation in the Adriatic basin, including kilometre-scale atmosphere-ocean approach, non-hydrostatic atmospheric models, fine vertical resolutions in both atmosphere and ocean and the location and forcing of the open boundary conditions, are thus discussed in the context of the different simulations. In 20 conclusion, a 31-year long run of the fine-resolution Adriatic Sea climate model is found to be able to outperform most aspects of the reanalysis product, the short-term hindcast and the data assimilated simulation, in reproducing the dense water dynamics in the Adriatic Sea.


Introduction
The focus of this study is the Adriatic Sea -an elongated semi-enclosed basin located in the northern Mediterranean Sea. The 25 main geomorphological features of the Adriatic Sea (Fig. 1a) include a shallow bathymetry of the northern Adriatic shelf, gradually increasing in depth towards the 280 m deep Jabuka Pit. The middle Adriatic is separated from the ~1200 m deep Southern Adriatic Pit (SAP) by the Palagruža Sill, whereas in the very south, the Otranto Strait connects the Adriatic with the northern Ionian Sea. The Adriatic region is also characterized by an extremely complex eastern coastline with many islands and large mountain ranges along the entire basin. 30 3 In the Adriatic, past numerical studies have shown the importance of several factors influencing the modelling of the boradriven dense water formation. At first, the atmospheric forcing used as a source of forcing for the ocean models were not capable to properly reproduce the extreme bora events driving this process (Bergamasco et al., 1999;Vested et al., 1998;Beg-Paklar et al., 2001;Zavatarelli et al., 2002;Oddo and Guarnieri, 2011). Indeed, the resolution of the atmospheric model has 75 been found to be one of the most important characteristics known to impact the bora wind speeds, due to an improved reproduction of the orography and the enhancement of jet flows in finer grids (Belušić et al., 2017). Also, the importance of the ocean model resolution has been demonstrated through many studies that used kilometre-scale limited-area models to simulate ocean processes driven by extreme conditions in the Adriatic Sea (e.g., Cavaleri et al., 2010Cavaleri et al., , 2018Ricchi et al., 2016;Carniel et al., 2016;Denamiel et al., 2020b). Further, the influence of the freshwater forcing in the ocean models was foun d 80 to be crucial in modelling the dense water formation. In particular, the river runoff climatology used in previous studies (Raicich, 1994) overestimated real river discharges along the eastern Adriatic coast (Zavatarelli and Pinardi, 2002;Chiggiato and Oddo, 2008) and has been replaced by a new climatology which was based on up-to-date observations (Janeković et al., 2014). That significantly improved the reproduction of the dense water formation, in particular at its secondary source location in the Kvarner Bay (Vilibić et al, 2016Mihanović et al., 2018). Other factors such as the choice of open boundary 85 conditions, the parametrization of vertical mixing and diffusion, etc., were also found to be important (Benetazzo et al., 2014;Janeković et al., 2014).
The proper representation of the bora-driven dense water dynamics in the Adriatic Sea is still nowadays challenging, whether it is for research purpose as hindcast simulations and reanalysis products or for operational purpose as forecast simulations . This is why, recently, data assimilation was explored as an avenue to improve free model simulations, including the dense 90 water dynamics in the Adriatic Sea (Yaremchuk et al., 2016;Janeković et al., 2020). More particularly, the Four -Dimensional Variational scheme (4D-Var; Courtier et al., 1994;Janeković et al., 2013;Iermano et al., 2015;Sperrevik et al., 2017) was used during the 2014-15 period when a large number of in situ salinity, temperature and current observations were available (Janeković et al., 2020). Further, a 31-year (1987Further, a 31-year ( -2017 evaluation simulation of the Adriatic Sea and Coast (AdriSC; Denamiel et al., 2019) climate model using kilometre-scale atmosphere-ocean models over the Adriatic basin has also been 95 recently completed and evaluated (Pranić et al., 2021;Denamiel et al., 2021b). This kind of simulations, also referred as "Control Run" in the climate community, produce several-decade-long results forced by reanalysis products (without data assimilation) and are mainly used for evaluation purpose in climate studies. As a free run (i.e., dynamically consistent over decades contrary to reanalysis products which depend on the availability of the observations; Thorne and Vose, 2010), the AdriSC evaluation simulation has already provided invaluable information about the, till now unknown, kilometre-scale 100 present trends and variability of the Adriatic Sea (Tojčić et al., 2023). The aim of this study is thus to compare the currently available modelling strategies used to represent the bora -driven dense water dynamics in the Adriatic Sea. The approaches considered in the study are: (1) the newest high-resolution physical reanalysis product for the Mediterranean Sea (Escudier et al., 2020(Escudier et al., , 2021, hereafter referred as MEDSEA, which is generated by a numerical system composed of the Nucleus for European Modelling of the Ocean (NEMO) V3.6 model (Madec et al., 2017) and a three-dimensional variational (3D-Var) data assimilation scheme OceanVAR (Dobricic and Pinardi, 2008), and 110 forced by the ERA5 reanalysis (Hersbach et al., 2020), (2) the year-long simulations of an atmosphere-ocean Adriatic Sea forecast model (Janeković et al., 2020) composed of ROMS (Regional Ocean Modelling System; Shchepetkin and McWilliams, 2009) and ALADIN/HR (Aire Limitée Adaptation dynamique Développement InterNational; Tudor et al., 2013Tudor et al., , 2015 models used in both hindcast mode (hereafter referred as ROMS-hind) and with a 4D-Var data assimilation procedure (hereafter referred as ROMS-full), and (3) the recently evaluated 31-year simulation of the fine-resolution atmosphere-ocean 115 AdriSC climate model which is based on a modified version of the Coupled Ocean-Atmosphere-Wave-Sediment-Transport (COAWST V3.3) modelling system (Warner et al., 2010). AdriSC model is composed of the ROMS model (hereafter referred as AdriSC-ROMS), and the Weather Research and Forecasting (WRF v3.9.1.1; Skamarock et al., 2005) model (hereafter referred as AdriSC-WRF).
In the following section, the numerical models as well as the methods used to perform this study are described. The results o f 120 the analyses are presented in Sect. 3 and discussed in detail in Sect. 4. Finally, the main conclusions of the study are summarized in Sect. 5.

Northern Adriatic Experiment (NAdEx)
The time period investigated in this study includes the Northern Adriatic Experiment (NAdEx) campaign which took place 125 between late autumn 2014 and summer 2015. The aim of the NAdEx campaign was to study the dense water generation and transport within and off the Kvarner Bay. It consisted in collecting temperature, salinity and current data using several instruments and observing platforms. To measure the currents, acoustic Doppler current profilers (ADCPs) were deployed at 9 locations between late November 2014 and mid-August 2015, while conductivity-temperature-depth (CTD) probes measured salinity and temperature at 5 of the ADCP locations. Additionally, vertical profiles of temperature and salinity wer e 130 acquired at 19 CTD stations during two cruises between 3 and 6 December 2014 and between 26 and 29 May 2015. An ocean glider was operated off the Kvarner Bay in a campaign lasting only for 3 daysfrom 24 to 27 February 2015, while an Arvor-C profiling float was deployed on 19 February 2015 in the northern part of Kvarner Bay and was recovered on 15 March 2015 off the Istria coast. The full description of the NAdEx campaign is provided in Vilibić et al. (2018). During the campaign, three severe bora episodes with gusts above 50 m s −1 in the Velebit channel occurred: between 28 December 2014 and 1 135 January 2015, between 4 and 7 February 2015 and between the 3 and 6 March 2015. The NAdEx campaign was thus able to partially observe the dense water generation within the Kvarner Bay during the 2014-15 period.
Due to the unique dataset collected during the NAdEx campaignwhich has already been used either in data assimilation experiments (Janeković et al., 2020) or in evaluation studies Pranić et al., 2021) the 2015 bora events present a unique opportunity to compare the capacity of different models (e.g., reanalysis, hindcast, assimilated simulations, 140 climate simulations) to reproduce the dense water dynamics in the Adriatic basin.

Methods
The main features of the numerical models/products used in this study are presented in Table 1. Additional information and a more detailed description of the models is provided as Supplementary Material (section S2).
In order to compare different simulations, model results with horizontal grid resolution coarser than 1 km are interpolated t o 145 the AdriSC-ROMS 1 km grid using the regridding routines based on the Earth System Modeling Framework (ESMF) software (http://earthsystemmodeling.org/, last access: 2 March 2023). More specifically, the results of the ocean models (MEDSEA, ROMS-hind and ROMS-full) and atmospheric models (ERA5, ALADIN/HR-hind, ALADIN/HR-full and AdriSC-WRF) are all regridded to a horizontal resolution of 1 km. However, model outputs were not interpolated in the vertical.
Hereafter, the bottom potential density anomalies (PDAs) are calculated using the function available within the NCAR 150 Command Language (NCL) library (Levitus et al., 1994a(Levitus et al., , 1994bDukowicz, 2000;https://www.ncl.ucar.edu/, last access: 14 November 2022) and the upward turbulent heat fluxes are computed as the sum of the latent and sensible heat fluxes based on the formulas provided in Denamiel et al. (2020aDenamiel et al. ( , 2021a. To be noted, heat fluxes from the ALADIN/HR-full results are also modified by the 4D-Var data assimilation process. Further, in this study "dense waters" are defined as having PDAs equal or larger than 29.2 kg•m -3 based on previous research dealing with dense waters in the Adriatic (e.g., Janeković et al., 2014, 155 Vilibić et al., 2016. A comparative evaluation of the simulations for the 2014-15 period is carried out against in-situ temperature and salinity observations extracted from the NAdEx campaign , the Palagruža Sill long-term monitoring transect (Institute of Oceanography and Fisheries, Croatia) as well as the database published by Vilibić (2021) anddescribed in Prani ć et al. (2021). The assessment is presented in the form of a Taylor diagram (Taylor, 2001) using multiple statistical parameters 160 (normalized standard deviations and correlations) and probability density functions of the biases (PDFs). Normalized standard deviations are calculated as the ratio of the standard deviation of model results and standard deviation of observations, while the of the simulations and the available observations (i.e., they are daily instantaneous bias errors). That is, the model result s are extracted at the location (i.e., near neighbour grid point), depth (i.e., linear interpolation from model depths to observati on depth) and timing (i.e., approximated to daily average) of the observations. The biases are then obtained as the differenc e between model results and observations at each point in time, depth and space. The probability density functions are derived with a kernel-smoothing method (Bowman and Azzalini, 1997) which calculates the probability density estimate based on a 175 normal kernel function, and is evaluated at 100 equally-spaced points. Also, for each model, the median and the median absolute deviation (MAD) of the biases are calculated.
Additionally, to determine the minimum horizontal resolution necessary to resolve the processes occurring in the Adriatic Sea, the baroclinic Rossby radii are calculated for the whole AdriSC-ROMS 1 km domain. The results are presented as spatial distributions of median and MAD of the Rossby radius as well as in the form of a time series of the Rossby radius for the four 180 subdomains (northern Adriatic, Kvarner Bay, Jabuka Pit and deep Adriatic). In this study, the potential density method described in Chelton et al. (1998) is used to calculate the Brunt-Väisälä frequency and, hence, the first mode of the Rossby radii. This method can result in an underestimation of the Rossby radii if the vertical spacing is not fine enough. As the AdriSC-ROMS model is providing results on 35 vertical sigma layers for the entire Adriatic Sea, this underestimation can thus only occur within the SAP area. 185 Further, to quantify the general conditions in the ocean and atmosphere throughout the whole 2014-15 period, an analysis of the extremes is performed. For each of the four simulations, the results are presented as spatial distributions of extremes accompanied with spatial distributions of their associated timing (in days). This includes the spatial distributions of the maximum wind stresses and the maximum upward turbulent heat fluxes in the atmosphere (at the sea surface), and of the maximum PDAs, minimum temperatures and maximum salinities in the ocean. 190 In order to better capture the dense water dynamics, two different temporal analyses of the results are also performed. First , time series of daily surface wind stresses and upward turbulent heat fluxes in the atmosphere, and of daily bottom PDAs, temperatures and salinities in the ocean are presented as the spatial average over different subdomains selected in areas where dense waters are known to be either generated or accumulated. These subdomains are the northern Adriatic and the Kvarner Bay for both atmosphere and ocean as well as the Jabuka Pit and the deep Adriatic for the ocean only (Fig. 1a). In addition, 195 the daily bottom PDA time series are presented without the mean and the seasonal signal (yearly and half-yearly) which are removed from the series at each point. More specifically, after subtracting the mean and detrending the time series, the seasonal signal is calculated using the least-squares method of a sine function and subtracted from the series. The final time series without the seasonality are obtained by adding the trend. Lastly, the time evolution of the spatial distributions of the bottom PDAs is presented at selected dates -1 March, 1 April, 1 May and 1 June 2015and for the whole 2014-15 period as a movie. 200 Deleted: ¶ Due to the unique dataset collected during the NAdEx campaign which has already been used either in data assimilation experimen (Janeković et al., 2020) or in evaluation studies (Vilibić et al., 201 Pranić et al., 2021) the 2015 bora events present a unique 205 opportunity to compare the capacity of different models (e.g., reanalysis, hindcast, assimilated simulations, climate simulations reproduce the dense water dynamics in the Adriatic basin. ¶

Methods ¶
The main features of the numerical models/products used in this 210 study are presented in Table 1. Additional information and a mor detailed description of the models is provided as Supplementary Material (section S2). ¶ In order to compare different simulations, model results with horizontal grid resolution coarser than 1 km are interpolated to th 215 AdriSC-ROMS 1 km grid using the regridding routines based on Earth System Modeling Framework (ESMF) software (http://earthsystemmodeling.org/, last access: 2 March 2023). Mo specifically, the results of the ocean models (MEDSEA, ROMS-h and ROMS-full) and atmospheric models (ERA5, ALADIN/HR-h 220 ALADIN/HR-full and AdriSC-WRF) are all regridded to a horizo resolution of 1 km. Model outputs were not interpolated in the vertical. ¶ Hereafter, the bottom potential density anomalies (PDAs) are calculated using the function available within the NCAR Comma 225 Language (NCL) library (Levitus et al., 1994a(Levitus et al., , 1994bDukowicz 2000; https://www.ncl.ucar.edu/, last access: 14 November 2022) the downward upward turbulent heat fluxes are computed as the s of the latent and sensible heat fluxes based on the formulas provid in Denamiel et al. (2020aDenamiel et al. ( , 2021a. To be noted, heat fluxes from t 230 ALADIN/HR-full results are also modified by the 4D-Var data assimilation process. Further, in this study "dense waters" are def for as having PDAs equal or larger than 29.2 kg•m -3 based on previous research dealing with dense waters in the Adriatic (e.g., Janeković et al., 2014, Vilibić et al., 2016. ¶ 235 Moved up [1]: A comparative evaluation of the simulations fo the 2014-15 period is carried out against in-situ temperature and salinity observations extracted from the NAdEx campaign (Vilibi al., 2018), the Palagruža Sill long-term monitoring transect (Instit of Oceanography and Fisheries, Croatia) as well as the database 240 published by Vilibić (2021) and described in Pranić et al. (2021). assessment is presented in the form of a Taylor diagram (Taylor, 2001) using multiple statistical parameters (normalized standard deviations and correlations) and probability density functions of t biases (PDFs). Normalized standard 245 Deleted: deviations are calculated as the square root of the ratio the variance of model results and variance of observations, while correlations are the Pearson's linear correlation coefficients. The biases are calculated as differences between the daily results of th simulations and the available observations (i.e., they are daily 250 instantaneous bias errors). ConsequentlyThat is, the model results extracted at the location (i.e., near neighbour grid point), depth (i. linear interpolation from model depths to observation depth) and timing (i.e., approximated to daily average) of the observations. T biases are then obtained as the difference The final analysis quantifies the total daily volume transport of the outflowing dense waters across four transects (T1-T4; Fig.   1a) for all depths. The outward transport is calculated as a double integral of velocities normal to the transect over the area of 260 the vertical plane of the transect.

Comparative evaluation during the 2014-15 period
A brief comparative evaluation of the four simulations is performed in order to quantify the skills of the ocean models against 18987 CTD measurements (Fig. 2c). The number of observations depending on the depth is: (1)  For the temperature biases ( Fig. 2b), ROMS-hind distribution has a median of -0.37 associated with a large peak and a MAD of ± 0.33 °C. ROMS-full distribution has a smaller peak and a median of -0.29 ± 0.31 °C. MEDSEA distribution has a median of 0.00 ± 0.84 °C with a heavier tail of positive biases up to 4.5 °C. AdriSC-ROMS distribution has the lowest peak and a median of -0.04 ± 0.61 °C. Therefore, the ROMS simulations systematically underestimate the sea temperature but the assimilation reduces the biases. Median temperature bias is smaller in AdriSC-ROMS and MEDSEA, but they have the largest 280 MAD due to an overestimation of temperatures in MEDSEA (up to 4.5 °C) and both an over-and underestimation of the temperatures between -2 and +2 °C in AdriSC-ROMS.
For the salinity biases ( Fig. 2d), ROMS-hind distribution has the lowest peak and a median of -0.16 ± 0.12. ROMS-full distribution has a larger peak and a median of -0.09 ± 0.09. MEDSEA distribution has a median of 0.00 ± 0.36, a tail of negative biases down to -2.0 and a heavy tail of positive biases with a secondary peak at approximately 1.0. AdriSC-ROMS 285 distribution has a slightly larger peak than MEDSEA and a median of 0.02 ± 0.16 with very low probabilities for negative biases below -0.2 but a heavy tail up to around 1.0 and a secondary peak around 0.4. Hence, the ROMS-full and hind 9 simulations both underestimate the observed salinities but the assimilation reduces the biases. The AdriSC-ROMS model tends to overestimate the salinity while the MEDSEA results display the largest over-and underestimations of salinities. 315 Lastly, the comparison of the performance of models with different resolutions may be affected by the double-penalty effect (Crocker et al., 2020) meaning that in pointwise comparison with observations the finer resolution models tend to be penalised more than the models with coarser resolution and therefore they can verify worse. When a model has sufficient resolution to reproduce a small-scale feature but it simulates it incorrectly, it is penalised twice: once for not simulating the feature where it should have been and once for simulating it where it hasn't been observed. Contrarily, if a model resolution is not sufficient 320 to reproduce a feature, it will be penalised only once for not reproducing the feature. This might partly explain why the AdriSC-ROMS model presents a larger bias variability in both temperature and salinity than the ROMS-hind and ROMS-full models.

Analysis of the extremes
To analyse how the different models capture the extremes during the 2014-15 period, the spatial distributions of daily maximum wind stresses, daily maximum upward turbulent heat fluxes and their associated timing are presented in Figures 3  325 and 4, while the spatial distributions of daily maximum bottom PDAs and their associated timing are presented in Figure 5.
Additionally, the spatial distributions of minimum temperatures and maximum salinities are provided and described in Supplementary Material ( Fig. S1 and S2).

Wind stresses and upward turbulent heat fluxes
It should be noted that ERA5, which is forcing the MEDSEA reanalysis, produces very small wind stresses over the whole 330 basin ( Fig. 3a), barely reaching 0.4 N m -2 in the northern Adriatic, while ALADIN/HR-hind and ALADIN/HR-full wind stress results (Fig. 3c,3e) are extremely similar despite the variational scheme of the assimilation changing the wind stresses (i.e., the differences between the ALADIN/HR-hind and full wind stresses are at least an order of magnitude smaller than their differences with the other atmospheric models). Further, AdriSC-WRF, which is the only kilometre-scale atmospheric model used in this comparison, overall generates the largest extremes for both wind stresses (> 1.5 N m -2 ; Fig. 3g) and upward 335 turbulent heat fluxes (> 1100 W m -2 ; Fig. 4g). However, for the upward turbulent fluxes, ERA5 produces maximum heat losses comparable to AdriSC-WRF (Fig. 4a, 4g) while ALADIN/HR-full maximum heat losses are at least twice smaller than in ALADIN/HR-hind (Fig. 4c, 4e). In fact, ALADIN/HR-full has the smallest maximum heat losses of all simulations and shows a patchy spatial distribution with the smallest values over the middle of northern Adriatic, barely reaching 750 W m -2 in February-March. Consequently, both MEDSEA and ALADIN/HR-full are strongly influenced by the assimilation (e.g., sea 340 surface temperature coming from remote sensing products or variational changes of the heat flux forcing, respectively).
Another important point is that the turbulent heat fluxes are strongly influenced by the sea surface temperature and the relative humidity, which are in turn influenced by the solar radiation. The maximum heat losses are thus more likely to be found in December 2014-January 2015due to a difference in air-sea temperatures of about 3-4 °C having a larger contribution in the Deleted: re upward turbulent heat flux calculation than the intensity of the wind stresses (Fairall et al., 1996) than in early February/March 2015 when the temperature differences are smaller.
In the northern Adriatic, the Trieste Jet is seen by ALADIN/HR-hind/full and AdriSC-WRF models with wind stress maxima 360 reaching 0.8 N m -2 and 1.3 N m -2 , respectively. It is also important to highlight that the Trieste Jet produced by ALADIN/HRhind/full is further extended offshore than in the AdriSC-WRF simulation. For the upward turbulent fluxes, ERA5 and ALADIN/HR-hind/full give small intensity along the Trieste Jet (less than 600 W m -2 ) between January and March 2015, while AdriSC-WRF reaches 850 W m -2 in December 2014.
The largest values of maximum wind stresses are found in the Kvarner Bay and along the Senj Jet for all simulations including 365 ERA5. They reach up to 1.3 N m -2 for ALADIN/HR-hind/full and more than 1.5 N m -2 for AdriSC-WRF over a far wider area than the other models. In this region, for the upward turbulent heat fluxes, the maximum values are reached by ERA5 (i.e., 900 W m -2 despite not reproducing the bora jets) and AdriSC-WRF (larger than 1100 W m -2 ), while ALADIN/HR-hind and full reach 850 W m -2 and barely 750 W m -2 , respectively.
In the middle Adriatic, high wind stresses up to 1.2 N m -2 for ALADIN/HR-hind/full and AdriSC-WRF models are produced 370 along the Karlobag/Sukošan Jets. However, AdriSC-WRF extends the Karlobag Jet to the middle of the Adriatic with values (up to 1.5 N m -2 ) several times larger than achieved with the ALADIN/HR-hind/full simulations. It also produces some strong wind stresses up to 1.3 N m -2 along the Dalmatian coast where other bora jets are known to be located. In terms of upward turbulent fluxes, the maximum values are in average 900 W m -2 for ERA5, 800 W m -2 for ALADIN/HR-hind, larger than 400 W m -2 for ALADIN/HR-full, and larger than -1100 W m -2 along the eastern Adriatic coast for AdriSC-WRF. 375 In the southern Adriatic, maximum wind stresses in ALADIN/HR-hind/full reach up to 0.7 N m -2 but are smaller along the coastline. In the AdriSC-WRF simulation, the wind stresses remain relatively small in the southern Adriatic (less than 0.5 N m -2 ), aside from a small patch of larger values off the southern Montenegrin coast. For the upward turbulent fluxes, the results obtained with ERA5 and AdriSC-WRF are quite similar with strong intensities along the eastern coast (in average 900 W m -2 and 1000 W m -2 , respectively) and values less than 700 W m -2 offshore. 380 Overall, for all models, maxima of wind stresses are associated with bora events, while upward turbulent heat fluxes seem to be influenced by the seasonal variations of the sea surface temperature (SST) more than the wind stresses. In other words, the largest input to the upward turbulent heat fluxes is coming from the bora wind, yet a small fractionwhich is found to influence maxima of the heat fluxesis coming from SST. That is the reason why maxima of heat fluxes occur mostly during bora episodes in late December/Early February (Fig. 4), whereas the maxima of wind stresses occur mostly during bora episodes in 385 early February/early March (Fig. 3). Additionally, the AdriSC-WRF model generates the strongest dynamics with, in average, the strongest wind stresses and the maximum heat cooling, while ERA5 has the weakest wind stresses and ALADIN/HR-full the smallest heat loss. 415

Potential density anomalies
In the northern Adriatic, all simulations produce the highest maximum PDA values during late winter (February-March 2015;Fig. 5b, d, f, h). They reach up to 29.4 kg m -3 on the shelf for MEDSEA, up to 29.6 kg m -3 along the coast but below 29.3 kg m -3 on the shelf for ROMS-hind, up to 29.8 kg m -3 along the coast and in average 29.5 kg m -3 on the shelf for ROMS-full, and finally, above 29.8 kg m -3 along the coast and in average 29.7 kg m -3 for AdriSC-ROMS (Fig. 5a, c, e, g). 420 In the Kvarner Bay, both MEDSEA and ROMS-hind have extremely low maximum PDAs (below 29.0 kg m -3 ), indicating no dense water formation in this area. In contrast, both ROMS-full and AdriSC-ROMS give high maximum PDAs (up to 29.6 kg m -3 and 29.7 kg m -3 respectively). However, ROMS-full presents patch-like PDA distributions with maxima occurring partly during winter and partly during spring, while AdriSC-ROMS has more homogeneous values over the whole Kvarner Bay with maxima occurring mostly in the winter but also in September in a few very small areas. Further, off the Kvarner Bay, ROMS-425 full produces a large patch of extremely dense waters (> 29.8 kg m -3 ) which does not seem to be smooth and continuous with the previous data assimilation cycle spatial PDA distributions over the rest of the Adriatic domain. In this case, data assimilation is correcting for the initial state of the ocean model at the start of the assimilation cycle, as the most cost-effective mechanism for correcting suboptimal atmospheric (hydrostatic and coarser) forcing and ocean model vertical and horizontal resolution constraints. This patch occurred in February and is located just southwest from the glider data assimilated in the 430 model, which is the strongest contributor to the data assimilation cost function at that period (Janeković et al., 2020).
In the middle Adriatic, ROMS-hind shows relatively low maximum PDAs (below 29.1 kg m -3 ) but the other models present some interesting spatial variations. In the Jabuka Pit, which is a known dense water reservoir, maximum PDAs reach up to In the northern Adriatic, all models present three prominent peaks of wind stresses (Fig. 6a), capturing the three severe bor a events that occur during the NAdEx campaign: 28 December 2014 -1 January 2015, 3-7 February 2015 and 3-6 March 2015.
These dominant wind stress events are also associated with peaks of upward turbulent heat fluxes in all models (Fig. 6c).
However, the intensities of the ERA5 wind stress peaks (0.15, 0.3 and 0.2 N m -2 ) are half those in ALADIN/HR-full/hind and AdriSC-WRF which are all similar (peaks at 0.3, 0.6 and 0.5 N m -2 ). Further, the intensity of the upward turbulent heat flux 470 peaks is often smaller and more spread or shifted over time in ALADIN/HR-full (peaks at 300, 450 and 300 W m -2 ) than in the other models, due to the variational scheme used in the assimilation. It should be noted that the strongest peaks in upward turbulent heat fluxes are always reached by ERA5 and/or AdriSC-WRF (peaks at 600, 400 and 350 W m -2 ), while ALADIN/HR-hind produced slightly smaller intensities in general (peaks at 500, 350 and 300 W m -2 ). Concerning the associated bottom PDA time variations (Fig. 7a), it should be first noted that the AdriSC-ROMS PDAs are systematically 475 higher than in the other models by 0.2-0.8 kg m -3 due to higher salinity (differences of about 0.3-0.6; Fig. S4a). Second, for all simulations, the maximum values are obtained between February and March 2015, when the dense water generation is found to occur . Further, in February 2015, a large increase of bottom PDAsprobably driven by the assimilation of the Arvor-C, towed CTD and glider data which influenced two 4-day cyclesis seen in ROMS-full, which reaches values Deleted: lower nearly as high as in AdriSC-ROMS. The PDAs without seasonality show that the peaks of density due to the bora-driven dense water formation are reproduced in all models (Fig. 8a). The highest increases in density during these peaks are always reached 495 by ROMS-full (0.4, 0.35 and 0.3 kg m -3 ) and the lowest by MEDSEA (below 0.1 kg m -3 for the 3 peaks). However, the MEDSEA and AdriSC-ROMS densities already increased before the first bora event by 0.2 kg m -3 which means that in fact the highest peak is reached by AdriSC-ROMS after the first bora event and that MEDSEA densities are close to AdriSC-ROMS values. The PDAs without seasonality also clearly show, for all models, a decrease in density during spring and summer when the denser waters are transported from the northern Adriatic towards the south. 500 In the Kvarner Bay, the three bora peaks of wind stresses (Fig. 6c) and the associated upward turbulent heat fluxes (Fig. 6d) are also seen by all models. However, ERA5 computed wind stresses are always extremely low (below 0.2 N m -2 ) while AdriSC-WRF produces stronger wind stresses (peaks at 0.5, 0.6 and 1.5 N m -2 ) than ALADIN/HR-full/hind (peaks at 0.25, 0.4 and 0.9 N m -2 ). The intensity of the upward turbulent heat flux peaks is again always less and more spread or shifted over time in ALADIN/HR-full (peaks at 400, 400 and 300 W m -2 ) than in the other models (peaks as large as 800, 500 and 600 N 505 m -2 ). Also, AdriSC-WRF models produce eight wind stress peaks above 0.25 N m -2 between December 2014 and April 2015, while ALADIN/HR-hind/full only surpasses this threshold for the three main bora events. That is, the non-hydrostatic kilometre-scale AdriSC-WRF model (at 3 km resolution) is capable to reproduce much higher wind stresses than the hydrostatic ALADIN/HR model (at 8 km resolution dynamically downscaled to 2 km for the winds only) due to the impact of the highly non-linear orographic processes on the dynamics of the bora-driven flows (e.g., Grubišić, 2004;Kuzmić et al., 510 2015). Next, the upward turbulent heat fluxes are less intense in ERA5 and ALADIN/HR-hind than in AdriSC-WRF, indicating that the cooling rates are smaller which thus should lead to less generation of dense waters. In terms of bottom PDA analysis ( Fig. 7b), similarly to the northern Adriatic subdomain, the AdriSC-ROMS model produces the highest values, while MEDSEA and ROMS-hind generally have the lowest values with differences up to 0.6 kg m -3 in February-March 2015. This difference is again mostly driven by salinity, which is the lowest in MEDSEA and again the highest in AdriSC-ROMS (Fig. S4). However, 515 salinity is much higher in ROMS-full than in ROMS-hind starting in December 2014, when near bottom salinity measurements were available continuously in the Kvarner Bay through the NAdEx campaign. Convincingly, these measurements moved the ROMS-full run from ROMS-hind towards the higher measured salinities and closer to the AdriSC-ROMS results. As for the northern Adriatic subdomain, the PDAs without seasonality show three main peaks linked to bora-driven dense water formation in all the models (Fig. 8b). However, the timing of the ROMS-full peaks as well as their intensity is generally 520 different than for the other models (which all behave quite similarly), particularly after the second and third bora events. This shows the impact of the assimilation of the NAdEx campaign observations within the ROMS-full model.
In the Jabuka Pit ( Fig. 7c and 8c), bottom PDAs (with and without seasonality) from the two free model runs (AdriSC-ROMS and ROMS-hind) increase from February 2015, when newly generated denser waters from the northern Adriatic start to fill the pit, and peak in late April 2015. However, AdriSC-ROMS PDAs are higher than ROMS-hind both in mean values (more 525 than 29.1 kg m -3 vs. less than 29.0 kg m -3 ) but particularly in increase rates (0.2 kg m -3 in 2 months vs. less than 0.1 kg m -3 in 2 months) during the known arrival time of dense waters in the Jabuka Pit This shows that no dense water arrival in the Jabuka Pit is seen by MEDSEA during spring 2015.
In the deep Adriatic (Fig. 7d), bottom PDA values are similar in all models with slightly higher values in ROMS-hind/full and smaller values in MEDSEA and AdriSC-ROMS. Further, temporal changes in PDAs are higher in ROMS-full and MEDSEA as they assimilate deep observations (e.g., by Argo profilers up to 700-800 m) which were available during the whole 2014-15 period (Kokkini et al., 2020), as can be clearly seen in the PDAs without seasonality (Fig. 8d). 550 Overall, the analysis of the time series spatially averaged over the subdomains where dense waters are either generated (i.e.,

Time evolution of the bottom PDA spatial distributions
To better visualize the evolution in time and space of the dense waters, the spatial distributions of the daily bottom PDAs a re 560 analysed both at specific datesi.e., 1 March (Fig. 9), 1 April (Fig. 10), 1 May (Fig. 11) and 1 June 2015 (Fig. 12)  Adriatic in ROMS-full and AdriSC-ROMS. However, due to the availability of assimilated measurements, ROMS-full first generates dense waters off the Kvarner Bay and then within. In contrast, AdriSC-ROMS clearly transports the dense waters generated within the Kvarner Bay towards the west along the bora jets. On the 1 March 2015 (Fig. 9), dense waters are starting to be collected within the Jabuka Pit in both ROMS-full and AdriSC-ROMS, while no dense water has been transported that far south in MEDSEA and  Between the third bora event and 1 April 2015, for ROMS-full and AdriSC-ROMS, after an initial increase along the bora jets, dense waters (above 29.5 kg m -3 ) are transported along the western coast from the northern Adriatic and the Kvarner Bay towards the south, and partially collected in the Jabuka Pit. ROMS-hind also shows some dense water transport (with PDAs barely reaching 29.2 kg m -3 ) from the northern Adriatic towards the Jabuka Pit. However, in MEDSEA, the dense waters generated in the northern shelf (up to 29.45 kg m -3 ) seem to slowly dissipate without being transported. On 1 April (Fig. 10) Between the 1 April and the 1 May 2015, in ROMS-full and AdriSC-ROMS, continuous transport towards the south results in 600 a larger amount of dense waters being collected in the Jabuka Pit from where they start to cascade towards the SAP via the deepest parts of the Palagruža Sill (Rubino et al., 2012). It should be noted that the cascading occurs along a narrower and more western path in AdriSC-ROMS than in ROMS-full. On the 1 May 2015 (Fig. 11), no dense water is present in the Overall, AdriSC-ROMS generates a larger amount of dense waters than the other models because of the strongest atmospheric forcing, while MEDSEA and ROMS-hind do not properly reproduce the dense water dynamics in the Adriatic basin. However, ROMS-full collects a larger amount of dense waters in the Jabuka Pit than all the other models. It can be concluded that 620 AdriSC-ROMS is probably too dissipative during the transport of the dense waters from the northern Adriatic and the Kvarner Bay towards the south. Further, in ROMS-full, the patchy distribution of very dense waters during winter and spring can be explained by the assimilation of data in 4-day cycles for which CTD measurementscollected at some given sites and for some specific days -took a significant role in adjusting the Adriatic dynamical solutions (Janeković et al., 2020). This demonstrates the importance of the coverage and the long-term availability of the assimilated data. A better representation of 625 the dense water dynamics within the Adriatic basin in ROMS-hind can thus be envisioned (and is possible as demonstrated by the results of the AdriSC model) before performing the data assimilation which, for the moment, is incapable to fully compensate the cumulated weaknesses of the ALADIN/HR and ROMS-hind models.

Daily volume transports along selected transects
To quantify the dense water outflow across different sections of the northern and middle Adriatic, the volume transports of 630 dense water defined by the PDA threshold of 29.2 kg m -3 through four transects (T1 to T4) are presented in Figure 13. The transport is defined as positive towards northwest (transects T1, T3 and T4) or northeast (transect T2). In general, MEDSEA and ROMS-hind transports are the lowest for all transects, which is expected as their overall PDA values are found to be the lowest of all simulations. With the same argument, the AdriSC-ROMS transport is the highest for all transects, except for T4, where the ROMS-full transport prevails (Fig. 13d). 635 The transports produced by MEDSEA at T1 are mostly very low, peaking at -0.03 Sv in February (Fig. 13a) ROMS transports indicate that the ratio between dense water originating from the northern Adriatic and the Kvarner Bay is roughly 60:40, which is similar to the transport ratio derived for the massive dense water generation in winter 2012 (Janeković et al., 2014).
For T3 and T4, both MEDSEA and ROMS-hind transports are null throughout the whole period. ROMS-full transports at T3 vary around -0.05 Sv from the middle of February to the end of May (Fig. 13c), being the highest in the second half of March 650 and reaching -0.20 Sv. Furthermore, the results show some similarities in the behaviour of the ROMS-full and AdriSC-ROMS transports. Interestingly, the dense water transports at T3 are lagged for about two to three weeks (depending on the simulation) after the transports at T1 and T2, from which an estimate of bottom density current may be computed (approximately 0.10-0.17 m s -1 ).
Lastly, ROMS-full transports are extremely large at T4, much larger than in AdriSC-ROMS, peaking during March-April with 655 values reaching -0.90 Sv and -0.70 Sv, respectively. For the rest of the time, the transports are relatively small, questioning if these outbursts of dense water are driven by the assimilated data or by an outflow of dense waters with high densities that are produced by ROMS-full northwest from transect T4, in the Jabuka Pit.

Discussion
The multi-model analysis performed in this study has demonstrated that reproducing the dense water dynamics within the 660 Adriatic basin is extremely complex as the presented models produced different or even divergent results despite all being thoroughly evaluated in previous studies (Escudier et al., 2021;Janeković et al., 2104;Vilibić et al., 2018;Pranić et al., 2021; Deleted: down to  Denamiel et al., 2021bDenamiel et al., , 2022. However, it is important to keep in mind that the presented results belong to different model categories: MEDSEA is a reanalysis product covering the full Mediterranean Sea for the 1987-2019 period, ALADIN/HR-ROMS does not cover the full Adriatic Sea and is used, in this study, either in hindcast mode (hind) or fully assimilated mode with 4-day cycles (full) for the 2014-15 period, and finally, AdriSC is the evaluation run of a climate model covering the full Adriatic for the 1987-2017 period. This implies that numerical schemes (e.g., discretization, parametrization) and set-up (e.g., 675 physics, resolution, forcing) used in these models as well as the type of simulation performed (free run vs. assimilated run) strongly influence the quality of the presented results. As this study only compares state-of-the-art models (ERA5, WRF, ALADIN in the atmosphere and NEMO, ROMS in the ocean), the differences in numerical schemes will not be discussed hereafter because it is difficult to quantify how they impact the dense water dynamics as they vary from model to model.
Though, the different model set-ups will be analysed with the aim to better understand their impact on the bora-driven dense 680 water dynamics in the Adriatic basin.

Impact of the resolution and the physics on the bora dynamics
First, the ERA5 reanalysis at 25 km resolution has been demonstrated to be incapable to capture the bora dynamics (Denamiel et al., 2021a). Consequently, in this study, ERA5 wind stresses are two to three times smaller than the AdriSC-WRF and ALADIN/HR results. However, both in the northern Adriatic and in the Kvarner Bay, heat losses calculated from the ERA5-685 MEDSEA modelvia bulk formulae using sea surface temperature assimilating remote sensing productsare comparable to the ALADIN/HR-ROMS-hind model (Fig. 6). These heat losses are still underestimated compared to the AdriSC model, particularly within the Kvarner Bay and the Gulf of Trieste as well as along all the bora jets (Fig. 3).
Second, the hydrostatic ALADIN/HR model at 8 km resolutionwith the wind fields further dynamically downscaled to 2 kmhas already been demonstrated to reproduce the basic bora dynamics (Horvath et al., 2009). However, in the Kvarner 690 Bay region, our results show that the ALADIN/HR wind stresses are not as intense and not covering as wide an area as the non-hydrostatic AdriSC-WRF model. Indeed, the bora cross-flow variability in the Kvarner Bay might occur at a kilometre scale, in particular during deep bora events (Kuzmić et al., 2015), while bora pulsations have a strong sub-kilometre spatial component, posing a challenge for proper quantification in any kilometre-scale atmospheric model. Nevertheless, Denamiel et al. (2021a) have demonstrated that, during 22 bora events including two in 2015, the AdriSC-WRF 3 km model reproduced 695 very well the wind speed observations at Pula, Rijeka, Ogulin, Zavižan, Gospić and Knin stations (all located in the Kvarner Bay region) above 20 m/s despite over predicting them by up to 5 m/s below this threshold. Further, the ALADIN/HR-ROMShind heat losses are always smaller than those computed from ERA5-MEDSEA and AdriSC models. It is documented that hydrostatic atmospheric models are not capable to capture all the details of the bora jets (Klemp and Durran, 1987;Blockley and Lyons, 1994;Grisogono and Belušić, 2009). Consequently, the hydrostatic approximation used in ALADIN/HR constrains 700 its ability to reproduce the finer-scale details of the bora flow (Horvath et al., 2009). Therefore, heat losses in ALADIN/HR-ROMS (hind and full) mostly occur along the Senj Jet but are still weaker than in AdriSC (Fig. 3). Further, quite surprisingly, the 4D-Var data assimilation scheme used in the ROMS-full assimilation is reducing the intensity of the turbulent heat fluxes and thus creating a dynamical imbalance between the wind stresses (which are similar in comparison to the differences between the different atmospheric models) and the heat losses forcing the ocean model. 705 Finally, the evaluations of AdriSC-WRF model performed both for the climate run over a 31-year period (Denamiel et al., 2021b) and during extreme bora events (Denamiel et al., 2020a(Denamiel et al., , 2020b(Denamiel et al., , 2021a have demonstrated that a 3 km resolution is appropriate to represent the atmospheric dynamics within the Adriatic basin. Further, the results of the AdriSC-WRF model at 3 km resolution (particularly the intensity of the winds) have been shown to converge towards the results obtained with the higher-resolution AdriSC-WRF-1.5 km model during bora events (Denamiel et al., 2021a). However, only sub-kilometre-scale 710 atmospheric models can properly capture the highly non-linear dynamics of the bora flows (Kuzmić et al., 2015) and thus using a 3 km non-hydrostatic model is still a compromise between accuracy and efficiency. This compromise is particularly important when running multi-year/climate simulations having a tremendous computational cost. This is also highlighted by Vodopivec et al. (2022), who conducted a sensitivity study over a 16-year period using different runoff configurations and different sources of atmospheric forcing and concluded that the atmospheric forcing has a substantial impact on the hydrology 715 and circulation of the Adriatic Sea.

Impact of the resolution and the bathymetry on the dense water dynamics
In the ocean models, the resolution is mostly going to impact the representation of the many islands located along the easter n Adriatic coast but more importantly, of the reservoirs collecting the dense waters within the Adriatic basin (i.e., Kvarner Bay, Jabuka Pit and SAP). To better understand the necessary horizontal resolution needed to reproduce the Adriatic Sea dynamics, 720 the spatial distributions of the median and MAD of the Rossby radii calculated from the AdriSC-ROMS results are presented for the entire model domain on Fig. 14a and 14b, respectively. In general, the median Rossby radius is decreasing from open seas towards shallower coastal areas. The largest values are found to be around 10.0 ± 2.0 km in the open northern Ionian Sea.
Median Rossby radii are slightly smaller in the SAP with values around 7.5 ± 1.3 km while sharply decreasing on the edges of the pit to around 5.0 ± 1.2 km. In the Jabuka Pit, the radii reach around 4.0 ± 1.2 km whereas in the rest of the middle Adriatic 725 around 2.5 ± 1.2 km. The deeper part of the Kvarner Bay presents high variability and Rossby radii around 2 ± 1.5 km. The smallest median Rossby radii as well as the smallest MAD are calculated for the northern Adriatic around 1.0 ± 0.4 km. Further, the time series of the Rossby radius are presented for the northern Adriatic and Kvarner Bay (Fig. 3c) as well as for the Jabuka Pit and deep Adriatic (Fig. 14d) Kurkin et al. (2020) in their study dedicated to analysing the first 745 Rossby Radii in European semi-enclosed basins. Overall, the baroclinic Rossby radii present large variability in the Adriatic Sea and the results suggest that even sub-kilometre scale ocean models are needed to simulate full range of processes in the Adriatic Sea, particularly the dense water dynamics. However, for climate simulations horizontal resolution finer than 1 km is not feasible yet.
Further, different Digital Terrain Models (DTMs) have been used to generate the bathymetries of the presented models. In 750 order to evaluate the joint impact of resolution and bathymetry, MEDSEA and ROMS-hind/full bathymetries are compared to the AdriSC-ROMS model at 1 km resolution (Fig. 1b, c). The MEDSEA model is clearly shallower than AdriSC-ROMS within the Kvarner Bay and the Jabuka Pit (by 60-80 m) but also in the middle of the SAP (by more than 100 m). Consequently, the capacity of the MEDSEA model to naturally collect the dense waters within the known Adriatic reservoirs is decreased compared to the AdriSC-ROMS model and thus relies heavily on the assimilation of the available data. In the ROMS-hind/full 755 model, the bathymetry is also generally shallower than in AdriSC-ROMS within the Kvarner Bay and along the canyon system between the Jabuka Pit and the SAP (by 20-40 m). This is particularly important as it might explain the differences in paths seen between ROMS-full and AdriSC-ROMS when the dense waters are transported from the Jabuka Pit towards the SAP.
However, concerning the Jabuka Pit and the SAP, the alternated positive and negative differences in bathymetry between ROMS-full/hind and AdriSC-ROMS clearly show some shifts in locations. Whether and how these shifts in location impact 760 the dense water dynamics is not clear with the results presented in this study.
Finally, it is important to highlight that the AdriSC-ROMS model uses 35 vertical sigma layers while the ROMS-full/hind model only has 20 of them. As the bora-driven dense water dynamics requires to properly resolve both the surface (for the sea temperature cooling) and the bottom (for the dense water transport) layers, the finer vertical resolution used in AdriSC-ROMS may play a major a role in the overall performance of the model. 765

Impact of the salinity forcing on the dense water generation
Dense water generation is highly sensitive to the background salinity content provided either through the open boundaries or the direct river outflows imposed on the ocean models.
First, in ROMS-hind, Janeković et al. (2014) quantified an underestimation of salinity by 0.2-0.5 for a simulation of the massive dense water formation in 2012. After updating the old river climatologies that largely overestimated the discharges, Vilibić et 770 al. (2016) confirmed that even the simulations using the most realistic river representation underestimate the background salinity content within the Adriatic basin. As the AREG model (forcing ROMS-full) is set-up with the old river climatologies and has a low salinity content over the entire Adriatic basin, far too much fresh water is inputted through the ROMS-hind open lateral boundary located in the southern Adriatic. Consequently, the ROMS-hind results presented in this study for the 2014-15 period have low basin-wide salinities and therefore generate dense waters with lower bottom PDA values. 775 Next, the AdriSC-ROMS model has been thoroughly evaluated over a 31-year period in Pranić et al. (2021). First, in the 780 northern Adriatic, despite a lack of accuracy for salinities under 36, due to the Po River misrepresentation, the AdriSC-ROMS model has been shown to perform well in reproducing dense water masses. Second, in the Kvarner Bay, AdriSC-ROMS salinities have been demonstrated to be too high, which could lead to a general overestimation of the dense water bottom PDAs in this region. And finally, in the SAP, the evaluation revealed that the salinities and the densest waters are captured rela tively well by the AdriSC-ROMS model. 785 Finally, salinities in MEDSEA are closer to the AdriSC-ROMS results in the southern Adriatic (i.e., Jabuka Pit and deep Adriatic subdomains) and to the ROMS-hind results in the northern Adriatic (i.e., northern Adriatic and Kvarner Bay subdomains) during the entire 2014-15 period (Fig. S4). It can thus be safely assumed that the old river climatologies used in MEDSEA are resulting in low salinities over the northern part of the Adriatic basin and hence lower bottom PDAs during the bora-driven dense water generation events. 790

Impact of the assimilation on the ocean dynamics
First, in ROMS-full, the 4DVar data assimilation is applied in 4-day cycles which means that the ocean dynamical properties are not continuously smooth in time between the cycles as the ROMS-full model adjusts the initial state at the beginning of each cycle. Consequently, despite the large improvement of the ocean fields used to minimize the cost function of the assimilation, the dense water generation and transport as a continuous process in time is not properly reproduced in ROMS-795 full. For example, as the salinity is generally underestimated in ROMS-hind, the data assimilation performed in ROMS-full is constantly trying to adjust salinities (and therefore bottom PDAs) to higher values. However, the data availability is highly variable during the investigated period and, for example, is more concentrated in the Kvarner Bay during the February-March 2015 period or along a northern Adriatic transect (Po-Rovinj) surveyed with a monthly or bimonthly frequency. This thus leads to having extremely high bottom PDAs present off the Kvarner Bay before the actual generation of the dense waters within 800 the Kvarner Bay or along the Trieste Jet in the ROMS-full model. Second, MEDSEA, contrarily to ROMS-full, uses a 3D-Var assimilation approach which is known to lose the temporal information contain in the observations through averaging (Janeković et al, 2020). In general, during the 2014/15 period, MEDSEA assimilates less data than ROMS-full which benefited from the observations collected during the NAdEx campaign.
Consequently, MEDSEA is incapable to adjust its solution in order to capture the proper dense water dynamics. For example, 805 in the Jabuka Pit, MEDSEA provides a constant decrease in bottom PDAs from autumn 2014 to winter 2015 opposite to all the other models and probably driven by the availability of the assimilated observations (e.g., Argo data). However, ROMSfull is likely to have assimilated the same observations within the Jabuka Pit but has also been assimilating Arvor -C and drifter data obtained off the Kvarner Bay during the NAdEx campaign. Further, during the winter, when bora episodes occur, only a small number of SST cloud free scenes are available for assimilation in ERA5. As a result, MEDSEA, contrary to 810 is mostly incapable to generate the bora-driven dense waters and hence to transport and collect them within the Jabuka Pit.

Conclusions
The aim of this study was to enhance our understanding of the bora-driven dense water dynamics in the Adriatic Sea using and analysing different state-of-the-art modelling approaches. The main findings of the study can be summarized as follows: • In the northern Adriatic and Kvarner Bay, dense water generation is better captured in ROMS-full and AdriSC-ROMS 815 than in MEDSEA and ROMS-hind which are producing lower volumes of dense waters. The AdriSC model generates the strongest dynamics of all the models during the bora events with the largest intensities in wind stresses, upward turbulent heat fluxes and bottom PDAs. Also, extreme dense waters are generated continuously in time and over the entire northern Adriatic in AdriSC-ROMS, while they appear as patches in ROMS-full in which a maximum is found off the southern tip of Istria, along the Senj Jet. This is linked to a combination of parameters including the 4-day 820 cycles of the 4D-Var data assimilation method used in ROMS-full and the use of atmosphere-ocean kilometre-scale models in AdriSC. Further, in the AdriSC simulation, due to the higher spatial resolution, the densest waters are collected within the Kvarner Bay where they stay for the longest amount of time.
• The transport of dense waters along the western coast is not quantitatively captured by MEDSEA and ROMS-hind.
Whereas, in the Jabuka Pit, ROMS-full collects a larger amount of dense waters than all the other models, indicating 825 that AdriSC-ROMS is probably far too dissipative. Lastly, in the SAP, the results show that the northern Adriatic dense waters did not reach the bottom of the SAP by the end of any simulation, classifying the winter of 2015 as a moderate in dense water formation over the northern Adriatic shelf.
• Impact of resolution of the atmospheric models is best seen in the ERA5 results which strongly underestimate the wind stresses. However, the heat losses are comparable between the models, but generally underestimated compared 830 to AdriSC-WRF. Concerning the hydrostatic approximation, the non-hydrostatic model AdriSC-WRF reproduces more intense wind stresses with larger spatial coverage and stronger heat losses than the hydrostatic ALADIN/HR model.
• As the Rossby radius of deformation is lower than 2 km in most of the domains during winter and spring when dense waters are generated and spreading, the differences in resolution of the ocean models and bathymetry clearly influence 835 the path and deposition of dense waters. However, it is not clear how this impacts the dense water dynamics.
• The ocean models are highly sensitive to the salinity input which plays an important role in the dense water generation.
In particular, the usage of old river climatologies causes lower salinities in ROMS-hind and MEDSEA, hence lower bottom PDAs, while AdriSC-ROMS reproduces higher salinities and PDAs.

Deleted: il
• Compared to ROMS-hind, the data assimilation in ROMS-full tends to increase the bottom PDA values in all the subdomains but particularly in the Kvarner Bay and the Jabuka Pit. Although assimilation made a large improvement of the ocean fields, the fields are reflecting initial state adjustments at the beginning of each assimilation cycle hence not producing long temporal smooth transitions. In addition, the lack of vertical resolution in the ROMS-full model 845 probably contributes to the improper representation of the dense water dynamics.
In summary, the reproduction of the dense water dynamics in the Adriatic Sea requires the use of (1) kilometre-or finerresolution atmosphere-ocean models and non-hydrostatic atmospheric models, (2) fine vertical resolutions in both atmosphere and ocean, (3) proper forcing of the open boundaries of the models, and, finally, (4) appropriate representation of the air-sea interactions (e.g., formulation of the surface wind drag). This study reveals that, if these conditions are fulfilled, models running 850 at the long temporal scales can outperform coarse resolution reanalysis products and assimilated simulations. Nevertheless, in addition to these prerequisites, 4D-Var data assimilation could be used to solve other model problemssuch as sea-surface temperature drifts, high mixing of the dense waters, etc.often found in long-term hindcasts and short-term forecasts.
However, such approach would be extremely expensive in terms of the required numerical and observational resources needed to achieve it. This study thus paved the way to a new generation of Adriatic circulation models which now should optimize 855 the accuracy of the results and the usage of the numerical resources.

Code availability
The code of the COAWST model as well as the ecFlow pre-processing scripts and the input data needed to re-run the AdriSC climate model in evaluation mode can be obtained under the Open Science Framework (OSF) data repository (Denamiel, 2021) under the MIT license. 860

Data availability
A major part of the observational data set used in this study can be obtained under the Zenodo data repository (Vilibić, 2021) under the Creative Commons by Attribution 4.0 International license. The remaining part of the observational data set is not publicly available as the data were collected within projects in which they were not publicly disseminated. The data were given for research purposes by the Institute of Oceanography and Fisheries (Croatia) upon request. 865 The model results used in this study can be obtained under the OSF data repository (Pranić, 2022)            AdriSC-ROMS simulations.