Implementation and assessment of a carbonate system model (Eco3M-CarbOx v1.1) in a highly dynamic Mediterranean coastal site (Bay of Marseille, France)

A carbonate chemistry balance module was implemented into a biogeochemical model of the planktonic food web. The model, named Eco3M-CarbOx, includes 22 state variables that are dispatched into 5 compartments: phytoplankton, heterotrophic bacteria, detrital particulate organic matter, labile dissolved organic, and inorganic matter. This model is applied to and evaluated in the Bay of Marseille (BoM, France), which is a coastal zone impacted by the urbanized and industrialized Aix–Marseille Metropolis, and subject to significant increases in anthropogenic emissions of CO₂.

The model was evaluated over the year 2017, for which in situ data of the carbonate system are available in the study site. The biogeochemical state variables of the model only change with time, to represent the time evolution of a sea surface water cell in response to the implemented realistic forcing conditions. The model correctly simulates the value ranges and seasonal dynamics of most of the variables of the carbonate system except for the total alkalinity. Several numerical experiments were conducted to test the response of carbonate system to (i) a seawater temperature increase, (ii) wind events, (iii) Rhône River plume intrusions, and (iv) different levels of atmospheric CO₂ contents. This set of numerical experiments shows that the Eco3M-CarbOx model provides expected responses in the alteration of the marine carbonate balance regarding each of the considered perturbation. When the seawater temperature changes quickly, the behavior of the BoM waters alters within a few days from a source of CO₂ to the atmosphere to a sink into the ocean. Moreover, the higher the wind speed is, the higher the air–sea CO₂ gas exchange fluxes are. The river intrusions with nitrate supplies lead to a decrease in the pCO₂ value, favoring the conditions of a sink for atmospheric CO₂ into the BoM. A scenario of high atmospheric concentrations of CO₂ also favors the conditions of a sink for atmospheric CO₂ into the waters of the BoM. Thus the model results suggest that external forcings have an important impact on the carbonate equilibrium in this coastal area.

1 Introduction

