Projected Effects of Climate-Induced Changes in Hydrodynamics on the Biogeochemistry of the Mediterranean Sea Under the RCP 8.5 Regional Climate Scenario

1. Introduction

Since the beginning of the industrial era, a vast amount of greenhouse gases have been released into the atmosphere through human activities. For example, the average atmospheric CO₂ concentration, which was around 280 ppm in the pre-industrial era, had reached 408 ppm in 2019 (Dlugokencky and Tans, 2019). These high concentrations of greenhouse gas in the atmosphere have accentuated the natural greenhouse effect and induced a global rise of the temperature of 1°C in 2018 (Allen et al., 2018), leading to major climate changes. However, in some regions, including the Mediterranean region, the temperature rise is expected to be higher than the global increase (Seneviratne et al., 2016) due to local characteristics such as land-use and/or land-cover (Pitman et al., 2009), urban development (Wilby, 2008) or aerosols (Levy et al., 2013).

The Mediterranean region has already been demonstrated to be particularly exposed to climate change by Giorgi (2006). The observed trend for the rate of climate change in the Mediterranean region is indeed higher than the global tendency (Cramer et al., 2018). For example, the average temperature is 1.4°C above the pre-industrial value and heat waves as well as droughts occur more frequently (Cramer et al., 2018). Furthermore, the Mediterranean region is increasingly impacted by human activities, with a population in North Africa that rose 4-fold between 1960 and 2015, while urbanization increased by a factor of 2 (Cramer et al., 2018). As a semi-enclosed sea surrounded by land and numerous industrialized countries, the Mediterranean sea (MS) is highly sensitive to climate change. Because of the specific response to climate change of the Mediterranean region (and consequently of the MS), there is a need to model climate scenarios at regional scale in order to fully cover the nature and the magnitude of the expected changes in this area.

In the MS, most hydrodynamic studies have predicted a global increase in temperature and salinity during the next century as well as a decline in deep water formation and a shallowing of the Mixed Layer (Somot et al., 2006; Durrieu de Madron et al., 2011; Adloff et al., 2015; Darmaraki et al., 2019). However, it has been shown recently by Soto-Navarro et al. (2020), who analyzed a set of 25 simulations, that the predicted warming is not homogeneous throughout the MS. All the models in this study showed a high rate of warming in the Levantine Basin, the north of the Aegean Sea, the Adriatic Sea, and the Balearic Sea. They also highlighted that despite the overall salinity increase throughout the MS, the surface of the Western Mediterranean Basin (Western Basin) will be mainly impacted by the characteristics of the Atlantic waters inflowing into the Western Mediterranean Basin due to the higher water input of north-eastern Atlantic waters predicted by most general circulation models (GCMs). Finally, Soto-Navarro et al. (2020) pointed out that these changes will probably lead to a reduction in deep water formation in the Gulf of Lions.

The dynamics of the Mixed Layer Depth (MLD) directly affects the planktonic food web through the availability of their essential nutrients [Nitrate (NO₃) and Phosphate (PO₄)] in the euphotic layer. In the open ocean, the availability of nutrients in the euphotic layer is mainly driven by the exchanges between the nutrient-depleted upper layers and the nutrient-rich deep layers through vertical mixing (Broecker and Peng, 1982; Van Den Broeck et al., 2004). In the mid-latitude areas (like the MS), the seasonal dynamics of the MLD is mainly responsible for these exchanges (Williams, 2003; Auger et al., 2014). In the MS, between May and November, the MLD is shallow and strong stratification separates the upper layers from the deep layers (Lionello et al., 2006). From December to April, because of strong winds and low temperatures that induce the vertical mixing of the water column, the MLD becomes deeper (Houpert et al., 2015). This process allows the homogenization of the water mass above the MLD and injects nutrients into the upper layers.

To our knowledge, only a few studies have investigated the effects of climate change on the biogeochemistry of the whole MS with a 3D coupled hydrodynamic/biogeochemical model. Most of those studies, except Moullec et al. (2019), have predicted a decline in primary production for the Western Basin and nutrient inputs to the upper layers due to a reduction in vertical mixing (Macias et al., 2015; Richon et al., 2019). Some studies also highlight a reduction of the phytoplankton biomass (mostly for diatoms) (Lazzari et al., 2014; Richon et al., 2019), while Moullec et al. (2019) show an increase in phytoplankton biomass. Finally Richon et al. (2019) show a deepening of the deep chlorophyll maximum (DCM) while Macias et al. (2015) show that it will become slightly shallower.

This study is intended to be a further contribution to present knowledge in this field. All these previous studies, except that of Lazzari et al. (2014), use a fixed stochiometry for some or all of the unicellular plankton functional types in the model. Moreover, some of these studies use discontinuous (i.e. not all the period has been simulated) simulations (Lazzari et al., 2014) and/or old SRES climate change scenarios (Lazzari et al., 2014; Richon et al., 2019). Furthermore, due to the scarcity of available data, all of them used rivers' nutrient inputs that are not consistent with the scenario used, and this is also the case for the present study.

Here we present the first study based on a biogeochemical model that considers a flexible stoichiometry for all the Plankton Functional Types (PFTs) (with the exception of copepods), and that also calculates the abundance of each PFT.

Furthermore, apart from considerations regarding the structure of the biogeochemical model, the biogeochemistry of the MS in the next century will be driven by the hydrodynamics of the MS, as well as by the input of nutrients at the different boundaries (i.e., inputs of nutrients by rivers, through the Gibraltar Strait, or by atmospheric deposition). In this study, our aim is to focus on the effects of the variations (expected according to the RCP 8.5 scenario) in hydrodynamics alone, apart from the effects induced by future changes in the inputs at the different boundaries. The reasons for this choice are related to our aim of exploring the relative weight of each of the internal and external contributions to the future changes in the MS. This choice also has the advantage of avoiding the uncertainties associated with the evolution of the external forcings (see the Discussion section for more details).

For this purpose, high-resolution simulations over the 1985–2100 period have been run using the Eco3M-Med biogeochemical model forced by the RCP 8.5 emission scenario for the twenty-first century provided by the high-resolution and fully-coupled regional climate system model CNRM-RCSM4 (Sevault et al., 2014; Darmaraki et al., 2019), while nutrient inputs by rivers and through Gibraltar Strait were kept unchanged.

2. Materials and Methods

2.1. The Biogeochemical Model: Eco3M-Med

The biogeochemical model used in this study is the Eco3M-Med configuration (Alekseenko et al., 2014; Guyennon et al., 2015; Pagès et al., 2020). It is coupled off-line with the hydrodynamical model NEMO (Nucleus for European Modeling of the Ocean, Madec and NEMO-Team, 2016). It is embedded in the Ecological Mechanistic and Modular Modeling (Eco3M) tool for biogeochemical modeling (Baklouti et al., 2006a) which is based, as far as possible, on mechanistic descriptions of biogeochemical fluxes (Baklouti et al., 2006b; Mauriac et al., 2011; Alekseenko et al., 2014; Guyennon et al., 2015).

