Carbonate system variability in the Mediterranean Sea: a modelling study

1 Introduction

It is now recognized that ocean is a significant sink of CO₂ (e.g., Sabine et al., 2004; Gruber et al., 2019), providing a buffer for anthropogenic carbon emissions and regulating global climate (Archer et al., 2009). A consequence of the CO₂ oceanic uptake is the resulting decrease of pH and carbonate ions concentration, commonly known as ocean acidification (Caldeira and Wickett, 2003), which may have implications (e.g., survival, growth, abundance) over a broad range of calcifying marine organisms in the open ocean (e.g., coccoliths, foraminifera) and coastal areas (e.g., corals, mollusks) (e.g., Kroeker et al., 2013; Leung et al., 2022).

The Mediterranean is a semi-enclosed oligotrophic basin, characterized by a well-defined eastward decreasing productivity (e.g., Moutin and Raimbault, 2002; Bosc et al., 2004; Kalaroni et al., 2020b). It is characterized by relatively higher alkalinity, as compared to the global ocean, due to high evaporation and lateral inputs from rivers and the Black Sea (Schneider et al., 2007). Its more alkaline waters, as well as the deep-water formation and fast overturning circulation, favor the absorption of CO₂ and its efficient transfer to deeper layers (Schneider et al., 2010; Alvarez et al., 2014; Palmiéri et al., 2015; Krasakopoulou et al., 2017; Hassoun et al., 2022). Given the important complexity of the Mediterranean ecosystem, long-term in situ measurements and dedicated modelling studies are required, in order to be able to assess and predict CO₂ storage and acidification trends in the basin.

Understanding carbonate system dynamics is fundamental for optimizing observational strategies and shaping climate change policies with more precise predictions. Given the limited spatio-temporal coverage of in situ data, coupled hydrodynamic/biogeochemical models provide useful tools to gain a better understanding on the carbonate system seasonal and inter-annual variability, quantify the role of different processes (e.g., Carroll et al., 2022) and perform future projections under changing climate conditions (e.g., Reale et al., 2022; Solidoro et al., 2022). Over the last decade, a few modelling studies have investigated different aspects of the carbonate system dynamics in the Mediterranean. Cossarini et al. (2015a) provided a basin-scale analysis on the spatio-temporal variability of alkalinity in the Mediterranean, while Cossarini et al. (2015b) investigated the mechanisms of CO₂ sequestration in the N. Adriatic Sea. Cossarini et al. (2021) provided aspects of the carbonate system variability and inter-annual trends, based on the Mediterranean re-analysis over the last two decades (1999-2019), as part of the Copernicus Marine Environment Monitoring Service (CMEMS) model product. Solidoro et al. (2022) investigated the changes on the Mediterranean plankton productivity, carbon pathways and carbonate system, under a high CO₂ emission (SRES-A2) future climate scenario. Reale et al. (2022) used a relatively higher (1/16°) spatial resolution model to provide a detailed assessment of future climate projected changes in the Mediterranean ecosystem (nutrients, oxygen, plankton production/biomass, dissolved inorganic carbon, pH), under two CO₂ emission scenarios (RCP4.5, RCP8.5). Butenschön et al. (2021) investigated the projected acidification in the Mediterranean, under a baseline future scenario (RCP4.5), along with additional alkalization scenarios. Finally, Palmiéri et al. (2015), used a perturbation approach, not explicitly describing biogeochemical processes, to simulate the anthropogenic carbon storage and resulted acidification in the Mediterranean, influenced by air-sea fluxes and exchanges with the Atlantic Ocean.

A major constraint for the successful setup and implementation of carbonate system models is the limited spatio-temporal coverage of in situ data, which are necessary for the model initial and boundary conditions, as well as to validate the simulated output. Even though the amount of available in situ data on carbonate parameters in the Mediterranean has significantly increased over the last decade, previous modelling efforts have not explicitly investigated the dynamics of the carbonate system variability, while model validation was mostly limited to basin-scale spatial variability. In the present study, a basin-scale Mediterranean carbonate system model has been setup, building on the POSEIDON operational biogeochemical model (Kalaroni et al., 2020a; Kalaroni et al., 2020b). The model output from a 13-year (2010-2022) simulation was evaluated against CARIMED (CARbon IN the MEDiterranean Sea, Alvarez et al., 2019) basin-scale dataset and time-series from three sites in the Western (DYFAMED, Villefranche-PointB) and Eastern (HCB) Mediterranean. Three additional sensitivity simulations, de-activating the effect of biology, freshwater and air-sea CO₂ fluxes, were performed in order to gain a better understanding on the carbonate system dynamics and the processes regulating its variability. The simulated CO₂ air-sea flux spatial variability and calculated inter-annual trends are also discussed.

2 Methods

2.1 Model description

2.1.1 Hydrodynamic/biogeochemical model

A three-dimensional coupled hydrodynamic/biogeochemical model is implemented at Mediterranean basin scale (1/10° resolution ~10km), based on the currently operational model within POSEIDON forecasting system (www.poseidon.hcmr.gr; Korres et al., 2007; Kalaroni et al., 2020a, b). The hydrodynamic model is based on the Princeton Ocean Model (POM, Blumberg and Mellor, 1983), which is a three-dimensional, sigma-coordinate and free surface primitive equation model. The model prognostic variables are temperature, salinity, velocity, sea surface height and turbulent kinetic energy. The vertical eddy viscosity and diffusivity coefficients are computed using the Mellor-Yamada 2.5 turbulence closure scheme (Mellor and Yamada, 1982) that solves the equations for turbulent kinetic energy and turbulence macroscale, considering the wind stirring and the stability induced by stratification. Horizontal diffusion is calculated along sigma-levels following a Smagorinsky formulation (Smagorinsky, 1963).

The biogeochemical model is based on the European Regional Seas Ecosystem Model (ERSEM; Baretta et al., 1995) that follows a “functional” group approach, including four phytoplankton groups (diatoms, nanoplankton, picoplankton, dinoflagellates), three zooplankton groups (heterotrophic nanoflagellates, microzooplankton, mesozooplankton) and bacteria. The pelagic model variables include also particulate and dissolved organic matter, along with dissolved inorganic nutrients (nitrate, ammonium, phosphate, silicate). The benthic-pelagic coupling is described by a simple first-order benthic return module, describing the settling of organic detritus and diffusional nutrient fluxes out of the sediment. For more details on the model setup and implementation, the interested reader may refer to Kalaroni et al. (2020a); Kalaroni et al. (2020b), where the model has been extensively validated against available satellite and in situ data.