Current climate change mostly originates from the carbon dioxide (CO₂) increase in the atmosphere at a high annual rate (+2.63 ppm from May 2018 to May 2019, https://www.esrl.noaa.gov/gmd/ccgg/trends/global.html, last access: September 2019). This atmospheric CO₂ increase impacts the carbonate chemistry equilibrium of the oceanic water column (Allen et al., 2009; Matthews et al., 2009). Oceans are known to act as a sink for anthropogenic CO₂, i.e., 30 % of emissions, which leads to a marine acidification (Gruber et al., 2019; Orr et al., 2005; Le Quéré et al., 2018). CO₂ is a key molecule in the biogeochemical functioning of the marine ecosystem. Photo-autotrophic organisms, mainly phytoplankton and macro-algae, fix this gas through photosynthesis in the euphotic zone and, in turn, produce organic matter and dissolved oxygen. Heterotrophic organisms, mainly heterotrophic protists and metazoans, consume organic matter and dissolved oxygen by aerobic respiration and, in turn, produce CO₂. In the ocean, the main processes regulating CO₂ exchanges between the atmosphere and sea are the solubility pump and the biological pump at different timescales. Overall, the thermohaline gradients drive the solubility pump, while the metabolic processes of gross primary production and respiration set the intensity of the biological pump (Raven and Falkowski, 1999).

The coastal zones, despite their small surface area and volume compared to those of the open ocean, have a large influence upon carbon dynamics and represent 14 % to 30 % of the oceanic primary production (Gattuso et al., 1998). At the interface between open ocean and continents, these zones receive large inputs of nutrients and organic matter from rivers, groundwater discharge, and from atmospheric depositions (Cloern et al., 2014; Gattuso et al., 1998). On coasts, shorelines are subject to an increasing density of population and associated urbanization (Small and Nicholls, 2003). This rapid alteration of shorelines all over the world accelerates the emissions of greenhouses gases near the coastal ocean, and it also involves large discharges of material into the seawater by wastewater runoff and/or rivers (Cloern, 2001). These anthropogenic forcings alter the biogeochemical functioning of these zones and could lead to a growing eutrophication (Cloern, 2001). Moreover, these forcings could affect the carbonate chemistry dynamics and amplify or attenuate the acidification in coastal zones. This alteration of the marine environment may provoke further changes in the structure of the plankton community, including, in fine, consequences on the populations with high trophic levels, such as teleosts (Esbaugh et al., 2012). At the global scale, coastal zones are considered to be a significant sink for atmospheric CO₂, with an estimated flux converging to 0.2 Pg C yr⁻¹ (Roobaert et al., 2019). However, some studies highlight that the status of these areas as a net sink or source still remains uncertain due to the complexity of the interactions between biological and physical processes, and also due to the lack of in situ measurements (Borges and Abril, 2011; Chen et al., 2013; Chen and Borges, 2009). Moreover, the capacity for coastal zones to absorb atmospheric CO₂ resulting from the increasing human pressure also remains poorly known. There are few works which highlight, under future atmospheric CO₂ levels, whether coastal zones will become a net sink or a reduced source of CO₂ (Andersson and Mackenzie, 2012; Cai, 2011).

The current increase in the atmospheric CO₂ partial pressure (pCO₂) is slowly shifting the marine carbonate chemistry equilibrium towards increases in the seawater pCO₂ and bicarbonate ions (${\mathrm{HCO}}_{\mathrm{3}}^{-}$) and decreases in pH and carbonate ions (${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$) (Hoegh-Guldberg et al., 2018). These trends were already described in several coastal and open-ocean locations worldwide (Cai et al., 2011). In a coastal northwestern Mediterranean site, a 10-year time series of in situ measurements highlights a trend of pH decrease and pCO₂ increase (Kapsenberg et al., 2017). Low pH values can inhibit the ability of many marine organisms to form the calcium carbonate (CaCO₃) used in the making of skeletons and shells (Gattuso et al., 2015). In an extreme case, this shift may promote dissolution of CaCO₃ because the water will become under-saturated with respect to CaCO₃ minerals (Doney et al., 2009).

The present study is dedicated to the implementation of a carbonate system module into a preexisting biogeochemical model of the planktonic food web. This new model, named Eco3M-CarbOx (v1.1), is then evaluated in a highly dynamic coastal site, i.e., the Bay of Marseille (BoM) in the northwestern Mediterranean Sea. This evaluation is performed on the seasonal dynamics of biogeochemical and carbonate modeled variables against that of the corresponding in situ data available over the year 2017. This study is extended by a fine analysis of the variability of the marine carbonate system (stocks, fluxes) in relation to physical (e.g., wind events, river intrusions, temperature increases, changes in the atmospheric pCO₂ levels) and biogeochemical processes (gross primary production (GPP) and respiration, R) in the study site. The BoM is suitable for this kind of study because this coastal area is subject to high emissions of atmospheric CO₂ from the nearby urban area, and it also receives effluents from the Aix–Marseille metropolis. In addition, strong wind events (the mistral) regularly occur, which could lead to (i) strong latent heat losses at the surface (Herrmann et al., 2011) and upwelling along the coast with a common consequence of a cooling effect and (ii) Rhône River plume intrusion under specific wind conditions (Fraysse et al., 2013, 2014). In this regional context, many anthropogenic forcings can interact with the dynamics of the carbonate systems. Natural determinants of the composition of the marine planktonic community can also play a crucial role in these dynamics.

Figure 1 Schematic diagram of the biogeochemical model Eco3M-CarbOx. Explicit state variables of the model are indicated in continuous-line box or circles except the implicit variable for zooplankton (dotted line box). State variables written in orange are added variables compared to the preexisting biogeochemical model of Fraysse et al. (2013). Arrows represent processes between two state variables. TA: total alkalinity; DIC: dissolved inorganic carbon; CO₂: dissolved carbon dioxide; O₂: dissolved oxygen; CaCO₃: calcium carbonate.

2 Materials and methods

2.1 Numerical model description

The Eco3M-CarbOx biogeochemical model was developed to represent the dynamics of the seawater carbonate system and plankton food web in the BoM. The model was implemented using the Eco3M (Ecological Mechanistic and Modular Modelling) platform (Baklouti et al., 2006). The model structure used is based on an existing model of the plankton ecosystem (Fraysse et al., 2013), including a description of carbon (C), nitrogen (N), and phosphorus (P) marine biogeochemical cycles. The Eco3M-CarbOx model includes 22 prognostic state variables that are split into several compartments: phytoplankton, heterotrophic bacteria, detrital particulate organic matter, labile dissolved organic matter, nutrients (ammonia, nitrate, and phosphate), dissolved oxygen, and carbonate system variables (Fig. 1). In this study, the state variables of the Eco3M-CarbOx model only change over time (i.e., usually termed “model 0D”); they are representative of the time evolution of a sea surface water cell, but this biogeochemical model is not coupled with a hydrodynamic model.

The model presented in this study includes a set of new developments and improvements in the realism of the plankton web structure and process formulations compared to the model of Fraysse et al. (2013). In order to improve the representation of chlorophyll concentration in the Bay of Marseille the phytoplankton is divided into two groups: one with some ecological and physiological traits of the Synechococcus cyanobacteria, which is one of the major constitutive members of pico-autotrophs in the Mediterranean Sea (Mella-Flores et al., 2011), and another with traits of large diatoms, which are generally observed during spring blooms at mid-latitudes (Margalef, 1978). For both of the phytoplankton, there is a diagnostic chlorophyll a variable related to the phytoplankton C biomass, the phytoplankton N-to-C ratio, and the limiting internal ratio $f_Q^N$ (Faure et al., 2010; Smith and Tett, 2000; Table B2, Appendix B). The functional response of primary production was modified using another formulation of the temperature limitation function which takes into account the optimal temperature of growth for each phytoplankton group. The exudation of phytoplankton was modified taking into account the intracellular phytoplankton ratio. For the uptake of matter by bacteria and the remineralization processes the dependence on the intracellular bacteria ratio was modified. A temperature dependence of all biogeochemical processes was added to take into account the effects of rapid and strong variations of seawater temperature on plankton during episodes of upwelling, for instance, that are usually observed in the BoM. Also certain parameters in some formulations were modified owing to the alterations of some formulations (Table B4, Appendix B).

Additionally, a carbonate system module was developed and three state variables were added: dissolved inorganic carbon (DIC), total alkalinity (TA) and the calcium carbonate (CaCO₃) implicitly representing calcifying organisms. The knowledge of DIC and TA allows the calculation of the pCO₂ and pH (total pH scale) diagnostic variables, necessary for resolving all the equations of the carbonate system. These equations use apparent equilibrium constants, which depend on temperature, pressure, and salinity (Dickson, 1990a, b; Dickson and Riley, 1979; Lueker et al., 2000; Millero, 1995; Morris and Riley, 1966; Mucci, 1983; Riley, 1965; Riley and Tongudai, 1967; Uppström, 1974; Weiss, 1974). The details of the resolution of carbonate system module are given in Appendix A. For this module three processes were also added: the precipitation and dissolution of calcium carbonate and the gas exchange of pCO₂ with the atmosphere. Based on the review of Middelburg (2019), it is considered that (i) TA decreases by 2 moles for each mole of CaCO₃ precipitated, by 1 mole for each mole of ammonium nitrified, by 1 mole for each mole of ammonium assimilated by phytoplankton, and TA increases by 2 moles for each mole of CaCO₃ dissolved, and by 1 mole for each mole of organic matter mineralized by bacteria in ammonium (Table B2, Appendix B); (ii) DIC is consumed during the photosynthesis and calcification processes and is produced by respiration (of phytoplankton, zooplankton, and bacteria) and the CaCO₃ dissolution processes. Moreover, the dynamics of DIC are altered by the CO₂ exchanges with the atmosphere (Table B2, Appendix B). The air–sea CO₂ fluxes are calculated from the pCO₂ gradient across the air–sea interface and the gas transfer velocity (Table B3, Appendix B) estimated from the wind speed and using the parametrization of Wanninkhof (1992).

In the Eco3M-CarbOx model, zooplankton is considered as an implicit variable. However, a closure term based on the assumption that all of the matter grazed by the zooplankton and higher trophic levels return as either organic or inorganic matter by excretion, egestion, and mortality processes is taken into account (Fraysse et al., 2013). The model considers a “non-Redfieldian” stoichiometry for phytoplankton and bacteria. All the biogeochemical model formulations, equations, and associated parameter values are detailed in Appendix B.

Figure 2 Map of study area: the Region Sud (a), Aix–Marseille Metropolis (area in the red rectangle) (b), the Bay of Marseille (c). CAV = Cinq Avenues Station (urban site); OHP: Observatoire de Haute Provence station (non-urban site); Carry, Solemio, Planier: marine study sites in the Bay of Marseille.

2.2 Study area

The BoM is located in the eastern part of the Gulf of Lion, in the northwestern Mediterranean Sea (Fig. 2). The city of Marseille, located on the coast of the BoM, is the second largest city of France, with a population of ca. 1 million. The Rhône River, which flows into the Gulf of Lion, is the greatest source of freshwater and nutrients for the Mediterranean Sea, with a river mean flow of 1800 m³ s⁻¹ (Pont et al., 2002). Several studies highlight the eastward intrusion events from the Rhône River plume in the BoM under east and southeasterly wind conditions, which favor biological productivity (Fraysse et al., 2014; Gatti et al., 2006; Para et al., 2010). The biogeochemistry of the BoM is complex and highly driven by hydrodynamics. For example, north–northwesterly winds induce upwelling events which bring cold and nutrient-rich waters upward (Fraysse et al., 2013). Moreover, the oligotrophic Northern Current occasionally intrudes into the BoM (Petrenko, 2003; Ross et al., 2016).

Despite the presence of several marine protected areas around the BoM (the Regional Park of Camargue, the Marine protected area Côte Bleue, and the National Park of Calanques), it is strongly impacted by diverse anthropogenic forcing, because industrialized and urbanized areas are located all along the coast. From the land, the BoM receives nutrients and organic matter from the urban area of the Aix–Marseille metropolis (Millet et al., 2018), the industrialized area of Fos-sur-Mer (one of the biggest oil-based industry areas in Europe), and the Berre Lagoon, which is eutrophized (Gouze et al., 2008; Fig. 2c). From the atmosphere, the BoM is subject to fine-particle deposition and greenhouse gas emissions (including CO₂) from the nearby urban area, and it also receives effluents from the Aix–Marseille metropolis.

2.3 Dataset

The modeled variables of the carbonate system (DIC, TA, pH, and pCO₂) and chlorophyll a are hereafter compared to observations collected at the SOLEMIO station (Figs. 2c and 3), which is a component of the French national monitoring network (Service d'Observation en Milieu Littoral – SOMLIT, http://somlit.epoc.u-bordeaux1.fr/fr/, last access: January 2020). Major biogeochemical parameters have been recorded since 1994. Carbonate chemistry variables (pH, pCO₂, DIC, and TA), sampled every 2 weeks, have been available since 2016.

Table 1 Forcing of the different scenarios (S) simulated with the model. See Sect. 2.4 for details of the scenarios.

2.4 Design of numerical experiments

In the present work, the Eco3M-CarbOx model was run for the whole year of 2017. This year was chosen because in situ data of carbonate systems (DIC, TA, pH, and pCO₂) are available for the whole year at the SOLEMIO station (Fig. 2c). The biogeochemical variables were initialized using in situ data from winter conditions (Table B1, Appendix B). The model was forced by time series of sea surface temperature and salinity, wind (at 10 m), light, and atmospheric CO₂ concentrations. The sea temperature time series is from in situ hourly data recorded at the Planier station (Fig. 2c). For salinity, hourly in situ data from the SOLEMIO station and from the CARRY buoy were used (Fig. 2c). Wind and light hourly time series were extracted from the WRF meteorological model at the SOLEMIO station (Yohia, 2017). Finally, we used hourly atmospheric CO₂ values from in situ measurements recorded at the Cinq Avenues station (CAV station; Fig. 2b) by the AtmoSud Regional Atmospheric Survey Network, France (https://www.atmosud.org, last access: August 2019). This simulation is the reference simulation (noted S0). As highlighted previously, Rhône River plume intrusions (due to wind-specific conditions) have an impact on the dynamics of primary production (Fraysse et al., 2014; Ross et al., 2016) and then on the seawater carbonate system. Moreover, the seawater temperature and atmospheric CO₂ variations control the seawater CO₂ dynamics via the solubility equilibrium and gas exchange with the atmosphere (Middelburg, 2019). In order to quantify the impact of different forcing, several simulations (hereafter noted S), which are summarized in Table 1, were conducted:

• Impact of seawater temperature increase, S1. The forcing time series of in situ temperatures was shifted by +1.5 ^∘C (Cocco et al., 2013).

• Impact of wind events. A first simulation S2 was run with a constant wind intensity of 7 m s⁻¹ (2017 annual average wind speed) throughout the year and a second one (S3) with two three-day periods of strong wind speed (20 m s⁻¹) representative of short bursts of the mistral (data not shown) starting on 15 May and 15 August, and a constant value of 7 m s⁻¹ the rest of the year.

• Impact of nutrient supply (nitrate) during a Rhône River plume intrusion (S4). A threshold of 37 has been chosen to identify the presence of low-salinity waters from the Rhône River plume in the forcing file of salinity. Here, concentrations in nitrate supplied by the river depend on the salinity. A relationship between these two variables was then established for the SOLEMIO point from the MARS3D-RHOMA coupled physical and biogeochemical model (Fraysse et al., 2013; Pairaud et al., 2011). This relationship has already been used successfully to reproduce realistic observed conditions in the studies of Fraysse et al. (2014) and Ross et al. (2016): NO_3intrusion (mmol m⁻³) = $-\mathrm{1.70}\times S+\mathrm{65}$.

• Non-urban atmospheric CO₂ concentrations (S5). This simulation takes into account the forcing of atmospheric CO₂ values measured at the Observatoire de Haute Provence station (OHP; Fig. 2b), located outside of the Aix–Marseille metropolis, by the ICOS National Network, France (http://www.obs-hp.fr/ICOS/Plaquette-ICOS-201407_lite.pdf, last access: August 2019).

In this work, we calculated the daily mean values of state variables, statistical parameters, and mean fluxes of modeled processes throughout the year and over two main hydrological periods: the stratified and mixed water column periods. The stratified water column (SWC) is defined with a temperature difference between the surface and bottom of more than 0.5 ^∘C (Monterey and Levitus, 1997). For the simulated year (2017), the SWC period lasts from 10 May to 20 October. The mixed water column (MWC) period corresponds to the rest of the year.

3 Results

3.1 Model skills

Following the recommendations of Rykiel (1996), three criteria were considered to evaluate the performance of our model:

• Does the model reproduce the timing of the observed variations of carbonate system at the seasonal timescale?

• Does the model reproduce the observed pCO₂ and pH ranges at the seasonal timescale?

• Are the observations correctly represented according to the Willmott skill score (WSS)? This index is an objective measurement of the degree of agreement between the modeled results and the observed data. A correct representation of observations by the model is achieved when this index is higher than 0.70 (Willmott, 1982).

Figure 3 Comparison of model results (red) and in situ data (blue) at the surface of the SOLEMIO station. (a) Chlorophyll a concentrations (mg m⁻³), (b) pCO₂ (µatm), (c) pH, (d) DIC (µmol kg⁻¹), (e) TA (µmol kg⁻¹). The value of each state variable represents the mean around ± 5 d of the seawater sampling date. The shaded area and dotted black line delimit the SWC and MWC periods.

Table 2 Statistical evaluation of observations vs. model for 2017: observed and simulated minimum and maximum values. WSS = Willmott skill score; N = number of measurements. Units of bias are those of modeled variables: chlorophyll a (Chl a, mg m⁻³), seawater partial pressure of CO₂ (seawater pCO₂, µatm), pH, dissolved inorganic carbon (DIC, µmol kg⁻¹), and total alkalinity (TA, µmol kg⁻¹).

^* Significant value of WSS (> 0.70).

Over most of the studied period, the model simulates lower chlorophyll a concentrations than the in situ observations, especially during the MWC period (Fig. 3a). Two maxima of chlorophyll a concentrations are observed in situ: the first one at ca. 1.71 mg m⁻³ in March and the second one at ca. 0.68 mg m⁻³ in May. They are both linked to Rhône River plume intrusions. Several in situ maxima between 0.50 and 0.70 mg m⁻³ are observed between March and April (at the end of the MWC period), and they signaled the spring bloom event (Table 2 and Fig. 3a). The biogeochemical model quantitatively reproduces the spring bloom observed at the end of the MWC period (Fig. 3a) with a maximum value of ca. 0.69 mg m⁻³. The model does not catch the two aforementioned maxima of chlorophyll, and it contains a low WSS and a strong bias (0.37 and +0.22 mg m⁻³, respectively – Table 2).

On the whole, the seasonal variations of the seawater pCO₂ are correctly simulated by the biogeochemical model (Fig. 3b), even if the values are rather overestimated during the MWC period. From January to February, the model reproduces the slight decrease in the observed pCO₂, and from February to March it reproduces the increase in pCO₂ even if the amplitude of this increase remains smaller. In mid-April, during the simulated spring bloom period, the observed drop in pCO₂ and increase in pH are also spotted in the model (Fig. 3b and c). The model especially succeeds in reproducing the observed increase in relation to high temperatures during the SWC period. The reduction of the CO₂ solubility due to thermal effects mostly explains the increase in pCO₂ during the SWC period. The strong SD of modeled values during the SWC period can be explained by the rapid changes in temperature, probably due to upwelling usually occurring at this time of the year (Millot, 1990). The range of modeled pCO₂ values (345–503 µatm) encompasses the range of observed values (358–471 µatm; Table 2). The statistical analysis provides a mean bias of +23 µatm, and a WSS of 0.69 (Table 2).

The seasonal dynamic of pH is mostly reproduced by the model, and in particular, the decrease during the SWC period (Fig. 3c). However, the modeled pH is generally underestimated throughout the year, except during the SWC period, with a mean bias of −0.015 (Table 2). The seasonal range is captured by the model with a minimum value during the SWC period (7.994 vs. 8.014 for observations; Table 2) and a maximum one during the MWC period (8.137 vs. 8.114 for observations; Table 2). The statistical analysis highlights an index of agreement between the in situ data and the model outputs higher than 0.70 (Table 2).

The seasonal variations of DIC show the highest values during the MWC period and a decrease (increase) during the beginning (the end) of the SWC period (Fig. 3d). The lowest values are observed during September. The Eco3M-CarbOx model closely matches the seasonal dynamic by reproducing the range of extreme observed values (Table 2). The mean bias is also small (−8.48 µmol kg⁻¹; Table 2). More than 70 % (0.73; Table 2) of modeled DIC concentrations are in statistical agreement with the corresponding observations.

The seasonal cycle of measured TA does not show a clear pattern (Fig. 3e). Large variations of values ranging between 2561 and 2624 µmol kg⁻¹ (Table 2) are observed, whatever the hydrological season being considered. The biogeochemical model provides almost constant values of around 2570 µmol kg⁻¹ all throughout the year, which is lower than in situ data. With a low WSS index of agreement and a large mean bias (Table 2), the model is not able to confidently reproduce the observed variations of TA (Fig. 3e and Table 2).

Figure 4 (a–f) 2017. (g–l) Temporal focus between 1 April and 31 July 2017. In situ daily average of (a, g) temperature (^∘C, black line) and salinity (blue line) at the SOLEMIO station (at the surface). Modeled daily average (b, h) chlorophyll a concentrations (mg m⁻³, black line) (c, i) DIC (µmol kg⁻¹, black line) and pH (blue line), (d, j) seawater pCO₂ (µatm, black line) and atmosphere pCO₂ from OHP (µatm, red line), (e, k) air–sea CO₂ fluxes ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$), (f, l) gross primary production ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, green line), total respiration ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, red line) and net ecosystem production ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, black line). The shaded areas and dotted black lines delimit the SWC and MWC periods.

Table 3 Mean flux values ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) and the contribution of each process to the DIC variations for the reference simulation over the year and SWC/MWC periods.

