Norwegian Sea net community production estimated from O₂ and prototype CO₂ optode measurements on a Seaglider

We report on a pilot study using a CO₂ optode deployed on a Seaglider in the Norwegian Sea from March to October 2014. The optode measurements required drift and lag correction and in situ calibration using discrete water samples collected in the vicinity. We found that the optode signal correlated better with the concentration of CO₂, c(CO₂), than with its partial pressure, p(CO₂). Using the calibrated c(CO₂) and a regional parameterisation of total alkalinity (A_T) as a function of temperature and salinity, we calculated total dissolved inorganic carbon content, c(DIC), which had a standard deviation of 11 µmol kg⁻¹ compared with in situ measurements. The glider was also equipped with an oxygen (O₂) optode. The O₂ optode was drift corrected and calibrated using a c(O₂) climatology for deep samples. The calibrated data enabled the calculation of DIC- and O₂-based net community production, N(DIC) and N(O₂). To derive N, DIC and O₂ inventory changes over time were combined with estimates of air–sea gas exchange, diapycnal mixing and entrainment of deeper waters. Glider-based observations captured two periods of increased Chl a inventory in late spring (May) and a second one in summer (June). For the May period, we found N(DIC) = (21±5) mmol m⁻² d⁻¹, N(O₂) = (94±16) mmol m⁻² d⁻¹ and an (uncalibrated) Chl a peak concentration of c_raw(Chl a) = 3 mg m⁻³. During the June period, c_raw(Chl a) increased to a summer maximum of 4 mg m⁻³, associated with N(DIC) = (85±5) mmol m⁻² d⁻¹ and N(O₂) = (126±25) mmol m⁻² d⁻¹. The high-resolution dataset allowed for quantification of the changes in N before, during and after the periods of increased Chl a inventory. After the May period, the remineralisation of the material produced during the period of increased Chl a inventory decreased N(DIC) to $(-\mathrm{3}\pm \mathrm{5})$ mmol m⁻² d⁻¹ and N(O₂) to (0±2) mmol m⁻² d⁻¹. The survey area was a source of O₂ and a sink of CO₂ for most of the summer. The deployment captured two different surface waters influenced by the Norwegian Atlantic Current (NwAC) and the Norwegian Coastal Current (NCC). The NCC was characterised by lower c(O₂) and c(DIC) than the NwAC, as well as lower N(O₂) and c_raw(Chl a) but higher N(DIC). Our results show the potential of glider data to simultaneously capture time- and depth-resolved variability in DIC and O₂ concentrations.

1 Introduction

Climate models project an increase in the atmospheric CO₂ mole fraction driven by anthropogenic emissions from a pre-industrial value of 280 µmol mol⁻¹ (Neftel et al., 1982) to 538–936 µmol mol⁻¹ by 2100 (Pachauri and Reisinger, 2007). The ocean is known to be a major CO₂ sink (Sabine et al., 2004; Le Quéré et al., 2009; Sutton et al., 2014); in fact, it has taken up approximately 25 % of this anthropogenic CO₂ with a rate of (2.5 ± 0.6) Gt a⁻¹ (in C equivalents) (Friedlingstein et al., 2019). This uptake alters the carbonate system of seawater and is causing a decrease in seawater pH, a process known as ocean acidification (Gattuso and Hansson, 2011). The processes affecting the marine carbonate system include air–sea gas exchange, photosynthesis and respiration, advection and vertical mixing, and CaCO₃ formation and dissolution. For that reason, it is important to develop precise, accurate and cost-effective tools to observe CO₂ trends, variability and related processes in the ocean. Provided that suitable sensors are available, autonomous ocean glider measurements may help resolve these processes.

To quantify the marine carbonate system, four variables are commonly measured: total dissolved inorganic carbon concentration, c(DIC); total alkalinity, A_T; the fugacity of CO₂, f(CO₂); and pH. At thermodynamic equilibrium, knowledge of two of the four variables is sufficient to calculate the other two. Marine carbonate system variables are primarily measured on research ships, commercial ships of opportunity, moorings, buoys and floats (Hardman-Mountford et al., 2008; Monteiro et al., 2009; Takahashi et al., 2009; Olsen et al., 2016; Bushinsky et al., 2019). Moorings equipped with submersible sensors often provide limited vertical and horizontal but good long-term temporal resolution (Hemsley, 2003). In contrast, ship-based surveys have higher vertical and spatial resolution than moorings but limited repetition frequency because of the expense of ship operations. Ocean gliders have the potential to replace some ship surveys because they are much cheaper to operate and will increase our coastal and regional observational capacity. However, the slow glider speed of 1–2 km h⁻¹ only allows a smaller spatial coverage than ship surveys, and the sensors require careful calibration to match the quality of data provided by ship-based sampling.

Carbonate system sensors suitable for autonomous deployment have been developed in the past decades, in particular pH sensors (Seidel et al., 2008; Martz et al., 2010; Rérolle et al., 2013) and p(CO₂) sensors (Atamanchuk et al., 2014; Bittig et al., 2012; Degrandpre, 1993; Goyet et al., 1992; Körtzinger et al., 1996). One of these sensors is the CO₂ optode (Atamanchuk et al., 2014), which has been successfully deployed to monitor an artificial CO₂ leak on the Scottish west coast (Atamanchuk, et al., 2015b), on a cabled underwater observatory (Atamanchuk et al., 2015a), to measure lake metabolism (Peeters et al., 2016), for fish transportation (Thomas et al., 2017) and on a moored profiler (Chu et al., 2020).

c(DIC) and c(O₂) measurements can be used to calculate net community production (NCP), which is defined as the difference between gross primary production (GPP) and community respiration (CR). At steady state, NCP is equal to the rate of organic carbon export and transfer from the surface into the mesopelagic and deep waters (Lockwood et al., 2012). NCP is derived by vertical integration to a specific depth that is commonly defined relative to the mixed layer depth (z_mix) or the bottom of the euphotic zone (Plant et al., 2016). A system is defined as autotrophic when GPP is larger than CR (i.e. NCP is positive) and as heterotrophic when CR is larger than GPP (i.e. NCP is negative) (Ducklow and Doney, 2013).

NCP can be quantified using bottle incubations or in situ biogeochemical budgets (Sharples et al., 2006; Quay, et al, 2012; Seguro et al., 2019). Bottle incubations involve measuring production and respiration in vitro under dark and light conditions. Biogeochemical budgets combine O₂ and DIC inventory changes with estimates of air–sea gas exchange, entrainment, advection and vertical mixing (Neuer et al., 2007; Alkire et al., 2014; Binetti et al., 2020).

The Norwegian Sea is a complex environment due to the interaction between the Norwegian Atlantic Current (NwAC) entering from the south-west, Arctic Water coming from the north and the Norwegian Coastal Current (NCC) flowing along the Norwegian coast (Nilsen and Falck, 2006). In particular, Atlantic Water enters the Norwegian Sea through the Faroe–Shetland Channel and Iceland–Faroe Ridge (Hansen and Østerhus, 2000) with salinity S between 35.1 and 35.3 and temperatures (θ) warmer than 6 ^∘C (Swift, 1986). The NCC water differs from the NwAC with a surface S<35 (Saetre and Ljoen, 1972) and a seasonal θ signal (Nilsen and Falck, 2006).

Biological production in the Norwegian Sea varies during the year and five different periods can be discerned (Rey, 2004): (1) winter with the smallest productivity and phytoplankton biomass; (2) a pre-bloom period; (3) the spring bloom when productivity increases and phytoplankton biomass reaches the annual maximum; (4) a post-bloom period with productivity mostly based on regenerated nutrients; (5) autumn with smaller blooms than in summer. Previous estimates of the DIC-based net community production (N(DIC)) were based on discrete c(DIC) samples (Falck and Anderson, 2005) or calculated from c(O₂) measurements and converted to C equivalents assuming Redfield stoichiometry of production or respiration (Falck and Gade, 1999; Skjelvan et al., 2001; Kivimäe, 2007). Glider measurements have been used to estimate NCP in other ocean regions (Nicholson et al., 2008; Alkire et al., 2014; Haskell et al., 2019; Binetti et al., 2020); however, as far as we know, this is the first study of net community production in the Norwegian Sea using a high-resolution glider dataset (>10⁶ data points; 40 s time resolution) and the first anywhere estimating NCP from a glider-mounted sensor directly measuring the marine carbonate system.

2 Materials and methods

2.1 Glider sampling

Kongsberg Seaglider 564 was deployed in the Norwegian Sea on 16 March 2014 at 63.00^∘ N, 3.86^∘ E, and recovered on 30 October 2014 at 62.99^∘ N, 3.89^∘ E. The Seaglider was equipped with a prototype Aanderaa 4797 CO₂ optode, an Aanderaa 4330F oxygen optode (Tengberg et al., 2006), a Sea-Bird conductivity–temperature–depth profiler (CTD) and a combined backscatter–chlorophyll a fluorescence sensor (Wetlabs Eco Puck BB2FLVMT). The mean sampling intervals for each sensor varied with depth (Table 1).