The Eco3M-Med configuration has been designed to represent the main characteristics of the MS. It includes six different plankton functional types (PFTs), two phytoplankton compartments, one large phytoplankton (PHYL) (i.e., phytoplankton with size above 10 μm) and one small phytoplankton (PHYS) (i.e., phytoplankton with size below 10 μm), three zooplankton compartments, the copepods (COP), the ciliates (CIL) and the heterotrophic nanoflagellates (HNF) as well as a compartment of heterotrophic bacteria (BAC). Three pools of inorganic matter are also represented, namely nitrate (NO₃), ammonium (NH₄), and phosphate (PO₄) as well as a pool of dissolved organic matter (DOM) and two compartments of detrital particulate matter, one with small particles (DETS) and the other with large particles (DETL). More information about the Eco3M-Med configurations is available in Guyennon et al. (2015) and Pagès et al. (2020).

As in the previous versions of the Eco3M-Med model, the intracellular ratios [Q_XY, where X and Y either stand for nitrogen (N), phosphorus (P) or carbon (C)] in organisms, as well as their cell quotas Q_X (intracellular quantity of X per cell) are flexible and dynamically calculated by the model. These intracellular ratios and quotas are used to regulate the different biogeochemical fluxes represented by the model. As already mentioned in several previous papers, the use of abundances in addition to biomass to characterize PFTs offers, among other advantages, the opportunity to calculate intracellular quotas which are absolute values giving access to the actual degree of limitation of organisms by each nutrient. More details on these aspects of the model can be found in Alekseenko et al. (2014), Guyennon et al. (2015), and Pagès et al. (2020).

Another important feature of this version of Eco3M-MED is that phytoplankton can directly uptake dissolved organic phosphorus (DOP). In this way, we implicitly consider the DOP hydrolysis by extracellular enzymes (i.e., alkaline phosphatase, nucleotidase, polyphosphatase, or phospho-diesterase) produced by heterotrophic bacteria in depleted conditions (Perry, 1972, 1976; Vidal et al., 2003), thereby making the hydrolyzed DOP also available for phytoplankton.

Finally, in the Eco3M-MED model, there is no regulation of kinetics by temperature. The use of PFTs enables us not only to include implicitly marine diversity in a simplified way (Leblanc et al., 2012), but also to implicitly consider that each external temperature will correspond to the optimum temperature of a given species inside each PFT. However, some minor changes since the version used in Pagès et al. (2020) have been made, since the previous version failed to reproduce the well-marked Deep Chlorophyll Maximum (DCM) in the Eastern Basin. To overcome this issue, the maximum CHL:C ratio has been correlated with available light and more specifically with photosynthetically available radiation (PAR), thereby allowing more efficient photo-acclimation in the model, as given by Equation (1):

As expected, the DCM intensities were stronger, especially in the Eastern Basin, which resulted in a significantly larger quantity of large DETL particles produced at the DCM depth. Since the latter sink at a rate of 50 m per day, it was necessary to reduce their turnover time (from 15 to 1 day) to avoid DETL accumulation in deep waters and too rapid changes in nutrient concentrations at depth.

2.2. The Physical Forcing: CNRM-RCSM4

In this study, the physical forcings of the biogeochemical Eco3M-Med model come from a 1950–2100 scenario simulation performed with the fully-coupled high-resolution Regional Climate System Model CNRM-RCSM4 dedicated to the study of the Mediterranean climate and Mediterranean Sea (Sevault et al., 2014).

The CNRM-RCSM4 simulates the main components of the Mediterranean regional climate system and their interactions. It includes four different components: (I) The atmospheric regional model ALADIN-Climate (Radu et al., 2008; Colin et al., 2010; Herrmann et al., 2011) characterized by a 50 km horizontal resolution, 31 vertical levels, and a time step of 1800 s, (II) the ISBA (Interaction between Soil Biosphere and Atmosphere) land-surface model (Noilhan and Mahfouf, 1996) at a 50 km horizontal resolution, (III) the TRIP (Total Runoff Integrating Pathways) river routine model (Oki and Sud, 1998), which is used to convert the runoff simulated by ISBA into rivers at a 50 km horizontal resolution (Decharme et al., 2010; Szczypta et al., 2012; Voldoire et al., 2013), and (IV) the Ocean general circulation model NEMO (Nucleus for European Modeling of the Ocean, Madec and NEMO-Team, 2016) in its NEMO-MED8 regional configuration (Beuvier et al., 2010). NEMO-MED8 is characterized by a horizontal resolution of 1/8° (grid cells size from 6 to 12 km), a vertical resolution of 43 vertical levels (cell height ranging from 6 to 200 m), and a time step of 1200 s. It covers the MS and a buffer zone (from 11 to 6°W) including part of the Atlantic ocean where boundary conditions are applied for temperature and salinity. Note that the rivers runoffs are fully coupled owing to the TRIP component. This means that freshwater inputs to NEMO-MED8 vary during the simulation according to precipitation and evapotranspiration over land. More details about the CNRM-RCSM4 model can be found in Sevault et al. (2014).

The physical simulation used in the current study is a projection simulation performed with CNRM-RCSM4 under the socio-economic scenario RCP 8.5 and described in Darmaraki et al. (2019). The simulation was run after 130 years of spin-up allowing us to stabilize the model state. The initial conditions came from MEDATLAS-II (Rixen et al., 2005) and the lateral boundary conditions for the atmosphere and ocean regional models came from a RCP 8.5 scenario simulation performed with the CNRM-CM5 General Circulation Model (Voldoire et al., 2013) participating in the CMIP5 multi-model exercise. The physical simulation is divided into two periods: a historical period from 1950 to 2005 and a scenario period from 2006 to 2100, here using the RCP 8.5 IPCC emission scenario. In order to ensure a present-climate state closer to observations, we applied two bias corrections constant in time during the CNRM-RCSM4 simulation: (1) to the atmosphere heat fluxes provided to the ocean by the atmosphere component of the regional coupled model and (2) to the near-Atlantic temperature and salinity 3D oceanic conditions provided by the driving GCM CNRM-CM5 (see Supplementary Material of Soto-Navarro et al., 2020 for more information). In addition we performed a control simulation with CNRM-RCSM4 parallel to the projection run (not shown) in order to verify the model stability. This simulation has already been used by Darmaraki et al. (2019) to study to occurrence of marine heat waves during the next century, by Moullec et al. (2019) to investigate the effect of climate change on the trophic web of the MS (end to end model) and by Soto-Navarro et al. (2020) to study the evolution of the seawater properties during the next century.

2.3. Characteristics of the Simulations

In this study, we ran two simulations with Eco3M-Med. The first simulation is a control covering only the 2006–2100 period. The second simulation includes a historical period between 1985 and 2005 and a scenario period between 2006 and 2100. The control simulation was run in order to ensure that the biogeochemical model does not significantly drift over the scenario period. The configuration of the control simulation is similar to the scenario simulation except for the hydrodynamic forcing, which has been built from the random draw of years among the historical period (1985–2005). To cover the 2006–2100 period we proceeded to several random samplings of 20 years without replacement into the historical period in order to remove any trend present in the historical period. The goal is to simulate a permanent and random present-climate state.

