Author responses to Reviewer#1 comments for the manuscript: “Acidification, deoxygenation, nutrient and biomasses decline in a warming Mediterranean Sea”

The authors projected the climate change-related impacts in the marine ecosystems of the Mediterranean Sea in the middle and at the end of the 21st century using an offline coupling model combining the physical model MFS16 and the transport-reaction model OGSTM-BFM, under emission scenarios RCP4.5 and RCP8.5, focusing on the middle and the end of 21 st century. Projected changes are presented for temperature, salinity, dissolved nutrients and oxygen, net primary production, respiration, organic matter, plankton and bacterial biomass, particulate organic matter, and biogeochemical parameters (DIC, pH). The paper provides interesting projections in a changing Mediterranean Sea that is already under multiple pressures. We thank the Reviewer #1 for their positive feedback and for providing detailed comments and suggestions, which will be considered to improve the manuscript. Reviewer’s comments are in bold, authors' responses are in normal font, italicized where they quote the proposed changes to the manuscript. phytoplankton biomass and net primary production by the of the in the they are sub-basins Mediterranean Sea. Abstract:

We thank the Reviewer for raising this point that was not thoroughly discussed in the manuscript. Setting the biogeochemical BCs at the Gibraltar Strait is not straightforward because no-verifiable and anomalous simulated tendencies outside the Strait can result in spurious trends within the basin.
We considered the temporal evolution of the carbonate system and dissolved nutrients concentration in the Atlantic area outside the Gibraltar Strait during the 21st century in the CMCC-CESM global simulation (Vichi et al., 2011). However, after a preliminary analysis, the projected changes simulated at the Gibraltar Strait for nutrients have not been used as boundary conditions due to the anomalous nitrate/phosphate ratio observed. Because of that we prefer to keep the profile corresponding to the present conditions. Moreover, the analysis of the temporal evolution during the 21st century of DIC and total alkalinity in the Atlantic area outside the Gibraltar Strait, showed a continuous increase in the DIC concentration (more than 4% with respect to the present) and stable value of total alkalinity (the variation is around zero in the 21st century). We considered these signals reasonable and, because of that, we kept the CMCC-CESM simulated vertical profiles of DIC at the Gibraltar Strait while the total alkalinity was kept constant throughout the simulation. We acknowledge that we should have used two different BCs for the BGC scenarios. However, besides the lack of availability of a consistent RCP4.5 (w.r.t. the used RCP8.5) global scenario, we decided to test the impact of the different scenario atmospheric forcings (i.e., atmospheric pCO2 and physical dynamics) keeping the BC constant. Thus, our analysis revealed the impact of different atmospheric forcing but, possibly, could have underestimated a potential difference between the two scenarios. The manuscript will be modified taking into account the Reviewer's concern. In particular we propose to write in section 2.3:

