Oxygen budget of the north-western Mediterranean deep- convection region

The north-western Mediterranean deep convection plays a crucial role in the general circulation and biogeochemical cycles of the Mediterranean Sea. The DEWEX (DEnse Water EXperiment) project aimed to better understand this role through an intensive observation platform combined with a modelling framework. We developed a three-dimensional coupled physical and biogeochemical model to estimate the cycling and budget of dissolved oxygen in the entire north-western Mediterranean deep-convection area over the period September 2012 to September 2013. After showing that the simulated dissolved oxygen concentrations are in a good agreement with the in situ data collected from research cruises and Argo floats, we analyse the seasonal cycle of the air–sea oxygen exchanges, as well as physical and biogeochemical oxygen fluxes, and we estimate an annual oxygen budget. Our study indicates that the annual air-to-sea fluxes in the deep-convection area amounted to 20 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. A total of 88 % of the annual uptake of atmospheric oxygen, i.e. 18 mol m⁻², occurred during the intense vertical mixing period. The model shows that an amount of 27 mol m⁻² of oxygen, injected at the sea surface and produced through photosynthesis, was transferred under the euphotic layer, mainly during deep convection. An amount of 20 mol m⁻² of oxygen was then gradually exported in the aphotic layers to the south and west of the western basin, notably, through the spreading of dense waters recently formed. The decline in the deep-convection intensity in this region predicted by the end of the century in recent projections may have important consequences on the overall uptake of atmospheric oxygen in the Mediterranean Sea and on the oxygen exchanges with the Atlantic Ocean, which appear necessary to better quantify in the context of the expansion of low-oxygen zones.

1 Introduction

Deep convection is a key process leading to a massive transfer of oxygen from the atmosphere to the ocean interior (Körtzinger et al., 2004, 2008b; Fröb et al., 2016; Wolf et al., 2018). Its weakening in some regions (de Lavergne et al., 2014; Brodeau and Koenigk, 2016), induced by enhanced stratification, is one of the primary factors, along with changing ventilation at intermediate depths, slowdown of the overturning circulation, warming-induced decrease in solubility modulated by salinity changes and changes in C:N utilisation ratios, that may explain the ongoing decline in the open-ocean oxygen inventory, or deoxygenation, observed and modelled since the middle of the 20th century (Bopp et al., 2002; Keeling and Garcia, 2002; Plattner et al., 2002; Joos et al., 2003; Keeling et al., 2010; Helm et al., 2011; Andrews et al., 2017; Ito et al., 2017; Schmidtko et al., 2017; Breitburg et al., 2018). The oxygen decline leads to an increase in the volume of hypoxic or even anoxic waters and to the expansion of oxygen minimum zones (OMZs), which substantially affect life and habitats of marine ecosystems and have implications for biogeochemical cycles (Ingall et al., 1994; Levin, 2003; Diaz and Rosenberg, 2008; Breitburg et al., 2009; Naqvi et al., 2010; Stramma et al., 2010; Scholz et al., 2014; Bristow et al., 2017). It is crucial to gain understanding of the actual ventilation occurring in deep-convection areas and to continue developing models to predict its future evolution under climate change.

The north-western (NW) Mediterranean Sea (Fig. 1a and b) is one of the few regions of the world where deep convection takes place (Schott et al., 1996). In autumn, during the preconditioning phase, a cyclonic gyre formed by the Northern Current and the Balearic front leads to the doming of the isopycnals and the rising of high-salinity intermediate waters, the Levantine Intermediate Waters (LIW), close to the surface. In winter, cold and dry northerly winds (Mistral, Tramontane) produce the cooling, evaporation and subsequent density increase in surface waters. The instability of the water column induces convective mixing of surface waters with deeper waters, and, when the process is intense, in the formation of new deep waters that spread into the western Mediterranean Sea, such as observed in 2004–2006 (Schroeder et al., 2008a). The depth and horizontal extension of convection in the NW region show strong interannual variability, driven by both the variability of the winter buoyancy loss and the stratification magnitude prior to the convection period (Mertens and Schott, 1998; Béthoux et al., 2002; Houpert et al., 2016; Somot et al., 2016). The deepest convection takes place in the centre of the Gulf of Lion where the yearly maximum of the mixed-layer depth varies from a few hundred metres to 2500 m when the bottom is reached (MEDOC Group, 1970).

Figure 1 Domain and bathymetry (m) of (a) the forcing coupled NEMO-Eco3MS Mediterranean model (red contour) and of (b) the coupled SYMPHONIE-Eco3MS western sub-basin model. The blue contour in (a) indicates the limits of the western sub-basin model. The red contoured area in (b) corresponds to the deep-convection area. (c) Scheme of the downscaling strategy from the Mediterranean Sea to the western sub-basin.

At the Mediterranean basin scale, the deep convection occurring in the north-western region is one of the major processes responsible for an enrichment of the euphotic layer with nutrients, compared to Atlantic influx as well as terrestrial and atmospheric inputs (Severin et al., 2014; Ulses et al., 2016; Kessouri et al., 2017). The replenishment of the surface layer with nutrients during deep convection is followed by an intense bloom in spring when vertical mixing weakens (Bernadello et al., 2012; Lavigne et al., 2013; Auger et al., 2014; Ulses et al., 2016; Kessouri et al., 2018). The formation of deep waters is also at the origin of a huge ventilation of the western Mediterranean Sea (Minas and Bonin, 1988; Copin-Montégut and Bégovic, 2002; Schroeder et al., 2008a; Schneider et al., 2014; Stöven and Tanhua, 2015; Touratier et al., 2016; Coppola et al., 2017; 2018; Li and Tanhua, 2020). The study of Mavropoulou et al. (2020), based on in situ observations over the period 1960–2011, indicated that its variability is one of the main drivers of the interannual variability of the dissolved oxygen (O₂) concentration in the deep waters of the western Mediterranean Sea. The oxygenation induced by recurrent intermediate and deep convection together with a relatively low primary production make the Mediterranean Sea a well-oxygenated basin (Tanhua et al., 2013). In the western Mediterranean open sea, the oxygen minimum layer (OML) is located in the LIW and shows a minimum oxygen concentration of 170–185 µmol kg⁻¹, above ≈70 % of the saturation levels (Tanhua et al., 2013; Coppola et al., 2018). Thus the OML in this region is clearly less pronounced than the OMZs in the open oceans or deep basins of other seas, such as the adjacent Black Sea, where hypoxic conditions (oxygen concentration <2 mL O₂ L⁻¹ or <61 µmol O₂ kg⁻¹; Diaz and Rosenberg, 2008; Breitburg et al., 2018) are encountered. However the semi-enclosed Mediterranean Sea with a fast warming was identified as one of the most vulnerable marine regions to climate change (Giorgi, 2006). Recently, regional ocean models of the Mediterranean Sea converged to predict a weakening of NW deep-convection intensity under climate change scenarios by the end of the 21st century (Soto-Navarro et al., 2020). Yet, Coppola et al. (2018), by analysing the evolution of observed oxygen profiles in the Ligurian Sea over a 20 year period, suggested that hypoxic conditions may be reached in water masses at intermediate depths after a period of 25 years without deep-convection events (presuming bacterial respiration remains the same).

One of the objectives of the DEWEX (DEnse Water EXperiment) project carried out in 2012/13 was to investigate the deep-convection process, the formation of north-western Mediterranean Deep Waters and the impact of deep convection on biogeochemical fluxes (Conan et al., 2018). Three cruises and the deployment of autonomous platforms (glider, Argo floats) provided an unprecedented intensive observation of this region before, during and after a deep-convection event and completed the observation effort during the stratified period operated since 2010 in the framework of the MOOSE-GE (Mediterranean Ocean Observing System for the Environment-Grande Échelle) programme (Estournel et al., 2016b). The 2012/13 event was identified by observational and modelling studies as one of the five most intense deep-convection events over the period 1980–2013 (Somot et al., 2016; Herrmann et al., 2017; Coppola et al., 2018) due to extremely strong buoyancy loss (Somot et al., 2016). Regarding the oxygen dynamics, DEWEX winter observations showed a strong increase in the O₂ inventory of the entire water column, which was concomitant to the deepening of the mixed layer and was attributed to a rapid intake of atmospheric dissolved oxygen (Coppola et al., 2017). However, these observations remain limited in time and space. Up to date, no high-resolution modelling of the oxygen dynamics in the NW deep-convection region that could complete the monitoring effort and provide quantification for the whole area has yet been proposed.

In this study, we take advantage of the DEWEX project to implement and constrain with in situ observations a 3D coupled physical–biogeochemical model representing the dynamics of dissolved oxygen and to gain understanding in the variability of the oxygen inventory in the whole NW Mediterranean deep-convection area, for the period between September 2012 and September 2013. In this framework, we investigate the seasonal cycle of the oxygen inventory and estimate its annual budget, and we analyse and quantify the relative contribution of air–sea exchanges, as well as of physical and biogeochemical processes in the budget. The following document is organised as follows: in Sect. 2, we describe the numerical model, its implementation and the observations used for its assessment. In Sect. 3, we compare our model results with in situ observations. In Sect. 4, we describe the seasonal cycle of atmospheric and physical conditions. In Sect. 5, we examine the seasonal cycle of oxygen inventory and fluxes, as well as the annual oxygen budget. We discuss our results in Sect. 6 and conclude in Sect. 7.

2 Material and methods

2.1 The numerical model

2.1.1 The hydrodynamic model