The initial conditions for the biogeochemical model come from the outputs of the hindcast simulation of Pagès et al. (2020) obtained in 1985. This simulation used almost the same biogeochemical model as the present one but with a different spatial resolution (1/12° instead of 1/8°). An additional 10-year spin-up has then been run (years 1985–1990 have been run twice).

Concerning nutrients input at the Gibraltar Strait (the only open boundary), the same procedure as in Pagès et al. (2020) has been applied: climatological data (monthly vertical profiles) were used in the buffer zone (previously defined) for each state variable of the biogeochemical model. The profiles for PO₄, NO₃, NH₄ and for the chlorophyll a concentration (CHL) were extracted from the World Ocean Atlas database while the profiles for the remaining variables were derived from the CHL profiles using conversion factors derived from published works (see Alekseenko et al., 2014 for more details). As already mentioned, nutrients concentrations in the buffer zone remained equal to their climatological monthly values during all runs.

For the rivers inputs, monthly values of NO₃ and PO₄ have been provided for the 29 main rivers all around the MS by the study of Ludwig et al. (2009). These nutrient inputs are available up to the year 2000. Between 2001 and 2100, the amount (in Kt.a⁻¹) of nutrient river inputs of the year 2000 are repeatedly used as already done in Macias et al. (2015). This implies that at rivers mouths, nutrients concentration can slightly vary, but in the end, the input of nutrients remains constant at the scale of the MED.

Finally, three 20-year periods will be distinguished in this paper to analyse the results: the 1986–2005 period denoted HIST, the 2031–2050 period denoted MID, and the 2081–2100 period denoted END.

2.4. Simulation Analyses

All the variables discussed in the following sections are calculated from monthly values through three steps:

(i) For each grid cell and each year, the yearly average is calculated.

(ii) The yearly values are spatially averaged over the Western and the Eastern Basins, excluding coastal areas where the depth is <200 m (see Figure 1) as done for example in Richon et al. (2018) (cf section 3.1).

(iii) These yearly values are averaged over the three 20-year periods, namely: the HIST period (1986–2005), the MID period (2031–2050), and the END period (2081–2100). The standard deviations (SD) provided in Table 1 are computed on the same 20-years periods.

Some variables, however, are averaged over specific seasons, to allow a more in-depth analysis. This is the case for the Mixed Layer Depth (MLD) which is an on-line and high-frequency diagnostic variable computed directly within the physical simulation (Darmaraki et al., 2019; Soto-Navarro et al., 2020). Starting from the surface, it corresponds to the first depth at which the potential density starts to differ from the density at 10 m by more than 0.011 kg m⁻³. In the present study, we focus on the winter maximum of the MLD instead of annual means. The winter maximum of MLD indeed directly constrains the amount of nutrients that will be brought to the surface layers during winter mixing, and therefore will sustain the biological production in the MED (Broecker and Peng, 1982; Williams, 2003; Van Den Broeck et al., 2004; Auger et al., 2014).

In the same way, the depth of the top nitracline (Ncline) and phosphocline (Pcline) are computed in summer (July to September) when stratification and depletion of nutrients in surface layers are the strongest. The top of the Ncline is defined as the first depth, starting from the surface, at which NO₃ equals the threshold value of 0.05 μmol·l⁻¹ which corresponds to the classical quantification limit of nitrate measurement using an autoanalyzer (Aminot and Kérouel, 2007). The top of the Pcline is based on the threshold value of 3 nmol·l⁻¹ which is the threshold for the top of the Ncline divided by 16, namely the Redfield N:P ratio.

The DCM is also calculated during summer (July to September), i.e., when it is present in most of the MED. Finally, the Stratification Index (SI) is computed in summer (July to September) from the surface to a given depth using the Brunt-Vaisala frequency which is an indicator of the strength of the vertical stratification. The latter increases with the SI. It has been used in many Mediterranean studies (Lascaratos and Nittis, 1998; Herrmann et al., 2008; Somot et al., 2016).

In order to assess whether the changes between the HIST, MID, and END periods are significant, we applied a Mann-Whitney U test (a non-parametric version of the Student t-test). If the p-values are lower than 0.05, the variations are considered as statistically significant, and if the p-values are above 0.05, they are not (see the gray cells in Table 1).

3. Results

3.1. Model Skill Assessment

The biogeochemical model configuration (Eco3M-Med) is almost the same as the one used in Pagès et al. (2020) in which its capacity to successfully reproduce the main patterns of the biogeochemistry of the MS (i.e., surface CHL, nutrients distribution, primary production, etc.) has been demonstrated with a physical forcing at higher resolution (1/12° instead of 1/8°). For this reason, only an additional comparison is provided here (Figure 2). This figure shows the modeled (mCHL) and the observed (oCHL), i.e., satellite-derived, mean surface chlorophyll a concentration over the 1999–2006 period. In the Western Basin, both maps exhibit high chlorophyll concentration in the Gulf of Lions, along the Algerian coast and in the Alboran Sea. However, the modeled mean concentration (0.12 μg·l⁻¹ over the Western Basin) is a little bit lower than the observed one (0.16 μg·l⁻¹). Regarding the Eastern Basin, the model succeeds in representing the oligotrophic offshore area, but the mCHL is too low in several areas: (i) in the Gulf of Gabes where several factors may explain this discrepancy between model outputs and satellite data, for example the shallow depth of this area which may mean that the satellite-based data is polluted by the reflection from benthic vegetation (Jaquet et al., 1999). Another complementary explanation is that phosphate inputs in this industrial area are not taken into account by the model, though some recent work managed to simulate relatively high phytoplankton biomass in this region without additional phosphate inputs (Macias et al., 2019), (ii) in the vicinity of the Nile mouth where the nutrient inputs are very poorly known, as well as throughout the Western Adriatic Sea probably due to a lack of phosphate in the Pô River inputs, (iii) at the scale of the Eastern Basin, the mean values of mCHL and oCHL are respectively equal to 0.07 and 0.12 μg·l⁻¹. However, if we exclude the coastal area (i.e., where the depth is <200 m) in order to avoid case II water (Bosc et al., 2004; Richon et al., 2018), the mean values of mCHL and oCHL are respectively equal to 0.055 and 0.06 μg·l⁻¹. Furthermore, to rule out the assumption that the changes in the biogeochemistry cannot be ascribed to model drift, a control simulation has been run over the forecast period, namely from 2006 to 2100 (see section 2.4 for details). This control simulation did not exhibit any trend regarding the state variables, thereby enabling us to conclude that there is no drift inherent to the biogeochemical model.

On that basis, and accounting for the numerous previous comparisons with data performed for the historical period (Guyennon et al., 2015; Pagès et al., 2020), it is considered that the present model can be used to investigate the impact of climate change on the basin-scale biogeochemistry over the next century.

3.2. The Mixed Layer Depth

