Simultaneous assessment of oxygen- and nitrate-based net community production in a temperate shelf sea from a single ocean glider

The continental shelf seas are important at a global scale for ecosystem services. These highly dynamic regions are under a wide range of stresses, and as such future management requires appropriate monitoring measures. A key metric to understanding and predicting future change are the rates of biological production. We present here the use of an autonomous underwater glider with an oxygen (O₂) and a wet-chemical microfluidic total oxidised nitrogen (${{\mathrm{NO}}_x}^{-}={{\mathrm{NO}}_{\mathrm{3}}}^{-}+{{\mathrm{NO}}_{\mathrm{2}}}^{-}$) sensor during a spring bloom as part of a 2019 pilot autonomous shelf sea monitoring study. We find exceptionally high rates of net community production using both O₂ and ${{\mathrm{NO}}_x}^{-}$ water column inventory changes, corrected for air–sea gas exchange in case of O₂. We compare these rates with 2007 and 2008 mooring observations finding similar rates of ${{\mathrm{NO}}_x}^{-}$ consumption. With these complementary methods we determine the O₂:N amount ratio of the newly produced organic matter (7.8 ± 0.4) and the overall O₂:N ratio for the total water column (5.7 ± 0.4). The former is close to the canonical Redfield O₂:N ratio of 8.6 ± 1.0, whereas the latter may be explained by a combination of new organic matter production and preferential remineralisation of more reduced organic matter at a higher O₂:N ratio below the euphotic zone.

1 Introduction

The coastal shelf seas are a vitally important human resource for numerous ecosystem services, including food, carbon storage, biodiversity, energy, and livelihoods (Halpern et al., 2015). They have a disproportionately large impact, relative to their surface area, on global carbon cycling (Thomas, 2004) as shelf seas provide 10 %–30 % of all marine primary production while comprising less than 10 % of the ocean surface (Harris et al., 2014; Sharples et al., 2019). Given their global role in carbon cycling, understanding the mechanisms driving primary production is vital for predicting how shelf systems will respond to and influence climate change (Palevsky et al., 2013; Legge et al., 2020).

The North Sea is a semi-enclosed region of the northwest European continental shelf between the island of Great Britain in the west and Norway, Denmark, Germany, and the Netherlands in the east and south. The northern boundary is facing the North Atlantic. It can be divided into two energetically distinct regimes: a seasonally stratifying northern part and a fully mixed southern and coastal part (Emeis et al., 2015). The shallower fully mixed region is believed to act as a net sink for CO₂ during the spring bloom periods and a source of CO₂ during the rest of the year (Emeis et al., 2015; Hull et al., 2016). Thomas et al. (2005) suggest that the seasonally stratified regions are a sink for carbon throughout the year. A more recent synthesis of Kitidis et al. (2019) indicates that the seasonally stratified regions were a source of CO₂ during December and January in 2015.

The North Sea is one of the most well-studied marine environments. However, despite the availability of large volumes of data, making robust statements on the long-term changes, natural variability, or predicted process responses to external forcing continues to be challenging (Emeis et al., 2015). As is typical for shelf regions, the North Sea displays high temporal and spatial variability, which makes the analysis of trends more challenging compared with the open ocean (Bozec et al., 2006; Bauer et al., 2013). At present, resolving the nutrient and oxygen dynamics in shelf sea regions remains a challenge for numerical biogeochemical models (Bauer et al., 2013; Emeis et al., 2015). Further development of model parametrisations for these processes, in addition to validation of the model outputs over appropriate spatial and temporal scales, is required.

The North Sea is changing, and the processes controlling the carbon, oxygen, and nutrient dynamics are highly complex and interconnected (Wakelin et al., 2020). Key among those processes is marine production. Estimates of net community production (NCP) can quantify the uptake of CO₂ from the atmosphere into the sea (Alkire et al., 2012; Hull et al., 2016; Plant et al., 2016). Changes to community production can affect how organic material is exported to and respired within the non-ventilated waters below the thermocline (Große et al., 2016). NCP is therefore a useful metric for understanding carbon cycling and an important control on bottom water oxygen concentration.

1.1 The Dogger Bank

The Dogger Bank is a large sandbank within the central North Sea (Fig. 1). Comprising predominantly fine sand and mud, this bank rises 20 m above the surrounding sea floor and is 15 m below sea level in its shallowest parts. The Dogger Bank sits south of a tidal front, which is located in an arc across the central North Sea from Yorkshire in the United Kingdom to the Frisian coast in the Netherlands. This hydrological divide separates well-mixed waters to the south from the seasonally stratified waters of the central North Sea to the north (Emeis et al., 2015). The combination of tidal, wind, and buoyancy forcing tends to develop a general anti-clockwise circulation within the North Sea (Mathis et al., 2015). Variability in the circulation is determined mainly by variation in local winds, which are predominantly from the west. Tidal forcings are dominated by the M2 constituent, with mean tidal excursion in the region of 1.5 km (Fig. 1). Prevailing currents move from the west along the tidal front towards the Dutch coast. These currents are supplied with water from the Scottish coastal current, itself primarily fed from the North Atlantic inflow through the Orkney–Shetland gap (Hill et al., 2008). The Dogger Bank experiences minimal riverine influence, such that maximal nutrient concentrations typically reflect those of the in-flowing Atlantic waters (Greenwood et al., 2010). The changes in N:P ratio observed in more coastal regions of the North Sea are not seen at Dogger Bank, and nitrogen remains the primary limiting nutrient (Burson et al., 2016).