Table 1 Average sampling interval of Sea-Bird CTD, Aanderaa 4330F oxygen optode, Aanderaa 4797 CO₂ optode and a combined backscatter/chlorophyll a fluorescence sensor (Wetlabs Eco Puck BB2FLVMT) in the top 100, from 100 to 500 and from 500 to 1000 m.

The deployment followed the Svinøy trench from the open sea towards the Norwegian coast. The glider covered a 536 km long transect eight times (four times in each direction) for a total of 703 dives (Fig. 1).

Figure 1 Map of the glider deployment and the main currents. The black dots are the glider dives; the green and the red dots are the water samples collected along the glider section and at Ocean Weather Station M (OWSM), respectively. The three main water masses (Skjelvan et al., 2008) are the Norwegian Coastal Current (yellow), the Norwegian Atlantic Current (NwAC, orange) and Arctic Water (green).

2.2 Discrete sampling

During the glider deployment, 70 discrete water samples from various depths (5, 10, 20, 30, 50, 100, 300, 500 and 1000 m) were collected on five different cruises on the R/V Haakon Mosby along the southern half of the glider transect on 18 March, 5 May, 6 and 14 June, and 30 October 2014. Samples for c(DIC) and A_T were collected from 10 L Niskin bottles following standard operational procedure (SOP) 1 of Dickson et al. (2007). The c(DIC) and A_T samples were preserved with saturated HgCl₂ solution (final HgCl₂ concentration: 15 mg dm⁻³) and analysed within 14 d after the collection. Nutrient samples from the same Niskin bottles were preserved with chloroform (Hagebo and Rey, 1984). c(DIC) and A_T were analysed on shore according to SOPs 2 and 3b (Dickson et al., 2007) using a VINDTA 3D (Marianda) with a CM5011 coulometer (UIC instruments) and a VINDTA 3S (Marianda), respectively. The precision of the samples' c(DIC) and A_T values was 1 µmol kg⁻¹ for both, based on duplicate samples and running Certified Reference Material (CRM) batch numbers 118 and 138 provided by professor A. Dickson, Scripps Institution of Oceanography, San Diego, USA (Dickson et al., 2003). Nutrients were analysed on shore using an Alpkem AutoAnalyzer. In addition, 43 water samples were collected at Ocean Weather Station M (OWSM) on five different cruises: 22 March on R/V Haakon Mosby, 9 May on R/V G.O. Sars, 14 June on R/V Haakon Mosby, 2 August and 13 November 2014 on R/V Johan Hjort from 10, 30, 50, 100, 200, 500, 800 and 1000 m depths. The OWSM samples were preserved and analysed for A_T and c(DIC) as the Svinøy samples. No phosphate and silicate samples were collected at OWSM. Temperature (θ) and salinity (S) profiles were measured at each station using a Sea-Bird 911 plus CTD. pH and f(CO₂) were calculated using the MATLAB toolbox CO2SYS (Van Heuven et al., 2011), with the following constants: K₁ and K₂ carbonic acid dissociation constants of Lueker et al. (2000), K(HSO${}_{\mathrm{4}}^{-}$/SO${}_{\mathrm{4}}^{\mathrm{2}-}$) bisulfate dissociation constant of Dickson (1990) and borate to chlorinity ratio of Lee et al. (2010). The precision of A_T and c(DIC) led to an uncertainty in the calculated c(CO₂) of 0.28 µmol kg⁻¹. For the OWSM calculations, we used nutrient concentrations from the Svinøy section at a time as close as possible to the OWSM sampling as input. In the case of the glider, we derived a parameterisation for phosphate and silicate concentration as a function of sample depth and time. This parameterisation had an uncertainty of 1.3 and 0.13 µmol kg⁻¹ and a R² of 0.6 and 0.4, for silicate and phosphate concentrations, respectively. The uncertainty was calculated as the root mean square difference between measured and parameterised concentrations. This nutrient concentration uncertainty contributed an uncertainty of 0.04 µmol kg⁻¹ in the calculation of c(CO₂), which is negligible and smaller than the uncertainty caused by A_T and c(DIC).

2.3 Oxygen optode calibration

The last oxygen optode calibration before the deployment was performed in 2012 as a two-point calibration at 9.91 ^∘C in air-saturated water and at 20.37 ^∘C in anoxic Na₂SO₃ solution. Oxygen optodes are known to be affected by drift (Bittig and Körtzinger, 2015), which is even worse for the fast-response foils used in the 4330F optode for glider deployments. It has been suggested that it is necessary to calibrate and drift correct the optode using discrete samples or in-air measurements (Nicholson and Feen, 2017). Unfortunately, no discrete samples were collected at glider deployment or recovery.

To overcome this problem, we used archived data to correct for oxygen optode drift. These archived concentration data (designated c_C(O₂)) were collected at OWSM between 2001 and 2007 (downloaded from ICES database) and in the glider deployment region between 2000 and 2018 (extracted from GLODAPv2; Olsen et al., 2016). To apply the correction, we used the oxygen samples corresponding to a potential density σ₀ > 1028 kg m⁻³ (corresponding to depths between 427 and 1000 m), because waters of these potential densities were always well below the mixed layer and therefore subject to limited seasonal and interannual variability, as evidenced by the salinity S and potential temperature θ of these samples: S varied from 34.88 to 34.96, with a mean of 34.90 ± 0.01; θ varied from 0.45 to −0.76 ^∘C, with a mean of (−0.15 ± 0.36) ^∘C.

Figure 2 shows that the glider oxygen concentration (c_G(O₂)) corresponding to σ₀ > 1028 kg m⁻³ was characterised by two different water masses separated at a latitude of about 64^∘ N. We used the samples collected north of 64^∘ N to derive the glider optode correction because this reflects the largest area covered by the glider. We did not use the southern region because the archived samples from there covered only 5 d. For each day of the year with archived samples, we calculated the median concentration of the glider and the archived samples. Figure 3 shows a plot of the ratio between $c_{\mathrm{C}}\left({\mathrm{O}}_{\mathrm{2}}\right)/c_{\mathrm{G}}\left({\mathrm{O}}_{\mathrm{2}}\right)$ against the day of the year and a linear fit, which is used to calibrate c_G(O₂) and correct for drift.

No lag correction was applied because the O₂ optode had a fast response foil and showed no detectable lag (< 10 s), based on a comparison between descent and ascent profiles.

Figure 2 Glider oxygen concentration, c_G(O₂), for σ₀ > 1028 kg m⁻³ coloured by latitude.

Figure 3 A linear fit of the ratio between the daily median of the discrete oxygen samples (c_C(O₂)) and glider oxygen data (c_G(O₂)) for σ₀ > 1028 kg m⁻³ was used to derive the c_G(O₂) drift and initial offset at deployment. The time difference Δt is calculated with respect to the deployment day on 16 March.

2.4 CO₂ optode measurement principle

The CO₂ optode consists of an optical and a temperature sensor incorporated into a pressure housing. The optical sensor has a sensing foil comprising two fluorescence indicators (luminophores), one of which is sensitive to pH changes, and the other is not and thus used as a reference. The excitation and emission spectra of the two fluorescence indicators overlap, but the reference indicator has a longer fluorescence lifetime than the pH indicator. These two fluorescence lifetimes are combined using an approach known as dual lifetime referencing (DLR) (Klimant et al., 2001; von Bültzingslöwen et al., 2002). From the phase shift (φ), the partial pressure of CO₂, p(CO₂), is parameterised as an eight-degree polynomial (Atamanchuk et al., 2014):

$$\begin{array}{}\text{(1)}& \log \left[p\right({\mathrm{CO}}_{\mathrm{2}})/\mathrm{\mu }\mathrm{atm}]=C_{\mathrm{0}}+C_{\mathrm{1}}\mathit{\phi }+\mathrm{\dots }+C_{\mathrm{8}}{\mathit{\phi }}^{\mathrm{8}},\end{array}$$

where C₀ to C₈ are temperature-dependent coefficients.

The partial pressure of CO₂ is linked to the CO₂ concentration, c(CO₂) and the fugacity of CO₂, f(CO₂), via the following relationship:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}c\left({\mathrm{CO}}_{\mathrm{2}}\right)& =p\left({\mathrm{CO}}_{\mathrm{2}}\right)/[\mathrm{1}-p({\mathrm{H}}_{\mathrm{2}}\mathrm{O})/p]F\left({\mathrm{CO}}_{\mathrm{2}}\right)\\ & =K_{\mathrm{0}}\left({\mathrm{CO}}_{\mathrm{2}}\right)f\left({\mathrm{CO}}_{\mathrm{2}}\right),\end{array}\end{array}$$

where F(CO₂) is the solubility function (Weiss and Price, 1980), p(H₂O) is the water vapour pressure, p is the total gas tension (assumed to be near 1 atm), and K₀(CO₂) is the solubility coefficient. F and K₀ vary according to temperature and salinity.

2.5 CO₂ optode lag and drift correction and calibration

