Modeling Carbon Budgets and Acidification in the Mediterranean Sea Ecosystem Under Contemporary and Future Climate

We simulate and analyze the effects of a high CO2 emission scenario on the Mediterranean Sea biogeochemical state at the end of the XXI century, with a focus on carbon cycling, budgets and fluxes, within and between the Mediterranean sub-basins, and on ocean acidification. As a result of the overall warming of surface water and exchanges at the boundaries, the model results project an increment in both the plankton primary production and the system total respiration. However, productivity increases less than respiration, so these changes yield to a decreament in the concentrations of total living carbon, chlorophyll, particulate organic carbon and oxygen in the epipelagic layer, and to an increment in the DIC pool all over the basin. In terms of mass budgets, the large increment in the dissolution of atmospheric CO2 results in an increment of most carbon fluxes, including the horizontal exchanges between eastern and western sub-basins, in a reduction of the organic carbon component, and in an increament of the inorganic one. The eastern sub-basin accumulates more than 85% of the absorbed atmospheric CO2. A clear ocean acidification signal is observed all over the basin, quantitatively similar to those projected in most oceans, and well detectable also down to the mesopelagic and bathypelagic layers.


INTRODUCTION
The Mediterranean Sea (MS) is the largest semi-enclosed European sea, and one of the most valued seas in the world. While mainly oligotrophic, it is rich in biological diversity, one of the world biodiversity hotspots, and comprises a variety of coastal and marine ecosystems that deliver benefits to all its coastal countries (Myers et al., 2000;Cuttelod et al., 2009). Currently, anthropogenic pressures on MS ecosystem are increasing because of tourism (about a third of international tourist arrivals; UNWTO-World Tourism Organization, 2015) and population growth, economic growth of both European and North African countries, sea-based activities such as fisheries, aquaculture, and ship traffic (Micheli et al., 2013). Furthermore, the Mediterranean region has been recognized as a "climate change hotspot" (Giorgi, 2006;Darmaraki et al., 2019b), i.e., an area particularly affected by climate change. Therefore, it is important to assess to which extent climate change will affect the dynamics of Mediterranea marine ecosystems (Lejeusne et al., 2010), and-eventually-its capability to provide ecosystem services (Costanza et al., 1997).
Models projections suggest that climate change will cause an increase in MS surface water temperature (Adloff et al., 2015;Darmaraki et al., 2019a;Soto-Navarro et al., 2020) and an intensification of extreme events (Galli et al., 2017), as well as changes in Mediterranean circulation and in vertical mixing/stratification processes, with less dense water formation (at least in the Western basin), a weaker thermohaline circulation and a stronger vertical stratification (Thorpe and Bigg, 2000;Somot et al., 2006;Adloff et al., 2015;Soto-Navarro et al., 2020).
However, climate change impacts on the MS biogeochemical properties (i.e., nutrients concentration and plankton productivity) at the basin-scale are less explored in scientific literature. Lazzari et al. (2014) performed several simulations to assess the impact of a family of Special Report on Emissions Scenarios (SRES) A1B climate scenarios, in conjunction with different scenarios of land use and related nutrient loads in rivers nutrients, suggesting that in all scenarios the increment in temperature will cause an increase in both plankton productivity and respiration, with the net effect of reducing plankton biomass and the amount of biomass available to higher trophic level organisms (with possible consequences on fishery compartments) and of increasing the labile dissolved organic matter. In that study, however, biogeochemical simulations described only the first and last 10 years of the 21st century, without a continuous simulation of the whole period, potentially introducing a significant bias and unknown uncertainties in future projections. Macias et al. (2015) projected a small increment in primary productivity and significant changes in spatial distribution of plankton biomass under both RCP 4.5 and 8.5 warming scenarios, likely related to changes in the stratification regime of the basin. However, they used a model not validated against experimental observations, and did not discuss possible existence of biogeochemical biases and /or numerical drifts. Richon et al. (2019) explored the impact of SRES A2 scenarios with a model setup similar to the one used in the present study, suggesting that the unbalanced nitrogento-phosphorus input from riverine discharge and fluxes via the Strait of Gibraltar used in their scenario lead to an expansion of phosphorus-limited regions across the MS and a significant reduction of plankton net productivity. More recently, Pagès et al. (2020) produced a continuous RCP8.5 run, that assumed no effect of temperature on the kinetics of the biogeochemical processes and no changes in nutrient loads from rivers, and showed that the climate change related variations in (only) hydrodynamics suffice in causing changes in biogeochemistry (increased oligotrophication, reduced production, decreased phytoplankton abundance, nitrate decline in the surface layer) but mainly in the second half of the 21st century, and only in the western basin. Other projections refer to simulations performed within recent EU projects (Boero et al., 2016), and in a variety of downscaling studies at the regional (Holt et al., 2012;Herrmann et al., 2014), subregional (Lamon et al., 2014) or local Salon et al., 2008;Canu et al., 2010) scales.
Yet, and as far as we know, no paper quantified the expected changes in ocean acidification in the MS, nor described the expected changes in carbon budgets and fluxes, rather than variables concentrations. Therefore, here we performed a continuous 100 years simulation under a high greenhouse gases emission scenario, in order to assess the potential impact of climate changes related effects on carbonate and biogeochemical systems, with a particular focus on changes in carbon cycling, budgets, and ocean acidification. Furthermore, and notably, the impact of climate change has been assessed by contrasting the climate model simulation against a century-long baseline simulation, used as a control run to assess multidecadal trends not related to climate change, but potentially present in long simulations, and which need to be removed to isolate climate change effects and reduce the uncertainty in climate scenarios impact studies. The use of this methodology, widely recognized as appropriate (Sarmiento et al., 1998;Teutschbein and Seibert, 2012), but rarely implemented in biogeochemical studies, increases the robustness of our conclusion and represents an additional element of novelty in our study.