The SYMPHONIE model used in this study is a 3D primitive equation model, with a free surface and generalised sigma vertical coordinate, as described in Marsaleix et al. (2008). The vertical diffusion is parameterised with a prognostic equation for the turbulent kinetic energy and a diagnostic equation for the mixing and dissipation lengths, following Gaspar et al. (1990). Atmospheric forcing (turbulent fluxes) is calculated using the bulk formulae described by Large and Yeager (2004). This model was previously used in the Mediterranean Sea to simulate open-sea convection (Herrmann et al., 2008; Estournel et al., 2016a; Ulses et al., 2016), shelf dense-water cascading (Estournel et al., 2005; Ulses et al., 2008b) and continental shelf circulation on the Gulf of Lion shelf (Estournel et al., 2001; 2003; Ulses et al., 2008a).

2.1.2 The biogeochemical model

The biogeochemical model Eco3M-S is a multi-nutrient and multi-plankton functional type model that simulates the dynamics of the biogeochemical decoupled cycles of several biogenic elements (carbon, nitrogen, phosphorus, silicon) and of non-Redfieldian plankton groups (Ulses et al., 2016). The model was previously used to study the biogeochemical processes on the Gulf of Lion shelf (Auger et al., 2011) and in the NW Mediterranean deep-convection area (Herrmann et al., 2013, 2017; Auger et al., 2014; Ulses et al., 2016; Kessouri et al., 2017, 2018). In this study, the model was extended to describe the dynamics of dissolved oxygen in the ocean interior and the air–sea exchanges of oxygen. Here we only describe the rate of change of the new state variable, the dissolved oxygen concentration and the parameterisation of the air–sea flux of oxygen, which were included in the model version described in detail by Auger et al. (2011). The rate of change of dissolved oxygen concentration due to biogeochemistry in the water column is governed by the following equation:

$$\begin{array}{ll}\text{(1)}& {\displaystyle \frac{\mathrm{d}{\mathrm{DO}}_x}{\mathrm{d}t}}=& {\displaystyle }\sum _{i=\mathrm{1}}^{\mathrm{3}}\left({\text{GPP}}_i-{\text{RespPhy}}_i\right){\mathit{\gamma }}_{\mathrm{C}/{\mathrm{DO}}_x}-{\displaystyle }& {\displaystyle }\sum _{i=\mathrm{1}}^{\mathrm{3}}\left({\text{RespZoo}}_i+{\text{RespZoo}}_i^{\text{add}}\right){\mathit{\gamma }}_{\mathrm{C}/{\mathrm{DO}}_x}-\\ {\displaystyle }& {\displaystyle }\text{RespBac}{\mathit{\gamma }}_{\mathrm{C}/{\mathrm{DO}}_x}-\text{Nitrif}{\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}/{\mathrm{DO}}_x},\end{array}$$

where DO_x is the dissolved oxygen concentration; GPP_i and RespPhy_i are gross primary production and respiration, respectively, for phytoplankton group i; RespZoo_i and ${\text{RespZoo}}_i^{\text{add}}$ are basal respiration and additional respiration fluxes to maintain constant N:C and P:C internal ratios, respectively, for zooplankton group i; RespBac is bacterial respiration; and Nitrif is nitrification. ${\mathit{\gamma }}_{\mathrm{C}/{\mathrm{DO}}_x}$ and ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}/{\mathrm{DO}}_x}$, equal to 1 and 2, respectively, are the mole of DO_x, used per mole of C in respiration and needed to oxidise 1 mol of ammonium in nitrification as described in Grégoire et al. (2008).

The flux of dissolved oxygen at the air–sea interface, DO_xFlux, is computed from

$$\begin{array}{}\text{(2)}& {\mathrm{DO}}_x\text{Flux}=K_w\left({{\mathrm{DO}}_x}_{\text{sat}}-{{\mathrm{DO}}_x}_{\text{surf}}\right),\end{array}$$

where ${{\mathrm{DO}}_x}_{\text{sat}}$ and ${{\mathrm{DO}}_x}_{\text{surf}}$ (mmol m⁻³) are the concentration of dissolved oxygen at saturation level and at the surface of the ocean, respectively, and K_w(m s⁻¹) is the gas transfer velocity. The oxygen solubility (or dissolved oxygen at saturation level) is determined using the equation of Garcia and Gordon (1992). The oxygen saturation anomaly (noted ΔO₂) is defined as $\mathrm{\Delta }{\mathrm{O}}_{\mathrm{2}}=({\mathrm{DO}}_x-{{\mathrm{DO}}_x}_{\text{sat}})/{{\mathrm{DO}}_x}_{\text{sat}}\times \mathrm{100}\mathit{\%}$. We computed here the gas transfer velocity using the parameterisation of Wanninkhof and McGillis (1999) with a cubic dependency to the wind speed, following the study in the convective Labrador Sea by Körtzinger et al. (2008b), who found that this parameterisation was one of those that gave the best results and recommended a stronger-than-quadratic wind speed dependency for the high wind speed range. In addition, sensitivity analyses using eight various parameterisations of the gas transfer velocity were performed to estimate uncertainties of air–sea exchanges and are discussed in Sect. 6.1. For these sensitivity tests, we used quadratic (Wanninkhof, 1992, 2014) and hybrid (Nightingale et al., 2000; Wanninkhof et al., 2009) wind speed dependency parameterisations, as well as parameterisations including air–sea fluxes due to bubble formation (Woolf, 1997; Stanley et al., 2009; Liang et al., 2013; Bushinsky and Emerson, 2018).

2.1.3 Implementation

The implementation of the hydrodynamic and biogeochemical simulations used in this study was described by Estournel et al. (2016a) and Kessouri et al. (2017; 2018). The numerical domain (Fig. 1a, b) covers most of the western Mediterranean basin, using a curvilinear grid with variable horizontal resolution (Bentsen et al., 1999). The mesh size ranges from 0.8 km in the north to 1.4 km in the south. The grid has 40 vertical levels with closer spacing near the surface (15 levels in the first 100 m in the centre of the convection area characterised by depths of ≈2500 m). As explained in Estournel et al. (2016a), the size of the grid is not small enough to explicitly represent convective plumes, which thus need to be parameterised. In our case, to prevent the development of static instabilities at the surface resulting in noise at the scale of the mesh, the heat and water fluxes are distributed over the whole mixed layer whose thickness is given by the depth at which the vertical density gradient becomes negative.

The biogeochemical model is forced offline by daily outputs of the hydrodynamic model. The advection and diffusion of the biogeochemical variables were calculated using the QUICKEST (QUICK with Estimated Streaming Terms) scheme (Leonard, 1979) on the horizontal and with a centred scheme on the vertical.

A strategy of downscaling from the Mediterranean basin to the western sub-basin scale was implemented in three stages (Fig. 1a and c) as described by Kessouri et al (2017). In a first step, the SYMPHONIE hydrodynamic model was initialised and forced at its lateral boundaries with daily analyses of the configuration PSY2V4R4, based on the NEMO ocean model at a resolution of 1/12^∘ over the Atlantic and the Mediterranean Sea by the Mercator Ocean International operational system (Lellouche et al., 2013). Second, the biogeochemical model was forced at the Mediterranean basin scale by the outputs of the same NEMO simulation. In a third step, the daily outputs of the two previous simulations were used to initialise and force the Eco3M-S biogeochemical model over the western Mediterranean Sea. This nesting protocol ensures the coherence of the physical and biogeochemical fields at the open boundaries. The basin configuration of the biogeochemical model was initialised in summer 2011, with climatological fields of in situ nutrient concentrations from the oligotrophic period in the MEDAR/MEDATLAS database (Manca et al., 2004) and according to oxygen observations from the Meteor M84/3 cruise carried out in April 2011 (Tanhua et al., 2013) and DYFAMED station observations in August 2011 (Coppola et al., 2018). The regional biogeochemical simulation started in August 2012. Due to strong vertical diffusivities in the basin-scale model, we corrected the initial oxygen concentration for the north-western region using DYFAMED observations carried out in the summer of 2012 (Coppola et al., 2018) and for the south-western region according to Meteor M84/3 observations.