"Boundary conditions are adopted to represent the external supply of biogeochemical tracers and properties from the Gibraltar Strait and the Mediterranean rivers into the basin. The exchanges of nutrients and other biogeochemical tracers
and properties in the Gibraltar Strait are achieved by relaxing the 3D fields in the Atlantic zone ( Fig. 1) to average vertical profiles which, for dissolved oxygen, phosphate, nitrate and silicate, refer to Salon et al. (2019), while total alkalinity is based on what was described in Cossarini et al. (2015). These profiles do not consider a seasonal cycle or a future temporal evolution, with DIC as the only exception, which is prescribed from a global ocean-climate simulation under RCP8.5 emission scenario performed within the framework of the CMIP5 project (Coupled Model Intercomparison Project Phase 5;Taylor et al., 2012) and based on the CMCC-CESM modeling system (CMCC-Coupled Earth System Model;Vichi et al., 2011). The reasons for these choices rely on: (i) anomalous values observed in N:P ratio under the RCP8.5 emission scenario, (ii) negligible variation, under emission scenario RCP8.5, of the total alkalinity along the 21 st century, (iii) lack of a consistent RCP4.5 scenario, (iv) the possibility, using the same conditions at the Atlantic boundary, to test the impacts of the different atmospheric and ocean forcings." Moreover, in the Discussion and Conclusions section we will acknowledge that: "Moreover, the use of the same Atlantic boundary conditions for the two scenarios (section 2.3) could have led to an underestimation of a potential difference between the two scenarios in the areas most influenced by the boundary (e.g.,

Alboran and Southern Western Mediterranean)"
Finally, we acknowledge that the Gibraltar strait is an interesting area to be investigated. However, the paper provides evidence of basin wide signals without focusing on specific areas due to the paper length limit and some plots of evolution on some specific areas, including Alboran Sea, are already in Fig. S4-7. If not further required by the Reviewer, we would prefer not to include other figures and tables due to the already high amount of materials discussed in the manuscript and just to direct the potential Readers to the provided references. We thank the Reviewer for highlighting this point, since we realize that also this information was missing in the text. We are aware that, in the CMIP5 simulations protocol, there is an historical simulation which covers the period  and then a scenario simulation which spans over the period 2006-2100. The choice to consider the period 2005-2020 as present climate for the validation relies on the idea to use more advanced datasets such as the CMEMS reanalysis (Teruzzi et al., 2019;Cossarini et al., 2021) and satellite Chl-a dataset (Colella et al., 2016) to evaluate the performance of the modeling system. Being that both the datasets cover a period spanning from 1999 to 2020, in order to avoid the overlapping between historical and scenario we limited the evaluation to the period 2005-2020. We agree with the Reviewer that the choice of the period could introduce a potential error in our analysis due to the fact that we are already in "scenario mode". However, we tested the statistical significance of the differences between two scenarios in the period 2005-2020 for temperature, salinity and current speed fields and we found that the differences are not statistically significant over most of the basin (with only exception of small areas in the southern Ionian, Adriatic Sea and Levantine basin) and anyhow lower than the climate change signal. We do not expect that the choice of the period of different length as PRESENT and MID-FUTURE / FAR-FUTURE significantly impacts our findings. However, in order to address the Reviewer's comment, we propose to add the following sentence:

Periods selection PRESENT MID and FAR
"The period 2005-2020 has been chosen as reference (also in the forthcoming validation) due to: (i) the availability, after 2000, of more advanced satellite and assimilated datasets to evaluate the biogeochemistry of the basin, (ii) to avoid the overlapping between historical and scenario part of the simulations (with the latter starting in 2005), (iii) the observed differences during this period between the two scenarios have been found not statistically significant over most of the basin (not shown)."

Statistical significance:
There is no indication in the text about the significance of the differences obtained for comparison between the MID and FAR periods. The word 'significant' is used but without statistics. In order to evaluate whether the numeric differences are substantial, it is necessary to calculate some parameters such as the t-student (See line 517 as an example) and to indicate if the variations are significant in the text and in the figures.
We thank the Reviewer for raising this point. We propose to redraw all the figures following their suggestion. In particular, we propose to assess the statistical significance of the differences in each point of the domain between PRESENT and MID-FUTURE and between PRESENT and FAR-FUTURE using a Mann-Whitney test with p<0.05. Just as an example here we report the proposed "New" Figure 5 with the relative caption.  PRESENT (2005-2020, a,b,c andd), and relative climate change signal (with respect to the PRESENT) in the MID-FUTURE (2040-2059, e,f,g andh) and FAR-FUTURE (2080-2099, i,j,k andl)  Some temporal correlations are indicated over the manuscript, for example between oxygen concentration and decrease of MLD, but any values are given (r and p-value) to assess the results. Please add some statistics.
We thank the Reviewer for the suggestion. In order to accommodate with the Reviewer's comment we propose to reformulate the sentence as follows: "Under RCP8.5 (Fig. 9 a,b), the progressive decline of oxygen concentration is timely corresponding to the progressive decrease in the maximum mixed layer depth (Fig. S4) and weakening of the zonal stream function (Fig.4) discussed in Section 3.2. In particular, in the North-Western Mediterranean the correlation coefficient between the average dissolved oxygen concentration in the first 100 m and the maximum mixed layer depth has been found equal to 0.64 (statistically significant with p<0.05)."