The shallow waters above Dogger Bank are not seasonally stratified as wind stress and tidal currents regularly homogenise the shallowest areas, which replenishes nutrients into the photic zone throughout the year (Riegman et al., 1990). In this paper we focus on the region north of the Dogger Bank where mixed waters converge with the seasonally thermally stratified waters. The exchange of nutrients and phytoplankton at this transition promotes enhanced primary production along the northern edge of the bank, which acts as a hotspot for marine life and supports several important fisheries (Plumeridge and Roberts, 2017). Since 2017, the bank has been designated as a Special Area of Conservation (SAC) by the UK government, signifying it as a region containing habitats and species which are considered most in need of conservation and legal protection. The area on and around Dogger Bank is thought to exhibit year-round phytoplankton production (Nielsen et al., 1993). Due to the shallow depth of Dogger Bank, the spring phytoplankton bloom can be initiated there months before stratification triggers the bloom within the more northern stratified regions of the North Sea (Nielsen et al., 1993). The region north of Dogger Bank typically sees a spring phytoplankton bloom commencing at the start of April.

1.2 The AlterEco project

The glider-based oxygen and total oxidised nitrogen (${{\mathrm{NO}}_x}^{-}={{\mathrm{NO}}_{\mathrm{3}}}^{-}+{{\mathrm{NO}}_{\mathrm{2}}}^{-}$) concentration measurements presented within this paper form part of the larger AlterEco (“An alternative framework to assess marine ecosystem functioning in shelf seas”; https://altereco.ac.uk, last access: 29 September 2021) project. AlterEco is a pilot study of a novel monitoring framework to deliver improved spatiotemporal understanding of key shelf sea ecosystem drivers though the use of autonomous systems, primarily underwater gliders, as part of the Marine Integrated Autonomous Observing Systems programme of the Natural Environment Research Council (NERC). Buoyancy gliders are capable of directed, long-term continuous monitoring over a variable spatial area while covering the full water column within a shelf sea (Alkire et al., 2012; Wihsgott et al., 2019). Gliders have previously been fitted with UV-fluorescence-based ${{\mathrm{NO}}_x}^{-}$ sensors (Johnson et al., 2006). In the central North Sea potentially low concentrations of ${{\mathrm{NO}}_x}^{-}$ together with interferences from organic matter are a challenge for these instruments, which are typically more suitable for observing high concentrations in non-optically complex waters. Recent advances in wet-chemical microfluidic instruments present an opportunity to gather high-resolution and high-quality glider-based measurements. We use these to resolve the ${{\mathrm{NO}}_x}^{-}$ and oxygen fluxes and determine the net community production associated with the spring bloom simultaneously with two complementary methods.

Figure 1 Seaglider path for the AlterEco AE7 2019 mission and the surrounding study area. Timestamps indicate when the glider returned to the western side of the transect. Mean tidal ellipse shown in red (not to scale). Red triangle marks the location of the Centre for Environment, Fisheries and Aquaculture Science (Cefas) North Dogger SmartBuoy, present between 2007 and 2008. Bathymetry data from GEBCO (2019) on a 15 arcsec grid.

2 Data collection and method

2.1 Study area

Observations were made along an east–west transect (Fig. 1) north of the Dogger Bank between 29 March and 25 April 2019 with an underwater glider (Seaglider sg602 “Scapa”). This period was the seventh group of AlterEco glider deployments and is therefore referred to by the mission code “AE7”. The Seaglider took 5 d to move from one end of the transect to the other and completed six transits. The glider was fitted with a standard non-pumped Seabird CT sail, an Aanderaa 4831 oxygen optode (Bittig et al., 2018), and a National Oceanographic Centre (NOC) lab-on-a-chip (LoC) microfluidic ${{\mathrm{NO}}_x}^{-}$ analyser (Beaton et al., 2012; Vincent et al., 2018).

The Seaglider conductivity, temperature (CT), and oxygen data were processed and quality-assured using the University of East Anglia (UEA) Seaglider toolbox (Queste, 2013). This comprised optimising the Seaglider flight parameters to determine speed though water, which was used to apply a cell thermal mass correction for the CT sensor. This is an iterative process as determining glider flight requires determining the water density, which is calculated from the measured salinity (Garau et al., 2011). Some salinity spikes still persist after performing the cell thermal mass correction due to strong vertical temperature gradients towards the end of the deployment (>0.8 ^∘C m⁻¹); these were manually quality-flagged (marked as bad) on a per-dive basis and not used in the analysis. As the Seaglider CT sensor relies on the forward movement of the vehicle to push water through the measurement cell, salinity and temperature values were also removed from the analysis where the glider speed was considered too slow to provide adequate flow through the cell (<0.1 m s⁻¹).

2.2 NO_x⁻ observations

The LoC ${{\mathrm{NO}}_x}^{-}$ sensor is a miniaturised wet-chemical colorimetric analyser offering high-precision in situ total oxidised nitrogen measurement (Beaton et al., 2012). The analyser takes approximately 20 s to collect a sample and a further 7.5 min to perform an analysis. The concentration of the on-board artificial seawater nutrient standard was determined by laboratory continuous gas segmented flow analysis (Beaton et al., 2012; Vincent et al., 2018). Accuracy was determined with reference to two autoclaved natural seawater standards (KANSO Co. LTD, Japan) and summarised in Table A1. ${{\mathrm{NO}}_x}^{-}$ values below the limit of detection are treated as 0 mmol m⁻³ in our calculations. The limit of detection is 0.025 mmol m⁻³. Vincent et al. (2018) demonstrated excellent stability and agreement between data-generated LoC analysers deployed on an ocean glider and samples simultaneously collected by traditional shipboard techniques and analysed via standard gas segmented flow analysis. This Celtic Sea study region of Vincent et al. (2018) represents a similar shelf sea environment which exhibits a similar range of ${{\mathrm{NO}}_x}^{-}$ concentrations to the North Sea.

As the maximum depth along the AlterEco transect is less than 85 m, the glider was configured to perform periodic loiter dives to enable multiple samples to be collected and so maximise the vertical resolution of the ${{\mathrm{NO}}_x}^{-}$ samples (Vincent et al., 2018). With the standard Seaglider flight mode the glider pitch and buoyancy are adjusted to maintain a uniform speed during both the descent and ascent phase. A loiter dive consists of a standard descent phase, during which a blank and standard measurement are performed, followed by a deliberately slow ascent which provides time for the ${{\mathrm{NO}}_x}^{-}$ sensor to take more samples (n between 3 and 14). Along the transect standard dives took 30 ± 3 min, while loiter dives took 110 ± 5 min. The reduced speed of the glider during loiter ascent results in inadequate flushing of the CT cell, with impaired salinity data quality. We therefore extrapolate the salinity from the descent to the ascent ${{\mathrm{NO}}_x}^{-}$ data assuming the same pressure–salinity relationship.

Photosynthetic fixation of atmospheric carbon into new organic material by phytoplankton is associated with a corresponding consumption of inorganic nutrients (Anderson and Sarmiento, 1994). Where ${{\mathrm{NO}}_x}^{-}$ is the main source of nitrogen, the ${{\mathrm{NO}}_x}^{-}$ uptake represents net community production ($J\left({{\mathrm{NO}}_x}^{-}\right)$) (Plant et al., 2016). As the number and position in the water column of the samples is variable, the inventories are calculated on a per-dive basis after linear interpolation onto a 1 m vertical grid, the sum of which is then calculated to provide a full water column inventory. The standard uncertainty in the inventory is calculated from the squared sum of the relative ${{\mathrm{NO}}_x}^{-}$ concentration uncertainties (± 0.13 %) and a 0.5 m uncertainty in water column height. We exclude three dives where we are unable to calculate a full water column inventory as either the vertical spacing between samples near the thermocline was too large (>10 m) or samples are unavailable from each of the three vertical layers (Fig. 2d). For direct comparison with the Vincent et al. (2018) central Celtic Sea study we also calculate ${{\mathrm{NO}}_x}^{-}$ drawdown within the upper section of the water column down to a 40 m depth. This coincides with the mean euphotic depth for our study region in April 2018 as observed by another AlterEco Seaglider ($z_{\mathrm{EU}}=\mathrm{40}\pm \mathrm{8}$ m) (Hull and Kaiser, 2020). The euphotic zone depth is defined as the depth to which photosynthetically active radiation (PAR) falls to 1 % of the level measured at the surface.

2.3 O₂ observations

During typical (non-loiter) dives the Seaglider's vertical speed was typically 0.1 m s⁻¹, and the oxygen sampling period was 10 s. The optode was multi-point-calibrated (10 concentrations at 4 temperatures each, 40 points total) by Aanderaa in January 2019 and referenced to in situ Winkler titrations taken during deployment and recovery. We adopted the recommendation of Bittig et al. (2018) that in the absence of reference data spanning a wide range of oxygen concentrations a correction factor (gain) should be applied to calibrate the oxygen concentration, i.e. where $C_{\mathrm{corrected}}=a+Cb$, a=0 and b>1. This is in contrast to an offset, a combination of offset and factor, or a correction applied to the raw phase reading of the optode (Uchida et al., 2008). Reference to the discrete samples remained challenging as the deployment location was within a region of significant horizontal heterogeneity. The reference location was as close as possible to the study region but was legally constrained to within 60 nm from the coast due to the small coastal vessel used. Drift in the optode response to oxygen was assumed to be linear with respect to time and corrected based on the deployment and recovery samples (−0.025 % d⁻¹) (Tengberg and Hovdenes, 2014; Bittig et al., 2018).

Biological oxygen production can also be used as a proxy for net community production (J(O₂)), with a net production of oxygen indicative of a net autotrophic system (Alkire et al., 2014; Hull et al., 2016). We implement an oxygen mass balance model to separate the physically and biologically driven components of the observed oxygen inventory changes. As with the ${{\mathrm{NO}}_x}^{-}$ mass balance we calculate values for the euphotic zone (z_I=40 m) and the full water column (z_I=82 m). Full details of this model and its implementation are found in the Appendix Sect. A1.

2.4 Supporting data

Local wind speed and atmospheric pressure were obtained from ECMWF ERA5 reanalysis (Copernicus, 2018). These hourly data were matched to the spatiotemporally closest Seaglider dive. For ocean currents at the AlterEco site we use the UK Met Office Operational Suite Atlantic Margin Model (FOAM) 1.5 km configuration (run 18 May 2020) (Graham et al., 2018). Bathymetry was taken from GEBCO (2019) on a 15 arcsec grid. Between 2007 and 2008 an autonomous monitoring buoy (SmartBuoy; Greenwood, 2016) was located to the south of the Seaglider transect (Fig. 1). This platform collected near-surface ${{\mathrm{NO}}_x}^{-}$ 0.5 and mid-water (25 m) O₂ concentration measurements using similar sensors to the Seaglider and presents an opportunity for determining interannual variability. Full details of this platform are provided in the Appendix A2.

3 Results

3.1 Overview

The Seaglider reached the transect on 29 March 2019 at 04:00 UTC (Fig. 1), at which point the water column was weakly stratified, with the base of the surface mixed layer at 30 m and a 0.03 ^∘C temperature difference across the thermocline (Fig. 2b). Horizontal heterogeneity was seen across the 64 km transect with the eastern end persistently cooler (0.5 to 1 ^∘C) than the western end. Figure 2a, b, and d demonstrate this gradient to be non-linear and differing between the upper and lower layers. This can also be observed in the salinity time series where a more saline region was seen between 1.8 and 2.4^∘ E, with fresher water at both ends of the transect. While the bottom mixed layer (BML) oxygen concentrations on the eastern side of the transect were 3 to 5 mmol m⁻³ higher than elsewhere, the saturation is similar throughout, which indicates that the increased oxygen concentration was mostly due to higher solubility at lower temperatures.

Mean wind speed was 6.7 m s⁻¹, with a peak of 12.6 m s⁻¹. Modelled residual surface currents were to the north-west (counter to the typical annual residual flow), with surface waters moving into the transect area from across the bank. Surface mean current speed was 0.17 m s⁻¹ with a maximum of 0.39 m s⁻¹. We calculate a cumulative total surface advection at the centre of the transect of 52 km over the 27 d occupation of the transect, with a mean residual flow of 0.02 m s⁻¹. Tidal advection was typically 5 km over 25 h. The tidal ellipse is shown in Fig. 1. Bottom water residual flow was small (0.004 m s⁻¹), with a period of flow to the north until 12 April 2019 and then to the south. As a result, bottom water is estimated to have moved less than 9 km while the glider was on the transect. Given the 10 d transit time for the glider to return to the same end of the transect, surface water may have been advected 19 km during that time.

Figure 2 Seaglider location along the transect (a), with associated temperature (b) and mixed layer depth (green line), oxygen concentration (c), salinity (d), and nitrate (e) profiles. Temperature, oxygen, and salinity data are vertically aggregated into 3 m bins for each dive. Dark-grey areas indicate data removed though quality control procedures. Samples which are excluded from the ${{\mathrm{NO}}_x}^{-}$ inventory are signified with a green triangle marker.

There is disagreement regarding the importance of horizontal advection to mass-balance-based estimates in this region (Queste et al., 2013; Rovelli et al., 2016; Große et al., 2016). Given the slow speed of the glider relative to the rates of change in the inventories and uncertainties due to the lack of observations to the north and south, it is not possible to fully resolve any persistent or evolving horizontal ${{\mathrm{NO}}_x}^{-}$ (or O₂) gradient for a 1D mass balance. The tidally driven advection of the glider is largely inconsequential to the data it collects as the glider drifts with the ebb and flow of the tidal currents. Using the modelled residual surface currents and taking the difference in near-surface ${{\mathrm{NO}}_x}^{-}$ concentrations observed across the first transect as indicative of the horizontal ${{\mathrm{NO}}_x}^{-}$ gradient ($\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ mmol m⁻⁴) suggests a horizontal advective ${{\mathrm{NO}}_x}^{-}$ flux of 0.05 mmol m⁻² d⁻¹, which can be neglected for the purposes of the present paper. The effect for O₂ is similarly small.

The thickness of the thermocline varies substantially over the deployment, ranging from 6 to 53 m (see Figs. 2b and A1). Vertical fluxes across the thermocline are only relevant to the mass balance when the base of the thermocline is below the 40 m integration depth (seen in Fig. 2b) and is never relevant for the full water column inventory. Diapycnal eddy diffusivity has been shown to vary over several orders of magnitude in the North Sea (Rovelli et al., 2016). If we assume a typical value for shelf sea diapycnal mixing of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹ and use the largest observed ${{\mathrm{NO}}_x}^{-}$ gradient across the thermocline of 0.77 mmol m⁻⁴, we estimate an upper bound for the diapycnal ${{\mathrm{NO}}_x}^{-}$ flux to be of the order of 0.7 mmol m⁻² d⁻¹. For O₂ we see a negative flux, with a maximal gradient of 4 mmol m⁻⁴ seen briefly in mid-April, resulting in a diapycnal O₂ flux of the order of −4 mmol m⁻² d⁻¹.

3.2 NO_x-based net community production (J(NO_x⁻))

As we focus our analysis on the determination of new production we restrict the period of analysis to be that between peak ${{\mathrm{NO}}_x}^{-}$ concentration (6 April 2019) and the time when the ${{\mathrm{NO}}_x}^{-}$ inventory stops decreasing (18 April 2019). We observe an average surface layer (0–40 m) with an initial inventory of 214 ± 7 mmol m⁻² drawn down to zero over 12 d. A linear least squares fit for this time series, which integrates over three complete transect lengths, finds a nitrate consumption rate as shown in Fig. 3 and reported in Table 1. In each of these analyses the water column depth is considered to be 82 ± 0.5 m based on average depth over the transect length calculated from both GEBCO bathymetry (GEBCO, 2019) and confirmed with the Seaglider altimeter. Converting from integrated to volumetric units for comparison, we find ${{\mathrm{NO}}_x}^{-}$ consumption rates of 0.5 ± 0.04 mmol m⁻³ d⁻¹ for the surface mixed layer (SML) and 0.4 ± 0.03 mmol m⁻³ d⁻¹ across the whole water column.

Figure 3 Water-column-integrated nitrate inventory. Points show inventory values with error bars indicating standard uncertainty. Lines indicate a linear model fit, with the shaded grey area signifying a 95 % confidence interval.

Figure 4 Oxygen mass balance. (a) Air–sea oxygen gas exchange with a negative sign signifying a flux out of the ocean. (b) The change in depth-integrated oxygen inventory. (c) Cumulative net community production. Lines indicate linear model fit to the period associated with nitrate drawdown.

3.3 Oxygen-based net community production (J(O₂))

At the start of the time series the water column was slightly supersaturated with oxygen (103 %–104 % in the SML and 102 % in the BML). The moderate winds associated with this period (8 to 9 m s⁻¹) would be expected to induce no more than 1 % supersaturation (Liang et al., 2013), which suggests a system which is already productive. There was some degree of supersaturation throughout the mission in both the SML and BML, with maximum supersaturation exceeding 120 % between 19 and 20 April 2019. The water column oxygen inventory changes are shown in Fig. 4b. There is net outgassing of oxygen throughout the mission (Fig. 4a), with cumulatively 2.3 ± 0.1 mol m⁻² released to the atmosphere over 27 d. This outgassing is primarily driven by supersaturation through biological production rather than due to warming of the surface waters. Solubility changes account for only 1 % of the observed air–sea oxygen flux during this period. The water column was net productive throughout the mission (Fig. 4c). J(O₂) is seen to be increasing until surface ${{\mathrm{NO}}_x}^{-}$ is depleted (19 April 2019), when a marked drop in production followed. As such J(O₂) follows a pattern consistent with the initiation of the spring bloom followed by ${{\mathrm{NO}}_x}^{-}$ limitation. The water column was likely still slightly net autotrophic at the end of April (Fig. 4c). We estimate cumulative net community oxygen production over the observation period of 3.3 ± 0.1 mol m⁻² (Fig. 4). As with the ${{\mathrm{NO}}_x}^{-}$ mass balance, we compare both the upper 40 m and the full water column with the results tabulated in Table 1.

3.4 O₂ : N ratio

Comparing J(O₂) and $J\left({{\mathrm{NO}}_x}^{-}\right)$ over the bloom observation period we determine the O₂:N ratios, which are summarised in Table 2. We note that the uncertainty associated with these ratios is small compared to previous studies using simultaneous fluxes owing to a more tightly constrained oxygen air–sea gas exchange flux (Alkire et al., 2012). The O₂:N ratio (Table 2) is substantially higher in the surface mixed layer relative to the rest of the water column. This is in keeping with observations from the Alkire et al. (2012) 3 d North Atlantic spring bloom study, which found a SML O₂:N ratio of 14 ± 4.8 and a photic zone (0 to 53 m) ratio of 10 ± 6 for the main bloom period.

Table 1 Summary of ${{\mathrm{NO}}_x}^{-}$ inventory changes. Inventories are calculated using either a depth consistent with the mean euphotic depth (z_I=40 m) or the full water column (z_I=82 m). $J\left({{\mathrm{NO}}_x}^{-}\right)$ is net community production in ${{\mathrm{NO}}_x}^{-}$ equivalents. J(O₂) is net community production in O₂ equivalents.

4 Discussion

4.1 Metabolic state

Our observations demonstrate a spring bloom which is highly productive and short in duration. A short duration makes it challenging to observe with remote sensing as much of the production is missed entirely due to cloud cover of the study region, typical of coastal seas at temperate latitudes. The water column is already net productive by the time the glider is occupying the transect, with weak stratification seen and partitioning of the oxygen and ${{\mathrm{NO}}_x}^{-}$ across the thermocline. While the rates of change in the ${{\mathrm{NO}}_x}^{-}$ concentration, relative to the speed of the glider, make the determination of horizontal gradients difficult, we can determine that in the BML ${{\mathrm{NO}}_x}^{-}$ concentrations are lower towards the eastern end as seen in Fig. 2. Comparing those concentrations observed near 2^∘ E demonstrates a decline in BML concentrations from 4 to 3 mmol m⁻³. The deepening of the thermocline seen between 11 and 16 April 2019 coincides with the glider moving west, and the subsequent shoaling coincides with the glider moving east. This makes it difficult to determine if this deepening is localised to the western end of the transect. However, differences in wind stress and solar heating are negligible across the transect such that we believe the deepening occurred elsewhere along the transect simultaneously. We therefore suggest that while ${{\mathrm{NO}}_x}^{-}$ continues to be consumed, the decline seen in BML ${{\mathrm{NO}}_x}^{-}$ concentrations is primarily through dilution and expansion of the BML to incorporate the nutrient-depleted waters from above.

4.2 Stoichiometry of new production

The Redfield ratio considers the remineralisation of an average marine plankton to have a 138:16 (8.6) O₂:N ratio. These ratios are known to vary systematically across the global oceans (Li and Peng, 2002). This “traditional” ratio appears applicable to the North Atlantic, with remineralisation ratios for the deep North Atlantic yielding 8.6 ± 1 (Li and Peng, 2002).

The ratios we determine for the central North Sea bloom are lower than those observed in the North Atlantic spring bloom (Table 2). ${{\mathrm{NO}}_x}^{-}$ consumption appears to exceed the nitrogen requirements of the average phytoplankton given the observed oxygen production and therefore carbon fixation. However, Redfield ratios are generally not applicable for the whole water column as the simultaneous production of organic matter together with consumption processes modifies the apparent O₂:N ratio. This is demonstrated in Fig. 5, where we can see that our observed O₂:N ratios could be explained by any number of simultaneous remineralisation processes. Hickman et al. (2012) found wide variations in the $\mathrm{C}:{{\mathrm{NO}}_{\mathrm{3}}}^{-}$ uptake ratio associated with changes in the local light climate during summer production in the Celtic Sea. Despite these variations, the particulate organic C:N ratio was observed to be close to the canonical Redfield ratio.

Figure 5 Conceptual diagram which shows how the observed O₂:N ratios (shown with coloured lines) can be explained by the presence of any number of simultaneous remineralisation processes (dashed lines).

Vincent et al. (2018)Alkire et al. (2012)Weston et al. (2005)

Table 2 Voluminal comparison with recent studies.

4.3 Comparison with prior regional estimates

Our mean volumetric drawdown rates are also summarised in Table 2 and are comparable to similar prior studies. A previous deployment of the ${{\mathrm{NO}}_x}^{-}$ sensor integrated with a Seaglider by Vincent et al. (2018) showed a decline in the surface mixed layer (top 40 m) ${{\mathrm{NO}}_x}^{-}$ of −4.3 mmol m⁻³ over 21 d (4 April to 4–25 April 2015; see Table 2) in the central Celtic Sea during the spring bloom. The 2019 central North Sea spring bloom therefore gave 25 %–150 % faster ${{\mathrm{NO}}_x}^{-}$ drawdown rates than the 2015 central Celtic Sea bloom. Direct comparisons of the total inventory change are not possible as the Celtic Sea Seaglider did not sample the stratifying surface waters for ${{\mathrm{NO}}_x}^{-}$ during the onset of the 2015 bloom.

The spring bloom in the seasonally stratified regions of the North Sea is thought to represent 40 % of the typical regional annual production (Richardson and Pedersen, 1998). Several previous studies have derived local production estimates for the region north of the Dogger Bank, but these have focused on the summer months (Weston et al., 2005; Bozec et al., 2005). There are relatively few contemporary NCP observations for this region, and most previous studies aggregate their production estimates over very large spatial or temporal scales. Such a broad approach provides limited information to describe production at the spatial and temporal scales comparable to modern biogeochemical models, which operate on <10 km resolution. Where disparities arise between model and observation, it is difficult to determine which regions or processes are inadequately captured. Therefore, accurate advice from the ocean observing community on how models might be improved is heavily constrained.

Richardson and Pedersen (1998) determined the winter nitrate maximum in the central stratified North Sea to be 8 mmol m⁻³ based on 1993 Convention for the Protection of the Marine Environment of the North-East Atlantic (OSPAR) data, and the 2017 Interim OSPAR Report determined it to be 6 to 8 mmol m⁻³ between 2006 and 2014. The deployment of a similar glider also fitted with an ${{\mathrm{NO}}_x}^{-}$ analyser between November 2017 and January 2018 observed a maximum of 7 ± 0.2 mmol m⁻³ (Birchill et al., 2021). Using 8 mmol m⁻³ as the winter nutrient concentration, taking the typical thermocline depth in a given year to be 30 m with zero vertical mixing, and with a C:N ratio of 5.68, they calculated total annual new production for the stratified North Sea to be 1360 ± 680 mmol m⁻² (C). By contrast our ${{\mathrm{NO}}_x}^{-}$ inventories, with a lower winter nitrate maximum, indicate a total spring bloom new production of 1043 ± 122 mmol m⁻² (C) using the above C:N ratio.

Bozec et al. (2006) calculated dissolved-organic-carbon-based NCP for the International Council for the Exploration of the Sea (ICES) regions, which contain the AlterEco transect (no. 4 and no. 14). For comparison we convert these values of NCP to oxygen equivalents using an O₂:C ratio of 1.4. This conversion provides a surface mixed layer NCP of 0.93 mmol m⁻³ d⁻¹ and a bottom water average NCP of −0.14 mmol m⁻³ d⁻¹ for April (Bozec et al., 2006, do not provide uncertainties). The large temporal and spatial integration range likely fails to resolve much of the local short-term production.

Figure 6 ${{\mathrm{NO}}_x}^{-}$ and oxygen concentration observations from the 2007–2008 North Dogger Cefas SmartBuoy mooring together with the 2019 AlterEco Seaglider (this study). ${{\mathrm{NO}}_x}^{-}$ and temperature observations are from 1 m depth, while oxygen concentrations are from 25 m in 2007 and 35 m in 2018. For 2019 all data between these depths are shown.

4.4 Interannual variability

Winter ${{\mathrm{NO}}_x}^{-}$ concentrations are predicted to reduce in the North Sea due to reduced exchange with the Atlantic Ocean (Gröger et al., 2013). Near-surface ${{\mathrm{NO}}_x}^{-}$ and mid-water oxygen from the 2007–2008 Cefas SmartBuoy mooring (Fig. 1) are shown in Fig. 6 together with the comparable Seaglider observations from 2019. We see that the peak ${{\mathrm{NO}}_x}^{-}$ concentrations prior to bloom onset are consistently 7 ± 0.5 mmol m⁻³ across all observed years within this region. The North Sea biogeochemical climatology (Hinrichs et al., 2017) has a single point for comparison: a discrete near-surface ${{\mathrm{NO}}_x}^{-}$ sample from February 2000 (7.9 mmol m⁻³). We therefore do not see any strong evidence for a reduction in peak winter ${{\mathrm{NO}}_x}^{-}$ with these data.

Extrapolating observations at the surface buoy to an unknown mixed depth to determine an inventory introduces substantial uncertainties. A more complete knowledge of the water column vertical structure is required to accurately determine NCP rates from these observations. We can use a subset of glider observations to demonstrate this; using only the near-surface value to calculate a 40 m inventory results in a 23 % larger calculated ${{\mathrm{NO}}_x}^{-}$ drawdown. However, the trend in Fig. 6 suggests similar maximal rates of ${{\mathrm{NO}}_x}^{-}$ consumption.

The duration and intensity of the ${{\mathrm{NO}}_x}^{-}$ drawdown demonstrate marked interannual variability. While there is an expected inverse correlation between pre-bloom temperature and oxygen concentration, oxygen saturation does not appear to be solely controlled by differing sea temperatures or bubble injection as 2007 is under-saturated in mid-March (93 %), 2008 is very close to saturation (99 %), and 2019 is 2 %–3 % supersaturated. We note that while there is a general trend of increasing sea surface temperatures in the North Sea (Capuzzo et al., 2018), 2019 was cooler than 2007. Große et al. (2016) suggested that year-to-year variations in the oxygen concentrations are mostly caused by variations in local primary production. Warmer years are thought to induce an earlier spring bloom (Capuzzo et al., 2018). This is reflected in these data, with late-March near-surface sea temperatures coinciding with the observed onset of the bloom.

5 Conclusions

We present well-resolved production estimates using two complementary methods from the central North Sea using a single autonomous vehicle. These results represent the most complete examination of the spring bloom formation within this region and demonstrate its highly episodic nature, a feature which is typically hard to observe with typical ship surveys or remote sensing. These observations highlight the strength of autonomous vehicles as a key part of future monitoring strategies, complementing remote sensing and biogeochemical modelling efforts. By comparing the rates of NCP through the two methods we determine the O₂:N ratio of new production within the water column, finding O₂:N ratios in the near surface close to the canonical Redfield North Atlantic ratio, while remineralisation processes in the deeper parts of the water column tend to lower the apparent ratio.

Appendix A

A1 Oxygen NCP model

We can describe the evolution of oxygen over time as the sum of various fluxes:

$$\begin{array}{}\text{(A1)}& z_I{\displaystyle \frac{\mathrm{d}{C}_I}{\mathrm{d}t}}=G+M+E-A+J,\end{array}$$

where C_I is the average oxygen concentration within integration depth z_I, G is the air–sea gas exchange with a flux into the ocean being defined as positive, M is mixing comprising horizontal (M_A) and vertical (M_z) mixing, E is entrainment, A is horizontal advection, and J is net community production.

Air–sea gas exchange is calculated as

$$\begin{array}{}\text{(A2)}& G=k_w\left(C_{\mathrm{sat}}(\mathrm{1}+B){\displaystyle \frac{P}{P_{\mathrm{0}}}}-C_{\mathrm{s}}\right),\end{array}$$

where k_w is the air–sea transfer velocity parameterisation of Nightingale et al. (2000) derived using the airsea R library (Hull and Johnson, 2015), using the Schmidt number parameterisation of Johnson (2010). C_s is the near-surface oxygen concentration. C_sat is the concentration of oxygen if it were in equilibrium with the atmosphere at a pressure of 1013.25 hPa (Garcia and Gordon, 1992; using the fit to the Benson and Krause data). P is local air pressure and P₀ standard atmospheric pressure (1013.25 hPa). B is the wind speed and temperature-dependent bubble supersaturation term from Liang et al. (2013). This bubble term has a relatively large effect on net community production estimates when C_s is close to C_sat (Emerson and Bushinsky, 2016; Hull et al., 2016; Liang et al., 2017). Gas exchange is only relevant within the ventilated partition of the water column, which is equivalent to the actively mixed surface layer (known as mixing layer; Brainerd and Gregg, 1995) (z_mix). If z_I<z_mix, then G needs to be scaled by the factor $z_I/z_{\mathrm{mix}}$.

Horizontal advection can be described by Eq. (A3).

$$\begin{array}{}\text{(A3)}& A=\left(u{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}x}}+v{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}v}}\right)z_I,\end{array}$$

