A new vision of the Adriatic Dense Water future under extreme warming

We use the Adriatic Sea and Coast (AdriSC) kilometer-scale atmosphere–ocean model to assess the impact of a far-future extreme-warming scenario on the formation, spreading, and accumulation of North Adriatic Dense Water (NAddW) over the entire basin, including the Jabuka Pit accumulation site, and Adriatic Deep Water (AdDW) over the Southern Adriatic Pit (SAP). Our key findings differ from previous studies that used coarser Mediterranean climate models and did not update the thresholds for dense-water and deep-water definitions to account for the far-future background density changes caused by warmer sea surface temperatures. We show that surface buoyancy losses at NAddW generation sites, driven by evaporation, are expected to increase by 15 % under extreme warming, despite a 25 % reduction in the intensity and spatial extent of Bora winds. As a result, future NAddW formation will remain similar to present conditions. However, the volume of dense water in the Jabuka Pit will decrease due to the increased far-future stratification. Additionally, dense-water transport between the Jabuka Pit and the deepest part of the SAP will stop, as future NAddW will be lighter than the AdDW. Regarding Ionian–Adriatic exchanges, extreme warming will not affect the impact of the bimodal oscillation system on the Adriatic salinity variability, but future AdDW dynamics will be determined by density changes in the northern Ionian Sea. Our findings highlight the complexity of climate change impacts on Adriatic atmosphere–ocean processes and the importance of high-resolution models for more accurate far-future projections in the Adriatic Sea.

1 Introduction

Dense waters, generated by extreme air–sea buoyancy losses, play a crucial role in the health and functioning of the oceans worldwide. These waters drive local and basin-wide thermohaline circulation (Broecker, 1991; Rahmstorf, 2002), ventilate deep ocean layers to support marine life, and facilitate the global carbon cycle (Emerson et al., 2001; Gruber, 2011). They transport essential nutrients (e.g., nitrogen, phosphorus, and iron) from the surface to deeper ocean layers, supporting primary production and marine ecosystems and promoting the growth of phytoplankton and other organisms at the base of the marine food web (Martin and Fitzwater, 1988; Boyd et al., 2007). Additionally, they drive vertical mixing and upwelling, which enhances biological productivity and biodiversity in surface waters (Vélez-Belchi et al., 2018; Doney et al., 2012). These processes influence regional and global climate patterns by transporting heat, moisture, and carbon dioxide across ocean basins (Rahmstorf et al., 2015; IPCC, 2019). However, with ongoing and future global warming, increased ocean stratification will inhibit the transport of heat, oxygen, and carbon dioxide from the surface to deeper layers, intensifying ocean acidification and impacting the marine food chain (Li et al., 2020).

In the Mediterranean Sea, the densest waters (observed potential density anomalies up to 30.6 kg m⁻³; Raicich et al., 2013; Mihanović et al., 2013) are formed in the northern Adriatic Sea (Fig. 1) during extreme winter windstorms known as Bora events, which produce hurricane-strength gusts up to 50 m s⁻¹ (Belušić and Klaić, 2004) and lead to significant sea surface cooling (e.g., Ličer et al., 2016). The dynamical properties of North Adriatic Dense Water (NAddW; Zore-Armanda, 1963) in the present climate have been extensively studied over the last 60 years, as summarized by Vilibić et al. (2023). NAddW formation occurs within the shallow northern Adriatic shelf (Fig. 1, NA subdomain) and Kvarner Bay (Fig. 1, KB subdomain) during strong surface heat and freshwater losses between December and March. The deepest part of Kvarner Bay (Fig. 1, DKB subdomain) also acts as a dense-water collector. NAddW then spreads southward along the western Adriatic coast and is partially collected within the Jabuka Pit (Fig. 1, JP subdomain) and the Southern Adriatic Pit (SAP; Fig. 1). In the SAP, Adriatic Deep Water (AdDW) is generated through open-ocean convection and strongly preconditioned by the presence of a permanent cyclonic gyre. Finally, the remaining NAddW and AdDW exit the Adriatic basin through the Strait of Otranto towards the northern Ionian Sea.

Figure 1 Topo-bathymetry of the AdriSC climate model with the locations of the five subdomains (colored polygons) and six transects (dotted black and white lines) used in the study. MSL stands for mean sea level.

NAddW and AdDW are the main sources of Eastern Mediterranean Deep Water (EMDW; Pollack, 1951; Malanotte-Rizzoli et al., 1997), playing a significant role in sustaining the Mediterranean overturning circulation (Li and Tanhua, 2020) and shaping the biogeochemical processes and ecosystem dynamics of the eastern Mediterranean Sea (Herut et al., 2016; Thingstad et al., 2005; Rahav and Herut, 2016). Nonetheless, the impact of climate change on their dynamical properties has not been thoroughly assessed. Soto-Navarro et al. (2020) analyzed the future evolution of deep-water formation in the Adriatic Sea with the Med-CORDEX ensemble of fully coupled regional climate models in the Mediterranean Sea, while Parras-Berrocal et al. (2023) studied the impact of climate change on dense-water formation in the eastern Mediterranean with one of the Med-CORDEX models. However, the results from both studies were averaged over the entire Adriatic Sea, and the Med-CORDEX regional climate system models (RCSMs) have coarse resolutions – 25 km in the atmosphere and about 15 km in the ocean – insufficient to represent the known NAddW dynamics accurately. Indeed, Denamiel et al. (2021b) and Pranić et al. (2023) have shown that only nonhydrostatic kilometer-scale atmospheric models and ocean models with at least 1 km resolution can properly reproduce the dense-water dynamics within the Adriatic basin.

The kilometer-scale Adriatic Sea and Coast (AdriSC) climate model (Denamiel et al., 2019) is thus used in this study. The abilities of the AdriSC model to simulate both extreme Bora events in the atmosphere and dense-water dynamics within the Adriatic basin have been assessed in the present climate, with many studies demonstrating the added value of such a kilometer-scale atmosphere–ocean climate approach (Denamiel et al., 2020a, b, 2021a, b, 2022; Pranić et al., 2021, 2023, 2024; Tojčić et al., 2023, 2024). Consequently, the present study focuses on understanding and analyzing in detail the far-future impacts of an extreme-warming scenario on the atmosphere–ocean processes driving NAddW and AdDW dynamics. The article is structured as follows. The AdriSC model and the methods used for the analyses are described in Sect. 2, while the impacts of climate change on the Bora events, Adriatic Dense Water dynamics, and Ionian–Adriatic exchanges are assessed and discussed in Sect. 3. Finally, conclusions about the main findings of the study are presented in Sect. 4.

2 Model and methods

2.1 Adriatic Sea and Coast (AdriSC) model

2.1.1 AdriSC model setup

The kilometer-scale Adriatic Sea and Coast (AdriSC) climate model (Denamiel et al., 2019) has been developed to represent the atmospheric and oceanic circulation over the Adriatic basin with greater accuracy than the available Mediterranean regional climate system models (RCSMs). It is based on the Coupled Ocean–Atmosphere–Wave–Sediment Transport (COAWST) modeling system (Warner et al., 2010), which dynamically couples the Weather Research and Forecasting (WRF; Skamarock et al., 2005) atmospheric model and the Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams, 2009). As illustrated in Fig. 2a, two nested grids of 15 km and 3 km resolution are used in the WRF model, while two nested grids of 3 km and 1 km resolution are used in ROMS. Vertical, terrain-following coordinates are used with 58 levels refined in the surface layer for the atmosphere (Laprise, 1992) and 35 levels refined both near the sea surface and the seafloor for the ocean (Shchepetkin and McWilliams, 2009).