GPP: gross primary production, R_A: autotrophic respiration, R_H: heterotrophic respiration, NEP: net ecosystem production.

3.2 Carbon fluxes and budgets

For the year 2017, the values of temperature vary between 13.3 and 25.9 ^∘C (Fig. 4a). The DIC variations closely match those of temperature (correlation coefficient −0.75). For example, the spring increase in temperature leads to a decrease in DIC concentrations (Fig. 4a and c), and the minimum values are reached at the end of SWC period. Over the simulated period, the air–sea CO₂ fluxes (F_aera) vary between −14 and 17 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, with a weakly positive annual budget of +6 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ (or +0.017 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, Table 3). Then, the BoM waters would act as a net source of CO₂ to the atmosphere on an annual basis. However, on a seasonal basis, the BoM waters would change from a net sink during the MWC period (F_aera < 0; Table 3) to a net source during the SWC one (F_aera > 0; Table 3).

On an annual basis, the gross primary production (GPP) and total respiration (R) are balanced, leading to a null average net ecosystem production (NEP; NEP = GPP−R) (Fig. 4f and Table 3). The intensity of autotrophic respiration (R_a) is lower than that of primary production (annual mean of 0.065 vs. −0.413 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, respectively – Table 3), while the zooplankton and bacterial respiration account for an average of 0.348 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ (Table 3). On a seasonal basis, the model highlights an ecosystem dominated by autotrophy during the MWC period (NEP > 0; Table 3) and heterotrophy during the SWC period with higher flux values (NEP < 0; Table 3). The biogeochemical fluxes show the strongest variations along the SWC period, following those of temperature (Fig. 4f). The maximum GPP occurs in April and is correlated with the maximum chlorophyll concentration. At this time, the ecosystem is autotrophic (NEP > 0; Fig. 4b and f) and is a net sink for atmospheric CO₂, which explains the DIC and seawater pCO₂ decreases during the bloom period (Fig. 4c–e)