Meteorological parameters including downward radiative fluxes were taken from the ECMWF (European Centre for Medium-Range Weather Forecasts) operational forecasts at 1/8^∘ horizontal resolution and 3 h temporal resolution based on daily analyses. River runoffs were considered based on realistic daily values for French rivers (data provided by Banque Hydro, http://www.hydro.eaufrance.fr/) and Ebro (data provided by SAIH Ebro, http://www.saihebro.com), and mean annual values for the other rivers. At the Rhône River mouth, nitrate, ammonium, phosphate, silicate and dissolved organic carbon concentrations were prescribed using in situ daily data (P. Raimbault, MISTRALS database). These data, combined with those of Moutin et al. (1998) and Sempéré et al. (2000), were used to estimate dissolved organic phosphorus and nitrogen and particulate organic matter concentrations as described in Auger et al. (2011). At the other river mouths, climatological values were prescribed according to Ludwig et al. (2010). Dissolved oxygen concentration at the river mouths was set to values at saturation. The deposition of organic and inorganic matter from the atmosphere was neglected in this study. Fluxes of inorganic nutrients and oxygen at the sediment–sea interface were considered by coupling the pelagic model with a simplified version of the meta-model described by Soetaert et al. (2001). The parameters of the latter model were set following the modelling study performed by Pastor et al. (2011) for the Gulf of Lion shelf.

2.1.4 Study area

For analyses and budget purposes, we defined the deep-convection area as the area where the daily averaged mixed-layer depth exceeded 1000 m at least once during wintertime (red contoured area in Fig. 1b), according to Kessouri et al. (2017, 2018). It covered an area of 61 720 km² in 2013. The mixed-layer depth is defined as the depth where the potential density exceeds its value at 10 m depth by 0.01 kg m⁻³ (Coppola et al., 2017). Heat fluxes and physical and biogeochemical parameters and fluxes presented in the following sections correspond to values averaged over all model grid points included in this area. The budget of oxygen inventory was computed in two layers based on biogeochemistry processes: in the upper layer (from the surface to 150 m) including the euphotic layer where photosynthesis influences the dynamics of oxygen and in the underlying aphotic layer (from 150 m to the bottom) where only respiration and nitrification processes are taken into account in the model. The maximum depth of the base of the euphotic layer was defined at 150 m, based on the regional minimum value of diffuse attenuation coefficient of light at 490 nm derived from satellite observations (http://marine.copernicus.eu/, products: OCEANCOLOUR_MED_OPTICS_L3_REP_OBSERVATIONS_009_095), and following the studies by Lazzari et al. (2012) and Kessouri et al. (2018).

2.2 Observations used for the model assessment

2.2.1 Cruise observations

To assess the horizontal and vertical distribution of the simulated dissolved oxygen concentration, we use 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 one, DEWEX Leg1, was carried out during the active phase of deep convection, in February 2013 (Testor, 2013), and the second one, DEWEX Leg2, was carried out during the following spring bloom, in April 2013 (Conan, 2013). In addition, we use 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). The dissolved oxygen measurements were performed during the DEWEX (Leg1: 74 stations, Leg2: 99 stations) and MOOSE-GE (74 stations) cruises, using a Sea-Bird SBE43 sensor. The calibration and quality control of the measurements were described by Coppola et al. (2017). The accuracy of the measurements was estimated at 2 % of oxygen saturation, i.e. 4 µmol kg⁻¹. A Winkler analysis performed on board was used to adjust the SBE43 raw data, as specified by the GO-SHIP programme (http://www.go-ship.org/).

We also compare our model results with high-frequency measurements of wind at 10 m and of ocean surface temperature, salinity (thermosalinograph and conductivity–temperature–depth (CTD)) and dissolved oxygen concentration (optode) at 3 m depth using the sea surface water continuous acquisition system (SACES) (Dugenne, 2017) during the two DEWEX cruises.

2.2.2 Argo floats

To evaluate the temporal evolution of the modelled oxygen inventory, we use data of three Argo-O₂ floats (floats 6901467, 6901470, 6901487) deployed in the NW Mediterranean Sea during the preconditioning phase (late November 2012) and the active phase (late January 2013) of dense water formation and operational until the end of the study period (Coppola et al., 2017; Argo, 2020). Dissolved oxygen measurements were made with a standard CTD sensor, equipped with an oxygen optode with fast time response (Aanderaa 4330). Calibrations of optodes were performed before the float deployment and also during the deployment using CTD profiles and seawater samples (Niskin bottles). Details on float deployment strategy and calibration are given by Coppola et al. (2017). We calculated the oxygen inventory from 1800 m to the surface for floats 6901467 and 6901470 and only from 1000 m to the surface for float 6901487 due to poor quality salinity data below this depth.

2.3 Statistical analysis

In order to quantify the performance of the model in its ability to represent the dynamics of dissolved oxygen for the study period, we computed four complementary metrics following the recommendations of Allen et al. (2007): (1) the SD ratio ($r_{\mathit{\sigma }}=\frac{{\mathit{\sigma }}_o}{{\mathit{\sigma }}_m}$ where σ_m and σ_o are the SD of model outputs and observations, respectively); (2) the Pearson correlation coefficient,

$$R={\displaystyle \frac{\frac{\mathrm{1}}{K}\sum _{k=\mathrm{1}}^K\left(y_k^m-\stackrel{\mathrm{‾}}{y^m}\right)\left(y_k^o-\stackrel{\mathrm{‾}}{y^o}\right)}{{\mathit{\sigma }}_m{\mathit{\sigma }}_o}},$$

where K is the number of observations, $y_k^o$ is the observation k, $y_k^m$ is the corresponding model output k, and $\stackrel{\mathrm{‾}}{y^o}$ and $\stackrel{\mathrm{‾}}{y^m}$ are the mean of observations and model outputs, respectively; (3) the normalised root-mean-square error,

$$\text{NRMSE}={\displaystyle \frac{\sqrt{\frac{\mathrm{1}}{K}\sum _{k=\mathrm{1}}^K{\left(y_k^o-y_k^m\right)}^{\mathrm{2}}}}{\stackrel{\mathrm{‾}}{y^o}}}\times \mathrm{100}\mathit{\%};$$

and (4) the percentage bias ($\text{PB}=\frac{\stackrel{\mathrm{‾}}{y^m}-\stackrel{\mathrm{‾}}{y^o}}{\stackrel{\mathrm{‾}}{y^o}}\times \mathrm{100}\mathit{\%}$). The model results are compared with the observations at the same dates and positions.

3 Evaluation of the model

The accurate representation of the winter mixing of water masses is an essential point for the simulation of the dissolved oxygen dynamics in this region, marked by a strong ventilation of the deep waters that plays a crucial role in its seasonal cycle (Copin-Montégut et al., 2002; Touratier et al., 2016; Coppola et al., 2017; 2018). A validation of the hydrodynamic part of the simulation is described by Estournel et al. (2016a), who showed similar spatial distribution of the modelled water column stratification in the entire deep-convection area, as well as modelled time evolution of the temperature profile in the centre of the Gulf of Lion open sea during the winter that is close to the observations.

Furthermore, an assessment of the biogeochemical part of the coupled model is presented in Kessouri et al. (2017; 2018). These studies showed that the model is able to accurately reproduce the timing and magnitude of the surface chlorophyll increase during the spring and autumnal blooms, as well as the concentrations of nutrients and depths of nutriclines and the dynamics and depth of the deep chlorophyll maximum during the stratified, oligotrophic period.

In this study, we focus the evaluation of the coupled model on its ability to realistically represent the dynamics of dissolved oxygen in the deep regions of the NW Mediterranean Sea. For this purpose, we first compare the model results to in situ observations from DEWEX and MOOSE-GE cruises conducted at three key periods: the winter mixing period, the phytoplankton bloom period and the stratified summer period. We then compare the model outputs to Argo data deployed in the area in terms of time evolution of oxygen inventory.

3.1 Comparisons to cruise observations

The comparisons of modelled wind velocity and ocean model outputs with in situ observations from the high-frequency SACES are shown in Fig. 2. Modelled wind provided by ECMWF and used to force the hydrodynamic model and to calculate the air–sea oxygen flux is highly correlated with the observations (R=0.96, p value <0.01). The low values of NRMSE (13.9 %) and percentage bias (−0.5 %) show the accuracy of this variable, found for all ranges of value. Regarding the surface ocean variables, we obtain statistically significant correlations equal to 0.64, 0.83 and 0.83 (p value <0.01), between observed and modelled values of, respectively, surface temperature, salinity and dissolved oxygen concentration. The NRMSEs are equal to 2.0 %, 0.3 % and 5.2 %, respectively. The percentage biases remain negligible for temperature (−0.7 %), salinity (0.002 %) and dissolved oxygen concentration (−1.2 %).

Figure 2 Time evolution during DEWEX Leg1 (in February 2013, a, c, d and e) and Leg2 (in April 2013, b, f) cruises of observed (red) and modelled (blue) (a and b) wind velocity (m s⁻¹), (c) surface temperature (^∘C), (d) surface salinity, and (e, f) surface dissolved oxygen concentration (µmol kg⁻¹). Trajectories of the measurements during DEWEX Leg1 and Leg2 cruises are indicated on inserted maps. Modelled wind velocity was provided by ECMWF. No surface temperature and salinity data are available over the period of DEWEX Leg2. The metrics indicated for the modelled wind velocity and surface oxygen concentration were calculated for both DEWEX Leg1 and Leg2 periods.

Figures 3 and 4 compare the observed and modelled dissolved oxygen concentration for the stations sampled during the DEWEX and MOOSE-GE cruises, respectively, at the surface (between 5 and 10 m) and along the south–north transect passing across the convection area (stations encircled in black in Fig. 3). Overall, the simulation correctly reproduces the spatial and temporal variability of the oxygen concentration observed at the surface and in the water column during and between the three cruises. During wintertime, the model simulates low surface oxygen concentrations (<220 µmol kg⁻¹) in the open sea of the Gulf of Lion and the Ligurian Sea, areas that coincide with the deep vertical mixing regions (Estournel et al., 2016a; Kessouri et al., 2017) (Fig. 3a and b). Figure 4a shows the oxygen homogenisation of the whole water column between 41.5 and 42.3^∘ N, the core of the deep-convection area. Concentrations above 240 µmol kg⁻¹ are modelled in the surface layer on the shelf and in the south at the Balearic front, in accordance with the observations (Figs. 3a, b and 4a). The model also agrees with observations showing a layer of low oxygen concentration (minimum concentration <185 µmol kg⁻¹ at depths around 500 m) located between 150 and 1500 m, mainly in the LIW (300–800 m), outside the deep-convection area (Fig. 4a). The metrics confirm the good agreement between model outputs and observations with a significant spatial correlation of 0.81 and 0.61 (p value <0.01), an NRMSE of 5.3 % and 15.7 %, and a negligible percentage bias of −1.1 % and 0.01 %, respectively, at the surface and along the south–north transect.

Figure 3 Surface dissolved oxygen concentration (µmol kg⁻¹) observed (a, c and e) and modelled (b, d and f) over the (a and b) DEWEX Leg1 (1–21 February 2013), (c and d) DEWEX Leg2 (5–24 April 2013) and (e and f) MOOSE-GE (11 June–9 July 2013) cruise periods. The black-outlined dots correspond to the measurement stations shown in Fig. 4.

Figure 4 Comparison between model outputs and observations on a transect crossing the deep-convection area (stations are circled in black in Fig. 3a, c and e) over the (a and b) DEWEX (Leg1: 10–12 February 2013, Leg2: 8–10 April) and (c) MOOSE-GE (27 June–5 July 2013) cruise periods. The left column shows observed and modelled profiles at 42 and 40.3^∘ N. The right column shows vertical section of dissolved oxygen concentration (µmol kg⁻¹) along the transect; the model is represented by background colours and observations are indicated in coloured circles.

During the spring cruise period, the model represents high dissolved oxygen values (>240 µmol kg⁻¹) at the surface throughout the region, as observed (Fig. 3c and d). The increase in modelled oxygen concentration in the surface layer between both campaigns is in agreement with observations (Figs. 3a–d and 4a, b). A zone of low oxygen concentration in the intermediate waters is present in the convection area in both datasets (Fig. 4b). However, it is worth noting that this zone of low oxygen concentration is heterogeneous in its magnitude and thickness both in model outputs and in observations, and it is not similarly distributed in space in the model compared to the measurements. At the surface and along the transect, the spatial correlation coefficients between modelled and observed dissolved oxygen are equal to 0.59 and 0.30 (p value <0.01), respectively, the NRMSE to 4.6 % and 19.3 %, respectively, and the percentage biases to −1.2 % and −4.4 %, respectively.

The north–south gradient, with lower surface concentrations in the south of the deep-convection area, observed during the stratified period (i.e. MOOSE-GE cruise period in June–July), is then well reproduced by the model (Fig. 3e and f). The minimum zone is more established than in spring in both in situ data and model results (Fig. 4c). Both sets of data represent a maximum in the subsurface at depths around 50 m, close to the deep chlorophyll maximum (shown in Fig. 5 in Kessouri et al., 2018), although an underestimation of its magnitude is visible between 41.5 and 42^∘ N in the model (Fig. 4c). We find a spatial correlation coefficient of 0.64 and 0.96 (p value <0.01), an NRMSE of 3.2 % and 3.5 %, and a negligible percentage bias (absolute values ≤0.4 %) between model outputs and observations at the surface and along the north–south transect, respectively.

Figure 5 The left column shows oxygen inventory integrated from the surface to 1800 m depth (a and b) or 1000 m (c) (mol m⁻²) in Argo float measurements (cyan) and model outputs (black) along Argo (a) 6901467, (b) 6901470 and (c) 6901487 float trajectories. The right column shows trajectories of corresponding Argo float.

The metrics computed using all station data from the three cruises are given in Table 1. The modelled dissolved oxygen concentration is significantly correlated with the observed concentration (R≥0.81, p<0.01), in particular for the winter period when the pattern of the oxygen distribution appears to be primarily shaped by deep-convection processes, shown to be accurately represented by Estournel et al. (2016a). The model results show low percentage biases (PB <1 %), low NRMSE (<8 %), and SD ratios ranging between 1.13 and 1.35, which indicate a larger variability in the observations than in the model outputs.

Table 1 Statistical analysis of model results: Pearson correlation coefficient (R), normalised root-mean-square error (NRMSE), percentage bias (PB) and SD ratio, calculated between modelled dissolved oxygen concentrations and observations from DEWEX winter and spring cruises and the MOOSE-GE summer cruise and from Argo-O₂ platforms.

3.2 Comparison to Argo float data

The model accurately reproduces the magnitude of oxygen inventory in the water column and its time evolution observed using Argo floats during the study period (Fig. 5). The model simulates the increase observed between early December and late February. This increase is estimated at ≈20 mol m⁻² over a layer from the surface to 1800 m, along the trajectory of the float 6901467 (Fig. 5a), and at ≈10 mol m⁻² over a layer from the surface to 1000 m, along the trajectory of the float 6901487 (Fig. 5c), both floats being located in the Gulf of Lion at that period. The oxygen inventory remains high during the month of March and then decreases significantly from early April to early June, in model outputs and Argo observations. In both datasets, the decrease reaches up to 20 mol m⁻² over 1800 m along the path of the Argo float 6901470 in the Gulf of Lion (Fig. 5b) and is less pronounced (≈10 mol m⁻²) along the trajectory of the float 6901467 in the Balearic Sea (Fig. 5a). More moderate decreases are then simulated and observed until September along all float trajectories. The statistical analysis shows that, in terms of oxygen inventory, significant correlation coefficients are obtained between the model outputs and the three float observations ($\mathrm{0.64}<R<\mathrm{0.83}$, p value <0.01), NRMSEs are smaller than or equal to 2.4 %, and the absolute values of percentage bias are smaller than or equal to 2 % (Fig. 5). In terms of dissolved oxygen concentration in the water column, we obtain significant correlation coefficients ($\mathrm{0.56}<R<\mathrm{0.93}$, p value <0.01), NRMSEs smaller than 10.5 %, percentage biases smaller than 1 %, and SD ratios close to 1 for floats 6901470 and 6901487 and of 1.37 for float 6901467 (Table 1).

4 Atmospheric and hydrodynamic conditions

In the NW Mediterranean Sea, deep convection takes place every winter but shows strong interannual variability in its magnitude and spatial extent. This interannual variability is partly related to the variability of heat fluxes (Somot et al., 2016). Over the study period, the convection area was marked by severe heat loss episodes from late October 2012 to mid-March 2013 (Fig. 6a). In particular, there was a first short but intense heat loss event (mean heat flux $<-\mathrm{1000}$ W m⁻²) at the end of October, followed by several long northerly wind episodes when heat loss peaks reached 500 W m⁻², during the months of December to February (late November to mid-December, mid-January, early February and late February). Finally, a last strong heat loss episode occurred in mid-March after a period of positive heat flux. The wind velocity averaged over the convection period (15 January–8 March, 15–24 March) was maximum in the centre of the Gulf of Lion, where it reached 10 m s⁻¹ (Fig. 7a). From April onwards, the convection region was mainly characterised by heat gains.

Figure 6 Time series of modelled (a) total heat fluxes (W m⁻²), (b) mixed-layer depth (m), (c) surface oxygen and oxygen solubility (µmol kg⁻¹), (d) air-to-sea oxygen fluxes ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), (e) downward oxygen transport at 150 m (dark blue) and lateral oxygen transport towards the convection area (light blue) ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), and (f) biogeochemical oxygen production (see Eq. 1) ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), spatially averaged over the convection area (spatial mean in solid line and shaded area for SD). Sources: ECMWF for heat fluxes, SYMPHONIE/Eco3M-S for the other parameters and fluxes. The blue shaded area corresponds to the deep-convection period (period when spatially averaged MLD >100 m). Note that the range of the y axis varies for the different oxygen fluxes, and due to higher values SD for vertical and lateral transport is not shown.

