Seasonal dynamics and annual budget of dissolved inorganic carbon in the northwestern Mediterranean deep-convection region

Deep convection plays a key role in the circulation, thermodynamics, and biogeochemical cycles in the Mediterranean Sea, which is considered to be a hotspot of biodiversity and climate change. In the framework of the DEWEX (Dense Water Experiment) project, the seasonal and annual budgets of dissolved inorganic carbon in the deep-convection area of the northwestern Mediterranean Sea are investigated over the period September 2012–September 2013 using a 3D coupled physical–biogeochemical–chemical modeling approach. At the annual scale, we estimate that the northwestern Mediterranean Sea's deep-convection region was a moderate sink of 0.5 mol C m⁻² yr⁻¹ of CO₂ for the atmosphere. The model results show the reduction of oceanic CO₂ uptake during deep convection and its increase during the abrupt spring phytoplankton bloom following the deep-convection events. We highlight the major roles in the annual dissolved inorganic carbon budget of both the biogeochemical and physical fluxes, which amount to −3.7 and 3.3 mol C m⁻² yr⁻¹, respectively, and are 1 order of magnitude higher than the air–sea CO₂ flux. The upper layer (from the surface to 150 m depth) of the northwestern deep-convection region gained dissolved inorganic carbon through vertical physical transport and, to a lesser extent, oceanic CO₂ uptake, and it lost dissolved inorganic carbon through lateral transport and biogeochemical fluxes. The region, covering 2.5 % of the Mediterranean, acted as a source of dissolved inorganic carbon for the surface and intermediate water masses of the Balearic Sea and southwestern Mediterranean Sea and could represent up to 22 % and 11 %, respectively, of the CO₂ exchanges with the Atlantic Ocean at the Strait of Gibraltar.

1 Introduction

Quantifying the ocean carbon pump and its evolution under ongoing global warming and a rising atmospheric CO₂ inventory is a challenging issue. Exchanges of CO₂ at the air–sea interface result from a complex interplay of chemical, biogeochemical, and physical processes in the ocean. Physical mechanisms can quantitatively play a comparable role as biogeochemical processes in the air–sea CO₂ flux at regional and global scales (Ayers and Lozier, 2012; Lévy et al., 2013; Stukel and Ducklow, 2017). In particular, deep-convection regions, such as the Labrador Sea, located at high latitudes in the Atlantic Ocean, are considered to be large sinks for atmospheric CO₂ due to strong cooling and high primary production leading to long, or even persistent, periods of deficit compared to the atmosphere (Takahashi et al., 2002). In these regions, large amounts of atmospheric CO₂ taken up at the ocean surface and biologically fixed carbon are transferred to the deep ocean during the intense vertical mixing periods (DeGranpre et al., 2006; Körtzinger et al., 2008a). Furthermore, respired organic carbon remaining above the winter mixing depth can be ventilated back to the surface during the following winter (DeGranpre et al., 2006; Körtzinger et al., 2008a; Palevsky and Nicholson, 2018). Lateral transport, often associated with the restratification of the water column, the dispersion of the newly formed dense water, and/or the exchanges with boundary currents, also greatly contributes to the budget of water masses and their biogeochemical contents (Wolf et al., 2018; Koelling et al., 2022).

Figure 1 Bathymetry (m) of the study area. The red contour indicates the limit of the deep-convection area defined for the budget calculation (see Sect. 2.1.3). The location of the BOUSSOLE buoy and EMSO-DYFAMED mooring sites in the Ligurian Sea are indicated with a green and black star, respectively. The insert representing the Mediterranean Sea indicates the limits of the coupled model used for this study in blue and of the forcing simulation in orange (see Sect. 2.1.2).

The northwestern Mediterranean Sea (Gulf of Lion and Ligurian Sea, Fig. 1), along with the South Adriatic, is one of the regions where deep convection occurs (Ovchinnikov et al., 1985; Mertens and Schott, 1998; Manca and Bregant, 1998; Béthoux et al., 2002). Few studies have investigated the dynamics of dissolved inorganic carbon (DIC hereafter) in this region, where the Western Mediterranean Deep Water is formed and which plays a crucial role in the circulation and ventilation of the Mediterranean Sea (Schroeder et al., 2016; Li and Tanhua, 2020; Mavropoulou et al., 2020). The objective of this study is to gain insights into the annual cycle of DIC by examining and quantifying the biogeochemical, physical, and air–sea fluxes.

In the northwestern Mediterranean region, a basin-scale cyclonic gyre is associated with a doming of isopycnals. The density increase, induced in winter in surface waters by cold and dry northerly winds, produces instabilities in the water column, leading to convective mixing of surface waters with deeper waters. With regards to the biogeochemical processes, the region is characterized at the sea surface by a moderate phytoplankton bloom in fall, interrupted by deep winter mixing and an abrupt phytoplankton bloom following deep winter mixing, which supplies inorganic nutrients to the euphotic layer (Severin et al., 2014; Bernardello et al., 2012; Lavigne et al., 2013; Ulses et al., 2016; Kessouri et al., 2017). At the annual scale, the net community production (NCP, defined as the gross primary production minus the community respiration) was found to be positive, leading to an autotrophic status for the area (Ulses et al., 2016; Coppola et al., 2018). The downward export of organic carbon and its interannual variability have been related to the intensity of the deep convection and the phytoplankton bloom (Heimbürger et al., 2013; Herrmann et al., 2013; Ulses et al., 2016).

Previous observational and modeling studies that have documented the dynamics of the CO₂ system in this region mostly focused on the Ligurian Sea at the EMSO-DYFAMED (European Multidisciplinary Seafloor and water column Observatory-Dynamique des Flux Atmospheriques en MEDiterranee) and BOUSSOLE (BOUée pour l'acquiSition d'une Série Optique à Long termE) mooring sites (Hood and Merlivat, 2001; Copin-Montégut and Bégovic, 2002; Bégovic and Copin-Montégut, 2002; Mémery et al., 2002; Copin-Montégut et al., 2004; Touratier and Goyet, 2009; Merlivat et al., 2018; Coppola et al., 2020), where the intensity of convection generally remains moderate compared to the Gulf of Lion (Fig. 1). These 1D studies showed a pronounced seasonal cycle of pCO₂ mostly controlled by the sea surface temperature. The thermal effect is counterbalanced in spring by the impact of phytoplankton growth, which leads to DIC drawdown, and in winter by intense mixing events, which bring DIC-rich water to the surface (Hood and Merlivat, 2001; Mémery et al., 2002; Copin-Montégut et al., 2004). On an annual timescale, the Ligurian Sea was found to be a medium to minor sink for atmospheric CO₂ (Hood and Merlivat, 2001; Mémery et al., 2002; Copin-Montégut et al., 2004; Merlivat et al., 2008). Based on cruise data, Touratier et al. (2016) complemented those mooring observations by describing the distribution of the carbonate system properties in the central region of the deep-convection region during two winter periods, during and just after the deep-convection event. The authors showed a rapid transfer of anthropogenic CO₂ to the ocean interior during the convection event and found an excess in CO₂ related to the atmosphere. Finally, D'Ortenzio et al. (2008) and Cossarini et al. (2021), based on 1D models and a 3D model, respectively, found that the whole deep-convection region is a major sink of atmospheric CO₂ in the open Mediterranean Sea. In the previous studies, the 3D dynamics of the CO₂ system over an annual cycle has never been specifically explored for the whole northwestern deep-convection region, and a complete DIC budget is still lacking for this region.

The aim of the DEWEX (Dense Water Experiment) project (Conan et al., 2018; Testor et al., 2018) was to investigate the deep-convection process and its impact on biogeochemical fluxes based on observational platforms and numerical models. In this framework, two research cruises were carried out in winter and spring of the year 2013, completing the MOOSE-GE (Mediterranean Ocean Observing System for the Environment-Grande Echelle) observational effort that has been performed each year during stratified periods since 2010. Due to extremely strong buoyancy loss, the 2012–2013 winter was characterized by intense deep-convection events and is considered to be one of the five most intense events over the period 1980–2013 (Somot et al., 2016; Herrmann et al., 2017; Coppola et al., 2018). Using a 3D coupled physical–biogeochemical simulation, Kessouri et al. (2018) estimated that the deep-convection region was characterized over the period 2012–2013 by a positive net community production and showed higher rates of export of organic carbon below the euphotic layer compared to the surroundings. They suggested that, due to high spatial and interannual variability and the dispersion of newly formed dense water to the southern Mediterranean Sea, a fraction of the exported carbon escapes a return to the surface layer during the following winter. Here, we took advantage of the in situ measurements from the DEWEX project, the MOOSE-GE program, and the BOUSSOLE and EMSO-DYFAMED mooring sites to implement and constrain a model of the dynamics of the CO₂ system and to complete the 3D coupled physical–biogeochemical modeling study on organic carbon by Kessouri et al. (2018). We (i) examined the seasonal cycle of DIC, (ii) estimated an annual carbon budget, and (iii) analyzed the contributions from air–sea CO₂ exchanges and biogeochemical and physical processes to the carbon budget.

2 Material and methods

2.1 The numerical model

2.1.1 The coupled hydrodynamic–biogeochemical–chemical model

The biogeochemical model Eco3M-S (Auger et al., 2011; Ulses et al., 2021) was forced offline by daily outputs (current velocities, turbulent diffusion coefficient, temperature, and salinity) of the 3D hydrodynamic model SYMPHONIE (Marsaleix et al., 2008). The SYMPHONIE model, a 3D primitive equation model, with a free surface and generalized sigma vertical coordinate, has been used to investigate open-sea convection (Herrmann et al., 2008; Estournel et al., 2016; Damien et al., 2017) and circulation in the northwestern Mediterranean Sea (Estournel et al., 2003; Ulses et al., 2008; Bouffard et al., 2008).

The biogeochemical model Eco3M-S is a multi-nutrient and multi-plankton functional-type model that simulates the dynamics of the pelagic planktonic ecosystem and the cycles of several biogenic elements (carbon, nitrogen, phosphorus, silicon, oxygen) (Auger et al., 2011; Many et al., 2021; Ulses et al., 2021). Particulate organic detritus and micro-phytoplankton have a constant settling velocity (1 m d⁻¹ for slow-sinking detritus and micro-phytoplankton and 90 m d⁻¹ for fast-sinking detritus based on the modeling studies of Levy et al. (1998), Lacroix and Grégoire (2002), and Raick et al. (2005) and on calibration tests (Auger et al., 2011; Kessouri et al., 2017)). The model has been used to study biogeochemical processes in the NW (northwestern) Mediterranean deep-convection area (Herrmann et al., 2013; Auger et al., 2014; Ulses et al., 2016, 2021; Kessouri et al., 2017, 2018) and in the whole Mediterranean Sea (Kessouri, 2015).

Figure 2 Scheme of the upgraded biogeochemical model Eco3M-S (redrawn from Ulses et al., 2021).