When looking in detail at the 2017 temperature and salinity time series (Fig. 4a), several crucial events can be seen occurring, including freshwater intrusions (e.g., 15 March and 6 May) into the BoM and large variations of temperature in relation to upwelling events or latent heat losses due to wind bursts. The largest freshwater intrusion from the Rhône River plume occurs in mid-March, with a minimum observed salinity of ca. 32.5 at the SOLEMIO station (Fig. 4a). During this event, the seawater pCO₂ decreases and pH increases concomitantly (Fig. 4c and d). Then, seawater appears to be temporarily under-saturated in CO₂, and the BoM waters thus act as a sink for atmospheric CO₂ at the time of intrusion (Fig. 4e).

During the SWC period, upwelling events quickly cool the surface seawater. In two days, from 25 to 27 July, the water temperature drops from 24.7 to 16.9 ^∘C (Fig. 4g). The decrease in temperature corresponds to the increase in DIC concentrations (Fig. 4i). Concomitantly, the value of seawater pCO₂ decreases from 497 to 352 µatm, and pH increases from 7.99 to 8.12 (Fig. 4i and j). This event quickly changes the BoM waters from a source to a sink for atmospheric CO₂ (from +17 to −14 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$; Fig. 4k), and also from a net heterotrophic to a net autotrophic ecosystem (Fig. 4l).

Figure 5 Modeled daily average chlorophyll a concentrations (mg m⁻³) (a), light limitation (b), temperature limitation, and a close-up from 15 April to 5 May of temperature (c) and nutrient limitation (d) for picophytoplankton and the same set for diatoms (f–i). Modeled daily average NEP ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, e) and air–sea CO₂ fluxes ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, j). Reference simulation (S0, blue line) and temperature-shifted simulation by 1.5 ^∘C (S2, red line). The shaded area and dotted black lines delimit the SWC and MWC periods.

3.3 Impact of external forcing on the dynamics of carbonate system

3.3.1 Temperature increase

Here we compare the reference simulation S0 with the S1 simulation (seawater temperature elevation of 1.5 ^∘C – Fig. 5). During the year, there are few changes on the carbonate system variables such as the pCO₂ and pH (data not shown). The main alterations occur during the blooms of phytoplankton. The simulated bloom of phytoplankton occurs later, at beginning of May, for both diatoms and picophytoplankton, with maximum values of chlorophyll at 1.4 and 0.4 mg m⁻³, respectively (Fig. 5a and f).

As both the limitations due to light and nutrients remain about the same during the simulations S0 and S1, this counterintuitive occurrence of a bloom relative to changes in temperature is mainly explained by the temperature limiting function, which depends on the optimal temperature of growth (T_opt). For the picophytoplankton, from January to April, the increase of 1.5 ^∘C drastically reduces the limitation by temperature (Fig. 5c), because the temperature is closer to the optimal temperature (T_opt = 16 ^∘C, Table A4) during S1 than S0. In the S0 simulation, the temperature reaches T_opt ca. 15 April and it induces the bloom, while at the same time in S1 the temperature moves slightly away from T_opt and it does not enable the triggering of a bloom. At the time of the bloom in S1, the opposite configuration occurs. In S0, the ambient temperature is again far from T_opt, explaining the absence of a bloom, while in the S1 the ambient temperature is closer to T_opt, enabling the occurrence of a bloom. The picophytoplankton bloom then occurs later in the warm simulation S1 than in the reference simulation S0 (Fig. 5a). The duration and termination of a bloom is controlled by both the nutrient availability and the temperature (Fig. 5c and d). Inversely, from January to April, the diatoms' growth limitation by temperature is strengthened in the warm simulation S1 (Fig. 5h), because the resulting ambient temperature is further from the optimum temperature (T_opt = 13 ^∘C, Table A4) than that in the reference simulation S0. This induces a slower growth of diatoms and a delay of the maximum concentration (Fig. 5f). Afterwards the photosynthesis is mainly limited by temperature (Fig. 5h).

The ecosystem is dominated by autotrophs at the time of blooms whatever the simulation considered (NEP > 0; Fig. 5e), and the quantity of DIC (not shown) fixed through autotrophic processes is larger than that released by heterotrophic processes. During the short period of a bloom, the seawater pCO₂ decreases, leading to some negative air–sea fluxes (i.e., an oceanic sink for atmospheric CO₂). In the warm simulation, the later occurrence of a bloom enables the period of the spring sink to extend by ca. 3 weeks over May relative to the reference simulation (Fig. 5j).

Figure 6 Temporal evolution for May (a–c) and August (d–f) 2017 of the wind speed (m s⁻¹, a, d); air–sea CO₂ fluxes ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, b, e); seawater pCO₂ (µatm, c, f). Constant wind scenario (S2, blue line) and wind event scenario (S3, red line). On (c) and (f), the dashed line represents the atmosphere pCO₂ (µatm) at the CAV station.

3.3.2 Wind speed

The Bay of Marseille is periodically under the influence of strong wind events (Millot, 1990). Here we compare two simulations: one with a constant wind value (S2) and the other one with two wind events that occur in May and August (S3) (Fig. 6a and d). The result of this numerical experiment shows that the stronger the wind speed is, the higher the air–sea fluxes are, mainly owing to the increase in gas transfer velocity. Depending on the gradient of CO₂ between seawater and the atmosphere, strong wind speeds will favor either the emission or uptake of CO₂ (Fig. 6b and e). In May, with the air–sea CO₂ flux being positive, the outgassing of CO₂ to the atmosphere is enhanced, leading to a decrease in seawater pCO₂ (Fig. 6c). On the contrary, in August the oceanic sink of atmospheric CO₂ is amplified, which leads to an increase in the seawater pCO₂ value (Fig. 6f).

Figure 7 (a–f) 2017. (g–l) Temporal focus between 1 May and 1 July 2017. (a, g) In situ daily average of salinity. Modeled daily average (b, h) nitrate concentrations (mmol m⁻³); (c, i) chlorophyll a concentrations (mg m⁻³); (d, j) DIC (µmol kg⁻¹); (e, k) seawater pCO₂ (µatm); and (f, l) pH. Reference simulation (S0, blue line) and nitrate supply simulation (S4, red line). On (e) and (k), the dashed line represents the atmosphere pCO₂ (µatm) at the CAV station. The shaded area and dotted black lines delimit the SWC and MWC periods.

3.3.3 Supply in nitrate by river inputs

According to the model results (Fig. 7), the occasional inputs of nitrate (S4) that are linked to Rhône River plume intrusions favor primary production, and they led to increased chlorophyll concentrations (Fig. 7b and c) five times during the SWC period. These blooms, as seen previously, lead to a decrease (increase) in the seawater pCO₂ (pH) (Fig. 7e and f). It can be noted that with the strongest river supply at mid-March (Fig. 7a and b) the occurrence of the spring bloom is earlier (Fig. 7c) than that occurring in the reference simulation (S0). The time lag between river nutrient supply and bloom is due to the temperature limitation (Fig. 4c). During blooms occurring within the SWC period following intrusions, the DIC concentrations are generally lower than those of the reference simulation, as in the case of the bloom of mid-May (decrease by ca. 15 µmol kg⁻¹; Fig. 7j), due to the autotrophic processes dominating the heterotrophic ones. In turn, the seawater pCO₂ drops by ca. 30 µatm (Fig. 7k) and pH increases by ca. 0.030 (Fig. 7l). Nitrate inputs, favoring primary production, reduce the source of CO₂ to the atmosphere and intensify the sink of atmospheric CO₂ into the waters of BoM (Fig. 7e and k).

Figure 8 (a) Temporal evolution for the year 2017 of the observed pCO₂ (µatm) in the atmosphere at the CAV station, called the “urban scenario” (S0, blue line), and at the OHP station, called the “non-urban scenario” (S6, red line). (b) Temporal evolution for the year 2017 of the modeled seawater pCO₂ (µatm) with forcings from the urban (S0, blue line) and non-urban (S6, red line) scenarios. The shaded area and dotted black lines delimit the SWC and MWC periods.

3.3.4 Urban air CO₂ concentrations