Figure 2 Spatial coverage and horizontal resolution of the different grids used in the AdriSC climate model (a) and pseudo-global-warming temperature forcing imposed in the AdriSC extreme-warming simulation (b).

The AdriSC modeling suite is installed and fully tested on the European Centre for Medium-Range Weather Forecasts (ECMWF) high-performance computing facilities. More details on the AdriSC setup can be found in Denamiel et al. (2019, 2021a) and Pranić et al. (2021).

2.1.2 Pseudo-global-warming (PGW) approach

In this study, the impact of climate change is assessed with two 31-year-long AdriSC climate simulations: a historical run for the 1987–2017 period and a far-future extreme-warming run (2070–2100 period) based on the Representative Concentration Pathway (RCP) 8.5 (hereafter RCP 8.5 simulation). As a rapid equilibrium is reached within the AdriSC ocean models (Pranić et al., 2021), a 2-month spin-up period allowing the atmosphere–ocean models to reach a steady state is used in both simulations. For the historical run, the initial and boundary conditions are provided to the WRF 15 km model by the 6 hourly ERA-Interim reanalysis fields at 0.75° resolution (Dee et al., 2011) and to the ROMS 3 km model by the Mediterranean Forecasting System (MFS) MEDSEA reanalysis at $\mathrm{1}/\mathrm{16}\mathit{°}$ resolution (Simoncelli et al., 2019). The AdriSC historical climate run has already been successfully evaluated (Denamiel et al., 2021a; Pranić et al., 2021) and proven to reproduce the known Adriatic multi-decadal dense-water dynamics (Pranić et al., 2024).

For the RCP 8.5 run, the pseudo-global-warming (PGW) methodology is used to address both the relative slowness of the AdriSC model (i.e., a month of results produced per day) and the low temporal and spatial resolutions (i.e., few vertical levels for daily or monthly results) of the Med-CORDEX RCSMs available to force the AdriSC WRF 15 km and ROMS 3 km models. The principle of the PGW simulations (Schär et al., 1996; Denamiel et al., 2020b) is to impose an additional climatological change on the reanalysis used to force the historical run. Here, the results of the LMDZ4-NEMOMED8 RCSM (Hourdin et al., 2006; Beuvier et al., 2010) are used to produce the PGW forcing; see Denamiel et al. (2020b) for a detailed description. For the atmosphere, the ERA-Interim air temperature, relative humidity, and horizontal wind velocities, defined on 37 atmospheric pressure levels, are modified between 1000 and 70 hPa with the climatological changes ΔT, ΔRH, ΔU, and ΔV, respectively. These changes are derived from the RCP 8.5 scenario of the LMDZ4-NEMOMED8 RCSM by subtracting the atmospheric results of the 1987–2017 period from those of the 2070–2100 period, producing 6-hourly three-dimensional climatological changes for the 366 d of the year. These new forcings are then used to provide the boundary and initial conditions for the WRF 15 km model in the PGW simulation. For the ocean, the MEDSEA ocean temperature, salinity, and currents, defined on 72 unevenly spaced vertical levels, are modified with the climatological changes ΔT ocean, ΔS ocean, ΔU ocean, and ΔV ocean, respectively. These changes are also derived from the RCP 8.5 scenario of the LMDZ4-NEMOMED8 RCSM to produce three-dimensional daily climatological changes for the 366 d of the year. These forcings are then used to provide the boundary and initial conditions for the ROMS 3 km model in the PGW simulation. In other words, the same climatological changes are used to modify the boundary conditions for each simulated year of the reanalysis period, and the PGW simulations “inherit” the synoptic environment and weather and ocean conditions from the atmosphere–ocean reanalyses at the lateral boundaries. As a result, the main limitation of this methodology, compared to traditional downscaling techniques (Brogli et al., 2023), is that potential changes in intra-annual and interannual variability may be missed in the PGW projections. Additionally, in the presented RCP 8.5 simulation, due to the location of the AdriSC ROMS 3 km boundary conditions, the northern Ionian ocean dynamics may be more influenced by the MEDSEA reanalysis than by the projected climatic changes.

As illustrated in Fig. 2b, the PGW temperature forcing imposed in the AdriSC RCP 8.5 simulation is about 1 °C warmer for the air than the sea at the surface. It is also below 0.5 °C in the ocean for all depths below 1000 m but can reach up to 3.5 °C between the surface and 200 m depth for the RCP 8.5 scenario. These strong vertical gradients of temperature imposed on the ocean reanalysis are thus expected to impact the density of the Adriatic Sea, which will be far lower in the shallow areas of the basin (e.g., NA and KB subdomains) than in its deepest part (e.g., DA subdomain) in the RCP 8.5 simulation.

2.2 Methods

2.2.1 Bora events

To understand the impact of climate change on the air–sea interactions driving NAddW formation, the atmospheric results – derived from the AdriSC WRF 3 km daily fields – are only examined over the northern Adriatic Sea (for latitudes above 43° N) during extreme Bora events, defined as wind speeds at 10 m greater than 13 m s⁻¹ (i.e., gale-force winds; Belušić and Klaić, 2004). First, the validity of this simple criterion is demonstrated by analyzing the median monthly wind speeds at 10 m ($\ge \mathrm{13}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) during the 31 years of the historical simulation and comparing the obtained results with the known Bora jet dynamics (Fig. 4). Second, the impact of climate change on the selected Bora winds is assessed with spatial plots of the climate adjustments (in percent, Fig. 5) defined, during the 31 years of the simulations, as the difference between RCP 8.5 and historical median monthly wind speeds divided by the historical median monthly wind speeds. Following this, monthly climatologies are presented as time series of the median and 25th and 75th percentiles of the historical and RCP 8.5 results for eight different variables (Figs. 6 and 7): horizontal wind transport at 10 m; accumulated surface buoyancy loss; total, sensible, and latent heat fluxes; air-minus-sea-saturation-specific humidity (SAT); relative humidity at 2 m; and freshwater fluxes (see Appendix A for the mathematical definition of the variables). Finally, the results are summarized with a box plot (Fig. 7) presenting the climate adjustments (in percent) for the eight variables, which are defined as the difference between the RCP 8.5 and historical monthly results divided by the historical monthly results during the December to March (DJFM) period when NAddW is known to be formed.

2.2.2 Dense-water dynamics

In this study, all ocean variables are derived from the daily AdriSC ROMS 1 km fields, and the potential density anomalies (PDAs, σ) and the thermal expansion and haline contraction coefficients are calculated with the equation of state introduced by McDougall, Wright, Jackett, and Feistel (MWJF; Levitus and Boyer, 1994; Levitus et al., 1994; Dukowicz, 2001). Under the present climate, NAddW is characterized by densities $\mathit{\sigma }\ge \mathrm{29.2}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ (Mantziafou and Lascaratos, 2008). However, NAddW is formed in the shallowest part of the Adriatic Sea, where a strong change in background density is imposed by the PGW forcing (Fig. 2b). Therefore, this threshold cannot be used to analyze the RCP 8.5 simulation. NAddW is known to exit the Adriatic basin along the shallow western shelf of the Strait of Otranto (Fig. 1, transect T6). In Fig. 3, the historical and RCP 8.5 PDAs are presented as spatial plots of median (over the 31 years of the simulations) along the T6 transect and probability density functions – calculated with a kernel-smoothing method (Bowman and Azzalini, 1997) and evaluated for 100 equally spaced points – derived along the Strait of Otranto at the bottom of the western shelf (Fig. 3a, black box). This analysis reveals that a density of $\mathit{\sigma }=\mathrm{29.2}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ is obtained for the 97th percentile of the historical PDAs, which corresponds to $\mathit{\sigma }=\mathrm{28.4}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ in the RCP 8.5 simulation (Fig. 3b) and defines the criterion used to identify the far-future NAddW. The SAP, which is the deepest area of the Adriatic Sea, is also excluded from the following analyses.