In previous versions of the model, particulate and dissolved organic carbon were considered, but the dynamics of the dissolved inorganic carbon was not described. To investigate the dynamics of the CO₂ system, the model was extended by implementing the carbonate chemistry model developed and described in detail by Soetaert et al. (2007) and applied by Raick-Blum (2005) in 1D in the northwestern Mediterranean Sea. The food web structure of the upgraded model and the biogeochemical processes interacting between compartments are schematically represented in Fig. 2. Two state variables were added in the upgraded version of the coupled model. The first added variable is the dissolved inorganic carbon concentration, the sum of the concentrations of the four carbon dioxide forms, namely dissolved carbonate dioxide, bicarbonate, carbonate ion, and carbonic acid. The rate of change of the concentration of DIC due to biogeochemical processes is governed by the following equation:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\displaystyle \frac{\partial \text{DIC}}{\partial t}}\mathrm{|}\text{bio}=\sum _{i=\mathrm{1}}^{\mathrm{3}}(-{\text{GPP}}_i+{\text{RespPhy}}_i)+{\sum }_{j=\mathrm{1}}^{\mathrm{3}}\\ & {\text{RespZoo}}_j+\text{RespBac},\end{array}\end{array}$$

where GPP_i and RespPhy_i are gross primary production and respiration, respectively, for phytoplankton size class i (size classes 1, 2, and 3 are pico-, nano-, and micro-phytoplankton, respectively, Fig. 2); RespZoo_j is the respiration flux for zooplankton size class j (size classes 1, 2, and 3 are nano-, micro-, and meso-zooplankton, respectively, Fig. 2), and RespBac is bacterial respiration. The second added state variable is the excess negative charge (denoted $\sum \left[-\right])$, which is the moles of negative charges over positive charges of the acid–base system (Table 2 in Soetaert et al., 2007). As in Soetaert et al. (2007), we use this excess negative charge instead of the total alkalinity, commonly measured for proton balance. Here, we assume that the uptake of ions is compensated for by the uptake or release of protons (electroneutrality) and that $\sum \left[-\right]$ is not impacted by changes in the concentrations of nitrate, phosphate, or ammonium, which is not the case for total alkalinity. The total alkalinity (TA) is then deduced from $\sum \left[-\right]:$

$$\begin{array}{}\text{(2)}& \text{TA}=\sum \left[-\right]+\sum {\mathrm{NH}}_{\mathrm{3}}-\sum {\text{NO}}_{\mathrm{3}}-\sum {\text{PO}}_{\mathrm{4}}.\end{array}$$

In this first study on DIC dynamics, we neglected calcium carbonate (CaCO₃) precipitation and dissolution. Schneider et al. (2007) indicated that the Mediterranean Sea is supersaturated with respect to calcite and aragonite throughout the water column and that calcium carbonate dissolution is thus not favored thermodynamically. Regarding the CaCO₃ precipitation, we are aware that future refinements will have to take this into account since neglecting it could lead notably to an overestimation of air–sea CO₂ flux. Sensitivity tests of this term were performed (see Sect. 2.1.4) and are presented in Sect. 5.

In this study, in the carbonate chemistry model, the dissociation equilibriums of carbonates, water, ammonium, phosphate, silicate, and borate were taken into account. The thermodynamic equilibrium constants of the carbonate system were calculated as a function of temperature, salinity, and pressure as in Millero (1995), with typographical correction from the CO2SYS program (Lewis and Wallace, 1998). In particular, carbonic acid dissociation constants are calculated in the same way as Mehrbach et al. (1973) constants, as refit by Dickson and Millero (1987).

The flux of CO₂ at the air–sea interface, CO₂flux, was calculated using the following equation:

$$\begin{array}{}\text{(3)}& {\text{CO}}_{\mathrm{2}}\text{flux}=\mathit{\rho }{\text{K}}_{\mathrm{0}}{\text{K}}_{\text{w}}\left(p{\text{CO}}_{\mathrm{2},\text{atm}}-p{\text{CO}}_{\mathrm{2},\text{sea}}\right),\end{array}$$

where pCO_2,atm and pCO_2,sea (µatm) are the atmospheric and sea surface partial pressure of CO₂, respectively; K₀ (mol kg⁻¹ atm⁻¹) is the solubility coefficient; K_w (m s⁻¹) is the gas transfer velocity; and ρ is the sea surface density (kg m⁻³). We calculated the solubility coefficient according to Weiss (1974) and the gas transfer velocity using the most often used parameterization of Wanninkhof et al. (1992), with a quadratic dependency on the wind speed 10 m above the sea. In addition, we performed sensitivity analyses using eight various parameterizations of the gas transfer velocity to estimate uncertainties of air–sea exchanges (see Sect. 2.1.4).

2.1.2 Model setup

The numerical domain covers most of the western Mediterranean Sea (blue contour on the insert in Fig. 1) using a curvilinear grid (Bentsen et al., 1999) with a horizontal resolution varying from 0.8 km in the north to 1.4 km in the south and with 40 vertical levels (Ulses et al., 2021). The implementation of the hydrodynamic simulation and the strategy of downscaling from the Mediterranean Basin to the western sub-basin scale in three stages (Fig. S1 in the Supplement) have been described in detail in Estournel et al. (2016) and Kessouri et al. (2017) and will be summarized here:

• The hydrodynamic simulation was implemented over the western Mediterranean sub-basin (delimited by blue lines in the inset of Fig. 1) and was initialized and forced at its lateral boundaries with daily hydrodynamic fields of the configuration PSY2V4R4 based on the NEMO ocean model at a resolution of 1/12^∘ over the Mediterranean Basin (delimited by orange lines in the inset of Fig. 1) by the Mercator Ocean International operational system (Lellouche et al., 2013). This simulation (step 1a, Fig. S1) was performed from 1 August 2012 to 31 October 2013.

• In parallel (step 1b, Fig. S1), the biogeochemical model was computed, in offline mode, at the Mediterranean Basin scale on the same $\mathrm{1}/{\mathrm{12}}^{\circ }$ NEMO grid using the same NEMO hydrodynamic fields as those used by the SYMPHONIE simulation in step 1a. This simulation was performed from 15 June 2011 to 15 November 2013. The carbonate system module in this configuration was initialized in June 2011 using mean values of dissolved inorganic carbon and total alkalinity observations carried out in 2011 from the Meteor M84/3 (Alvarez et al., 2014), CASCADE (CAscading, Surge, Convection, Advection and Downwelling Events; Touratier et al., 2016), and MOOSE-GE cruises (Testor et al., 2010), as well as from the EMSO-DYFAMED mooring (Coppola et al., 2021) and BOUSSOLE buoy (Golbol et al., 2020) sites, over bio-regions defined in Kessouri (2015) based on Lavezza et al. (2011). To deduce the excess negative charge from total alkalinity (Eq. 2), we also used the nutrient concentration data from the MEDAR/MEDATLAS database as in Kessouri et al. (2017). Recently, Davis and Goyet (2021) described a method based upon the property variability to precisely quantify the uncertainties at any point of an interpolated data field. This approach could be used in the near future to improve both the at-sea sampling strategy (Guglielmi et al., 2022a, b) and the accuracy of model initialization.

• In a second step (Fig. S1), the Eco3M-S biogeochemical model was implemented over the western Mediterranean sub-basin using the grid and the hydrodynamics fields of the aforementioned SYMPHONIE simulation (step 1a) in offline mode. This simulation was performed from 15 August 2012 to 30 September 2013. The initial state and lateral boundary conditions of the biogeochemical fields are provided by the biogeochemical simulation of the Mediterranean Basin of step 1b.

