Net community production in the northwestern Mediterranean Sea from glider and buoy measurements

. The Mediterranean Sea comprises just 0.8 % of the global oceanic surface, yet considering its size, it is regarded as a disproportionately large sink for anthropogenic carbon due to its physical and biogeochemical characteristics. An underwater glider mission was carried out in March–April 2016 close to the BOUSSOLE and DyFAMed time series moorings in the northwestern Mediterranean Sea. The glider deployment served as a test of a prototype ion-sensitive ﬁeld-effect transistor pH sensor. Dissolved oxygen (O 2 ) concentrations and optical backscatter were also observed by the glider and increased between 5 19 March and 1 April, along with pH. These changes indicated the start of a phytoplankton spring bloom, following a period of intense mixing. Concurrent measurements of CO 2 fugacity and O 2 concentrations at the BOUSSOLE mooring buoy showed ﬂuctuations, in qualitative agreement with the pattern of glider measurements. Mean net community production rates ( N ) were estimated from glider and buoy measurements of dissolved O 2 and inorganic carbon (DIC) concentrations

the interior ocean is key for quantifying the effects of a future warmer climate (Bauer et al., 2013). Whether a location is predominantly autotrophic (dominated by photosynthesis) or heterotrophic (dominated by respiration) determines the sign of net community production, which is defined as gross primary production (by phytoplankton) minus total respiration (by phytoplankton, zooplankton and bacteria) (Alkire et al., 2012).
The surface of the Mediterranean Sea (2.5 × 10 6 km 2 ) represents just 0.8 % of oceans globally, but relative to its size, it is regarded as an important sink for anthropogenic carbon dioxide emissions due to higher levels of anthropogenic carbon than in the Atlantic and Pacific oceans (Lee et al., 2011;Schneider et al., 2010). This is due to a low Revelle factor (Broecker et al., 1979) related to relatively warm, salty and high-alkalinity waters, encouraging a net flux of carbon dioxide from the atmosphere to the ocean. Carbon dioxide dissolved in water (CO 2 (aq) and H 2 CO 3 ) dissociates to bicarbonate (HCO − 3 ) and carbonate (CO 2− 3 ), releasing H + ions (Zeebe and Wolf-Gladrow, 2001). CO 2 (aq), H 2 CO 3 , HCO − 3 and CO 2− 3 make up total dissolved inorganic carbon (DIC), with HCO − 3 accounting for 90 % of DIC. Carbon dioxide absorbed by the ocean is thought to reach the interior via deep water formation and biological processes (Álvarez et al., 2014;Arrigo et al., 2008).
The northwestern Mediterranean Sea displays strong seasonal variability. At the surface, temperatures remain at around 13 to 14 • C during winter, increasing to maxima of 26 to 28 • C during summer. Wind-driven vertical mixing occurs during autumn and winter, whilst surface stratification is common during summer as a result of solar heating (Copin-Montégut et al., 2004;D'Ortenzio et al., 2008). Vertical mixing can transport nutrients from greater depths to oligotrophic surface waters (Marty and Chiavérini, 2002;de Fommervault et al., 2015). A combination of nutrient availability and increased stability driven by surface warming of 0.2 • C can trigger phytoplankton blooms (Yao et al., 2016;Copin-Montégut et al., 2004). The onset of the spring bloom in the northwestern Mediterranean Sea varies from March to April, as observed between 2013 and 2015 at the BOUS-SOLE buoy (Bouée pour l'acquisition de Séries Optiques à Long Terme, http://www.obs-vlfr.fr/Boussole/, last access: 20 August 2022) close to the DyFAMed (Dynamique des Flux Atmospheriques en Mediterraneé) site in the Ligurian-Provençal basin (Merlivat et al., 2018). Significant increases in particulate and dissolved organic carbon (POC and DOC) concentrations have been observed during bloom events (e.g. Carlson et al., 1998). POC and DOC export from the ocean surface constitutes the "biological carbon pump" (Carlson et al., 1998;Van Der Loeff et al., 1997;Alkire et al., 2014).
Quantifying these processes in detail requires sufficient data coverage in space and time. Few DIC time series have been maintained continuously, among them the DyFAMed mooring, which is complemented by monthly ship hydro-casts (Copin-Montégut et al., 2004;Antoine et al., 2008;Taillandier et al., 2012). The DyFAMed site is considered an open-ocean location as it is roughly 52 km from the coast in > 2000 m deep water. The mooring is useful for studying processes occurring at specific depth levels at one location (Merlivat et al., 2018;Copin-Montégut et al., 2004;Hood and Merlivat, 2001), but a lack of vertical and horizontal spatial information is a limiting factor when quantifying mass and energy budgets. Autonomous underwater gliders have been used to survey the northwestern Mediterranean Sea since 2005 (Niewiadomska et al., 2008;Cyr et al., 2017). They are useful platforms for a range of physical and biogeochemical sensors and can operate autonomously for many months in up to 1000 m deep water using battery power (Eriksen et al., 2001;Piterbarg et al., 2014;Queste et al., 2012). The deployment of autonomous platforms, such as underwater gliders, complements fixed-depth time series by enabling observations of biogeochemical and physical horizontal and vertical gradients. This paper aims to estimate net community production (N) at the DyFAMed site using in situ continuous measurements from a mooring and an underwater glider deployed in March-April 2016. The additional glider data help overcome limitations of the single-depth mooring data in terms of both vertical data coverage and the contribution of horizontal advection. Furthermore, the glider mission served as a test for a prototype ion-sensitive field-effect transistor (ISFET) pH sensor, which complemented the standard temperature, salinity and c(O 2 ) sensors, to provide both O 2and DIC-based net community production.
2 Data collection and quality

Ship measurements
Ship CTD and water sample profiles were collected by RV Téthys II on 7 March and 16 April 2016 at the DyFAMed site ( Fig. 1). Ship hydrocast profiles of temperature and salinity (Sea-Bird Scientific SBE 9 CTD), and c(O 2 ) (Sea-Bird Scientific SBE 43 sensor), were supplemented by discrete Niskin bottle samples for c(O 2 ), c(DIC), total alkalinity (A T ) and concentrations of total oxidised nitrogen (NO − x , i.e. the total of NO − 3 and NO − 2 ), silicate (Si(OH) 4 ) and phosphate (PO 3− 4 ); see Appendix A for details. Implausible outliers in the CTD profiles (< 1 % of values) were flagged and excluded from further analysis.

BOUSSOLE and meteorological buoy measurements
At BOUSSOLE, CO 2 fugacity f (CO 2 ) (Wanninkhof and Thoning, 1993) at 10 m depth was measured spectrophotometrically via a CARIOCA sensor using the thymol blue pH indicator. Inside an exchanger cell, dissolved CO 2 equilibrates with the pH indicator across a silicon membrane. The  change in the optical absorbance of the pH indicator is converted to hourly readings of f (CO 2 ) (Hood and Merlivat, 2001), with an accuracy of 3 µatm (Copin-Montégut et al., 2004). The CARIOCA sensor is replaced roughly every 6 months with a newly calibrated instrument (Merlivat et al., 2018). To remove temperature effects, we show f 13 (CO 2 ), normalised to a temperature of 13 • C (Takahashi et al., 1993). Temperature and salinity were measured hourly using two Sea-Bird Scientific SBE 37-SM MicroCAT sensors at 3 and 10 m depth. The 10 m MicroCAT failed on 15 March 2016. Up to that point, the mean temperature difference between 10 and 3 m was (−0.01 ± 0.03) • C. The corresponding mean practical salinity S(PSS-78) difference was −0.069 ± 0.004. The 10 m MicroCAT had the same salinity offset to the ship CTD cast on 7 March. The 3 m MicroCAT agreed to within 0.007 • C with temperature and to within 0.006 with salinity of the ship CTD cast on 7 March. A total of 3 d before the ship CTD cast on 16 April, the 3 m MicroCAT also failed. However, its final salinity reading was within 0.01 of the ship value. We therefore used the 3 m MicroCAT temperature and salinity readings for calculating buoy-related results.
c(O 2 ) was measured at BOUSSOLE using Aanderaa oxygen optodes at 3 and 10 m. Only nighttime optode measurements were used because the daytime readings were affected by light-induced spikes (Binetti et al., 2020). Winkler titration samples taken during the ship visits on 7 March and 16 April were used to calibrate optode c(O 2 ) (Coppola et al., 2018).
The meteorological buoy "Côte D'Azur" maintained by Météo-France is located close to the BOUSSOLE buoy ( Fig. 1). This meteorological buoy measures wind speed at 3.8 m height above sea level (which is extrapolated to 10 m by adding 10 % of the value), wind direction, air temperature, atmospheric pressure, relative humidity, precipitation and sea surface salinity. Data are archived in the database "SErvice de DOnnées de l'OMP (SEDOO) Mistrals" (http: //mistrals.sedoo.fr/Data-Download, last access: 20 August 2022). Wind speed and sea level pressure were used in this study.

Glider measurements
An iRobot Seaglider model 1KA (sg537; named "Fin") with an ogive fairing was deployed close to the DyFAMed mooring site. A total of 147 dives (294 profiles) were completed by the glider between 7 March and 5 April 2016, covering a diamond-shaped pattern at 43.22-43.50 • N, 7.64-8.00 • E (Fig. 1). The glider completed seven circumnavigations of the survey pattern, each in approximately 4 d, with each of the four sides of the pattern taking approximately 1 d to complete. The glider was equipped with a non-pumped SBE model CTD sensor, an Aanderaa model 4330F oxygen optode sensor, a WET Labs Eco Puck sensor measuring optical backscatter at two wavelengths (470 and 700 nm) and two paired prototype ISFET pH and p(CO 2 ) sensors (Shitashima et al., 2013;Hemming et al., 2017).
Conductivity, temperature and pressure were used to derive potential temperature (θ ) and practical salinity (S). Glider c(O 2 ) was calibrated against ship-based c(O 2 ) profiles ( Fig. 2a, Appendix B1). Outliers (< 1 % of values) outside specified standard deviation ranges (3.5 standard deviations for depths > 400 m; 10 standard deviations for depths < 400 m) were flagged and discarded from further analysis.
The deployment offered a second opportunity to trial prototype ISFET sensors, previously tested on a glider during the REP14-MED experiment in the Sardinian Sea Onken et al., 2018). ISFET pH was corrected ( Fig. 2b) for drift and pressure effects, similar to the steps undertaken by Hemming et al. (2017). The drift-and pressure-corrected ISFET pH was calibrated by linear regression against ship pH on 7 March and 16 April 2016 (Appendix B2).

Development of a spring bloom
The spring bloom at the DyFAMed site is characterised by a decrease in surface f 13 (CO 2 ) due to photosynthesis. Its onset varies from year to year (e.g. April 2013, March 2014) and usually occurs after a period of deep mixing (Merlivat et al., 2018). The glider deployment in March-April 2016 provided horizontal and vertical context to the fixed-depth near-surface mooring sensors when the spring bloom was expected to occur.
The higher-f 13 (CO 2 ), highersalinity waters were likely the results of wind-induced deep mixing and increased convection (Fig. 4a). These waters originate from 50-150 m depth where high-salinity (Fig. 5b), high-f 13 (CO 2 ), low-c(O 2 ) (Fig. 5c) Levantine Intermediate Water (LIW) exists (Knoll et al., 2017). SST increased by up to 0.6 • C between 19 and 20 March and continued to increase intermittently to a maximum of 14.3 • C on 5 April (total increase of 1 • C over the deployment period).
Whilst SST increased between 19 March and 1 April, f 13 (CO 2 ) decreased by 85 µatm. This was the start of the spring bloom, associated with increased surface c(O 2 ), pH T and optical backscatter (Figs. 4b,5c,d,e). O 2 is produced by photosynthesis. The pH T increase reflects biological CO 2 uptake altering the carbonate equilibria (Cornwall et al., 2013;Copin-Montégut and Bégovic, 2002), and increases in temperature. Optical backscatter relates to predominantly POC (Stramski et al., 2004), with minor contributions from mineral particles and gas bubbles at this site. Therefore, all of these observations suggested an increase in biological production.
A clear relationship between potential temperature (θ ), c(O 2 ) and pH T existed within the euphotic layer (mean depth: z eu = 46 ± 7 m) (Figs. 3, 5a, c, e), with higher potential temperature associated with higher pH T and c(O 2 ). The stronger surface stratification resulting from calmer meteorological conditions (e.g. reduced wind speed after 19 March, Fig. 4a) enhanced light supply and stability and caused productivity to increase (Sverdrup, 1953;Pingree et al., 1977). As a result of increased photosynthesis, surface waters became O 2 -supersaturated by 27 March (Figs. 4b, 5c).
Nutrient concentrations obtained by ship hydrocasts for depths of 10-90 m on 7 March were generally highest at 90 m where there is increased remineralisation and lowest at the surface where there is increased usage by phytoplankton (Fig. 6a). Clear differences can be seen in nutrient samples obtained before and after the spring bloom (Fig. 6b). NO − x concentrations decreased by 7 mmol m −3 , Si(OH) 4 by 2.5 mmol m −3 and PO 3− 4 by 0.22 mmol m −3 at 30 m depth, associated with biological production during the spring bloom period. A deep chlorophyll a maximum is often found around this depth (Estrada, 1996) later in the season but was not apparent in our data. At 70 and 90 m depth, nutrient concentrations did not vary much between March and April, while there was a significant decrease in c(O 2 ) and a small change in potential temperature (Fig. 6c). This corresponded to increased stratification with a stronger potential temperature gradient within the top 70 m.

Estimating net community production from O 2 and DIC mass budgets
Net community production N is estimated based on mass budgets using the method described by Alkire et al. (2014): This equation applies to N (O 2 ). For DIC, the N (O 2 ) term is replaced with −N (DIC) because during photosynthesis,  DIC is consumed and O 2 is produced, and vice versa for respiration.
The following diagnostic equations are therefore used to calculate N based on O 2 and DIC mass budgets.
I t is the integrated column inventory change over time, F adv is horizontal advection, F ase is the air-sea flux (positive for the direction ocean to atmosphere), F e is the entrainment flux due to mixed-layer deepening and F v is the diapycnal eddy diffusion flux. The time period chosen to derive the fluxes required to calculate N spanned 25 d between 10 March and 3 April and included 285 individual profiles (glider dives 10 to 147) spanning the entire survey domain. The first nine dives were omitted because they were used for glider flight trimming. Table 1 gives the uncertainties considered for each budget term. Mean values over the 25 d period are shown in Table 2.

Inventory change I / t
Daily mean column inventories were integrated between surface and z lim for all glider profiles within a 4 d moving window centred on each date. As integration depth z lim , we used the mean euphotic depth of 46 m, derived from satellite measurements using the method described by Lee et al. (2007). Inventory changes ( I t ) were calculated from the day-to-day inventory differences.

Advection F adv
The glider measured within an Eulerian framework, sampling geographic locations over time. In contrast to a Lagrangian framework, estimates of F adv are needed to close the budget. Advection was calculated following Alkire et al.
(2014) using zonal and meridional mean horizontal concen- tration gradients and current velocities within z lim using a moving time window of 8 d. This longer time window of 8 d was chosen to reduce the error in the mean gradients ( Fig. 7). Mean horizontal concentration gradients A c (z) and B c (z) were calculated for each 5 m depth bin using robust plane fits of all concentration measurements in the 8 d time window: where d c (z) is a constant (not used), and x and y are Cartesian coordinates calculated using zonal and meridional distances, respectively. Bisquare weights were applied when determining the fit coefficients (Gross, 1977) to limit the effect of outliers. The zonal (east-west) concentration gradients were small for both O 2 and DIC (Fig. 8). Opposing meridional (northsouth) gradients for O 2 and DIC indicate biological patchiness. However, interestingly the positive gradient for O 2 (corresponding to higher concentrations in the north) is opposite to the pattern seen in the satellite-derived ocean colour (Fig. 1b), which shows lower chlorophyll a concentrations in the north (implying that the bloom propagates from south to north). This highlights that primary production and net community production are not always directly correlated.
To estimate the absolute geostrophic currents within the survey domain, the dynamic height anomaly ( ) was calculated relative to the surface using glider salinity, temperature and pressure (Roquet et al., 2015), and planes were fitted for a moving time window of 8 d: Using these plane fits, the meridional and zonal components of geostrophic shear were derived and referenced to the 8 d running means of the meridional and zonal components of OSCAR velocities (Bonjean and Lagerloef, 2002). Referencing the geostrophic shear with glider dive-averaged currents (DACs) was initially explored, but when compared with satellite data products (Fig. 8c, d) it was clear that the meridional velocities (v) were unrealistic (see Sect. 7).
The advection flux (F adv ) was calculated from the horizontal concentration gradients and the absolute geostrophic velocities: The mean gradients using measurements binned every 5 m between the surface and z lim were averaged at each time step and used for estimating advection (Fig. 8a, b).

Air-sea gas exchange F ase
Air-sea gas exchange (negative for fluxes from atmosphere to ocean) was estimated for O 2 and DIC using with the equilibrium bubble correction ∆ bubble = 0.01[w 10 /(9 m s −2 )] 2 (Woolf and Thorpe, 1991) and   The gas transfer velocity k is calculated according to Wanninkhof (2014). c surf is the mean concentration in the top 10 m over a moving 4 d time window. The mixed-layer depth z mix was estimated using the potential density criteria of the algorithm of Holte and Talley (2009), with a threshold of 0.03 kg m −3 and gradient of 0.0005 kg m −3 dbar −1 . These criteria are sensitive enough to detect the actual mixedlayer depth relevant for gas exchange fluxes (Brainerd and Gregg, 1995). c sat (O 2 ) is the O 2 air-saturation concentration of García and Gordon (1992), corrected to actual sealevel pressure, p baro . c sat (CO 2 ) = χ (CO 2 )p baro F (CO 2 ) is the CO 2 air-saturation concentration, calculated using the dry mole fraction of atmospheric CO 2 , χ (CO 2 ), at the nearest NOAA Carbon Cycle network station Lampedusa, Italy, in March 2016 (ftp://aftp.cmdl.noaa.gov/data/trace_gases/co2/ flask/surface, last access: 20 August 2022, site code: LMP, Dlugokencky et al., 2021). F (CO 2 ) is the CO 2 solubility function in mol dm −3 atm −1 (Weiss and Price, 1980). The last term of the equations is a scaling factor that apportions the air-sea exchange flux to the depth layer of interest (z lim = 46 m) when z mix is deeper than z lim .

Entrainment F e
When the mixed layer deepens to depths greater than z lim , water from below with different concentrations mixes with water within the layer. As a result, there is either a reduction or an increase in mean concentration within z lim depending on the vertical concentration distribution. To account for this, F e was estimated following Binetti et al. (2020): where I 1 (z mix (t 2 )) is the inventory at time t 1 within the mixed-layer depth z mix at time t 2 , I 1 (z lim ) is the inventory at t 1 within z lim and t 2 − t 1 represents the time difference. Mean values of z mix and c averaged over a 4 d moving time window were used. The mean F e (O 2 ) was (−65 ± 1) mmol m −2 d −1 and the mean F e (DIC) was (−74 ± 1) mmol m −2 d −1 for the 7 d in the 25 d period when entrainment (mixed-layer deepening) occurred and z mix > z lim .

Vertical diapycnal eddy diffusion F v
Vertical diapycnal eddy diffusion across the base of z lim was estimated using the measured concentration-density gradients (Copin-Montégut, 2000) and an estimate of the turbulence energy dissipation rate : where the concentration gradient was evaluated using centred finite differences at z lim ± 5 m. = 1.5 × 10 −9 m 2 s −3 was taken from measurements in the Mediterranean bottom pycnoline (Cuypers et al., 2012). Copin-Montégut (2000) guessed a 33 times higher value of 5 × 10 −8 m 2 s −3 in her study of net community production at DyFAMed, but this value was probably an overestimate because even in the presence of eddies, Cuypers et al. (2012) estimate to be only 8.5 × 10 −9 m 2 s −3 . Concentration gradients were calculated for each profile within a 4 d moving window and then averaged to get one gradient per time step. The mean F v (O 2 ) was (−0.8±0.2) mmol m −2 d −1 and the mean F v (DIC) was (−0.1±0.1) mmol m −2 d −1 over the 25 d period. With the higher value of 8.5 × 10 −9 m 2 s −3 found in eddies, these fluxes would be 5.7 times higher, but still negligible with respect to the other four terms in the mass budget.

Buoy measurements
We calculate N from c(O 2 ) and c(DIC) at the BOUSSOLE buoy (N b ) for comparison with N g estimates. N b was estimated similarly to N g , in that we follow the budget approach defined in Sect. 4, and equations defined thereafter. Inventory changes were calculated from surface observations multiplied by z lim because depth profiles are not measured at the buoy. When z mix < z lim , this most likely results in an overestimate of the actual inventory change because N decreases with depth. Also, F ase was scaled in the same way as the glider-based estimates when z mix > z lim . Finally, F adv , F e and F v could not be derived from the single-depth measurements at BOUSSOLE. Instead, the glider-derived fluxes were also used to compute N b .

Net community production N
Prior to 21 March, N g and N b were negative or close to zero most of the time (Fig. 9e). This could be a sign of local net heterotrophic conditions, but could also be due to underestimated entrainment (F e ) or vertical eddy diffusion fluxes (F v ). This observation cannot be explained by physical undersaturation driven by recent cooling because temperatures decreased by only 0.2 • C between 1 and 21 March, explaining at most 0.4 % undersaturation. After 21 March, N exceeds zero until the end of the deployment period. The approximate start of the spring bloom was 19 March as evidenced by the buoy c(DIC) (Fig. 4), which is also the time when N b (DIC) becomes positive, lasting until the end of the deployment period. The contributions of F adv and F e are significant at times. F e (O 2 ) < −80 mmol m −2 d −1 on days when wind is high and the surface mixed layer is deepening (Fig. 9), and there are periods when absolute F adv (O 2 ) > 50 mmol m −2 d −1 . Absolute F adv (DIC) is on average 34 mmol m −2 d −1 higher than F adv (O 2 ), and on some days absolute F adv (DIC) > 100 mmol m −2 d −1 . F e (DIC) is generally of similar magni- tude to absolute F e (O 2 ). Surface waters are undersaturated with O 2 prior to 29 March (negative F ase ) and then oversaturated until the end of the deployment period (positive F ase ). Surface waters are undersaturated in CO 2 throughout the deployment period (negative F ase ).

Stoichiometric relationship
Dividing N(O 2 ) by N(DIC) provides the photosynthetic quotient Q P , previously estimated as 1.45 ± 0.15 (Laws, 1991;Anderson and Sarmiento, 1995;Anderson, 1995). The average Q P in this study using N g (O 2 ) and N g (DIC) between 10 March and 3 April was 0.14 ± 0.81, excluding periods when N g (DIC) < 30 mmol m −2 d −1 (i.e. close to zero within its uncertainty). When using N b (O 2 ) and N b (DIC) over the same time period, average Q P was 0.49 ± 0.60. Using linear regression of N(O 2 ) against N (DIC) provides another way to derive Q P . The resulting slopes are 0.25 ± 0.06 for the glider and 0.54 ± 0.06 for the buoy N estimates. These Q P values therefore do not match the canonical values of 1.45 ± 0.15, but similarly low values of 0.63 ± 0.30 were observed by Copin-Montégut (2000) in early May 1995 at the same site (Table 2).
The change of nutrient concentrations provides an alternative to estimating net community production. Assuming the observed NO − to (128±90) mmol m −2 d −1 , indicating nutrient consumption in line with the assumed stoichiometric ratio. The discrepancy between Redfield and observed stoichiometry is larger for the O 2 -based net community production, especially for the glider-based values. Non-Redfield ratios have been documented before, including at the same site. For example, Copin-Montégut (2000) found non-Redfield ratios for O 2 and DIC-derived N at DyFAMed. Hull et al. (2021) found non-Redfield ratios for O 2 -and nitrate-derived N in the North Sea. We do not have a physiological explanation for such large deviations from the canonical stoichiometric values. It is unlikely to be due to a calibration error in the O 2 concentration because it would have to be of the order of 20 mmol m −3 , and our observations at depths > 400 m agree to within 3 mmol m −3 with the record of historic measurements.

Discussion
Quantifying net community production helps us to understand the role of biological production and consumption on carbon export from the surface to deep waters. If the rate of carbon export is known, the rate of atmospheric CO 2 drawdown into the ocean can be better quantified, and the accuracy of future climate projections might improve. Using underwater gliders as tools to observe the water column on timescales of less than a month over a wider area allows us to estimate physical processes (e.g. mixing, advection) that affect biogeochemical tracer concentrations and estimates of net community production derived therefrom. Table 2. Comparison of I / t , F adv , F ase , F v and N estimates using c(O 2 ) and c(DIC) and the corresponding flux terms used in this study with others in the literature. The N estimates of Copin-Montégut (2000) are based on I / t, F ase and F v (calculated with a turbulence energy dissipation rate of = 5.0 × 10 −8 m 2 s −3 ).

Time period
Variable During the spring bloom (19 March to 3 April) we estimate N between (9 ± 36) and (128 ± 90) mmol m −2 d −1 on average, with maxima > 300 mmol m −2 d −1 at times. Few studies have estimated N at the DyFAMed site, and there are currently no N estimates that incorporate F adv using concentrations over a larger area surrounding DyFAMed.
We can compare our study with other studies elsewhere that estimate N using glider measurements. Alkire et al. (2014)  Uncertainties for N (DIC) are generally higher than those for N (O 2 ), and the highest uncertainties are associated with F adv . F adv uncertainties are high because the uncertainties associated with the concentration plane fits are high due to c(DIC) variability within the time-centred window. For example, c(DIC) within z lim varied over 33 mmol m −3 for the time-centred window on 16 March, leading to low standard errors (Fig. 8b), while c(DIC) varied over 40 mmol m −3 for the time-centred window on 26 March, leading to high standard errors. For this latter window centred on 26 March, the decrease in c(DIC) over time between 25 and 30 March caused this high uncertainty. We use an 8 d time-centred moving window to reduce the potential effect of biologically related concentrations on the mean gradients, typically seen on short timescales. However, as demonstrated here, a tradeoff is the higher uncertainty. Additional gliders deployed in future missions could potentially decrease the uncertainty associated with F adv by combination of alternative sampling strategies, leading to more measurements in horizontal space and thus improving the robustness of the plane fits.
We initially explored using the glider DACs to reference geostrophic shear, required to estimate F adv . However, DAC v velocities were positive throughout most of the deployment period in an area where negative v velocity is expected, as seen in satellite data products (Fig. 8). We think the erroneous DAC v velocities relate to the glider's roll during flight, which was consistently > 100 • starboard. The disadvantage of using OSCAR velocities to derive velocity profiles relates to resolution; OSCAR velocities are derived every 8 d on a grid with 0.25 • resolution.
The N b estimates incorporate measurements at a depth of approximately 10 m that we extrapolate over the average euphotic depth (z lim = 46 m), while N g estimates account for the vertical variability between the surface and z lim . This means that the glider-based N estimates should in principle be more accurate. It also explains why the buoy-based inventory changes are larger than for the glider, due to concentrations varying more near the surface than at depth.
We chose z lim = 46 m as the integration depth for estimating N. We investigated the sensitivity of the integration depth to N estimates by estimating N over the top 36 m and top 56 m of the water column. We found that N estimates are not very sensitive to the integration depth; e.g. using z lim = 36 m, we obtain N values between (11 ± 28) and (113±75) mmol m −2 d −1 for the period 19 March to 3 April. Using z lim = 56 m, we obtain N values between (3 ± 44) and (147 ± 109) mmol m −2 d −1 .

Conclusions
The spring bloom started in this area of the northwestern Mediterranean Sea around 19 March 2016. We identified the spring bloom through a decrease in the DIC concentration (inferred from a pH increase) and an increase in the oxygen concentration. The spring bloom followed a period of high winds and a well-mixed higher nutrient water column. The glider also detected increased optical backscatter with spikes on 19 and 20 March, an indication of increased particle concentrations. The spring bloom start date fell into the time frame of phytoplankton blooms in previous years according to buoy-based sensors. However, this was the first time that high-resolution glider profiles were available, providing insights into the biogeochemical and physical processes at depth and in waters horizontally surrounding the DyFAMed mooring.
This study demonstrates the potential of using autonomous glider measurements for estimating net community production (N). It is the first study of this kind in the northwestern Mediterranean Sea. Both oxygen and dissolved inorganic carbon concentrations were used to derive mass budgets, taking into account N , air-sea gas exchange, horizontal advection, vertical entrainment and diapycnal eddy diffusion. Estimating N using different data sets was challenging, which is reflected in the range of N estimates and their uncertainties presented here. We found that day-to-day variations in derived N estimates were dominated by net inventory changes, air-sea exchange and horizontal advection. Entrainment only came into play on 7 out of 25 d, i.e. during periods when the mixed-layer depth deepened beyond the integration limit of the mean euphotic zone depth of 46 m. Diapycnal eddy diffusion was negligible at all times. Over the duration of the deployment from 10 March to 3 April, we estimate N over the mean euphotic zone depth to be (−17 ± 36) mmol m −2 d −1 based on glider O 2 , (44 ± 94) mmol m −2 d −1 based on glider DIC, (−17 ± 37) mmol m −2 d −1 based on buoy O 2 and (49 ± 86) mmol m −2 d −1 based on buoy DIC. Within their uncertainties these values agree and show net metabolic balance. Unlike previous studies in this region of the northwestern Mediterranean Sea, our results and their uncertainties explicitly consider the role of horizontal advection.
We showed that advection can potentially be significant in the area of our study, bearing in mind that the associated uncertainty was high. Daily mean N values estimated from point measurements typical of fixed moorings were affected by advection on some days. Advection may also be relevant on timescales longer than 1 month, but the current study was too short to explore this. However, when the advection flux was averaged over the 25 d period, there was less of an effect on N, suggesting that at least in this region, its importance might diminish for longer integration periods. To better resolve the influence of advection on shorter timescales and reduce its uncertainty, similar future studies should increase the number of gliders circumnavigating the central mooring so that averaging intervals can be reduced and the time resolution can be enhanced.
This was the second test deployment of the prototype IS-FET sensors on gliders, following our first test in 2014 . There were still considerable problems with drift and a need for calibration. Since then, Saba et al. (2018) and Takeshita et al. (2021) had more success with a Sea-Bird Scientific ISFET sensor and a Deep-Sea-DuraFET sensor, respectively. However, these sensors are not yet commercially available for gliders.

Appendix A: Discrete water samples
Water samples used to measure c(O 2 ) were collected at the DyFAMed site on 7 March (n = 10) and 16 April 2016 (n = 8) using 12 L Niskin bottles (General Oceanics 1010X) from the top 1000 m of the water column; four samples each were from the top 100 m (Fig. B1b). Reagents needed for the fixation of oxygen were added to the samples at the time of water sample collection on board the ship. An automated Winkler titration method with endpoint detection was used after each cruise in the laboratory at the Observatoire Océanologique de Villefranche-sur-Mer, France, to determine c(O 2 ). Replicates were obtained to determine precision. c(O 2 ) measured by the rosette-mounted SBE43 sensor within the top 150 m was 15 mmol m −3 higher than the Winkler measurements (Fig. B1b). SBE 43 sensor c(O 2 ) was corrected by regressing against the Winkler c(O 2 ) values. The resulting regression coefficients were 0.892 ± 0.009 for slope and (18.7 ± 2.0) mmol m −3 for intercept on 7 March (r 2 = 0.9993; σ = 0.9 mmol m −3 ). The corresponding values on 16 April were 0.920 ± 0.008 and (11.8 ± 1.6) mmol m −3 (r 2 = 0.9997; σ = 0.7 mmol m −3 ).
A Marianda Versatile INstrument for the Determination of Titration Alkalinity (VINDTA 3C; https://www.marianda. com, last access: 20 August 2022) was used to measure c(DIC) and A T at 10 depth levels. A total of 19 bottles of certified reference material (CRM) supplied by Scripps Institution of Oceanography (San Diego, California, USA) were run to calibrate the instrument. Coulometry following standard operating procedure SOP 2 was used to measure c(DIC) (Dickson et al., 2007;Johnson et al., 1985), and potentiometric titration following SOP 3b was used to measure A T (Dickson et al., 2007).
Nutrients at DyFAMed were measured using water collected in the same Niskin bottles as those used for c(O 2 ), c(DIC) and A T . Samples were analysed via a standard automated colorimetry system (Seal Analytical continuous flow AutoAnalyser III) at Observatoire Océanologique de Villefranche-sur-Mer for NO − x (nitrate NO − 3 and nitrite NO − 2 combined) (Bendschneider and Robinson, 1952), Si(OH) 4 (Murphy and Riley, 1962) and PO 3− 4 (Strickland and Parsons, 1972), with detection limits of 0.01, 0.02 and 0.02 mmol m −3 , respectively (de Fommervault et al., 2015).

Appendix B: Glider calibrations B1 O 2 measurements
Glider c(O 2 ) values were calibrated to take into account the response time (τ ) of the sensor, which is dependent on the thickness, and usage of the sensor foil (McNeil and D'Asaro, 2014), as well as temperature, and to account for the difference between the glider measurements and c(O 2 ) obtained by the ship (Fig. B1b) (Binetti et al., 2020;Hemming et al., 2017). To correct c(O 2 ) measured by the glider sensor, the sensor-related oxygen engineering parameters TCPhase and CalPhase were used. A mean τ of 8 s was applied to correct measurements for the sensor time lag. This mean τ value was determined from the lowest root mean square difference between descending and ascending TCPhase profiles shifted in time by values of τ ranging from 0 to 50 s. After lag correction of the glider TCPhase, the relationship between the ship sensor pseudo-CalPhase measured on 7 March and on 16 April and the glider TCPhase measured on 8 March and on 4 April was determined (Fig. B1a). The ship sensor pseudo-CalPhase was calculated from the calibrated ship sensor c(O 2 ) by inverting the set of equations normally used to obtain c(O 2 ) from glider optode TCPhase and CalPhase. The slope and offset coefficients obtained from the relationship between the glider TCPhase over the entire deployment time period and ship sensor pseudo-CalPhase on 7 March and 16 April (Fig. B1a) were linearly interpolated over the duration of the deployment. The interpolated coefficients were matched with glider measurements in time and used to correct all glider CalPhase profiles obtained during the deployment, allowing a calibration of glider c(O 2 ) (Fig. B1b). After calibration, the glider c(O 2 ) agreed with buoy c(O 2 ) (Fig. B1c), as would be expected because the same Winkler samples were used to calibrate ship and buoy oxygen sensors.

B2 ISFET pH measurements
The deployment offered a second opportunity to trial two prototype ISFET pH-p(CO 2 ) sensor pairs, previously tested on an underwater glider in the Sardinian Sea during the REP14-MED experiment Onken et al., 2018). The sensors are custom-built non-commercially by Kiminori Shitashima's group at Tokyo University of Marine Science and Technology, Japan. ISFET sensors measure pH using the interface potential between a reference chlorine electrode (Cl-ISE) and a semiconducting ion sensing transistor .
One pH-p(CO 2 ) sensor pair was stand-alone, meaning measurements were logged and stored by the sensor and retrieved after the deployment. Another sensor pair was integrated into the glider electronics, allowing measurements to be sent remotely by satellite in near-real time. Both the stand-alone and integrated sensors were positioned on the underside of the glider to limit the effect of sunlight on measurements, and backup batteries were provided to supply power in between sampling . Unfortunately, the stand-alone sensor ceased operating after less than 3 d around 07:00 CET on 10 March because of an issue with its power supply. For this reason, pH T obtained by the stand-alone sensor was not used. The sensors were placed in a bucket of locally collected coastal surface seawater for a period of 10 h before deployment. Several weeks of predeployment conditioning have been recommended (Bresnahan et al., 2014;Takeshita et al., 2014), but this was not possible due to time constraints.
A comparison of sensor measurements with pH T calculated from discrete c(DIC) and A T measurements both during the deployment and with archived data from between 1998 and 2013 at the DyFAMed site (http://mistrals.sedoo. fr/Data-Download, last access: 20 August 2022) indicated problems in accuracy and stability (Fig. B2). The pH T of the integrated sensor drifted by 0.8 over the course of the 29 d deployment, with a depth-dependent range in pH T of 1. pH T measurements were corrected for drift using the mean difference (offset) between dive 10 pH T measured between 350 and 950 m and pH T measured at the same depths during each subsequent dive. This method assumes that pH T variations over the 29 d deployment between 350 and 950 m are negligible with respect to the sensor precision. The offsets were added to the full-depth glider pH T profiles, and the range of pH T was reduced 3-fold to approximately 0.3 (Fig. B2). After correcting pH T for drift, pH T was corrected for pressure using a pressure coefficient of +5.6×10 −5 dbar −1 . After adding this correction term, the range of pH T was reduced to approximately 0.2 (Fig. B2).
Finally, integrated ISFET pH T values were corrected relative to "true" discrete pH T . Regression coefficients calculated using glider drift-and pressure-corrected pH T and ship pH T collected on 7 March and on 16 April were used to correct glider pH T . The range of glider pH T was now matched with that of historical discrete samples (Fig. B2).

Appendix C: Deriving parameters using CO2SYS
The CO2SYS software package (Van Heuven et al., 2011;Orr et al., 2015) was used to derive discrete pH T for correcting glider ISFET pH T , using discrete measurements of c(DIC), A T , in situ temperature, salinity, pressure, and Si(OH) 4 and PO 3− 4 concentrations. Equilibrium constants (Mehrbach et al., 1973;Dickson and Millero, 1987), sulfate acid dissociation constant (Dickson, 1990) and total borate concentration (Uppström, 1974) were used, as recommended by previous Mediterranean-based studies (Álvarez et al., 2014;Key et al., 2010). c(DIC) and A T had an uncertainty of 3.5 and 3.2 µmol kg −1 , respectively, determined from CRM measurements, representing a combined uncertainty of 0.008 in derived pH T .
Calibrated glider pH T and A T were used to derive glider f (CO 2 ) and c(DIC) using CO2SYS. The mean uncertainties of derived f (CO 2 ) and c(DIC) were 5.3 µatm and 4.1 µmol kg −1 , respectively, calculated as the mean of the absolute differences of f (CO 2 ) and c(DIC) derived using single values of A T and pH (including the margin of error associated with the A T -S relationship and the pH T pressure correction, respectively).
To obtain buoy c(DIC), observed f (CO 2 ) at 10 m was used with A T calculated from S. The uncertainty of derived buoy c(DIC) was 6.6 µmol kg −1 , calculated as the absolute difference between the ship c(DIC) sample collected on 7 March at 10 m depth and average c(DIC) measured by the buoy on the same day. Figure B2. A comparison of glider pH T after various corrections have been applied. Raw pH T (blue) is corrected for drift (turquoise) and subsequently for pressure effects (yellow). Finally pH T is compared with "true" ship pH T and corrected (purple). Historical pH T anomalies (black) collected in March and April between 1998 and 2013 are shown for comparison, alongside the ship pH T on 7 March (green squares) and 16 April (pink diamonds).
Author contributions. This work resulted from MPH's PhD project at the University of East Anglia under the supervision of JK, KJH, DCEB and JB. GAL, MCG and KS implemented the ISFET sensors on the glider and deployed it. JK led the UEA piloting team. DA acquired funding for the Bio-optics and Carbon Experiment (BIOCAREX) that supports the CO 2 measurements at BOUSSOLE and led ship operations and ancillary data collection. The data were analysed by MPH and JK, who also wrote the manuscript. All other authors contributed to further editing and writing. from the Agence Nationale de la Recherche (ANR, Paris) to the Bio-optics and Carbon Experiment (BIOCAREX).
Review statement. This paper was edited by Vanessa Cardin and reviewed by three anonymous referees.