The CO₂ optode was fully functional between dives 31 (on 21 March 2014) and 400 (on 24 July 2014). After dive 400, the CO₂ optode stopped sampling in the top 150 m. Figure 4 shows the outcome of each calibration step: (0) uncalibrated optode output (blue dots), (1) drift correction (red dots), (2) lag correction (green dots) and (3) calibration using discrete water samples (black dots).

Figure 4 Panel (a) shows the calibrated p(CO₂) (p_cal(CO₂)) in black and the discrete samples in azure. (b) Plot of p(CO₂) versus depth where the continuous vertical lines are the mean every 50 m and the error bars represent the standard deviation. Blue colour shows p_u(CO₂) without any correction; red shows p_d(CO₂) corrected for drift; green represents p_c(CO₂) corrected for drift and lag; black shows p_cal(CO₂) calibrated against water samples (azure dots) collected during the deployment (Sect. 2.5). p_cal(CO₂) had a mean standard deviation of 22 µatm and a mean bias of 1.8 µatm compared with the discrete samples.

In order to correct for the drift occurring during the glider mission, we selected the CO₂ optode measurements in water with σ₀>1028 kg m⁻³ (just as for O₂; Sect. 2.3). We calculated the median of the raw optode phase shift data (“CalPhase” φ_cal) for each Seaglider dive. Then, we calculated a drift coefficient (m_i) as the ratio between the median φ_cal for a given dive divided by the median φ_cal of dive 31. Drift-corrected φ_cal,d values were calculated by dividing the raw φ_cal by the specific m_i for each dive.

The CO₂ optode was also affected by lag (Atamanchuk et al., 2014) caused by the slow response of the optode to ambient c(CO₂) changes in time and depth. The lag created a discrepancy between the depth profiles obtained during glider ascents and descents. To correct for this lag, we applied the method of Miloshevich (2004), which was previously used by Fiedler et al. (2013) and Atamanchuk et al. (2015b) to correct the lag of the Contros HydroC CO₂ sensor (Fiedler et al., 2013; Saderne et al., 2013). This CO₂ sensor has a different measurement principle (infrared absorption) than the CO₂ optode, but both rely on the diffusion of CO₂ through a gas-permeable membrane.

To apply the lag correction, the sampling interval (Δt) needs to be sufficiently small compared to the sensor response time (τ) and the ambient variability (Miloshevich, 2004). Before the lag correction, φ_cal,d was smoothed to remove any outliers and “kinks” in the profile using the Matlab function rLOWESS. The smoothing function applies a local regression every nine points using a weighted robust linear least-squares fit. Subsequently, τ was determined such that the following lag-correction equation (Miloshevich, 2004) minimised the φ_cal,d difference between each glider ascent and the following descent:

$$\begin{array}{}\text{(3)}& p_{\mathrm{c}}\left({\mathrm{CO}}_{\mathrm{2}}{t}_{\mathrm{1}}\right)={\displaystyle \frac{p_{\mathrm{d}}({\mathrm{CO}}_{\mathrm{2}},t_{\mathrm{1}})-p_{\mathrm{d}}({\mathrm{CO}}_{\mathrm{2}},t_{\mathrm{0}})e^{-\mathrm{\Delta }t/\mathit{\tau }}}{\mathrm{1}-e^{-\mathrm{\Delta }t/\mathit{\tau }}}},\end{array}$$

where p_d(CO₂,t₀) is the drift-corrected value measured by the optode at time t₀, p_d(CO₂,t₁) is the measured value at time t₁, Δt is the time between t₀ and t₁, τ is the response time, and p_c(CO₂,t₁) is the lag-corrected value at t₁. We calculated a τ value for each glider dive and used the median of τ (1384 s, 25th quartile: 1101 s; 75th quartile: 1799 s) (Fig. 5), which was larger than Δt(258 s) and therefore met the requirement to apply the Miloshevich (2004) method. To apply the lag correction, the glider needs to sample same water mass during the ascent and descent. The difference between the ascent and descent was minimal: it was (0.13 ± 0.33) ^∘C for θ and 0.02 ± 0.04 for S. This lag correction reduced the average difference between glider ascent and descent from (71 ± 30) to (21 ± 26) µatm.

Figure 5 The histogram shows the distribution of the τ calculated from glider dives 31 to 400 to correct the CO₂ optode drift using the algorithm of Miloshevich (2004).

The CO₂ optode output was calibrated using the discrete samples collected throughout the mission. Using the discrete sample time and potential density σ₀, we selected the closest CO₂ optode output. A linear regression between optode output and c(CO₂) from the discrete samples (c_WS(CO₂) was used to calibrate the optode output p_c(CO₂) in terms of c(CO₂). c(CO₂) had a better correlation than p(CO₂) (R²=0.77 vs. R²=0.02).

Plotting the regression residuals (c_r(CO₂), calculated as the difference between c_WS(CO₂) and the value predicted by the regression) revealed a quadratic relation between the regression residuals and water temperature (θ). We have therefore included θ and θ² in the optode calibration (Fig. 6a). This second calibration increased the correlation coefficient R² from 0.77 to 0.90 and decreased the standard deviation of the regression residuals from 1.3 to 0.8 µmol kg⁻¹. Even with the explicit inclusion of temperature in the calibration, the CO₂ optode response remained more closely related to c(CO₂) than p(CO₂) (Fig. 6b).

Figure 6 Regression (black lines, reg1) of the CO₂ optode output p_c(CO₂) against (a) co-located concentration c_WS(CO₂) that has an uncertainty of 0.28 µmol kg⁻¹ (b) and partial pressure p_WS(CO₂) of CO₂ in discrete water samples (black dots). Also shown are the values predicted by including θ and θ² in the regression used for optode calibration (red dots, reg2). The regression equations are (a) reg1: c_WS(CO₂)/(µmol kg⁻¹) = (0.033 ± 0.003)p_c(CO₂)/µatm – 1.8 ± 1.6 (R²=0.77); reg2: c_WS(CO₂)/(µmol kg⁻¹) = (0.12 ± 0.14)θ/^∘C – (0.071 ± 0.011)(θ/^∘C)² + (0.0094 ± 0.0048)p_c(CO₂)/µatm + 16 ± 4 (R²=0.90). (b) reg1: p_WS(CO₂)/µatm = (0.05 ± 0.05)p_c(CO₂)/µatm + 344 ± 33 (R²=0.02); reg2: p_WS(CO₂)/µatm = (21 ± 3)θ/^∘C – (1.9 ± 0.2)(θ/^∘C)² + (0.2 ± 0.1)p_c(CO₂)/µatm + 209 ± 76 (R²=0.60).

2.6 Regional algorithm to estimate A_T

To calculate c(DIC)_, we used two variables: (1) glider c(CO₂) derived as described in Sect. 2.5 and (2) A_T derived using a regional algorithm based on S and θ depths of less than 1000 m. The algorithm followed the approach of Lee et al. (2006) and was derived using 663 water samples collected at OWSM from 2004 to 2014 and GLODAPv2 (Olsen et al., 2016) data from the year 2000 in the deployment region. Discrete samples with S<33 were removed because these values were lower than the minimum S measured by the glider. The derived A_T parameterisation is

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{A}_{\mathrm{T},\mathrm{reg}}/\left(\mathrm{\mu }\mathrm{mol}{\mathrm{kg}}^{-\mathrm{1}}\right)& =\mathrm{2317.03}\pm \mathrm{12}+\mathrm{33.12}\pm \mathrm{6.21}(S-\mathrm{35})\\ & +\mathrm{7.94}\pm \mathrm{11.38}(S-\mathrm{35}{)}^{\mathrm{2}}\\ & +\mathrm{0.96}\pm \mathrm{1.79}(\mathit{\theta }/{}^{\circ }\mathrm{C}-\mathrm{20})\\ & +\mathrm{0.01}\pm \mathrm{0.06}(\mathit{\theta }/{}^{\circ }\mathrm{C}-\mathrm{20}{)}^{\mathrm{2}}.\end{array}\end{array}$$

The parameterisation has an uncertainty of 8.2 µmol kg⁻¹ calculated as the standard deviation of the residual difference between actual and parameterised A_T.

To test this parameterisation, we compared the predicted A_T,reg values with discrete measurements (A_T,WS) collected close in terms of time, potential density (σ₀) and distance to the glider transect (n=60). These discrete samples and the glider had mean temperature and salinity differences of (0.17 ± 0.68) ^∘C and 0.03 ± 0.013, respectively. The mean difference between A_T,WS and A_T,reg was (2.1 ± 6.5) µmol kg⁻¹.

This A_T parameterisation was used in CO2SYS (Van Heuven et al., 2011) to calculate c(DIC) from A_T,reg and the calibrated c(CO₂), c_G,cal(CO₂). These calculated c_G,cal(DIC) values were compared with c_WS(DIC) of the same set of discrete samples used to calibrate c_G,cal(CO₂), the only difference being that instead of the actual total alkalinity of the water sample (A_T,WS), we used A_T,reg. The mean difference between c_G,cal(DIC) and c_WS(DIC) was (3 ± 11) µmol kg⁻¹, with the non-zero bias and the standard deviation due to the uncertainties in the A_T,reg parameterisation and the c_G,cal(CO₂) calibration.

2.7 Quality control of other measurement variables

The thermal lag of the glider conductivity sensor was corrected using the method of Gourcuff (2014). Single-point outliers in conductivity were removed and replaced by linear interpolation. The glider CTD salinity was affected by presumed particulate matter stuck in the conductivity cell (Medeot et al., 2011) during dives 147, 234, 244, 251, 272, 279, 303, 320 and 397, and sensor malfunction caused a poor match between glider ascent and descent during dives 214, 215, 235 and 243. These dives were removed from the subsequent analysis.

Glider-reported chlorophyll concentrations, c_raw(Chl a), were computed using the factory coefficients. c_raw(Chl a) was affected by photochemical quenching during the daytime dives. To correct for quenching, we used the method of Hemsley et al. (2015) based on the nighttime relationship between fluorescence and optical backscatter. This relationship was established in the top 60 m and the nighttime values were selected between sunset and sunrise. We calculated a linear fit between c_raw(Chl a) measured at night, c_N(Chl a), and the backscatter signal measured at night (b_N). The slope and the intercept were then used to derive corrected daytime c_D(Chl a). The glider-reported chlorophyll concentration has not been calibrated against in situ samples and is not expected to be accurate, even after correction for quenching. However, it should give an indication of the depth of the deep chlorophyll concentration maximum (z_DCM) and the direction of chlorophyll concentration change (up/down). The 8 d means of c_raw(Chl a) were compared with satellite 8 d composite chlorophyll concentration (Fig. 7) from Ocean Colour CCI (https://esa-oceancolour-cci.org/, last access: 7 May 2020) and gave a mean difference of (0.12 ± 0.08) mg m⁻³.

Figure 7 Comparison between the 8 d glider c(Chl a) (c_G(Chl a)) mean and the 8 d satellite c(Chl a) (c_S(Chl a)) download from Ocean Colour CCI (https://esa-oceancolour-cci.org/, last access: 7 May 2020) as a time series (a) and scatter plot (b).

2.8 Calculation of oxygen-based net community production N(O₂)

Calculating net community production N from glider data is challenging because the glider continuously moves through different water masses. For that reason, we subdivided the transect by binning the data into 0.1^∘ latitude intervals to derive O₂ concentration changes every two transects. The changes were calculated between transects in the same direction of glider travel (e.g. transects 1 and 3, both in the N–S direction) to have approximately the same time difference (40–58 d) at every latitude. If instead we had used two consecutive transects, this would lead to a highly variable time difference of near 0 to about 50 d along the transect.

We calculated N(O₂) (in mmol m⁻² d⁻¹) from the oxygen inventory changes ($\mathrm{\Delta }I\left({\mathrm{O}}_{\mathrm{2}}\right)/\mathrm{\Delta }t$) corrected for air–sea exchange Φ(O₂), normalised to z_mix when z_mix was deeper than the integration depth of z_lim, entrainment E(O₂) and diapycnal eddy diffusion F_v(O₂):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}N\left({\mathrm{O}}_{\mathrm{2}}\right)& ={\displaystyle \frac{\mathrm{\Delta }I\left({\mathrm{O}}_{\mathrm{2}}\right)}{\mathrm{\Delta }t}}+\mathit{\Phi }\left({\mathrm{O}}_{\mathrm{2}}\right){\displaystyle \frac{\min \left(z_{\lim },z_{\mathrm{mix}}\right)}{z_{\mathrm{mix}}}}\\ & -E\left({\mathrm{O}}_{\mathrm{2}}\right)-F_{\mathrm{v}}\left({\mathrm{O}}_{\mathrm{2}}\right).\end{array}\end{array}$$

The inventory changes were calculated as the difference between two transects of the integrated oxygen concentration C(O₂). C(O₂) (in mmol m⁻³) was derived from the oxygen content c(O₂) (in µmol kg⁻¹) by multiplication with the water density (about 1027 kg m⁻³, but we used the actual values). A default integration depth of 45 m was chosen to capture the deepest extent of the deep chlorophyll maximum (z_DCM) found during the deployment, which likely represents the extent of the euphotic zone.

The inventory changes for every latitude bin were calculated using the following equation:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{\Delta }I\left({\mathrm{O}}_{\mathrm{2}}\right)}{\mathrm{\Delta }t}}={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{45}\mathrm{m}}{C}_{n+\mathrm{1}}({\mathrm{O}}_{\mathrm{2}},z)\mathrm{d}z-{\int }_{\mathrm{0}}^{\mathrm{45}\mathrm{m}}{C}_n({\mathrm{O}}_{\mathrm{2}},z)\mathrm{d}z}{t_{n+\mathrm{1}}-t_n}},\end{array}$$

