Projected climate oligotrophication of the Adriatic marine ecosystems

Mentaschi, Lorenzo; Lovato, Tomas; Butenschön, Momme; Alessandri, Jacopo; Aragão, Leonardo; Verri, Giorgia; Guerra, Roberta; Coppini, Giovanni; Pinardi, Nadia

doi:10.3389/fclim.2024.1338374

ORIGINAL RESEARCH article

Front. Clim., 02 February 2024

Sec. Predictions and Projections

Volume 6 - 2024 | https://doi.org/10.3389/fclim.2024.1338374

This article is part of the Research Topic Coastal Climate Monitoring and Downscaling for Adaptation Planning in the Adriatic Sea View all 11 articles

Projected climate oligotrophication of the Adriatic marine ecosystems

$\r\nLorenzo Mentaschi$ Lorenzo Mentaschi¹^*

Giorgia Verri³^†

Nadia Pinardi¹^†

¹Department of Physics and Astronomy “Augusto Righi” (DIFA), University of Bologna, Bologna, Italy
²Fondazione Centro Euro-Mediterraneo sui Cambiamenti Climatici (CMCC), Bologna, Italy
³Fondazione Centro Euro-Mediterraneo sui Cambiamenti Climatici (CMCC), Lecce, Italy

The Adriatic Sea hosts diverse marine ecosystems, characterized by rich biodiversity and unique ecological dynamics. Its intricate coastal habitats and open waters support a range of species and contribute to the region's ecological and economic significance. Unraveling the consequences of the ongoing climate changes on this delicate environment is essential to ensure the future safeguard of this basin. To tackle this problem, we developed a biogeochemical model for the entire basin, with a horizontal resolution of about 2 km and 120 vertical levels, forced by the projections of atmosphere, hydrology and ocean circulation between 1992 and 2050, under emission scenario RCP8.5. The changes projected between 2031–2050 and 1992–2011 were evaluated on ecoregions characterized by different trophic conditions, identified using a k-medoid classification technique. The results point toward a generalized oligotrophication of the basin, especially intense in the northern estuarine areas, driven by a substantial decrease in river discharge projected for the rivers of the Po Plain. This scenario of unproductive and declining resources, together with the ongoing warming, salinization, and acidification of marine waters, cast doubt on the long-term resilience of the Northern Adriatic food web structure, which has evolved to thrive in high trophic conditions. The outcome of this study provides the stakeholders with a tool to understand how potential long-term decreases in the regimes of the Northern Adriatic Rivers could affect the marine ecosystem and its goods and services in the future.

Introduction

Despite its relatively small size, the Adriatic Sea hosts a variety of marine ecosystems, that support a range of habitats and species. The Northern Adriatic Rivers carry large nutrient loads toward the sea, resulting in a nutrient-rich environment (e.g., Solidoro et al., 2009) on the western side while the north-eastern regions are dominated by oligotrophy (Gilmartin et al., 1990; Polimene et al., 2006, 2007; Mozetič et al., 2010). The middle-eastern coasts are characterized by relatively small karstic rivers, and nutrient scarcity, which induce prevalent oligotrophic conditions (e.g., Vilicic et al., 2002). Finally, the southern part is characterized by relevant estuaries and upwelling in the East, producing a shelf nutrient-rich area (Marini et al., 2010), and by oligotrophic conditions at the surface layers of the deep waters (e.g., Zavatarelli et al., 2000; Polimene et al., 2006), though convection events can create phytoplankton blooms at the surface during winter (e.g., Civitarese et al., 2005; Bensi et al., 2013; Najdek et al., 2014; Korlević et al., 2015; Batistić et al., 2019).

Such variety allows for the presence of a range of different biomes. The north-western coastal areas host rich and productive habitats, with abundant and diversified biocenosis. For example, seagrass meadows and associated grazers (e.g., Facca et al., 2014), or sandy, nutrient-rich benthic areas, important for burrowing animals, such as worms and bivalves (e.g., Fedra et al., 1976; Scroccaro et al., 2022). The pelagic zone is productive and rich in biomass, harboring a rich ichthyofauna, and making the Adriatic Sea one of the richest areas in the Mediterranean Sea (Giani et al., 2012; Lipej et al., 2022). Such conditions reflect flourishing economic activities in the area such as fishing (Bastardie et al., 2017), shellfish farming (e.g., Brigolin et al., 2017) and harvesting (e.g., Viaroli et al., 2005).

The more oligotrophic environment of the eastern and southern areas results in pristine and nutrient-poor waters, where biological communities adapted to nutrient scarcity can thrive. Such conditions are often regarded as limiting the risk of Harmful Algal Blooms (HAB), which can have damaging impacts on ecology and coastal economic activities (e.g., Sellner et al., 2003). Indeed, the north-eastern part of the Adriatic Sea, which is prevalently oligotrophic, is not subject to recurrent eutrophication events since the 1980s (e.g., Faganeli et al., 1985; Giani et al., 2012; Kralj et al., 2019; Grilli et al., 2020), unlike the more eutrophic north-western portion of the basin (e.g., Alvisi and Cozzi, 2016; Brush et al., 2020). It is noteworthy that HABs are extreme events whose frequency cannot be inferred from the sole knowledge of the average conditions. For example Totti et al. (2019) suggest that more frequent HAB episodes are possible amid mean tendencies toward oligotrophication. On the other hand, the low levels of nutrients reduce the overall productivity of the marine ecosystem and can hinder the sustainability of high levels in specific economic sectors, such as fishing or clam harvesting. Moreover, the widespread practice of aquaculture in these areas can have negative consequences on the environment, such as pollution by the release of excess nutrients, or the aggregation of communities of wild fish in the proximity of aquaculture sites (e.g., Tičina et al., 2020).