where $\frac{\mathrm{d}C}{\mathrm{d}x}$ is the zonal oxygen gradient, and u is the zonal velocity, with equivalent terms for the meridional (v) direction.

Table A1 Summary of ${{\mathrm{NO}}_x}^{-}$ concentrations in calibration standard.

Figure A1 Example vertical temperature profiles recorded by the Seaglider.

We cannot determine A from our glider data as gliders are semi-Lagrangian; they are displaced by tidal currents, and as such the method of Hull et al. (2020) is not applicable. Gliders are also relativity slow compared with the observed rates of oxygen concentration change and the residual currents along the transect.

We treat the water column as a three-layer system, with an SML and BML separated by a transitional thermocline. Both the SML and BML are considered to be homogeneous, maintained by wind and tidal mixing, respectively. The extent of these layers is typically calculated from the vertical profile of temperature (or density), using the depth at which the temperature is sufficiently different from a reference depth (within the homogeneous layer) as the base of the layer. A reference depth of 10 m, which is commonly used for determining the surface mixed layer depth (z_mix) in open-ocean mass balance studies (de Boyer Montégut, 2004; Kara et al., 2000), is inappropriately deep for our data (Fig. 2b). Data within the top 3 m are of poor quality due to the Seaglider surface manoeuvres and wave action. A reference depth of 5 m is therefore chosen as deep enough to avoid biases from poor surface data quality while shallow enough to determine the air–sea oxygen gradient (see Sect. A1). We adopt a temperature-based threshold for z_mix of 0.2 ^∘C (as per Castro-Morales and Kaiser, 2012), which visually agrees with oxygen profiles (Figs. 2b–d and A1); z_mix is therefore defined here as the shallowest depth at which the temperature differs (warmer) from that at 5 m by more than 0.2 ^∘C. We reverse the procedure for the bottom mixed layer depth and use the deepest temperature recorded as the reference threshold. The choice of SML or BML definition has minimal impact on our analysis as we do not use the SML as our integration depth for the ${{\mathrm{NO}}_x}^{-}$ and O₂ mass balance.

