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 CO2. 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 CO2 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 CO2 to the atmosphere to a sink into the ocean. Moreover, the higher the wind speed is, the higher the air– sea CO2 gas exchange fluxes are. The river intrusions with nitrate supplies lead to a decrease in the pCO2 value, favoring the conditions of a sink for atmospheric CO2 into the BoM. A scenario of high atmospheric concentrations of CO2 also favors the conditions of a sink for atmospheric CO2 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.


Abstract.
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 2 .
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 2 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 pertur-bation. When the seawater temperature changes quickly, the behavior of the BoM waters alters within a few days from a source of CO 2 to the atmosphere to a sink into the ocean. Moreover, the higher the wind speed is, the higher the airsea CO 2 gas exchange fluxes are. The river intrusions with nitrate supplies lead to a decrease in the pCO 2 value, favoring the conditions of a sink for atmospheric CO 2 into the BoM. A scenario of high atmospheric concentrations of CO 2 also favors the conditions of a sink for atmospheric CO 2 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.

Introduction
Current climate change mostly originates from the carbon dioxide (CO 2 ) 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 2 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 2 , 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).
Published by Copernicus Publications on behalf of the European Geosciences Union. 296 K. Lajaunie-Salla et al.: A carbonate system model in a Mediterranean coastal site CO 2 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 2 . In the ocean, the main processes regulating CO 2 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 2 , with an estimated flux converging to 0.2 Pg C yr −1 (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 2 resulting from the increasing human pressure also remains poorly known. There are few works which highlight, under future atmospheric CO 2 levels, whether coastal zones will become a net sink or a reduced source of CO 2 (Andersson and Mackenzie, 2012;Cai, 2011).
The current increase in the atmospheric CO 2 partial pressure (pCO 2 ) is slowly shifting the marine carbonate chemistry equilibrium towards increases in the seawater pCO 2 and bicarbonate ions (HCO − 3 ) and decreases in pH and carbon-ate ions (CO 2− 3 ) (Hoegh-Guldberg et al., 2018). These trends were already described in several coastal and open-ocean locations worldwide . In a coastal northwestern Mediterranean site, a 10-year time series of in situ measurements highlights a trend of pH decrease and pCO 2 increase (Kapsenberg et al., 2017). Low pH values can inhibit the ability of many marine organisms to form the calcium carbonate (CaCO 3 ) used in the making of skeletons and shells (Gattuso et al., 2015). In an extreme case, this shift may promote dissolution of CaCO 3 because the water will become under-saturated with respect to CaCO 3 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 2 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 2 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(Fraysse et al., , 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.

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  Fraysse et al. (2013). Arrows represent processes between two state variables. TA: total alkalinity; DIC: dissolved inorganic carbon; CO 2 : dissolved carbon dioxide; O 2 : dissolved oxygen; CaCO 3 : calcium carbonate. 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 N Q (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 3 ) implicitly representing calcifying organisms. The knowledge of DIC and TA allows the calculation of the pCO 2 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 2 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 3 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 3 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 3 dissolution processes. Moreover, the dynamics of DIC are altered by the CO 2 exchanges with the atmosphere (Table B2, Appendix B). The air-sea CO 2 fluxes are calculated from the pCO 2 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.

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 3 s −1 (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, northnorthwesterly 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 2 ) from the nearby urban area, and it also receives effluents from the Aix-Marseille metropolis.

Dataset
The modeled variables of the carbonate system (DIC, TA, pH, and pCO 2 ) 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 2 , DIC, and TA), sampled every 2 weeks, have been available since 2016.

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 2 ) 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 2 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 sta- tion (Yohia, 2017). Finally, we used hourly atmospheric CO 2 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 2 variations control the seawater CO 2 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 −1 (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 −1 ) 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 −1 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 −3 ) = −1.70 × S + 65.
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.

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 2 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).
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 −3 in March and the second one at ca.  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 −3 ), seawater partial pressure of CO 2 (seawater pCO 2 , µatm), pH, dissolved inorganic carbon (DIC, µmol kg −1 ), and total alkalinity (TA, µmol kg −1 ). 0.68 mg m −3 in May. They are both linked to Rhône River plume intrusions. Several in situ maxima between 0.50 and 0.70 mg m −3 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 quan-titatively reproduces the spring bloom observed at the end of the MWC period ( Fig. 3a) with a maximum value of ca. 0.69 mg m −3 . 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 −3 , respectively -Table 2). On the whole, the seasonal variations of the seawater pCO 2 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 2 , and from February to March it reproduces the increase in pCO 2 even if the amplitude of this increase remains smaller. In mid-April, during the simulated spring bloom period, the observed drop in pCO 2 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 2 solubility due to thermal effects mostly explains the increase in pCO 2 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 2 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 −1 ; 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 −1 (Table 2) are observed, whatever the hydrological season being considered. The biogeochemical model provides almost constant values of around 2570 µmol kg −1 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).

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 2 fluxes (F aera ) vary between −14 and 17 mmol m −3 d −1 , with a weakly positive annual budget of +6 mmol m −3 yr −1 (or +0.017 mmol m −3 d −1 , Table 3). Then, the BoM waters would act as a net source of CO 2 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 mmol m −3 d −1 , respectively - Table 3), while the zooplankton and bacterial respiration account for an average of 0.348 mmol m −3 d −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 2 , which explains the DIC and seawater pCO 2 decreases during the bloom period  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 2 decreases and pH increases concomitantly ( Fig. 4c and d). Then, seawater appears to be temporarily under-saturated in CO 2 , and the BoM waters thus act as a sink for atmospheric CO 2 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 2 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 2 (from +17 to −14 mmol m −3 d −1 ; Fig. 4k), and also from a net heterotrophic to a net autotrophic ecosystem (Fig. 4l).
3.3 Impact of external forcing on the dynamics of carbonate system

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 2 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 −3 , 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 2 decreases, leading to some negative airsea fluxes (i.e., an oceanic sink for atmospheric CO 2 ). 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).

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 airsea fluxes are, mainly owing to the increase in gas transfer velocity. Depending on the gradient of CO 2 between seawater and the atmosphere, strong wind speeds will favor either the emission or uptake of CO 2 (Fig. 6b and e). In May, with the air-sea CO 2 flux being positive, the outgassing of CO 2 to the atmosphere is enhanced, leading to a decrease in seawater pCO 2 (Fig. 6c). On the contrary, in August the oceanic sink of atmospheric CO 2 is amplified, which leads to an increase in the seawater pCO 2 value (Fig. 6f).

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 2 (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 −1 ; Fig. 7j), due to the autotrophic processes dominating the heterotrophic ones. In turn, the seawater pCO 2 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 2 to the atmosphere and intensify the sink of atmospheric CO 2 into the waters of BoM ( Fig. 7e and k).