Figure 3 Spatial plots of the vertical transect along the Strait of Otranto for the median of the PDAs over the 31-year historical and RCP 8.5 simulations (a) and historical and RCP 8.5 PDA probability density functions at the bottom of the western shelf (black box) and for depths greater than 800 m (white box) along the Strait of Otranto (b). It should be noted that historical and RCP 8.5 results are presented with different extrema in density in the color plots for a better visualization but have identical ranges (0.3 kg m⁻³) to emphasize the increased stratification.

The impact of climate change is first assessed for the isopycnal depth (ID; see Appendix A for the mathematical definition) using spatial plots of the median (over the 31 years of the simulations) of both the historical monthly maximums and the climate adjustments (in percent), defined as the difference between RCP 8.5 and historical monthly maximums divided by the historical monthly maximums (Figs. 8 and 9). Following this, daily climatologies are presented as the median and 25th and 75th percentiles of the historical and RCP 8.5 results for three different variables (Figs. 10 to 12): dense-water volume (DWV) and stratification index (SI) over four different subdomains (NA, KB, DKB, and JP; Fig. 1) and outward (i.e., exiting the Adriatic basin) NAddW mass transport along five different transects (T1 to T5; Fig. 1) defined along the known dense-water pathways (see Appendix A for the mathematical definition of the three variables). The NA and KB subdomains are geographically defined. They cover the northern Adriatic shelf (with depths below 50 m) and Kvarner Bay (with depths ranging from 0 to 100 m) and are previously identified dense-water formation sites (e.g., Zore-Armanda, 1963; Pranić et al., 2024). Transects T1 and T2 are defined along the open boundary of these subdomains. The DKB and JP subdomains are defined for depths above 80 and 200 m, respectively, and are accumulation sites. The dense waters generated in Kvarner Bay, which is much deeper than the adjacent northern Adriatic shelf, are gravitationally attracted in the DKB, while the JP is a well-researched dense-water accumulation site (e.g., Zore-Armanda, 1963; Pranić et al., 2024). Transects T3 and T4 are located north and south of the JP subdomain with the aim of properly quantifying and discriminating NAddW transported southward from the that accumulated in the Jabuka Pit. Transect T5 is located north of the deepest part of the Adriatic (SAP) to quantify how much NAddW is reaching the middle Adriatic. Finally, all of the above results are summarized with box plots (Fig. 13) of the climate adjustments (in percent), defined as the difference between the RCP 8.5 and historical daily results divided by the historical daily results during DJFM for DWV and mass transport and in December for SI, and further analyzed with historical and RCP 8.5 PDA pycnoclines along the T3 to T5 transects and within the JP subdomain. An animation of the isopycnal depth (ID) over the Adriatic Sea is also provided for the RCP 8.5 simulation (Movie S1 in the Supplement).

2.2.3 Ionian–Adriatic exchanges

Due to the strong density gradients between the shallow and deep areas of the Adriatic Sea under the RCP 8.5 scenario, most of the AdDW exchanges within the SAP are expected to occur with the northern Ionian Sea, a deep basin with depths greater than 3500 m. An analysis of the PDAs for depths below 800 m along the Strait of Otranto (Fig. 3a, white box) reveals that the present climate deep-water density threshold $\mathit{\sigma }=\mathrm{29.2}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ (Gačić et al., 2001) is obtained for the 15th percentile in the historical simulation, which corresponds to $\mathit{\sigma }=\mathrm{29.09}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ in the RCP 8.5 simulation (Fig. 3b). This defines the criterion used to identify the far-future AdDW.

Empirical orthogonal functions (EOFs) are used to compare, in space and time, the most important variability patterns in the Adriatic and northern Ionian seas for both historical and RCP 8.5 simulations. Denamiel et al. (2022) demonstrated that the long-term variability of the AdriSC model is well described by the change in sign of the main EOF components derived from sea surface height (SSH) in the northern Ionian Sea. The main modes of variability of the Ionian–Adriatic exchanges are thus derived from the AdriSC ROMS 3 km monthly northern Ionian SSH over the 31-year period of the simulations, while their impact on the Adriatic Sea is extracted from the AdriSC ROMS 1 km results. All presented spatial EOFs are obtained via a covariance matrix and are normalized. The time series of the amplitudes associated with each eigenvalue in the EOF are derived via the dot product of the data and the EOF spatial patterns, with the mean subtracted from each component time series.

The Ionian–Adriatic exchanges are also characterized with time series of both inward and outward deep-water mass transports along the Strait of Otranto (T6 transect) for the RCP 8.5 simulation and from the deep-water volume and SI within the DA subdomain for both historical ($\mathit{\sigma }\ge \mathrm{29.2}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ criterion) and RCP 8.5 ($\mathit{\sigma }\ge \mathrm{29.09}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ criterion) simulations. The DA subdomain is defined for depths above 1000 m and encompasses the SAP identically to the study of Pranić et al. (2024). Additionally, an animation of the isopycnal depth for deep water in the southern Adriatic Sea for the RCP 8.5 simulation is provided (Movie S1).

3 Results

3.1 Bora events

3.1.1 Spatial extent and intensity

The spatial extent and intensity of the selected windstorms and their associated climate adjustments are first presented in Figs. 4 and 5. Between November and March in the historical simulation (Fig. 4a–c, k, and l), the horizontal wind speeds vary from north to south along the eastern Adriatic coast, with intense jets above 16 m s⁻¹ separated by weaker speeds. This behavior is characteristic of the known Bora wake and gap jet dynamics (e.g., Jiang and Doyle, 2005; Gohm et al., 2008; Alpers et al., 2009; Signell et al., 2010), and the historical simulation can reproduce the known Trieste, Senj, Karlobag, and Sukošan main Bora jets (see Fig. 4). For the rest of the year, in accordance with the known literature, some Bora jets are still present, but their intensity is on average decreased ($\le \mathrm{14}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$). Consequently, using wind speeds $\ge \mathrm{13}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ to identify Bora winds is a simple but efficient criterion, as windstorms are dominated by these events over the northern Adriatic Sea.

Figure 4 Historical monthly climatology of the spatial extent, intensity, and direction (as vectors) of the selected windstorms ($\ge \mathrm{13}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) defined as the monthly median wind speed at 10 m over the 31 years of the historical simulation.

Figure 5 Monthly climatology of the climate adjustments (in percent) associated with the selected Bora events (wind speeds at 10 m $\ge \mathrm{13}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) and defined, over the 31 years of the simulations, as the difference between RCP 8.5 and historical monthly medians divided by historical monthly medians.

In terms of climate adjustments (Fig. 5), the far-future intensity of the Bora jets during DJFM is mostly reduced by about 5 % within Kvarner Bay but increased by up to 15 % (less than 5 % on average) along the Trieste jet. Additionally, in October and December, the intensity of the selected Bora winds is overall increased by 5 % to 15 %. Finally, for most months, an alternation of strong reduction in intensity (up to 15 % and 10 % on average) and moderate increase (up to 15 % but 5 % on average) along the known Bora jet locations can be seen in Fig. 5. Consequently, in the RCP 8.5 simulation, the locations of the main Bora jets are most probably shifted in space while their intensity is reduced overall.