Depending on the choice of integration depth, i.e. the volume of water we calculate NCP over, the fluxes to consider change. For example Barone et al. (2019) eliminated the influence of entrainment by defining their integration depth corresponding to an isopycnal, which is always deeper than the base of the surface mixed layer. Often the integration is restricted to the productive part of the water column (euphotic depth; e.g. Bushinsky and Emerson, 2015, or Binetti et al., 2020). However, the shallow nature of the study area presents an opportunity to reduce the overall uncertainty by calculating the change in oxygen inventory over the whole water column, the depth of which is known with a high level of certainty. There is little variation in depth across the E–W transect, and as such we calculate our full water column mass balance in terms of the average total water column depth (z_I= (82 ± 0.5) m).

As diurnal warming can reduce the ventilated water column to a region of water above the persistent seasonal thermocline, z_mix is not necessarily equivalent to the thermocline depth. It is this volume which is influenced by air–sea gas exchange directly. Typically the dynamic nature of this layer and the difficulty determining its extent from density profiles can introduce significant uncertainties to the magnitude of the air–sea gas flux. The 0.2 ^∘C threshold we use to determine the SML visually agrees well with the oxygen profiles used for determining the value of C_s.

By integrating either the entire water column or to a depth which encompasses the BML, our mass balance differential equation is reduced to Eq. (A4): we simply need to estimate the oxygen concentration within the mixing layer for use in determining the gas exchange flux. This also means entrainment and diapycnal mixing fluxes can be ignored as these do not influence the full water column inventory. Rovelli et al. (2016) observed that within the central North Sea these vertical fluxes can be substantial while being both difficult to quantify and temporally variable. Formulating the mass balance to avoid the need for these terms is thus highly advantageous.

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}}=G+J\end{array}$$