Climate change poses a threat to the delicate physical and biogeochemical equilibria underpinning this ecological diversity. Rising temperatures can alter the ecophysiology of marine organisms, potentially leading to alterations in the nutrient cycles and to a decrease in nutrients availability. Warmer waters come with decreased oxygen solubility, which can expand the oxygen minimum zone and cause further stress on marine organisms (Keeling et al., 2010). Past studies suggest an imbalanced response to warming of photosynthesis and respiration, with respiration rates raising more rapidly than primary productivity, potentially leading to a decrease in biomass (e.g., O'Connor et al., 2009; Boscolo-Galazzo et al., 2018; Kralj et al., 2019). Increasing concentration of CO₂ in the atmosphere causes a decline in ocean pH (ocean acidification, e.g., Gregor and Gruber, 2021; Ma et al., 2023), with potential wide-range impacts on marine organisms (Doney et al., 2009). The magnitude of these impacts differs between taxa. For example, corals, coccolithophores, and mollusks exhibit the highest declines in calcification, with significant survival reduction in mollusks (Kroeker et al., 2013). Climate projections suggest that the global ocean will be affected by increasing stratification, with consequent oligotrophication and reduction in productivity (e.g., Fu et al., 2016). Very important for the Adriatic Sea, shifting precipitation and evapotranspiration patterns and depletion of alpine glaciers could alter the hydrological regimes of the river catchments, with substantial changes in the discharge of freshwater and nutrient loads (e.g., Vezzoli et al., 2015). A major threat to the Adriatic ecosystem is human-induced changes in river discharges, connected to agricultural practices and climate precipitation trends. Alterations in river patterns are known to significantly impact nutrient loads and marine biogeochemistry in regions affected by river estuaries (e.g., Cozzi et al., 2018). Nevertheless, the examination of such effects in the context of future climate changes is a new endeavor for the Adriatic Sea. This study marks one of the initial efforts, incorporating river discharge projections and employing a climate downscaling model in this exploration.

From a point of view of physical modeling, high-resolution is needed to account for the complex orography/bathymetry, and the model must be set up for the correct reproduction of mesoscale/submesoscale dynamics (e.g., Zavatarelli and Pinardi, 2003). Furthermore, a careful representation of the riverine release plays a pivotal role in both physics and biogeochemistry in the Adriatic Basin (Vilibić et al., 2016; Verri et al., 2018, 2024). From the perspective of biogeochemistry, the Northern Adriatic is extremely sensitive to the balance of freshwater and nutrients, as well as to the southward transport of the constituents by the western coastal currents (e.g., Spillman et al., 2007). A further challenge involves the need to simulate vastly different conditions within a single biogeochemical model, notably the contrast between nutrient-rich conditions in the northwest and oligotrophic conditions in other areas.

The last decades saw the development of 3D models coupling transport and biogeochemistry at the Adriatic or sub-Adriatic scale, with increasing degree of accuracy (e.g., Zavatarelli et al., 2000; Polimene et al., 2007; Spillman et al., 2007; Scroccaro et al., 2022). However, up to now, predictions/forecasts of future changes have only been derived from (a) pioneering research using 1D models at certain sites (Vichi et al., 2003), and (b) larger-scale 3D models that lack the resolution needed to accurately describe the intricate circulation patterns and biogeochemical dynamics of the Adriatic Sea (e.g., Teruzzi et al., 2014; Macias et al., 2018; Richon et al., 2019; Butenschön et al., 2021; Reale et al., 2022; Tojčić et al., 2023).

To fill in some of the gaps, within the AdriaClim project we developed a high-resolution 3D biogeochemical model at the Adriatic scale, forced by atmospheric, circulation and hydrological data of a limited area climate model (Verri et al., 2024) downscaled from a Med-CORDEX (Ruti et al., 2016) simulation under Representative Concentration Pathway (RCP) scenario RCP 8.5 (van Vuuren et al., 2011). We then examined the model's results for two time periods, the past (1992–2011) and the future (2031–2050), considering a set of Essential Ocean Variables (EOVs, Muller-Karger et al., 2018), in particular temperature, oxygen concentration, concentration of chlorophyll-a, inorganic nitrate and phosphate, phytoplankton biomass, plus pH and Net Primary Productivity (NPP). In the reminder of the manuscript the model is described, together with its ability to accurately reproduce the biogeochemical dynamics, and the projected changes in relevant Adriatic ecoregions are discussed.

Materials and methods

Model setup

The AdriaClim biogeochemical modeling system consists of the Biogeochemical Fluxes Model (BFM) (Vichi et al., 2007), version 5.2 (Vichi et al., 2020), coupled with the offline transport module Tracer in Ocean Paradigm (TOP) of the Nucleus for European Modeling of the Ocean (NEMO), version 3.6 (Madec and NEMO team, 2017; Aumont et al., 2018).

BFM is a well-established open-source community model, based on the living and non-living Functional Groups approach, and with variable stoichiometry, capable of reproducing the main features of the pelagic and benthic lower trophic layers and biogeochemical cycles. BFM operates by describing the fluxes of key Chemical Functional Families (Carbon, Nitrogen, Phosphorus, Silica, Chlorophyll-a) between different living and non-living Functional Groups (FG) as well as with dissolved inorganic nutrients (Nitrate ${NO}_{3}^{-}$ , Ammonium ${NH}_{4}^{+}$ , Phosphate ${PO}_{4}^{3 -}$ , Silicate ${SiO}_{4}^{4 -}$ ). Our implementation uses a layout similar to Lazzari et al. (2016) and Cossarini et al. (2021), with 4 living FGs of phytoplankton (diatoms, autotrophic nanoflagellates, picophytoplankton, large phytoplankton), 4 living FGs describing mesozooplankton (omnivorous and carnivorous mesozooplankton) and microzooplankton (microzooplankton and heterotrophic nanoflagellates), one bacterial FG, and 4 non-living FGs for the organic matter (labile and semi-labile Dissolved Organic Matter – DOM, Particulate Organic Matter – POM, and inorganic matter). Altogether, the biogeochemical dynamics of each pelagic cell is described employing 51 state variables. Our model adopts the approach suggested by Vichi et al. (2004) for nutrient limitation and bacterial biology, and Lazzari et al. (2012) for chlorophyll synthesis. Furthermore, as a large portion of the Adriatic Sea is characterized by shallow depths, we included a benthic return closure, able to provide an approximate balance between deposition of organic material and remineralization at the bottom. As BFM treats each numerical model grid cell as a separate biogeochemical reactor, it is one-way coupled with the transport module TOP that uses physical transport fields and turbulent diffusion coefficients from the output of an Ocean General Circulation Model-OGCM (Lovato et al., 2020) (Figure 1).

Figure 1

Figure 1. Diagram of the AdriaClim offline biogeochemical model architecture: BFM is coupled with the transport module of the NEMO modeling framework and fed with physical input provided by the AdriaClim modeling system and chemical input from Ludwig et al. (2009) and Butenschön et al. (2021).

The OGCM was implemented on a domain consisting of a regular grid covering the whole Adriatic Sea with a horizontal resolution of ~2 km and 120 vertical levels. The physical input (circulation data, surface water balance, river discharge, winds and shortwave radiation) was provided by the AdriaClim modeling system, developed to downscale the Med-CORDEX climate model LMDz-NEMOMED (L'Hévéder et al., 2013) at Adriatic scale, under scenario RCP8.5. The system includes an atmosphere model (WRF) with a resolution of about 6 km covering the Central Mediterranean area (Fedele et al., 2024), a hydrological model (WRFHydro) implemented over the same domain as WRF which includes the NOAH land surface submodel and a double shallow water system prognostically solving 145 catchments ending into the Adriatic Sea at about 600 m resolution, a marine circulation model (NEMO) and a wave model (Moulin et al., 2024) implemented on the same domain as the biogeochemical model i.e., the Adriatic basin with 2 km horizontal resolution.

The OGCM output useful for biogeochemistry consists of the horizontal (U, V) and vertical (W) components of the velocity, temperature (T), salinity (S), and vertical diffusivity (K) for the whole water column. Wind speed, short-wave radiation, and fresh-water balance (evaporation - precipitation) are provided at the surface. The river water discharge is provided as a surface boundary condition at the cells corresponding with the river mouths. The river chemical inputs consist of river loads of inorganic nutrients (nitrate, phosphate, silicate from Ludwig et al., 2009) and carbonate system parameters (Dissolved Inorganic Carbon, DIC, from Ludwig et al., 2009, and Total Alkalinity, TA, from Butenschön et al., 2021). As past studies found a strong linear correlation between runoff and nutrient transport in the Northern Adriatic (Cozzi and Giani, 2011; Ciglenečki et al., 2020) here we assumed a proportionality between river discharge and nutrient load. The atmospheric concentration of CO₂ was provided by the RCP8.5 specification (van Vuuren et al., 2011). The concentrations of oxygen, DIC, nitrate, phosphate, silicate, and TA at the lateral open boundary in the Ionian Sea (not shown) were obtained from Butenschön et al. (2021) and their exchanges are handled using the Orlanski Normal Propagation of Oblique radiation scheme (Orlanski, 1976; Marchesiello et al., 2001). A zero-gradient Neumann scheme was adopted for the open boundary conditions of the other biogeochemical state variables.

The model simulations comprise two different time windows: 1992–2020 for the baseline period, 2021–2050 for the future scenario. As marine biogeochemistry contains long-memory process, and the model takes long before reaching an equilibrium state, the model was run the first time with initial conditions estimated from the Copernicus Marine Service (CMS) reanalysis (Cossarini et al., 2021) for oxygen, DIC, nitrate, phosphate, silicate, and using constant values for the other constituents. Simulations were subsequently executed running the 10 years 1992–2011 recursively until reaching dynamic stability, and then using the state at the end of 2011 to re-initialize the model at the beginning 1992 for the Baseline Scenario Simulation (BSS). The state of BSS at the end of 2020 was then used to initialize the Future Scenario Simulation (FSS) for the period 2021–2050. To better characterize the ongoing tendencies, and the differences between the biogeochemical state of the Adriatic Sea at the beginning and the middle of the twenty-first century, we considered the baseline scenario simulation over the 20 years 1992–2011 (BSS_{1992 − 2011}), and the future scenario simulation over the 20 years 2031–2050 (FSS_{2031 − 2050}).

Classification of marine ecoregions

To identify the ecoregions of the Adriatic Basin, a classification was carried out using the K-Medoid methodology (Kaufman and Rousseeuw, 1987). Similar to K-Means, K-Medoid selects the cluster's centroid as the data point that has the lowest sum of distances to all other points in the cluster. The centroid is thus computed as the cluster's “Medoid,” rather than its barycenter as in K-Means, and bound to be an element of the analyzed sample. This comes with advantages in terms of robustness and flexibility, most noticeably, a reduced sensitivity to outliers, and the possibility of adopting customized, non-Euclidean metrics.

In this study, we designed a non-dimensional metric able to consider both the physical and chemical properties of the environment. Given 2 geographical points P₁ and P₂, the distance ξ is defined as Equation (1):

\begin{array}{l} ξ (P_{1}, P_{2}) = {[\int {(T_{N} (P_{2}, z) - T_{N} (P_{1}, z))}^{2} d z d t]}^{^{\frac{1}{2}}} \\ + {[\int {({[C h l]}_{N} (P_{2}, z) - {[C h l]}_{N} (P_{1}, z))}^{2} d z d t]}^{^{\frac{1}{2}}} \\ + {[\int {({[O 2]}_{N} (P_{2}, z) - {[O 2]}_{N} (P_{1}, z))}^{2} d z d t]}^{^{\frac{1}{2}}} \\ + κ_{d p t} \log {\frac{max [depth (P1), depth (P2)]}{min [depth (P1), depth (P2)]}} & (1) \end{array}

where T_N, [Chla]_N and [_O₂]N are the time-varying vertical profiles of temperature and of the concentrations of chlorophyll-a, and oxygen, normalized with the field mean and standard deviation, according to the formula X_N = [X−mean(X)]/σ_X. The last term is a “depth-penalty” introduced to increase the distance between points with different depths. Given 2 points P1 and P2, this term is equal to 0 when depth (P1) = depth (P2), while it increases logarithmically with the ratio between the depths. If a shallow point, e.g., with a depth of about 10 m, would be compared with an abyssal point, with a depth of about 1000 m, the term would be roughly equal to 5k_dpt. In general, we found good results by setting k_dpt to the low value of 0.02. The results of clustering techniques such as K-Means and K-Medoid depend on the algorithm seeds, i.e., on the initial set of first-guess centroids. In this study we selected 6 locations to diversify the geographic position, depth, and dynamic/biogeochemical regime (Figure 2A). The number 6 was chosen after several attempts, and allows for the generation of contiguous clusters, with a good interpretability compared with subdivisions existing in literature. The algorithm was applied to a monthly climatology of temperature, chlorophyll-a and oxygen computed from the baseline simulations BSS_{1992 − 2011}.

Figure 2

Figure 2. Classification of marine ecoregions based on the baseline scenario (A). Mean satellite and baseline scenario concentrations of chlorophyll-a, for winter (B, E), summer (C, F), and annual means (D, G). Sat, satellite observations.

Comparison with observations and reanalysis

To verify the ability of the model to provide realistic data, the output of BSS_{1992 − 2011} was compared, for each cluster, with observational climatologies for winter (January, February, March) and summer (July, August, September). These seasons were selected as the overall most and least productive ones, due to the combined effect of nutrient supply and mixing (e.g., Zavatarelli et al., 1998). To show the overall model's skill, we also compared simulation results and observations averaged on the year. The concentration of chlorophyll-a integrated over the water column down to the Secchi depth was compared with climatological means of the CMS color-based satellite observations for the Mediterranean Sea (Volpe et al., 2019). The climatological winter, summer and year profiles of chlorophyll-a, nitrate, phosphate and oxygen were compared with those provided by the EMODnet chemistry climatology from 1969 to today (Giorgetti et al., 2020) and by the CMS reanalysis from 1999 to present (Cossarini et al., 2021).

Assessment of future changes

Future changes were assessed by comparing the concentrations of chlorophyll-a, oxygen, phosphate, nitrate, Net Primary Productivity (NPP), phytoplankton biomass and the values of temperature and pH on the baseline period (1992–2011) vs. the future period (2031–2050). The values were averaged for each cluster, over the winter, summer, and the whole year.

Results

Ecoregions of the Adriatic Sea

The clustering technique identifies regions of the Adriatic Sea characterized by different biogeochemical dynamics, which were labeled based on their geographic location and trophic regime (Figure 2A). The regions are:

• Northern Estuarine (NE). This area is defined by the estuaries of the Northern Adriatic river catchments, in particular the Po Delta, and by high loads of nutrients. It stands out as the most nutrient-rich region in the Adriatic Sea. The prevalently cyclonic circulation of the Northern Adriatic and its Western Adriatic Coastal Current intensification carries the constituents southward, justifying the extension of the first cluster along the central Italian coasts. This area overlaps with the Shallow Northern Adriatic area identified by Zavatarelli et al. (1998).

• Italy Center-South (ICS). This region comprises the coastal areas of center and southern Italy, and is defined by the western middle and southern Adriatic currents (W-MAd and W-SAd) that transport southward temperature and constituents, and by lower trophic levels with respect to NE.

• Eastern Coastal (EC). Coastal area of Southern Croatia, Montenegro and Albania. Characterized by the deltas of some important rivers, among the others the river system of Buna-Bojana and Drin, Vjosa, and Neretva. Its circulation is also influenced by the Levantine Intermediate Water (LIW) entering the Adriatic Sea through the Otranto Channel.

• Northern Mesotrophic (NM). Comprises a large portion of the offshore Northern Adriatic and the northern portion of the Croatia coasts, and shows more oligotrophic conditions with respect to NE, though still influenced by the hydrological basins of the Northern Adriatic. It partially coincides with the Deep Northern Adriatic area identified by Zavatarelli et al. (1998).

• Intermediate-Depth Center (IDC). Area in the central Adriatic with a median depth above 100 m and a maximum depth of about 270 m, which includes the Middle Adriatic Pit, and with prevalently oligotrophic conditions. This cluster stretches southward parallel to the ICS region, in response to the prevalent W-MAd currents.

• Southern Deep (SD). This is the southern abyssal portion of the Adriatic, with a median depth of about 800 m, characterized by the Southern Adriatic gyre and affected by the LIW entering the Adriatic Sea through the Otranto Channel.

Skill vs. historical data

The model reproduces spatial and temporal patterns of chlorophyll-a qualitatively similar to those provided by color-based satellite observations, both on an annual and seasonal basis (Figures 2B–G). The north-western Adriatic basin stands out as highly productive in model results and satellite observations, especially during winter, and stretches southward along the Italian coast, consistent with the prevalent currents. Furthermore, our model also captures relatively high chlorophyll-a concentrations along the eastern coasts in proximity to important rivers' mouths, especially during winter, though it shows limited ability to reproduce the nutrient-rich coastal area at the Albanian coasts. While our model reproduces annual chlorophyll-a concentrations that are lower than those from satellite data in nutrient-rich areas, chlorophyll-a concentrations in oligotrophic areas are higher than those from CMS satellite imagery during Winter.

The fit with observations noticeably improves when we compare the model's vertical profiles of chlorophyll-a with the EMODnet climatology for each marine ecoregion identified in the classification, especially during winter (Figures 3, 4). In the NE region, the bias of median concentration of chlorophyll-a vs. EMODnet is 0.09 μg L⁻¹, or little more than 10%, while the Root Mean Squared Difference (RMSD) is 0.1, <20% of the median values (Figure 3A). During winter the largest differences with the observations are in Eastern Coastal (Figure 3R) and in the Intermediate-Depth Center (Figure 4I), where model output shows a higher concentration of chlorophyll than EMODnet, and in the Northern Mesotrophic area, where the model reproduces conditions significantly more oligotrophic than EMODnet (Figures 4A–H). Notably, during winter, the AdriaClim model reproduces generally higher concentrations of chlorophyll-a than the CMS reanalysis and is usually more in line with the EMODnet estimations. On the other hand, the satellite data shows concentrations of chlorophyll-a substantially higher than those from EMODnet in NE and ICS ecoregions during winter. During summer the model shows the presence of Deep Chlorophyll Maxima (DCM) more pronounced than both EMODnet and CMS reanalysis, with correspondingly lower concentrations of chlorophyll-a in the proximity of the surface (Figures 3E, N, V, 4E, N, V). The depth of the DCM is generally well reproduced (e.g., Figures 3V, 4N, V).

Figure 3

Figure 3. Vertical profiles for winter (A-D, I-M, R-U) and summer (E-H, N-Q, V-Y) for the concentrations of chlorophyll-a (A, E, I, N, R, V), oxygen (B, F, K, O, S, W), phosphate (C, G, L, P, T, X), nitrate (D, H, M, Q, U, Y), for the ecoregions Northern Estuarine (NE, A–H), Italy Center-South (ICS, L–Q), Eastern Coastal (EC, R–Y). The thick blue line is the model median, the cyan area is the model interquartile range, the red dots/lines are the EMODnet climatology median/interquartile range, the hatched areas are the profiles from the CMS reanalysis, the green dot/line is the CMS satellite observations median/interquartile range at the surface. The model vs. EMODnet bias and Root Mean Squared Difference (RMSD) are also reported.

Figure 4

Figure 4. As in Figure 3, for the ecoregions Northern Mesotrophic (NM, A–H), Intermediate-Depth Center (IDC, I–Q), and Southern Deep (SD, R–Y).

The vertical profiles of dissolved oxygen concentration reproduced by the model show a very good agreement with those of the EMODnet climatology and the CMS reanalysis, in all ecoregions (Figures 3C, G, L, P, T, X, 4C, G, L, P, T, X). During winter/summer, the bias of the median profile in the NE vs. EMODnet is −10.81/11.13 μM equivalent to 3.94/3.70% of the predicted value averaged over the water column. In the notation, the symbol μM represents a micro-molar concentration (IUPAC, 2019). The RMSD is 11.13/8.15 μM, equivalent to 4.05/3.79% of the predicted value averaged on the water column. The other ecoregions all exhibit lower biases, and values of RMSD are lower in winter than in summer, except for NM, where the winter/summer RMSD is 5.33/3.48 μM.

The concentrations of inorganic nutrients such as phosphate and nitrate exhibit higher uncertainty, with substantial differences between BSS_{1992 − 2011} and EMODnet, BSS_{1992 − 2011} and the CMS reanalysis, and the EMODnet and the CMS reanalysis. In NE, the concentration of phosphate shows vertical profiles significantly higher than both EMODnet and the CMS reanalysis (Figures 3C, G). In all the other marine ecoregions the BSS_{1992 − 2011} concentrations of phosphate are similar to those from the CMS reanalysis, but substantially lower than those from the EMODnet climatology (Figures 3L, P, T, X, 4C, G, L, P, T, X). For nitrate the picture is more diversified: in some areas/seasons our model's concentration profiles are closer to those from EMODnet (ICS summer, NM winter, Figures 3M, 4D), in others to those from the CMS reanalysis (NE summer and winter, NM summer, SD summer and winter, Figures 3D, H, 4H, U, Y), in others EMODnet and the CMS reanalysis are closer to each other than in BSS_{1992 − 2011} (ICS winter, EC and IDC winter and summer, Figures 3Q, U, Y, 4M, Q). Similar to phosphate and nitrate, the vertical profiles of silicate ${SiO}_{4}^{4 -}$ concentration exhibit significant differences with respect to the ones from EMODnet (Supplementary Figures S1, S2). In particular, the concentration of silicate from our model is higher than EMODnet in NE (Supplementary Figures S1A, D) and lower in the deeper layers of SD (Supplementary Figures S2C, F), while the ecoregions where our vertical profiles most closely align with EMODnet are ICS (Supplementary Figures S1B, E) and NM (Supplementary Figures S2A, D).

Projected changes

According to our projections, the Adriatic Sea will experience strong oligotrophication along with warming, deoxygenation, and acidification in the next decades. This trend is observable for all the EOVs considered in this work within the identified marine ecoregions (Figure 5). The temperature is projected to rise by <1°C between 1992–2011 and 2031–2050 within all ecoregions, with distinct differences between winter and summer, especially in the shallow areas of the Northern Adriatic. The NE ecoregion displays the smallest mean annual increase in winter and summer temperatures (0.28°C and 1.33°C), the NM ecoregion the largest increase (0.76 and 1.41°C, respectively), with the SD reporting smaller differences with a mean increase corresponding to 0.91°C in both seasons (Figures 5A–C).

Figure 5

Figure 5. Projected changes of temperature (A-C), oxygen (D-F), pH (G-I), Net Primary Productivity (NPP, J–L), chlorophyll-a (M–O), phytoplankton biomass (P–R), inorganic phosphate (S-U) and nitrate (V-X), in each ecoregion of the Adriatic Sea: Northern Estuarine (NE), Italy Center South (ICS), Eastern Coastal (EC), Northern Mesotrophic (NM), Intermediate-Depth Center (IDC), Southern Deep (SD), for winter (A, D, G, J, M, P, S, V), summer (B, E, H, K, N, Q, T, W), and on the year (C, F, I, L, O, R, U, X). At the base of the histogram is reported the projected change for each variable/area, in absolute terms for T, pH and oxygen, as a percentage for the remaining variables. The symbol μM C, used for NPP and phytoplankton biomass, indicates a micromolar concentration of carbon equivalents.

The oxygen concentration is also projected to consistently decline within all marine ecoregions (Figures 5D–F). On an annual basis, the ICS and the NE are the ecoregions displaying the largest and smallest decrease in oxygen concentrations (5.46 μM and 2.15 μM year on year, respectively), and those with the maximum and minimum increase in temperature. The pH is also projected to decrease in all marine ecoregions, with relatively small seasonal variability and annual differences future vs baseline >0.07 in EC, IDC and SD ecoregions, with the NE showing the lowest difference (0.03, Figures 5G–I).

According to our model, the nutrient-rich NE and the ICS and NM marine ecoregion, which are heavily influenced by the regimes of the Northern Adriatic hydrological catchment, are projected to experience the greatest decline in biological activity-related variables, such as NPP, chlorophyll-a and phytoplankton carbon. In these three ecoregions, chlorophyll-a levels will drop on average by 49.4%, 14.2%, and 39.6% respectively, while phytoplankton carbon levels will decrease by 43.5%, 24.6%, and 40.7% year on year, with respect to the baseline period. The projected change is less pronounced for the other ecoregions (EC, IDC, SD), with an average decline in chlorophyll-a lower than 10%, and in phytoplankton C lower than 15% on the year. Furthermore, these ecoregions exhibit a rise in NPP (Figures 5J–R).

The concentrations of nutrients such as inorganic phosphate and nitrate are projected to decrease substantially across all marine ecoregions in the next decades, especially during winter. When the annual average is considered, the decrease in phosphate concentration will be >35% in all ecoregions. In NE, NM and IDC, the phosphate concentration is projected to decrease by 48.6%, 66.3% and 56.1% respectively. In the case of nitrate, the decrease in concentrations is projected to reach 34.6% in NE, 62.8% in NM and 40.2% in IDC.

Discussion

Our classification of the marine ecoregions based on the physical and biogeochemical outputs of the BSS_{1992 − 2011} shows a realistic subdivision of the Adriatic Sea. The NE ecoregion closely resembles the Shallow Northern Adriatic (SNA) area identified in the pioneering work by Zavatarelli et al. (1998), who also discussed the dominant role of river nutrient input in that area. The NM ecoregion partially coincides with a superimposition of the Deep Northen Adriatic (DNA) and the northern part of the Eastern Coastal Region described in Zavatarelli et al. (1998). The other areas differ significantly from the more traditional subdivision of the basin into “stripes” in a North-West to South-East direction adopted in previous studies on physics (e.g., Artegiani et al., 1997; Vilibić and Supić, 2005; Mantziafou and Lascaratos, 2008; Janeković et al., 2014) or biogeochemistry (Gacic et al., 1999; Polimene et al., 2007; the same Zavatarelli et al., 1998 for other areas), often due to limitations in spatial and temporal coverage related to field-collected observations. Except for the NM cluster, which combines mixing mesotrophic coastal and offshore features, the classification provided in this study results in a separation between sub-areas with predominant coastal regimes (in clusters NE, ICS, EC) and areas located offshore (in clusters IDC and SD).

This technique offers a different perspective on how the Adriatic Sea can be subdivided, in that it can identify marine ecoregions consistent from the point of view of their physical and trophic properties. The model capacity to capture the primary biogeochemical features/variables of the Adriatic ecosystems on the baseline scenario is demonstrated by the soundness of the subdivision and the overall correspondence of the vertical profiles with observations and reanalysis in each of the marine ecoregions. In particular, the model reproduces nutrient-rich conditions in the northern portion of the basin, and oligo-mesotrophic median conditions elsewhere, and correctly detects winter as the most productive season, in general agreement with the literature (e.g., Zavatarelli et al., 1998; Polimene et al., 2007; Fiori et al., 2016).

The vertical profiles of EOVs such as chlorophyll-a, oxygen, inorganic phosphate, and nitrate show different degrees of similarity with those from the EMODnet climatology and the CMS reanalysis. For oxygen the resemblance is remarkable in all the ecoregions and in both seasons, especially when compared with EMODnet, showing that the BFM model is capable of accurately capturing the dynamics of this important parameter (e.g., Di Biagio et al., 2022).

In comparison to the EMODnet and the CMS reanalysis, the AdriaClim model shows a systematic tendency to produce a more pronounced DCM during summer. DCMs tend to develop in the presence of low vertical mixing, otherwise an intense mixing would smooth the vertical gradients (e.g., Huisman et al., 2006; Weiss and Provenzale, 2008; Liccardo et al., 2013). In shallow areas such as the NE and the NM, the difference between AdriaClim and CMS can be explained by the lower velocity-sheer-dominated vertical eddy diffusivity, which comes with a reduced mixing at the surface (Supplementary Figures S2A, D), together with the fact that unlike CMS, our model implements a benthic return closure that withholds the biogeochemical constituents in the bottom layers of the water column and releases them back to the pelagic system. The vertical eddy diffusivity in the ICS is lower than CMS both at the surface (Supplementary Figure S2B), and below 40 m depth, where the chlorophyll-a concentration is at its maximum (Supplementary Figure S3B). The marked DCM in the SD ecoregion can be associated with vertical eddy diffusivity that is significantly lower than in CMS below a depth of 30 m (Supplementary Figure S3F). In the other 2 areas (EC and IDC) the difference between AdriaClim and CMS seem less connected to the vertical mixing, as in both cases the CMS vertical eddy diffusivity is higher at the surface (Supplementary Figures S2C, E) but not in the deeper layers, where the profiles are both characterized by very low mixing between 30 and 80 m depth (Supplementary Figures S3C, E). A hypothesis is that in these 2 areas the horizontal circulation from neighboring areas may play an important role.

The vertical profiles of chlorophyll-a reproduced by BSS_{1992 − 2011} during the winter months often resemble those from EMODnet more closely, with concentrations higher than the CMS reanalysis everywhere. This can be related to the more accurate representation of the rivers in our model with respect to the CMS reanalysis. Indeed in our modeling system a subset of 71 rivers among the 145 solved by WRFHydro model enter the marine components NEMO and BFM, following criteria described in Verri et al. (2024). On the other hand, the CMS reanalysis considers only 15 rivers flowing into the Adriatic basin with no rivers on the western coastline south of Po River mouths. The difference between the surface concentration of chlorophyll-a determined from satellites, EMODnet, and the models in the NE and ICS is noteworthy, which suggests a possible overestimation of color-based satellite estimates in estuarine areas characterized by high sediment load. Noteworthy, Mozetič et al. (2010) and Kovač et al. (2018) reported similar results when comparing long-term time-series data of in-situ measured chlorophyll with satellited-derived time series in the North Adriatic Sea.

The discrepancies between observed and modeled concentrations of phosphate and nitrate have many fold causes. One of these is sparsity and uncertainty in the observations, which hamper a thorough assessment of the modeled profiles (Guyennon et al., 2015). Despite the development of Biogeochemical Argo floats in recent years could change the scenario for large-scale models (Teruzzi et al., 2014; Salon et al., 2019), these drifters are not frequently available in a shallow marginal basins such as the Adriatic Sea. Additionally, according to Zhu and Ma (2020) using currently available analytical techniques to quantify phosphate in aquatic environments still presents several challenges. These include issues with precision in trace-level determination, compromised analytical functionality in challenging environmental circumstances, and the instability of in situ ${PO}_{4}^{3 -}$ techniques over extended periods of deployment.

Finally, the fact that these datasets are examined over different periods (1992–2011 for BSS_{1992 − 2011}, 1978–2020 for EMODnet and 1999–2011 for CMS) represents a limitation in the comparison vs. EMODnet and the CMS reanalysis. On one hand, our model's goal is to reproduce a reliable historical climate baseline, which justifies the comparison with climatological fields from different time frames. On the other hand, the possible presence of underlying trends could contribute to amplifying the differences between CMS and EMODnet.

Oligotrophication trend

The projected tendency confirms trends that have been observed in historical data on warming (e.g., Shaltout and Omstedt, 2014; Mohamed et al., 2019; Bonacci et al., 2021), oligotrophication (e.g., Mozetič et al., 2010; Coppini et al., 2013; Cerino et al., 2019; Novak et al., 2019; Ciglenečki et al., 2020; Penezić et al., 2022), acidification (Luchetta et al., 2010; Hassoun et al., 2015, 2022; Giani et al., 2023), and corroborates the results of recent climate projections at Mediterranean scale (Reale et al., 2022). Our study provides a closer look at the diversification of the projected changes within the Adriatic basin.

The oligotrophication trend is particularly strong in the northern portion of the Adriatic basin, especially in NE, where the annual NPP and the concentrations of phytoplankton biomass and chlorophyll-a are projected to decrease respectively by 37.8%, 46.3%, and 54.5% by the 20 years 2030–2050. This sharp change is mostly due to a decrease in discharge of the Northern Adriatic rivers projected by the hydrological model (Verri et al., 2024), and the related reduction in nutrient load. Furthermore, we rely on the hypothesis, consistent with scenario RCP8.5, that no policy will be enforced to counteract this tendency. This decrease amounts to about 35% on the year for the rivers with discharge location in the Northern Estuarine area, rising to about 40% if the sole Po Plain is considered (Figure 6, histograms NE and NE_PV). In this regard it is worth mentioning that an early signal of a regime shift toward oligotrophication was revealed at the beginning of the 1980s by long time-series of water column nutrients, chlorophyll a, and phytoplankton biomass (Marić et al., 2012; Marini and Grilli, 2023), and became more evident at the turn of the 2000s in the Gulf of Trieste (Cozzi et al., 2020) and in the North Eastern Adriatic (Djakovac et al., 2012), evidencing main the role of the decline in Po River discharge (Cozzi et al., 2018) in the North Adriatic basin, and that of minor rivers in the coastal areas of the Central Adriatic (Ricci et al., 2022).

Figure 6

Figure 6. Projected changes of river discharge, in the ecoregions located along the coasts: Northern Estuarine (NE), Italy Center South (ICS), Eastern Coastal (EC), Northern Mesotrophic (NM), and in the NE region restricted to the Po Plain (NE_PV), for winter (A), summer (B), and on the year (C). The projected percentage change of discharge is reported at the base of the histogram.

The NM ecoregion, strongly influenced by the dynamics of the northern estuaries via horizontal transport and mesoscale filaments, also shows a significant decrease in annual biological activity: −12.1% in NPP, −33.2% in chlorophyll-a, −33% in phytoplankton biomass year on year. These changes appear aligned with the recent occurrence of extremely low annual discharges (e.g., Gervasio et al., 2023; Montanari et al., 2023; Aragão et al., 2024), and the low levels of nutrients and chlorophyll-a concentrations measured in the water column of the NE ecoregion (e.g., ARPAE, 2021; Gašparović et al., 2023). Although it is still too soon to interpret these data as harbingers of future mean changes, they surely represent an unprecedented event in our records (Bonaldo et al., 2023; Montanari et al., 2023). A climate study conducted at high resolution on the Po catchment reports a decrease in discharge comparable to what is presented in Figure 6 by the end of the century under scenario RCP8.5 (Vezzoli et al., 2015). On the other hand, our results are only in partial agreement with past climate studies based on ensembles of lower-resolution CORDEX models, that point toward a significant increase in the intensity and frequency of droughts in the Po Plain, at 1.5° and 2° of global warming levels, albeit without a significant change in mean discharge, and even an increase in the high extremes (Cammalleri et al., 2020; Mentaschi et al., 2020; Naumann et al., 2021; Zittis et al., 2021; Dottori et al., 2023). These studies reveal a substantial variability in the projected changes of hydrological regimes within the Adriatic Sea catchments, dependent on the specific model and scenario employed. At any rate, the findings of this study offer insights into the potential repercussions of prolonged decreases in the discharge of Northern Adriatic Rivers. This scenario gains more significance in light of the recent unprecedented low flow observed in the rivers of the Po Plain. If the trends reproduced by our model will be confirmed by future observations, these changes, together with warming and acidification, can potentially threaten the resilience of the nutrient-rich ecosystems of the Northern Adriatic, and the sustainability of related human activities.

The growth in NPP in all the remaining areas mirrors an increase in plankton metabolic rates related to the rise in temperature. This outcome, which corroborates a result by Reale et al. (2022) for the Mediterranean Sea, does not contradict past climate studies that found a decrease in NPP at a global scale related to global warming (e.g., Fu et al., 2016; Tagliabue et al., 2021). In general, such reduction is associated with an increase in ocean stratification at a global scale, while according to our model, in the Adriatic Sea the stability of the water column is projected to decrease in response to the reduction in river discharge and the associated increase in surface salinity (Verri et al., 2024). The growth in NPP occurs together with a generalized decrease in chlorophyll-a and in phytoplankton biomass. This could be due to the fact that amid warming conditions the heterotrophic respiration rates tend to rise faster than photosynthetic efficiency (e.g., O'Connor et al., 2009; Boscolo-Galazzo et al., 2018; Kralj et al., 2019). A temperature-related intensification of the biogeochemical cycles could also justify the sharp decrease in inorganic nutrients, especially phosphate, together with the reduction in river discharge and nutrient load (Figure 6), along with a decrease in nitrate and phosphate entering through the open boundary condition that will mainly affect the southern areas. Interestingly, the decrease in river discharge is especially strong in the EC ecoregion but with less dramatic effects (Figure 6), than in the NE, as the EC is a less estuarine and deeper area, which is influenced by the abyssal dynamics of SD, as well as by the southern boundary conditions. Overall, marine ecoregions except the NE and the NM display a noteworthy decline in phytoplankton biomass, phosphate and nitrate concentrations, and a lower reduction in chlorophyll-a, thus confirming a predominant oligotrophication trend in the future.

The projected changes in mean oxygen concentration, all negative, are mostly explained by the warming-related reduction in solubility (e.g., Matear and Hirst, 2003; Schmidtko et al., 2017). Another possible cause mentioned in the literature, such as the increased ocean stratification (Helm et al., 2011), may indirectly affect the southern part of the Adriatic Sea through the open boundary (Verri et al., 2024). The acidification trend obtained within the different regions of the Adriatic Sea is of the same order of changes detected at the global scale (e.g., Doney et al., 2009) as well as in the Mediterranean Sea (Cossarini et al., 2021; Hassoun et al., 2022; Reale et al., 2022; Solidoro et al., 2022; Kristiansen et al., 2024). The Northern Estuarine area stands out as more acidic than the others at present day, and at the same time less affected by acidification, hinting at the important role that rivers can play as sources of carbon and alkalinity, already observed in past studies (Müller et al., 2016; Butenschön et al., 2021; Gomez et al., 2021; Reale et al., 2022), and object of further research in future projects on the Adriatic Sea. Indeed, in the literature this behavior is associated with the fact that in the Northern Adriatic Sea the TA is in the upper range of that measured in the Mediterranean Sea, primarily because riverine inputs transport carbonates dissolved from Alpine dolomites and karstic watersheds (Brush et al., 2020). Acidification poses a serious threat to marine life, as it increases the undersaturation of aragonite. Ocean acidification has adverse effects on marine ecosystems, and compromises their capacity to provide functions and services, the provision of ecosystem services including mariculture and fisheries, with important consequences on local economies within the Mediterranean Sea (Partelow et al., 2023). Multiple lines of evidence, ranging from simulated experiments and laboratory studies have demonstrated the negative effects of global warming and ocean acidification on phytoplankton community structure (Dutkiewicz et al., 2015), fish eco-physiology and behavioral performances (Cattano et al., 2018), as well as on calcification, larval production and midstage growth in farmed bivalves (Barton et al., 2012; Gazeau et al., 2014), which play an important role in Northern Adriatic's ecology and economy (Lacoue-Labarthe et al., 2016).

Final remarks

The Adriatic Sea is a unique basin due to its semi-enclosed geomorphology, variable depth, and diversified trophic state, with nutrient-rich conditions prevailing along the north-western coasts, and oligotrophic conditions elsewhere. According to our projections, the basin will be affected by oligotrophication, acidification and deoxygenation in the next decades, confirming findings from earlier studies, both on historical and climate data. Our high-resolution model reveals that projected trends across the Adriatic Sea are diverse and related to the various trophic conditions in its sub-regions as well as to many underlying causes. For instance, alterations in hydrological regimes are particularly relevant in the shallow Northern Adriatic, whereas shifts in stratification will have a more pronounced impact on mesotrophic and oligotrophic regions of the central-southern Adriatic. In this respect, the classification of the marine ecoregion employed in this study served as a useful tool to identify local dynamics and influencing factors.

The oligotrophication is projected to be particularly intense in the shallow regions of the North Adriatic basin by 2031–2050. This result comes with the clear limitation that it was obtained downscaling a single Med-CORDEX model, and is associated with low river discharge, which was most severe in the Po Plain. If confirmed by future research on climate models and by an evident trend in current data, this tendency, along with the ongoing acidification and warming, casts doubt on the resilience and long-term sustainability of the nutrient-rich marine ecosystem of the Adriatic Sea. Oligotrophy is often regarded as a healthy state, and even a desirable outcome of policies devised to counter past conditions of nutrient pollution and eutrophication (e.g., Justić, 1987; Tsiamis et al., 2013; Fiori et al., 2016; Giovanardi et al., 2018; Le Fur et al., 2019). On the other hand, excessive oligotrophication can come with unproductive and declining biological resources, especially in ecosystems long evolved to be supported by high trophic conditions (Stockner et al., 2000). In this sense, the findings of this study provide insight into how the projected decline in the hydrological regimes of the Northern Adriatic could negatively affect the future of the marine food web structure. The ongoing and projected oligotrophication of the Adriatic Sea warrants deeper investigation, in view of the potential harm it poses to important sectors of the economy such as fisheries. This aligns with the concern, recently expressed by the European Commission, on the resilience and sustainability of healthy ecosystems, and with the objectives stated by the Sustainable Development Goal (SDG) 14 of the United Nations, “life below water”.

Data availability statement

The raw data supporting the conclusions of this article are available at the CMCC ERDDAP portal, https://erddap.cmcc-opa.eu/erddap/index.html.

Author contributions

LM: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Supervision, Validation, Visualization, Writing—original draft. TL: Investigation, Software, Writing—review & editing. MB: Investigation, Writing—review & editing. JA: Writing—review & editing. LA: Writing—review & editing. GV: Writing—review & editing. RG: Investigation, Writing—review & editing. GC: Funding acquisition, Writing—review & editing. NP: Supervision, Writing—review & editing.

Funding

This research was funded by the Interreg IT-HR AdriaClim, project id 10252001. MB received funding from the ACTNOW project within the European Union's Horizon Europe Research and Innovation Programme under grant agreement No. 101060072.

Acknowledgments

The authors acknowledge prof. Marcello Vichi for his support and suggestions.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fclim.2024.1338374/full#supplementary-material

References

Alvisi, F., and Cozzi, S. (2016). Seasonal dynamics and long-term trend of hypoxia in the coastal zone of Emilia Romagna (NW Adriatic Sea, Italy). Sci. Total Environ. 541, 1448–1462. doi: 10.1016/j.scitotenv.2015.10.011

PubMed Abstract | Crossref Full Text | Google Scholar

Aragão, L., Mentaschi, L., Pinardi, N., Verri, G., Senatore, A., Sabatino, D., et al. (2024). The freshwater discharge into the Adriatic Sea revisited. Front. Clim. (under review).

Google Scholar

ARPAE (2021). Qualità Ambientale Delle Acque Marine in Emilia-Romagna. Rapporto Annuale 2021. Washington, DC: ARPAE.

Google Scholar

Artegiani, A., Paschini, E., Russo, A., Bregant, D., Raicich, F., Pinardi, N., et al. (1997). The Adriatic Sea general circulation. part i: air–sea interactions and water mass structure. J. Phys. Oceanogr. 27, 1492–1514. doi: 10.1175/1520-0485(1997)027<1492:TASGCP>2.0.CO;2