Model validations:
Please add more quantitative validations with for example a Taylor diagram.
Following also the request from the Reviewer#1, we propose to enrich the validation part by including a more quantitative description of the biases. The use of the Taylor diagram was a possibility discussed during the preparation of the manuscript to provide a synthetic summary of validation. However, it was discarded because the standard deviation of the modelled and observed data was quite different below 200m, possibly misleading in this way the main messages of the validation. In fact, we preferred to focus on a more qualitative but explicative comparison where showing the biases is the main objective. Indeed, as also suggested by Reviewer#1, we will revise the section adding more explanations about model performance in reproducing the mean present state.

The model is validated based on 2005-2020 averages made with the CTRL simulation forced with RCP scenario.
As mentioned before, it is not well suited to use that for model validations because the scenarios already have an impact. Could you discuss the potential implications?
As we reported before, the choice to consider the period 2005-2020 as present climate for the validation relies on the idea to use more advanced datasets such as the CMEMS reanalysis (Cossarini et al., 2021, Teruzzi et al., 2019 and satellite chl-a dataset (Colella et al., 2016) to evaluate the performance of the modeling systems and what we wanted to consider as present-day condition. Being that both the datasets cover a period spanning from 1999 to 2020, in order to avoid the overlapping between historical and scenario we limited the evaluation to the period 2005-2020. We agree with the Reviewer that the choice of the period could introduce a potential error in our analysis due to the fact that we are already in "scenario mode". However, we tested the statistical significance of the differences between two scenarios in the period 2005-2020 for temperature, salinity and current speed fields and we found that the differences are not significant over most of the basin (with only exception of small areas in the southern Ionian, Adriatic Sea and Levantine basin) and lower than the climate change signal. Additionally, in the method section, we clearly specify that the climate impacts are evaluated with respect to the 2005-2020 period that we introduced as PRESENT.

Discussion limitations:
The discussion is interesting overall but there is a lack of links between the paragraphs and some parts could be more detailed.
We thank the Reviewer for this observation. The discussion and conclusions will be revised in order to address their concerns as detailed below. We propose to restructure the discussion and conclusions sections, leaving the stage to a more detailed discussion on the comparison between our projected changes and those discussed in the scientific literature. First we will clearly state that:  (Hermann et al., 2014;Lazzari et al., 2014;Macias et al., 2015;Moullec et al., 2019;Richon et al., 2019;Pagès et al., 2020;Kwiatkowski et al., 2020;Solidoro et al., 2021)." Then we propose to discuss separately the projections related to nutrients, net integrated primary production and respiration, deoxygenation and dissolved inorganic carbon and pH. Finally, we will discuss the physical mechanisms driving the projected changes and the related source of uncertainties:

"Our projections for the biogeochemical tracers, properties and processes at the end of the 21st century show several signals (the decrease in dissolved nutrients in the euphotic layer of the basin and in the intermediate layer of the central
part of the Mediterranean Sea, the increases in the net primary production and respiration, the decline of the stocks of particulate carbon biomass, an uniform surface and subsurface deoxygenation of the water column, acidity of the water column; table SP1) that are mostly in agreement with previous studies (Hermann et al., 2014;Lazzari et al., 2014;Macias et al., 2015;Moullec et al., 2019;Richon et al., 2019;Pagès et al., 2020;Kwiatkowski et al., 2020;Solidoro et al., 2021).
The decline in the dissolved nutrients at the surface is comparable with that observed in Richon et al. (2019). However, they project an overall increase in the concentration of both nutrients at the surface after 2050, which is ascribed by the authors to river and Gibraltar inputs that are not constant over time (as in our case) but are based on a global climate scenario simulation. As highlighted by Richon et al. (2019), the sensitivity of the biogeochemical fluxes at the river loads and Gibraltar exchanges is of paramount importance, and surely worthy of further investigation. Nevertheless, the

