The simulations have been performed by using a state-of-the-art transport-biogeochemical model, purposely developed and validated for describing plankton productivity and carbon and nutrients biogeochemical cycles in the MS (Lazzari et al., 2012, 2016). The model is the product of a long-term research activity on MS ecosystem modeling (Crise et al., 1998; Crispi et al., 2001; Solidoro et al., 2003; Lazzari et al., 2012), it is now embedded in the European Copernicus Marine Environment Monitoring Service (CMEMS) infrastructure, and routinely used to produce short-term forecasts of the MS biogeochemistry (Salon et al., 2019; Cossarini et al., 2021). The capability of the model to reproduce (without any data assimilation) the main features of MS biogeochemistry has been assessed by contrasting model results against satellite estimates of surface chlorophyll (Lazzari et al., 2012), available information on primary productivity (Lazzari et al., 2012), nutrients (Lazzari et al., 2016), dissolved inorganic carbon and alkalinity data (Cossarini et al., 2015). It has been used in a number of applications, both at the basin wide scale (e.g., Lazzari et al., 2014; Melaku Canu et al., 2015; Galli et al., 2017) and for smaller subregions (Lamon et al., 2014; Fraysse et al., 2013).

The physical dynamics were pre-computed by the ocean general circulation model NEMOMED8 at 1/8° horizontal resolution (Sevault et al., 2009), then used in Beuvier et al. (2010), Herrmann et al. (2010), Albouy et al. (2013), Adloff et al. (2015), and Richon et al. (2019) for climate scenarios studies. NEMOMED8 supplied the temporal evolution at daily frequency of the 3D fields of horizontal and vertical current velocities, vertical eddy diffusivity, temperature, salinity, in addition to surface data for solar shortwave irradiance and wind stress. The meshgrid has 1/8° × 1/8°cos(latitude) horizontal resolution ranging from 9 to 12 km from the North to the South of the basin with square meshes and with a stretched and tilted grid at the Gibraltar Strait that increase locally the resolution up to 6 km (160 × 394 grid points in latitude and longitude directions), and 43 vertical z-levels (non-terrain following).

Details on the NEMOMED8 model and the protocol of the simulation scenario A2 can be found in Albouy et al. (2013, see their supplementary material) and Adloff et al. (2015). The same configuration for the SRES A2 scenario was adopted by Richon et al. (2019).

The biogeochemical dynamics are simulated using the Mediterranean implementation of the Biogeochemical Flux Model (BFM; Lazzari et al., 2012), where chlorophyll and carbon dynamics are based on the parameterization of chlorophyll synthesis proposed by Geider et al. (1997), and on a Plankton Functional Type (PTF, Anderson, 2005) representation of the planktonic food web that includes nine plankton pools, subdivided in photosynthetic producers (phytoplankton), consumers (zooplankton), and decomposers (bacteria). These broad functional groups are further subdivided into functional subgroups, so to have four phytoplankton groups (representative of diatoms and of not-diatoms of different size), heterotrophic bacteria, two microzooplankton groups preying on bacteria and smaller phytoplankton groups, and two mesozooplankton groups feeding of microzoo and larger phytoplankton groups. Each PFT is described in terms of its internal concentration of carbon, nitrogen, phosphorus and—for diatoms—silica, and contribute to the interlaced biogeochemical cycles of these elements. The model also includes the descriptions of the labile, semi-refractory, and refractory dissolved organic matter and of the particulate matter, as well as the dynamics of the microbial loop, which is known to play an important role in MS (Thingstad and Rassoulzadegan, 1995).

The carbonate system consists of two prognostic state variables: alkalinity (ALK) and dissolved inorganic carbon (DIC), from which pH is also computed (Cossarini et al., 2015). DIC exchange at the air-sea is resolved by computing the seawater pH, pCO₂ and gas transfer formula (OCMIP I model, Orr, 1999; Melaku Canu et al., 2015).

