The Three Rs: Resolving Respiration Robotically in Shelf Seas

Ocean deoxygenation threatens ocean productivity, carbon cycling and marine ecosystems. Shelf seas are highly dynamic regions, which contributes to their high productivity and also makes monitoring and constraining their oxygen status a challenge. Here, using the temperate Celtic shelf sea (April and July 2015) as a case study, we present high‐resolution ocean glider observations of turbulence and biogeochemical parameters, demonstrating the potential of these autonomous platforms for environmental monitoring. We estimate vertical turbulent oxygen fluxes be 25% higher in summer than in spring, due to the presence of subsurface chlorophyll and associated oxygen maxima at the seasonal thermocline. We demonstrate that glider‐based estimates were able to constrain similar bottom layer respiration rates as those derived from traditional ship‐based measurements. We suggest ocean gliders are useful monitoring tools that can aid sustainable management of shelf sea ecosystems.

Future shelf sea scenarios indicate prolonged, earlier onset of stratification and stronger thermocline stability (Lowe et al., 2009;Meire et al., 2013;Sharples et al., 2013), all of which affects BML ventilation and O 2 dynamics (Wakelin et al., 2020). However, despite the prominence of gathering O 2 data in order to assess the health of shelf seas, these measurements are either seasonally or spatially limited (Große et al., 2016). It is apparent that there is an immediate need for high spatio-temporal resolution O 2 measurements alongside estimates of physical mixing to elucidate the physical and biogeochemical processes that govern the development of O 2 deficiency, enabling effective ecosystem management of shelf seas.
Here, we present a novel data set of autonomously sampled high-resolution O 2 and turbulent dissipation measurements over a 40-day period in the temperate, seasonally stratified Celtic Sea. The Celtic Sea supports a large phytoplankton bloom in spring, strong stratification in summer and a second bloom event in autumn at the breakdown of stratification (Holligan et al., 1984;Pingree et al., 1976Pingree et al., , 1978. The impact of vertical exchange across the thermocline on nutrient transfer is already locally well characterized (Davis et al., 2014;Sharples et al., 2001;, thus it is a good study site to assess the role of diapycnal mixing on O 2 fluxes and distribution. Using an ocean glider we were able to quantify and assess the diapycnal exchange of O 2 between the SML and BML. Finally, we use these measurements, combined with estimates of the horizontal dispersion and advection, to calculate respiration during spring and summer 2015.

Glider Measurements
Ocean gliders were deployed to conduct "virtual mooring" profiles at a study site in the seasonally stratified central Celtic Sea (station CCS, 49° 24' N, 8° 36' W; see Figure 1) during spring 2015 (6th-28th April, decimal day 95-117) and summer 2015 (15th July-2nd August, decimal day 195-213). The integrated approach adopted in this study, combining ship based and glider measurements enabled estimates of spatial gradients while also minimising tidal aliasing that would likely be introduced by long spatial transects with the glider. A Slocum (Teledyne Webb Research, Falmouth, USA) Ocean Microstructure Glider (OMG, see Palmer et al., 2015) was equipped with a MicroRider microstructure package (Rockland Scientific International) to measure the microstructure of velocity shear, a Seabird SBE42 CTD sensor to measure temperature, salinity and pressure, and an Aanderaa 4831 oxygen optode to measure O 2 (precision 0.2 μmol kg −1 ). Measurements were taken within 5 m of the bed and 2 m of the surface on most dives, with each yo-yo profile taking approximately 20 min. The glider AA4831 optode is known to experience severe lag across strong oxygen gradients, and therefore oxygen data was corrected where possible for optode membrane lag following Bittig et al. (2014). In comparison to other oxygen optodes, the AA4831 has been documented by various scientific studies as being an extremely stable optode with low detectable drift (<0.5% yr −1 ) and high precision of <0.2 μmol kg -1 (Champenois & Borges, 2012;Johnson, 2010;Kortzinger et al., 2004;Nicholson et al., 2008). Optode drift was calculated in this study by comparing discrete Winkler-analyzed samples taken at deployment and recovery of the gliders, identifying a downward drift of 0.001% d −1 , in close agreement with quoted manufacturer values.
Glider sensors (temperature, salinity and 2 ) were calibrated against nearby ship CTD profiles (CTD calibrated 1 month prior to cruise, SBE 43 precision = 2% of 2 saturation) and discrete water samples collected within 3 hr and 2 km of glider deployment and recovery times and glider position, respectively, as part of the Shelf Sea Biogeochemistry program (RRS Discovery, DY029 and DY033). Error estimates for the total change in 2 BML (μmol kg −1 ) were calculated as the sum of the optode precision (0.2 μmol kg −1 ) and drift over the entire respective deployments (<0.1 μmol kg −1 ). Currents were monitored throughout the glider deployments by a mooring at the CCS study site, which was equipped with an acoustic current profiler (ADCP), salinometer and thermistors that provided near-continuous data (

Oxygen Budget Calculations
The change in 2 BML can be represented by the following equation: where C( 2) represents the concentration of 2 (μmol kg −1 ). The term on the LHS refers to the change observed in 2 BML over the sampling periods. The first three terms on the RHS represents the advection of water with respect to the x, y (horizontal), and z (vertical) coordinates. The fourth and fifth terms represent the vertical and horizontal turbulent flux of oxygen into the BML respectively, where K z refers to the vertical eddy diffusivity, K h refers to the horizontal eddy diffusivity and BML represents the thickness of the BML (m). The final term, Δ ( 2) , represents the net effect of biological processes on oxygen. In the BML, this will be the depletion of 2 due to nitrification and bacterial respiration.

Vertical Diffusivity
Vertical turbulent diffusion (K z , m 2 s −1 ), was calculated using profiles of the rate of dissipation of turbulent kinetic energy (TKE; ε, m 2 s −3 ) and the buoyancy frequency (N 2 , s −2 ) following Equation 2 (e.g., Osborn, 1980): Here Γ, the mixing efficiency, is considered constant at 0.2 for the stratified water column (Osborn, 1980;Tweddle et al., 2013;Williams, Sharples, Green, et al., 2013). While there is ongoing discussion on the assumption of a constant mixing efficiency in stratified fluids, the vast majority of studies have not conclusively arrived at a suitable improvement to that proposed by Osborn (1980), so this simple solution has been employed here as current best practice (Gregg et al., 2018).
The buoyancy frequency was calculated using the density profiles (Equation 3): where is the acceleration due to gravity (9.81 m s −2 ), ρ is density, and z is the vertical coordinate (m, positive upward).
The base of the pycnocline represents the interface for diapycnal mixing and 2 fluxes and was defined as the depth that marked the top of the BML. We chose this to be where the density decreased by 15% of the total water column density change. This allowed us to standardize our method over various densities during spring and summer.

Vertical O 2 Flux
Vertical turbulent driven fluxes of 2 from the base of the pycnocline into the BML may be calculated following Equation 4 (e.g., Sharples et al., 2001) where ( 2 ) represents the vertical gradient in oxygen (mmol m −4 ) and K z is the eddy diffusivity defined in Equation 2. Confidence limits (95%) for 2 were calculated using Efron Gong Bootstrap resampling method (Efron & Gong, 1983).

Horizontal Advection of O 2
The Celtic Sea is bordered by the Northeast Atlantic to the west, St Georges Channel and the Irish Sea to the north, and the Bristol and English Channels to the northeast and east, respectively. Fronts form at each of these boundaries which have the potential to generate eddies that transfer properties (e.g., salt, heat and 2 ) and momentum laterally across the shelf (Huthnance et al., 2009;Proctor et al., 2003) before being eroded by wind and tidal mixing (Badin et al., 2010). In the Celtic Sea, internal tides can be large and non-linear resulting in solitons that 5 of 11 can drive on-shelf exchange of oceanic water (Huthnance et al., 2009). In April 2015 and August 2015, oceanic waters moved on-shelf in the bottom layer at 1.4 and 1.5 km d −1 respectively, due to wind-driven Ekman transport, cross-shelf pressure gradients and/or internal tidal wave Stokes drift (Ruiz-Castillo et al., 2019). Based on these estimates, shelf edge waters would take 110 and 150 days to travel from the shelf edge to our study site respectively. It is important to note that there are likely large uncertainties in these estimations, however the estimations by Ruiz-Castillo et al. (2019) are still the best measurements of horizontal fluxes in this area. During this period, the bottom layer horizontal oxygen gradient indicated that 2 concentrations were lowest at the shelf edge and increased with distance on shelf to the highest concentrations at our study site ( Figure S1 in Supporting Information S1). This negative on-shelf gradient suggests that water with comparatively lower 2 concentration would be slowly advected to our study site.

Horizontal Dispersion of O 2
Observations of horizontal dispersion coefficients (K h ) in shelf seas are typically between 10 and 600 m 2 s −1 (Houghton et al., 2009;Sanders & Garvine, 2001). In order to calculate the horizontal dispersion of 2 the derivative of the lateral gradient in 2 is needed (Equation 1).

Water Column Characteristics
The onset of thermal stratification at our site in the central Celtic Sea occurred on the 6th April (decimal day 95) coinciding with a period of reduced wind mixing and increased solar heating (Wihsgott et al., 2019). This supported the development of an intense spring bloom occurring between days 95 and 114 (Poulton et al., 2019), where 2 was relatively high throughout the water column (280-310 μmol kg −1 ) as a result of increased primary productivity. There was less than 1°C difference between SML and BML temperature, separated by a 30 m thick thermocline (Figure 1a). Tidal currents were 0.2-0.4 ms −1 with the fastest currents in the SML ( Figures S3 and  S4 in Supporting Information S1). The BML temperature at CCS was 2°C warmer than that at the shelf edge ( Figure S1 in Supporting Information S1). O 2 BML at CCS was 17.7 μmol kg −1 larger than O 2 BML at the shelf edge ( Figure S1 in Supporting Information S1) during spring 2015 therefore on-shelf advection in the BML (e.g., Ruiz-Castillo et al., 2019) would have diluted O 2 BML concentrations at our study site.
Stratification was more pronounced in July with the SML being 5°C warmer than the BML (Figure 1c). The vertical gradient in O 2 was also different from spring, with the lowest concentrations in the SML (249 μmol kg −1 ), then the BML (253 μmol kg −1 ), and the highest concentrations within the 15 m thick seasonal thermocline (273 μmol kg −1 , Figure 1c). Following the depletion of nutrients in the SML, phytoplankton thrive at the base of the pycnocline, where nutrients (mixed up from the BML) and light are available (Sharples et al., 2001). Phytoplankton photosynthesize at the base of the pycnocline, and thus a SCM as well as O 2 maxima, is observed here (Figure 1d). The thermocline and depth of the base of the pycnocline were stronger and deeper, respectively, in summer than in spring, with Z BML decreasing from 95 m in spring to 85 m in summer (Figures 1b and 1c). The vertical oxygen gradient (dO 2 /dz) at the base of the pycnocline was more than twice as strong in summer than in spring (Table 1). O 2 BML at CCS was 7 μmol kg −1 larger than O 2 BML at the shelf edge (Table 1), indicating a decrease in the across shelf horizontal O 2 BML gradient in summer.
The strongest mixing and highest values of K z at the base of the pycnocline were observed on day 107 (6.2 × 10 −3 m 2 s −1 , Figure 2d), corresponding with near-gale conditions on day 107 (wind speeds >17 m s −1 , see Figure S2 in Supporting Information S1) occurring immediately prior to spring tides (day 108-110). Persistent strong winds on days 107-113 ( Figure S2 in Supporting Information S1) prolonged strong mixing and a relatively WILLIAMS ET AL.
10.1029/2021GL096921 6 of 11 high oxygen flux beyond the spring-tide period (Figure 2). During summer, high levels of diapycnal K z exceeding 1.8 × 10 −3 m 2 s −1 were also observed during spring tides (days 196-198, 210-212; Figure 2d). The lowest values of K z at the base of the pycnocline were at molecular levels (O (1e-9) m 2 s −1 ) observed immediately after neap tides, and highest values were observed to coincide with spring tides (Figure 2d).

Vertical Mixing of O 2
The diapycnal 2 flux ( 2 ) ranged from O (1e-7) to 6.49 × 10 −4 mmol m −2 s −1 , with the highest 2 (6.49 × 10 −4 mmol m −2 s −1 ) being observed during summer (day 211) coinciding with K z spikes at spring tides. The lowest 2 fluxes were observed immediately after neap tides concurring with low K z values. Average 2 fluxes across the base of the seasonal pycnocline were higher in July by ∼25%, this was due to the vertical 2 gradient being twice as strong during July (Table 1). It is important to note that the calculation of K z (Equation 2) considers a constant mixing efficiency and thus does not take into account the dampening effect of stronger stratification on diapycnal mixing in summer compared to spring, which has been shown in previous studies to potentially affect the vertical flux (Schultze et al., 2017). Fluxes were almost a factor of 10 larger during spring tides than neap tides in both spring and summer driven by enhanced K z .

Horizontal Advection and Diffusion of O 2
The horizontal advection of 2 to our study site was calculated using the values of horizontal advection (Ruiz-Castillo et al., 2019) together with the measured horizontal gradient (d 2∕ ) in the BML for both spring and summer ( Figure S1 in Supporting Information S1). The horizontal 2 gradient from the shelf edge to CCS in the BML was −15.5 × 10 −5 mmol m −4 during spring, the strength of this gradient reduced in summer to −6.1 × 10 −5 mmol m −4 ( Figure S1 in Supporting Information S1). The horizontal advection of on-shelf waters acted to dilute O 2 at our study site during both spring (−0.2 mmol m −2 d −1 ) and summer (−0.1 mmol m −2 d −1 ).
Using the upper and lower estimates of horizontal dispersion across the shelf of 10-600 (m 2 s −1 ; Houghton et al., 2009;Sanders & Garvine, 2001), together with the derivative of the lateral gradient (d 2 2 /dx 2 ) in the BML from J4 to the shelf edge (Table 1), we were able to quantify the horizontal dispersion of 2 to our study site. The horizontal dispersion acted as a source of 0.1-6.6 mmol m −2 d −1 of 2 to CCS during spring. During July however, the horizontal dispersion acted to dilute the 2 at our study site by −0.1 to −1.5 mmol m −2 d −1 .

Estimates of Respiration
The net effect of biological processes (Δ ) on 2 BML during the sampling period can be estimated from the supply rate of 2 to the BML via the vertical and horizontal fluxes and the in situ 2 BML concentrations measured by the gliders (rearranging Equation 1). Once the BML was isolated from direct air-sea gas exchange by stratification, there was no change in 2 concentration due to reduced gas solubility associated with the increase in BML temperature of 0.003°C d −1 (Figure 2a). During the two separate glider sampling periods, the depth-integrated 2 BML was observed to decrease at a rate of −28.5 ± 0.10 mmol m −2 d −1 in spring (decimal days 95-114), and −8.5 ± 0.09 mmol m −2 d −1 in summer (decimal days 193-212), with 2 saturation in the BML reaching 92% by decimal day 234 (Figure 2). Rearranging Equation 1, the depth-integrated Δ was estimated to be in the range 23.9-32.9 mmol m −2 d −1 in spring and 8.6-15.4 mmol m −2 d −1 in summer, indicating net consumption of 2 in the BML during both periods.

Discussion and Conclusions
Oxygen decline is evident in several shelf seas, with large regions of the Northwest European continental shelf seas having the potential to become seasonally deficient in 2 in late summer under current conditions (Ciavatta et al., 2016) and more extensively under future predicted climate scenarios (Wakelin et al., 2020). In this study, using a novel high-resolution data set collected using ocean gliders, we accurately balanced the 2 inventory in the central region of a broad, stratified shelf sea.
Our high resolution glider-based 2 measurements in the Celtic Sea were in good agreement with nearby shipbased CTD measurements. For example, mean BML 2 concentrations from the glider were 289 ± 2 μmol kg −1 in spring and 261 ± 2 μmol kg −1 in summer, respectively, with equivalent values from the CTD of 289 ± 1 and 261 ± 1 μmol kg −1 , respectively. Error estimates for the total change in 2 BML (μmol kg −1 ) were calculated as the sum of the optode precision (0.2 μmol kg −1 ) and drift (<0.1 μmol kg −1 ) over the entire respective deployments (±0.3 μmol kg −1 ). Good agreement between our K z values at the base of the seasonal pycnocline and those previously reported for the Celtic Sea  and other temperate shelf seas (e.g., Mackinnon & Gregg, 2003) suggests that we obtained accurate K z measurements at high spatial and temporal resolution during our study. A caveat of our method is that the calculation of K z (Equation 2) considers a constant mixing efficiency and thus does not take into account the dampening effect of stronger stratification on diapycnal mixing in summer compared to spring, which has been shown in previous studies to potentially affect the vertical flux (Schultze et al., 2017). There is no consensus for a constant mixing efficiency there is agreement that below a certain buoyancy Reynolds number (Re b ) turbulence is not mixing diapycnally (Schultze et al., 2017), thus Equation 2 would be valid only in developed turbulence. If we compute the minimum diffusivity (K z min ) of an isotropic turbulence patch using Re b > 7 to (Schultze et al., 2017), this gives a K z min value of ∼1.4 × 10 −6 m 2 s −1 (Figure 2d, dashed line). Only a small number of diffusivities fall below this (Figure 2d), indicating that a constant mixing efficiency of 0.2 is suffice for the majority of our study, but in low energy environments our method of calculating the vertical 2 flux might not hold.
Diapycnal mixing in temperate shelf seas maintains the ecosystem in a variety of ways; (a) supplying nutrients to the euphotic zone (Sharples et al., 2001, 2007, Williams, Sharples, Green, et al., 2013, (b) supplying organic matter to the BML (Davis et al., 2014; (c) supplying 2 to the BML (Rovelli et al., 2016). Downward 2 fluxes were largest during spring tides, with fluxes increasing 8-fold in spring and 13-fold in summer (Table 1) relative to neap periods. Despite similar K z values, 2 fluxes were more 25% larger during summer than spring due to a stronger ( 2 ) gradient associated with the oxygen maximum (Figures 1b and 1c; Table 1). During high wind conditions (>14 ms −1 ) spring tide 2 fluxes were still a factor of 10 larger than during neap tide under similar wind conditions (e.g., 50 × 10 −5 mmol m −2 s −1 and 0.4 × 10 −5 mmol m −2 s −1 , respectively, on days 107 and 204, Figure 2d, Figure S1 in Supporting Information S1). These 2 fluxes were significantly smaller than values measured over a 3 days study in the North Sea by Rovelli et al. (2016) (∼54 mmol m −2 d −1 ), which were comparable to the biological consumption of 2 in the BML (∼60 mmol m −2 d −1 ). Conversely, in our study, it appears that 2 provides a steady leakage of 2 into the BML which is smaller than Δ indicating that 2 is not enough to sustain the biological 2 demand but does act to slow the rate of 2 BML depletion in the BML. A recent study concluded that the vertical downward mixing of 2 was not sufficient to maintain the hypoxic transition zone of the Baltic Sea (Holtermann et al., 2020). In contrast to our study, horizontal intrusions were an important 2 source in the Baltic.
Water of lower 2 concentration was advected from the shelf edge to our study site and was twice as strong in spring than summer (0.2 and 0.1 mmol m −2 d −1 respectively). This seasonal difference was due to the horizontal gradient 9 of 11 across the shelf in 2 being stronger in spring than summer (−15.5 × 10 −5 mmol m −4 and −6.1 × 10 −5 mmol m −4 respectively). The greater lateral 2 gradient during spring was partly due to a supersaturated BML at CCS following the spring bloom ( Figure S1b in Supporting Information S1). Our measurements indicated that 2 consumption by biological processes occurs in the BML in both spring and summer, with the estimated Δ being 3-fold higher in spring compared to summer (68 ± 0.1 mmol m −2 d −1 and 20 ± 0.1 mmol m −2 d −1 , respectively, Table 1). Elevated 2 consumption in the BML during spring likely reflect the higher availability of organic matter for remineralization via downward sinking particles and fresher dissolved organic matter in spring relative to summer Garcia-Martin et al., 2018). Our estimates of net 2 consumption are in good agreement with community respiration rates reported by Garcia-Martin et al. (2018) at the same study site during spring and summer (36-133 mmol m −2 d −1 and 9-26 mmol m −2 d −1 , respectively; Table 1). Our findings suggest that BML respiration rates can be accurately estimated using this method of combined glider and ship-based observations, providing insight into the importance of various consumption processes on BML 2 . The results we have presented of high resolution shelf sea 2 distribution and vertical diffusivity have only been possible by using ocean gliders, capable of measuring over long periods and during varying sea states that would be problematic from a research vessel. Our results show that diapycnal mixing supplies 2 to the BML over the pycnocline, varying relative to tides and wind events and modulated by the strength of vertical O 2 stratification. Modifications to these controls has the potential to have important consequences for marine ecosystem health in shelf seas. Future IPCC climate scenarios predict an increase in water temperature, stronger and more prolonged stratification (Lowe et al., 2009), and reduced 2 concentrations in the BML in shelf seas (Schmidtko et al., 2017;Wakelin et al., 2020). Diapycnal 2 fluxes observed in this study may become a critical mechanism aiding the ventilation of the BML, helping mitigate 2 deficiency on continental shelf seas and coastal regions.