During the HIST period, the average winter MLD (Figure 3A) is deeper in the area of deep water formation (i.e., the Gulf of Lion, South Adriatic Sea, and the South Aegean Sea) and shallower in most of the other parts of the MED. Between the HIST and MID periods (Figure 3B) the winter MLD exhibits a slight decrease of almost 5.7% in the Western Basin (from 87 m for the HIST period to 82 m for the MID period) and of 2.8% in the Eastern Basin (from 122 m for the HIST period to 118 m for the MID period). However, those decreases are not statistically significant. Between the HIST and END periods, the winter MLD shows (Figure 3C) a substantial decrease of 21% in the Western Basin (from 87 m for the HIST period to 68 m for the END period) and of 8% in the Eastern Basin (from 122 m for the HIST period to 112 m for the END period). The Gulf of Lions region, and more generally the Western Basin, is the most affected, showing the strongest reduction in the maximum MLD, especially after 2050.

Regarding the stratification index during the HIST period (Figure 3D), the area of deep water formation shows the weakest values. Furthermore, the Western Basin shows an SI stronger than the Eastern Basin with respectively 1.3 and 1 m²·s ⁻². The variation of SI between the HIST and MID periods (Figure 3E) shows an increase of 11% (from 1.3 for the HIST period to 1.45 m²·s ⁻² for the MID period) in the Western Basin and of 5% in the Eastern Basin (not statistically significant). Between HIST and END periods (Figure 3F) the increase is even more substantial with 29% (from 1.3 for the HIST period to 1.68 m²·s ⁻² for the END period) in the Western Basin and of 16% in the Eastern Basin (from 1 for the HIST period to 1.24 m²·s ⁻² for the END period). For SI as for the MLD, the Western Basin is more impacted than the Eastern Basin.

The variations between the three different periods mostly occur between the MID and the END periods, i.e., between 2050 and 2100. The details of variations between each period are given in Table 1. Further information on the variations in the MLD in the CNRM-8.5 simulation can be found in Soto-Navarro et al. (2020).

3.3. Scenario Results Associated With the Biogeochemistry of the MS

3.3.1. Nutrient Concentrations

The annual average surface concentration of NO₃ and PO₄ in the Western and Eastern Basins (Figure 4) do not show any clear trend up to the year 2050.

After 2050 in the Western Basin, the surface NO₃ concentration (Figure 4A) progressively declines from 0.42 μmol·l⁻¹ in the MID period to 0.22 μmol·l⁻¹ for the END period, corresponding to a decrease of 48.3% between MID and END (50% between HIST and END, Figure 5).