The physical and biogeochemical models are coupled through a modified version of the OPA transport model (Foujols et al., 2000), named OGSTM-BFM (OGS Transport Model—Biogeochemical Flux Model; Lazzari et al., 2012), which was already used in many studies on MS biogeochemistry (Lazzari et al., 2012, 2014, 2016; Cossarini et al., 2015; Melaku Canu et al., 2015) and it is currently employed as a core component of the CMEMS forecast system for biogeochemical properties of the MS (Salon et al., 2019), integrated with a data assimilation scheme (Teruzzi et al., 2014, 2018, 2019).

The scenario simulation carried out covers the period between January 2000 to December 2099 over the MS domain (Figure 1), and assuming the atmosphere forcing are in agreement with a high CO₂ scenario, representing a future socio-economic system in which economic activities and industrial processes continue to grow in line with the current situation. In particular, we adopted a scenario derived from the A2 emission scenario, therefore at the higher end of the SRES family related to a story line characterized by heterogeneity in the socio-economic system and a continuous population growth that reaches over 10 billion by 2050. According to Integrated Assessment Models this scenario has an expected CO₂ concentration of about 575 and 870 ppm for the middle and end of the 21st century, respectively. Therefore, it is close and—for all of our purposes equivalent—to the RCP8.5 scenario, which assumes similar CO₂ concentration regardless of the underpinning socio-economic dynamics (IPCC, 2014; Harrison et al., 2019). The time course of CO₂ atmospheric concentration used in the simulation is illustrated in the plot within Figure 2.

We chose to focus on this scenario because it might trigger signals large enough to be detected. At the same time, it is not unrealistic and, actually, analysis of historical data over the last decades indicates that it corresponds to the current trajectory (2005 to present) of emissions (Schwalm et al., 2020).

The boundary condition at Gibraltar and the biogeochemical initial conditions of the control run were set up in agreement to Lazzari et al. (2012), and therefore based on the MEDA-MEDATLAS database, while nutrient inputs from major rivers and Black Sea were based on Ludwig et al. (2009). Regarding the carbonate system variables (DIC and alkalinity), the initial conditions and inputs from rivers and Black Sea are as described in Cossarini et al. (2015), while the boundary condition at the Gibraltar were based on the CMIP5 global models. In the scenario simulation, rivers inputs were kept constant as in present-day condition, the Gibraltar boundary condition evolve consistently with the CMIP five global model output. The control run spin-up period was 50 years.

While the coupled model considers about 50 state variables and an even larger number of processes, it might be interesting to synthetize the major processes driving the carbon cycle in a conceptual scheme, useful also to highlight the relationships between the nutrients (N,P,Si) and carbon cycling and environmental forcing. The scheme is illustrated in Figure 3.