2.1.2 Carbonate system model

The hydrodynamic/biogeochemical model POM-ERSEM has been coupled to a carbonate system sub-model (Blackford and Gilbert, 2007), which uses the HALTAFALL (Ingri et al., 1967) carbonate speciation code to calculate the CO₂ partial pressure (pCO₂), pH and carbonate system components (H₂CO₃, HCO₃, CO₃) in the water column, using dissolved inorganic carbon (DIC) and total alkalinity (TA) as input. The carbonic acid dissociation constants are based on Hanson and Mehrbach data, as refitted by Dickson and Millero (1987) and dissociation constant for HSO₄ from Dickson (1990). pH is calculated in the sea water scale (SWS). The carbonate system model has been setup in the Mediterranean considering both DIC and TA as model state variables (submitted to advection, diffusion etc.). The seawater pCO₂ is influenced by biogeochemical (i.e., photosynthesis, respiration) and physical (solubility/temperature) processes, as well as by the air-sea CO₂ exchange, based on the pCO₂ difference with the atmosphere. The CO₂ flux across the air-sea interface is calculated as a function of temperature, wind speed and the difference between the seawater (pCO_2water) and atmospheric (pCO_2air) partial pressure of CO₂:

where H is the solubility of CO₂ in seawater, based on Weiss (1974) and k is the gas transfer coefficient, which is calculated following Nightingale et al. (2000) as:

where U₁₀ is the wind speed 10m above water (m/s) and Sc is the temperature depending Schmidt number.

The variability of TA and DIC is evaluated in the model based on the effect from advection/diffusion (phys) and biological (bio) processes:

At the sea surface, these are also affected by freshwater (E-P) fluxes, which result in their concentration (evaporation) or dilution (precipitation), as in the case of salinity. DIC is also modified based on the CO₂ flux (CO₂flux, Equation 1) at the air-sea interface, which can be positive or negative, depending on the pCO_2air and pCO_2water values. With regard to biological processes, TA is increased (decreased), with the same amount of moles, by the uptake (release) of nitrates and phosphates or the release (uptake) of ammonia (Wolf-Gladrow et al., 2007). Similarly, DIC is decreased by the CO₂ uptake from phytoplankton photosynthesis and increased by CO₂ production from plankton community respiration. The effect of calcium carbonate formation/dissolution on TA and DIC is currently not considered, as calcifying organisms are not explicitly described in the model and information on their abundance and distribution in the Mediterranean is still poorly quantified (e.g., Siokou-Frangou et al., 2010).

2.2 Model setup