In response to the autumnal heat loss events, the mixed layer (ML) began to deepen below 50 m at the end of November (Fig. 6b). Its deepening was strongly enhanced in winter over four periods that coincided with the four episodes of intense northerly wind associated with heat loss mentioned above (Fig. 6a and b). Deep convection reached the bottom layer (≈2000 m) in the core of the convection zone (latitude ≈42^∘ N, 4^∘ E < longitude <5^∘ E) in early February, and the spatially averaged mixed layer reached a maximum depth of about 1500 m at the end of February (Fig. 6b). At the end of the main convection event, end of February–early March, the spatially averaged mixed layer abruptly decreased to less than 100 m (Fig. 6b). Finally, during the secondary convection event from 15 to 24 March, it reached almost 800 m. Figure 7b shows the modelled mixed-layer depth (MLD) averaged over the convection periods. It exceeded 1000 m in a central area of the Gulf of Lion, between 41.5 and 42.5^∘ N and 3.5 and 7^∘ E and was smaller than 500 m in the Ligurian Sea. From mid-April to the end of the period, the mixed layer was shallow (depth <50 m) and its depth remained above the nutriclines (Kessouri et al., 2017) and the deep chlorophyll maximum (Kessouri et al., 2018).

Figure 7 Modelled (a) wind velocity (m s⁻¹), (b) mixed-layer depth (m) (dark grey lines represent 500, 1000 and 1500 isocontours and light grey line the contour of the deep-convection area), (c) oxygen saturation anomaly (%) at the surface and (d) air-to-sea oxygen flux ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), averaged over the 2013 deep-convection period (15 January–8 March; 15–24 March).

5 Results

5.1 Seasonal cycle of dissolved oxygen