Atmospheric CO₂ dissolves in sea water through the air-sea exchange process (arrow #1 in figure), proportional to the difference between concentration of CO₂ in the atmosphere and in the surface ocean, and occurring at a rate depending on sea temperature T and other factors affecting the interface, such as wind intensity. Since the major driver of the dissolution/degassing process is the CO₂ gradient between atmosphere and surface sea, all factors able to change the CO₂ sea surface concentration indirectly affect this process. The factors acting to reduce sea surface CO₂ will favor CO₂ atmospheric dissolution, those acting to increase it will favor CO₂ degassing. Most of these factors depend on sea temperature, in turn a function of climate warming and ocean transport, but many are also influenced by biogeochemical processes.

Once dissolved in the sea, CO₂ chemically reacts with water to form carbonic acid (H₂CO₃), which partially dissociates in bicarbonate (HCO3−) and bicarbonate (CO3=) ions, while releasing hydrogen ions H⁺. The concentration of H⁺ defines the ocean acidity, measured by pH. Other substances present in the ocean that partially dissociate in hydrogen ions and weak chemical bases impacts on the carbonate system equilibria by bonding to (and therefore removing) some of the hydrogen ions. The final equilibrium concentrations of the carbon dioxide system (carbonate system) and of hydrogen ions (pH) therefore result from multiple equilibria, and depend on the total amount of dissolved inorganic carbon (DIC, the sum of CO₂ hco3 and co3), the total alkalinity (excess of proton acceptors over proton donors), and temperature (Zeebe and Wolf-Gladrow, 2001; Middelburg et al., 2020).

The photosynthesis (arrow #2 in Figure 3) reduces DIC concentrations by transforming DIC in living particulate organic matter. Beside temperature, photosynthesis rate depends on light availability, in turn a function of water transparency and surface irradiance level. The opposite process is respiration (arrow #3 in figure), a process continuously performed by all living organisms and during which organic carbon is burned and transformed back to DIC. During photosynthesis, the photosynthetic organisms also need to assimilate nitrogen and phosphorus, that are uptaken from the sea. Therefore, a shortage of dissolved inorganic nutrients causes a limitation to photosynthetic activity, and hence in carbon cycling. The death of living organism implies a flux from living POC to non-living POC (arrow #4), which then sinks (arrow #5) and is dissolved and mineralized to DOC (arrow #6) and—thanks to bacterial activity- DIC (arrow #9). During mineralization, nutrients are released in their dissolved inorganic forms too. Finally, several biogeochemical processes (exudation, assimilation, microbial food web activity) drive the fluxes that support exchange of carbon from DOC to POC (arrows #7 and 8). Those processes also depend on N and P concentrations.

It is therefore evident that the cycle of carbon in sea, as well as the exchange of carbon dioxide though the sea-air interface, depends on many biogeochemical processes, and is deeply interlaced with the cycling of nitrogen and phosphorus.

Coupled transport-biogeochemical models are based on the numerical integration of a system of a large number of partial differential equations, subjected to prescribed initial and boundary conditions. Given the complexity of the system and the large number of uncertainty sources present in any code, it is only fair to assume that any model, when run over a period as long as 100 years, would show some multi-year trend in some of its state variables. These trends might be representative of a real dynamic of the system, moving from its initial conditions in agreement with boundary conditions, and/or the result of some numerical drifts, or the cumulative effect of the propagation of some uncertainties. While to reduce the uncertainty sources in a model simulation surely is a meritorious activity, it is true that any complex model will always present some degree of dysfunctionality, approximations and uncertainty sources (Oreskes et al., 1994; Stow et al., 2009; Hipsey et al., 2020), so that—while keeping on improving model structures and parameterizations—it is important to be able to extract the largest fraction of useful information from currently available, though still imperfect, simulations.

A viable way to detect the effect of a specific forcing is to run pairs of simulations, differing only for that specific forcing, so that one simulation can be used as a “control,” and any difference between the two simulations can be related to the single changed factor.

In our case, this implies to get a “control” run by performing a 100-years long simulation forced by present day atmospheric conditions (i.e., by repeating for 10 times the atmospheric conditions used for the period 2000–2010) in addition to the 100-years long “scenario” simulations, in which the atmospheric forcing changes in agreement with the climate scenarios, all the rest being equal.

In reality, in our application, as a preliminary study, we performed not a single run but a number of “control” simulations, in order to fine tune a few highly uncertain model parameters, and to adjust the long-term dynamic to a somehow reasonable state. The model set-up has then been frozen, but for atmospheric forcing which have been replaced by those corresponding to the future scenarios. Finally, we computed the multidecadal signal in the “control” simulation, removed it (if present) from the “climate” scenario simulations, and analyzed the corrected results.

Figure 2 exemplifies what above described: the time series of the basin-scale averaged surface temperature and dissolved inorganic carbon concentration (taken as examples) are shown for the control run and the scenario simulation. Changes in the control run (that cannot be ascribed to climate change) are then estimated as the difference between the decadal means 2090–2099 and 2000–2009 in the control run (S_CTRL^FUT and S_CTRL^PRES, respectively); this bias is then removed from the difference in the scenario simulation between the decadal means 2090–2099 and 2000–2009 (S_CLIM^FUT and S_CLIM^PRES, respectively) to highlight the climate change signal ΔS:

where < . > denotes the decadal means.

To demonstrate the model capability to reproduce the most important properties of MS biogeochemistry, we compare the model results in the historical period 2000–2010 with a series of observational evidences.

The model correctly simulates the mean 3D structure of the chlorophyll field, in agreement with previous studies using models (e.g., Lazzari et al., 2012) and observations (e.g., Mignot et al., 2014; Lavigne et al., 2015). The mean surface concentration of the control run computed in the period 2001–2010 (Figure 4, top) clearly shows the east-west gradient and the highly productive areas in north-western MS and Alboran Sea. The model also properly reproduces the mean seasonal cycle of the averaged vertical chlorophyll profile in western and eastern sub-basins (Figure 4, bottom), featuring the winter bloom, mainly limited to the upper 100 m layer and with highest values in the western areas, and the summer stratification, characterized by the deep chlorophyll maximum (DCM) at depths increasing from 75 m in the western areas to 100 m in the eastern ones.

Model nutrients (phosphate and nitrate), dissolved oxygen, dissolved inorganic carbon, and alkalinity concentrations computed for the historical period are compared with a reference observational climatology derived from available moorings and vessels provided by the OGS National Oceanographic Data Center (NODC-OGS), spanning the period 2000–2011 (reported in Cossarini et al., 2015; Lazzari et al., 2016). Standard skill metrics (RMSD, bias, and correlation) are computed following a basin-aggregated approach for eight different vertical layers (Table 1) and compared to typical variability of the variables (i.e., standard deviation computed on spatial and temporal variability of the MS).

Carbonate system variables are characterized by a very high level of accuracy compared to their typical variability. Phosphate shows some inaccuracy in the surface layers (i.e., RMSD higher than the typical variability of the basin), however, it must be noted that in absolute terms the surface phosphate values (i.e., mean value is 0.02 ± 0.02 mmol/m³) and small absolute errors are amplified. The accuracy of oxygen, which is compared with its typical temporal variability, is very good in all layers except the deepest ones where the model overestimates of about 30 mmol/m³. Nitrate is the least accurate variable, with errors exceeding its typical variability in several layers. However, given the fact that phosphate is the limiting nutrient in the MS (Lazzari et al., 2016), it is reasonable to consider nitrate inaccuracy impacting less on the scenario reliability.

In general, correlation is larger than 0.7, confirming the ability of the model to correctly reproduce the basin-wide gradients along the water column, with the exception of dissolved inorganic carbon below 150 m. In this case, it is worth to remind that the spatial gradient of carbonate system variables is pretty low in the deep layers (Cossarini et al., 2015).

The long-term dynamics of the biogeochemical properties are assessed by analyzing the scenario simulations in 2090–2099 after removing the trend in the control simulation, if present (see Section 2). The analysis shows that at the end of the century, the concentration of nutrients (both nitrate and phosphate) increases in both epi- and mesopelagic layers, with larger increments in the surface layers and in the western part (around +17% with respect to the present conditions, and in partial agreement with Richon et al., 2019, specifically for nitrate), whereas in the bathypelagic layer we observe a very modest increase of nitrate in the western basin, a modest decrease in the eastern part, and no variation in phosphate (Table 2).

As a results of changes in temperature and nutrients, gross primary production increases in both western and eastern surface waters (as already suggested in Lazzari et al., 2014), with higher absolute variation in the west, but its increase is more than compensated by an increase in total respiration (absolute values), so that the epipelagic layer concentration of total carbon associated to living component (plankton) decreases (around −5% at west, −6 % at east), similarly to chlorophyll (−5% at west, −3% at east) and oxygen concentration (−6% at west, −4% at east).

In the mesopelagic layer, where primary production is not occurring, the decrease in living carbon component is even higher, reaching 15 and 20% in the western and eastern sub-basin, respectively. Changes are up to 30 and 70% in the deep layer, tough there the absolute values are two order of magnitude lower than in the mesopelagic layer.

DIC concentration increases everywhere, as a result of increased respiration, but also because of the increase in atmospheric CO₂ and of the changes in CO₂ solubilization driven by changes in temperature and salinity (see also Richon et al., 2019), with higher changes in the epipelagic layer and variation in the eastern part larger than that in the western one.

Similarly, pH decreases everywhere, but acidity variations are more marked in the epipelagic layer (−3.5%) and in the eastern basin, where significant changes are projected also in the bathypelagic part of the basin. The basin-scale reduction in pH results of about 0.3 in a century, and therefore similar to what is projected in the global ocean (Bopp et al., 2013; Jiang et al., 2019; and recently reviewed by Kwiatkowski et al., 2020), even if the MS alkalinity is higher than in the ocean. In fact, the relevant processes acting to change alkalinity level (namely evaporation/dilution and -to some extent- rivers input) also cause a corresponding change in DIC concentration, and the two factors have opposite effects on pH. It is evident that those change are a significant hazard for the Mediterranean marine ecosystems (Zunino et al., 2019, 2021).

The analysis of the alkalinity time series (not shown) reveals the presence of a small negative trend in the upper layers (values of around 0.6 mmol/m³/y Table 2), which is induced by changes in the Gibraltar boundary condition derived by the global ocean model, possibly resulting from a freshening of the Northern Atlantic induced by Greenland and Artic ice melting, or from changes in the ocean general circulation pattern close to the area (also see Reale et al., 2020).

The differential behavior of eastern and western basin may be also related to the weakening of the thermohaline circulation in the MS (as already observed by Adloff et al., 2015, who analyzed the same physical forcing adopted in this work), with the eastern basin being more sensitive to the strength reduction of the circulation deep cell.

The analysis of spatial distribution of changes at the regional and sub-basins level is not the major focus of the present manuscript, but it might be worth highlighting a few major features in the surface layer (Figure 5).

While the changes in pH are homogenous all over the basin, some heterogeneity is present in the alkalinity and DIC, which show variations of opposite sign and are generally spatially correlated. A larger degree of spatial variability is observed in nutrient concentrations, plankton primary production and chlorophyll density.

Indeed, from visual inspection of Figure 5 it is possible to recognize widespread positive anomaly regions in phosphate concentrations in the north-western part of the basin, the Adriatic Sea, large part of the Ionian sea, and the Aegean sea, which are contrasted by negative regions in the Levantine basin from the Rhodus gyre till Cyprus.

The primary production anomaly map is well correlated to the phosphate one, suggesting the presence of phosphorus limitation in plankton growth in most of the above-mentioned area, with the notable exception of the Aegean Sea, where plankton production decreases even if the phosphate concentration increases, possibly suggesting nitrogen limitation. The correlation between phosphate and primary production anomalies is observable also at smaller scales, such as the production hotspot in the Sicily channel or close to the Balearic islands.

The spatial distribution of total chlorophyll concentration resembles the primary production one, while nitrate concentration shows a general decrement, suggesting that nitrogen is not the limiting nutrient in most of the basin. Conversely, the pattern in the spatial distribution of nitrate anomaly might be co- shaped also by the uptake of nitrogen related to phosphorus-limited plankton growth.

It is clear, anyway, that different choices of future scenarios of land use and land changes might produce significant changes in the spatial pattern of nutrients, and as a consequence in plankton productivity ones. As an additional lateral comment, our model results (not shown here) point to the presence of different responses to climate scenarios from different plankton functional groups, with cascading consequences on different zooplankton groups. This result, however, requires a deeper analysis to be developed further and cannot be expanded in the present manuscript.

The effect of the increasing CO₂ atmospheric concentration is well captured also in Figure 6, representing the depth distribution along a west-east transect of the DIC anomaly, obtained as the difference between future (2090–2099 decadal means) and present (2000–2009) DIC concentration. The increment in atmospheric concentration forces higher CO₂ dissolution, resulting in higher concentration of surface DIC, which in turn drives significantly higher downward fluxes of DIC. The increment in DIC fluxes and concentrations is well recognizable down to 1,500 m, and more marked in the eastern sub-basin than in the eastern one. Carbon vertical fluxes can be related to the superposition and the cumulated effect of changes in solubility (solubility pump), vertical mixing and advection (physical pump), fluxes related to sinking of particulate organic carbon (biological pump), but also to the fluxes related to dense water formation and cascading along the continental shelves (continental shelf pump). Indeed, the strong increment all over the first 1,500 m between 20°E and 25°E is related to the spreading of the DIC rich dense water along the Adriatic and the Aegean shelves.

The analysis of our simulation also enables us to assess changes in the budgets of carbon stocks and fluxes for present and end of century conditions. Figure 7 summarizes the present situation (Figures 7A,B) and the main changes projected for the considered climate scenario (Figures 7C,D), both in terms of variable concentration and of variable total mass. In fact, while the budgets are assessed in terms of the distribution of the total mass of carbon through the different compartment, the physical and biogeochemical processes that support the mass transfer between compartments are often driven by differences in variable concentrations.

The simulation of the present day conditions shows the well-known east–west gradient in surface water POC concentration (Figure 7B), resulting from the superposition of the estuarine inverse circulation (defined by the inflow of surface water at Gibraltar Strait and a west to east advection driven flux, and the opposite east to west intermediate and deep fluxes) and the sinking of particulate organic carbon (the biological pump), which transfers carbon from the surface layer to the intermediate one, and while doing so reduces the entity of the eastward carbon flux (Crise et al., 1999; Crispi et al., 2001).

Model projections for future conditions indicate an increase in external inputs from Gibraltar and atmospheric fluxes (Figures 7C,D), which combined with changes in biological activity and internal circulations results in a reversal of the net global Atlantic-Mediterranean exchange, which changes from a slightly negative balance (export carbon to the Atlantic) to a slightly positive one (import carbon from the Atlantic). Globally, the system is slightly richer in carbon, with slightly higher horizontal exchanges among sub-basins, and higher total eastern carbon fluxes.

By contrast, the stock of organic particulate carbon, POC, and of surface dissolved organic carbon will decrease all over the basin (8%), as will the stock of bottom POC (Figures 7A,C), possibly because of the reduction in sinking from surface water and of the increase in respiration.

At bottom, the scenario projects an increase of the pool of dissolved organic carbon at east (1%) and a decrease at west (−3%) as a response of the thermohaline transport weakening, as indicated by the flux of DOC at the Sicily Strait.

The largest increment of the carbon pool is observed in the stock of dissolved inorganic carbon at bottom (+2,759 TgC) whereas the increase at surface is of 671 TgC (Figures 7A,C) even if the increase of the concentration is about 5% at surface and 2% at surface and bottom, respectively, as a result of the slow deepward invasion of absorbed atmospheric CO₂ (+360 and 400% in the eastern and western sub-basins, respectively. The eastern sub-basin accumulated more than 85% of the newly assimilated carbon, because its volume is 60% greater and its rate of accumulation in the deeper layer per unit of surface is 50% higher than the western one (Figures 7B,D). The higher carbon accumulation efficiency in accumulation is a direct consequence of the lowering of the intensity of the thermohaline circulation.

The analysis of the metabolic fluxes (arrows in Figure 7) reveals that, in spite of the increase of atmospheric CO₂ input and the increase in temperature which will fuel ecosystem metabolism, the increment in the respiration processes appear to be more relevant than the increment in photosynthesis, with the final results of smaller concentrations in POC and total organic matter, and higher DIC concentration and ocean acidity.