where n is the transect number, t is the day of the year, and C(O₂,z) is the vertical O₂ concentration profile.

The air–sea flux of oxygen, Φ(O₂), was calculated for each glider dive using the median C(O₂), θ and S in the top 10 m. We followed the method of Woolf and Thorpe (1991) that includes the effect of bubble equilibrium supersaturation in the calculations:

$$\begin{array}{}\text{(7)}& \mathit{\Phi }\left({\mathrm{O}}_{\mathrm{2}}\right)=k_{\mathrm{w}}\left({\mathrm{O}}_{\mathrm{2}}\right)\mathit{\left\{}\right(C\left({\mathrm{O}}_{\mathrm{2}}\right)-[\mathrm{1}+{\mathit{\Delta }}_{\mathrm{bub}}({\mathrm{O}}_{\mathrm{2}}\left)\right]C_{\mathrm{sat}}\left({\mathrm{O}}_{\mathrm{2}}\right)\mathit{\}},\end{array}$$

where k_w(O₂) is the gas transfer coefficient, Δ_bub(O₂) is the increase of equilibrium saturation due to bubble injection, and C_sat(O₂) is the oxygen saturation. C_sat(O₂) was calculated from S and θ using the solubility coefficients of Benson and Krause (1984), as fitted by Garcia and Gordon (1992). Δ_bub(O₂) was calculated from the following equation:

$$\begin{array}{}\text{(8)}& {\mathit{\Delta }}_{\mathrm{bub}}\left({\mathrm{O}}_{\mathrm{2}}\right)=\mathrm{0.01}{\left({\displaystyle \frac{U}{U_{\mathrm{0}}}}\right)}^{\mathrm{2}},\end{array}$$

where U is 10 m wind speed with 1 h resolution (ECMWF ERA5; https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5, last access: 1 February 2021) and U₀ represents the wind speed when the oxygen concentration is 1 % supersaturated and has a value of 9 m s⁻¹ (Woolf and Thorpe, 1991). U has a spatial resolution of 0.25^∘ latitude and 0.25^∘ longitude and was interpolated to the glider position at the beginning of the dive.

The transfer velocity k_w(O₂) was calculated based on Wanninkhof (2014):

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{k_{\mathrm{w}}\left({\mathrm{O}}_{\mathrm{2}}\right)}{\mathrm{cm}{\mathrm{h}}^{-\mathrm{1}}}}=\mathrm{0.251}{\left({\displaystyle \frac{\mathit{Sc}\left({\mathrm{O}}_{\mathrm{2}}\right)}{\mathrm{660}}}\right)}^{-\mathrm{0.5}}{\left({\displaystyle \frac{U}{\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}}}\right)}^{\mathrm{2}}.\end{array}$$

The Schmidt number, Sc(O₂), was calculated using the parameterisation of Wanninkhof (2014). To account for wind speed variability, k_w(O₂) applied to calculate N(O₂) was a weighted mean. This value was calculated using the varying daily-mean wind speed U in the time interval between t_n and t_n+1(Δt) (50 d) using a five-point-median z_mix (Sect. 3.2) (Reuer et al., 2007). The time interval is the same as the one used to calculate $\frac{\mathrm{\Delta }I\left({\mathrm{O}}_{\mathrm{2}}\right)}{\mathrm{\Delta }t}$.

The entrainment flux, E(O₂), was calculated as the oxygen flux when the mixed layer depth deepens in time and is greater than z_lim at time t₂:

$$\begin{array}{}\text{(10)}& E\left({\mathrm{O}}_{\mathrm{2}}\right)={\displaystyle \frac{I\left({\mathrm{O}}_{\mathrm{2}},t_{\mathrm{1}},z_{\mathrm{mix}}\left(t_{\mathrm{2}}\right)\right)\frac{z_{\lim }}{z_{\mathrm{mix}\left(t_{\mathrm{2}}\right)}}-I({\mathrm{O}}_{\mathrm{2}},t_{\mathrm{1}},z_{\lim })}{t_{\mathrm{2}}-t_{\mathrm{1}}}},\end{array}$$

where t₂−t₁ represents the change in time, z_mix is the mixed layer depth, I(O${}_{\mathrm{2}},t_{\mathrm{1}},z_{\mathrm{mix}}\left(t_{\mathrm{2}}\right)$) is the expected inventory that would result from a mixed layer deepening to z_mix(t₂) between t₂ and t₁, and I(O${}_{\mathrm{2}},t_{\mathrm{1}},z_{\lim }$) is the original inventory at t₁.