The good agreement found between model results and in situ measurements (Sect. 3) gave us confidence in the model that we use here to analyse the evolution of the oxygen inventory in the deep-convection area and to quantify the relative contribution of each oxygen flux in its variation: exchanges at the air–sea interface, as well as physical and biogeochemical fluxes in the ocean interior. Based on the evolution of vertical mixing and the phytoplankton growth in the study area, Kessouri et al. (2017) divided the study period into four sub-periods. The first period from September to the end of November, which we will refer to as the autumn period, is characterised by a stratified water column (mean MLD <50 m) and respiration dominating primary production (Kessouri et al., 2018). The second period, from the end of November to the end of March, referred to here as the winter period, is characterised by a sustained vertical mixing (mean MLD >50 m). The third period, called spring, ran from late March to early June. It corresponds to the period of restratification of the water column (Estournel et al., 2016a) and of the peak of the phytoplankton bloom at the sea surface followed by the formation of a deep chlorophyll maximum (Kessouri et al., 2018). The last period, summer, from early June to September, is characterised by a strong stratification (mean MLD <20 m) and the permanent presence of a deep chlorophyll maximum below 40 m depth (Kessouri et al., 2018). In the following, we will analyse the dynamics of dissolved oxygen for these four periods. The time evolution of daily oxygen budget terms is shown in Fig. 6d–f, while the time evolution of cumulative oxygen fluxes and the resulting variation in oxygen inventory for the upper (surface to 150 m) and deeper (150 m to bottom) layers is presented in Fig. 8. The biogeochemical term of the budget is defined as the sum of oxygen production through photosynthesis and of oxygen consumption through respiration by phytoplankton, zooplankton and bacteria and through oxidation of ammonium (nitrification) (see Eq. 1). The physical term is decomposed into two modes of transport: a net lateral transport due to advection (positive values correspond to an input for the deep-convection area) and a net vertical downward transport at the interface between the two layers, at 150 m depth, due to advection and turbulent mixing. Finally, the time evolution of the dissolved oxygen concentration and the oxygen saturation anomaly, ΔO₂, averaged over the convection area is shown in Fig. 9.

Figure 8 Time series from September 2012 to September 2013 of the oxygen inventory (black line) and cumulative air–sea flux (red line), downward transport (dark blue), lateral transport (positive values: input for the convection area, light blue), and biogeochemical flux (green line) in the (a) upper (surface to 150 m) and (b) deeper (150 m to bottom) layer. Unit: mol m⁻².

Figure 9 Time evolution of (a) the dissolved oxygen concentration (µmol kg⁻¹) and (b) the oxygen saturation anomaly (%), with mixed-layer depth (m) indicated by the black line, all horizontally averaged over the deep-convection area.

Autumn

From September to the end of November 2012 (91 d), depth-integrated respiration exceeded depth-integrated primary production in the upper layer (Fig. 6f). The result of biogeochemical processes in the water column was a net consumption of oxygen and a decrease of 1.8 mol m⁻² in oxygen inventory (Fig. 8). Lateral transport was low for autumn and yielded a slight decrease of 0.9 mol m⁻² in oxygen inventory (Fig.  8). The heat loss and vertical mixing caused by the northerly wind gust at the end of October 2012 led to a decrease in surface temperature and consequently to an increase in oxygen solubility (Fig. 6c). In addition, the vertical mixing reached the depth of the oxygen maximum present in the subsurface (Fig. 9). This caused its erosion and an increase in the surface oxygen concentration which was, however, lower than the oxygen solubility (Figs. 9 and 6c). From this event, the NW deep-convection area became undersaturated at the surface (Fig. 9b) and the sea began to absorb atmospheric oxygen (flux towards the ocean of 80 mmol d⁻¹ on 29 October, Fig. 6d). Over the autumnal period, the cumulative air–sea oxygen flux amounted to 0.3 mol m⁻² (Fig. 8a). Globally, the convection area was characterised by a decrease in oxygen inventory of 2.4 mol m⁻², more than two-thirds of which occurred in the upper layer.

Winter

The winter period was defined from late November 2012 to late March 2013 but can be further divided into two sub-periods based on the intensity of the vertical mixing (Kessouri et al., 2017). During the first sub-period, from the end of November to mid-January (44 d), the mixing intensified, but remained moderate: the ML averaged over the deep-convection area remained above the depth of the maximum euphotic layer (150 m; see Sect. 2.1.3) (Fig. 9). The vertical mixing induced a supply of inorganic nutrients in the upper layer that supported primary production. Kessouri et al (2018) identified the beginning of this period as the beginning of a first bloom. From mid-December, the net biogeochemical production of oxygen became positive in the upper layer (Fig. 6f). However, over this sub-period, the influence of biogeochemical processes on the oxygen inventory remained low (−0.3 mol m⁻², Fig. 8). Air-to-sea oxygen flux was marked by several peaks, greater than 250 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Fig. 6d), coinciding with cold gales from the north. Its contribution to the oxygen inventory over this sub-period amounted to 3.6 mol m⁻². Regarding the lateral oxygen export, it contributed to a loss of 1.3 mol m⁻². The sum of the contributions of the different processes in the water column and at the air–sea interface yielded an increase in O₂ inventory of 2.0 mol m⁻² in the water column. A total of 90 % of this increase occurred in the upper layer, from which 0.7 mol m⁻² of O₂ was exported toward the deeper layers.

The second winter sub-period, from mid-January to late March (69 d), corresponds to the period of deep convection. From the middle to the end of January, the surface water masses previously enriched with oxygen, due to primary production and air–sea exchanges, were mixed with the intermediate water masses characterised by a minimum of oxygen (Fig. 9). From the beginning of February, the vertical mixing intensified, causing a net oxygen transport towards deeper layers (depth >800 m, Figs. 6e, 8 and 9). O₂ concentration decreased significantly at the surface and the difference between surface oxygen concentration and oxygen solubility deepened further, with the oxygen saturation anomaly reaching −15 % until the end of the convection period (Fig. 9). Over this sub-period, the whole NW convection area was undersaturated at −10 % to −15 % (Fig. 7c). Strong undersaturation and wind intensity led to very high air–sea fluxes. Several peaks reaching 800 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ are modelled until mid-March (Fig. 6d). The contribution of air–sea fluxes over this period amounted to 18.0 mol m⁻² (Fig. 8a). Over the deep-convection period, air–sea oxygen exchanges are characterised by high spatial variability (Fig. 7d) with a SD of 38 %. The air-to-sea oxygen flux averaged over the deep-convection period varied between 300 and 460 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the heart of the convection area and between 65 and 200 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the Ligurian Sea.

With regard to biogeochemical processes, as shown in previous studies (Auger et al., 2014; Kessouri et al., 2018), zooplankton growth was largely reduced by the deep-convection process due to a dilution-induced decoupling of prey and predators. In the upper layer, oxygen production through primary production exceeded oxygen consumption processes (respiration, nitrification) (Fig. 6f). In parallel, the export of organic matter into the intermediate and deep layers during deep convection (Kessouri et al., 2018) led to an increase in remineralisation processes (Fig. 6f) and consequently a decrease in oxygen inventory in these aphotic layers. The sum of biogeochemical fluxes over the entire water column resulted in a small increase in oxygen inventory of 0.4 mol m⁻², negligible compared to that induced by air–sea fluxes, in consistency with the previous study of Minas and Bonin (1988).

Over this period, the lateral export of dissolved oxygen had high values, reaching 220 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Fig 6e). In the upper layer, the total lateral transport over the period was low (0.5 mol m⁻²), while it is estimated that in the deeper layers 6.7 mol m⁻² was exported horizontally from the convection area between mid-February and the end of the convection period (Fig. 8b). The downward transport at the base of the upper layer showed strong peaks reaching 500 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Fig. 6e), concomitant with the peaks of the air-to-sea fluxes and the deepening of the ML.

The model results indicate that atmospheric oxygen injected at the surface and, to a lesser extent, produced by phytoplankton or horizontally advected in the upper layer, was massively transported to the intermediate and deep layers (20.1 mol m⁻²). It is worth noting that vertical fluxes showed high spatial variability within the convection area. Over this period, the lateral transport from the aphotic layer outside the convection area represents 33 % of the amount of downward transport. Globally, the different contributions led to an increase in the water column oxygen inventory of 12.3 mol m⁻².

Spring (late March to early June, 74 d)

In spring, net biogeochemical production of O₂ remained high in the upper layer until the bloom peak in mid-April; afterwards it decreased but generally remained positive until the end of that period (Fig. 6f). Oxygen consumption through heterotrophic respiration in the deeper layers also remained relatively high. The result of biogeochemical contributions was a small increase of 0.3 mol m⁻² in the O₂ inventory of the water column.

During this period, primary production led to a sharp increase in surface oxygen concentration from 220 to 280 µmol kg⁻¹ at the peak of the phytoplankton bloom (Fig. 6c), and the latter became above saturation in early April, when the convection area became a source of oxygen for the atmosphere (Fig. 6c and d). This oversaturation situation at the surface then persisted until the end of the period. The model simulates significant outgassing during the bloom peak (235 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ on 18 April 2013, Fig. 6d) when the mean saturation anomaly reached a maximum value of 15 % (Figs. 6c and 9). Overall, the convection area released 0.8 mol m⁻² of oxygen to the atmosphere during spring. During this restratification phase, a moderate oxygen export to the deep layers is found (3.2 mol m⁻², Fig. 8). Lateral export to regions surrounding the convection area continued at a high rate with a cumulative value of 5.1 mol m⁻². Finally, over this period, the water column in the convection area was subjected to a 5.7 mol m⁻² decrease in its oxygen inventory, due to the lateral export of oxygen via the spreading of dense waters in the deeper layers and a slight outgassing to the atmosphere.