For the surface PO₄ concentrations in the Western Basin (Figure 4C) after 2050, we observe a decrease from 0.027 μmol·l⁻¹ for the MID period to 0.017 μmol·l⁻¹ for the END period. This represents a 38% decrease between MID and END (and a 43% decrease between HIST and END (Figure 5).

In the Eastern Basin, we observe the same trends for both NO₃ and PO₄ as for the Western Basin (Figures 4B,D). However, the decrease between the MID and the END periods is less than in the Western Basin and equal to 37.6% for NO₃ and 37% for PO₄. Between HIST and END the decrease reaches 42% for NO₃ and 54% for PO₄ (Figure 5).

Overall, the decline in nutrient concentrations is stronger in the Western Basin than in the Eastern Basin between MID and END. We can also observe that NO₃ concentrations decrease more than those of PO₄ between MID and END (respectively by 45% and by 37%) in the Western Basin, while in the Eastern Basin between MID and END the decrease is almost the same 37% for both nutrients. Between HIST and END, NO₃ concentrations still decrease more than those of PO₄ in the Western Basin (respectively by 50% and by 43%). But in the Eastern Basin, NO₃ concentrations decrease less than those of PO₄ (respectively by 42% and by 54%). All the values are available in Table 1.

The model also predicts a deepening of the top of the nitracline (Ncline) and the phosphocline (Pcline) after 2050 (Figure 5). The top of the Ncline remains constant in both basins until 2050, with an average depth of 53 and 64 m in the Western and Eastern Basins, respectively. Thenceforth, in 2050 Ncline starts to deepen to end up in the END period at 56 m in the Western Basin and 67 m in the Eastern Basin a deepening of respectively 6.9 and 5% between HIST and END (Figure 5). Regarding the Pcline, the same trends are shown by the model as for the Ncline, with an average depth of 51 m (68 m) for the Western (Eastern) Basin, for the MID period, and 50 m (70 m) for the END period, corresponding to a deepening of 0.7% (not statistically significant) in the Western Basin and to 3.2% in the Eastern Basin. Between HIST and END, the Pcline deepens of 4.5 and 8%, respectively in the Western and Eastern Basins.

3.3.2. Chlorophyll Concentration

The annual average of CHL exhibits a strong decline after 2050 in both basins. In the Western Basin (Supplementary Figure 1A), CHL quickly drops from 0.12 (0.11) μg·l⁻¹ in the HIST (MID) periods down to 0.08 μg·l⁻¹ in the END period, representing a decrease of ~31% (25.4%) (see Figure 5). In the Eastern Basin (Supplementary Figure 1B), CHL exhibits the same trends as in the Western Basin, the decrease is equal to 21% (16%) between HIST (MID) and END, CHL varies from 0.058 (0.054) to 0.045 μg·l⁻¹. Looking at the depth of the Deep Chlorophyll Maximum (DCM) during summer (July to September) in the Western Basin (Supplementary Figure 1C), it can be observed that, the DCM deepens from 64 m (66 m) in the HIST (MID) period to 69 m in the END period, a deepening of 7.3% (6%) (see Figure 5). In the same way, in the Eastern Basin, no trend could be observed for the DCM between HIST and MID period, with an average DCM depth of 90 m. After 2050, the DCM deepens from 90 m in the MID period to 94 m in the END period, a deepening of 5.5% (Figure 5).

3.3.3. Phytoplankton Community

According to the model, the vertically-integrated cellular abundance of phytoplankton remains stable over the first half of the twenty-first century almost everywhere in the Mediterranean. After 2050, the vertically-integrated abundance of small phytoplankton (${\text{PHYS}}_{int}^{cell}$) remains stable over both basins (changes are not statistically significant), as shown in Supplementary Figure 2 (left). The vertically-integrated cellular abundance of large phytoplankton (${\text{PHYL}}_{int}^{cell}$) shown in Supplementary Figure 2 (right), decreases by 26% in the Western Basin and by 35% in the Eastern Basin between the MID and the END periods.

3.3.4. Primary Production Rates

In both basins, the annual vertically-integrated net primary production (NPP_int) remains stable (Supplementary Figures 3A,B) between the HIST and the MID periods (variations are not statistically significant, p-value > 0.05). In the Western Basin, NPP_int declines from 32 (30.6) gC·m⁻²·a⁻¹ during the HIST (MID) period to 27.6 gC·m⁻²·a⁻¹ at the END period, a reduction of 13.8% (10%) (see Figure 5). In the Eastern Basin, NPP_int declines from 20.3 to 17.3 gC·m⁻²·a⁻¹ between the HIST and the END periods (i.e., a reduction of 14%).

Regarding the annual vertically-integrated gross primary production (GPP_int), it follows the same trend as NPP_int but with higher values. In the Western Basin, GPP_int declines from 62.9 gC·m⁻²·a⁻¹ during the HIST period to 53.6 gC·m⁻²·a⁻¹ for the END period (i.e., a reduction of 14.7%). In the Eastern Basin, NPP_int declines from 43.3 to 38.1 gC·m⁻²·a⁻¹ between the HIST and the END periods (i.e., a reduction of 12%).

3.3.5. Copepod Abundance and Biomass

The mean annual vertically-integrated carbon biomass (${\text{COP}}_{int}^C$) of copepods decreases in both basins after 2050 (Supplementary Figure 3, top). ${\text{COP}}_{int}^C$ declines between the MID and END periods by 7% in the Western Basin and 6% in the Eastern Basin. Regarding the vertically-integrated abundance (${\text{COP}}_{int}^{cell}$) of COP (Supplementary Figure 3, bottom), we observe a decline from 159 to 147 ind·m⁻² between the MID and the END periods, representing a 8% decrease in the Western Basin. In the Eastern Basin, between MID and END, ${\text{COP}}_{int}^{cell}$ decreases from 147 to 139 ind·m⁻², which corresponds to a decline of 5.4%.

3.3.6. Particulate and Dissolved Organic Carbon

The model also highlights changes (Figure 5) in Detrital Organic Carbon (DETC). In both basins, DETC integrated concentrations remain stable between HIST and MID (values are not statistically significant), with an average value of 123 and 75 mmolC·m⁻², respectively for the Western Basin and the Eastern Basin. Then they strongly decline in both basins from 124 and 76.4 mmolC·m⁻², respectively for the Western Basin and the Eastern Basin during HIST period to end up at 103.8 and 64 mmolC·m⁻² at the END period, a decrease of 16% in both basins. Concerning the Dissolved Organic Carbon (DOC) the model did not exhibit any trend despite the slight increase in PHYS abundance (not shown).

3.3.7. Export of Organic Carbon

The annual average of particulate and dissolved organic carbon export at 1000 m (${\text{EXP}}_{1000m}^{org}$) remains steady up to 2050 in both basins with average values of 0.32 gC·m⁻²·a⁻¹ for the Western Basin (Supplementary Figure 5A) and 0.16 gC·m⁻²·a⁻¹ for the Eastern Basin (Supplementary Figure 5B). After 2050, the ${\text{EXP}}_{1000m}^{org}$ starts to decline to end up at 0.25 gC·m⁻²·a⁻¹ in the Western Basin and 0.13 gC·m⁻²·a⁻¹ in the Eastern Basin, which corresponds to a decrease of respectively 20 and 22% between the HIST and the END periods, respectively (Figure 5).

3.3.8. Nutrients Availability After Winter Mixing

Figure 6 shows the mean nutrient profiles in autumn (September to November) and the position of the maximum MLD, averaged over the Western Basin, for the HIST, MID, and END periods. It should be noted that for better legibility, NO₃ and PO₄ concentrations in Figure 6 are shown on a different horizontal axis. The slope of the Ncline between its top and 100 m depth is weaker than the slope of the Pcline, and the annual maximum of MLD becomes shallower from HIST to END. Furthermore, the intersection between the deeper MLD and autumn nutrient profiles provides a rough estimation of the nutrient concentrations mixed with upper surface water. For the HIST period, winter mixing reaches a conceptual value of 3.42 (0.2) μmol·l⁻¹ of NO₃ (PO₄). These values drop to 2.98 (0.17) and 1.11 (0.08) μmol·l⁻¹ for NO₃ (PO₄), respectively during the MID and the END periods. This represents a 67% decrease in the NO₃ concentrations induced by winter mixing between HIST and END, and a 60% decrease in PO₄. In the Eastern Basin, the MLD variation (not shown) is statistically significant only between the HIST and the END periods. For the HIST period, winter mixing reaches a conceptual value of 2.59 (0.13) μmol·l⁻¹ of NO₃ (PO₄). These values decrease to 1.66 (0.08) μmol·l⁻¹ for NO₃ (PO₄), for the END period. This represents a 35% decrease between HIST and END in the NO₃ concentration potentially reached by winter mixing, and a 38% decrease for that of PO₄.

3.3.9. Growth-Limiting Nutrients

The most limiting nutrient between 0 and 150 m among N and P during summer (which corresponds to the period where nutrients most limit phytoplankton growth) is shown for PHYS and PHYL in Figure 7. The most limiting nutrient X corresponds to the element for which Q_X (in mol X cell⁻¹) is minimum, and ${\tilde{Q}}_X$ is the mean normalized (i.e., comprised between 0 and 1) intracellular quota in element X. Furthermore, we consider that there is co-limitation when the difference between ${\tilde{Q}}_N$ and ${\tilde{Q}}_P$ does not exceed 5%. During the HIST period, PHYS (Figure 7A) is mainly limited by phosphorous (P) in the Western Basin (30% of its area) or co-limited by N and P (52% of its area). Nitrogen (N) limitation only occurs on 18% of the Western Basin area. During the END period, PHYS (Figure 7B) is mainly limited by N (42% of its area) or co-limited (41% of its area). P limitation only occurs on 17% of the Western Basin area. In the Eastern Basin during the HIST period, PHYS is mainly limited by P (75% of the area). N and N-P colimitation only occur on respectively 15 and 8% of the Eastern Basin area. During the END period, the limitation remains almost unchanged as compared to the situation during the HIST period.

For PHYL (Figure 7C), during the HIST period, 46% of the area of the Western Basin is limited by P, 45% is co-limited and only 9% is limited by N. During the END period (Figure 7D), the limitation in N increases to reach 33% of the area of the Western Basin, the P limitation decreases down to 19%, and the co-limited area remains almost unchanged. In the Eastern Basin, for the HIST and the END periods, the limitation remains driven by phosphorus for more than 80% of its area.

4. Discussion

4.1. Relevance of the Present Scenario

This scenario explores the effects of the expected changes in hydrodynamics (according to the RCP 8.5 scenario) on the biogeochemistry of the MS. The other potential sources of variability of the biogeochemistry, namely the inputs of nutrients at the different boundaries (i.e., inputs of nutrients by rivers, at the Gibraltar Strait, or by atmospheric deposition) have deliberately not been considered at this stage. There are several reasons for this: (i) first, it is very informative to know the relative weight of each of the internal and external contributions to the future changes in the MS and one of the main strengths of modeling is precisely that it enables us to disentangle the role of each process and forcing. Aware of the potential role that could be played by river inputs throughout the MS, and not only near the river mouths (Macias et al., 2014; Moon et al., 2016; Pagès et al., 2020), we deliberately decided to maintain nutrient river outputs constant (the year 2000 is repeatedly used after 2000). This study has therefore enabled us to explore the changes ascribed to hydrodynamics alone, independently of the other sources of variability. It also provides a reference for future studies focused on the understanding of the role of nutrients boundary inputs on the biogeochemistry of the MS, (ii) unlike scenarios for hydrodynamics, the only available scenario for nutrient inputs by rivers provided by Ludwig et al. (2010) is very rudimentary since it provides two values (one for 2030 and another for 2050) for each of four scenarios of the Millenium Ecosystem Assessment. Such nutrients input scenarios are therefore liable to introduce a considerable additional uncertainty in the model outputs, all the more so as they are not consistent with the RCP 8.5 climate scenario used in the present study. Concerning atmospheric deposition, to our knowledge there is no transient scenario available so far describing the evolution of atmospheric deposition over the MS, (iii) by contrast, several coupled physical-biogeochemical simulations run at global scale and forced by the RCP 8.5 scenario are available. However, they provide very different (and even divergent) dynamics for nutrient concentrations in the next century. Moreover, some of them only provide the nutrients concentrations at the sea surface. In this context, and to consider the effects of hydrodynamics alone, the nutrients boundary concentrations at the Gibraltar Strait were kept constant (describing an annual climatology scenario). However, since the physical simulation used to force the biogeochemical model predicts an increase in the surface water inflow at Gibraltar straits of 0.03 Sv (Darmaraki et al., 2019), corresponding to an increase of 3.7% in the surface inflow over the next century (Soto-Navarro et al., 2020), this increase in the Gibraltar inflow implies an increase in the inflow of nutrients in the surface layer at the Gibraltar Strait and acts against nutrients decrease, and probably mitigates it in our simulation.

4.2. Analysis of the Scenario Results

The simulation shows a significant decrease in surface nutrients over the second half of the twenty-first century in the MS. Between the MID (2031–2050) and END (2081–2100) periods, the model shows a decrease of 48.3 and 37.6% of NO₃ and 38 and 37% of PO₄ (Figure 4) in the Western and the Eastern Basins, respectively. As already mentioned, these changes in nutrient concentrations and in the other state variables of the biogeochemical model can be ascribed to the effects of variations in hydrodynamics induced by climate change, independently of other sources of variability. Furthermore, considering that our model does not include the effect of temperature on the kinetics of the biogeochemical processes, and that currents are expected to remain almost unchanged (except locally in some areas), the main factor that could explain the decrease in nutrient concentrations is the vertical mixing and the associated MLD. The MLD is expected to progressively become shallower, mostly after 2050, as shown in Figure 3 and by Soto-Navarro et al. (2020). This decrease is directly linked to the increase of the stratification shown in Figure 3. The changes in water temperature and salinity over the next century induces a stronger stratification. Due to that stronger stratification of the water column, more energy is needed in winter to break it, and therefore, the winter mixed layer is shallower. This strong decrease in the winter MLD (nearly 20% in the Western Basin between HIST and END and 17% between MID and END), is probably responsible for the decline in surface nutrients after 2050 that mostly impacts the Western Basin (Figure 4). The decline in surface nutrients over the next century caused by an enhanced stratification has already been predicted in most global scale studies for the sub-tropical areas (Bopp et al., 2005, 2013; Marinov et al., 2010; Dutkiewicz et al., 2013, 2019; Bindoff et al., 2019). The same conclusion has been drawn for the MS by Lazzari et al. (2014) (SRES-A1B) or for the north-Western MS by Herrmann et al. (2014) (SRES-A2), though they both use discontinuous simulations and different climate scenarios. This conclusion is also that suggested by Macias et al. (2015) using a continuous simulation under RCP 8.5 climate condition and constant river inputs, but only for the Western Basin, whereas they predict an increase in nutrient concentrations in the Eastern Basin (due to the intensification of vertical mixing because of the increase in evaporation). This result is not in contradiction with our work since the nutrient decrease in the Eastern Basin is significantly weaker than that in the Western Basin. The study by Richon et al. (2019) obtained with a SRES-A2 scenario indicates that nitrate will be accumulating in the surface layers of the MS over the twenty-first century, while there is no trend in surface phosphate concentrations. However, according to these authors, this is the result of the effects of river and Gibraltar inputs since the reduction of vertical mixing makes surface water more sensitive to external nutrient inputs.

The decrease in nutrient concentrations in the surface layer is concomitant with the deepening of the top of the Ncline and the top of the Pcline, as shown by the model in both basins in summer [the depth of the top of the Ncline increases by 6.9% (5%) and that of the Pcline by 4.5% (8%) in the Western Basin (Eastern Basin)]. The deepening of nutriclines results in less nutrients being brought to the surface layer during winter, which results in the oligotrophication of the surface layer. In the model, the top of the nutriclines under stratified conditions is driven by the equilibrium between the production flux (mineralization through bacterial activity), and the consumption flux of nutrients. Those processes are extensively explained in Pagès et al. (2020). Since nutrient quantities brought to the upper layer during vertical mixing progressively decrease, phytoplankton has to go deeper to find available nutrients in summer to sustain its growth and it progressively erodes the top of the nutriclines. The deepening of phytoplankton during summer is indeed predicted by our simulation through the deepening of the DCM (Supplementary Figures 1C,D). The depth of the DCM increases by 7.3% in the Western Basin and by 5.4% in the Eastern Basin, following the same trend as nutricline depths. A slight deepening of the DCM has already been shown at regional scale by Richon et al. (2019), but Macias et al. (2015) predict a small reduction of the mean depth of the DCMs (mean reduction of 1.2 m over the MS).

Another effect of nutrient decline in the surface layer is the reduction of the surface CHL concentration after 2050 (Supplementary Figure 1). Our simulation do indeed predict a 31% decrease in surface CHL in the Western Basin and 21% in the Eastern Basin between HIST and END, however most of this decrease occurs between MID and END with a reduction of 25% for the Western Basin and of 16% for the Eastern Basin. This is in accordance with the study of Dutkiewicz et al. (2019) at global scale which shows a reduction of CHL in several regions of the ocean after 2030 (0.5% a⁻¹ for the North Atlantic). Richon et al. (2019) also highlighted an almost 9% decline in the integrated CHL over the MS (while this study shows a 16% decline in the vertically-integrated CHL in the Western Basin and of 3.4% in the Eastern Basin). Finally, the study by Herrmann et al. (2014) shows an opposite trend with a 8% increase in the integrated CHL over the north-western MS. However, primary production is positively correlated with temperature in their model through an exponential law (Eppley, 1972), which increases the photosynthesis rate as temperature increases.

The aforementioned drop in surface CHL concentrations is related to the reduction of PHYL abundance (${\text{PHYL}}_{int}^{cell}$) (Supplementary Figure 2) and carbon biomass (${\text{PHYL}}_{int}^c$) simulated mostly after 2050 (not shown), as a consequence of the increasing oligotrophy. By contrast, PHYS which is better adapted to oligotrophic conditions, slightly increases in the Eastern Basin, but not enough to compensate the lack of CHL due to PHYL decline (PHYL is indeed responsible for a large part of the CHL production, mostly during the spring bloom period). The study by Flombaum et al. (2013) also predicts at global scale, that small phytoplankton (Prochlorococcus and Synechococcus) abundance will be enhanced (respectively by 29 and 14%) in the next century. Our study is also in agreement with the study of Herrmann et al. (2014) who showed that the planktonic community structure should evolve toward larger biomasses of small-size groups [picophytoplankton (+25%), bacteria (+18%) and nanozooplankton (+17%)]. Richon et al. (2019) showed that diatoms were expected to be more sensitive to climate change as their simulated biomass decreased more sharply than that of nanophytoplankton. Overall, it has been shown by most studies at global scale that large phytoplankton are more disadvantaged when external conditions become more nutrient-depleted (Bopp et al., 2005; Marinov et al., 2010; Dutkiewicz et al., 2013).

As a result of the drop in ${\text{PHYL}}_{int}^{cell}$ and in surface and DCM CHL (Supplementary Figure 1), the integrated primary-production decreases by 12.7 and 14%, respectively for NPP_int and GPP_int in the Western Basin and by 9.3 and 11%, respectively for NPP_int and GPP_int in the Eastern Basin between MID and END periods. The historical value of GPP_int, namely 63 and 43 gC·m⁻²·a⁻¹, respectively for the Western and the Eastern basins are in good agreement with the literature, although at the lower end of the range. For example if we focus on modeling studies, Guyennon et al. (2015) calculate a GPP of 82 gC·m⁻²·a⁻¹ for the whole Mediterranean basin, while Lazzari et al. (2012) obtain an average of 98 gC·m⁻²·a⁻¹ and Crispi et al. (2002) of 88 gC·m⁻²·a⁻¹. The satellite values are more sparse from 68 gC·m⁻²·a⁻¹ (Uitz et al., 2012) to 136 gC·m⁻²·a⁻¹ (Bosc et al., 2004). Regarding in situ values, Siokou-Frangou et al. (2010) show ranges between 59 and 150 gC·m⁻²·a⁻¹ and according to Ryther (1969), oligotrophy is characterized by a primary production smaller than 50 gC·m⁻²·a⁻¹.

Most of the climate scenarios at global scale suggest a decrease in NPP, though in some specific regions, an increase has been predicted (Bopp et al., 2005, 2013; Dutkiewicz et al., 2013, 2019; Yool et al., 2013). Bopp et al. (2013) and Yool et al. (2013) have predicted a 8.6 and 6.3% reduction of the NPP_int at global scale by the end of the century, respectively. Finally, Bindoff et al. (2019) show that tropical NPP will decline by 7–16% by 2100. At a regional scale, the conclusions of previous studies are more sparse. Lazzari et al. (2014) did not predict any trend in the NPP (but with a discontinuous scenario only spanning the 1990–2000 and 2090–2100 periods), Macias et al. (2015) pointed out a decrease of the primary production rate in the Western Basin and an increase in the Eastern Basin. Richon et al. (2019) show a 10% reduction of the NPP_int. Finally, some studies also predict an increase in primary production, with Herrmann et al. (2014) showing a 19% increase in the Gross Primary Production (GPP) for the north-western MS and Moullec et al. (2019) reporting an increase in the NPP stronger in the Eastern Basin than in the Western Basin.

The decline in ${\text{PHYL}}_{int}^{cell}$ and ${\text{PHYL}}_{int}^c$ also directly impacts its main predator in the model, namely the mesozooplankton (Supplementary Figure 4). Our model shows that ${\text{COP}}_{int}^{cell}$ (${\text{COP}}_{int}^c$) decreases by 7% in the Western Basin and 5% in the Eastern Basin. This result is consistent with the fact that large organisms are disadvantaged in oligotrophic areas. In addition, Stock et al. (2014) show a 7.9% reduction of mesozooplankton using a RCP 8.5 scenario. In our model, the fact that ${\text{COP}}_{int}^{cell}$ decreases, may have potentially significant effects on small pelagic fishes such as sardines and anchovy, as pointed out by Woodworth-Jefcoats et al. (2017) for the tropical regions and suggested in Pagès et al. (2020) for the MS.

The decrease in ${\text{PHYL}}_{int}^{cell}$ and ${\text{COP}}_{int}^{cell}$ also induces a 16% (14.7%) and 15.6% (13.2%) reduction of DETC_int in the Western Basin and Eastern Basin, respectively between HIST and END (MID and END). In the model, the DETC compartment is only fueled by PHYL and COP through the processes of mortality or fecal pellets production. The simulated dissolved organic carbon (DOC) did not exhibit any trend despite the increase in ${\text{PHYS}}_{int}^{cell}$. The decrease in DETC_int directly affects the export of organic matter at 1000 m (${\text{EXP}}_{1000m}^{org}$), which decreases by 22% (18%) and 20% (15%) in the Western Basin and the Eastern Basin, respectively (Supplementary Figure 5) between HIST and END (MID and END). At global scale, Bopp et al. (2001) show a 6% reduction of carbon export at 100 m, while Dutkiewicz et al. (2013) show a 10% decrease in organic export below 100 m. Bopp et al. (2013) also highlight a decrease in carbon export at 100 m at global scale due to reduced primary production combined with a change in the plankton community in favor of small organisms. Recently, it has been shown for the global ocean by Bindoff et al. (2019) that the downward carbon flux at 1000 m should decrease by 9 to 16% by 2100 under a RCP 8.5 scenario. At regional scale however, Herrmann et al. (2014) and Moullec et al. (2019) show that the organic export remains stable or even increases.

4.3. The Strong Decrease in NO₃ Concentration in the Western Basin and Its Effect on Plankton Community

In the previous section, we highlighted that nutrient concentrations in the surface layers are expected to decline in both basins. However, in the Western Basin, the drop is higher for NO₃ (50%) than PO₄ (43%). By contrast, in the Eastern Basin, the decrease in PO₄ (54%) is stronger than the decline in NO₃ (42%) between HIST and END. Since the MLD is, in the present simulation, the main factor underlying the predicted trends in nutrient concentrations (see section 4.2), we can consider that the maximum depth of the MLD reached during winter drives the quantity of nutrients that will be available in the surface layer to sustain the biological production for the rest of the year. This quantity of nutrients is also strongly linked to the shape of nutrients vertical profiles in the water column. When looking at the profiles shape for the HIST, MID and END periods in the Western Basin (Figure 6), the slope of NO₃ profiles is weaker than that of PO₄, which means that NO₃ concentrations decline more rapidly with depth. For the HIST period, the average of the maximum MLD values is 84 m (see section 3 and Table 1), corresponding to values of 3.42 and 0.2 μmol·l⁻¹, for NO₃ and PO₄, respectively, that can be achieved by the vertical mixing. During the MID period, the values drop, respectively to 2.98 and 0.17 μmol·l⁻¹ for NO₃ and PO₄, which means a decrease of 20% for both nutrients. For the END period, the values are now equal to 1.11 (0.08) μmol·l⁻¹ for NO₃ (PO₄), corresponding to a 62.7% (52.9%) decline between MID and END periods. In other words, the slight difference in the shapes of the Ncline and the Pcline that arose between the MID and the END periods led to a 9% deficit in NO₃ inputs by winter vertical mixing as compared to those of PO₄. This provides a possible explanation for the stronger deepening of the Ncline (5.3%) than the Pcline (2%) in the Western Basin shown by the model. Since 2050, the NO₃:PO₄ ratio at the winter MLD depth regularly decreases (not shown). This leads to a drop in the N:P ratio in DET and DOM (not shown) in the upper layer of the Western Basin. Due to that decrease, the available amount of N for mineralization is lower than the amount of P at depth. This implies a stronger deepening of the Ncline than the Pcline in the Western Basin which progressively increases the shift between the NO₃ and PO₄ surface concentrations.

If we focus now on the Eastern Basin (not shown), between HIST and END the decrease of nutrients concentration at the MLD depth reaches quite similar values for NO₃ and PO₄, namely 35 and 38%. Those similar values are due to the MLD being deeper in the Eastern Basin than in the Western Basin and to the fact that both nutriclines at this depth are steeper and their slopes are closer.

To explain the 12% gap in the surface layer between PO₄ and NO₃, we need to focus on the rivers' nutrient inputs. In the simulation, in the Eastern Basin, PO₄ inputs strongly decrease between 1985 and 2000 (Ludwig et al., 2009). This induces a progressive deepening of the phosphocline (while the nitracline remains at the same depth) in the Eastern Basin, as shown in a previous paper (Pagès et al., 2020). Due to that deepening of the phosphocline the amount of PO₄ brought to the surface by winter vertical mixing decreases. That leads to a strong PO₄ reduction in the surface layer lasting in the model at least until 2010 despite the fact that rivers nutrients inputs after the year 2000 are set constant in this simulation. This reduction in rivers' nutrient inputs could explain why the Eastern Basin shows a stronger decrease in PO4 than in NO₃ between HIST and END. For the Western Basin, the changes in rivers' nutrient inputs during the 1985–2000 period are far less significant and do not impact the nutriclines depth (Pagès et al., 2020).

Since NO₃ concentration declines more than that of PO₄ in the surface layer of the Western Basin, changes in factors controlling plankton growth during summer occur. A significant part of the 0–150 m depth area where PHYS growth was mostly limited by P, progressively becomes limited by N or co-limited by N and P in the Western Basin (Figure 7). We observe a shift from a 30% P-limited and 18% N-limited in the Western Basin for the HIST period, to a 17% P-limited and 42% N-limited in the Western Basin by the end of the century. By contrast, there is no significant change in limiting nutrients in the Eastern Basin since NO₃ and PO₄ decline in a similar way in this basin.

For PHYL, we observe the same trend as for PHYS toward an expansion of the N-limitation in the Western Basin.

5. Conclusions

This study offers new insights regarding the possible effects of variations in hydrodynamic forcing induced by climate change on the biogeochemistry of the Mediterranean Sea and enhances the findings of other studies. It is the first RCP 8.5 continuous scenario run for the MS between 1985 and 2100 with a fully flexible-stoichiometry model, here the Eco3M-MED model having the distinctive characteristics of considering the abundance of the different PFTs and not directly taking into account the temperature effects on the biogeochemical processes. Moreover, the physical forcing used here is innovative because it is derived from a coupled model (CNRM-RCSM4) with a high degree of consistency between oceanic and atmospheric forcing.

This study highlights that most of the changes in the biogeochemistry of the MS could occur (under the RCP 8.5 scenario) after 2050. Before that, the changes remain sparse for most variables. The model clearly shows a strong decline in surface nutrients throughout the MS, which is in good accordance with regional and global studies. To our knowledge, this study is the first to show a discrepancy between the respective decrease of NO₃ and of PO₄ surface concentrations in the Western Basin after 2050 and to provide an explanation for it: because the slope of the average autumn NO₃ profile is flatter than the slope of the average autumn PO₄ profile, when the MLD becomes deeper during winter the amount of NO₃ brought into the upper layers by vertical mixing is weaker than the amount of PO₄. This discrepancy induces significant changes in the Western Basin, such as a deepening of the Ncline under the Pcline as well as significant changes in the limiting nutrients for phytoplankton and bacteria.

We also show in this work a strong decline in surface CHL concentrations (as well as at the DCM depth), in integrated net primary production and in biomass and abundance of large organisms (PHYL and COP) that could have a strong impact on the small pelagic fish community.

We also pointed out a strong decrease in organic carbon exports at 1000 m in both basins mostly due to the reduction of the abundance and biomass of PHYL and COP. We show in addition that variations in hydrodynamic forcing induced by climate change alone under the RCP 8.5 emission scenario will lead to a strong increase in the oligotrophy of the MS, and that the Western Basin seems to be more impacted than the Eastern Basin.

Finally, our study shows that the projected effects of climate-induced changes in hydrodynamics on the biogeochemistry of the Mediterranean Sea under a RCP 8.5 regional climate scenario are far from negligible.

As such, this scenario was highly instructive and has enabled us to disentangle the effects of hydrodynamics alone. These changes could be even more important if nutrient inputs by rivers and at the Gibraltar Strait were to decrease, or mitigated if they were to increase. The relative contributions of each of the internal and external sources of variability should be calculated in the future when all the forcing scenarios would be available.

Further work is necessary to provide a consistent prediction of the biogeochemistry of the MS by the end of the century. This study is indeed intended to be a first step toward an integrated prediction of land and sea responses to global change in the Mediterranean Basin using an integrated process-based model which will allow the use of consistent inputs of nutrients by rivers.

To conclude, this work and the forthcoming one will not be enough to explore all the different sources of uncertainties that could impact future projections (choice of scenarios, choice of models, internal variability, external nutrients inputs) and for this purpose, coordinated studies should be carried out in the future, as already done for physical projections.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

This study has been made possible thanks to a Ph.D. thesis funded by the LaSeR-Med project. This project was part of the Labex OT-Med (no. ANR-11-LABX-0061) funded by the French Government Investissements d'Avenir program of the French National Research Agency (ANR) through the A*MIDEX project (no ANR-11-IDEX-0001-02). This work was performed using HPC resources from GENCI-IDRIS (Grant 2018-0227). This work was also a contribution to the MISTRALS/MerMex program. The physical simulation performed with CNRM-RCSM4 has been carried out on the Meteo-France computing facilities. It was a contribution to the Med-CORDEX multi-model initiative (Ruti et al., 2016) and can be downloaded from the Med-CORDEX database (www.medcordex.eu).

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.

Supplementary Material

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

Supplementary Figure 1. Annual average of surface chlorophyll concentration μg.l⁻¹ (blue line on the top) for the Western Basin (A) and the Eastern Basin (B). The black dotted line shows the control simulation. The summer (July to September) average of the deep chlorophyll maximum depth (green) and of the chlorophyll concentration at the Deep Chlorophyll Maximum depth (red) are also plotted for the Western Basin (C) and the Eastern Basin (D).

Supplementary Figure 2. Map of the percentage of changes of the average vertically-integrated (surface to bottom) abundance (left for PHYS and right for PHYL) between the MID period (2031–2050) and the END period (2081–2100).

Supplementary Figure 3. Vertically-integrated net primary production averaged over the Western Basin (A) and the Eastern Basin (B) and vertically-integrated gross primary production averaged over the Western Basin (C) and the Eastern Basin (D) in gC·m⁻²·a⁻¹. The black dotted line shows the control simulation.

Supplementary Figure 4. Annual mean of the vertically-integrated (surface to bottom) biomass in mmolC·m⁻² (green) and abundance in cell·m⁻² (blue) of Copepods averaged over the Western Basin (A,C) and the Eastern Basin (B,D). The black dotted line shows the control simulation.

Supplementary Figure 5. Annual mean of particular and dissolved organic carbon export in gC·m⁻²·a⁻¹ at 1000 m averaged over the Western Basin (A) and the Eastern Basin (B). The black dotted line shows the control simulation.

References