The effect of diapycnal eddy diffusion (F_v) was calculated at z_mix when it was deeper than z_lim and at z_lim when z_mix was shallower than z_lim, using the following equation:

$$\begin{array}{}\text{(11)}& F_{\mathrm{v}}\left({\mathrm{O}}_{\mathrm{2}}\right)=K_z{\displaystyle \frac{\partial C\left({\mathrm{O}}_{\mathrm{2}}\right)}{\partial z}},\end{array}$$

for a vertical eddy diffusivity (K_z) of 10⁻⁵ m² s⁻¹ (Naveira Garabato et al., 2004). The effect of F_v(O₂) on N(O₂) was negligible (Fig. A2b) with a median of (−0.1 ± 0.5) mmol m⁻² d⁻¹.

2.9 Calculation of dissolved inorganic carbon-based net community production, N(DIC)

N(DIC) was expressed in mmol m⁻² d⁻¹ and was calculated from DIC inventory changes ($\mathrm{\Delta }I\left(\mathrm{DIC}\right)/\mathrm{\Delta }t$), air–sea flux of CO₂, Φ(CO₂), entrainment E(DIC) and diapycnal diffusion F_v(DIC):

$$\begin{array}{}\text{(12)}& \begin{array}{rl}N\left(\mathrm{DIC}\right)& =-{\displaystyle \frac{\mathrm{\Delta }\left(\mathrm{DIC}\right)}{\mathrm{\Delta }t}}-\mathit{\Phi }\left({\mathrm{CO}}_{\mathrm{2}}\right){\displaystyle \frac{\min \left(z_{\lim },z_{\mathrm{mix}}\right)}{z_{\mathrm{mix}}}}\\ & +E\left(\mathrm{DIC}\right)+F_{\mathrm{v}}\left(\mathrm{DIC}\right).\end{array}\end{array}$$

Firstly, Φ(CO₂) was calculated using the 10 m wind speed with 1 h resolution downloaded from ECMWF ERA5. As for oxygen, we selected the closest wind speed data point at the beginning of each glider dive. We used the monthly mean atmospheric CO₂ dry mole fraction (x(CO₂)) downloaded from the Greenhouse Gases Reference Network Site (https://www.esrl.noaa.gov/gmd/ccgg/ggrn.php, last access: 28 February 2019) closest to the deployment at Mace Head, County Galway, Ireland (Dlugokencky et al., 2015). Using x(CO₂), we calculated the air-saturation concentration C_atm(CO₂):

$$\begin{array}{}\text{(13)}& C_{\mathrm{atm}}\left({\mathrm{CO}}_{\mathrm{2}}\right)=x\left({\mathrm{CO}}_{\mathrm{2}}\right)p_{\mathrm{baro}}F\left({\mathrm{CO}}_{\mathrm{2}}\right),\end{array}$$

where p_baro is the mean sea level pressure and F(CO₂) is the CO₂ solubility function (in mol dm⁻³ atm⁻¹) calculated from surface θ and S (Weiss and Price, 1980).

The seawater C(CO₂) at the surface was calculated using the median in the top 10 m between the glider ascent and descent of the following dive C(CO₂). From this, Φ(CO₂) was calculated:

$$\begin{array}{}\text{(14)}& \mathit{\Phi }\left({\mathrm{CO}}_{\mathrm{2}}\right)=k\left({\mathrm{CO}}_{\mathrm{2}}\right)\left[C\right({\mathrm{CO}}_{\mathrm{2}})-C_{\mathrm{atm}}({\mathrm{CO}}_{\mathrm{2}}\left)\right].\end{array}$$

k(CO₂) was calculated using the parameterisation of Wanninkhof (2014):

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{k_{\mathrm{w}}\left({\mathrm{CO}}_{\mathrm{2}}\right)}{\mathrm{cm}{\mathrm{h}}^{-\mathrm{1}}}}=\mathrm{0.251}{\left({\displaystyle \frac{\mathit{Sc}\left({\mathrm{CO}}_{\mathrm{2}}\right)}{\mathrm{660}}}\right)}^{-\mathrm{0.5}}{\left({\displaystyle \frac{U}{\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}}}\right)}^{\mathrm{2}}.\end{array}$$

Sc(CO₂) is the dimensionless Schmidt number at the seawater temperature (Wanninkhof, 2014). To account for wind speed variability, k_w(CO₂) applied to calculate N(DIC) was a weighted mean based on the varying daily-mean wind speed U in the time interval between t_n and t_n+1(Δt) used to calculate $\frac{\mathrm{\Delta }I\left(\mathrm{DIC}\right)}{\mathrm{\Delta }t}$ and for 40–50 d to calculate Φ(CO₂) (Sect. 3.2) (Reuer et al., 2007).

The DIC inventory changes were calculated in the top 45 m with the following equation:

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\mathrm{\Delta }I\left(\mathrm{DIC}\right)}{\mathrm{\Delta }t}}={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{45}\mathrm{m}}{C}_{n+\mathrm{1}}(\mathrm{DIC},z)\mathrm{d}z-{\int }_{\mathrm{0}}^{\mathrm{45}\mathrm{m}}{C}_n(\mathrm{DIC},z)\mathrm{d}z}{t_{n+\mathrm{1}}-t_n}}.\end{array}$$

Just as for C(O₂), C(DIC) (in mmol m⁻³) was derived from the DIC content c(DIC) (in µmol kg⁻¹) by multiplication with the water density (about 1027 kg m⁻³, but we used the actual values).

The entrainment flux, E(DIC), was calculated as the DIC flux when the mixed layer depth deepens in time and is greater than z_lim at time t₂:

$$\begin{array}{}\text{(17)}& E\left(\mathrm{DIC}\right)={\displaystyle \frac{I\left(\mathrm{DIC},t_{\mathrm{1}},z_{\mathrm{mix}}\left(t_{\mathrm{2}}\right)\right)\frac{z_{\lim }}{z_{\mathrm{mix}\left(t_{\mathrm{2}}\right)}}-I(\mathrm{DIC},t_{\mathrm{1}},z_{\lim })}{t_{\mathrm{2}}-t_{\mathrm{1}}}}.\end{array}$$

As for oxygen, the effect of diapycnal eddy diffusion (F_v) was calculated at z_mix when it was deeper than z_lim and at z_lim when z_mix was shallower than z_lim, using the following equation:

$$\begin{array}{}\text{(18)}& F_{\mathrm{v}}\left(\mathrm{DIC}\right)=K_z{\displaystyle \frac{\partial C\left(\mathrm{DIC}\right)}{\partial z}},\end{array}$$

for a K_z of 10⁻⁵ m² s⁻¹ (Naveira Garabato et al., 2004). The effect of F_v(DIC) was negligible (Fig. A2a) with a median of (0.1 ± 0.3) mmol m⁻² d⁻¹.

The contribution of horizontal advection to N(DIC) was considered minimal over the timescales where we calculated inventory changes because previous studies have shown that changes in C(DIC) during summer are mainly controlled by biology and air–sea interactions (Gislefoss et al., 1998). For that reason, previous studies that estimated N in the Norwegian Sea have also neglected advective fluxes (Falck and Anderson, 2005; Falck and Gade, 1999; Kivimäe, 2007; Skjelvan et al., 2001).

Uncertainties in N(DIC) and N(O₂) were evaluated with a Monte Carlo approach. The uncertainties of the input variables are shown in Table 2; we repeated the analysis 1000 times. The total uncertainty in N was calculated as the standard deviation of the 1000 Monte Carlo simulations.

Table 2 Uncertainty associated with N(DIC) and N(O₂) input variables calculated by a Monte Carlo approach.

3 Results

The uncorrected p(CO₂) values presented in Fig. 8 were analysed up to dive 400 (24 July 2014). For the following dives, the CO₂ optode stopped sampling in the first 150 m (Fig. 8d). Instead, the uncorrected temperature θ, salinity S, c(O₂) and c_raw(Chl a) were analysed for all the dives (30 October 2014). The raw optode c(O₂) data were calibrated and drift corrected, and c(CO₂) was drift corrected, lag corrected and recalibrated, then used to quantify the temporal and spatial changes in N and Φ together with the quenching corrected c_raw(Chl a) to evaluate net community production changes.

Figure 8 Raw glider data for all 703 dives with latitude of the glider trajectory at the top (black: NwAC; red: NCC, separated by a S of 35). (a) Temperature θ, (b) salinity S, (c) oxygen concentration c_G(O₂), (d) uncorrected CO₂ optode output p_u(CO₂) and (e) chlorophyll a concentration c_raw(Chl a). The white space means that the sensors did not measure any data. The pink line is z_mix calculated using a threshold criterion of Δθ=0.5 ^∘C to the median θ in the top 5 m (Obata et al., 1996; Monterey and Levitus, 1997; Foltz et al., 2003). The dotted black line designates z_lim, used as a depth limit to calculate N. Contoured black lines represent isopycnals.

3.1 O₂ and CO₂ optode calibration