Urban air CO 2 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 2 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 2 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 2 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 2 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 2 values from urban or non-urban sites can lead to significant differences in the values of the seawater pCO 2 , especially during the MWC period. The air-sea gradient of pCO 2 is higher when using a forcing derived from the CO 2 concentrations originating from an urban area than from a non-urban area, which strengthens the sink of atmospheric CO 2 into the waters of BoM. The seawater pCO 2 is then lower with nonurban area pressure (S5) than with urban area pressure (S0), because of lower CO 2 solubility in the BoM (Fig. 8b).

Model performance
The evaluation of model skill vs. in situ data highlights that the modeled pH, pCO 2 , 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 2 , 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.

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 2 variations and associated fluxes would be mostly driven by the seawater temperature dynamics. Moreover, the seasonal variations of the air-sea CO 2 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 2 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 2 in seawater due to the decrease in the CO 2 solubility (Middelburg, 2019) and, at the same time, the fixation of DIC by phytoplankton is favored, leading to a de-crease in the pCO 2 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 2 during the bloom of phytoplankton ( Fig. 5e and j).

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 2 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 2 values (Fig. 4j) and, in fine, an alteration of the CO 2 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 eu- photic layer and would, in turn, contribute to the seawater pCO 2 decrease. North and northwestern winds through latent heat losses and/or upwelling events could then enhance the sink for atmospheric CO 2 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 −1 ) amplified the gaseous exchange of CO 2 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 sig-nificant impact on the seawaterpCO 2 values during a longer period of ca. 15 d (Fig. 6). A combination of high atmospheric pCO 2 values and high wind speeds would then favor the sink for CO 2 into the waters of the BoM. The aeration process also depends on the choice of the formulation of the gas transfer velocity (k 600 ). 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 2 owing to the nutrient supply favoring the phytoplankton development (Fig. 7). In this scenario, the oceanic sink of atmospheric CO 2 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 2 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 2 levels in the model. This drop of pCO 2 is due to the decrease in the CO 2 solubility when salinity decreases (Middelburg, 2019).
In the scenario of forcing the model by using urban atmospheric pCO 2 time series, the air-sea gradient increases, and then it enhances the status of the BoM as a sink for atmospheric CO 2 . 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 2 , the CO 2 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 2 sink and a potential threat to coastal marine biodiversity in coastal areas owing to high atmospheric CO 2 levels.

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 2 ). 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 2 and a dominantly autotrophic ecosystem during the MWC period, and a net source of CO 2 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 2 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 2 with the atmosphere. A Rhône River plume intrusion with input of nitrate favors pCO 2 decreases, and the sink of atmospheric CO 2 into the BoM waters is enhanced. The higher atmospheric pCO 2 values from the urban area intensify the oceanic sink of atmospheric CO 2 .
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 2 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 2 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. The concentrations in conservative elements (relative to salinity) are calculated as follows: -Total fluoride (TF) concentrations from Riley (1965) in mol kg −1 : -Total sulfate (TS) concentration from Morris and Riley (1966) in mol kg −1 : TS = 0.14 96.062 · S 1.80655 .
-Calcium ion concentration from Riley and Tongudai (1967) in mol kg −1 : -Total boron (TB) concentration from Uppström (1974) in mol kg −1 : -Ionic strength (IonS) from Millero (1982): 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: -From free pH scale to total pH scale: FREEtoTOT = 1 + TS K S .