The Aix–Marseille metropolis is strongly subject to urban emissions to the atmosphere (Xueref-Remy et al., 2018a). The seasonal variability of atmospheric CO₂ concentrations at the urban site (CAV station; Fig. 2) is much higher than that observed in a non-urban area (OHP station; Fig. 2), especially during the MWC period (Fig. 8a): CO₂ concentrations vary between 379 and 547 µatm at the CAV station and between 381 and 429 µatm at the OHP station. Moreover, in winter the atmospheric pCO₂ is higher in the urban area than the non-urban area, whereas in summer those of both areas are quite close. These differences in the seasonal pattern and between areas are usually explained by (i) the thinner atmospheric boundary layer, (ii) the decreased fixation of CO₂ by terrestrial vegetation, and (iii) the greater influence of anthropogenic activities by emissions from heating (Xueref-Remy et al., 2018b). Forcing the model by atmospheric pCO₂ values from urban or non-urban sites can lead to significant differences in the values of the seawater pCO₂, especially during the MWC period. The air–sea gradient of pCO₂ is higher when using a forcing derived from the CO₂ concentrations originating from an urban area than from a non-urban area, which strengthens the sink of atmospheric CO₂ into the waters of BoM. The seawater pCO₂ is then lower with non-urban area pressure (S5) than with urban area pressure (S0), because of lower CO₂ solubility in the BoM (Fig. 8b).

4 Discussion

4.1 Model performance

The evaluation of model skill vs. in situ data highlights that the modeled pH, pCO₂, and DIC are in acceptable agreement with observations (Fig. 3). The seasonal variations observed for the different variables are captured by the model, including for example the seasonal decrease in DIC and pH during the SWC period, in relation to the increase in pCO₂, and the inverse scenario during the MWC period. The chlorophyll content variability is not well reproduced, especially during spring (Fig. 3a), even taking into account the nitrate supply from the Rhône River plume intrusion (Fig. 7c). This is due to the multiple origins of chlorophyll, organic matter, and nutrients in the BoM that are not accounted for in the Eco3M-CarbOx model: autochthonous marine production, and allochthonous origins from the Rhône and Huveaune River plumes (Fraysse et al., 2013). The observed variations and levels of TA are not correctly simulated by the model (Fig. 3f). The study of Soetaert et al. (2007) highlights that the main variations of TA in the marine coastal zones are linked to freshwater supplies and marine sediments. The present study does not take into account the inputs of TA from the Rhône River and the water–sediment interface, and it may explain why the TA variable is not correctly predicted by our model.

4.2 Contribution of physical and biogeochemical processes to the variability of carbonate system

The contribution of each biogeochemical process to the DIC variability can be assessed using the presented model: the aeration process contributes to 78 % of the DIC variations, and all biogeochemical processes together contribute to 22 % (Table 3). As mentioned by Wimart-Rousseau et al. (2020), the model suggests that the seawater pCO₂ variations and associated fluxes would be mostly driven by the seawater temperature dynamics. Moreover, the seasonal variations of the air–sea CO₂ flux are in agreement with some previous field studies (De Carlo et al., 2013; Wimart-Rousseau et al., 2020), which measured a weak oceanic sink for atmospheric CO₂ during winter and a weak source to the atmosphere during summer.

The model results reveal that temperature would play a crucial role in controlling two counterbalanced processes: (1) the carbonate system equilibrium and (2) the phytoplankton growth. The increase in temperature during SWC leads to a higher pCO₂ in seawater due to the decrease in the CO₂ solubility (Middelburg, 2019) and, at the same time, the fixation of DIC by phytoplankton is favored, leading to a decrease in the pCO₂ level. The imbalance between the latter two processes leads to a change in the ecosystem status (autotrophic or heterotrophic) and the corresponding behavior as a sink or source to the atmosphere. In case of a 1.5 ^∘C rise over the whole year, the temperature variation has a very small impact on the carbonate system dynamics. However, it favors the autotrophic processes and strengthens the oceanic sink of atmospheric CO₂ during the bloom of phytoplankton (Fig. 5e and j).

4.3 Contribution of the external forcing to the variability of carbonate system

In line with several previous works on the northwestern Mediterranean Sea (De Carlo et al., 2013; Copin-Montégut et al., 2004; Wimart-Rousseau et al., 2020), the model also suggests that the status of the Bay of Marseille regarding sink or source for CO₂ could change at high temporal frequency (i.e., hours to days). Bursts of north–northwestern winds lead to sudden and sharp decreases in seawater temperature (< 2 d; Fig. 4g) either directly by latent heat loss through evaporation at the surface (Herrmann et al., 2011) or indirectly by creating upwelling (Millot, 1990), with the consequences of a decrease in the seawater pCO₂ values (Fig. 4j) and, in fine, an alteration of the CO₂ air–sea fluxes. Model results suggest that the fast variations of temperature could lead to rapid changes of the sink vs. source status in this coastal zone (Fig. 4k). Moreover, Fraysse et al. (2013) highlight that upwelling in the BoM favors ephemeral blooms of phytoplankton by nutrient supplies up to the euphotic layer and would, in turn, contribute to the seawater pCO₂ decrease. North and northwestern winds through latent heat losses and/or upwelling events could then enhance the sink for atmospheric CO₂ due to the temperature drop and nutrients inputs. However, these results remain preliminary because in our experimental design only the cooling effect of upwelling on the carbonate balance is taken into account. But concomitantly, upwelling usually bring nutrients and DIC at the surface and these supplies could also perturb the balance of the carbonate system. A further coupling of the Eco3M-CarbOx model with a tridimensional hydrodynamic model would certainly enable the multiple effects of upwelling on the dynamics of the carbonate system in this area to be embraced and the results presented in this study to be refined.

High wind speeds (> 7 m s⁻¹) amplified the gaseous exchange of CO₂ considerably (De Carlo et al., 2013; Copin-Montégut et al., 2004; Wimart-Rousseau et al., 2020). The model highlights that a strong wind event of 3 d has a significant impact on the seawaterpCO₂ values during a longer period of ca. 15 d (Fig. 6). A combination of high atmospheric pCO₂ values and high wind speeds would then favor the sink for CO₂ into the waters of the BoM. The aeration process also depends on the choice of the formulation of the gas transfer velocity (k₆₀₀). In this study, the formulation of Wanninkhof (1992) is used and depends of the wind speed at 10 m above the water surface. However, the current velocity could favor the gas exchange, and suspended matter concentration could limit the gas exchange (Abril et al., 2009; Upstill-Goddard, 2006; Zappa et al., 2003). Due to the important heterogeneity of physical and biogeochemical forcings in coastal zones, other factors that control the air–sea gas exchange should certainly be taken into account.

The simulation with intrusions of the Rhône River plume shows that inputs of nitrate cause a drop of seawater pCO₂ owing to the nutrient supply favoring the phytoplankton development (Fig. 7). In this scenario, the oceanic sink of atmospheric CO₂ is enhanced. But rivers also supply TA (e.g., Gemayel et al., 2015; Schneider et al., 2007) and DIC (e.g., Sempéré et al., 2000) that shift the carbonate system equilibrium toward a pCO₂ decrease and a DIC increase (Middelburg, 2019). Taking into account these further supplies may sensibly modify the modeled carbonate balance in the BoM. A next step to the present work will be to design more realistic numerical experiments to refine the results obtained in this preliminary study. The intrusions of the Rhône River plume also induce a salinity decrease in the BoM waters, which leads to a drop in the pCO₂ levels in the model. This drop of pCO₂ is due to the decrease in the CO₂ solubility when salinity decreases (Middelburg, 2019).

In the scenario of forcing the model by using urban atmospheric pCO₂ time series, the air–sea gradient increases, and then it enhances the status of the BoM as a sink for atmospheric CO₂. As suggested by the in situ study of Wimart-Rousseau et al. (2020), the Eco3M-Carbox model highlights the crucial role of the coastal ocean in an urbanized area, and with an increase in atmospheric CO₂, the CO₂ uptake by the coastal ocean may increase. This result is in line with the studies of Andersson and Mackenzie (2004) and Cai (2011), which predict an increase in the intensity of the CO₂ sink and a potential threat to coastal marine biodiversity in coastal areas owing to high atmospheric CO₂ levels.

5 Conclusions

A marine carbonate chemistry module was implemented in the Eco3M-CarbOx biogeochemical model and evaluated against in situ data available in the Bay of Marseille (northwestern Mediterranean Sea) over the year 2017. The model correctly simulates the value ranges and seasonal dynamics of most of the variables of the carbonate system except for the total alkalinity. Several numerical experiments were also conducted to test the sensitivity of the carbon balance to physical processes (temperature and salinity), biogeochemical processes (GPP and respiration processes), and external forcing (wind, river intrusion, and atmospheric CO₂). This set of numerical experiments shows that the Eco3M-CarbOx model provides expected responses in the alteration of the marine carbonate balance regarding each of the considered perturbation.