The O₂ optode drift caused a continuous and unexpected decrease of the uncorrected c_G(O₂) from 290 to 282 µmol kg⁻¹ for σ₀ > 1028 kg m⁻³ (Fig. 8c). The ratio $c_{\mathrm{C}}\left({\mathrm{O}}_{\mathrm{2}}\right)/c_{\mathrm{G}}\left({\mathrm{O}}_{\mathrm{2}}\right)$ against day of the year used for the drift correction had a good correlation with time (R²=0.90), showing a continuous increase of 0.0004 d⁻¹ (Fig. 3), equivalent to a decrease in the measured glider O₂ concentration of 0.11 µmol kg⁻¹ d⁻¹. It was possible to apply the correction because c_C(O₂) had low temporal variability for the chosen potential density σ₀ > 1028 kg m⁻³. The c_C(O₂) values from OWSM and GLODAPv2 had a mean of (305 ± 3) µmol kg⁻¹, varying from 294 to 315 µmol kg⁻¹ (Fig. A1). The drift correction reduced the variability of c_G(O₂) in the selected potential density range from a standard deviation of 7.3 µmol kg⁻¹ to a standard deviation of 2.4 µmol kg⁻¹ (Fig. 9a).

Figure 9 (a) c(O₂); (b) s(O₂) = $c\left({\mathrm{O}}_{\mathrm{2}}\right)/c_{\mathrm{sat}}\left({\mathrm{O}}_{\mathrm{2}}\right)$ with z_DCM (red line), z_mix (pink line) five-point median z_mix (dotted pink line). The black line indicates σ₀=1028 kg m⁻³. The top panel indicates glider latitude (black: NwAC; red: NCC).

Following drift, lag and scale corrections, glider fugacity f_G(CO₂) derived from Eq. (2) had a mean difference of (2 ± 22) µatm to the discrete samples (n=55; not shown) and c(DIC) had a mean difference of (3 ± 11) µmol kg⁻¹ (Fig. 10). p(CO₂) and f(CO₂) are almost identical, but f(CO₂) takes into account the non-ideal nature of the gas phase. The optode was able to capture the temporal and spatial variability showing that NCC had a lower DIC concentration than NwAC. Restricting the f(CO₂) comparison to the discrete samples in the top 10 m gave a mean difference of (19 ± 31) µatm (n=6). We also compared glider f_G(CO₂) with the Surface Ocean CO₂ Atlas (SOCAT) f(CO₂) (Bakker et al., 2016) data in the region during the deployment (Fig. 11). During the whole deployment, there was general agreement between f_G(CO₂) and f_SOCAT(CO₂). f_G(CO₂) varied between 204 and 391 µatm, while f_SOCAT(CO₂) varied between 202 and 428 µatm (Fig. 11).

Our results are in agreement with Jeansson et al. (2011), who found the surface NCC was the region with the lowest c(DIC) values (2083 µmol kg⁻¹) in the Norwegian Sea. This was confirmed during our deployment because c(DIC) was (2081 ± 39) µmol kg⁻¹ in the NCC region and (2146 ± 27) µmol kg⁻¹ in the NwAC region (Fig. 10) and c(O₂) was > 300 µmol kg⁻¹ in the NwAC and < 280 µmol kg⁻¹ in the NCC.

Figure 10 c(DIC) contour plot with z_DCM (red line), z_mix (pink line) five-point median z_mix (dotted pink line). The black line indicates σ₀ = 1028 kg m⁻³. The top panel indicates glider latitude (black: NwAC; red: NCC).

Figure 11 Comparison between surface f(CO₂) from 2014 SOCAT and CO₂ optode on the glider. (a) The black lines are the median glider f(CO₂) in the top 10 m, with f_c(CO₂) (dotted line) corresponding to regression 1 (Fig. 6a) and f_t(CO₂) (continuous line) to regression 2 (Fig. 6a). Discrete samples collected during the deployment are shown as red dots, with the other coloured dots representing cruises in the SOCAT database (Bakker et al., 2016). (b) Glider and SOCAT data positions (same colours as in panel a).

3.2 Air–sea exchange

The surface water was supersaturated with oxygen all summer (Fig. 12). From May, this supersaturation drove a continuous O₂ flux from the sea to the atmosphere. However, the flux varied throughout the deployment, having a median of 25 mmol m⁻² d⁻¹ (5th percentile: −31 mmol m⁻² d⁻¹; 95th percentile: 88 mmol m⁻² d⁻¹). Prior to the spring period of increased Chl a inventory, the supersaturation varied between 0 to 10 µmol kg⁻¹. Φ(O₂) had a median of −1.4 mmol m⁻² d⁻¹ (5th percentile: −49 mmol m⁻² d⁻¹; 95th percentile: 23 mmol m⁻² d⁻¹). Then, during the spring period of increased Chl a inventory, the surface concentration increased by over 35 µmol kg⁻¹, causing a peak in Φ(O₂) of 140 mmol m⁻² d⁻¹. A second period of increased Chl a inventory was encountered in June and had a larger Φ(O₂) up to 118 mmol m⁻² d⁻¹, driven by supersaturation of 68 µmol kg⁻¹. The fluxes were smaller than during spring and were associated with an increase of c_raw(Chl a) from 2.5 mg m⁻³ to the summer maximum of 4.0 mg m⁻³. However, prior to the increased spring Chl a inventory, Φ(O₂) showed a few days of influx into seawater caused by a decrease of θ from 7.6 to 5.9 ^∘C that increased C_sat(O₂). The influx at the beginning of the deployment is partly due to the Δ_bub(O₂) correction that resulted in $[\mathrm{1}+{\mathit{\Delta }}_{\mathrm{bub}}({\mathrm{O}}_{\mathrm{2}}\left)\right]c_{\mathrm{sat}}\left({\mathrm{O}}_{\mathrm{2}}\right)>c\left({\mathrm{O}}_{\mathrm{2}}\right)$ for U > 10 m s⁻¹. In August, the surface supersaturation decreased to 2.3 µmol kg⁻¹ and Φ(O₂) decreased to a monthly minimum of −7.6 mmol m⁻² d⁻¹. In the second half of September, the surface water became undersaturated by −2.6 µmol kg⁻¹, causing O₂ uptake with a median flux of −13 mmol m⁻² d⁻¹ (5th percentile: −39 mmol m⁻² d⁻¹; 95th percentile: 10 mmol m⁻² d⁻¹).

Figure 12 Air–sea flux of O₂ and CO₂ during spring and summer for CO₂ and during spring, summer and autumn for O₂, (a) c_sat(O₂) in blue and c(O₂) in red, (b) c_sat(CO₂) in blue and c(CO₂) in red, (c) Δc(O₂) = c(O₂)−c_sat(O₂), (d) Δc(CO₂) = c(CO₂)−c_sat(CO₂), (e) sea surface temperature, (f) k_w(O₂) in blue and k_w(CO₂) in red normalised back to 50 d (Reuer et al., 2007), (g) oxygen air–sea flux Φ(O₂) and (h) CO₂ air–sea flux Φ(CO₂). The flux from sea to air is positive, while that from air to sea is negative.

The CO₂ flux from March to July was always from the air to the sea (Fig. 12), with a median of −5.2 mmol m⁻² d⁻¹ (5th percentile: −14 mmol m⁻² d⁻¹; 95th percentile: −1.5 mmol m⁻² d⁻¹). An opposite flux direction is expected for Φ(O₂) and Φ(CO₂) during the productive season when net community production is the main driver of concentration changes. After the summer period of increased Chl a inventory, the flux had a median of −11 mmol m⁻² d⁻¹ (5th percentile: −16 mmol m⁻² d⁻¹; 95th percentile: −6.8 mmol m⁻² d⁻¹), in agreement with previous studies that classified the Norwegian Sea as a CO₂ sink (Takahashi et al., 2002; Skjelvan et al., 2005). Φ(CO₂) for the discrete samples from 18 March to 14 June (n=13) varied from 0.1 to −13 mmol m⁻² d⁻¹.

3.3 N(O₂)

To capture the entire euphotic zone, we calculated N(O₂) and N(DIC) using an integration depth of z_lim=45 m because the mean deep chlorophyll maximum (DCM) depth was z_DCM = (20 ± 18 m) (Fig. 9). For comparison, the mixed layer depth was deeper, varied more strongly and had a mean value of z_mix = (68 ± 78) m, using a threshold criterion of Δθ=0.5 ^∘C to the median θ value in the top 5 m of the glider profile (Obata et al., 1996; Monterey and Levitus, 1997; Foltz et al., 2003).

The two N values were calculated as the difference in inventory changes between two transects when the glider moved in the same direction.

During the deployment, we sampled two periods of increased Chl a inventory, the first one in May and a second one in June. The chlorophyll a inventory ($I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$) was calculated integrating c_raw(Chl a) to z_lim. To remove outliers, we used a five-point moving mean of $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$.