A 13-year simulation was performed over the 2010-2022 period. The initial conditions for the hydrodynamic (e.g., temperature, salinity, currents) and biogeochemical (e.g., dissolved inorganic nutrients, plankton groups biomass) models were obtained from Kalaroni et al. (2020a, b) simulation. The initial fields for DIC and TA (Figure 1) were obtained from GLODAP.v2.2016b global gridded fields (1° x 1° resolution, available in Copernicus database, Lauvset et al., 2016), which are based on in situ data from 1972-2013 period. Prior to the 13-year simulation, a 2.5-year spin-up (2010-2011 repeated) was performed starting from June, assuming that the initial DIC and TA fields represent spring period, which is close to the yearly mean, considering the DIC/TA seasonal variability. This 2.5-year spin-up was found sufficient for the purposes of the present study, focusing on the near surface carbonate system variability. At the model open boundary, west of the Gibraltar straits, DIC and TA during water inflow were obtained from GLODAP annual climatology fields (Figure 1). Deeper layer DIC and TA are restored to their initial fields, applying a weak relaxation (timescale:1 year) in order to prevent a systematic model drift over the long-term simulation. Average river inputs of TA and DIC for major Mediterranean rivers were obtained from GEMS/GLORI database (Schneider et al., 2007, http://www.gemswater.org) and available data from the literature (Skoulikidis, 1993, for North Aegean rivers). Average TA and DIC concentration of the inflowing Black Sea Water (BSW) was adopted from Copin-Montegut (1993). The concentration of atmospheric CO₂ was assumed to linearly increase from ~375.2 µatm during 2002 to ~417 µatm during 2022 (~2.1 µatm/year), based on in situ data from ICOS-FINOKALIA station (Crete Isl., Greece, Gialesakis et al., 2023). This is quite similar to the trend (2.08 µatm/year) based on in situ data from Plateau Rosa, Italy (Kapsenberg et al., 2017).

Figure 1 Initial fields for (A) total alkalinity (µmol/kg) and (B) dissolved inorganic carbon (µmol/kg), interpolated in the model grid from Glodap.v2.2016b global gridded fields (1° x 1° resolution, available in Copernicus database). DYFAMED, Point-B and HCB stations are indicated with cyan dots, as well as Gulf of Lions and Levantine basin areas, over which model outputs are averaged in Figures 6–8.

The atmospheric forcing was obtained from the POSEIDON operational weather forecast (Papadopoulos and Katsafados, 2009), using the bulk formulae set from Korres and Lascaratos (2003) for the calculation of momentum, heat and freshwater fluxes at the air-sea interface.

Along with the 13-year reference simulation, three additional sensitivity simulations were performed over the 2010-2014 period, de-activating the effect of biology, freshwater and air-sea CO₂ fluxes on DIC and TA (Equations 3, 4), in order to gain a better understanding on the processes regulating the model simulated carbonate system variability. To evaluate the overall effect of different processes (see section 3.3), the difference between the mean (2010-2014) DIC, TA, pCO₂ and air-sea CO₂ flux, obtained from the reference and sensitivity simulations was calculated in two representative areas (Gulf of Lions: increased mixing and primary production, Levantine basin: stratification, low primary production), along with the basin-scale average. A monthly mean (2010-2014) climatology from the reference and sensitivity simulations was also calculated to investigate the seasonal variability in the contribution of different processes in the two areas (Gulf of Lions, Levantine basin) (see section 3.3).

We should note that seasons in the present study are shifted by one month, as compared to meteorological seasons, to account for the thermal inertia of the ocean. Thus, winter, throughout the text, refers to January-March period, spring to April-June etc.

2.3 Validation datasets

The model simulated DIC and TA spatial variability was evaluated against CARIMED dataset, which includes carbonate in situ data (TA, DIC) from several cruises in the Mediterranean from 1976 until 2018 (Alvarez et al., 2019). From CARIMED data (DIC, TA, Temperature, Salinity), pCO₂ data were calculated using CO2SYS software (Lewis and Wallace, 1998), with carbonic acid dissociation constants from Mehrbach et al. (1973), as refitted by Dickson and Millero (1987) and dissociation constant for HSO₄ from Dickson (1990), as recommended by Alvarez et al. (2014).

The model simulated carbonate system temporal variability was also validated against available in situ near surface data from DYFAMED (7.9° E, 43.3° N, Ligurian Sea), Villefranche-PointB (7.3°E, 43.5°N, Ligurian Sea, 85m bottom depth) and Heraklion Coastal Buoy (HCB, 35.43°N, 25.07°E, Crete Isl., 180 m bottom depth) sites. DYFAMED in situ near surface data (DIC, TA, pH) over 2010-2019 period were obtained from DYFAMED database (Marty, 2002; Coppola et al., 2021). Villefranche-PointB in situ surface (1m) data (DIC, TA, pCO₂) over 2010-2019 period were obtained from PANGAEA database (doi:10.1594/PANGAEA.727120, Kapsenberg et al., 2017). HCB data included high-frequency pH (3-hour, December 2020 - July 2022) and pCO₂ (6-hour, December 2020 - November 2022) data, obtained from two sensors, deployed at near-surface (~1.5m) water, as well as monthly surface (~2m) in situ DIC and TA data, from Niskin bottle water samples. Additional pCO₂ data were also derived from pH and TA (bottle) data, using CO2SYS (version 3.0, Pierrot et al., 2021) software (see Frangoulis et al., 2024, this issue, for more details on HCB data).

Model simulated sea surface height and temperature were compared against satellite altimetry (~12 km resolution) and sea surface temperature (~4 km resolution) data, obtained from the European Copernicus data base (http://marine.copernicus.eu/), while simulated Chl-a was compared against GLOBCOLOUR (www.globcolour.info) monthly merged chlorophyll-a data(~4 km resolution). Model simulated sea surface salinity was compared against the CMEMS Mediterranean re-analysis salinity product (http://marine.copernicus.eu/).

pCO₂ at DYFAMED and pH at Villefranche-PointB were calculated from TA, DIC and temperature/salinity in situ data, using CO2SYS software (Lewis and Wallace, 1998), with the carbonic acid dissociation constants mentioned above for CARIMED data. pH calculated at Villefranche-PointB, as well as pH in situ data at HCB and DYFAMED, are all expressed in total scale. HALTAFALL speciation code uses a different pH scale (SWS) and dissociation constants (see section 2.1.2). The introduced mean model bias (model-data) from using different dissociation constants and pH scales, as calculated with CO2SYS software, was found approximately -0.016 for pH at PointB, +16.9 µatm for pCO₂ at DYFAMED and +16.5 µatm for pCO₂ from CARIMED data. pCO₂ at constant temperature (pCO₂@20°C) was calculated, using the Takahashi et al. (2002) empirical formula on pCO₂ solubility temperature dependence, which largely drives the pCO₂ seasonal variation. The variability of this temperature normalized pCO₂ (or non-thermal pCO₂) may, thus, be more directly related with other physical (e.g., mixing, advection, air-sea CO₂flux) and biological processes (see discussion in section 3.3). Moreover, it can be used to distinguish the effect of increasing air pCO₂ from that of increasing water temperature on pCO₂ long-term trends (see discussion in section 3.2).

3 Results

3.1 Model validation

In Figure 2 (see also Supplementary Figure S.1), the simulated mean (2010-2022) surface DIC, TA and pCO₂ fields are compared (same locations and seasons) with CARIMED data (pCO₂ calculated from TA and DIC). The observed spatial variability is reasonably captured (pearson correlation: r~0.6 for DIC, r~0.8 for TA, r~0.5 for pCO2, see Supplementary Figure S.2). This is to a point expected, since the model initial DIC and TA fields from the GLODAP climatology are based on a dataset (1972-2013) that largely overlaps with CARIMED (1976-2018). However, the GLODAP gridded fields still present some deviations from the CARIMED actual data, due to smoothing/extrapolation in areas with scarce data availability and also due to the additional data in CARIMED over the 2013-2018 period. In any case, such a comparison ensures that the model does not have a significant bias/drift, preserving the dynamics in different areas of the Mediterranean (see also Supplementary Figures S.3, S.4). A model deviation (also seen in the GLODAP initial fields) is the overestimation of DIC and TA along the pathway of Atlantic water in the Algerian basin and Sicily channel. Another model deviation, which is also seen in the initial fields and is partly removed in the model simulation, is the underestimation of DIC and TA in the North Western Mediterranean (Gulf of Lions and Ligurian Sea). One may also notice the simulated increased TA in the Northern Adriatic, in agreement with the data, which is due to higher alkalinity riverine inputs. Finally, the simulated TA appears slightly overestimated (~15 µmol/kg on average) in the Eastern Levantine (see also Supplementary Figure S.4) and Ionian basins. This might be partly attributed to a model overestimation of evaporation fluxes, as the adopted basin-average mean (2010-2022) evaporation-precipitation flux (~0.88 m/yr) is slightly higher, as compared to existing estimates (0.7 m/yr, Béranger et al., 2005) and also to a small bias of the initial fields, as compared to the data (see Figures 2A, C; Supplementary Figures S.1B, E). The simulated pCO₂ appears slightly underestimated in the Gulf of Lions, similar with DIC, mentioned above and in the Levantine basin (see also Supplementary Figure S.4), which is due to the overestimation of TA, mentioned above. In the Alboran Sea, an overestimation of pCO₂ may be also seen, which might be due to the relatively lower TA, as compared to DIC (see also Supplementary Figure S.4), although the few DIC/pCO₂ data points in this area may be not sufficient for such a comparison. In any case, as the Alboran Sea is directly influenced by the Atlantic water inflow, model deviations in this area suggest that the adopted boundary conditions for TA and DIC may be slightly biased. This bias may be related to the small amount of in situ data, used to derive the GLODAP annual climatology, not properly resolving the seasonal variability of the Atlantic water TA and DIC.

Figure 2 Difference between mean (2010-2022) simulated surface (A) total alkalinity (µmol/kg, TA), (B) dissolved inorganic carbon (µmol/kg, DIC) and (E) pCO₂ (µatm), against CARIMED data. Difference between GLODAP initial model TA and DIC and CARIMED data are also shown in (C, D). Model values are extracted on the same locations and seasons with the data.

In Figure 3, the simulated carbonate system variables (DIC, TA, pCO₂, pH) are shown against available insitu data at HCB, DYFAMED and Villefranche-PointB stations. The simulated DIC seasonal variability is found in reasonable agreement with insitu data (see also Supplementary Figures S.5, S.6), particularly at DYFAMED and Villefranche-PointB stations. A time-lag (1-2 months) of the simulated DIC winter peak, as compared with the data, may be noticed during some years (e.g. 2012, 2013, 2016, 2019) at Villefranche-PointB station (see discussion in section 4). At HCB station, the simulated DIC is found within the range of observations, which does not exhibit a clear seasonal pattern. The simulated pCO₂ and pH show an overall reasonable agreement with insitu data (see also Supplementary Figures S.5, S.6). At all sites, a generally stronger seasonal variability of pCO₂ and pH may be seen, as compared with DIC, which is largely driven by the temperature dependence of pCO₂ solubility. The simulated pH appears slightly underestimated during autumn and overestimated during late winter-spring, which is mostly noticable at HCB (mean model-data, autumn: -0.016, winter-spring: +0.029) and during some years at Villefranche Point-B station (mean model-data, autumn: -0.012, winter-spring: +0.019). This model deviation of pH at Villefranche PointB is associated with the time lag mentioned above for DIC (see discussion section 4). The opposite pattern may be seen for pCO₂ at Villefranche Point-B (i.e., overestimation during autumn and underestimation during late winter-spring). At HCB station, the simulated pCO₂ appears slightly underestimated, particularly during spring period. The simulated TA presents a similar variability with DIC (pearson correlation DYFAMED and Villefranche PointB: r=0.95, HCB: r=0.79), showing a peak during winter due to increased vertical mixing, while a secondary peak, related with evaporation, is simulated during late summer/autumn. This TA secondary peak is more clearly visible at HCB, where simulated TA appears reasonably close to the observations.

Figure 3 Model simulated (red line) surface (A-C) total alkalinity (µmol/kg, TA), (D-F) dissolved inorganic carbon (µmol/kg, DIC), (G-I) pH and (J-L) water pCO₂ (µatm), at DYFAMED (7.9° E, 43.3° N), Villefranche PointB (7.3° E, 43.5° N) and HCB (25.07° E, 35.43° N) respectively, compared with available insitu data (blue dots). For HCB station blue dots refer to data from the pH and pCO₂ sensors, while green open circles refer to bottle data (see section 2.3).

As insitu data for TA do not show a clear seasonal pattern in all sites, a monthly climatology was derived from all available data (Figure 4), in order to gain a better understanding on the seasonality of different variables. At Villefranche-PointB station, the climatology of TA in situ data is characterized by a double peak during winter and later summer-autumn, which appears similar to the seasonality of salinity, being mostly controlled by mixing and evaporation. A less clear variability may be seen for TA at DYFAMED. The simulated salinity shows a similar pattern with the observed seasonality, peaking at winter and late summer. However, the simulated TA shows a stronger decrease during spring and fails to reproduce the observed late summer-autumn peak (see discussion in section 4). At both Villefranche-Point-B and DYFAMED, a much higher variability is simulated, as compared with the data. The resulting poor model skill in reproducing the observed TA seasonal variability is indicated by the low correlation coefficients at both sites (DYFAMED: r~0.04, PointB: r~0.13, see Supplementary Figure S.6). At HCB station, the derived climatology from just two years is characterized by significant uncertainty, particularly for TA and DIC, showing data gaps for a couple a months. Therefore, the observed seasonality may not be properly resolved, providing only some limited indication. The in situ TA data variability suggest some correlation with salinity, which is particularly noticed during March-June period and a relatively higher TA during August-October period. Such a correlation between TA and salinity has been also suggested from other studies in the Aegean Sea (e.g., Hassoun et al., 2015; González-Dávila et al., 2016). The model TA at HCB shows a double peak (February-March and September-October).

Figure 4 Monthly mean (2010-2020) climatology of model simulated (red line) (A-C) surface salinity, (D-F) total alkalinity (µmol/kg, TA), (G, H, I) dissolved inorganic carbon (µmol/kg, DIC), (J-L) water pCO₂ (µatm) and (M-O) water pCO₂@20°C (µatm), at (A, D, G, J, M) DYFAMED (7.9° E, 43.3° N),(B, E, H, K, N) Villefranche Point-B (7.3° E, 43.5° N) and (C, F, I, L, O) HCB (25.07° E, 35.43° N), compared with available in situ data (blue line) monthly climatology.

The simulated DIC reasonably captures the observed seasonality at DYFAMED and Villefranche Point-B, showing a peak during winter. The DIC variability appears slightly overestimated at Villefranche PointB. At HCB station, observations do not exhibit a clear seasonal pattern, but the simulated DIC is found within the observed range. One may notice that, in contrast with the other sites, DIC at HCB station follows a similar seasonality with TA, also showing a secondary peak during late summer/autumn. This peak is related with evaporation, being relatively more significant in the Eastern Mediterranean (see later discussion in section 3.3). The simulated DIC winter peak in all cases is related to increased mixing, which results in the entrainment of deeper water with relatively higher DIC concentration. After this peak, DIC starts to decrease, following its consumption by primary production that shows a peak in the same period (see Supplementary Figure S.7). The DIC decrease during spring period appears overestimated in the model, which is reflected also in the simulated pCO₂ and particularly pCO₂@20°C (see later discussion in section 4). During summer period, there is a small effect from remineralization processes, as net CO₂ production (gross production-community respiration, see Supplementary Figure S.7) becomes sligtly negative (i.e. CO₂ net production), at least at surface. This CO₂ input from biological processes is however partly compensated by CO₂ outgassing at the air-sea interface during summer (see discussion in the next section). The relatively lower seasonal variability of DIC at HCB (~35 µmol/kg), as compared to DYFAMED (~58 µmol/kg) and Villefranche PointB (~68 µmol/kg) is mostly related to the relatively lower primary production at HCB (see Supplementary Figure S.7).

3.2 CO₂ air-Sea flux and annual trends

In Figure 5, the mean (2010-2022) air-sea CO₂ flux (CO_2flux), obtained from the 13-year simulation is shown. This is calculated from Equations 1, 2, based on the difference between water pCO₂ from CO₂ air concentration (pCO_2air). As no spatial variability is adopted in the latter, the CO_2flux largely depends on pCO₂ seasonal variability (Figure 6) and whether this is above or below pCO_2air. Net CO₂ sink (positive values) can be seen in most of the Mediterranean, except some areas in the Eastern Levantine and the Libyan Sea that show a weak (>-1 mmol/m2/day) release of CO₂. The simulated CO₂ air-sea flux appears relatively higher in the Western Mediterranean, particularly the Gulf of Lions and Alboran Sea, as well as in river influenced areas, such as the N. Adriatic and N. Aegean and along the pathway of the Black Sea water intrusion in the Aegean. Similarly, relatively lower pCO₂ (Figure 5) is simulated in the same areas.

Figure 5 Simulated mean (2010-2022) (A) air-sea CO₂ flux (mmol/m²/day), (B) surface water pCO₂ (µatm) and (C) sea surface temperature (°C). Positive air-sea flux values indicate net CO₂ ocean uptake. The contour line indicating zero air-sea CO₂ flux is also shown in (A).

Figure 6 Simulated annual mean (A) air-sea CO₂ flux (mmol/m²/day) (positive values indicate net CO₂ uptake), (C) water pCO₂ (µatm), (E) pH and (G) SST (°C), and their monthly climatology in (B, D, F, H) respectively, averaged over the Gulf of Lions (3-6°E, 40.8-43°N, red line), Levantine basin (22-36°E, 30-35°N, black line) and entire Mediterranean (blue line) (see Figure 1). The adopted air pCO₂ increase is also indicated in (C) (green line).

The spatial variability of pCO₂ and CO₂ air-sea flux (Figure 5, see also Supplementary Figure S.8) is mainly controlled by temperature and primary production. The latter shows an eastward decreasing gradient (see also Supplementary Figure S.9), opposite with temperature (i.e., increases eastward). A significant part of the pCO2 and CO₂ air-sea flux spatial variability is related to the dependence of pCO₂ solubility on the water temperature (i.e., decrease of solubility with increasing temperature). Thus, pCO₂ follows an eastward increasing gradient, being higher in the Eastern Mediterranan. Moreover, in open sea areas, lower temperature is indicative of deeper winter mixing (e.g., Gulf of Lions, Southern Adriatic, North Aegean), resulting in an overall higher primary production and CO₂ net biological uptake. River influenced areas, such as the Northern Adriatic and North Aegean, are also characterized by generally lower temperatures, along with high productivity, sustained by riverine nutrient inputs. The same applies for areas receiving lateral inputs from the Atlantic (Alboran Sea) and the Black Sea (N. Aegean). In these areas, however, pCO₂ is also strongly influenced by the DIC and TA of the inflowing water. Thus, the relatively low pCO₂ in the Alboran Sea is mainly related to the DIC/TA concentrations of the Atlantic inflowing water. Overall, at basin-scale, the concurrent variability of temperature and productivity results in a significant correlation between CO₂ air-sea flux and temperature (r=-0.72). A slightly lower correlation (r=0.57) is found between pCO₂ and SST.

In Figure 6, the mean (2010-2022) monthly climatology and inter-annual variability of air-sea CO₂ flux, pCO₂ and SST are shown in two representative areas: Gulf of Lions (increased mixing/low temperature, high primary production) and Levantine basin (stratification/high temperature, low primary production), along with the basin-scale average. As mentioned above, the water pCO₂ is largely controlled by its temperature depending solubility, showing a similar seasonal cycle to SST, while the CO₂flux depends on the difference between water and atmospheric pCO₂ (mean 2010-2022 pCO₂air ~404.5 μatm, see Figure 6C), following Equation 1. In the Gulf of Lions, the simulated CO_2flux is positive throughout most of the year, with a peak during February-March (~11 mmol/m²/day) and a minimum during August (~-4 mmol/m²/day), being slightly negative (net source of CO₂) during summer period (July-September). A similar pattern may be seen in the Levantine basin (maximum ~5 mmol/m²/day, minimum~-5 mmol/m²/day), but the negative air-sea flux extends over a longer period (June-November). The mean (2010-2022) basin-scale averaged air-sea annual mean CO₂ flux is +2.1 mmol/m²/day, indicating a net sink of CO₂. A much higher mean flux (+4.1 mmol/m²/day) is simulated in the Gulf of Lions and a close to zero flux (+0.35 mmol/m²/day) in the Levantine basin.

In Table 1, the inter-annual (linear) trends are shown for different variables. These were calculated based on their monthly mean time series in the three different areas (Mediterranean, Gulf of Lions, Levantine basin), after removing their monthly climatology (not shown) as in Kapsenberg et al. (2017). The resulting trends are quite similar to those calculated based on the yearly mean variability, shown in Figure 6, although the latter trends are in many cases not statistically significant, given the small amount of data in the yearly time-series. In some cases, particularly in the Levantine basin, the calculated trends are not statistically significant (p-value>0.05 in Table 1), which might be due to the relatively short time period (13 years) considered. In the Gulf of Lions, all variables show a statistically significant trend, except the air-sea CO₂ flux, which is probably related to the important trend of pCO₂ (+1.52 µatm/year), almost parallel with the pCO_2air trend (+2.1 µatm/year). A slightly lower trend is found for pCO₂@20°C (+0.55 µatm/year), with the remaining part explained by the increasing trend of SST (+0.06 °C/year). TA shows a slightly lower trend (+0.53 µmol/kg/year), as compared to DIC (+0.77 µmol/kg/year). On Mediterranean basin-scale, trends have the same sign, but are much weaker (~1/3), as compared with the Gulf of Lions, except the air-sea flux (+0.113 mmol/m²/day/year), due to the much lower trend of pCO₂ (+0.6 µatm/year) in relation to pCO_2air. One may notice that the stronger increasing trend of pCO₂ in the Gulf of Lions, which is due to relatively higher CO2 air-sea flux in this area (Figure 6, see also Figure 5), results in a faster acidification trend (i.e., pH decrease, Figure 6E), as compared to the Mediterranean basin-average. In the Levantine basin, most variables does not show a statistically significant trend, except the air-sea flux (+0.155 mmol/m²/day/year), DIC (+0.21 µmol/kg/year) and pCO₂@20°C (+0.35 µatm/year), again driven by the increasing pCO_2air.

Table 1 Calculated annual (linear) trends over 2010-2022 period for SST (°C/year), DIC (µmol/kg/year), TA (µmol/kg/year), pCO₂ (µatm/year), pCO2@20°C (µatm/year), air-sea CO₂ flux (mmol/m²/day/year) and pH (/year) in Gulf of Lions, Levantine Basin and Mediterranean (see Figure 6).

3.3 Effect of plankton productivity and freshwater fluxes

In Figure 7, the mean (2010-2014) difference between the reference and sensitivity simulations, without the effect of biology (NoBio), E-P fluxes (NoEvap) and air-sea CO2 flux (NoAir), are shown for TA, DIC, pCO₂ and CO₂ air-sea flux, averaged in the Gulf of Lions, the Levantine basin and the entire Mediterranean (see Figure 1). The mean effect of biochemical processes on TA, resulting from the uptake/release of dissolved inorganic nutrients (see section 2.1.2) is relatively small (-2 µmol/kg on basin-average), being slightly higher in the more productive Gulf of Lions (-3µmol/kg) (see also Supplementary Figure S.10). On the other hand, de-activating biochemical processes (i.e., uptake of CO₂ by photosynthesis and release by respiration) results in a significant increase of DIC (basin-average +18.5 µmol/kg, Gulf of Lions +25.7 µmol/kg, Levantine +9.1 µmol/kg) and pCO₂ (basin-average +36 µatm, Gulf of Lions +53 µatm, Levantine +18 µatm) (see also Supplementary Figures S.10, S.11). We should note that in this case the CO₂ air-sea flux becomes negative (CO₂ net source) in the entire basin, including the Gulf of Lions. De-activating the effect of E-P fluxes results in a quite similar decrease on DIC and TA (Basin-average TA -15.1 µmol/kg and DIC -15.0 µmol/kg), which is much stronger in the Levantine basin (TA -19.8 μmol/kg, DIC -20.7 μmol/kg), as compared to the Western Mediterranean (Gulf of Lions TA -10.5 μmol/kg, DIC -10.6 μmol/kg)). This similar decrease of DIC and TA results in a small but relevant decrease of pCO₂, particularly in the Levantine basin (basin-average -4.8 µatm, Gulf of Lions -3.5 µatm, Levantine -8.3 µatm), which is characterized by significant E-P losses. Finally, de-activating the CO₂ air-sea flux has a relatively smaller overall effect on DIC (basin-average -2.6 µmol/kg, Gulf of Lions -5.4 µmol/kg, Levantine +2.1 µmol/kg). We should note that de-activating CO₂ air-sea flux results in an increase of DIC in the Levantine area, where there is a net CO₂ release in the reference simulation over the 2010-2014 period. Overall, the relative changes mentioned above provide an indication of the importance of different processes, in relation to the simulated pCO₂ and the air-sea flux, as well as the potential model bias, when for example those processes are overestimated/underestimated in the model.

Figure 7 Difference (sensitivity-reference) between mean (2010-2014) surface (A) total alkalinity (µmol/kg, TA), (B) dissolved inorganic carbon (µmol/kg, DIC), (C) water pCO₂ (µatm) and (D) CO₂ air-sea flux (mmol/m²/day), in the reference simulation and the simulations without the effect of evaporation-precipitation fluxes (NoEvap, green bar), without the effect of biology (NoBio, blue bar) and without the air-sea CO₂ flux (NoAir, yellow bar), averaged in the Gulf of Lions (3-6°E, 40.8-43°N), the Levantine basin (22-36°E, 30-35°N) and the entire Mediterranean (see Figure 1). The air-sea CO₂ flux in the case of NoAir sensitivity simulation is calculated, but not included in the model.

As shown in Figure 8, the effect of CO₂ net uptake from photosynthesis/respiration on DIC and pCO₂ seasonal variability peaks during early spring period (see Figures 8E–H), which is slightly delayed as compared with the primary production peak (winter). This is because during winter period, DIC is primarily controlled by vertical mixing and thus all simulations are converging (see Figures 8A, B). In the Gulf of Lions, stronger vertical mixing and increased biological activity (see SupplementaryFigure S.12) results in a much higher seasonal variation of DIC (~42.5 µmol/kg), as compared to the Levantine basin (~26.8 µmol/kg). When biological processes are switched-off (simulation NoBio), the DIC seasonal variation is significantly decreased (Gulf of Lions -37%, Levantine -36%) in both areas (Gulf of Lions ~ 26.8 µmol/kg, Levantine ~17.4 µmol/kg). The effect of E-P fluxes on DIC has the strongest impact during late summer period, synchronously with the peak of net evaporation (see Figures 8E, F). In the Levantine basin, one may notice the stronger effect of evaporation, as compared to biological processes, given the low productivity and high evaporation fluxes in this area. The effect of CO₂ air-sea fluxes in the Gulf of Lions shows a peak in late spring, slightly delayed with regard to the pCO₂ minimum (late winter), as with the effect of biological processes. On the other hand, in the Levantine, de-activating the CO₂ air-sea fluxes has an opposite effect (DIC and pCO₂ increase, see Figures 8B, D), which peaks during summer period (see Figures 8F, H), when pCO₂ is higher (i.e., higher air-sea CO₂ export in the reference simulation). The effects of biology and air-sea fluxes on pCO₂@20°C appear quite similar with those on DIC (see Figures 8C, D, G-H). On the other hand, as E-P fluxes affect both DIC and TA, they have a much lower impact on pCO₂, as compared to DIC. This is why pCO₂ shows a slightly different seasonal variability from DIC in the reference simulation (see Figures 8A–D), being relatively lower during the late summer-autumn period, when E-P fluxes peak.

Figure 8 Monthly climatology (2010-2014) of simulated concentration of (A, B) DIC (µmol/kg) and (C, D) pCO₂@20 C° (µatm) in the reference simulation (Ref, black line) and the simulations without the effect of evaporation-precipitation fluxes (NoEvap, green line), without the effect of biology (NoBio, blue line) and without the air-sea CO₂ flux (NoAir, red line) and difference of (E, F) DIC (µmol/kg) and (G, H) pCO₂@20 C° (µatm) (Ref-NoEvap: green line, NoBio-Ref: black line, Ref-NoAir: magenta line), averaged in the Gulf of Lions (3-6°E, 40.8-43°N) (left) and Lenantine basin (22-36°E, 30-35°N) (right).

4 Discussion and concluding remarks

A basin-scale Mediterranean carbonate system model has been setup, building on the POSEIDON operational biogeochemical model (Kalaroni et al., 2020a, b). The spatial variability of simulated carbonate system variables (DIC, TA, pCO₂) from a 13-year simulation (2010-2022) was validated against CARIMED in situ data, showing reasonable agreement, with the exception of a DIC and TA overestimation along the Atlantic water pathway (Algerian basin and Sicily channel) and also a slight overestimation of TA in the Eastern Levantine basin. This might be related to an overestimation of evaporation fluxes in the model and/or a small bias of the initial TA fields, as suggested by Cossarini et al. (2021), who found a similar bias in their model results. The time-variability of the model simulated outputs (DIC, TA, pCO₂, pH) were also validated, against available time-series from Western (DYFAMED, Villefranche-PointB) and Eastern Mediterranean (HCB) sites. Data from HCB station provided the first time-series with measurements of pCO₂ in the Eastern Mediterranean (Adriatic Sea excluded), which was accompanied by pH, DIC and TA measurements (see Frangoulis et al., 2024, this issue). The simulated variability was found in good agreement with the data, particularly for pCO₂, pH and DIC. However, the model failed to reproduce the observed late summer/autumn peak of TA at PointB site, which is related with evaporation. This model deviation is related to the advection of lower alkalinity waters of Atlantic origin towards the Ligurian Sea during spring/summer period (see SupplementaryFigure S.13). A comparison of the near surface circulation, as depicted by the simulated sea surface height, against satellite altimetry data (see Supplementary Figure S.14) suggests a model deviation in the anti-cyclonic pattern west of Sardinia, in the Algerian basin This deviation results in a overestimated northward propagation of the Atlantic water. Such a northward pathway of Atlantic water through anti-cyclonic eddies and instabilities formed in the Algerian basin has been documented in the Western Mediterranean (see e.g., Béthoux, 1980; Millot and Taupier-Letage, 2005; Schroeder et al., 2008). However, the dominant pathway of Atlantic water carried by the Algerian current is eastward (e.g. Schroeder et al., 2008), bifurcating toward the Tyrrhenian basin and Eastern Mediterranean (through the straits of Sicily) and this northward pathway is probably exaggerated in the model simulation. This has also been confirmed with a comparison with surface salinity, obtained from CMEMS re-analysis (see Supplementary Figures S.15-S.16), which shows a relatively lower model simulated salinity toward the Ligurian Sea and relatively higher salinity in the Tyrrhenian basin. The overestimated advection of Atlantic water towards Ligurian Sea results in an underestimation of TA at DYFAMED and Villefranche-PointB, mostly noticed during spring/summer stratified period, as during winter this effect is largely compensated by vertical mixing. The same effect also explains the DIC and pCO₂ underestimation at DYFAMED and Villefranche-PointB during spring, given that the spring bloom net primary production is probably underestimated in the model, as indicated by a comparison of simulated with satellite Chl-a (see SupplementaryFigures S.17, S.18). The same comparison at HCB indicates an overestimation of primary production, due to overestimated winter vertical mixing in the Eastern Mediterranean (Kalaroni et al., 2020a, b), which explains the underestimation of pCO₂ at HCB during spring.

The relatively higher time-resolution of in situ data at Villefranche-PointB allowed a more detailed assessment of the simulated DIC, pCO₂ and pH seasonal variability. In particular, a time-lag (1-2 months) of the simulated winter DIC peak could be seen in some years (e.g., 2013, 2015, 2016, 2019), which was reflected also in the variability of pH and pCO₂ (Figure 3, see also Supplementary Figure S.19). This could be related to a different timing of phytoplankton bloom in the area, revealed in some cases from a comparison with satellite Chl-a (Supplementary Figure S.18). However, a similar time lag could be also seen in the sensitivity simulation with the effect of biology switched-off (see section 3.3). Moreover, the simulated SST variability was found very well correlated (r~0.99) with in situ SST, showing a similar winter minimum during February-March period (see Supplementary Figure S.19). Therefore, the simulated DIC winter time-lag in some years can be attributed to the inter-annual variability of near surface circulation in the area and the advection of relatively higher DIC water from the Tyrrhenian basin (not shown).

Along with the reference simulation, three sensitivity simulations were performed, de-activating the effect of biology, evaporation and CO₂ air-sea fluxes on DIC and TA, in order to gain insight into the processes regulating the model simulated carbonate system variability. The impact of biochemical processes on TA was relatively small, being slightly stronger in the more productive North Western Mediterranean, in agreement with Cossarini et al. (2015a) findings. On the other hand, biochemical processes had an important impact on DIC, particularly in the more productive North Western Mediterranean (annual mean difference ~26 µmol/kg), while evaporation had a stronger impact in the Levantine basin (annual mean difference ~21 µmol/kg). As evaporation fluxes affect equally DIC and TA, pCO₂ was still mostly affected by biological processes throughout the basin (mean difference Gulf of Lions +53 µatm, Levantine +18 µatm). Thus, when biological processes are switched-off the CO₂ air-sea flux becomes negative (CO₂ net source) in the entire basin, which is consistent with findings of Canu et al. (2015). An interesting finding in the present study is that when net evaporation that affects equally DIC and TA is switched-off, a small but relevant decrease of pCO₂ is simulated, particularly in the Levantine basin, which is characterized by significant E-P losses. Part of this effect is related to the slightly higher dependence of pCO₂ on DIC, as compared to TA (see Equation 3, 4 in Takahashi et al., 1993). Assuming a 1% decrease in both DIC (-22.8 µmol/kg) and TA (-27.2 µmol/kg), which is representative of the effect of E-P fluxes in the Levantine basin (see Supplementary Figure S.20), a decrease of 0.63% (-2.76 µatm) can be calculated for pCO₂, using CO2SYS and HALTAFALL codes (constant temperature=22°C and salinity =39). However, the simulated mean decrease of pCO₂ in the Levantine basin, over the 5-year (2010-2014) period, is relatively higher (~2%, -8.3 µatm, see Figure 7), as compared to this theoritical calculation. This is because the increase due to evaporation counterbalances the DIC removal from biological activity and, thus, DIC shows a slightly stronger decrease (-0.9%), as compared to TA (-0.76%), when E-P is de-activated (see Supplementary Figure S.20). This small difference between DIC and TA response results in the relatively stronger decrease of pCO2. This process might be important in areas with strong evaporation, such as the Eastern Mediterranean.

The simulated annual mean CO₂ air-sea flux was found higher in the Western Mediterranean, particularly in the Gulf of Lions and the Alboran Sea, as well as in river influenced areas, such as the N. Adriatic and along the pathway of the Black Sea Water in the Aegean. A weak release of CO₂ was found in the Eastern Levantine and Libyan Sea. The simulated air-sea flux spatial variability is consistent with existing estimates from other observational and modelling studies (DYFAMED: +3.8 mmol/m²/day, Copin-Montégut and Bégovic, 2002; N. Aegean +6.2/11.8 mmol/m²/day, Krasakopoulou et al., 2009; N. Adriatic: +2.9 mmol/m²/day, Cossarini et al., 2015a; South-Eastern Levantine: -2.3 mmol/m²/day, Sisma-Ventura et al., 2017; S. Aegean +0.5 mmol/m²/day, Frangoulis et al., 2024, this issue). It is also similar with results of Cossarini et al. (2021), except in the Alboran Sea, where air-sea CO₂ flux is relatively higher in the present simulation, due to relatively lower pCO₂, based on the adopted Atlantic water DIC and TA. The strong impact of the Atlantic water properties on the air-sea CO₂ flux in this area has been also noted by Canu et al. (2015). The simulated mean (2010-2022) CO₂ air-sea flux (+2.1 mmol/m²/day), averaged over the entire Mediterranean indicates a net CO₂ uptake. This is relatively higher (mean 2010-2019: +1.5 mmol/m²/day), as compared to the, also positive (net uptake), basin-average air-sea CO₂ flux (+0.35 mmol/m²/day), simulated by CMEMS re-analysis over 2010-2019 period (Cossarini et al., 2021), which may be partly attributed to the relatively higher Alboran Sea air-sea CO₂ flux, simulated in the present study.

The simulated basin average air-sea CO₂ flux presents a positive trend (+0.11 mmol/m²/day/year), indicating a gradual increase of net CO₂ uptake. This may be attributed to the lower increasing trend of water pCO₂ (+0.60 mmol/m²/day/year), as compared to the air CO₂ (+2.1 mmol/m²/day/year). Part of the water pCO₂ increasing trend is also related to the increasing water temperature (SST trend=+0.0125 °C/year, pCO₂@20°C trend=+0.34 µatm/year). A slightly higher increasing trend (+0.15 mmol/m²/day/year) was simulated by CMEMS re-analysis model output over 2010-2019 period (Cossarini et al., 2021). The in situ pCO₂ data at Villefrance-PointB show a much higher trend (+3.53 μatm/year, Kapsenberg et al., 2017), as compared with the model (+1.51 μatm/year). However, the CO₂ dynamics in such coastal areas are probably not fully represented, being complicated by additional processes (fresh water inputs, benthic processes etc.) (Cossarini et al., 2021). The simulated positive trends of DIC and TA in the North Western Mediterranean (DIC: +0.77 µmol/kg/year, TA: +0.53 µmol/kg/year) are consistent with the simulated trends in Cossarini et al. (2021) (Western Mediterranean DIC: +0.7 µmol/kg/year, TA: +0.6 µmol/kg/year). In the Levantine basin, a relatively smaller trend (+0.21 µmol/kg/year) was simulated for DIC. However, no significant trend was simulated for TA, in contrast to Cossarini et al. (2021), who reported stronger trends for both DIC and TA (Eastern basin DIC: +2 µmol/kg/year, TA: +1.5 µmol/kg/year), as compared to the Western basin (DIC: +0.7 µmol/kg/year, TA: +0.6 µmol/kg/year). These stronger increasing trends in the Eastern Mediterranean were found aligned with simulated increasing salinity (Escudier et al., 2021) and other observational evidence (e.g. Wimart-Rousseau et al., 2021). Such an increasing trend for Eastern Mediterranean salinity and TA was not simulated in the present study, which is consistent with E-P fluxes, also not showing any increasing trends to explain this increase. Despite that, there has been evidence for long-term trends of increasing salinity and freshwater fluxes in the Mediterranean (1950-2010, see e.g Skliris et al., 2018), but the simulated inter-annual variability was probably not long-enough to properly identify such long-term trends, as this was not the primary focus of this study.

Besides the relatively short simulation period, there are also other important model limitations and uncertainties of the simulated carbonate system. The most significant uncertainties are related with the initial DIC and TA fields from GLODAP climatology, which might be biased in areas with scarce data availability (e.g. Eastern Mediterranean), as well as the open boundary conditions (Atlantic and Black Sea water inflow) and riverine inputs of DIC and TA that are still based on fragmented information. Simulated CO₂ air-sea fluxes are also characterized by an important uncertainty, as water pCO₂ data are still scarse, while derived carbonate data often present significant discrepencies, as compared with the observed (Artioli et al., 2012). Indicative of this uncertainty is the difference between the estimate of CO₂ air-sea fluxes in the present study and other basin-scale model results (e.g. Cossarini et al., 2021). To reduce this model uncertainty, multi-model ensembles are currently becoming common practice, particularly for future climate projections (see e.g., Gooya et al., 2023). Another important limitation of the present study is that the effect of calcium carbonate formation/dissolution on TA and DIC is not considered in the model, as calcifying organisms are not explicitly described, due to the limited information on their abundance and distribution in the Mediterranean (e.g., Siokou-Frangou et al., 2010). The formation of 1 mole calcium carbonate results in the decrease of TA by 2 moles and DIC by 1 mole (Wolf-Gladrow et al., 2007). Given this differential effect on TA and DIC, calcium carbonate formation in the euphotic zone results in an increase of pCO₂ (i.e., a decrease of air-sea CO₂ flux). Missing this effect in the model, results in a potential underestimation of pCO₂ and an overestimation of DIC and particularly TA, as carbonate formation/sedimentation is considered a major sink for alkalinity in the Mediterranean (Schneider et al., 2007). Recent satellite data may provide information on the distribution of different phytoplankton functional types, including calcifying organisms (e.g., Di Cicco et al., 2017; El Hourany et al., 2019), which can be used to explicitly describe this process in future modelling efforts.

The present study is a first step in setting up a new carbonate system model in the Mediterranean, with significant skill in reproducing the observed spatio-temporal variability. Next steps that are expected to further improve the model skill include the parameterization of calcium carbonate formation/dissolution, the update of riverine and Black Sea water DIC/TA inputs and particularly the use of more realistic boundary conditions for the Atlantic water inflow (i.e. seasonally/inter-annually varying DIC/TA), which is an important driver for acidification trends in the Mediterranean (e.g., Palmiéri et al., 2015; Wimart-Rousseau et al., 2023).

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

KT: Conceptualization, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing, Formal Analysis. CF: Formal Analysis, Investigation, Methodology, Project administration, Writing – review & editing. NS: Formal Analysis, Investigation, Methodology, Software, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was partially supported by the JERICO-S3 project. The JERICO-S3 project is funded by the European Commission’s H2020 Framework Programme under grant agreement No. 871153.

Acknowledgments

We thank Dr. Louisa Giannoudi for providing bottle in situ data at HCB station. We also thank Professor Evangelia Krasakopoulou and Dr. George Petihakis for the valuable comments on the manuscript.

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.

Publisher’s note

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

Supplementary material

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

References