3.1.2 Monthly climatologies

For the selected Bora events, the impact of climate change on the air–sea dynamics is presented separately for the historical and RCP 8.5 simulations as monthly climatologies of horizontal wind transports; accumulated surface buoyancy losses; total, latent, and sensible heat fluxes (Fig. 6); and air-minus-sea-saturation-specific humidity (SAT), relative humidity, and freshwater fluxes (Fig. 7).

Figure 6 This figure shows the selected Bora events, monthly climatologies of the median and 25th and 75th percentiles of the horizontal wind transport at 10 m (a); the accumulated surface buoyancy loss (b); and the total (c), latent (d), and sensible (e) heat fluxes defined over the 31 years for the historical, RCP 8.5, and RCP 8.5 minus historical conditions.

Figure 7 This figure shows the selected Bora events, monthly climatologies of the median and 25th and 75th percentiles of the air-minus-sea-saturation-specific humidity (a), the relative humidity at 2 m (b), and the freshwater flux (c) defined over the 31 years of the historical, RCP 8.5, and RCP 8.5 minus historical conditions. Climate adjustments (in percent) for the eight variables used in the Bora event analyses are presented as box plots during DJFM (d).

For both historical and RCP 8.5 simulations, the strongest horizontal wind transports (median value above $\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{9}}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$) occur between November and March. Compared to the historical results, the horizontal wind transports are overall reduced in the RCP 8.5 simulation, i.e., between $\mathrm{0.02}\times {\mathrm{10}}^{\mathrm{9}}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$ in February and $\mathrm{0.28}\times {\mathrm{10}}^{\mathrm{9}}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$ in March, but increased in January and August by about 0.17×10⁹ and $\mathrm{0.14}\times {\mathrm{10}}^{\mathrm{9}}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$, respectively. In terms of the most extreme wind transports, defined as the 75th percentile, they are reduced in the RCP 8.5 simulation by up to $\mathrm{3.07}\times {\mathrm{10}}^{\mathrm{9}}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$ in January and $\mathrm{1.25}\times {\mathrm{10}}^{\mathrm{9}}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$ in February but are increased by up to $\mathrm{0.82}\times {\mathrm{10}}^{\mathrm{9}}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$ in September and $\mathrm{1.33}\times {\mathrm{10}}^{\mathrm{9}}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$ in December.

The strongest accumulated surface buoyancy losses occur between September and March in both historical and RCP 8.5 simulations and can reach a monthly median of more than 0.030 m² s⁻² in December. In contrast with the horizontal wind transports, the median RCP 8.5 accumulated surface buoyancy losses are overall increased by 0.004 m² s⁻² on average compared to the historical simulation. This increase varies between 0.001 m² s⁻² in December and 0.017 m² s⁻² in November. The extreme RCP 8.5 buoyancy losses, defined as the 75th percentile, are also increased all year long by 0.005 m² s⁻² on average and by a maximum of 0.015 m² s⁻² in November.

Regarding the total, latent, and sensible monthly heat fluxes, both RCP 8.5 and historical simulations reach their maximum losses (median value above 150 W m⁻²) between September and March. Overall, compared to the historical results, the RCP 8.5 total heat losses increase between 11 W m⁻² in February and 88 W m⁻² in November, with an average of 35 W m⁻² between August and March, while the RCP 8.5 total heat gain decreases by about 17 W m⁻² on average between May and July. In contrast, the RCP 8.5 latent heat losses increase all year long by at least 3 W m⁻² in April and up to 72 W m⁻² in November (an average of 33 W m⁻²), while the RCP 8.5 sensible losses decrease by 6 W m⁻² on average most of the year (except in March, October, and November, which have increased losses between 4 and 9 W m⁻²). In terms of extremes, defined as the 25th percentile, the monthly RCP 8.5 latent heat losses are increased by 40 W m⁻² on average and up to 42 W m⁻² in March and 52 W m⁻² in November.