The N(O₂) changes were dominated by Φ(O₂), which had an absolute median of 34 mmol m⁻² d⁻¹ (5th percentile: 4.3 mmol m⁻² d⁻¹; 95th percentile: 86 mmol m⁻² d⁻¹), followed by I(O₂), which had a median of 15 mmol m⁻² d⁻¹ (5th percentile: 2.3 mmol m⁻² d⁻¹; 95th percentile: 29 mmol m⁻² d⁻¹), F_v(O₂), which had an absolute median of 0.3 mmol m⁻² d⁻¹ (5th percentile: 0 mmol m⁻² d⁻¹; 95th percentile: 1.0 mmol m⁻² d⁻¹) and E(O₂), which had a median of 0 mmol m⁻² d⁻¹ (5th percentile: −1.2 mmol m⁻² d⁻¹; 95th percentile: 0 mmol m⁻² d⁻¹).

At the beginning of May, $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$ increased to 97 mg m⁻² and N(O₂) = (95 ± 16) mmol m⁻² d⁻¹. After this period, $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$ decreased to 49 mg m⁻² and N(O₂) = (−4.6 ± 1.6) mmol m⁻² d⁻¹. During the summer $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$ increased to 110 mg m⁻², which caused a sharp increase of N(O₂) to (126 ± 25) mmol m⁻² d⁻¹. $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$ remained higher than 50 mg m⁻² until the end of June when N(O₂) was (31 ± 9) mmol m⁻² d⁻¹. The passage of the glider from NwAC to NCC was accompanied by a drop of surface c(O₂) from 330 to 280 µmol kg⁻¹ (Fig. 9) that resulted in lower Φ(O₂) and N(O₂) values (Fig. 13). At the same time, $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$ decreased to 35 mg m⁻² showing that the decrease of N(O₂) depended on the passage to NCC and a decrease of biological production. After the beginning of August, $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$ decreased to 49 mg m⁻² and N(O₂) turned negative with a minimum of (−23 ± 25) mmol m⁻² d⁻¹. In October, during the last glider transect, $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$ continued decreasing to 27 mg m⁻² leading to the minimum N(O₂) of (−52 ± 11) mmol m⁻² d⁻¹.

Integrating N(O₂) from March to October gives a flux of 4.9 mol m⁻² a⁻¹ (Table 3; discussed in Sect. 4.2).

Figure 13 (a) Components of the N(O₂) calculation: $\mathrm{\Delta }I\left({\mathrm{O}}_{\mathrm{2}}\right)/\mathrm{\Delta }t$ (red), E(O₂) (blue), F_v(O₂) (violet), Φ(O₂) (green) with k_w(O₂) weighted over 50 d, N(O₂) (yellow). (b) Chl a inventory in the top 45 m, $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$ (violet). Chl a inventory for the whole water column, $I_{\mathrm{raw},z_{\mathrm{all}}}\left(\mathrm{Chl}a\right)$ (dotted violet line). The vertical black lines represent each glider transect. Between the two vertical red lines, the glider was in the NCC region.

3.4 N(DIC)

In the case of N(DIC), the main drivers were the inventory changes with an absolute median of 29 mmol m⁻² d⁻¹ (5th percentile: 1.3 mmol m⁻² d⁻¹; 95th percentile: 57 mmol m⁻² d⁻¹), followed by Φ(CO₂), which had an absolute median of 7.0 mmol m⁻² d⁻¹ (5th percentile: 0.8 mmol m⁻² d⁻¹; 95th percentile: 15 mmol m⁻² d⁻¹), F_v(DIC), which had an absolute median of 0.2 mmol m⁻² d⁻¹ (5th percentile: 0 mmol m⁻² d⁻¹; 95th percentile: 1.3 mmol m⁻² d⁻¹) and E(DIC), which had a median of 0 mmol m⁻² d⁻¹ (5th percentile: 0 mmol m⁻² d⁻¹; 95th percentile: 3.4 mmol m⁻² d⁻¹). During the period of increased Chl a inventory N(DIC) was (21 ± 4.5) mmol m⁻² d⁻¹. Later, $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$ decreased to 30 mg m⁻² driving N(DIC) to negative values with a minimum of (−2.7 ± 5.0) mmol m⁻² d⁻¹. In the next transect, the glider measured the maximum $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$ of 111 mg m⁻² that increased N(DIC) to (85 ± 4.5) mmol m⁻² d⁻¹. This maximum was reached during a transect when the glider moved in NCC, which had a c(DIC) of 2080 µmol kg⁻¹ at the surface compared with the 2150 µmol kg⁻¹ in NwAC and drove a continuous positive N(DIC), which had a minimum of (36 ± 7.4) mmol m⁻² d⁻¹ (Fig. 14).

Integrating N(DIC) from March to July gives a flux of 3.3 mol m⁻² a⁻¹ (Table 3; discussed in Sect. 4.2).

Table 3 Net community production (N) estimates in the Norwegian Sea (with integration depth z_lim). Falck and Anderson (2005) used year-round data from 1960 to 2000 between 62 and 70^∘ N and from 1991 to 1994 at OWSM. Skjelvan et al. (2001) used year-round data from 1957 to 1970 and from 1991 to 1998 between 67.5^∘ N, 9^∘ E, and 71.5^∘ N, 1^∘ E, and along 74.5^∘ N from 7 to 15^∘ E. Kivimäe (2007) used year-round data from 1955 to 2005, and Falck and Gade (1999) used year-round data from 1955 to 1988 in all of the Norwegian Sea. While the previous studies report annual N estimates, the present study derives N(O₂) between March and October and N(DIC) between March and July.

Figure 14 (a) Components of the N(DIC) calculation: $\mathrm{\Delta }I\left(\mathrm{DIC}\right)/\mathrm{\Delta }t$ (red), E(DIC) (blue), F_v(CO₂) (violet), Φ(CO₂) (green) with k_w(CO₂) weighted over 50 d, N(DIC) (yellow). (b) Chl a inventory in the top 45 m, $I_{\mathrm{raw},z_{\lim }}\left(\mathrm{Chl}a\right)$ (violet). Chl a inventory for the whole water column, $I_{\mathrm{raw},z_{\mathrm{all}}}\left(\mathrm{Chl}a\right)$ (dotted violet line). The vertical black lines represent each glider transect. Between the two vertical red lines, the glider was in the NCC region.

4 Discussion

4.1 Sensor performance

This study presents data from the first glider deployment with a CO₂ optode. The initial uncalibrated p_u(CO₂) measured by the CO₂ optode had a median of 604 µatm (5th percentile: 566 µatm; 95th percentile: 768 µatm), whereas the p(CO₂) of discrete samples varied from 302 to 421 µatm.

We applied corrections for drift (using deep-water samples as a reference point), sensor lag and calibrated the CO₂ optode against co-located discrete samples throughout the water column.

Atamanchuk et al. (2014) reported that the sensor was affected by a lag that varied from 45 to 264 s depending on temperature. These values were determined in an actively stirred beaker. However, in this study, the sensor was mounted on a glider and was not actively pumped, which increased the response time to 23 min (25th quartile: 18 min; 75th quartile: 30 min). Also, the optode was affected by a continuous drift from 637 to 5500 µatm that was larger than the drift found by Atamanchuk et al. (2015a) which increased by 75 µatm after 7 months.

In this study, the drift- and lag-corrected sensor output showed a better correlation with the CO₂ concentration c(CO₂) than with p(CO₂). The latter two quantities are related to each other by the solubility that varies with θ and S (Weiss, 1974) (Eq. 2). The better correlation with c(CO₂) was probably related to an inadequate temperature parameterisation of the sensor calibration function. Including both temperature and temperature squared in the calibration gave a better fit with c(CO₂) than with p(CO₂) but overall still a lower calibration residual for the former. The sensor output depends on the changes in pH that are directly related to the changes of c(CO₂) in the membrane and – indirectly – p(CO₂), via Henry's law. The calibration is supposed to correct for the temperature dependence of the sensor output (Atamanchuk et al., 2014). So the fact that the sensor output correlated better with c(CO₂) than p(CO₂) is perhaps due to a fortuitous cancellation of an inadequate temperature parameterisation and Henry's law relationship between c(CO₂) and p(CO₂).

The calibrated optode output captured the c(DIC) changes in space and time with a standard deviation of 11 µmol kg⁻¹ compared with the discrete samples. c(DIC) decreased from 2130 to 2000 µmol kg⁻¹ and increased with depth to 2170 µmol kg⁻¹. This shows the potential of the sensor for future studies that aim to analyse the carbon cycle using a high-resolution dataset.

The optode-derived CO₂ fugacity f_G(CO₂) had a mean bias of (1.8 ± 22) µatm compared with the discrete samples. These values are comparable with a previous study when the CO₂ optode was tested for 65 d on a wave-powered profiling cRAWLER (PRAWLER) from 3 to 80 m (Chu et al., 2020), which had an uncertainty between 35 and 72 µatm. The PRAWLER optode was affected by a continuous drift of 5.5 µatm d⁻¹ corrected using a regional empirical algorithm that uses c(O₂), θ, S and σ_o to estimate A_T and c(DIC).

4.2 Norwegian Sea net community production