increase in the concentration of nutrients in the intermediate layers of both the Western Mediterranean and Levantine
basin can be also traced back to the reduced vertical mixing resulting from the increase in the vertical stratification (Somot et al., 2006;Adloff et al., 2015;Richon et al., 2019).
An overall increase in the net primary production and respiration is projected in both scenarios although many recent studies in the Mediterranean region have shown a different response of integrated net primary production to climate change in both Western and Eastern basin (e.g. Macias et al., 2015;Moullec et al., 2019;Pagès et al., 2020). In fact, this response may vary according to the sensitivity of the assumptions (model equations) for primary production and recycling processes to changes in temperature (Moullec et al., 2019). In the BFM model temperature regulates most of the metabolic rates with a Q10 formulation (Vichi et al., 2015). The increase in net primary production is a consequence of such dependence. In other studies (Eco3M-Med model;Pages et al., 2020) (Bindoff et al., 2015).
An increase in the dissolved inorganic carbon content and acidity of the water column is found as observed also in (Solidoro et al., 2022). The overall accumulation of CO2 in the basin resulted in an acidification of the Mediterranean water with a decrease in pH of approximately 0.23 units, which is slightly lower than the 0.3 projected on a global scale (Kwiatkowski et al., 2020) and lower than the value in Goyet et al. (2016), who projected, using thermodynamic equations of the CO2/carbonate system chemical equilibrium in the seawater, a variation of 0.45 pH units in the basin under the worst SRES case scenario (0.25 units in the most optimistic scenario). However, this last estimate probably tends to overestimate the future acidification of the basin, as it does not consider the decrease in the exchanges and the penetration of CO2 across the ocean-atmosphere interface due to the warming of the water column (MedECC, 2020).This difference in the response to climate change between the Western and Eastern basins has been also observed for the dissolved inorganic carbon accumulation and reflects indeed different factors such as the different ventilation and residence time of water masses in the two basins as well as the exchanges in the Gibraltar Strait (e.g. 'Alvarez et al., 2014;Stöven and Tanhua, 2014;Cardin et al., 2015;Hassoun et al., 2019). Results show that the Western basin, while adsorbing greater quantities, accumulates only a half of the atmospheric carbon stored by the Eastern basin (1.85 PgC) because, in the former, the carbon is partly exported to the northern Atlantic Ocean, while in the latter, it is also affected by a more intense reduction of the thermohaline circulation and therefore in the vertical transport processes, the carbon is retained together with the atmospheric CO2 sink. Additionally, in our case, the use of a high resolution for the biogeochemical projections has shown that in many coastal areas the observed acidification is lower by approximately 8% with respect to the open ocean due to damping effects of total alkalinity input from the rivers (not shown here).
The magnitude of the projected changes has been shown to be scenario-dependent with the largest deviations from the present climate state observed in the RCP8.5 (worst-case) scenario. As shown in the previous results sections our simulations, by covering also the RCP4.5 scenario, highlights how an intermediate greenhouse emission scenario produces results that are not simply an average between the present condition and the RCP8.5, but (at least for some variables) something quantitatively different. For example, the temporal evolution of pH (Fig.15) is similar in two scenarios in the first 30 years of the 21st century. Conversely, after 2050, the pH undergoes a substantial decrease under RCP8.5 while it remains almost stable under RCP4.5 with a final projected variation lower than the half with respect to the worst-case scenario. This supports the idea -possibly based on the existence in a system like the Mediterranean Sea of a certain buffer capacity and renewal rate, that the implementation of policies of reducing CO2 emission could be, indeed, effective and could contribute to the foundation of ocean sustainability science and policies.
The decline in many biogeochemical tracers and properties in the euphotic layer begins in the 2030-2035 period, in correspondence to the weakening of the thermohaline circulation in the basin (Fig. 4), and it is particularly marked in the Eastern basin. This shows that the modification of the circulation resulting from future climate scenarios has substantial effects on the biogeochemical properties of the basin. Changes in the thermohaline circulation of the basin also explain the increase in the nutrient concentration in the intermediate layer of the Levantine basin, which is a result of the weakening of the westward transport of nutrients through the Strait of Sicily (Fig.S5)." The use of 1/16 resolution grid is one of the strengths of this study but I think this needs to be more discussed as said in the General comments section. This point has been highlighted in the introduction and is supposed to be the main improvement of this study. However, in the current state of the discussion, it seems that the same conclusions may have been reached with a lower resolution.
We thank the Reviewer for raising this issue that was missing in the text. The result section and discussion will be integrated with a discussion about the importance of high resolution in allowing the detection, for the first time, of spatial gradients in the same sub-basin or between open and coastal areas that involve signs and the statistical significance of the projected changes.
In particular, in the discussion and conclusion we will state: Could you accentuate the discussion around the difference observed between the RCP 8.5 and the RCP 4.5?
The comparison between RCP8.5 and RCP4.5 is shown in different paragraphs of the result sections. However, given that available results in the scientific literature are often focused on the worst emission scenario, we decided to focus our discussion and conclusion mainly on the RCP8.5 scenario. On the other hand, we agree with the Reviewer that a discussion between the two scenarios should be accentuated. Thus, we propose to include the following sentence in the discussion and conclusion section: "The magnitude of the projected changes has been shown to be scenario-dependent with the largest deviations from the present climate observed in the RCP8.5 (worst-case) scenario. As shown in the previous results section, our simulations, by covering also the RCP4.5 scenario, highlight how an intermediate greenhouse emission scenario produces results that are not simply an average between the present condition and the RCP8.5, but (at least for some variables) something quantitatively different. For example, the temporal evolution of pH (Fig.15)  Agreed. The mixed layer depth is already provided as output by the physical model and is diagnosed using the criterion based on a density difference with the surface lower than 0.01 kg.m −3 . In order to accommodate the Reviewer's comment we propose to modify the caption of Fig. S4

Vertical stratification affects the sinking velocity of particles: It is suggested that vertical stratification affects the sinking velocity of particles (lines 544-545). To my knowledge, even if stratification might impact the downward flux of particles, most of the models didn't take that into account. Could you please explain how it is taken into account in BFM?
We thank the Reviewer for raising this point because we realize that the formulation of the sentence was incorrect. We propose to modify the sentence as follows: "The overall increase in the respiration community has as consequence the decrease in the organic stock matter in the water column."

NPP:
The effects of the model equations on the NPP trends could be more discussed. This is an important difference between the models which needs to be addressed. Your model equations assume that the planktonic community will remain the same over the next century without adaptation to warmer temperatures, which is unlikely.

Furthermore, O'Connor et al., (2011) is cited to explain particulate organic matter decrease despite NPP increase.
This paper focuses on terrestrial plants and herbivores which are very different from phytoplankton. Please, add other more appropriate references.
We thank the Reviewer for raising this point. The present configuration of the BFM does not consider any physiological adaptation to average temperature over a certain time window. Given the lack of clear indication on entity and velocity of adaptation to be added to future scenario simulations, performing a sensitivity study, with numerical simulations based on different degrees of adaptation to environmental change, could be an interesting theme to explore in future works.
We report that the paper cited O' Connor et al., (2011) has a theoretical nature so it is of broad scope. Additionally, its results are applied to plankton data (see their Fig. 3 for example) so we think it is appropriate to cite it in the present manuscript.
We propose to add the following sentence: "In the BFM model, temperature regulates most of the metabolic rates with a Q10 formulation (Vichi et al., 2015). The increase in the net primary production is a consequence of such dependence. In other studies (Eco3M-Med model;Pages et al., 2020) organisms are always optimally adapted and no temperature dependence is accounted for in the physiology.

This different parameterization could be connected to the different results in terms of trends; in fact, the scenarios based
on Eco3M-Med model results in a reduction of net primary production. In this case surface nutrient reduction, rather than temperature, affects the net primary production trend producing a decrease. The relative impact of different drivers (nutrient supply versus organism's adaptation to average water temperature) could be explored with dedicated sensitivity experiments. "

Nutrients:
A nutrients peak is obtained with the RCP 4.5 simulation between 2055 and 2075 (figure S7) and described in lines 455-456. How do you explain this peak? At Gibraltar's Strait the nutrient concentrations are supposed to be fixed, so did the surface water flux change? Could you explain why?
In our computational domain, the Western boundary is located at ~8.8°W (see Figure 1), while the fluxes in Figure S7 are computed at 5°W (inner part of the domain with respect to the Gibraltar Strait). Thus, the shown fluxes account for the variability of water masses inflowing through the Gibraltar Strait combined with the changes in the concentration of nutrients determined by the underlying biogeochemical dynamics. These peaks are primarily driven by the variability of water masses entering the system through the lateral boundaries conditions, taken from the CMCC-CM global simulations (further detailed in the reply to the next comment). The associated signal then extends to a larger part of the Western Mediterranean (see Fig. 7), whereas the changes in nutrients concentration is rather small (±0.01 mmol m -3 for PO4 and ±0.1 mmol m -3 for NO3) when compared to their mean values (see Fig. 5 and Fig. 6).
We recognized that the text at lines 455-456 might be misleading, since we refer there to the inflow of nutrients at the Gibraltar Strait while Figure S7 illustrates the fluxes into the Alboran region (as reported in the related caption) for a longitudinal section located at 5°W that is far from the model open boundary.
The manuscript text is proposed to be modified to comply with the content of Fig. S7 as follows: "[...] to an increase in the nutrient inflow into the Alboran Sea (Fig. S7)." and, accordingly, the caption of Fig.S7 will report the exact location of the section " Fig.S7

Timeseries of the yearly total transport (in Mmol year-1) of Phosphate (left panels) and Nitrate (right panels) in the Alboran Sea in the period 2005-2099 (longitudinal section at 5°W), [...]".
To the best of my knowledge, the Atlantic inflow is hypothesized to increase a bit over the next century due to stronger evaporation in the eastern basin that will increase the SSH gradient. Might this strong peak (and the even higher peak in 2095) be related to model instability? Could you also explain why we didn't observe the peak in the RCP 8.5 simulation where the effect of climate change should be even stronger?This is concerning as this nutrient input affects the response of the WMED (see for example your figure 12 where a peak is visible for the organic matter).
The MFS16 simulation under RCP8.5 scenario projects an increase of the water flow into the Mediterranean basin of about 0.015 Sv (computed as the difference between 2080-2099 and 2005-2020), which is close to results obtained by Adloff et al. (2015) under similar A2 scenario conditions (see Table 2 therein). The simulation under RCP4.5 conditions is instead characterized by a smaller increase of the water flow between 2080-2099 and 2005-2020, equal to 0.005 Sv.
The variability of water inflows through the Gibraltar Strait in the RCP4.5 scenario is determined by the entrance of water masses with lower salinity through the lateral boundary conditions extracted from CMCC-CM global simulations. These events originate from episodic 'leakages' of Arctic waters toward the eastern north Atlantic and are a consequence of the CMCC-CM internal climate variability. In particular, the RCP8.5 scenario is characterized by a far stronger increase of temperature which does not determine the occurrence of similar Arctic waters 'leakages'.
The competing changes in temperature and salinity fields of water masses entering the Strait of Gibraltar is the main driver of the resulting changes in water transport, as previously reported also by Somot et al. (2006) and Adloff et al. (2015).
Such variability of water fluxes in the RCP4.5 simulation results in specific events of larger inflow of nutrients which extend over the Western Mediterranean and stimulate the microbial loop, with increased bacterial and organic matter concentrations (See Fig. 14), and to a lesser extent also the growth of phytoplankton and zooplankton (see Fig. 13). As noted in the previous reply, the changes in biogeochemical variables related to these episodic events is rather smaller when compared to their mean values (see previous reply).

DIC and pH:
Line 626: "Disentangle the temperature and DIC contributions …" how did you do this?
We thank the Reviewer for raising this point that is missing in the manuscript. In order to disentangle the temperature and DIC contributions to pCO2 we used the carbonate system solver of the BFM in an offline mode with the 2005-2099 simulated variables of each scenario. The offline calculation was made twice, keeping constant, alternatively, temperature and DIC concentration at their 2005-2020 values. By comparing the offline calculations with the actual scenario results, we computed the relative impact of the two factors. The sum of the impacts of two factors (computed with the offline method) differs from the actual scenario result of less than 0.2% (i.e., the offline calculations do not account for the synergic impact of the changes of both variable) which is reasonably low to consider our estimations of the relative impacts as robust. We propose to modify the sentence as follows: "In order to assess the temperature and DIC contributions to the pCO2 temporal evolution, the carbonate system We thank the Reviewer for the suggestions. Our proposition is to follow this point and to extensively redraw all the figures. We agree with the Reviewer about the high numbers of figures included in the manuscript. However, we think that all of them provide interesting information to a potential Reader about spatial patterns, temporal dynamics and statistical significance of projected changes and thus, if not specifically required by Reviewer, we would prefer to maintain all of them in the manuscript. As an example, here we show the new version of Figure 5 and relative caption:  (2005-2020, a,b,c and d), and relative climate change signal (with respect to the PRESENT) in the MID-FUTURE (2040-2059, e,f,g andh) and FAR-FUTURE (2080-2099, i,j,k andl)  In order to accommodate with the Reviewer's comment, we will state in the manuscript: " Figure 4 shows the temporal evolution of the Mediterranean thermohaline circulation during the 21st century using the zonal overturning stream function (or ZOF;Myers and Haines, 2002;Somot et al., 2006). The ZOF has been computed by the meridional integration from south to north and from the bottom to the top of the water column of the zonal velocity (see Adloff et al., 2015). The domain of the integration is the same as shown in Figure 1 with the exclusion of the Atlantic area, outside the Gibraltar strait." - Figure 16, there is no label for the year, and apparently, you give twice the WMED, d,e,f should be EMED The figure will be corrected - Figure S7, the y-axis label is too stretched, use scientific notation you will have more space.
The figure will be corrected - Figure S8 (S9) the size of both color bars seems to be differents The difference size of two colorbars is a consequence of the different numerical interval considered. Because of that the size of the two colorbar can be reduced but will not be the same. The figure will be redrawn following the Reviewer's comment. The buffer zone is essentially the area outside the Gibraltar strait and is already shown in Figure 1 (and not shown any more in other figures). In order to better address the Reviewer's request, the caption of Figure 1 will be modified as follows: Agreed. As reply to Reviewer#1, we propose to include and discuss the suggested references in the manuscript. In the introduction we will state:

"An assessment of the effects of climate change on the biogeochemistry and marine ecosystem dynamics of the Mediterranean Sea has been considered in a certain number of studies based on different emission scenarios".
Line 76: A1B climate change: please explain Agreed. We meant "emission scenario". The text will be changed accordingly.
Line 83: A2 emission scenario: please, relate to the RCP scenarios used here Agreed. The sentence will be reformulated as follows:  Sea (e.g. Hermann et al., 2008;Moutin and Prieur, 2012;Guyennon et al., 2015;Ramirez-Romero et al., 2020) We apologize with the Reviewer because the different number of vertical levels in the two models could lead to some confusion. The reason for referring to 72 levels in the MFS16 original configuration is that the domain of the physical model covers a wider area outside Gibraltar (i.e., the western boundary is at 20 o W, see Figure 1 of Lovato et al., 2013) and, above all, deeper with respect to our computational domain shown in Fig.1  The text will be modified accordingly Line 138: Same thing as before for OPA this time.
The text will be modified accordingly Line 145: Same for NEMO.
The text will be modified accordingly Line 153: Same for CMCC-CM. I will not continue to list those errors, please be sure that every acronym is defined as it is first used.
We thank the Reviewer for the suggestion. The manuscript will be extensively checked for that.

Line 188: chemical reactions
The text will be modified accordingly

Line 191-192: CaCO3
The text will be modified accordingly Lines 294-276: I think this is an interesting approach but this paragraph is difficult to read. Could you try to make it a bit more clear? I think that bringing the S1 figure here could be helpful for the reader.
We thank the Reviewer. The methodology part will be extensively written modifying the terminology as follows: Finally, the temporal evolution of the climate change signal (CCS) with respect to the present is given by:

CCS(k)SCEN = X_evo(k)SCEN -X_evoSCEN-PRESENT (2)
where X_evoSCEN-PRESENT is the average of X_evo(k)SCEN in the PRESENT. Hereafter, if not differently specified, all the shown time-series will be represented by CCSSCEN.
Horizontal spatial averages are computed considering the sub-basin defined in Fig. 1 forcing, such as riverine loads. Because of that, following the approach already adopted in previous works (Lazzari et al., 2012;2016;Di Biagio et al., 2019;Reale et al., 2020 a,b) they are not considered in the spatial averages related to WMED and EMED.
Temporal averages of the climate change signals are computed over two 20-year periods: 2040-2059, hereafter MID-FUTURE and 2080-2099 (Fig.S1) We agree that moving Figure S1 in the manuscript would improve the readability of this section. However, we are concerned about the already high numbers of figures. Thus, if not insisted by the Reviewer, we would like to maintain S1 in the supplementary materials to avoid the inflation of materials in the manuscript.
Line 297: Specify the periods of the satellite climatology.

Done.
Line 310 -315: The spin-up information needs to be completed and added before in the manuscript.
Please see our replies to Your major comment concerning the spin-up. We thank the Reviewer for the suggestion. The figure will be redrawn following the Reviewer's suggestion. The text will be modified accordingly

(unless I miss it) so it is confusing.
We thank the Reviewer for pointing out this error in the figures. All the figures will be redrawn following the Reviewer's suggestions. All the time series shown are corrected for the CTRL, except for the temperature and salinity where the CTRL run is not available. This is mentioned in the manuscript at the lines 261-262 where it is stated that: "Hereafter, if not differently specified, all the shown time series will be represented by CCSSCEN." Please note that the methodology part has been rewritten (see below). Because of that this information is not included in the captions which are already long enough. The caption also contains the information about the 10-year running mean which has been applied to reduce the noise in the figure. Because of that this filter is discussed only in the caption and not in the main manuscript. We thank the Reviewer because we forgot to include this information in the caption of Figure   7 that will be modified accordingly.
Line 338: Give more citations here.
Done. Adloff et al., (2015) will be included Line 348: remove " Fig. S2 in the supplementary material" by ( Fig. S2 and S3). The "in the supplementary material" is not necessary and you are not using it all the time, so, remove it everywhere in the manuscript.
Done. The text will be modified accordingly.
Line 360: replace "increased freshwater deficit" by "decreasing freshwater discharge" Agreed. The text will be modified accordingly. Figure 4 will be redrawn following the Reviewer's suggestion. Concerning the Reviewer's request, in the present manuscript, the zonal stream function (ZOF) has been computed following the approach described in Adloff et al., (2015) where it is computed through the meridional integration from south to the north of the zonal velocity which is then integrated from the bottom to the top of the water column. This means the domain of integration is the entire Mediterranean Sea shown in Fig.1. The Atlantic Ocean area, outside the Gibraltar strait, has been excluded from the computations.
In order to accommodate with the Reviewer's comment, we will state in the manuscript: " Figure 4 shows the temporal evolution of the Mediterranean thermohaline circulation during the 21st century using the zonal overturning stream function (or ZOF;Myers and Haines, 2002;Somot et al., 2006). The ZOF has been computed by the meridional integration from south to north and from the bottom to the top of the water column of the zonal velocity (see Adloff et al., 2015). The domain of the integration is the same as shown in Figure 1 with the exclusion of the Atlantic area outside the Gibraltar strait." We thank the Reviewer for the important comment about the figures. The captions of all the figures will be rewritten including some important information that were previously missed. Just as example we propose to state: "Phosphate concentration (in mmol m -3 ) in the layers 0-100m and 200-600m in the PRESENT (2005-2020, a,b,c andd), and relative climate change signal (with respect to the PRESENT) in the MID-FUTURE (2040-2059, e,f,g andh) and FAR-FUTURE (2080-2099, i,j,k andl)