Model System
The simulations have been performed by using a state-of-theart transport-biogeochemical model, purposely developed and validated for describing plankton productivity and carbon and nutrients biogeochemical cycles in the MS (Lazzari et al., 2012(Lazzari et al., , 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 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 . 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 . DIC exchange at the air-sea is resolved by computing the seawater pH, pCO 2 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(Lazzari et al., , 2016Cossarini 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 , integrated with a data assimilation scheme (Teruzzi et al., 2014(Teruzzi et al., , 2018.

The Scenario Numerical Experiment
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 2 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 2 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 2 concentration regardless of the underpinning socio-economic dynamics (IPCC, 2014;Harrison et al., 2019). The time course of CO 2 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.

Conceptual Framework, Major Processes and Relationships Among Nutrients and Carbon Cycling
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 2 dissolves in sea water through the airsea exchange process (arrow #1 in figure), proportional to the difference between concentration of CO 2 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 2 gradient between atmosphere and surface sea, all factors able to change the CO 2 sea surface concentration indirectly affect this process. The factors acting to reduce sea surface CO 2 will favor CO 2 atmospheric dissolution, those acting to increase it will favor CO 2 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 2 chemically reacts with water to form carbonic acid (H 2 CO 3 ), which partially dissociates in bicarbonate (HCO − 3 ) and bicarbonate (CO = 3 ) 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 2 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.

Reducing Uncertainties in Long-Term Dynamic
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.

Model Corroboration and Dynamics During the Historical Period
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 subbasins (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 3 ) 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 3 . 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 .

Long Term Dynamics
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  For each variable, absolute bias and RMSD are lower than 50% of the spatial standard deviation, except in the cases highlighted by symbols: (*) absolute bias larger than 50%, (**) larger than 100%; ( ∧ ) RMSD larger than 50%, ( ∧∧ ) larger than 100%.
Correlation is significant at p-level of 0.05, except in cases highlighted by ( ).
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, 2 | Mean present-day values for western and eastern basins and variations at the end of the scenario simulation (as in Eq. 1), for the epipelagic, mesopelagic and bathypelagic layers, defined, respectively, as the average between 0 and 250, 250, and 1,000 m, 1,000 and 5,000 m. Each cell reports the percentage ratio between the net variation estimated at the end of the century (2090-2099, numerator), and the present-day reference (2000-2009, denominator).

Nitrate (mmolN/m 3 )
Phosphate ( 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 2 and of the changes in CO 2 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(Zunino et al., , 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 3 /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 2 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)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)) DIC concentration. The increment in atmospheric concentration forces higher CO 2 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 subbasin 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.

Carbon Budgets
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 wellknown 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 2 (+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 2 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.

CONCLUSION
A set of Mediterranean Sea biogeochemical simulations forced by a high greenhouse gas emission scenario projects a basin-scale increment of the nutrient pools in the epipelagic and mesopelagic layers, an increment in surface water temperature, an increase in CO 2 fluxes from the atmosphere, an increase of DIC inflow from Gibraltar Strait. As a result there is an increment in DIC all over the basin and an increase in gross primary production, which, however, is more than compensated by the increase in the system total respiration, so that globally there is a decrement in the concentrations of total living carbon, chlorophyll and oxygen in the epipelagic layer. The POC fraction also decreases, because of the increased respiration in the upper layer and of the reduced vertical transport in terms of sinking, as will -even if to a smaller extent-the DOC fraction. The decrease in the concentration of total living carbon concentration is even more significant in the mesopelagic layer.
Nutrients concentration, primary production and respiration all increase more in the western sub-basin than in the eastern one, but the plankton primary production is large enough to deplete the nutrient stocks in the surface layers, giving rise (in this land use scenario) to phosphorus limitation to plankton growth in most areas. sub-basins. Because of the interplay between primary production and respiration, the decrease in the concentration of total living carbon will decrease everywhere, but more in the eastern sub-basin than in the western one. The increment in surface DIC related to the increasing CO 2 atmospheric concentration drives a downward flux of DIC supported by the solubility, physical and biological pumps, but also by the effects of the dense water formation and of the continental shelf pumps acting along the Adriatic and the Agean seas.
In terms of mass budgets, the increase in the DIC compartment at the bottom constitutes the largest contribute in the MS carbon budget at the end of the century, with the eastern sub-basin accumulating more than 85% of the absorbed atmospheric CO 2 , while the major changes in the exchanges through the boundaries occur at the air-water surface, where the increase in atmospheric CO 2 forces a large increase in CO 2 dissolution in surface water. The horizontal carbon exchange between the two sub-basins also increases.
Finally, even if alkalinity in the MS is higher than in the ocean, a clear ocean acidification signal is observed all over the basin. In fact, the acidification rate is similar to those observed in most oceans, the major processes acting to change alkalinity level also cause a corresponding change in DIC concentration, and the two processes have opposing effects on pH, while CO 2 invasion affects only DIC. The acidification signal is well detectable down to the mesopelagic and bathypelagic layers, and potentially able to impact Mediterranean surface and deepsea ecosystems.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
CS and SSa contributed to the conception and design of the study and wrote the first draft of the manuscript with the contributions from all authors. PL, SSo, GC, GG, and GB organized the database and performed the simulation. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This study has been co-sponsored by the MEDSEA EU project, the COCONET EU project, and the grant "Cambiamenti Climatici in Friuli Venezia Giulia" of the Regione Friuli Venezia Giulia (Italy).

ACKNOWLEDGMENTS
GG was supported by OGS and CINECA under HPC-TRES program award number 2015-05. This study has been conducted using data from the OGS National Oceanographic Data Center (NODC-OGS) and information from the E. U. Copernicus Marine Service. The physical NEMOMED8 simulation has been run on the Meteo-France supercomputer center (Toulouse, France) as a contribution to the Med-CORDEX initiative (www.medcordex.eu). We would like to thank F. Sevault for running the NEMOMED8 simulation and sharing the forcing dataset.