Calculation of the fugacity factor
We suppose that the pressure is at one atmosphere or close to it (Weiss, 1974): 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.
if (nbiter < 1) pH = 8 pHtol = 0.001 ! tolerance for iterations end deltapH = pHtol + 1 do while (abs(deltapH) > 0.0001) BAlk · H K B + H Slope = log 10 · Slope deltapH = Residual/Slope ! this is Newton's method do while (abs(deltapH) > 1) deltapH = deltapH 2 ! to keep the jump from being too big enddo pH = pH + deltapH ! Is on the same scale as K 1 and K 2 were calculated, i.e., total pH scale pCO 2 = DIC · H 2 H 2 + K 1 · H + K 1 · K 2 · 10 6 K 0 · FugFac ! in µatm CO 2 = DIC · 10 6       More rapid attenuation of polychromatic light near the sea surface 1.0 1.0 - Tett (1990) α Chl a Chlorophyll-specific light absorption coefficient 8 × 10 −6 5 × 10 −6 m 2 mol C (g Chl a J) −1 Leblanc et al. (2018) T opt Temperature optimal of growth 16.0 13.0 • C -T let Lethal temperature 11.0 9.0 • C b Shape factor for temperature curve 0.5 0.8 - Lacroix and Grégoire (2002)     Author contributions. KLS conceptualized this study, made code developments on the carbonate module, designed the numerical experiments, performed the final analysis to the model results (validation, statistics), developed MATLAB software to visualize the results of model and statistics tests, wrote and edited the initial draft, revised the manuscript, and answered reviewers. FD participated in the conceptualization of the present study, supervised this research, wrote some parts of the initial draft, helped to answer reviewer comments, revised the draft after the second review, and edited the final version. CWR participated in the data acquisition of the seawater pCO 2 and corresponding measurements at the laboratory and reviewed and helped with the initial draft. TW participated in the data acquisition of the seawater pCO 2 and corresponding measurements at the laboratory, helped to answer reviews, and helped with the final editing. DL acquired funding, conducted the project administration, and participated in the data acquisition of the seawater pCO 2 . CY helped to maintain computing resources and gave expertise on the Eco3M-CarbOx code and software developments and compiling. IX-R acquired funding, conducted the project administration, and provided data on the atmospheric pCO 2 at OHP. BN participated in the acquisition of pCO 2 data at OHP and reviewed and edited the first submitted draft. AA provided temporal data on the atmospheric pCO 2 at CAV. CP participated in the conceptualization of the present study, provided the initial Eco3M code (without the carbonate module), supervised this research, reviewed the initial draft, and helped to answer reviewer comments.