For the remaining variables (Fig. 7), while the RCP 8.5 median monthly relative humidity changes by less than ±1 % compared to the historical simulation, the median and extreme (represented by the 25th percentile) monthly losses of both air-minus-sea SAT and freshwater flux are expected to increase all year long by an average of 0.8 g kg⁻¹ and $\mathrm{1.15}\times {\mathrm{10}}^{-\mathrm{8}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, respectively, and up to 1.5 g kg⁻¹ in November and $\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{8}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ in October, respectively.

3.1.3 Discussion

To summarize the results of the previous section, the monthly climate adjustments are presented as box plots (Fig. 7d). These reveal that the median (and extreme, given by the 25th and 75th percentiles depending on the negative and positive sign of the median, respectively) differences between RCP 8.5 and historical results are −25 % (−52 %) for the horizontal wind transports at 10 m and −4 % (−9 %) for the sensible heat fluxes but +15 % (+30 %) for the buoyancy losses, +8 % (+20 %) for the total heat fluxes, +17 % (+24 %) for the latent heat fluxes, +20 % (+27 %) for the air-minus-sea SAT, +1 % (+2 %) for the relative humidity at 2 m, and finally +17 % (+24 %) for the freshwater fluxes.

Consequently, a strong reduction in the intensity and spatial extent of the winter Bora winds is projected in the AdriSC RCP 8.5 far-future simulation. This confirms the findings of other regional- and kilometer-scale atmospheric long-term models (e.g., Benetazzo et al., 2012; Androulidakis et al., 2015; Bonaldo et al., 2017; Belušić Vozila et al., 2019) and the 3 d long AdriSC climate simulations (Denamiel et al., 2020a, b). As previously seen in Denamiel et al. (2020a), the accumulated buoyancy losses, particularly the latent heat and freshwater losses, are strongly increased (by more than 15 %) in the RCP 8.5 simulation, leading to strong cooling at the air–sea interface. In contrast with what was previously hypothesized in Denamiel et al. (2020a), the changes in relative humidity at 2 m are minor and cannot explain this increase. However, as with the latent heat losses, the air-minus-sea SAT and the freshwater losses are projected to increase by at least 17 %, and the increase in the buoyancy losses under the RCP 8.5 conditions is mainly controlled by the increased evaporation and not by the decrease in Bora wind intensity and spatial extent. Given these results, in contrast with the findings of Parras-Berrocal et al. (2023), NAddW formation is expected to be similar in far-future and present climates.

3.2 Adriatic Dense Water dynamics

3.2.1 Spatial extent and intensity

The changes in the spatial extent and intensity of the dense-water formation, propagation, and accumulation under the far-future extreme warming are first presented as spatial plots of the monthly historical ID and their associated climate adjustments (Figs. 8 and 9). Importantly, the dynamical behavior of the SAP will be discussed in Sect. 3.3 and will not be analyzed here.

Figure 8 Historical monthly climatology of the median isopycnal depth (ID; defined for $\mathit{\sigma }\ge \mathrm{29.2}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) over the 31 years of the historical simulation.

Figure 9 Monthly climatology of the climate adjustments (in percent) associated with the ID and defined, over the 31 year of the simulations, as the difference between RCP 8.5 and historical monthly medians divided by historical monthly medians.

In the historical simulation, the ID ($\mathit{\sigma }\ge \mathrm{29.2}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ criterion) reaches a maximum within the northern Adriatic shelf (between 40 m in the shallow areas and 75 m in the deepest parts) and the Kvarner Bay (above 75 m) between December and April when NAddW is formed and fills the formation sites (NA and KB subdomains). Within the deepest parts of Kvarner Bay (i.e., DKB subdomain), the ID is still about 25 m in May and decreases to below 5 m in September when no dense water is left in this accumulation site before December. In the Jabuka Pit accumulation site, the ID peaks in February and March (above 160 m) but remains above 125 m all year long in the deepest parts of the middle and western areas of the pit. Along the Italian coast (i.e., the western side of the SAP), where the dense water is known to exit the Adriatic basin, the ID peaks in December and February (above 100 m) but varies between 5 and 75 m the rest of the year when exchanges of dense water occur between the Jabuka Pit and the SAP.

In terms of climate adjustments (Fig. 9), the RCP 8.5 ID presents changes smaller than ±10 % compared to the historical simulation all year long within the northern Adriatic, except between April and May when it decreases up to 60 % in the area where the dense waters are known to exit Kvarner Bay. However, during December the RCP 8.5 ID increases by up to 100 %. Within and off Kvarner Bay the RCP 8.5 ID decreases by up to 25 % in April but increases by up to 15 % between July and September and by up to more than 100 % in October when dense water will still be present within the DKB subdomain in the RCP 8.5 simulation. It is thus expected that less NAddW is transported from Kvarner Bay between April and June. However, the biggest changes in RCP 8.5 ID (up to ±100 %) occur all year long within the Jabuka Pit, where it decreases between 10 % in June and up to 100 % in February and March, and along the western side of the SAP, where it increases between 10 % in April and 100 % in July but decreases by up to 100 % between January and March.

3.2.2 Daily climatologies

The impact of climate change on dense-water dynamics is presented separately for the historical and RCP 8.5 simulations. These are shown as monthly climatologies of dense-water volume (DWV; Fig. 10) and stratification index (SI; Fig. 11) within the NA, KB, DKB, and JP subdomains (Fig. 1) and dense-water transports (Fig. 12) along the transects T1 to T5 (Fig. 1).

Figure 10 Daily climatologies of the median and 25th and 75th percentiles of the dense-water volume (DWV) defined over the NA (a), KB (b), DKB (c), and JP (d) subdomains for the 31 years of the historical (defined for $\mathit{\sigma }\ge \mathrm{29.2}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$), RCP 8.5 (defined for $\mathit{\sigma }\ge \mathrm{28.4}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$), and RCP 8.5 minus historical conditions.

Figure 11 Daily climatologies of the median and 25th and 75th percentiles of the stratification index (SI) defined over the NA (a), KB (b), DKB (c), and JP (d) subdomains for the 31 years of the historical, RCP 8.5, and RCP 8.5 minus historical conditions.

Figure 12 Daily climatologies of the median and 25th and 75th percentiles of the dense-water transport defined along the T1 to T5 (a–e) transects for the 31 years of the historical, RCP 8.5, and RCP 8.5 minus historical conditions.

For both historical and RCP 8.5 simulations (Fig. 10), the largest DWV (defined for $\mathit{\sigma }\ge \mathrm{29.2}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ and $\mathit{\sigma }\ge \mathrm{28.4}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, respectively) occurs between December and March (without much variability between the 25th and 75th percentiles) within the NAddW formation sites (up to 5 and 3.7×10¹¹ m³ for the NA and KB subdomains, respectively) and within the DKB accumulation site (up to 0.55×10¹¹ m³). However, for the JP subdomain, this occurs only between February and May (with variability reaching 4 and 2.5×10¹¹ m³ for the historical and RCP 8.5 simulations, respectively). Compared to the historical results, the RCP 8.5 DWV is overall identical within the NA, KB, and DKB subdomains (differences below 0.1×10¹¹ m³ for both the mean and extrema) but is reduced within the Jabuka Pit by an average of less than 0.2×10¹¹ m³ but reaching an extreme (represented by the 75th percentile) of up to 1.7×10¹¹ m³ in March.

The highest values of the SI (Fig. 11) occur across all subdomains during summer (JAS) when the sea surface temperature is at its maximum. During DJFM, when NAddW is formed, the SI is always below 0.1 m² s⁻², except for the JP subdomain where it varies between 0.2 and 0.7 m² s⁻² and 0.4 and 1.1 m² s⁻² in the historical and RCP 8.5 simulations, respectively. During this period, the RCP 8.5 SI median (and extreme, represented by the 75th percentile) gains compared to the historical simulation reach up to 0.01 (0.04) m² s⁻² for NA, 0.02 (0.04) m² s⁻² for KB and DKB, and 0.2 (0.3) m² s⁻² for JP.

In terms of the NAddW transports (Fig. 12), they mostly occur between December and May outward of the formation sites: up to $\mathrm{18.0}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}$ along T1 in March and also in January for the RCP 8.5 simulation; up to $\mathrm{7.5}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}$ along T2 in February and December for the historical and RCP 8.5 simulations, respectively; and up to nearly $\mathrm{35.0}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}$ along T3 in March. In both simulations the transports towards the Ionian Sea are overall reduced along T4 compared to T3 by up to $\mathrm{10.0}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}$ in March in the historical simulation and $\mathrm{15.0}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}$ in December in the RCP 8.5 simulation. Finally, along T5, the transports are overall reduced under the RCP 8.5 conditions compared to the historical simulation by $\mathrm{5.0}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}$ and up to $\mathrm{9.0}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}$ in March.

3.2.3 Discussion

The presented results contrast with the study by Parras-Berrocal et al. (2023), which used the same threshold to define NAddW in the present and future climates and considered NAddW and AdDW as deep water without distinction due to the coarser resolution of their RCSM. They demonstrate that under far-future RCP 8.5 conditions the median accumulated buoyancy losses are expected to increase by 15 % (Fig. 7d) and NAddW formation within the NA and KB subdomains and the accumulation within the DKB subdomain are expected to remain identical. Indeed, there is no major change in median DWV (Fig. 13a). It should be noted that despite the increase between 15 % and 32 %, the median RCP 8.5 SI remains really small in these areas during December (Fig. 11 and Fig. 13a, middle panel). Furthermore, for both historical and RCP 8.5 simulations, more NAddW is transported through T3 than through T1 and T2 combined, which means some NAddW is probably formed offshore of the NA and KB sites. Under the RCP 8.5 scenario, the offshore formation of NAddW is expected to increase as the median NAddW transports increase by 13 % along T3, decrease by 9 % along T2, and do not change along T1 (Fig. 13a, right panel).

Figure 13 Climate adjustments (in percent) for the three variables used in the NAddW dynamics analyses are presented as box plots during DJFM for DWV and mass transports and in December for SI (a). Variability of the historical and RCP 8.5 pycnoclines defined between the median and 99th percentile of the PDA along the T3 to T5 transects and within the Jabuka Pit during DJFM (b).

Within the Jabuka Pit accumulation site, the main cascading and accumulation of NAddW shifts from March in the historical simulation to December under the RCP 8.5 conditions (i.e., the maximum reduction in dense-water transports between T3 and T4 that frame the JP is obtained in March under the historical conditions and in December under the RCP 8.5 conditions; Fig. 12). However, despite the increase in transports by 13 % along T3 and their decrease by 19 % along T4 (Fig. 13a, right panel), the RCP 8.5 DWV within the Jabuka Pit is reduced by 5 % compared to the historical simulation (Fig. 13a, left panel). Comparing the pycnoclines along T3, JP, T4, and T5 (Fig. 13b) reveals that there are more occurrences of NAddW filling the full water column in the historical simulation than under the RCP 8.5 conditions; i.e., above 100 m depth, the area between 29.2 kg m⁻³ and the 99th percentile (in blue) is twice as large as the area between 28.4 kg m⁻³ and the 99th percentile (in red). Over the Jabuka Pit, the increase by 40 % of the RCP 8.5 SI is thus likely to hamper vertical mixing, causing a diminution of the RCP 8.5 DWV despite the presence of RCP 8.5 NAddW varying between 28.5 and 29.2 kg m⁻³ at the bottom of the pit below 200 m depth (Fig. 13b).