On the whole, the model results suggest that the carbonate system is mainly driven by the seawater temperature dynamics. At a seasonal scale, the BoM marine waters appear to be a net sink of atmospheric CO₂ and a dominantly autotrophic ecosystem during the MWC period, and a net source of CO₂ to the atmosphere during the SWC period, which is mainly characterized by a dominance of heterotrophic processes. However, the model results highlight that sharp seawater cooling observed within the SWC period, probably owing to upwelling events, cause the CO₂ status of the BoM marine waters to change from a source to the atmosphere to a sink into the ocean within a few days. External forcing as the temperature increases leads to a delay in the bloom of phytoplankton. Strong wind events enhance the gas exchange of CO₂ with the atmosphere. A Rhône River plume intrusion with input of nitrate favors pCO₂ decreases, and the sink of atmospheric CO₂ into the BoM waters is enhanced. The higher atmospheric pCO₂ values from the urban area intensify the oceanic sink of atmospheric CO₂.

The BoM biogeochemical functioning is mainly forced by wind-driven hydrodynamics (upwelling, downwelling), urban rivers, wastewater treatment plants, and atmospheric deposition (Fraysse et al., 2013). In addition, Northern Current and Rhône River plume intrusions frequently occurred (Fraysse et al., 2014; Ross et al., 2016). Moreover, the BoM harbors the second biggest metropolis of France (Marseille), which is impacted by many harbor activities. The next step of this study will be to couple the Eco3M-CarbOx biogeochemical model to a 3-D hydrodynamic model that will mirror the complexity of the BoM functioning. In this way, the contributions of hydrodynamic, atmospheric, anthropic, and biogeochemical processes to the DIC variability will be able to be determined with higher refinement and realism, and an overview of the air–sea CO₂ exchange could be made at the scale of the Bay of Marseille. The main results of our study could be transposed to other coastal sites that are also impacted by urban and anthropic pressures. Moreover, in this paper we highlighted that fast and strong variations of pCO₂ values occur, and thus it is essential to acquire more in situ values at high frequency (at least with an hourly resolution) to understand the rapid variations of the marine carbon system at these short spatial and temporal scales.

Appendix A: Details of resolution of carbonate system module

A1 Calculation of carbonate system constants

The concentrations in conservative elements (relative to salinity) are calculated as follows:

• Total fluoride (TF) concentrations from Riley (1965) in mol kg⁻¹:

  $$\text{TF}={\displaystyle \frac{\mathrm{0.000067}}{\mathrm{18.998}}}\cdot {\displaystyle \frac{S}{\mathrm{1.80655}}}.$$

• Total sulfate (TS) concentration from Morris and Riley (1966) in mol kg⁻¹:

  $$\text{TS}={\displaystyle \frac{\mathrm{0.14}}{\mathrm{96.062}}}\cdot {\displaystyle \frac{S}{\mathrm{1.80655}}}.$$

• Calcium ion concentration from Riley and Tongudai (1967) in mol kg⁻¹:

  $${\mathrm{Ca}}^{\mathrm{2}+}={\displaystyle \frac{\mathrm{0.02128}}{\mathrm{40.087}}}\cdot {\displaystyle \frac{S}{\mathrm{1.80655}}}.$$

• Total boron (TB) concentration from Uppström (1974) in mol kg⁻¹:

  $$\text{TB}={\displaystyle \frac{\mathrm{0.000416}\cdot S}{\mathrm{35}}}.$$

• Ionic strength (IonS) from Millero (1982):

  $$\text{IonS}={\displaystyle \frac{\mathrm{19.924}\cdot S}{\mathrm{1000}-\mathrm{1.005}\cdot S}}.$$

The constants are calculated on the total pH scale except for K_S on free pH scale. If necessary, pH scale conversion factors are as follows:

• From seawater pH scale (SWS) to total pH scale: $\text{SWStoTOT}=\frac{\mathrm{1}+\frac{\text{TS}}{K_S}}{\mathrm{1}+\frac{\text{TS}}{K_S}+\frac{\text{TF}}{K_F}}$;

• From free pH scale to total pH scale: $\text{FREEtoTOT}=\mathrm{1}+\frac{\text{TS}}{K_S}$.

Further carbonate system constants were calculated as follows:

• K_S equilibrium constant of dissociation of ${\mathrm{HSO}}_{\mathrm{4}}^{-}$ from Dickson (1990a) in mol kg⁻¹:

  $$\begin{array}{ll}{\displaystyle }{K}_S=& {\displaystyle }{\displaystyle \frac{-\mathrm{4276.1}}{T_{\left(K\right)}}}+\mathrm{141.328}-\mathrm{23.093}\cdot \log \left(T_{\left(K\right)}\right)\\ {\displaystyle }& {\displaystyle }+\left(\mathrm{324.57}-\mathrm{47.986}\cdot \log \left(T_{\left(K\right)}\right)-{\displaystyle \frac{\mathrm{13856}}{T_{\left(K\right)}}}\right)\\ {\displaystyle }& {\displaystyle }\cdot {\text{Ions}}^{\mathrm{2}},\\ {\displaystyle }{K}_S=& {\displaystyle }{K}_S+(-\mathrm{771.54}+\mathrm{114.723}\cdot \log (T_{\left(K\right)})\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{35}\mathrm{474}}{T_{\left(K\right)}}})\cdot \text{Ions}+{\displaystyle \frac{-\mathrm{2698}}{T_{\left(K\right)}}}\cdot {\text{Ions}}^{\frac{\mathrm{3}}{\mathrm{2}}}+{\displaystyle \frac{\mathrm{1776}}{T_{\left(K\right)}}}\\ {\displaystyle }& {\displaystyle }\cdot {\text{Ions}}^{\mathrm{2}},\\ {\displaystyle }{K}_S=& {\displaystyle }{e}^{K_S\cdot (\mathrm{1}-\mathrm{0.001005}\cdot S)};\end{array}$$

• K_F equilibrium constant of dissociation of hydrogen fluoride (HF) formation from Dickson and Riley (1979) in mol kg⁻¹:

  $$K_{\mathrm{F}}=e^{\frac{\mathrm{1590.2}}{T_{\left(K\right)}}-\mathrm{12.641}+\mathrm{1.525}\cdot {\text{Ions}}^{\frac{\mathrm{1}}{\mathrm{2}}}}\cdot \left(\mathrm{1}-\mathrm{0.001005}\cdot S\right);$$

• K_B equilibrium constant of dissociation of boric acid from Dickson (1990b) in mol kg⁻¹:

  $$\begin{array}{ll}{\displaystyle }{K}_{\mathrm{B}}=& {\displaystyle }(-\mathrm{8966.9}-\mathrm{2890.53}\cdot S^{\frac{\mathrm{1}}{\mathrm{2}}}-\mathrm{77.942}\cdot S\\ {\displaystyle }& {\displaystyle }+\mathrm{1.728}\cdot S^{\frac{\mathrm{3}}{\mathrm{2}}}-\mathrm{0.0996}\cdot S^{\mathrm{2}})/T_{\left(K\right)},\\ {\displaystyle }{K}_{\mathrm{B}}=& {\displaystyle }{K}_{\mathrm{B}}+\mathrm{148.0248}+\mathrm{137.1942}\cdot S^{\frac{\mathrm{1}}{\mathrm{2}}}+\mathrm{1.62142}\cdot S\\ {\displaystyle }& {\displaystyle }+(-\mathrm{24.4344}-\mathrm{25.085}\cdot S^{\frac{\mathrm{1}}{\mathrm{2}}}-\mathrm{0.2474}\cdot S)\\ {\displaystyle }& {\displaystyle }\cdot \log \left(T\right)+\mathrm{0.053105}\cdot S^{\frac{\mathrm{1}}{\mathrm{2}}}\cdot T;\end{array}$$