Summer

During the summer period (87 d), the surface oxygen concentration remained higher than the oxygen solubility (Fig. 6c). A supersaturated situation occurred in the deep chlorophyll maximum zone, due to primary production and a general stratification (Fig. 9). We estimate that the ocean released 1.4 mol O₂ m⁻² to the atmosphere over this period (Fig. 8), mainly during moderate northerly gales. Depth-integrated oxygen-consuming biogeochemical processes exceeded depth-integrated primary production on average over this period. The result of biogeochemical fluxes was responsible for a consumption of 0.8 mol m⁻² of oxygen. In addition, the convection area continued to export oxygen to the adjacent zone (1.3 mol m⁻²), but at a lower rate (15 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) than in the two previous periods (90 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ over the deep-convection period and 69 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in spring). Finally, the oxygen inventory decreased by 3.5 mol m⁻² in the whole water column of the deep-convection area (Fig. 8).

5.2 Annual oxygen budget

Figure 10 illustrates the oxygen budget of the NW Mediterranean convection area over the period September 2012 to September 2013. At the annual scale, the deep-convection area is a net sink of oxygen for the atmosphere, estimated at 20.0 mol O₂ m⁻². A total of 88 % (17.7 mol O₂ m⁻²) of this amount was injected into the ocean interior during the period when the deep-convection process took place.

Figure 10 Schematic showing the terms of the annual oxygen budget (mol O₂ m⁻²) for the north-western Mediterranean deep-convection area over the period from September 2012 to September 2013. F_A–S: air-to-sea flux; F_H: net horizontal transport; F_v 150: net downward transport at the base of the euphotic layer (150 m); F_B: net biogeochemical production; Δ∫O₂: variation in oxygen inventory. Positive fluxes are inputs for the deep-convection zone. The terms of the budget are estimated for the upper euphotic layer (surface to 150 m) and the deeper aphotic layers (150 m to bottom).

The annual net biogeochemical production of oxygen in the euphotic layer (0–150 m) is estimated at 1.6 mol O₂ m⁻². The net annual NCP (net community production, defined as gross primary production minus community respiration in the euphotic zone) is estimated at 3.9 mol O₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, yielding autotrophy in this area. In the deeper layers (150 m–bottom) an oxygen consumption of 3.8 mol O₂ m⁻² was associated with respiration of heterotrophic organisms by 70 % and oxidation of ammonium by 30 %. This led to an annual net biogeochemical consumption of 2.2 mol O₂ m⁻² over the whole water column.

The model indicates that 27.1 mol m⁻² of O₂ was exported from the upper layer to deeper layers. This net transport toward the bottom occurred during 68 % of the events of deep vertical mixing of oxygen-rich surface waters with oxygen-poor underlying waters. Finally, the budget shows that the deep-convection area appears as a net source for dissolved oxygen for the rest of the western Mediterranean Sea with an annual net horizontal transport of 15.0 mol O₂ m⁻². This transport breaks down into an input of 5.3 mol O₂ m⁻² in the upper layer and an export of 20.3 mol O₂ m⁻² in the deeper layer.

At the end of the annual cycle, a negligible decrease (0.3 mol m⁻², i.e. 0.05 %) in the oxygen inventory of the upper euphotic layer is found, while 3.1 mol m⁻² (i.e. 0.66 % of the inventory) was stored in the deeper water masses.

6 Discussion

6.1 Air–sea oxygen flux

Our model results indicate that the NW Mediterranean deep-convection area was a net sink for the atmospheric oxygen at a rate of 20.0 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ between September 2012 and September 2013 and at a rate of 280 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (17.7 mol m⁻² over 63 d) during the 2013 deep-convection period. Inside the area, the annual air–sea flux shows strong spatial heterogeneity, with a range extending from 2.7 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ at the periphery to 36.0 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in the centre. Considering its sea surface area (61 720 km²), the NW deep-convection zone received 1233 Gmol of oxygen from the atmosphere over the period September 2012 to September 2013, including 1090 Gmol during the winter 2013 intense vertical mixing period. We showed that the strong oxygen ingassing was essentially driven by a high undersaturation ($<-\mathrm{10}$ %) and intense northerly winds during the deep-convection period.

Nevertheless, uncertainties in the net uptake rate remain. First, uncertainties are linked to errors in modelled ocean surface variables (dissolved oxygen, temperature and salinity) and wind velocity used for the calculation of the air–sea flux. The comparisons of model results with in situ high-frequency measurements at the surface during the period of maximum flux (deep-convection period) indicate a bias of less than or close to 1 % and a NRMSE smaller than 14 % for the wind velocity, surface temperature, salinity and oxygen concentration (Sect. 3.1). A second source of uncertainty is linked to the parameterisation chosen for the calculation of the gas transfer velocity. In the standard run, we used the cubic dependence with wind speed parameterisation proposed by Wanninkhof and McGillis (1999). Sensitivity analyses were performed using eight other parameterisations for the calculation of air–sea flux (Wanninkhof, 1992, 2014; Woolf, 1997; Nightingale et al., 2000; Wanninkhof et al., 2009; Stanley et al., 2009; Liang et al., 2013; Bushinsky and Emerson, 2018; see Sect. 2.1.2). Estimates of annual air–sea flux, as well as flux and amount of atmospheric oxygen captured by the study area during the deep-convection period, calculated with all these parameterisations are gathered in Table 2. All estimates show a net sink for atmospheric oxygen for the study area. They range from 14.2 to 21.5 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ at the annual scale, with a mean value of 17.7±2.8 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, and from 188 to 285 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, with a mean value of 242±38 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, during the deep convection. Both estimates in the standard run are in the upper range of all estimates. Considering all estimates, we determine an uncertainty (SD) of 15 %–16 % for the annual and convection period air–sea flux. This uncertainty, associated with the parameterisation of the gas transfer velocity, propagates to the estimates of vertical and lateral transport of oxygen in the ocean interior. Depending on the gas transfer parameterisation used, at the annual scale, downward export below the euphotic zone ranges from 22.0 to 27.9 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (mean value: 25.1±2.2 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), and lateral transport ranges from 4.8 to 6.1 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (mean value: 5.6±0.5 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) in the euphotic layer and from −17.1 to −20.6 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (mean value: $-\mathrm{19.0}\pm \mathrm{1.3}$ $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) in the aphotic layer. During the deep-convection event, downward export below the euphotic zone ranges from 223 to 302 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (mean value: 265±30 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). The uncertainty on the transport terms of the annual budget thus remains smaller than 12 %. The values of the NRMSE between cruise observations and modelled dissolved oxygen from sensitivity tests are found very close to the NRMSE obtained for the standard run. Slightly smaller NRMSE values (≈10 %) are found only for the winter DEWEX-Leg1 cruise period using the parameterisations of Wanninkhof and McGillis (1999), Woolf (1997), Stanley et al. (2009), and Liang et al. (2013), which give a higher oxygen transfer than the other parameterisations. Previous studies on oxygen air–sea flux in deep-convection zones recommended the use of parameterisations with high transfer during periods of strong wind and convection (Copin-Montégut and Bégovic, 2002; Körtzinger et al., 2008b; Koelling et al., 2017; Atamanchuk et al., 2020). Atamanchuk et al. (2020), comparing flux estimates based on several parameterisations, found that these flux estimates may vary by an order of magnitude and warned of the possibility of a strong underestimation of air–sea oxygen flux in biogeochemical models that do not include bubble-mediated terms. In our study, the range of estimates obtained with both types of parameterisations, those that are only diffusive and those that include bubble-mediated terms, is similar. Although the parameterisation of Wanninkhof and McGillis (1999) used in our standard run does not include an explicit bubble-mediated transfer term, it provides estimates of air–sea fluxes close to those obtained with the bubble-flux-inclusive one of Stanley et al (2009), preferred by Atamanchuk et al. (2020) in their Labrador Sea study. The strong undersaturation obtained in the north-western Mediterranean during the convection period, between −10 % and −20 %, may explain a greater contribution of the diffusive flux compared to the air injection by bubbles. Moreover, winter conditions are less extreme than in the Labrador Sea where strong wind speeds of more than 13.8 m s⁻¹ were encountered for at least 40 d. In the NW Mediterranean Sea and in winter 2012/13, only 13 % of the convection area was characterised by a number of days with wind speeds >13.8 m s⁻¹ varying between 30 and 35 d. An experimental study of flux measurements in this region over a whole year would allow a better assessment of the contribution of air injection in the total air–sea flux and hence of the different parameterisations of gas transfer.

Table 2 Estimates of air-to-sea oxygen flux (F_A–S) for the period September 2012–September 2013 and during the deep convection (15 January–8 March, 15–24 March), using different parameterisations of gas transfer velocity.