Increases in N(O₂) and N(DIC) were associated with increases in depth-integrated c_raw(Chl a), designated as periods of increased Chl a inventory I_raw(Chl a), at the beginning of May and in June. During May, I_raw(Chl a) reached 135 mg m⁻². In June, I_raw(Chl a) reached again 135 mg m⁻². Between these two periods, N(DIC) briefly turned negative, indicating remineralisation of the high Chl a inventory material during this period. The period of increased Chl a inventory coincided with a surface temperature increase from 7 to 11 ^∘C and shoaling of the mixed layer from 200 to 20 m. c(O₂) reached a summer maximum of 340 µmol kg⁻¹ and c(DIC) decreased to a summer minimum of 1990 µmol kg⁻¹. In both cases, the main components of the N changes were the inventory and air–sea flux, while the smallest driver was the entrainment. Also, the glider sampled two different water masses characterised by different c(DIC) and c(O₂). This might be the cause of the smaller values N(O₂) and higher values N(DIC) in June and July in NCC compared to NwAC (Figs. 13 and 14). Another explanation might be a consumption of O₂ due to remineralisation and a delay in the response of the c(DIC) that was lowered during the two blooms. A fully functional CO₂ optode in the second part of the deployment would have helped to uncover the cause of the higher N(DIC) than of N(O₂).

Table 3 shows estimates of net community production (N) in the Norwegian Sea. All other studies used ships to gather observations. The estimated N in the four other studies varied from 2.6 to 11.1 mol m⁻² a⁻¹ for N(O₂) and was 3.4 for N(DIC). In our glider study, we obtained between March and July a N(DIC) of 3.3 mol m⁻² a⁻¹ and a N(O₂) of 4.2 mol m⁻² a⁻¹, in agreement with these studies. The ratio of N(O₂) and N(DIC) for an integration depth of 45 m gave a photosynthetic quotient (PQ) of 1.3, in agreement with the Redfield ratio of 1.45 ± 0.15 (Redfield, 1963; Anderson, 1995; Anderson and Sarmiento, 1994; Laws, 1991). The N(O₂) estimate is influenced primarily by the air–sea exchange flux Φ(O₂) (median: 34 mmol m⁻² d⁻¹), followed by the inventory change (15 mmol m⁻² d⁻¹). In contrast, N(DIC) is dominated by the inventory change (−29 mmol m⁻² d⁻¹), followed by Φ(CO₂) (−7.0 mmol m⁻² d⁻¹). This reflects the slower gas-exchange time constant of CO₂ compared with O₂, due to DIC buffering. To compare our results with previous studies, we also used z_lim=30 m (Falck and Gade, 1999) and 100 m (Falck and Anderson, 2005; Kivimäe, 2007). The calculated N(DIC; 30 m) was 3.1 mol m⁻² a⁻¹, N(DIC; 100 m) was 3.4 mol m⁻² a⁻¹, N(O₂; 30 m) was 4.1 mol m⁻² a⁻¹ and N(O₂; 100 m) was 3.7 mol m⁻² a⁻¹. The N(DIC; 100 m) value is in agreement with the value of 3.4 mol m⁻² a⁻¹ given by Falck and Anderson (2005). However, the latter estimate was for the entire year, whereas our estimate only covers the months from March to July. N(O₂) was similar for z_lim=30 and 45 m but lower for z_lim=100 m because of O₂ consumption during organic matter remineralisation below the euphotic zone. The PQ value at 30 was 1.3 and at 100 m decreased to 1.1. Extending N(O₂) to October increased N(O₂; 30 m) and N(O₂; 45 m) to 5.0 and 4.9 mol m⁻² a⁻¹, respectively. Instead, N(O₂; 100 m) decreased to 3.6 mol m⁻² a⁻¹, confirming the consumption of O₂ below the euphotic zone. The calculated N(O₂) until October was in agreement with the previous studies that varied between 2.6 and 11 mol m⁻² a⁻¹.

Some of the previous N(DIC) estimates derived c(DIC) from other variables such as c(O₂), c(PO${}_{\mathrm{4}}^{\mathrm{3}-}$), c(NO${}_{\mathrm{3}}^{-}$), assuming Redfield ratios P : N : C : O₂ 1 : 16 : 106 : −138 (Redfield, 1963). During photosynthesis c(PO${}_{\mathrm{4}}^{\mathrm{3}-}$) and c(NO${}_{\mathrm{3}}^{-}$) are taken up by phytoplankton to form organic matter and are released again after remineralisation of the organic matter giving an indication of NCP changes. Our N(DIC) estimate was 3.3 mol m⁻² a⁻¹ and is similar to 3.4 mol m⁻² a⁻¹ estimated by Falck and Anderson (2005) who used c(DIC) samples directly. The carbon $/$ nutrient ratios vary between water masses and during photosynthesis (Thomas et al., 1999; Copin-Montégut, 2000; Osterroht and Thomas, 2000; Körtzinger et al., 2001).

The difference of the annual N(O₂) and N(DIC) with the previous studies can also be caused by the yearly variability of N in the Norwegian Sea. In fact, Kivimäe (2007) saw an annual variability of N(O₂) from 1955 to 2005 between 4.7 and 18.3 mol m⁻² a⁻¹. In order to understand what is causing these interannual changes, it is important capture inventory and air–sea changes. Also, this study showed that the Norwegian Sea spring, summer and autumn N is strongly affected by the time of sampling. For that reason, N estimated from low-resolution datasets makes the result strongly dependent on the time of sampling. To quantify this interannual variability in N, more high-resolution studies are needed.

5 Conclusions

To the best of our knowledge, this study represents the first glider deployment of a CO₂ optode. The CO₂ optode together with a O₂ optode shows the potential of using these sensors on autonomous observing platforms like Seagliders to quantify the interactions between biogeochemical processes and the marine carbonate system at high spatiotemporal resolution. The deployment helped to uncover NCP and air–sea flux variability over a period of 8 months.

Despite all the problems (drift, lag and poor calibration), the CO₂ optode data could be used to quantify dissolved inorganic carbon concentration variations. The temporal sampling resolution was 106 s in the top 100 m (increasing to 381 s from 500 to 1000 m). This could be improved to less than 10 s, but this would reduce the length of the deployment due to the limited glider battery capacity. With better calibration and stability improvements, the CO₂ optode could be routinely used to measure the carbonate system on gliders, floats and surface vehicles. Glider deployments up to 8 months are possible thanks to the sensor's low power consumption of 8 mW at 5 s sampling intervals and 7 mW at 60 s sampling intervals (Atamanchuk et al., 2014). Combined with other novel sensors that measure another DIC-related quantity such as A_T or c(DIC), CO₂ optodes on gliders could help provide estimates of NCP, air–sea flux, respiration and remineralisation and aragonite saturation.

During our deployment, we calculated O₂ and DIC-based NCP over the spring and summer period. In the future, extended deployments could be used to estimate annual (full-year) NCP. To have an accurate estimate of annual NCP, at least one additional glider deployment is needed to have continuous coverage (Binetti et al., 2020). Similar deployments can be used in other areas of the globe to fill gaps in N(DIC) and N(O₂). In particular, glider deployments have potential in undersampled areas of the globe such as the Southern Ocean and the Arctic. Also, they can be used in well-studied areas such as North Sea and Mediterranean Sea to reduce monitoring costs and compare NCP estimates with previous studies that used other sampling strategies.

Appendix A

Figure A1 Discrete samples c_C(O₂) (yellow), raw glider oxygen c_G(O₂) (blue) and drift-corrected glider oxygen c_G,cal(O₂) (red) for a potential density > 1028 kg m⁻³ at depths less than 1000 m.

Figure A2 Diapycnal mixing (F_v) calculated for the glider descent and ascent for (a) c(DIC) and (b) O₂ at the mixed layer depth (z_mix) when deeper than 45 m (z_lim) and at z_lim when z_mix was shallower than 45 m. In the calculations, we used a vertical eddy diffusivity (K_z) of 10⁻⁵ m² s⁻¹ (Naveira Garabato et al., 2004).

Appendix B: List of symbols

Symbols	Quantity (unit)
AT	total alkalinity (µmol kg−1)
b	backscatter signal (engineering units)
c	amount content (µmol kg−1)
C	amount concentration (mmol m−3)
Chl a	chlorophyll a
DIC	dissolved inorganic carbon
E	entrainment flux (mmol m−2 d−1)
Fv	diapycnal eddy diffusion flux (mmol m−2 d−1)
f(CO2)	fugacity of CO2 (µatm)
I	inventory (mmol m−2)
Kz	diapycnal eddy diffusivity (m2 s−1)
N	net community production (mmol m−2 d−1)
p(CO2)	partial pressure of CO2 (µatm)
S	practical salinity ()
t	time (s)
U	wind speed (m s−1)
x	dry mole fraction (mol mol−1)
zDCM	depth of the deep chlorophyll maximum (m)
zlim	integration depth (m)
zmix	mixed layer depth (m)
Φ	air–sea flux (mmol m−2 d−1)
φ	CO2 optode CalPhase (∘)
σ0	potential density (kg m−3)
θ	Celsius temperature (∘C)
τ	response time (s)