Finally, the reduction in the RCP 8.5 transports by 13 % along T4 and 34 % along T5 compared to historical conditions (Fig. 13a, right panel) can have two explanations. First, the accumulated NAddW from the Jabuka Pit cannot cascade within the deepest part of the SAP where the densities are higher and, hence, in contrast with the historical results, no strong density current is present. Second, the decrease in densities of NAddW (from the Jabuka Pit to transect T5) due to the interaction with the ambient Adriatic waters is greater for the RCP 8.5 than the historical conditions, up to 0.3 and 0.2 kg m⁻³, respectively (Fig. 13b). This suggests that most NAddW exits the Adriatic Sea along the western side of the SAP under RCP 8.5 conditions, which also explains the increase in ID in this area between April and November (Fig. 9).

3.3 Ionian–Adriatic exchanges

In this section, the SSH EOFs over the northern Ionian Sea are used to define the main modes of the Ionian–Adriatic exchanges. First, for both historical and RCP 8.5 simulations, the first SSH EOFs are linked to the interannual variability and are not displayed here. Second, both an analysis and discussion of the results are presented together below.

3.3.1 Bimodal oscillation system (BiOS)

For the historical simulation, as described in Denamiel et al. (2022), the second SSH EOF, representing nearly 10 % of the total signal (Fig. 14a, left panels), is linked to the Ionian–Adriatic bimodal oscillation system or BiOS (Gačić et al., 2010). In the present climate, the BiOS connects the quasi-decadal reversals of the North Ionian Gyre (NIG) circulation to the salinity variability in the Adriatic Sea. During the anticyclonic phase of the NIG, the southern Adriatic Sea salinity decreases due to the advection of lower-salinity Atlantic Water. During the cyclonic phase of the NIG, the salinity increases due to the advection of high-salinity Levantine and eastern Mediterranean waters.

Figure 14 Normalized spatial EOF components and associated time series of amplitude for both the historical and RCP 8.5 AdriSC ROMS 3 km sea surface height (SSH) over the northern Ionian Sea (a). RCP 8.5 AdriSC ROMS 1 km salinity at 100 m and at the bottom of the Adriatic Sea (b). Time series of the daily vertical potential density anomaly (PDA) and dense-water volume (DWV) in the JP subdomain for the RCP 8.5 simulation (c).

In the RCP 8.5 simulation, the BiOS signal appears as the third SSH EOF and represents only 8 % of the total signal (Fig. 14a, right panels). From these results, the expected BiOS signal in the Ionian Sea for the RCP 8.5 scenario is similar in both spatial pattern and time series to the one obtained for the historical simulation. Furthermore, under RCP 8.5 conditions, the correlations between the BiOS signal and the first salinity EOFs at 100 m depth and the bottom of the Adriatic Sea – representing 74 % and 56 % of the total signal, respectively (Fig. 14b) – reach more than 60 % with a 2-year lag.

Consequently, as these results are similar to those found for the historical simulation by Denamiel et al. (2022), the BiOS remains the main driver of Adriatic salinity variability under the presented PGW extreme-warming scenario. Importantly, at the bottom of the Adriatic Sea, the BiOS does not affect the deepest part of the SAP (i.e., the DA subdomain).

Finally, for the RCP 8.5 simulation, the Adriatic BiOS-driven salinity phases strongly impact the renewal of the Jabuka Pit collector site (Fig. 14d): during the cyclonic phases (Fig. 14b, in blue in the EOF time series), both PDA and DWV increase (up to 29.2 kg m⁻³ and 5.0×10¹¹ m³, respectively), whereas during the anticyclonic phases (Fig. 14b, in red in the EOF time series), the PDA is largely decreased over the entire water column (down to below 28.2 kg m⁻³), while no (or very little, below 0.5×10¹¹ m³) dense water is present in the Jabuka Pit collector.

3.3.2 Deep-water exchanges

In the RCP 8.5 simulation, the second SSH EOF – representing about 15 % of the total signal over the Ionian Sea (Fig. 15a, right panels) – is a mode of Ionian–Adriatic exchanges that is not present in the historical simulation. The anti-correlation, without lag, between this new mode and the first EOF of the bottom PDA over the SAP – representing about 63 % of the total PDA signal for depths below 800 m – reaches more than 80 %. Consequently, this mode controls the deep-water content of the deepest part of the SAP (i.e., the DA subdomain) and hence the presence of AdDW under the PGW RCP 8.5 scenario. The first phase of this mode (Fig. 15a, in blue in the EOF time series of SSH over the Ionian Sea) is present for 7 years at the beginning of the RCP 8.5 simulation. The second phase (Fig. 15a, in red in the EOF time series of SSH over the Ionian Sea) lasts for 20 years, while during the last 4 years of the simulation, the mode reverts to the first phase.

Figure 15 Normalized spatial EOF components and associated time series of amplitude for the RCP 8.5 AdriSC ROMS 3 km sea surface height (SSH) over the northern Ionian Sea and the RCP 8.5 AdriSC ROMS 1 km bottom PDA for depths greater than 800 m in the Adriatic Sea (a). Time series of the daily deep-water transport along the T6 transect for the RCP 8.5 simulation (b). DWV (c) and SI (d) over the DA subdomain for the historical and RCP 8.5 simulations.

The switch between the first and second phase of this new mode corresponds to the year 1994 in the historical simulation, which marks the shift in dominant deep-water source in the northern Ionian Sea from the Adriatic Sea to the Aegean Sea. In the historical simulation, this event – known as the Eastern Mediterranean Transient (EMT) – is characterized by massive dense-water formation triggered by extreme heat losses and high salinity in the Aegean Sea during winter 1992–1993 (Roether et al., 1996, 2007; Klein et al., 1999; Velaoras et al., 2017). During the EMT, the northern Ionian Sea is filled with very dense water from the Aegean Sea, and the intrusion of water originating from the Adriatic into the Levantine basin is blocked (Akpinar et al., 2016; Li and Tanhua, 2020). As the PGW method uses the historical boundary forcing, the RCP 8.5 scenario presented in this study is also forced with the EMT signal modified with an extreme-warming climatological change.

In the RCP 8.5 simulation, for $\mathit{\sigma }\ge \mathrm{29.09}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, no deep water is present within the DA subdomain before the EMT. First, most NAddW (defined for $\mathit{\sigma }\ge \mathrm{28.4}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) is too light to cascade into the deepest part of the SAP (which has an ambient density of about 29.0 kg m⁻³ before the EMT). Second, over the DA subdomain, the RCP 8.5 stratification index (SI) is multiplied by at least 7 compared to the historical conditions (Fig. 15d), which greatly hampers the far-future deep-ocean convection. This strongly contrasts with the present climate conditions for which NAddW is known to partly transform into AdDW during deep-convection processes over the SAP. However, after 7 years of simulation, the EMT triggers new Ionian–Adriatic exchanges of deep water, and the DA subdomain is filled with deep water coming from the northern Ionian Sea – i.e., inward transports and DWV within the DA subdomain up to $\mathrm{10.0}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}$ and 5.0×10¹² m³, respectively (see Movie S1 and Fig. 15b).