• K₀ constant of CO₂ solubility from Weiss (1974) in $\mathrm{mol}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{atm}}^{-\mathrm{1}}$:

  $$\begin{array}{ll}{\displaystyle }{K}_{\mathrm{0}}=& {\displaystyle }\exp (-\mathrm{60.2409}+\mathrm{93.4517}\cdot {\displaystyle \frac{\mathrm{100}}{T_{\left(K\right)}}}+\mathrm{23.3585}\\ {\displaystyle }& {\displaystyle }\cdot \log \left({\displaystyle \frac{T_{\left(K\right)}}{\mathrm{100}}}\right)+S\cdot (\mathrm{0.023517}-\mathrm{0.023656}\cdot {\displaystyle \frac{T_{\left(K\right)}}{\mathrm{100}}}\\ {\displaystyle }& {\displaystyle }+\mathrm{0.0047036}\cdot {\left({\displaystyle \frac{T_{\left(K\right)}}{\mathrm{100}}}\right)}^{\mathrm{2}}\left)\right);\end{array}$$

• K_e dissociation constant of water from Millero (1995) in (mol kg⁻¹)²:

  $$\begin{array}{ll}{\displaystyle }{K}_e=& {\displaystyle }\exp ({\displaystyle \frac{-\mathrm{13847.26}}{T_{\left(K\right)}}}+\mathrm{148.9802}-\mathrm{23.6521}\\ {\displaystyle }& {\displaystyle }\cdot \log \left(T_{\left(K\right)}\right)+(-\mathrm{5.977}+{\displaystyle \frac{\mathrm{118.67}}{T_{\left(K\right)}}}+\mathrm{1.0495}\\ {\displaystyle }& {\displaystyle }\cdot \log \left(T_{\left(K\right)}\right))\cdot S^{\frac{\mathrm{1}}{\mathrm{2}}}-\mathrm{0.01615}\cdot S),\\ {\displaystyle }{K}_e=& {\displaystyle }{K}_e\cdot \text{SWStoTOT, on total }p\mathrm{H}\text{ scale in }\mathrm{mol}{\mathrm{kg}}^{-\mathrm{1}};\end{array}$$

• K₁ and K₂ from Lueker et al. (2000) in mol kg⁻¹:

  $$\begin{array}{ll}{\displaystyle }{K}_{\mathrm{1}}=& {\displaystyle }{\mathrm{10}}^{\left(\begin{array}{l}\frac{-\mathrm{3633.86}}{T_{\left(K\right)}}+\mathrm{61.2172}-\mathrm{9.6777}\cdot \log \left(T_{\left(K\right)}\right)\\ +\mathrm{0.011555}\cdot S-\mathrm{0.0001152}\cdot S^{\mathrm{2}}\end{array}\right)},\\ {\displaystyle }{K}_{\mathrm{2}}=& {\displaystyle }{\mathrm{10}}^{\left(\begin{array}{l}\frac{-\mathrm{471.78}}{T_{\left(K\right)}}+\mathrm{251.929}-\mathrm{3.16967}\cdot \log \left(T_{\left(K\right)}\right)\\ +\mathrm{0.01781}\cdot S-\mathrm{0.0001122}\cdot S^{\mathrm{2}}\end{array}\right)};\end{array}$$

• K_ca for calcite from Mucci (1983) in (mol kg⁻¹)²:

  $$K_{\text{ca}}={\mathrm{10}}^{\left(\begin{array}{l}-\mathrm{171.9065}-\mathrm{0.077993}\cdot T_{\left(K\right)}+\frac{\mathrm{2839.319}}{T_{\left(K\right)}}\\ +\mathrm{71.595}\cdot \log \mathrm{10}\left(T_{\left(K\right)}\right)+(-\mathrm{0.77712}\\ +\mathrm{0.0028426}\cdot T_{\left(K\right)})+\frac{\mathrm{178.34}}{T_{\left(K\right)}}\cdot S^{\frac{\mathrm{1}}{\mathrm{2}}}\\ -\mathrm{0.07711}\cdot S\\ +\mathrm{0.0041249}\cdot S^{\frac{\mathrm{3}}{\mathrm{2}}}\end{array}\right)}.$$

All the constants are corrected by the effect of hydrostatic pressure:

R=83.1451 in $\mathrm{ml}{\mathrm{bar}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{mol}}^{-\mathrm{1}}$;

$$\begin{array}{ll}{\displaystyle }\ln K_{\mathrm{1}}\text{fac}=& {\displaystyle }{\displaystyle \frac{\begin{array}{l}(\mathrm{25.5}-\mathrm{0.1271}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}+\mathrm{0.5}\\ \cdot \left(\frac{-\mathrm{3.08}+\mathrm{0.0877}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}}{\mathrm{1000}}\right)\cdot P_{\text{bar}})\cdot P_{\text{bar}}\end{array}}{R\cdot T_{\left(K\right)}}},\\ {\displaystyle }& {\displaystyle }{K}_{\mathrm{1}}=K_{\mathrm{1}}\cdot e^{\ln K_{\mathrm{1}}\text{fac}};\\ {\displaystyle }\ln K_{\mathrm{2}}\text{fac}=& {\displaystyle }{\displaystyle \frac{\begin{array}{l}(\mathrm{15.82}-\mathrm{0.0219}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}+\mathrm{0.5}\\ \cdot \left(\frac{\mathrm{1.13}+\mathrm{0.1475}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}}{\mathrm{1000}}\right)\cdot P_{\text{bar}})\cdot P_{\text{bar}}\end{array}}{R\cdot T_{\left(K\right)}}},\\ {\displaystyle }& {\displaystyle }{K}_{\mathrm{2}}=K_{\mathrm{2}}\cdot e^{\ln K_{\mathrm{2}}\text{fac}};\\ {\displaystyle }\ln K_{\mathrm{B}}\text{fac}=& {\displaystyle }{\displaystyle \frac{\begin{array}{l}(\mathrm{29.48}-\mathrm{0.1622}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}+\mathrm{0.002608}\\ \cdot T_{\left({}^{\circ }\mathrm{C}\right)}^{\mathrm{2}}+\mathrm{0.5}\cdot \left(\frac{-\mathrm{2.84}}{\mathrm{1000}}\right)\cdot P_{\text{bar}})\cdot P_{\text{bar}}\end{array}}{R\cdot T_{\left(K\right)}}},\\ {\displaystyle }& {\displaystyle }{K}_{\mathrm{B}}=K_{\mathrm{B}}\cdot e^{\ln K_{\mathrm{B}}\text{fac}};\\ {\displaystyle }\ln K_e\text{fac}=& {\displaystyle }{\displaystyle \frac{\begin{array}{l}(\mathrm{20.02}-\mathrm{0.1119}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}+\mathrm{0.001409}\\ \cdot T_{\left({}^{\circ }\mathrm{C}\right)}^{\mathrm{2}}+\mathrm{0.5}\cdot \left(\frac{-\mathrm{5.13}+\mathrm{0.0794}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}}{\mathrm{1000}}\right)\cdot P_{\text{bar}})\\ \cdot P_{\text{bar}}\end{array}}{R\cdot T_{\left(K\right)}}},\\ {\displaystyle }& {\displaystyle }{K}_e=K_e\cdot e^{\ln K_e\text{fac}};\\ {\displaystyle }\ln K_{\mathrm{F}}\text{fac}=& {\displaystyle }{\displaystyle \frac{\begin{array}{l}(\mathrm{9.78}-\mathrm{0.009}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}+\mathrm{0.0009429}\\ \cdot T_{\left({}^{\circ }\mathrm{C}\right)}^{\mathrm{2}}+\mathrm{0.5}\cdot \left(\frac{-\mathrm{3.91}+\mathrm{0.054}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}}{\mathrm{1000}}\right)\cdot P_{\text{bar}})\\ \cdot P_{\text{bar}}\end{array}}{R\cdot T_{\left(K\right)}}},\\ {\displaystyle }& {\displaystyle }{K}_{\mathrm{F}}=K_{\mathrm{F}}\cdot e^{\ln K_{\mathrm{F}}\text{fac}};\\ {\displaystyle }\ln K_S\text{fac}=& {\displaystyle }{\displaystyle \frac{\begin{array}{l}(\mathrm{18.03}-\mathrm{0.0466}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}+\mathrm{0.000316}\\ \cdot T_{\left({}^{\circ }\mathrm{C}\right)}^{\mathrm{2}}+\mathrm{0.5}\cdot \left(\frac{-\mathrm{4.53}+\mathrm{0.009}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}}{\mathrm{1000}}\right)\cdot P_{\text{bar}})\\ \cdot P_{\text{bar}}\end{array}}{R\cdot T_{\left(K\right)}}},\\ {\displaystyle }& {\displaystyle }{K}_S=K_S\cdot e^{\ln K_S\text{fac}}.\end{array}$$