Previous studies based on in situ observations have proposed estimates for the air–sea oxygen flux in the study area. Our modelled seasonal cycle of air–sea oxygen flux agrees with the results of Copin-Montégut and Bégovic (2002) and Coppola et al. (2018) in the Ligurian Sea, at the DYFAMED site, who observed an annual cycle with a net ingassing from December to March and net outgassing from April to November. In the Ligurian Sea the deep-convection process does not occur each winter. When occurring, it is generally shorter and shallower than in the centre of the Gulf of Lion, the core of dense water formation. Coppola et al. (2018) using temperature, salinity and oxygen monthly profiles and the gas transfer parameterisation of Ho et al. (2006) estimated for the period 1994–2014 a monthly air-to-sea flux varying from −15.1 to 14.8 mol O₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, with an annual mean value of −2.6 mol O₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. Over this 20-year period, the authors identified one winter, winter 2005/2006, with intense vertical mixing reaching the deep layers, and four winters (1999, 2000, 2006 and 2013) with moderate vertical mixing reaching intermediate depths. From the difference in O₂ inventory between December 2005 and April 2006, they deduced that 24 mol m⁻² of atmospheric O₂ was injected between 350 and 2000 m at a rate of 300 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. At the same site, Copin-Montégut and Bégovic (2002) estimated an air–sea ingassing of 5 and 2.6 mol O₂ m⁻², respectively, for the moderate cold winters 1999 (for 26 d, rate of 190 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and 2000 (for 23 d, rate of 110 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), respectively, using in situ surface measurements of oxygen in winter and the formulation of gas transfer velocity from Wanninkhof and McGillis (1999). Those estimates were twice as small as their observation of variation in the oxygen content in the first 600 m of depth, namely 11 and 15 mol m⁻² (in 1 month) for winters 1999 and 2000, respectively. Those authors suggested an underestimation in their estimates due to low-frequency measurements and an underestimation of the gas transfer coefficient. At the same location, for the study period, we found a net oxygen ingassing of 9.2 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ at the annual scale and 135 mmol O₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ during the period of deep-convection events (63 d) (Fig. 7d). Thus our calculation of atmospheric oxygen uptake in the Ligurian Sea is close to the ones of Copin-Montégut and Bégovic (2002) for moderate convective winters. Our estimates in the centre of the Gulf of Lion, where convection reached the deep waters, with values of 20–28 mol O₂ m⁻² during deep convection are also close to the estimate by Coppola et al. (2018) for the intense vertical mixing winter 2005/2006 in the Ligurian Sea.

Finally, our model calculation of air–sea oxygen flux for the NW Mediterranean is in the same range found for other worldwide deep-convection areas. At the centre of the Labrador Sea, Körtzinger et al. (2008b) found an annual air–sea ingassing of 10.0±3.1 mol O₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ over the period 2004/2005, using in situ observations at the K1 mooring site and the Wanninkhof (1992) parameterisation. By quantifying the relative contribution of biogeochemical and lateral fluxes, Koelling et al. (2017) estimated an oxygen ingassing of 29.1±3.8 mol m⁻² over winter 2014/2015 at the same mooring site K1. Wolf et al (2018) derived from Argo float observations in the Labrador Sea mean air–sea fluxes with a large range of values varying from 5.7 to 22.8 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ using various parameterisations. Using parameterisations including bubble-mediated fluxes (Liang et al., 2013; Yang et al., 2017), their estimates of atmospheric oxygen uptake ranged from 21.6 to 36.6 mol m⁻² for the convective winter 2003/2004. Based on measurements of oxygen from a moored profiler and Argo floats and on the Stanley et al. (2009) parameterisation, Atamanchuk et al. (2020) estimated an annual air–sea flux of oxygen of 19.3±3.4 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for the year 2016/2017. For the Irminger Sea, Maze et al. (2012) estimated an abiotic air–sea oxygen flux of 13±3 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for the years 2002, 2004 and 2006, using an optimisation method and observations from three surveys and Word Ocean Atlas 2009.

6.2 The role of the NW deep-convection area in the ventilation of the western Mediterranean Sea

Open-sea convection and shelf dense water cascading (Canals et al., 2006; Ulses et al., 2008b) in the NW Mediterranean are the main mechanisms for the ventilation of the entire western Mediterranean Sea. Over the past decades, several observational studies reported increases in O₂ concentration in deep water masses at several sites in the western Mediterranean where deep convection did not occur and where winter vertical mixing was limited to surface or intermediate levels. Coppola et al. (2018) associated the high concentrations of O₂ observed in deep layers of the Ligurian Sea in 1994 and 2005, when convection was limited to intermediate waters, with the arrival of deep water formed in the open sea of the Gulf of Lion or formed on the Gulf of Lion shelf and cascading down to the deep basin. Using measurements from five cruises, Schroeder et al. (2008a) documented an abrupt increase in heat, salt and O₂ inventory of deep waters in an extensive area of the western Mediterranean, occurring in 2005 and 2006. The authors attributed these changes, referred to as the western Mediterranean Transient (hereafter WMT; Zunino et al., 2012), to the propagation of the new dense waters formed in the NW deep-convection area during the winters 2004/2005 and 2005/2006, when severe weather conditions caused intense dense water formation (Lopez-Jurado et al., 2005; Schroeder et al., 2006). The study of Schroeder et al. (2008a) showed the presence of these new O₂-rich deep waters in the Balearic Sea, the Ligurian Sea and in the Algerian sub-basin in June 2005 and their propagation to the whole Algerian sub-basin and the west of the Alboran Sea in October 2006. New oxygenated waters were also observed in the entire deep layers of the Algerian sub-basin in 2011 (Schneider et al., 2014; Stöven and Tanhua, 2015), 2014 (Keraghel et al., 2020), 2016 and 2018 (Li and Tanhua, 2020). Moreover, the results of Li and Tanhua (2020) showed a ventilation of the deep waters of the Tyrrhenian Sea through an overflow of well-oxygenated water masses from the Algerian basin into the deep layer, between 2011 and 2016.

Somot et al. (2016) found that winter 2012/13 is one of the five winters over the 33-year period 1980–2013 showing high dense water formation rates (above 0.6 Sv), using the CNRM-RCSM4 model. According to their estimates, the cumulative volume of dense water formed over the winters 2011/2012 (0.45 Sv) and 2012/13 (0.7 Sv), amounting to 1.15 Sv, may be close to the volume of dense water formed in winter 2004/2005 of 1.2 Sv. As a result, these successive 2012 and 2013 deep-convection events could have been responsible for a similar ventilation as the one observed after the event of 2005 (Schroeder et al., 2008b; Schneider et al., 2004; Stoven and Tanhua, 2015), assuming similar air–sea exchanges.

Our modelling study indicates that, over the period September 2012 to September 2013, the upper layer of the NW deep-convection area captured 5.3 mol O₂ m⁻² from the surrounding regions, in addition to the 20.0 mol m⁻² of oxygen from the atmosphere, while the deeper layers released 20.3 mol m⁻² toward the adjacent seas (Sect. 5.2). Considering the deep-convection surface area of 61 720 km², the lateral transport led to a gain in the upper layer of 330 Gmol yr⁻¹ in the convection area and a loss of 1250 Gmol yr⁻¹ towards the adjacent deep areas. As a result, the NW convection area appears as a source of 920 Gmol yr⁻¹ of oxygen for the rest of the western basin for the period 2012/13. Lateral O₂ inputs in the upper layer occurred mainly from February to September with two peak periods, in early March, a calm period between two convective events, and in April, during restratification. These imports were mainly related to eddies produced by the baroclinic instability that was triggered at the periphery of the convection zone when strong wind ceased (Killworth, 1976; Testor et al., 2018). These inputs from the peripheral zone contributed to the vertical export of oxygen to the aphotic layer. First in the short term, the oxygen imported between two convection events was exported at depth by the following events. At longer timescales (April–September), the convection area was also fed by the peripheral zones and in turn produced a vertical export to the aphotic layer. These exchanges were of lower intensity and concerned shallower layers but are not negligible when integrated over the year. Regarding the lateral transport in the deeper layer, our model outputs show that the O₂-rich dense waters formed in the NW deep-convection area propagated towards the Balearic Sea, first at intermediate depths (150–800 m) from the beginning of the winter mixing period and then through deep layers (800 m to bottom) from mid-February (not shown). These water masses then mostly flowed towards the south of the western basin, while a smaller part was advected back in the convection area through mesoscale circulations counteracting the effect of the intrusions of low-oxygen LIW during the restratification period, in increasing the oxygen inventory of intermediate waters (not shown). A preferential pathway to the south of the basin was the one along the eastern coast of Minorca in the Algerian sub-basin, in agreement with previous observational and modelling studies that examined the spreading of waters formed in winter in the NW region (Pinot and Ganachaud, 1999; Schroeder et al 2008b; Beuvier et al., 2012). Our simulated circulation of oxygen in the western basin is also consistent with the study of Piñeiro et al. (2019), who reported the arrival of new dense water masses formed in the deep-convection area east of Minorca over the 2011–2013 period using temperature and salinity observations at the hydrographic stations RADMED. In our model outputs, the offshore Balearic Sea (bathymetry >1000 m, surface area: 19 700 km²) and Algerian sub-basin (bathymetry >1000 m, surface area: 171 610 km²) experienced an increase in their oxygen inventory, during and after the NW deep-convection events, receiving oxygen through lateral transport (271 and 1276 Gmol, respectively) while the amounts of oxygen captured at the air–sea interface during the period of intense vertical mixing were smaller in those areas than in the NW deep-convection area by a factor of 10 and 3, respectively (i.e. 104 and 385 Gmol versus 1090 Gmol for the NW deep-convection area). This suggests that an important part of the oxygen absorbed at the air–sea interface of the NW deep-convection area, exported first vertically towards its deeper layers and then horizontally towards the adjacent regions, was stored, at least temporarily, in the Algerian sub-basin.

Finally, our results demonstrate that the total oxygen supply by air–sea exchanges in the NW deep-convection region for the period 2012/13 (1233 Gmol yr⁻¹), which was then mainly released to adjacent seas in the aphotic layer, constitutes a major source of oxygen at the scale of the whole Mediterranean Sea. Indeed this supply is close to the biogeochemical oxygen consumption within the Mediterranean Sea estimated at 1545 Gmol yr⁻¹ by Huertas et al. (2009) using in situ measurements at the Strait of Gibraltar over the period 2005–2007.

The present study of the period 2012/13 constitutes a first step in our analysis and quantification of the oxygen budget for the western Mediterranean Sea. Previous observational studies (Coppola et al., 2018; Mavropoulou et al., 2020) over periods of 20 years or more showed that the mean oxygen concentration in the western Mediterranean and in particular in the NW deep-convection area is subjected to a strong interannual variability, mainly in response to the variability of deep convection, the latter being influenced by transient changes such as the WMT event. A deeper analysis of the physical processes involved in the vertical and horizontal transport in the convection zone as well as of the spreading of the oxygen-enriched dense waters, formed in the NW deep-convection area, in the western basin and toward the Atlantic Ocean through the Gibraltar Strait will be conducted in further studies using a numerical simulation with extended domain and period.

6.3 Net community production

Our budget calculation shows that in this region characterised by intense vertical mixing the biogeochemical terms remained very low compared to the air–sea oxygen flux over the period 2012/13. Our modelling results indicate that the net biogeochemical production of oxygen in the euphotic layer was positive from mid-December to the end of July and negative the rest of the year. It was maximum during the spring bloom from mid-March to mid-April. We estimate a net annual NCP (in the upper layer) of 46.8 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (3.9 molO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; see Eq. 1). This indicates a net autotrophy for the euphotic layer of the NW Mediterranean deep-convection area over the period September 2012–September 2013. If consumption of oxygen through nitrification is considered, net biogeochemical production amounted to 1.6 mol O₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. It is worth noting that nitrification, discussed in Kessouri et al. (2017), who estimated a nitrogen budget using the same coupled model, accounts only for 7 % of the total oxygen consumption but for 60 % of the NCP, suggesting that this process should be considered when estimating NCP from oxygen concentration.

Our value of NCP is higher than the net downward export of organic carbon at the base of the euphotic layer estimated at 35 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (25 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for particulate organic carbon and 10 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for dissolved organic carbon) by Kessouri et al. (2018) over the same period and using the same coupled model. Also using the same coupled model, Kessouri et al. (2017), who analysed the nitrogen cycle over the study period, obtained a new primary production varying from 65 to 77 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in the deep-convection zone. By analysing the carbon, nitrogen and oxygen cycles, the NW deep-convection region is always found to be an autotrophic ecosystem. The discrepancies in magnitude obtained depending on the element considered reflect different dynamics for these elements in the euphotic layer, possibly due to variable ${\mathrm{O}}_{\mathrm{2}}:\mathrm{C}:\mathrm{N}$ ratios in biogeochemical production and consumption processes as shown by Copin-Montégut (2000) in the Ligurian Sea using high-frequency measurements.

Our estimate of NCP is smaller than the estimate of 85.2 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in the Ligurian Sea at the DYFAMED site over the period 1994–2014 by Coppola et al. (2018) using monthly observations. It is close to the estimate by Ulses et al. (2016) of 42.8 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ over the period 2003–2008 using the same numerical model and considering an area extending to the whole offshore NW Mediterranean and a 100 m thick upper layer. It is also similar to the previous estimates of new primary production by Severin et al. (2014) varying from 46 to 63 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ over the period February–March 2011, based on in situ nutrient concentrations.

NCP is often used to estimate the strength of the biological pump and the potential capacity of a system to capture atmospheric CO₂. Although the NW deep-convection pelagic ecosystem appears as a net annual sink for atmospheric CO₂ from our modelling results and previous studies (Coppola et al., 2018; Ulses et al., 2016), the role of this region in terms of carbon sequestration remains highly uncertain. Deep convection generates a strong downward transport of organic carbon below the euphotic layer (Ulses et al., 2016; Kessouri et al., 2018). A large amount of organic carbon transferred below the euphotic zone is then consumed and demineralised after deep convection (Santinelli et al., 2010), leading to an increase in CO₂ inventory into the deeper reservoir that could be raised back in the euphotic layer during the following deep-convection events as shown in the Atlantic Ocean by Körtzinger et al. (2008a) and in the Pacific Ocean by Palevsky et al. (2016). Episodes of oversaturation of sea surface pCO₂ related to atmospheric pCO₂ were reported during short wind gusts and intense vertical mixing events in the Ligurian Sea (Copin-Montégut et al., 2004; Merlivat et al., 2018) and in the central Gulf of Lion open sea (Touratier et al., 2016). The authors explained those oversaturation episodes by the increase in CO₂ concentration at the ocean surface induced by the mixing of surface CO₂-poorer waters with deep CO₂-rich waters. Those punctual observations suggested short releases of CO₂ by the ocean induced by deep convection. On the other hand, using a 0.5^∘ resolution array of coupled 1D hydrodynamic–biogeochemical models of the upper layer in the Mediterranean Sea over the period 1998–2004, D'Ortenzio et al. (2008) estimated that the NW region is a sink for atmospheric CO₂ in winter and at the annual scale (between 12 and 24 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) and found that biogeochemical processes dominate air–sea exchanges and mixing processes in this region most of the year. In another region of deep convection, the Labrador Sea, DeGrandpre et al. (2006) and Körtzinger et al. (2008b) also found an annual uptake of atmospheric CO₂ that amounted to 55.2 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for the period 2000–2001 and 32.4±9.6 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for the period 2004–2005, using mooring observations (and also a 1D biogeochemical model for DeGrandpre et al., 2006).

In the study area, our results show that lateral transport dominated the budget of oxygen during the restratification period when deep dense waters spread in the western basin and LIW reintegrated the deep-convection zone. This suggests that coupled 3D biogeochemical–physical models, including a carbonate system module, could be useful tools to complete the previous 1D studies of dissolved inorganic carbon dynamics and integrate on an annual scale the exchanges at the air–sea interface by taking into account lateral transport and mesoscale structures influencing the spreading of water masses and their compounds during convection and restratification phases and impacting the budgets.

7 Conclusions

Our study is the first attempt to describe the seasonal cycle of dissolved oxygen and to estimate the oxygen budget over the whole NW Mediterranean deep-convection area, using a high-resolution coupled 3D physical–biogeochemical model. The assessment of the model results using in situ measurements from DEWEX and MOOSE-GE cruises and from Argo-O₂ floats shows the ability of the model to capture the main spatial and temporal variability of dissolved oxygen observed. From our modelling results, the following conclusions can be drawn for the period 2012/13.

• The seasonal cycle of surface dissolved oxygen in this area exhibited a winter period with strong undersaturation due to a decrease in temperature and surface dissolved oxygen concentration, induced by strong heat loss and vertical mixing of surface O₂-rich waters with the underlying low-O₂ waters. The undersaturation averaged over the whole area indeed reached −15 % during the deep-convection event. During the stratified period, an oversaturation situation occurred with a maximum surface value of 15 % during the peak of the spring bloom.

• The NW Mediterranean deep-convection area acted as a large sink for atmospheric oxygen. We estimate that the area captured 20 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ of atmospheric oxygen at the annual timescale. An uptake of 18 mol m⁻² of atmospheric oxygen, which equals 88 % of the annual uptake, took place during the deep-convection period. This uptake is characterised by high spatial variability, with a SD of 38 % in this area including the open sea of the Gulf of Lion and Ligurian Sea. The magnitude of the uptake is maximum inside a central zone of the Gulf of Lion where the average over the deep-convection period (63 d) reached a rate ranging between 300 and 460 mmol O₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$.

• The NW Mediterranean deep-convection area acts as a conveyor of atmospheric oxygen, as well as of oxygen produced in the upper layer, both locally and in the surrounding areas, towards the intermediate and deep layers of the western Mediterranean Sea. Based on the rate of dense water formation (Somot et al., 2016), the ventilation due to deep convection in the NW area in 2012/13 may represent half of the ventilation observed in 2004/2005 by Schroeder et al. (2008a). The magnitude of atmospheric O₂ uptake and lateral transport to the adjacent regions in the aphotic layers in 2012/13 is close to the magnitude of the oxygen consumption of the whole Mediterranean Sea estimated by Huertas et al. (2009).

• Sensitivity tests to the parameterisation of the gas transfer velocity yield an estimate of the budget term (air–sea exchanges and transport terms) uncertainty of 12 %–16 %.

• As expected for this very energetic region, the annual budget of oxygen is clearly dominated by air–sea exchanges and physical transport over convective years such as 2012/13. The net biogeochemical production in the euphotic zone is estimated to account for 10 %, i.e. 2 mol O₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, of the net atmospheric oxygen uptake. In deeper depths, heterotrophic organisms' respiration and nitrification resulted in an oxygen consumption of 4 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$.

• The NW Mediterranean deep-convection area is found to be an autotrophic ecosystem with an annual NCP (in the 150 m upper layer) estimated at 47 g C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$.

The high interannual variability of deep convection in the NW Mediterranean (Houpert et al., 2016; Somot et al., 2016) suggests a high variability of the oxygen budget. Further modelling at pluri-annual and multi-decadal scales is thus needed to investigate the interannual variability of the annual budget over the whole western basin, as well as the evolution of this budget under climate warming, the effects of which could have been masked for the time being by the significant impacts of climatic transient shifts such as WMT according to Mavropoulou et al. (2020).