Figure 16 Visual summary of the study. Adriatic Dense Water and Adriatic Deep Water far-future dynamics as seen by the kilometer-scale atmosphere–ocean AdriSC model under the pseudo-global-warming assumption. The numbers and letters correspond to the description of the different findings as described in Sect. 4.

In the far-future simulation, the amount of dense water within the DA subdomain is thus controlled by the Ionian–Adriatic exchanges and is far lower than under the historical conditions (Fig. 15c). Under present climate conditions, the ventilation of the deepest part of the SAP by NAddW is indeed known to occur regularly (Cardin et al., 2020). This is marked by strong peaks in DWV occurring every 3 to 5 years over the DA subdomain in the historical simulation and their absence in the RCP 8.5 simulation.

Interestingly, these results can be compared to the study of Soto-Navarro et al. (2020), who found that most Med-CORDEX models project a reduction in the intensity of the deep-convection events, while one model projects an intensification of the convection in the Aegean Sea similar to what happened during the EMT in the 1990s. However, in the Aegean, Soto-Navarro et al. (2020) found that most Med-CORDEX models project a reduction in the dense water formation, and hence the EMT-like situation seen in the AdriSC model under the PGW assumption is unlikely to occur.

4 Conclusions

In this study, an analysis of the dynamics of North Adriatic Dense Water (NAddW) and Adriatic Deep Water (AdDW) is presented using the kilometer-scale atmosphere–ocean AdriSC model under the pseudo-global-warming (PGW) assumption. Several findings differing from previous studies based on coarser Mediterranean climate models are revealed and summarized in Fig. 16 and as follows.

First, employing PGW forcing in the far-future simulation and thus imposing a strong vertical temperature gradient to the AdriSC ROMS 3 km initial and boundary conditions clearly emphasizes the necessity to update thresholds for defining dense and deep waters to account for background density changes. In fact, this result is aligned with the supplementary study done by Parras-Berrocal et al. (2023), which demonstrates that the choice of threshold significantly influences the results of dense-water formation within the Adriatic basin. However, this study reduces the NAddW threshold to 28.4 kg m⁻³, while the lowest threshold tested by Parras-Berrocal et al. (2023) is 28.8 kg m⁻³.

Second, analysis of air–sea interactions at NAddW generation sites demonstrates a 15 % increase in winter surface accumulated buoyancy losses under extreme warming. This finding contrasts with previous studies that did not reproduce the changes in coastal evaporation which compensate for the well-known reduction in intensity and spatial extent of far-future Bora winds (found to be, on average, 25 % in this study).

Third, as a consequence of the first two points, the major finding of this study is that far-future NAddW formation under extreme warming is expected to remain similar to present conditions. However, in terms of NAddW transport and accumulation, the volume of dense water in the Jabuka Pit is projected to decrease due to higher stratification hampering the vertical mixing, while transports between the Jabuka Pit and the deepest part of the Southern Adriatic Pit (SAP) are expected to stop because NAddW will be lighter than AdDW in the far future.

Fourth, the deepest part of the Adriatic basin is found to be mostly disconnected from the NAddW dynamics, and the far-future AdDW dynamics are expected to depend on density changes in the northern Ionian Sea. However, the presented RCP 8.5 AdDW results are strongly influenced by boundary conditions imposed on the AdriSC ROMS 3 km in the northern Ionian Sea and using the PGW methodology. Consequently, exchanges between the northern Ionian Sea and the deepest part of the SAP should be further investigated, for example, with kilometer-scale models capable of properly representing (and sampling) the Strait of Otranto and having open boundaries away from the Ionian Sea.

Finally, under the PGW assumption, this study finds that extreme warming is unlikely to affect the impact of the bimodal oscillation system (BiOS) on salinity variability in the Adriatic basin. However, similar to the previous point, the impact of extreme warming on the BiOS itself (i.e., on the reversal of the North Ionian Gyre) is likely not captured by the AdriSC modeling suite due to the strong influence of boundary conditions imposed in the northern Ionian Sea.

Beyond the presented results, dense-water formation in the Adriatic Sea plays a crucial role in sustaining a variety of marine species, ranging from deep-sea corals (Cushman-Roisin et al., 2001; Grubišić et al., 2014) to shallow-water mussels (Ballarin and Frizzo, 2004), sea urchins (Pais et al., 2012), and seagrasses (Boudouresque et al., 2009). Pelagic species such as European pilchards and anchovies also benefit indirectly from nutrient upwelling caused by dense-water formation, which increases plankton availability (Santojanni et al., 2006a, b). Currently, the impact of climate change on these species has not been comprehensively studied in the Adriatic Sea but has only been addressed in global assessments (e.g., Wernberg et al., 2011; Doney et al., 2012; Bopp et al., 2013). As demonstrated in this study, which provides new insights into far-future NAddW dynamics, climate change impacts on Adriatic atmosphere–ocean processes are highly complex and necessitate the use of high-resolution models. These processes also influence the biogeochemistry of the Adriatic basin, suggesting that this study may pave the way for new assessments of the impact of extreme warming on ecology and fisheries in the Adriatic Sea.

Appendix A

A1 Atmospheric variables

U10	horizontal wind speed at 2 m	[m s−1]
Ua	horizontal wind speed at 2 m	[m s−1]
Ta	air temperature at 2 m	[°C]
rh	relative humidity at 2 m	[%]
ρa	density of moist air at 2 m	[kg m−3]
ρw	density of freshwater	[kg m−3]
Pa	mean sea level pressure	[hPa]
Cpa=1004.67	specific heat capacity of the air	[$\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{kg}}^{-\mathrm{1}}$]

A2 Ocean variables

Ts	sea surface temperature	[°C]
TsK	sea surface temperature	[K]
Ss	sea surface salinity	[PSU]
g=9.81	gravitational acceleration	[m s−2]
ρ0=1025	reference density of seawater	[kg m−3]
Cp0=3991.87	specific heat capacity of seawater	[$\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{kg}}^{-\mathrm{1}}$]
α	thermal expansion coefficient	[K−1]
β	haline contraction coefficient	[PSU−1]
σ	potential density anomaly (PDA)	[kg m−3]
σT	PDA threshold for dense or deep waters	[kg m−3]
ρ	density of the seawater	[kg m−3]

A3 Horizontal wind transport

The horizontal wind transport is defined as the integration of gale-force winds (i.e., horizontal wind speeds at 10 m greater than 13 m s⁻¹) over the area where they occur in the northern Adriatic Sea (for latitudes above 43° N). In this study, the monthly median of the horizontal wind transports is used as a proxy to quantify the impact of climate change on the intensity and spatial extent of the extreme Bora events.

$T_{\text{wind}}=\iint U_{\mathrm{10}}\left(U_{\mathrm{10}}\ge \mathrm{13}\right)\mathrm{d}x\mathrm{d}y$

[m3 s−1]

A4 Total, sensible, latent heat, and freshwater fluxes

For comparison purposes, the air–sea fluxes are calculated in the same way as in Denamiel et al. (2020a, b). In this study, the monthly medians of the total, latent, and sensible heat fluxes; relative humidity at 2 m; air-minus-sea-saturation-specific humidity (SAT); and freshwater fluxes are used to quantify the impact of climate change on the air–sea interactions over the northern Adriatic during extreme Bora events (U₁₀≥13).