This nesting protocol ensures the coherence of the physical and biogeochemical fields at the open boundaries of the western Mediterranean model. For both biogeochemical simulations (steps 1b and 2), we prescribed at the river mouths the mean DIC concentration measured by Sempéré et al. (2000) for the Rhone River and climatological values according to Ludwig et al. (2010) and Schneider et al. (2007) at the other river mouths. To compute the gas transfer velocity, we used the 3 h wind speed, pressure, and humidity provided by the ECMWF model on a $\mathrm{1}/{\mathrm{8}}^{\circ }$ grid, consistently with the hydrodynamic simulation. The atmospheric pCO_2,atm was deduced from the flask-air measurements of mole fraction, measured monthly at the Lampedusa site (World Data Centre for Greenhouse Gases: https://gaw.kishou.go.jp/,last access: 5 August 2019, Lan et al., 2022). Fluxes of dissolved inorganic carbon at the sediment–sea interface were considered by coupling the pelagic model with a simplified version of the meta-model described in Soetaert et al. (2001).

2.1.3 Study area and computation of DIC balance

We computed DIC fluxes and the resulting variation in the DIC inventory for the whole deep-convection area. The deep-convection area was defined as the area that includes the model grid points where the mixed-layer depth exceeded 1000 m for at least 1 d of the study period based on Kessouri et al. (2017, 2018). This area covers 70 100 km². The budget was calculated for two vertical layers: the photic upper layer, where the photosynthesis process takes place, and the aphotic deeper layer. The base of the upper layer was set at 150 m based on the regional minimum value of the diffuse attenuation coefficient of light at 490 nm as derived from satellite observations (http://marine.copernicus.eu/, last access: 19 November 2017, product: OCEANCOLOUR_MED_OPTICS_L3_REP_OBSERVATIONS_009
_095) and following the studies by Lazzari et al. (2012) and Kessouri et al. (2018).

The variation of the DIC inventory in the upper layer between times t₁ and t₂ (Δ DIC I_upper) is equal to the sum of all DIC fluxes within the deep-convection area between t₁ and t₂:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& \mathrm{\Delta }\text{DIC}{I}_{\text{upper}}=\text{DIC}{I}_{\text{upper},t_{\mathrm{1}}}-\text{DIC}{I}_{\text{upper},t_{\mathrm{2}}}\\ & =\underset{t_{\mathrm{1}}}{\overset{t_{\mathrm{2}}}{\int }}\left(F_{\text{DIC,air–sea}}+F_{\text{DIC,lat}}+F_{\text{DIC,vert}}+F_{\text{DIC,bgc}}\right)\mathrm{d}t,\end{array}\end{array}$$

where F_DIC,lat and F_DIC,vert are the lateral and vertical fluxes at the boundaries of the deep-convection area, F_{DIC,air–sea} is the air–sea CO₂ flux, and F_DIC,bgc is the biogeochemical flux.

DIC I_upper,t was computed from the following:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\text{DIC}}_{\text{upper},t}=\int \int \underset{\left(x,y\right)\mathit{ϵ}\text{DCA}/z\mathit{ϵ}\text{upper layer}}{\int }\\ & \text{DIC}\left(x,y,z,t\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z,\end{array}\end{array}$$

where ($x,y,z$) belongs to the upper layer (150 m to the surface) of the DCA (deep-convection area).

The lateral exchange flux was computed from the following:

$$\begin{array}{}\text{(6)}& F_{\text{DIC,lat}}=\int \underset{(x,y,z)\mathit{ϵ}A}{\int }\text{DIC}\left(x,y,z,t\right)v_t\left(x,y,z,t\right)\mathrm{d}A,\end{array}$$

where v_t is the current velocity normal to the limit of the deep-convection area (m s⁻¹), and A (m²) is the area of the section from the base of the upper layer (150 m) to the surface of the deep-convection area.

The F_DIC,air-sea was computed from

$$\begin{array}{}\text{(7)}& F_{\text{DIC,air-sea}}=\int \underset{(x,y)\mathit{ϵ}\text{DCA}}{\int }{\text{CO}}_{\mathrm{2}}\text{flux}\left(x,y,t\right){\mathrm{10}}^{-\mathrm{3}}\mathrm{d}x\mathrm{d}y,\end{array}$$

where CO₂flux (µmol C m⁻² s⁻¹) is the air–sea flux given by Eq. (3).

F_DIC,bgc was computed from the following:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& F_{\text{DIC,bgc}}=\int \int \underset{\left(x,y\right)\mathit{ϵ}DCA/z\mathit{ϵ}\text{upper layer}}{\int }\\ & \text{BGCflux}\left(x,y,z,t\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z,\end{array}\end{array}$$

where BGCflux is the biogeochemical flux, i.e., the sum of DIC release through respiration by living organisms and of DIC consumption through photosynthesis.

Finally, the vertical transport flux, F_DIC,vert, was derived from all the other terms of Eq. (4). The computation of the DIC balance in the deeper layer is computed in a similar way, with the variation of the inventory being the sum of the lateral and vertical flux at the boundaries and of the biogeochemical flux. Here, the fluxes at the sea–sediment interface were taken into account but were negligible with respect to the other terms of the balance.

2.1.4 Sensitivity tests

We performed various sensitivity tests to estimate the uncertainties of the modeled air–sea CO₂ flux. A first set of tests was based on the parametrization of the gas transfer coefficient. For these tests, we used quadratic (Wanninkhof, 2014), cubic (Wanninkhof and McGillis, 1999), and hybrid (Liss and Merlivat, 1986; Nightingale et al., 2000; Wanninkhof et al., 2009) wind speed dependency parameterizations of diffusive flux, as well as parameterizations explicitly including air–sea fluxes due to bubble formation (Woolf, 1997; Stanley et al 2009; Liang et al., 2013). In the second set of sensitivity tests, we prescribed the atmospheric mole fraction by adding and subtracting an associated uncertainty of 3 ppm due to spatial variabilities (Keraghel et al., 2020). Finally, in a third set of sensitivity tests, we performed two simulations by adding simple estimates of the calcium carbonate production into the equations of DIC (Eq. 1) and excess negative charge. Following the study of Palevsky and Quay (2017), we first estimated it based on the PIC:POC ratio and the NCP. Miquel et al. (2011) estimated that the PIC:POC ratio at 200 m depth varied between 0.31 and 0.78, with a mean value of 0.5, based on sediment trap measurements at the EMSO-DYFAMED site. Besides, Kessouri et al. (2018) estimated that POC export represents ∼70 % of the total OC (TOC) export (the remaining 30 % being attributed to DOC export). Thus, by assuming the ratio of calcium carbonate production to NCP is close to the PIC:TOC ratio, we added into Eq. (1) a consumption term representing 36 % of NCP for the mean value of the PIC:POC ratio and 22% and 55% for the minimum and maximum ratio values, respectively. This term, multiplied by 2, was added into the equation of the rate of change of the excess negative charge (Middelburg, 2019). In a second sub-test, we added a CaCO₃ production term based on the parametrization used in the Gulf of Lion's shelf modeling study by Lajaunie-Salla et al. (2021) (their Table A4, Precip = $k_{\text{precip}}\frac{({\mathrm{\Omega }}_c-\mathrm{1})}{\mathrm{0.4}+({\mathrm{\Omega }}_c-\mathrm{1})}\sum _{i=\mathrm{1}}^{\mathrm{3}}({\text{GPP}}_i-{\text{RespPhy}}_i)$, where k_precip is the PIC:POC ratio, and Ω_c is the aragonite saturation, set at 3.5 based on Schneider et al. (2007)).

2.2 Data used for model skill assessment

2.2.1 BOUSSOLE buoy and EMSO-DYFAMED mooring site observations

To assess the time evolution of sea surface properties, we used high-frequency temperature, salinity, and pCO₂ data collected at 3 m depth at the BOUSSOLE mooring site (43^∘22^′ N, 7^∘54^′ E; depth: ∼2400 m; green star in Fig. 1), in the Ligurian Sea in 2013 (Antoine et al., 2006; Merlivat et al., 2018). Temperature and salinity were measured using a Sea-Bird SBE 37-SMP MicroCAT instrument. The sensors were cross-calibrated before and after each mooring deployment with the ship CTD by performing a high-temporal-resolution sampling cast with 30 min long time series at the fixed depths of 300 and 1000 m. This allows for high accuracy of 0.001 ^∘C in temperature and 0.005 in salinity (Houpert, 2013). fCO₂ measurements were monitored using a CARIOCA sensor whose accuracy is estimated at 2 µatm. A detailed description of these data is given in Merlivat et al. (2018).

We also used monthly vertical profiles of temperature, salinity, dissolved oxygen, dissolved inorganic carbon, and total alkalinity collected from September 2012 to September 2013 at the EMSO-DYFAMED site (43^∘25^′ N; 7^∘52^′ E; depth: 2350 m; black star in Fig. 1) (Coppola et al., 2020, 2023), located 5 km from the BOUSSOLE site. Note that temperature and salinity were collected using a Sea-Bird SBE911. Dissolved oxygen measurements were performed using Winkler titration at each CTD cast and were used to correct the SBE43 sensor data by adjusting the calibration coefficients (Coppola et al., 2018). DIC and total alkalinity were measured via potentiometric titration following the methods described by Edmond (1970) and DOE (1994) with an accuracy estimated to be between 1.5 and 3 µmol kg⁻¹. They were analyzed by the SNAPO-CO2 national service (Service National d'Analyse des Paramètres Océaniques du CO₂). pCO₂ and pH_T (pH at total scale) were deduced from total alkalinity and total inorganic carbon using the carbonate system CO2SYS program (Lewis et Wallace, 1998; van Heuven et al., 2011), as in the system carbonate module described in Sect. 2.1.1.

2.2.2 DEWEX and MOOSE-GE cruise observations

To assess the horizontal and vertical distribution of the simulated DIC concentration, we used in situ observations collected during two cruises carried out in the framework of the DEWEX project on board the RV Le Suroît. The first cruise, DEWEX Leg1, was carried out in February 2013 during the active phase of deep convection (Testor, 2013), and the second one, DEWEX Leg2, was carried out in April 2013 during the following spring bloom (Conan, 2013). In addition, we used observations from the 2013 MOOSE-GE cruise, conducted during the stratified, oligotrophic season in June–July 2013, on board RV Tethys II (Testor et al., 2013). During the three cruises, the total dissolved inorganic carbon measurements (DEWEX Leg1: 19 stations, DEWEX Leg2: 14 stations, MOOSE-GE: 20 stations) were collected into acid-washed 500 cm³ borosilicate glass bottles and poisoned with 100 mm³ of HgCl₂ following the recommendations of DOE (1994) and Dickson et al. (2007). Samples were stored in the dark at 4 ^∘C pending analysis. Analyses were also performed by the national service SNAPO-CO2. Following Wimart-Rousseau et al. (2021), the values of total dissolved inorganic carbon and total alkalinity below 500 m depth, outside the range defined by ±2 the standard deviation in relation to the mean value, were not considered. The accuracy of the measurements was estimated at 1.5–3 µmol kg⁻¹.

2.3 CANYON-MED neural networks

Since in situ observations of the carbonate system remain scarce, the comparison of model outputs versus in situ observations is completed with a comparison with dissolved inorganic carbon, total alkalinity, and pH_T derived from the CANYON-MED neural networks developed by Fourrier et al. (2020) for the Mediterranean Sea. pCO₂ was derived from CANYON-MED outputs using the carbonate system CO2SYS program. The neural networks were applied at the EMSO-DYFAMED and BOUSSOLE sites using as input parameters pressure, temperature, salinity, and dissolved oxygen measured between 3 and 14 m depth, as well as the geolocation and date of the sampling. The accuracies of the derived pH_T, dissolved inorganic carbon, and total alkalinity were estimated at 0.014, 12 and 13 µmol kg⁻¹, respectively, for the entire Mediterranean Sea area. However, the accuracy of the neural networks was greatly improved locally, as Fourrier et al. (2022) showed specifically for the Gulf of Lion and Ligurian Sea over the period 2013–2020.

Figure 3 Time series of (a) temperature, (b) salinity, (c) pCO₂, (d) DIC, (e) total alkalinity, and (f) pH at total scale, modeled at 3 m depth (line in black), observed (small red dots at BOUSSOLE site and pink points at EMSO-DYFAMED site between 3 and 14 m depth) and computed with CANYON-MED neural networks (small blue dots at BOUSSOLE at 3 m and blue squares at EMSO-DYFAMED site between 3 and 14 m depth; error bars are indicated in gray). Correlation coefficient, RMSE, and bias between model outputs and BOUSSOLE observations are indicated in (a), (b), and (c).

3 Assessment of the model skills

3.1 Comparison at the BOUSSOLE and EMSO-DYFAMED sites

Figure 3 shows the seasonal cycle of temperature, salinity, pCO₂, DIC, total alkalinity, and pH_T observed and modeled at the surface at the EMSO-DYFAMED and BOUSSOLE sites. The temperature is very well simulated, with a highly significant correlation of 0.997 (p value < 0.01), an RMSE (root mean square error) of 0.50 ^∘C, and a bias of −0.31 ^∘C compared to BOUSSOLE buoy observations (Fig. 3a). Regarding the salinity, the model is close to the observations, except from December to February, when it underestimates the BOUSSOLE observations (Fig. 3b). The correlation (R=0.59, p value < 0.01) remains significant compared to the observations at the buoy, where the model has a bias of −0.04 and an RMSE of 0.10 over the whole study period. The modeled pCO₂ is in good agreement with observations and values derived with CANYON-MED neural networks, with a significant correlation of 0.90 (p value < 0.01), a bias of 5.74 µatm, and an RMSE of 25.57 µatm, compared to the BOUSSOLE observations (Fig. 3c). The model simulates low values in winter, when temperatures were minimum, and in spring, during the phytoplankton bloom identified by Kessouri et al. (2018). The maximum values are modeled in summer due to warming, as in observations and in CANYON-MED results. The seasonal variation of modeled DIC is in agreement with that observed and deduced from CANYON-MED (Fig. 3d), showing an increase in winter until the end of the deep-mixing period and a drop in spring when the growth of phytoplankton was maximum. The seasonal dynamics of modeled alkalinity shows minimum values in November and December, an increase in winter, and low variations in spring (Fig. 3e). The increase in winter is also found in observations and CANYON-MED results. In summer, the model underestimates both datasets by ∼10–15 µmol kg⁻¹. The pH_T seasonal variation in observations, simulation, and CANYON-MED results all indicate a drop in summer, mostly due to the increase in temperature (Hagens and Middleburg, 2016). Finally, the model results for the variables of the carbonate system are also consistent with the seasonal variability derived by Coppola et al. (2020) from monthly mean DYFAMED observations over the period 1998–2016. The modeled variables fall within the range of the observed values gathered in this synthesis study: 300–570 µatm for pCO₂, 2200–2340 µmol kg⁻¹ for DIC, 2510–2600 µmol kg⁻¹ for alkalinity, and 7.9–8.2 for pH_T.

3.2 Comparison with DEWEX and MOOSE-GE cruise observations

Previous studies based on the present coupled model concluded that the model shows good performance in reproducing fall and winter mixing (Estournel et al., 2016), as well as the timing and intensity of the phytoplankton blooms (Kessouri et al., 2018) and the seasonal dynamics of the dissolved oxygen (Ulses et al., 2021) and inorganic nutrients (Kessouri et al., 2017) over the three cruise periods. Here, we focus the assessment of the model skills on the seasonal dynamics and spatial variability of the carbonate system, especially the DIC concentration.

Figure 4 Surface dissolved inorganic carbon (DIC) concentration (µmol kg⁻¹) observed (a, c, e) and modeled (b, d, f) over the (a, b) DEWEX Leg1 (1–21 February 2013), (c, d) DEWEX Leg2 (5–24 April 2013), and (e, f) MOOSE-GE (11 June–9 July 2013) cruise periods. The correlation coefficient (R), root mean square error (RMSE), and bias between surface observed and modeled DIC are indicated in (b, d, f).

A comparison of modeled surface (from 5 to 10 m depth) DIC concentrations from DEWEX Leg1, DEWEX Leg2, and MOOSE-GE cruise observations is shown in Fig. 4. Figure 5 shows the modeled and observed DIC vertical profiles in the deep-convection area (indicated in Fig. 1 and defined in Sect. 2.1.3) and south of this zone (latitude <41^∘ N) in the Balearic Front, where winter vertical mixing is shallower (Ulses et al., 2021), during the three cruise periods. Comparisons were performed by extracting model outputs at the same date and location as measurements. The statistical analysis for surface DIC concentrations indicates significant spatial correlations of 0.78, 0.67, and 0.54 (p value < 0.01); RMSEs of 18.71, 24.25, and 12.04 µmol kg⁻¹; and biases of 5.25, 15.27, and 1.94 µmol kg⁻¹ , compared, respectively, to DEWEX Leg1, DEWEX Leg2, and MOOSE-GE observations. The model correctly represents observed spatial variability. In winter, during the intense vertical mixing period, maximum values near the sea surface are found in the deep-convection zone (Fig. 4a–b), where the vertical profiles are almost homogeneous over the whole water column (Fig. 5a). In spring, the model represents the drops observed in the surface layer in both zones (Figs. 4c–d and 5b) due to phytoplankton growth (Kessouri et al., 2018). Finally, the model reproduces the low values observed at the surface in the deep-convection zone during the stratified period where the vertical profiles approach those observed in the southern zone (Figs. 4e–f and 5c). The statistical analysis based on the whole vertical profiles shows that the model is significantly correlated with the observations (R>0.7, p value < 0.01), has an RMSE smaller than 20 µmol kg⁻¹ and a standard deviation smaller than 25 µmol kg⁻¹, close to that of the observations (bottom panel in Fig. 5a–c).

Figure 5 Comparison between observed and modeled dissolved inorganic carbon (DIC) in the northwestern Mediterranean Sea over the (a) DEWEX-Leg1 (10–12 February 2013), (b) DEWEX-Leg2 (8–10 April 2013), and (c) MOOSE-GE (27 June–5 July 2013) cruise periods. (a–c) Observed (blue and red – mean is shown with dashed lines, and shaded areas are for standard deviation) and modeled (green and orange – mean is shown with solid lines, and shaded areas are for standard deviation) profiles in the deep-convection area and south of it (latitude < 41^∘ N); (d–f)  Taylor diagram summarizing the statistical comparisons between the entirety of the observations (noted O) collected during the three cruises and the corresponding model outputs (noted M): radius is standard deviation, angle is correlation coefficient, and distance from the origin is root mean square error (RMSE).

4 Results

4.1 Seasonal cycle of dissolved inorganic carbon

We analyze here the seasonal cycle of the modeled dissolved inorganic carbon over the period September 2012–September 2013. Figure 6 shows the time evolution of atmospheric and hydrodynamic conditions, as well as of surface pCO₂ and DIC fluxes, while Fig. 7 displays the cumulative DIC fluxes and the resulting change in DIC inventory for the upper (surface–150 m) and deeper (150 m–bottom) layers since 1 September 2012, along with schemes of seasonal budgets. Figure 8 shows maps of the seasonal mean pCO₂ difference between the atmosphere and surface seawater (pCO_2,atm–pCO_2,sea) and air-to-sea CO₂ flux. Finally, the time evolution of the DIC concentration profile averaged over the deep-convection area is shown in Fig. 9. The study year was divided into seasonal periods defined according to the timing of stratification and biogeochemical processes, specific to this year, according to the studies of Kessouri et al. (2017, 2018).

Figure 6 Time series of modeled (a) total surface heat fluxes (W m⁻²), (b) mixed-layer depth (m), (c) sea surface temperature (^∘C), (d) sea surface and atmospheric pCO₂ (µatm), (e) wind speed (m s⁻¹), (f) air-to-sea CO₂ flux (mmol C m⁻² d⁻¹), (g) DIC total (vertical plus lateral) transport in the upper (light blue) and deeper layer (dark blue) towards the deep-convection area (mmol C m⁻² d⁻¹), and (h) DIC biogeochemical production (see Eq. 1, mmol C m⁻² d⁻¹) in the upper (green) and deeper (brown) layer. All the parameters are spatially averaged over the defined deep-convection area (spatial mean is shown with the solid line, and the shaded area is for SD). Sources: ECMWF for air–sea heat flux and wind speed, SYMPHONIE/Eco3M-S for the other parameters and fluxes. The blue shaded area corresponds to the deep-convection period (period when spatially averaged mixed-layer depth >100 m). Note that the range of the y axis varies for the different carbon fluxes, and due to higher values, SD for transport is not shown.

4.1.1 Autumn (1 September–27 November, 88 d)

After a period of alternating heat gain and loss events, and with the intense northerly wind event that occurred at the end of October (spatial mean heat loss reached 1000 W m⁻²), the deep-convection area continuously transferred heat to the atmosphere (Fig. 6a and e). The sea surface heat loss induced drops of surface temperature (Fig. 6c) and vertical mixing with a mixed layer that, on average, was shallower than 50 m during the whole autumn period (Fig. 6b).

Regarding biogeochemical processes, the shallowing of the nutricline and vertical mixing events induced nutrient supplies into the upper layer that favored primary production and growth of phytoplankton near the surface from the end of October (Kessouri et al., 2017, 2018, and Fig. 9a). This led, notably, temporally, to a low DIC concentration in the mixed layer at the end of October and end of November (Fig. 9b). However, over the whole fall period, respiration dominated primary production (Kessouri et al., 2018) (Fig. 6h), yielding a cumulative flux of DIC of 0.56 mol C m⁻² and 0.30 mol C m⁻² in the upper (surface–150 m) and deeper (150 m–bottom) layers, respectively, over the 88 d period (Fig. 7).

Figure 7 Time series of variation in dissolved inorganic carbon (DIC) inventory since 1 September 2012 (black) and cumulative air–sea (red), physical (light and dark blue), and biogeochemical (bright and brown green) flux of dissolved inorganic carbon in the (a) upper (surface to 150 m) and (b) deeper (150 m to bottom) layers from September 2012 to September. The cumulative flux at a day d is the time-integrated flux over the period from the 1 September 2012 to day d (mol C m⁻²). Positive values of fluxes represent DIC inputs for the deep-convection area. The blue shaded area corresponds to the deep-convection period (period when spatially averaged mixed-layer depth >100 m). The DIC inventory on 1 September 2012 was 353 and 5560 mol C m⁻² in the upper and deeper layers, respectively. (c) Scheme of cumulative seasonal fluxes (in mol C m⁻²) over the respective periods (fall: 88 d, winter: 116 d, spring: 74 d, and summer: 87 d). Resp. stands for respiration, and GPP stands for gross primary production. The direction of the arrows indicates the direction of the fluxes, and positive values correspond to DIC inputs for the deep-convection area.

The net physical fluxes of DIC at the deep-convection area boundaries fluctuated between −150 and 150 mmol C m⁻² d⁻¹ in the upper and deeper layers (Fig. 6g). More specifically, our model results show a net cumulative upward transport of DIC of 35.08 mol C m⁻² into the upper layer of the area, while the cumulative lateral transport led to a loss of DIC of 34.59 mol C m⁻² in the upper layer and a gain of DIC of 33.70 mol C m⁻² in the deeper layer of the deep-convection area (Fig. 7c). This is probably induced by the upwelling and the surface divergence associated with the dynamics of the cyclonic gyre (Estournel et al., 2016). The lateral and vertical physical transfers at the boundaries resulted in a net increase in DIC inventory in the upper layer of 0.49 mol C m⁻² and a net decrease in DIC inventory in the deeper layer of −1.38 mol C m⁻².

Over the whole autumn period, sea surface pCO₂ decreased with temperature (Fig. 6c and d, coefficient correlation of 0.99 with p value < 0.01). The air–sea CO₂ flux displayed strong out-gassing peaks exceeding −20 mmol C m⁻² d⁻¹ during the northerly wind events that occurred in early and mid-September (Fig. 6f), when the pCO₂ difference was greater than 60 µatm (Fig. 6d). During the intense event of heat loss and cooling at the end of October, the deep-convection area became deficient compared to the atmosphere (sea surface pCO₂ smaller than atmospheric pCO₂) (Fig. 6d) and started to absorb atmospheric CO₂ (Fig. 6f). The air–sea flux displayed ingassing peaks smaller than 6 mmol C m⁻² d⁻¹ in November (Fig. 6f), characterized by pCO₂ differences smaller than 30 µatm (Fig. 6d) and moderate wind speeds (Fig. 6e). Considering the whole autumn period (88 d), the deep-convection area was a weak source of CO₂ for the atmosphere, with a cumulative air–sea flux that amounted to −0.19 mol C m⁻² (Fig. 7a, c). Characterized by low temperature, the deep-convection area, outside its southwestern region, showed small averaged pCO₂ differences (<25 µatm, Fig. 8a) and weak mean out-gassing fluxes compared to the surrounding open sea (Fig. 8e). In the southwestern region, intrusions of warm water from the Balearic Sea, characterized by higher pCO₂ and associated with strong wind speeds (not shown), favored maximum out-gassing fluxes.

Figure 8 Seasonal averages of modeled (a–d) pCO₂ difference (pCO_2,atm–pCO_2,sea, in µatm) and (e–h) air-to-sea CO₂ flux (mmol C m⁻² d⁻¹). Note that the periods of seasons here are defined in Sect. 4.1 according to mixed-layer depth and biogeochemical processes (fall: 1 September–27 November, winter: 28 November–23 March, spring: 24 March–5 June, summer: 6 June–31 August). Gray lines indicate in (a) to (d) pCO₂ difference isolines (−50, −25, −10, 0, 10, 25, 50 µatm) and in (e) to (h) CO₂ flux isolines (−10, −5, 0, 5, 10 mmol C m⁻² d⁻¹); the black line indicates the limit of the deep-convection area.

Finally, the inventory of DIC changed by 0.85 and −1.08 mol C m⁻² over the autumn period in the upper and deeper layers, respectively (Fig. 7). To sum up, the upper layer gained DIC through biogeochemical and physical processes and lost DIC through out-gassing to the atmosphere.

Figure 9 Time evolution of (a) Chlorophyll-a (mg Chl m⁻³) and (b) dissolved inorganic carbon (DIC, µmol C kg⁻¹) concentration profile, with mixed-layer depth (m) indicated by the black line, horizontally averaged over the deep-convection area.

4.1.2 Winter (28 November–23 March, 116 d)

The winter period can be further divided into two sub-periods based on the magnitude of vertical mixing (Kessouri et al., 2017).

• Winter sub-period 1. During the first winter period (end of November–mid-January), heat loss events induced an intensification of vertical mixing that remained moderate with the spatially averaged mixed-layer depth above the euphotic layer depth (150 m) (Figs. 6a, b and 9). Vertical mixing induced new supplies of inorganic nutrients into the upper layer, supporting primary production near the surface (Kessouri et al., 2017, 2018). From mid-December, a net consumption of DIC is modeled in the whole upper layer (Fig. 6h). The cumulative biogeochemical fluxes in the upper layer showed a progressive decrease during the first winter period, leading to a weak cumulative consumption of DIC of 0.05 mol C m⁻² over 49 d (Fig. 7a).

  The physical fluxes at the limit of the upper layer of the deep-convection area showed similar patterns to those during autumn, with a cumulative upward flux of DIC into the upper layer of 41.40 mol C m⁻² over a 2.5 month period, almost counterbalanced by a cumulative lateral export of DIC of 40.44 mol C m⁻² in the upper layer and a cumulative lateral input of DIC of 39.90 mol C m⁻² in the deeper layer. The net physical fluxes led to a gain of 0.97 mol C m⁻² in the DIC inventory of the upper layer and a loss of DIC of 1.51 mol C m⁻² in the deeper layer (Fig. 7).

  Until mid-January, sea surface pCO₂ continued to decrease with temperature, yielding a reinforcement of the pCO₂ difference from 30 to 40 µatm (Fig. 6c and d). The spatial mean air–sea CO₂ flux reached 15 mmol C m⁻² d⁻¹ during the northerly wind events (Fig. 6e and f). Over the first winter period (49 d) the deep-convection area absorbed 0.35 mol m⁻² of atmospheric CO₂ (Fig. 7a).

  In summary, the euphotic layer showed an increase in the DIC inventory of 1.27 mol C m⁻², resulting from a gain through air–sea fluxes and physical transport, and a weak net biogeochemical flux (Fig. 7a).

• Winter sub-period 2. The second winter period, which began in mid-January, was characterized by deep convection. The mixed layer deepened strongly during the intense heat loss events that occurred until the end of February (Fig. 6a–b). After a 16 d pause during which surface restratification caused the mixed layer to be shallower, a new northerly wind generated a secondary deep-convection event in late March. During the second winter period (mid-January–end of March), vertical mixing reached deep water masses (Fig. 6b). Surface temperature remained relatively constant at a value of ∼12.9 ^∘C, close to deep water temperature (Fig. 6c).

  Regarding the biogeochemical fluxes, the net consumption of DIC progressively increased in the upper layer with the decrease in heterotrophic respiration and the moderate increase in primary production rates (Kessouri et al., 2018). It accelerated when vertical mixing ceased in mid-March and remained high until the end of the period. The cumulative biogeochemical flux reached −1.49 mol C m⁻² over this sub-period of 67 d.

  Vertical and horizontal exchanges showed the same patterns as during the preconditioning period (fall and first period of winter), but the difference between vertical and horizontal transports became more pronounced, leading to net physical fluxes exceeding 110 mmol C m⁻² d⁻¹ in the upper layer and −160 mmol C m⁻² d⁻¹ in the deeper layer during the four northerly wind events. The physical fluxes integrated over the second winter period reached 3.48 and −5.36 mol C m⁻² in the upper and deeper layers, respectively.

  Despite the biogeochemical consumption of DIC, a progressive increase in DIC concentration in the upper layer is clearly visible in Fig. 9b due to vertical transport. The upward fluxes led to an increase in sea surface pCO₂, showing values that were, on average, close to equilibrium (Fig. 6d). The pCO₂ difference decreased; despite intense wind events, air–sea flux peaks remained lower than 12 mmol C m⁻² d⁻¹; and, finally, the cumulative air–sea flux reached 0.28 mol C m⁻² over the second winter sub-period of 67 d (a lower value and flux (3.1 versus 7.3 mmol C m⁻² d⁻¹) than over the first winter period).

  To summarize, the upper layer showed a gain in DIC inventory through vertical transport and, to a lesser extent, uptake of atmospheric CO₂, while it lost DIC through lateral transport and biogeochemical processes (Fig. 7a).

• All winter periods. The pCO₂ difference and air–sea fluxes averaged over the whole winter period (end November–end March; Fig. 8b and f, respectively) integrate various processes: (1) mistral and tramontane northerly winds blowing, on average, over a northwest–southeast axis over the Gulf of Lion (Fig. S2e, f) intensified air–sea fluxes on this axis; (2) low sea surface temperature in the deep-convection region favored an amplification of the pCO₂ difference and maximum air–sea fluxes during the first winter period, especially at the northern edge of the convection zone (Fig. S2a, c); and (3) high surface DIC concentrations in the regions of intense vertical mixing generated a reduction in the pCO₂ difference and air–sea fluxes, especially during the second winter period (Fig. S2b, d). The maxima of the lateral DIC transport in the upper layer of the water column averaged over the whole winter period are found in the general circulation, especially in the Northern Current, the Balearic Current, and the Balearic Front, separating the southern, less salty Atlantic waters from the deep-convection salty waters (Fig. 10a, c). The instabilities developing at the periphery of the deep-convection area favored the incorporation of saltier and DIC-rich waters in the general circulation through a bleeding effect, similarly to that described by Herrmann et al. (2008) for the export of newly formed dense waters from the deep-convection area (i) at the western boundaries of the deep-convection area towards the Balearic Sea and towards the Algerian basin through the southern extension of the Balearic Current and (ii) along the Balearic Front between the Minorca Balearic Island and Corsica, as illustrated in Fig. 10c, d. Finally, Fig. 10b shows that the vertical DIC supply into the upper layer during winter resulted from small-scale upward and downward vertical fluxes due to the absence of stratification.

Figure 10 (a) Winter horizontal flux of dissolved inorganic carbon (DIC, mol C m⁻² s⁻¹) vertically integrated over the upper layer (0–150 m), (b) winter vertical DIC flux (mol C m⁻² s⁻¹) at 150 m, (c) surface salinity, and (d) DIC concentration (µmol C kg⁻¹) on 4 March 2013. The red and black lines in (a) and (b), respectively, indicate the limit of the deep-convection area.

4.1.3 Spring (24 March–5 June, 74 d)

From the end of March, total surface heat flux remained mostly positive (Fig. 6a). The water column rapidly stratified with a mixed-layer thickness lower than 50 m (Fig. 6b). The cessation of deep vertical mixing favored the onset of a spring phytoplankton bloom, with a peak of primary production and phytoplankton concentration at the surface mid-April (Fig. 9a; Kessouri et al., 2018), leading to a sustained consumption of DIC (Figs. 6h, 7a, and 9b). From mid-April, the near-surface layer became depleted in nutrients, a deep chlorophyll maximum (DCM) formed and progressively deepened (Fig. 9a; Kessouri et al., 2018), and the consumption of DIC was slowed down in the upper layer (Figs. 6h and 7a). In the deeper layer, the rate of DIC production through heterotrophic remineralization of organic carbon, exported during the deep convection and the bloom, was maximum during this period (Fig. 6h). The contribution of biogeochemical processes to the DIC inventory over the spring period resulted in a loss of 2.19 mol C m⁻² in the upper layer and a gain of DIC of 0.91 mol C m⁻² in the deeper layer over the 74 d period (Fig. 7).

The net exchanges of DIC at the deep-convection area boundaries weakened in both the upper and deeper layers compared to the previous period (Fig. 6g). Over this restratification period, characterized by baroclinic instabilities (Jones and Marshall, 1997), the cumulative physical exchange flux was negative, leading to a loss of DIC of 0.80 mol C m⁻² in the upper layer and of 1.62 mol C m⁻² in the deeper layer over the 74 d period (Fig. 7).

Sea surface pCO₂ decreased in early spring (Fig. 6d), when DIC was consumed through strong primary production. The pCO₂ difference reached a maximum positive value of 75 µatm in mid-April, at the peak of the bloom. Afterwards, when the DCM was formed (Fig. 9a), sea surface pCO₂ varied with sea surface temperature, which increased until early May and remained around 15 ^∘C until the end of the spring period, when the pCO₂ difference was around 30 µatm (Fig. 6c and d). The air–sea flux showed positive values, with a maximum uptake of 30 mmol C m⁻² d⁻¹ in mid-April during the bloom. Over the whole spring period, the deep-convection area absorbed 0.45 mol C m⁻² of atmospheric CO₂ (Fig. 7a and c). The mean spring pCO₂ difference varied between 30 and 40 µatm over the deep-convection area (Fig. 8c). The mean air–sea CO₂ flux was the strongest in the Gulf of Lion (Fig. 8g), where wind speed was maximum (not shown).

To sum up, the upper layer gained DIC through air–sea flux and lost DIC through biogeochemical processes and physical fluxes (Fig. 7). The loss of DIC by physical fluxes resulted from a prevalence of the loss through lateral transport in comparison to the gain through vertical transport.

4.1.4 Summer (6 June–31 August, 87 d)

Two long episodes of heat gain in June and July (Fig. 6a) generated a strong stratification of the water column with a mixed layer shallower than 20 m (Figs. 6b and 9), along with increases in surface temperature (Fig. 6c). Early August surface temperature started to slowly decrease. From August onwards, the DIC drawdown due to biogeochemical processes decreased, the primary production rate becoming close to the respiration rate (Fig. 6h). We estimate the contribution of biogeochemical processes in summer to be −0.57 mol C m⁻² in the surface layer and 0.39 mol C m⁻² in the deeper layer over an 87 d period (Fig. 7). The cumulative physical fluxes were again negative in both layers (Fig. 7). On average, the lateral export of DIC from the upper layer prevailed over the vertical supply of DIC into this layer.

Sea surface pCO₂ shows a similar evolution to that of temperature (correlation coefficient of 0.99 (p value < 0.01), Fig. 6c and d). Very quickly at the beginning of the summer period, the deep-convection area came to be in excess of CO₂ relative to the atmosphere. The pCO₂ difference reached a maximum absolute value of 135 µatm at the end of July. The model outputs show peaks of CO₂ out-gassing varying between 18 and 28 mmol C m⁻² d⁻¹ during the northerly wind events that occurred from the end of July. The deep-convection area released 0.42 mol C m⁻² over the whole summer period (Fig. 7a, c). Similarly to the fall situation, the out-gassing fluxes were maximum in the western part of the delimited deep-convection area (Fig. 8h), influenced by both strong wind speeds and the arrival of warm waters from the Balearic Sea through anticyclonic circulations, which were characterized by higher pCO₂ differences (between 70 and 90 µatm, Fig. 8d).

In summary, in summer, the upper layer of the deep-convection area decreased its DIC inventory in response to physical and biogeochemical processes, as well as out-gassing towards the atmosphere (Fig. 7a, c).

Figure 11 Scheme of the annual carbon budget for the period September 2012 to September 2013 from the coupled model SYMPHONIE-Eco3M-S. Fluxes are indicated in mol C m⁻² yr⁻¹. The direction of the arrows indicates the direction of the fluxes, and positive values of fluxes represent DIC inputs for the deep-convection area (positive vertical fluxes represent inputs for the upper layer).

4.2 Annual carbon budget

Figure 11 shows a schematic of the annual budget of dissolved inorganic carbon in the deep-convection zone. Our model results show that the deep-convection area acted as a moderate CO₂ sink for the atmosphere on an annual scale over the period September 2012–September 2013. We estimate that it absorbed 0.5 mol C m⁻² yr⁻¹ of atmospheric CO₂. This uptake of atmospheric CO₂ displayed spatial variability (Fig. 12). It was greater than 1 mol C m⁻² yr⁻¹ in the northern edge of the area along the Northern Current flowing over the Gulf of Lion continental slope and became less than 0.25 mol C m⁻² yr⁻¹ in the western and eastern edge areas. One can notice that the annual rate remained lower than on the Gulf of Lion's shelf, which is beyond the scope of this study. Within the sea, biogeochemical processes induced an annual DIC consumption of 3.7 mol C m⁻² yr⁻¹ in the upper layer and a DIC gain of 2.3 mol C m⁻² yr⁻¹ in the deeper layers. The deep-convection area thus appears as a net autotrophic region from a biogeochemical point of view, with a DIC consumption of 1.5 mol C m⁻² yr⁻¹ considering the whole water column.

Our estimate of net physical fluxes (lateral plus vertical) is an input of 3.3 mol C m⁻² yr⁻¹ into the upper layer and an export of 11.0 mol C m⁻² yr⁻¹ from the deeper layer. Specifically, the model indicates a vertical DIC supply of 133.2 mol C m⁻² yr⁻¹ from the deeper layer to the upper layer, partly offset by a lateral outflow of 129.8 mol C m⁻² yr⁻¹ in the upper layer and an inflow of 122.2 mol C m⁻² yr⁻¹ in the deeper layer. The budget in the deep layer masks different signs of physical fluxes: if the deeper layer is subdivided into an intermediate layer (150–800 m) and the deeper most layer (800 m–bottom), we find that the former, the intermediate layer, gained 83.1 mol C m⁻² yr⁻¹ of DIC through vertical transport, while it lost 87.6 mol C m⁻² yr⁻¹ of DIC through lateral export. Finally, our model shows that the convection zone was a source of DIC of 8.7 mol C m⁻² yr⁻¹ for the rest of the western Mediterranean Sea. While the DIC inventory in the upper layer remained stable (decrease of 0.07 mol C m⁻² yr⁻¹), the DIC inventory in the deeper layer experienced a decrease of 8.7 mol C m⁻² yr⁻¹. This loss occurred mainly during deep convection and, to a lesser extent, during the preconditioning period (in autumn and early winter).

Finally, we complete the inorganic carbon budget with the labile organic carbon fluxes (refractory organic carbon is not considered in our model). We estimate that, during the studied period a lateral export of organic carbon of 1.1 and 0.3 mol C m⁻² yr⁻¹, took place in the upper and deeper layers, respectively. The modeled downward export of organic carbon amounted to 2.3 mol C m⁻² yr⁻¹.

5 Discussion

Based on high-resolution 3D modeling, we have estimated a DIC budget in the northwestern Mediterranean deep-convection zone over an annual period, September 2012–September 2013. Our results show that biogeochemical and physical processes, through their impacts on DIC concentration, both play a major role in the intensity and sign of the air–sea exchanges in the deep-convection area.

5.1 Assessment of the seasonal cycle of the pCO₂

The seasonal pattern of the simulated sea surface pCO₂ averaged over the deep-convection area, in deficit compared to the atmosphere from November 2012 to early June 2013 and in excess relative to the atmosphere the rest of the year, is similar to the one simulated at the EMSO-DYFAMED and BOUSSOLE sites (Fig. 3) and is in good agreement with those described in previous observational and modeling studies at these sites (Hood and Merlivat, 2001; Copin-Montégut and Bégovic, 2002; Bégovic and Copin-Montégut, 2002; Mémery et al., 2002; Copin-Montégut et al., 2004; Merlivat et al., 2018; Coppola et al., 2020). The high-frequency measurements at the CARIOCA buoy described by Hood and Merlivat (2001) and Merlivat et al. (2018) indicated an interannual variability of 4–5 weeks in terms of the date at which the pCO₂ difference changes sign, depending on the interannual variability of air–sea heat flux and of the bloom onset. They showed that, in autumn, the change of sign extends from mid-September to the end of October, and in spring, it extends from early May to mid June. The magnitude of the variation of the modeled sea surface pCO₂ (spatial mean: 193 µatm, between 322 and 515 µatm; at EMSO-DYFAMED site: 192 µatm) is in the range of those deduced from observations (120 µatm by Bégovic and Copin-Montégut (2002), 230 µatm by Hood and Merlivat (2001), and ∼200 µatm by Merlivat et al., 2018). More specifically, the impact of deep convection on sea surface pCO₂ through a large upward transport of CO₂-rich waters, leading to its increase, is consistent with previous studies. The measurements at the BOUSSOLE site gave evidence of brief windy periods marked by sea surface pCO₂ higher than atmospheric pCO₂ (Hood and Merlivat, 2001; Copin-Montégut and Bégovic, 2002; Copin-Montégut et al., 2004; Merlivat et al., 2018). The observations of the CASCADE cruise in March 2011 in the Gulf of Lion also showed high surface concentrations of DIC and sea surface pCO₂ higher than atmospheric pCO₂ in deep-convection cells (Touratier et al., 2016).

Figure 12 Modeled annual air-to-sea CO₂ flux (mol C m⁻² yr⁻¹) averaged over the period September 2012–September 2013. Gray lines indicate CO₂ flux isolines (0, 0.5, 1, 1.5 mol C m⁻² yr⁻¹), and the black line indicates the limit of the deep-convection area.

5.2 Estimate of the annual air–sea CO₂ flux and its uncertainties

On an annual scale, our results indicate that the deep-convection area was a sink for atmospheric CO₂. Our estimate of the annual air-to-sea flux is 0.47 mol C m⁻² yr⁻¹, which, considering the surface of the zone, corresponds to an uptake of atmospheric CO₂ of 0.4 Tg yr⁻¹. This estimate is associated with various sources of uncertainties related to the modeling of the different physical, biogeochemical, and air–sea exchange processes. Regarding the wind speed accuracy, Ulses et al. (2021) calculated a percentage bias of −0.5 % and a normalized RMSE of 13.9 % based on comparisons between ECMWF forcing fields and high-frequency measurements during the DEWEX cruises. The statistical analysis in Sect. 3 indicates that the model has low to moderate RMSE for surface temperature (0.50 ^∘C), salinity (0.10), pCO₂ (<26 µatm), and DIC (<24 µmol kg⁻¹) and low biases (respectively, −0.31 ^∘C, −0.04, 6 µatm, and <15 µmol kg⁻¹).

Figure 13 Sensitivity tests of the annual air-to-sea CO₂ flux estimate to the parameterization of gas transfer velocity, the variability of the mole fraction of CO₂ in the atmosphere, and the calcification processes. The black bar indicates the annual estimate in the reference simulation, and gray bars indicate the mean value for each of the three sets of sensitivity tests. For the sensitivity tests on the parametrization of gas transfer (from 2 to 9), relations with a quadratic (2), hybrid (3 to 5), or cubic (6) wind speed dependency are, respectively, in light pink, yellow, and orange, and relations that include explicit bubble parametrizations (7 to 9) are in red. For the test (14) of calcification processes, the bar indicates the result found for the mean PIC:POC ratio, while the black line indicates the range using the minimum and maximum PIC:POC ratios.

To assess the uncertainties linked to the calculation of the gas transfer coefficient, we performed sensitivity tests using eight other parameterizations of this coefficient (Sect. 2.1.4). The estimates of the annual air–sea flux using these parameterizations are displayed in Fig. 13 and Table S1. The results of these tests indicate that the deep-convection zone is found as a moderate CO₂ sink for the atmosphere using all parametrizations. The estimate in the reference run, based on the wind speed quadratic dependency relation established by Wanninkhof (1992), is close to the mean value of all the estimates at 0.43 (±0.12) mol C m⁻² yr⁻¹. The highest estimates were obtained using the relation from the cubic wind dependency parametrization of Wanninkhof and McGillis (1999) and the bubble-inclusive parametrizations of Woolf (1997) and Stanley et al. (2009). The lowest estimate, which is almost twice as small as the mean value, was obtained using the parametrization of Liss and Merlivat (1986). The second set of tests of sensitivity for atmospheric CO₂ forcing shows that the annual air–sea flux varies between 0.33 and 0.61 mol C m⁻² yr⁻¹ (SD of 0.20 mol C m⁻² yr⁻¹) if an uncertainty value of 3 ppm in relation to the atmospheric mole fraction is constantly subtracted and added, respectively. Finally, sensitivity tests taking into account supplementary consumption terms in the equations of DIC and of the excess of negative charges were performed to assess the potential influence of CaCO₃ precipitation on the air–sea CO₂ flux (Sect. 2.1.4). They show that not taking into account calcification processes could lead to an overestimation of the annual air–sea CO₂ uptake by 16 % to 57 % with estimates of 0.29 mol C m⁻² yr⁻¹ based on the mean PIC:POC ratio (Miquel et al., 2011) and the NCP (varying between 0.20 and 0.36 mol C m⁻² yr⁻¹ based on the measured maximum and minimum PIC:POC ratios (Miquel et al., 2011), respectively) and of 0.40 mol C m⁻² yr⁻¹ based on the parametrization used in Lajaunie-Salla et al. (2021). This demonstrates the need to better constrain this term in future studies on carbonate system dynamics.

5.3 Comparisons of air–sea CO₂ flux with previous studies in the Mediterranean Sea

Our estimates of the annual air–sea flux over the whole deep-convection area and at the DYFAMED site – 0.47 and 0.33 mol C m⁻² yr⁻¹, respectively – are close to those provided in previous observational and modeling studies at the DYFAMED site. Based on the parametrization of Liss and Merlivat (1986) and for the period 1995–1997, Hood and Merlivat (2001) found a value of 0.10–0.15 mol C m⁻² yr⁻¹ using hourly measurements, while Mémery et al. (2002) found a value of 0.15 ± 0.07 mol C m⁻² yr⁻¹ using a 1D model. We obtained an annual flux of 0.17 mol C m⁻² yr⁻¹ at the DYFAMED site using the same gas transfer relationship. The estimates of Bégovic (2001) and Copin-Montégut et al. (2004) varied between 0.42 and 0.68 mol C m⁻² yr⁻¹ for the period 1998–2000 using monthly CO₂ measurements and the parametrization proposed by Wanninkhof and McGillis (1999). Using the same parametrization, our estimate amounts to 0.40 mol C m⁻² yr⁻¹ at DYFAMED. Finally, Merlivat et al. (2018) estimated a close annual CO₂ air–sea flux of 0.45 mol C m⁻² yr⁻¹ using hourly measurements for the period 2013–2015.

Based on a 1D satellite data approach (Antoine and Morel, 1995) applied with a horizontal resolution of 0.5^∘ to the whole Mediterranean Sea over the period 1998–2004, D'Ortenzio et al. (2008) estimated an annual mean air–sea CO₂ flux ranging between 0 and 4 mol C m⁻² yr⁻¹ over the NW deep-convection area. The larger homogeneity in our estimates (varying between −0.1 and 1.2 mol C m⁻² yr⁻¹ inside the deep-convection area) could be partly ascribed to the horizontal diffusion and advection that were accounted for in our model. Using a 3D coupled physical–biogeochemical reanalysis of the Mediterranean Sea, von Schuckmann et al. (2018), over the period 1999–2016, and Cossarini et al. (2021), over the period 1999–2019, estimated a mean annual air–sea flux in the deep-convection zone ranging between 0 and 0.5 mol C m⁻² yr⁻¹, and 0 and 1 mol C m⁻² yr⁻¹, respectively. Our results in terms of spatial distribution, with minimum values in the western edge of the deep-convection zone and maximum values in the northern area of the Gulf of Lion, are also consistent with their results.

Our estimate is close to the annual flux estimated around 0.5 mol C m⁻² yr⁻¹ by Cossarini et al. (2021) in the South Adriatic Sea, another deep-convection area of the Mediterranean Sea.

Finally, it is also noteworthy that our estimate is found in the lower range of the annual flux estimated from experimental studies for the northern Adriatic and Aegean shelves, where dense-water formation also takes place and which are identified as sinks for atmospheric CO₂ most of the year and on an annual basis. With respect to the northern Adriatic shelf, our estimate is found to be close to the estimate of 0.4–0.5 mol C m⁻² yr⁻¹ for the period of 2014–2015 from Urbini et al. (2020) and between about 2- to 4-fold lower than the estimates of 0.8–0.9 mol C m⁻² yr⁻¹ from Urbini et al. (2020) over the period 2016–2017, of 1–1.1 mol C m⁻² yr⁻¹ from Catalano et al. (2014) and Cossarini et al. (2015), and of 2.2 mol C m⁻² yr⁻¹ from Cantoni et al. (2012) and Turk et al. (2010). Regarding the northern Aegean Sea, we found a lower winter flux than the one deduced from observations in February 2006 by Krasakopoulou et al. (2009) (8.6–14.7 mmol C m⁻² d⁻¹ versus 4.9 mmol C m⁻² d⁻¹ in our study). Our estimates are also lower than the CO₂ uptake exceeding 1 mol C m⁻² yr⁻¹ found for the northern shelves in the modeling studies of Cossarini et al. (2015, 2021). The higher fluxes over the continental shelves compared to our study area could be explained by a lower seawater temperature in winter, riverine nutrient inputs favoring intense primary production, and a transport of DIC associated with dense-water outflow towards the deep basin (Cantoni et al., 2016; Ingrosso et al., 2017).

5.4 The major influence of physical transport on the DIC budget of the deep-convection area

Our study confirms the importance of deep convection to the vertical DIC distribution and surface pCO₂ in the study area, as shown in previous observational studies (Copin-Montégut et al., 2004; Touratier et al., 2016), and highlights that physical transports play a crucial role in the DIC budget in this highly energetic region. Our 3D model results allowed us to distinguish the contribution of vertical and lateral transports to the net physical exchange flux. They both show a similar seasonal cycle, with greater magnitude (positive for the vertical transport and negative for the lateral transport with regard to the upper layer) in fall, the preconditioning phase, and in winter, the convection period, both being sea surface heat loss periods. During those periods, vertical supply overwhelmed lateral export in the upper layer. Conversely, during the stratification phase and stratified period, in spring and summer, lateral export prevailed over vertical supply. At the annual scale, we estimate that the vertical supplies amounted to 133 mol C m⁻² yr⁻¹. They were almost counterbalanced by a lateral transfer of 130 mol C m⁻² yr⁻¹ to adjacent upper-layer areas, which acted as major sinks of DIC for the deep-convection upper layer.

By estimating a water mass budget, we found that lateral and vertical DIC fluxes are highly significantly correlated with lateral (R=0.9998, p value < 0.01) and vertical (R=0.9998, p value < 0.01) water fluxes, respectively. Moreover, a detailed calculation of the water and DIC fluxes at the limits of the deep-convection area allowed us to evaluate the contribution of (1) the difference in inflowing and outflowing water fluxes at constant mean DIC concentration and (2) the difference in DIC concentrations at constant flux. We found that the first contribution is largely dominant compared to the second one, highlighting that the lateral transports of DIC are essentially related to the difference of inflows and outflows of water rather than to DIC concentration differences between deep-convection waters and surrounding waters. Strong mesoscale activities and instabilities within and on the edge of the mixed patch that characterized the convection zone shown in previous works (Marshall and Shott, 1999; Testor et al., 2018) could lead to this strong lateral transfer of water and associated DIC, as illustrated in Fig. 10. The findings of Waldman et al. (2018) showing water sinking in the general circulation suggest that DIC could be partly transferred back into deep waters in the boundary current. Studying the dissolved oxygen dynamics, Wolf et al. (2018) also found that lateral processes could play a major role in the biogeochemical annual budget in the deep convection located in the central Labrador Sea. We consider that the study of the contribution of the various lateral exchange mechanisms in the lateral DIC transfer, such as Ekman-driven transport, geostrophic advection, frontal processes, and submesoscale coherent vortices, is out of scope of this first work on the DIC budget, but further complementary works will be dedicated to this subject.

5.5 Net community production and air–sea flux relationships

Previous modeling studies (Herrmann et al., 2013; Ulses et al., 2016) showed that the northwestern Mediterranean deep-convection area acts as an autotrophic ecosystem with, on an annual timescale, gross primary production dominating over respiration and hence a positive NCP. The present modeling study displays an NCP in the upper layer of 3.7 mol C m⁻² yr⁻¹ and a DIC buildup of 2.3 mol C yr⁻¹ through respiration of heterotrophic organisms in the deeper layer. Our budget shows that, in the upper layer, the net biogeochemical loss term of DIC is counterbalanced by 88 % by physical gain fluxes and only by 13% by air–sea gain fluxes. It clearly appears that deep vertical mixing and advection significantly slowed down the atmospheric CO₂ uptake in winter by bringing into the upper layer remineralized organic carbon. We quantify here that the annual air–sea CO₂ flux is ∼8 times smaller than the annual NCP in the upper layer. These results are in line with the previous findings of Oschlies and Kähler (2004), Körtzinger et al. (2008a, b) and Palevsky and Nicholson (2018) and Palevsky and Quay (2017) regarding the reducing effect of winter ventilation on atmospheric CO₂ uptake in the northern Atlantic Ocean and in the Pacific Ocean, respectively.

In our simulation, the downward export of organic carbon (OC) at the base of the photic zone of the deep-convection area is estimated as 2.3 mol C m⁻² yr⁻¹. The results of Herrmann et al. (2013) and Ulses et al. (2016) showed, using a similar coupled modeling approach, that OC export is characterized by high interannual variability, with standard deviations between 24 % to 37 %, linked to the variability of the convection strength. The intensity of the winter vertical mixing has been shown to be subject itself to high interannual variability (Houpert et al., 2016; Margirier et al., 2020). Observations in the core of the convection zone by Margirier et al. (2020) evidenced that only intermediate convection events took place in the 4 years following the 2013 events. Thus, organic carbon transferred into the deeper layer could either have been stored in the deep-convection area until the next 2018 events (Fourrier et al., 2022), when it could have been reinjected in its remineralized form, or it could have been transferred, partly in remineralized form, towards the southwestern Mediterranean through the dispersion of the newly formed dense water (Schroeder et al., 2008; Beuvier et al., 2012). Here, we estimate that an amount of 0.3 mol C m⁻² yr⁻¹ of organic carbon was laterally exported from the deeper layer. We found that the vertical supply of DIC into the upper layer is 2 orders of magnitude higher than OC export or the upper layer and depth-integrated NCP. This is explained by the equilibrium role of the DIC lateral transfers towards the surrounding zone. This shows that a 1D approach would not be appropriate to take into account the complexity of the 3D mechanisms of exchanges for the DIC budget of this deep-convection area, which has to be integrated into a whole regional system, especially in the context of a changing atmosphere and ocean.

5.6 Contribution of the northwestern deep-convection region to the carbon budget of the Mediterranean Sea

Our results indicate that the NW Mediterranean deep-convection zone was a sink of carbon for the atmosphere and a source of carbon for the western Mediterranean Sea over the period of September 2012 to September 2013. More specifically, we found that the exchanges with the surrounding region were characterized by a net lateral input of total carbon into the deep layers of the deep-convection region, although organic carbon was exported towards the surrounding region, and a net lateral export of both organic and inorganic carbon in upper water masses (Fig. 11).

Previous studies investigating the air–sea CO₂ flux at the scale of the whole Mediterranean Sea showed that this sea acted as a moderate sink of atmospheric CO₂ over the past decades (Copin-Montégut, 1993; D'Ortenzio et al., 2008; Cossarini et al., 2021). According to those studies, the northern continental shelves and open seas (Gulf of Lion, Adriatic and Aegean seas) absorbed atmospheric CO₂, while the southeastern Mediterranean was in excess of CO₂ relative to the atmosphere and released CO₂ to the atmosphere. The water formation areas and, in particular, the northwestern Mediterranean deep-convection area were shown to be regions of relatively strong atmospheric CO₂ uptake (Cossarini et al., 2021). Estimates of air–sea flux for the whole Mediterranean Sea varied between 0.2 Tg C yr⁻¹ (D'Ortenzio et al., 2008) and 2.6 Tg C yr⁻¹ (Cossarini et al., 2021) if only the open seas were considered and between 4.2 Tg C yr⁻¹ (Copin-Montégut, 1993; Cossarini et al., 2021) and 12.6 Tg C yr⁻¹ (Solidoro et al., 2022) when including the continental shelves. Thus the NW Mediterranean deep-convection area, which represents 2.5% of the Mediterranean Sea surface and which we estimate here absorbed at the sea surface 0.4 Tg C yr⁻¹, could strongly contribute to the uptake of atmospheric CO₂ in the open Mediterranean Sea.

Our results show that DIC was transferred from the deep depths and, to a much lesser extent, from the atmosphere to the surface and intermediate water masses and was then transferred laterally to the neighboring sub-basins. The transfer of DIC into surface water masses, which is here estimated at 109 Tg C yr⁻¹, could also mitigate the air–sea CO₂ uptake in winter and spring in the surrounding western and southern seas. It could represent 21 % of the DIC export from the western to the eastern Mediterranean sub-basin estimated by Solidoro et al. (2022) to be 509 Tg C yr⁻¹ and between 8 % to 22 % of the Atlantic CO₂ surface inflow estimated to be between 660 to 1310 Tg C yr⁻¹ by Aït-Ameur and Goyet (2006) and 487 Tg C yr⁻¹ by Solidoro et al. (2022). Finally, the transfer of DIC into intermediate waters, estimated here to be 73 Tg C yr⁻¹, could represent up to 11% to the Mediterranean DIC export towards the Atlantic Ocean at the Gibraltar Strait, estimated to range between 680 and 1380 Tg C yr⁻¹ (Aït-Ameur and Goyet, 2006; Huertas et al., 2009), and up to 100 % of the net DIC outflow (difference between Atlantic surface inflow and Mediterranean outflow), estimated to be between 18 and 70 Tg C yr⁻¹ (Dafner et al., 2001; Aït-Ameur and Goyet, 2006; de la Paz et al., 2009; Huertas et al., 2009). The large discrepancy found here in the percentages of DIC fluxes at the Gibraltar Strait that DIC outflows from the deep-convection area represent is due to the large range in the DIC flux estimates at the Gibraltar Strait found by Aït-Ameur and Goyet (2006). This latter discrepancy is in turn explained by the large range in the estimates of water mass transport used in their computations, which can be attributed to the complex topography and strong tidal currents at the strait (Jordà et al., 2017). In their synthesis study on the exchanges at the Gibraltar Strait, Jordà et al. (2017) found a mean value of the water flux estimates based on observations close to that proposed by Baschek et al. (2001), from which Aït-Ameur and Goyet (2006) derived the lower values of DIC fluxes. This suggests that the higher percentage values that we calculated would be more reliable.

Our results for the northwestern deep-convection area could be compared to those obtained in one of the other major deep-water formation areas of the Mediterranean Sea, the Adriatic Sea. The latter has been shown to be a sink of atmospheric CO₂ (Cossarini et al., 2021) and a sequestration region of anthropogenic carbon (Krasakopoulou et al., 2011; Palmiéri et al., 2015; Hassoun et al., 2015; Ingrosso et al. 2017), as is the case with the study area (Touratier et al., 2016). In particular, experimental studies showed that the deep layer of the southern Adriatic Sea was occupied by dense water rich in DIC and anthropogenic carbon formed in the deep-convection regions of the southern Adriatic Pit and Pomo Pit, as well as on the northern shelf (Krasakopoulou et al., 2011; Cantoni et al, 2016; Ingrosso et al. 2017). The deep dense waters could then be transferred towards the Ionian Sea and the Mediterranean general deep circulation. Krasakopoulou et al. (2011) deduced from in situ measurements over February 1995 inorganic carbon fluxes crossing the Otranto Strait, which connects the Ionian Sea to the southern Adriatic Sea. They estimated that, on an annual basis, the Adriatic Sea could act as a sink of 314 Tg C yr⁻¹ of dissolved inorganic carbon for the Ionian Sea. This net flux resulted from an inflow of 1563 Tg C yr⁻¹, with 27 % in the Levantine Intermediate Water. and an outflow of 1249 Tg C yr⁻¹, with 21 % in the Adriatic Deep Water. Thus, the northwestern Mediterranean deep-convection region and the southern Adriatic, which includes shallower areas, could have opposite contributions in the deep and intermediate layers of the Mediterranean general circulation. However, our DIC budget assessment is limited to a single year and will need to be extended to a longer period to investigate, in particular, the question of carbon sequestration. The deep-convection process exhibits strong interannual variability related to that of air–sea heat flux intensity and the hydrological properties of water masses (Houpert et al., 2016; Somot et al., 2016; Estournel et al., 2016; Margirier et al., 2020). The interannual variability in the intensity and extent of deep convection is expected to directly impact the vertical supply of respired carbon into the upper layer. Furthermore, through variability in temperature and in the vertical supply of dissolved inorganic nutrients into the upper layer, it should also alter biogeochemical and air–sea fluxes of CO₂.

Along with interannual variability, the Mediterranean Sea is experiencing changes in atmospheric forcing and water mass properties in response to global warming and rising atmospheric CO₂ levels. The findings of Merlivat et al. (2018) based on pCO₂ measurements over the two periods 1995–1997 and 2013–2015 suggested an increasing trend in surface DIC concentration. Based on observations from EMSO mooring sites and MOOSE cruises, Coppola et al. (2020) confirmed this trend and also evidenced increasing trends in intermediate (300–800 m) and deep (>2000 m) waters of the Ligurian Sea over the period 1998–2016. The increasing trend in surface DIC concentration could be explained by the trend in atmospheric pCO₂ and is also expected to be influenced by biogeochemical changes in water masses flowing in at the Gibraltar Strait (Merlivat et al., 2018), as well as over the entire Mediterranean Basin, especially in the intermediate- and deep-water formation areas of the eastern basin (Wimart-Rousseau et al., 2021).

Finally, the reduction of winter mixing and the intensification of marine heat waves predicted by models in the second half of the 21st century (Darmaraki et al., 2019; Soto-Navarro et al., 2020) should clearly modify the contribution of the NW deep-convection zone in the Mediterranean. Based on coupled models over the entire Mediterranean Sea, Solidoro et al. (2022) and Reale et al. (2022) predicted, in response to the increase in atmospheric pCO₂, temperature, and stratification, an increase in atmospheric CO₂ uptake and in the DIC inventory in the whole Mediterranean Sea and modifications of the exchange fluxes between the eastern and western sub-basins. The 3D coupled models clearly constitute useful tools to gain insight into carbon budgets, and multi-model ensemble exercises on these issues, as performed by Friedland et al. (2021) regarding the influence of inorganic nutrient river inputs, could allow a refinement of the carbon budget terms and their evolution together with an assessment of their uncertainties.

6 Conclusions

We have estimated for the first time a CO₂ budget for the whole northwestern Mediterranean deep-convection zone over an annual period using a high-resolution 3D coupled hydrodynamic–biogeochemical–chemical model. An assessment of the model results through comparisons with DEWEX and MOOSE-GE cruise observations, as well as EMSO-DYFAMED mooring and BOUSSOLE buoy site observations and outputs of the CANYON-MED neural networks, shows the ability of the model to describe the seasonal cycle and spatial variability of the DIC dynamics in this region with good accuracy. Based on the present study over the period 2012–2013, we can draw the following conclusions for this key region in the ocean circulation and biogeochemical cycles in the Mediterranean Sea.

The CO₂ dynamics in the NW Mediterranean deep-convection area underwent large seasonal variation. The region was marked by a deficit of CO₂ compared to the atmosphere from November to early June, which led to a 7-month ingassing of atmospheric CO₂. The deficit situation, to a large extent controlled by temperature variability, was, on the one hand, reduced by vertical supply of DIC during the period of deep convection and, on the other hand, accentuated and extended by the spring phytoplankton bloom. This underlines the findings of Mémery et al. (2002) regarding the importance of data of sea surface DIC or pCO₂ data during deep convection for precise estimates of air–sea CO₂ flux in this area in addition to sea surface temperature observations and spring NCP estimates.

On an annual basis, the NW Mediterranean deep-convection area acted as a sink of atmospheric CO₂. We estimate an annual uptake of 0.47 mol C m⁻² yr⁻¹. The maximum fluxes (>1 mol C m⁻² yr⁻¹) occurred in the northern Gulf of Lion region, subject to strong northerly continental winds and located at the edge of the deep convection between the Northern Current and the core of the deep convection, while minimum values (close to null values) are found in the western zone where warm anticyclonic gyres developed. The sensitivity tests of the annual estimate to the parametrization of the gas transfer velocity indicate an uncertainty of 28 %. Calcification processes could lead to an overestimation of 16 % to 57 % of the annual uptake, highlighting the need for the refinement of the model in future studies.

The annual DIC budget in the upper layer of the deep-convection area was co-dominated by biogeochemical and physical fluxes, estimated here to be −3.7 and 3.3 mol C m⁻² yr⁻¹, respectively. The net physical flux resulted from a balance of a net upward transfer and a net lateral export, both exhibiting maximum intensity during the preconditioning and deep-convection period. The air–sea CO₂ flux only represents 13 % of the upper-layer NCP and 31 % of the NCP integrated over the whole water column. These results confirm that the DIC budget in this region should be addressed with a 3D approach that considers the complex physical mechanisms taking place.

The NW Mediterranean deep-convection area acted as a source of DIC for the surface and intermediate water masses flowing towards the southwestern Mediterranean. The transfers of DIC into the adjacent surface and intermediate water masses could also mitigate the atmospheric CO₂ uptake in the surrounding open sea of the sub-basin and represent up to 22 % and 11 %, respectively, of the DIC exchanges with the Atlantic Ocean at the Strait of Gibraltar.