In order to best include our systematic and measurement uncertainties we adopt a Bayesian Markov chain Monte Carlo approach to infer our parameter estimates. We implement the probabilistic model described in Eq. (A5), where the $\widehat{}$ operator signifies our observations of the unknown (latent) state parameter. C_s is the concentration of oxygen in the surface actively mixed layer (z<z_mix). I is the oxygen concentration integrated over our integration depth (z_I); ε_Cs is the uncertainty bias in the oxygen measurements, for this study considered to be ± 1 mmol m⁻³. This model equates to a 15 % error in k_w and a 50 % error in B; ε_I is calculated as the standard error in I; α, β, and σ are left as improper priors, as in a standard least squares regression. Implementation of this model uses the probabilistic programming language Stan (Carpenter et al., 2017). Output variables are presented with 95 % credible intervals, which are the intervals within which an unobserved parameter value falls with a particular probability. These are therefore analogous to frequentist confidence intervals.

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}\mathrm{\Sigma }J& \sim \mathcal{N}\left(\mathit{\alpha }+\mathit{\beta }t,\mathit{\sigma }\right)\\ \mathrm{\Sigma }{J}_t& =\left(\mathrm{\Sigma }{I}_t-I_{\mathrm{0}}\right)-\left(\mathrm{\Sigma }{G}_t\mathrm{\Delta }t\right)\\ G& =k_{\mathrm{w}}\left(C_{\mathrm{sat}}(\mathrm{1}+B){\displaystyle \frac{P}{P_{\mathrm{0}}}}-C_{\mathrm{s}}\right)\\ \widehat{I}& \sim \mathcal{N}(I,{\mathit{\epsilon }}_I)\\ \widehat{k_w}& \sim \mathcal{N}(k_w,\widehat{k_w}\mathrm{0.15})\\ \widehat{C_{\mathrm{s}}}& \sim \mathcal{N}\left(C_{\mathrm{s}},{\mathit{\epsilon }}_{C_{\mathrm{s}}}\right)\\ \widehat{B}& \sim \mathcal{N}(B,\widehat{B}\mathrm{0.5})\\ \widehat{C_{\mathrm{sat}}}& \sim \mathcal{N}\left(C_{\mathrm{sat}},\mathrm{1}\right)\end{array}\end{array}$$

A2 SmartBuoy mooring

The SmartBuoy monitoring station known as “North Dogger” was occupied between 25 February 2007 and 14 September 2008. The moorings consisted of a Cefas ESM2 data logger with an Aanderaa 3835 oxygen optode, Aanderaa 3319B (conductivity and temperature), and an Envirotech NAS-3X nutrient analyser. The NAS-3X performs the same wet-chemistry method as the LoC.