Qswn	net shortwave radiation	[W m−2]
Qlwd	downward longwave wave radiation	[W m−2]
ϵ=0.97	infrared emissivity
$\begin{array}{rl}{\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\text{Stef-Bolt}}\\ {\displaystyle }& {\displaystyle }=\mathrm{5.670374419}\times {\mathrm{10}}^{-\mathrm{8}}\end{array}$	Stefan–Boltzmann constant	[$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{4}}$]
$\begin{array}{rl}{\displaystyle }& {\displaystyle }{Q}_{\text{lwn}}\\ {\displaystyle }& {\displaystyle }=Q_{\text{lwd}}-\mathit{ϵ}{\mathit{\sigma }}_{\text{Stef-Bolt}}{T}_{\text{sK}}^{\mathrm{4}}\end{array}$	net longwave radiation	[W m−2]
esat(T)	saturation vapor pressure	[hPa]
$L\left(T\right)=\mathrm{2}\mathrm{501}\mathrm{000}-\mathrm{2370}T$	latent heat of vaporization	[J kg−1]
CH,CE	turbulent transfer coefficients

$\begin{array}{rl}{\displaystyle }& {\displaystyle }{q}_{\mathrm{a}}\\ {\displaystyle }& {\displaystyle }\approx {\displaystyle \frac{\mathrm{0.62197}\left(\mathrm{0.01}{r}_{\mathrm{h}}{e}_{\text{sat}}\left(T_{\mathrm{a}}\right)\right)}{p_{\mathrm{a}}}}\end{array}$	air-saturation-specific humidity at 2 m	[kg kg−1]
$\begin{array}{rl}{\displaystyle }& {\displaystyle }{q}_{\mathrm{s}}\\ {\displaystyle }& {\displaystyle }\approx {\displaystyle \frac{\mathrm{0.62197}\left(\mathrm{0.98}{e}_{\text{sat}}\left(T_{\mathrm{s}}\right)\right)}{p_{\mathrm{a}}}}\end{array}$	sea-surface-saturation-specific humidity	[kg kg−1]
$Q_{\mathrm{s}}={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{H}}{C}_{p\mathrm{a}}{U}_{\mathrm{a}}\left(T_{\mathrm{a}}-T_{\mathrm{s}}\right)$	sensible heat flux	[W m−2]
$\begin{array}{rl}{\displaystyle }& {\displaystyle }{Q}_{\mathrm{l}}\\ {\displaystyle }& {\displaystyle }={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{E}}{U}_{\mathrm{a}}L\left(T_{\mathrm{s}}\right)\left(q_{\mathrm{a}}-q_{\mathrm{s}}\right)\end{array}$	latent heat flux	[W m−2]
$E_{\mathrm{v}}=\frac{{\mathit{\rho }}_{\mathrm{a}}}{{\mathit{\rho }}_{\mathrm{w}}}C{}_{\mathrm{E}}{U}_{\mathrm{a}}\left(q_{\mathrm{a}}-q_{\mathrm{s}}\right)$	evaporation rate	[m s−1]
Pr	precipitation rate	[m s−1]
$\begin{array}{rl}{\displaystyle }& {\displaystyle }{Q}_{\text{Total}}\\ {\displaystyle }& {\displaystyle }=Q_{\text{swn}}+Q_{\text{lwn}}+Q_{\mathrm{s}}+Q_{\mathrm{l}}\end{array}$	total heat fluxes	[W m−2]
$\text{FWF}=P_{\mathrm{r}}-E_{\mathrm{v}}$	freshwater fluxes over the sea	[m s−1]

A5 Surface buoyancy fluxes and losses

In this study, the monthly accumulated surface buoyancy losses (BLs) are used to quantify the air–sea fluxes over the northern Adriatic Sea during extreme Bora events (U₁₀≥13). The surface buoyancy fluxes (BFs) are defined as in Parras-Berrocal et al. (2023), but the buoyancy losses (BLs) are calculated monthly instead of over the DJFM period.

$\text{BF}=g\left(\frac{\mathit{\alpha }{Q}_{\text{Total}}}{{\mathit{\rho }}_{\mathrm{0}}{C}_{\text{pw}}}+\mathit{\beta }{S}_{\mathrm{S}}\text{FWF}\right)$	surface buoyancy fluxes	[m2 s−3]
$\text{BL}=-\int \text{BF}\mathrm{d}t$	monthly surface buoyancy losses	[m2 s−2]

A6 Isopycnal depth and volume for dense and deep water

The isopycnal depth (ID) for dense and deep waters is calculated over the vertical for a specific isopycnic surface (σ_T) assuming that the water column is stable. It is used to derive the dense-water or deep-water volume (DWV) quantifying the amount of dense and deep water present within the specific subdomains chosen in this study.

${\mathit{\delta }}_{\mathit{\sigma }{\mathit{\sigma }}_{\mathrm{T}}}=\left\{\begin{array}{ll}\mathrm{0}& \text{if}\mathit{\sigma }<{\mathit{\sigma }}_{\mathrm{T}}\\ \mathrm{1}& \text{if}\mathit{\sigma }\ge {\mathit{\sigma }}_{\mathrm{T}}\end{array}\right.$	Kronecker delta
$\text{ID}=\int {\mathit{\delta }}_{\mathit{\sigma }{\mathit{\sigma }}_{\mathrm{T}}}\mathrm{d}z$	isopycnal depth	[m]
$\text{DWV}=\iint \text{ID}\mathrm{d}x\mathrm{d}y$	dense- or deep-water volume	[m3]

A7 Stratification index

The stratification index (SI; Turner, 1973) is used in this study to assess the daily water column stratification (i.e., low values indicate a weak stratification and vice versa). For comparison purposes, the same vertical integration is used as in Parras-Berrocal et al. (2023), but the SI is defined as the median value over the specific subdomains chosen in this study and not over the whole Adriatic Sea.

$N^{\mathrm{2}}=\frac{g}{{\mathit{\rho }}_{\mathrm{0}}}\frac{\partial \mathit{\rho }}{\partial z}$	with N the Brunt–Väisälä frequency	[s−2]
$\text{SI}={\int }_{\mathrm{0}}^h{N}^{\mathrm{2}}z\mathrm{d}z$	stratification index with h=650 m	[m2 s−2]

A8 Dense- or deep-water inward and outward transports along a vertical transect T

The dense- and deep-water transports are calculated along the transects selected in this study and can be outward transports (i.e., exiting the Adriatic basin) or inward transports (i.e., entering the Adriatic basin).

UN	ocean velocity normal to the transect T	[m s−1]
xT	distance along the transect T	[m]
$\begin{array}{rl}{\displaystyle }& {\displaystyle }{M}_{\mathrm{T}\mathit{\_}\text{outwards}}\\ {\displaystyle }& {\displaystyle }=\iint \mathit{\sigma }\left(\mathit{\sigma }\ge {\mathit{\sigma }}_{\mathrm{T}}\right)\\ {\displaystyle }& {\displaystyle }{U}_{\mathrm{N}}\left(U_{\mathrm{N}}\le \mathrm{0}\right)\mathrm{d}{x}_{\mathrm{T}}\mathrm{d}z\end{array}$		[kg s−1]
$\begin{array}{rl}{\displaystyle }& {\displaystyle }{M}_{\mathrm{T}\mathit{\_}\text{inwards}}\\ {\displaystyle }& {\displaystyle }=\iint \mathit{\sigma }\left(\mathit{\sigma }\ge {\mathit{\sigma }}_{\mathrm{T}}\right)U_{\mathrm{N}}\\ {\displaystyle }& {\displaystyle }\left(U_{\mathrm{N}}\ge \mathrm{0}\right)\mathrm{d}{x}_{\mathrm{T}}\mathrm{d}z\end{array}$		[kg s−1]