$$\begin{array}{ll}{\displaystyle }\ln K_{\text{ca}}\text{fac}=& {\displaystyle }{\displaystyle \frac{\begin{array}{l}(\mathrm{48.76}-\mathrm{0.5304}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}\\ +\mathrm{0.5}\cdot \left(\frac{-\mathrm{11.76}+\mathrm{0.3692}\cdot T_{\left({}^{\circ }\mathrm{C}\right)}}{\mathrm{1000}}\right)\cdot P_{\text{bar}})\cdot P_{\text{bar}}\end{array}}{R\cdot T_{\left(K\right)}}};\\ {\displaystyle }& {\displaystyle }{K}_{\text{ca}}=K_{\text{ca}}\cdot e^{\ln K_{\text{ca}}\text{fac}}\end{array}$$

Calculation of the fugacity factor

We suppose that the pressure is at one atmosphere or close to it (Weiss, 1974):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{P}_{\text{atm}}=\mathrm{1.01325}\mathrm{bar},\\ {\displaystyle }& {\displaystyle }\text{delta}=\mathrm{57.7}-\mathrm{0.118}\cdot T\text{ in }{\mathrm{cm}}^{\mathrm{3}}{\mathrm{mol}}^{-\mathrm{1}},\\ {\displaystyle }& {\displaystyle }b=-\mathrm{1636.75}+\mathrm{12.0408}\cdot T-\mathrm{0.0327957}\cdot T^{\mathrm{2}},\\ {\displaystyle }& {\displaystyle }+\mathrm{3.16528}\cdot \mathrm{0.00001}\cdot T^{\mathrm{3}}\text{ in }{\mathrm{cm}}^{\mathrm{3}}{\mathrm{mol}}^{-\mathrm{1}},\\ {\displaystyle }& {\displaystyle }\text{FugFac}={\exp }^{\frac{(b+\mathrm{2}\cdot \text{delta})\cdot P_{\text{atm}}}{R\cdot T}}.\end{array}$$

A2 Resolution of carbonate system

To resolve the carbonate system, we calculate the deltapH, which is the difference of pH between two iterations of the model. We initialize the run by imposing a pH value of 8.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{if}(\text{nbiter}<\mathrm{1})\text{pH}=\mathrm{8}\\ {\displaystyle }& {\displaystyle }\text{pHtol}=\mathrm{0.001}\text{! tolerance for iterations end}\\ {\displaystyle }& {\displaystyle }\text{deltapH}=\text{pHtol}+\mathrm{1}\\ {\displaystyle }& {\displaystyle }\text{do while (abs(deltapH)}>\mathrm{0.0001})\\ {\displaystyle }& {\displaystyle }H={\mathrm{10}}^{-\text{pH}}\\ {\displaystyle }& {\displaystyle }\text{Denom}=H^{\mathrm{2}}+K_{\mathrm{1}}\cdot H+K_{\mathrm{1}}\cdot K_{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }\text{CAlk}=\text{DIC}\cdot K_{\mathrm{1}}\cdot \left({\displaystyle \frac{H+\mathrm{2}\cdot K_{\mathrm{2}}}{\text{Denom}}}\right)\\ {\displaystyle }& {\displaystyle }\text{BAlk}={\displaystyle \frac{\text{TB}\cdot K_{\mathrm{B}}}{K_{\mathrm{B}}+H}}\\ {\displaystyle }& {\displaystyle }\text{OH}={\displaystyle \frac{K_e}{H}}\\ {\displaystyle }& {\displaystyle }\text{FreetoTot}=\mathrm{1}+{\displaystyle \frac{\text{TS}}{K_S}}\\ {\displaystyle }& {\displaystyle }\text{Hfree}={\displaystyle \frac{H}{\text{FreetoTot}}}\\ {\displaystyle }& {\displaystyle }{\mathrm{HSO}}_{\mathrm{4}}={\displaystyle \frac{\text{TS}}{\mathrm{1}+\frac{K_S}{\text{Hfree}}}}\\ {\displaystyle }& {\displaystyle }\text{HF}={\displaystyle \frac{\text{TF}}{\mathrm{1}+\frac{K_{\mathrm{F}}}{\text{Hfree}}}}\\ {\displaystyle }& {\displaystyle }\text{Residual}=\text{TA}-\text{CAlk}-\text{BAlk}-\mathrm{OH}+\text{Hfree}+{\mathrm{HSO}}_{\mathrm{4}}\\ {\displaystyle }& {\displaystyle }+\text{HF}\\ {\displaystyle }& {\displaystyle }\text{Slope}=\text{DIC}\cdot H\cdot K_{\mathrm{1}}\cdot (H^{\mathrm{2}}+K_{\mathrm{1}}\cdot K_{\mathrm{2}}+\mathrm{4}\cdot H\cdot K_{\mathrm{2}})\\ {\displaystyle }& {\displaystyle }\text{Slope}={\displaystyle \frac{\text{Slope}}{{\text{Denom}}^{\mathrm{2}}}}+\mathrm{OH}+H+{\displaystyle \frac{\text{BAlk}\cdot H}{K_{\mathrm{B}}+H}}\\ {\displaystyle }& {\displaystyle }\text{Slope}=\log \mathrm{10}\cdot \text{Slope}\\ {\displaystyle }& {\displaystyle }\text{deltapH}=\text{Residual}/\text{Slope}\text{! this is Newton's method}\\ {\displaystyle }& {\displaystyle }\text{do while }\left(\text{abs}\right(\text{deltapH})>\mathrm{1})\text{deltapH}={\displaystyle \frac{\text{deltapH}}{\mathrm{2}}}\\ {\displaystyle }& {\displaystyle }\text{! to keep the jump from being too big}\\ {\displaystyle }& {\displaystyle }\text{enddo}\\ {\displaystyle }& {\displaystyle }\text{pH}=\text{pH}+\text{deltapH}\text{! Is on the same scale as}{K}_{\mathrm{1}}\text{and}{K}_{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }\text{were calculated, i.e., total}p\mathrm{H}\text{scale}\\ {\displaystyle }& {\displaystyle }p{\mathrm{CO}}_{\mathrm{2}}=\left({\displaystyle \frac{\text{DIC}\cdot H^{\mathrm{2}}}{H^{\mathrm{2}}+K_{\mathrm{1}}\cdot H+K_{\mathrm{1}}\cdot K_{\mathrm{2}}}}\right)\cdot {\displaystyle \frac{{\mathrm{10}}^{\mathrm{6}}}{K_{\mathrm{0}}\cdot \text{FugFac}}}\\ {\displaystyle }& {\displaystyle }\text{! in }\mathrm{\mu }\mathrm{atm}\\ {\displaystyle }& {\displaystyle }{\mathrm{CO}}_{\mathrm{2}}={\displaystyle \frac{\text{DIC}\cdot {\mathrm{10}}^{\mathrm{6}}}{\mathrm{1}+\frac{K_{\mathrm{1}}}{H}+\frac{K_{\mathrm{1}}\cdot K_{\mathrm{2}}}{H^{\mathrm{2}}}}}\\ {\displaystyle }& {\displaystyle }{\mathrm{HCO}}_{\mathrm{3}}={\displaystyle \frac{K_{\mathrm{1}}\cdot {\mathrm{CO}}_{\mathrm{2}}}{H}}\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{CO}}_{\mathrm{3}}={\displaystyle \frac{K_{\mathrm{2}}\cdot {\mathrm{HCO}}_{\mathrm{3}}}{H}}\\ {\displaystyle }& {\displaystyle }\text{Omega}={\displaystyle \frac{\mathrm{Ca}\cdot {\mathrm{CO}}_{\mathrm{3}}\cdot {\mathrm{10}}^{-\mathrm{6}}}{K_{\text{ca}}}}\end{array}$$

Appendix B: Biogeochemical model variables and parameters

Table B1 Initial conditions of the state variables of the Eco3M-CarbOx model.

^* Diagnostic variables.

Table B2 Balance equations of the Eco3M-CarbOx model.

Table B3 Biogeochemical processes simulated by the Eco3M-CarbOx model.

Table B4 Value of parameters.

Appendix C: Short user manual

After uploading the whole archive on the Zenodo site (Lajaunie-Salla et al., 2020; https://doi.org/10.5281/zenodo.3757677), the exact version of the Eco3M-CarbOx code used in this study can be run as follows:

make !two executable files will be created: eco3M_ini.exe and eco3M.exe

Please note the following:

• The file config.ini allows the following to be defined: the time, time step, and save time of simulation variables biogeochemical process.

• Results files are stocked in “SORTIES” directory.

• Boundary conditions and forcings data are stocked in “DATA” directory.

• All subroutines of biogeochemical processes are stocked in “F_PROCESS” directory.

For further information, please contact Frédéric Diaz (frederic.diaz@univ-amu.fr) or Christel Pinazo (christel.pinazo@univ-amu.fr).