Technical note: Novel triple O₂ sensor aquatic eddy covariance instrument with improved time shift correction reveals central role of microphytobenthos for carbon cycling in coral reef sands

The aquatic eddy covariance technique stands out as a powerful method for benthic O₂ flux measurements in shelf environments because it integrates effects of naturally varying drivers of the flux such as current flow and light. In conventional eddy covariance instruments, the time shift caused by spatial separation of the measuring locations of flow and O₂ concentration can produce substantial flux errors that are difficult to correct. We here introduce a triple O₂ sensor eddy covariance instrument (3OEC) that by instrument design eliminates these errors. This is achieved by positioning three O₂ sensors around the flow measuring volume, which allows the O₂ concentration to be calculated at the point of the current flow measurements. The new instrument was tested in an energetic coastal environment with highly permeable coral reef sands colonised by microphytobenthos. Parallel deployments of the 3OEC and a conventional eddy covariance system (2OEC) demonstrate that the new instrument produces more consistent fluxes with lower error margin. 3OEC fluxes in general were lower than 2OEC fluxes, and the nighttime fluxes recorded by the two instruments were statistically different. We attribute this to the elimination of uncertainties associated with the time shift correction. The deployments at ∼ 10 m water depth revealed high day- and nighttime O₂ fluxes despite the relatively low organic content of the coarse sediment and overlying water. High light utilisation efficiency of the microphytobenthos and bottom currents increasing pore water exchange facilitated the high benthic production and coupled respiration. 3OEC measurements after sunset documented a gradual transfer of negative flux signals from the small turbulence generated at the sediment–water interface to the larger wave-dominated eddies of the overlying water column that still carried a positive flux signal, suggesting concurrent fluxes in opposite directions depending on eddy size and a memory effect of large eddies. The results demonstrate that the 3OEC can improve the precision of benthic flux measurements, including measurements in environments considered challenging for the eddy covariance technique, and thereby produce novel insights into the mechanisms that control flux. We consider the fluxes produced by this instrument for the permeable reef sands the most realistic achievable with present-day technology.

1 Introduction

This study introduces a new eddy covariance instrument and demonstrates its functionality through a measuring series addressing benthic oxygen flux in a dynamic backreef area covered by highly permeable carbonate sands. The aquatic eddy covariance technique is a powerful technique for quantifying fluxes at the seafloor as it measures over any type of substrate and integrates over a relatively large area (Berg et al., 2003; Lorrai et al., 2010; McGinnis et al., 2008). This non-invasive technique derives flux by averaging the turbulent vertical advective transport of O₂ (or other solutes) above the sediment over time (Berg et al., 2003; Lorrai et al., 2010; McGinnis et al., 2008). Since the vast majority of the recorded flux signals originates from the seafloor upstream of the instrument (termed “footprint”), the flux measurements include the effects of the natural flow and light fields, as well as the benthic sedimentary community (Berg et al., 2003; Lorrai et al., 2010; McGinnis et al., 2008; Huettel et al., 2020). The eddy covariance method thus can produce flux data in dynamic shelf ecosystems such as coral reefs that are strongly influenced by flow and light (Long et al., 2013; Long, 2021; Huettel et al., 2020). The technology so far has been adapted to measure temperature, salinity, oxygen, hydrogen, sulfide, and nitrate fluxes (McGinnis et al., 2011; Johnson et al., 2011; Long et al., 2015; Crusius et al., 2008; Weck and Lorke, 2017).

Since the fluxes are calculated from minute concentration changes measured at a high frequency required to account for all water movements that transport the solute, it is critical that the flow data are accurately aligned in time with the associated solute data. This presently is a potential source of error as conventional aquatic eddy covariance instruments cannot measure current velocities and solute concentrations in the same location. Typically, current velocities are recorded with an acoustic Doppler velocimeter (ADV) and solutes with a fast-responding electrochemical or optical sensor (Kuwae et al., 2006; Berg et al., 2003; Reimers et al., 2012; Attard et al., 2016; Donis et al., 2016; Glud et al., 2010; McGinnis et al., 2014; Lorke et al., 2013; Huettel et al., 2020). The tip of the solute sensor in these instruments is positioned at a few centimetres horizontal distance from the ADV's measuring volume to prevent disturbances of the flow and the acoustic ADV signal. This spatial separation of flow and solute measurements causes a misalignment between the two measurement time series, which requires a time shift correction of the data. In environments with dynamic currents, this misalignment changes continuously as direction and velocity of the turbulent flow varies. Algorithms were developed that shift the O₂ data in time such that they are synchronised with the velocity data (McGinnis et al., 2008; Berg et al., 2015; Reimers et al., 2016). A common procedure is to move a short sequence (e.g. 15 min) of solute data in time relative to the current flow data recorded at that time until a maximum in flux is reached. In steady unidirectional flow, this procedure largely can eliminate time shift errors, but it is difficult to apply an effective correction in dynamic settings (Donis et al., 2015; Reimers et al., 2016). Since the rapid changes in solute concentration and vertical flow velocity are relatively small and affected by signal noise, a distinct maximum in flux may not be found when time shifting the data, which can result in erroneous corrections and fluxes. Furthermore, wave orbital motion in shelf environments produces oscillating bottom currents that may change in magnitude and direction at the timescale of seconds, complicating a correct alignment of the data and producing further potential sources of uncertainty in the flux calculations. In the conventional eddy covariance instruments with one or two solute sensors, the cumulative effect of small errors in the time shift correction thus can lead to significant under- or overestimates of the flux, which in extreme cases can reverse the direction of the calculated flux relative to the true flux (Berg et al., 2015; Reimers et al., 2016). To remove this potential source of error, we designed a triple O₂ sensor eddy covariance instrument (3OEC) that eliminates the uncertainties caused by the spatial separation of flow and concentration measurements.

The instrument was tested in the Florida Keys at an exposed inner shelf site with carbonate sands, clear oligotrophic water, and substantial wave action, i.e. in an environment considered challenging for eddy covariance measurements due to the low particle concentrations in the water (ADV measurements rely on sound reflection from particles) and the dynamic flows (causing the data misalignments addressed in this study). We selected that site because carbonate sand beds are an integral part of reef environments and may play a central role in their carbon and nutrient cycles (Eyre et al., 2018; Santos et al., 2011; Cyronak et al., 2013). As warm water coral reefs grow in shallow, high-energy environments, these sediments typically are dominated by highly permeable coarse sands (Yahel et al., 2002; Harris et al., 2015) colonised by microphytobenthos (Werner et al., 2008; Jantzen et al., 2013). Owing to the rapid pore water exchange facilitated by the high permeability, biogeochemical processes in the sediment surface layers can respond almost instantly to changes in flow, organic matter input, and light (Huettel et al., 2014). This makes the eddy covariance technique the unrivalled method for measuring interfacial fluxes in this environment, provided that potential errors associated with time shift correction in dynamic flows can be removed.

The goals of this study therefore were to develop an eddy covariance instrument not affected by uncertainties caused by spatial separation of flow and solute measuring points and to demonstrate its functionality through measurements in a dynamic coral reef environment. The instrument measured O₂ flux, which is considered a good proxy for benthic production and respiration (Berg et al., 2013; Glud, 2008). Comparison with parallel measurements with a conventional eddy covariance instrument reveals the improvements achieved by the new instrument design.

Figure 1 The triple O₂ sensor eddy covariance instrument (3OEC). (a) The 3OEC deployed over carbonate sand at the study site in the Florida Keys. (b) Positioning of the ADV sensor head (1) and the three O₂ sensors (A, B, C). The average O₂ concentration calculated from the readings of the three O₂ sensors approximates the concentration at the geometric centre of the red triangle (side lengths 6.4 cm) defined by the three O₂ sensor tips. This geometric centre is also the centre of the flow measuring volume of the ADV, marked as yellow cylinder. (c) Vertical view of the positioning of the O₂ sensors around the ADV measuring volume (yellow circle) that is located 15 cm below the central sensor stem. The sensor tips are located at 3.0 cm horizontal distance from the edge of the ADV measuring volume. The red triangle is the vertical view of the red triangle shown in (b). (d) Proof of the equivalence of the average of the three sensor readings and the concentration at the centre of the equilateral triangle defined by the positions of the three optode sensor tips (red triangle in b and c). We postulate that the concentration gradients between the sensor tips are linear (see also text). Accordingly, the half-way point of a side of the triangle corresponds to the average concentration measured by the optodes at the two endpoints of that side. Applying the law of sines and triangle congruence criteria (transitive property of congruence, angle bisector theorem, converse of angle bisector theorem), the concentration at the centre point (M) of the equilateral triangle (ABC) equals the concentration at the vertex (M) of the right triangle defined by one vertex of the equilateral triangle (C in the above example) and the midpoint ((A + C) $/$ 2 in the above example) between that vertex and the second vertex (A in the above example) on that line. According to the congruence criteria, this is valid for analogous, congruent right triangles constructed on the other sides of the triangle. As these triangles are based on the average concentration of two vertices of the equilateral triangle, it follows that the centre point concentration is equivalent to the average concentrations measured by the three sensors. The heat map visualises the equivalence of the concentrations at the centre of the equilateral triangle calculated using this approach and the average of the three sensor signals. Photographs: Markus Huettel.

2 Methods

2.1 Triple O₂ sensor eddy covariance instrument (3OEC)

The new 3OEC instrument utilises a data averaging approach to remove misalignments in time caused by the spatial separation of flow and O₂ measurements and thereby eliminates potential errors caused by time shift corrections. The 3OEC measures simultaneously with three O₂ fibre optodes positioned at 120^∘ angular spacing in the same horizontal plane around the centre point of the water volume where current flow is measured by the ADV (Fig. 1). Assuming approximately linear concentration gradients within the 6.4 cm distance between the O₂ sensors – justifiable according to planar optode readings of water column oxygen distributions in turbulent flows (Glud et al., 2001; Larsen et al., 2011; Oguri et al., 2007) – the O₂ concentration at the location where the flow is measured can be calculated through averaging of the three simultaneous sensor signals. The three optode tips, positioned at the corners of an equilateral triangle, present three equidistant points on a circle with the measuring volume of the ADV at its centre. For linear concentration gradients between the three measuring points, it can be proven with the law of sines that the O₂ concentration at the centre of this circle corresponds to the average of the three sensor signals (Fig. 1d).

The optodes are ultra-high-speed O₂-needle sensors (PyroScience™ – OXR430-UHS; Table A1) with a response time of 200–300 ms (Merikhi et al., 2018). The three optodes are pointing downward at a 45^∘ angle, with their sensing tips positioned at 3.7 cm horizontal distance from the centre of the flow measuring volume. This placement, within the recommended distance of 10 Kolmogorov scale lengths from the ADV measuring volume (Lorrai et al., 2010), prevents any disturbance of the flow within that volume and potential interferences with the acoustic pulses of the ADV. The sensors are read by three FireStingO2-Mini O₂ meters (PyroScience™; Table A2). The ADV is a Nortek™ Vector acoustic Doppler velocimeter (Table A3) that measures the 3D velocity field within a cylindrical measuring volume (1.5 cm diam. × 1.5 cm, located 15 cm below the central acoustic transducer) at a sampling rate of 32 Hz. A DataQ DI-710-UH USB data logger (14-bit A/D conversion) records simultaneously the output of O₂ meters and the ADV at a rate of 64 Hz to prevent aliasing. O₂ meters and data logger are contained in an underwater housing (A.G.O Environmental Electronics) fitted with three PyroScience™ fibre feedthrough plugs for connecting the O₂ sensors, as well as Impulse micro inline plugs for connecting the ADV and external battery (4 × lithium-ion 12 V, 50 Wh). The ADV, O₂-meter housing, and battery pack are mounted on a stainless-steel tripod with 1.2 m side length and 1.2 m height (Berg and Huettel, 2008). In addition, the frame carries a miniDOT O₂ logger (PME) and an Odyssey PAR (photosynthetically active radiation) logger (Dataflow Systems) for collection of temperature and O₂ reference data (once per minute) and PAR data (once per 10 min), respectively.

2.2 Dual O₂ sensor eddy instrument (2OEC)

An eddy covariance instrument with conventional sensor configuration (2OEC) was deployed parallel to the 3OEC to analyse the potential flux error caused by the time shift and the effectiveness of standard data corrections. The 2OEC, described in detail in Huettel et al. (2020), measures simultaneously with two O₂ optodes positioned on one side of the ADV measuring volume with their measuring tips 1 cm horizontally apart. Deployments of this instrument at the same study site in the Florida Keys simultaneously with benthic advection chambers (Huettel and Gust, 1992; Janssen et al., 2005; Huettel et al., 2020) produced eddy covariance fluxes (3.7 ± 0.9 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) that were similar to those of the chamber fluxes (3.9 ± 3.0 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) during daytime and of similar order of magnitude during nighttime (2OEC: −2.5 ± 1.3 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; chambers: −3.4 ± 0.8 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). These fluxes obtained with an independent measuring technique corroborate the magnitude of the eddy covariance fluxes, but it should be noted that the chambers do not account for changes in flow and organic matter supply during the incubation, which both can have a significant influence on the flux. While the chambers under relatively steady conditions and short incubation periods can produce fluxes similar to those recorded by eddy covariance instruments, discrepancies between fluxes measured by the two techniques were observed in dynamic environments (Berg et al., 2013).

2.3 Data processing

Eddy covariance flux calculations are based on the assumption that the flux signal is transported by a bottom current with steady state mean flow and O₂ concentration that reaches the instrument unobstructed after passing the footprint area (Massman and Lee, 2002; Baldocchi, 2003; Kuwae et al., 2006; Berg et al., 2007). In coastal environments, such conditions rarely are met, requiring post processing of the flux data to correct for infringements of these assumptions, as well as errors caused by technical limitations (Holtappels et al., 2013; Reimers et al., 2016; Huettel et al., 2020). We applied the same routine corrections to 3OEC and 2OEC data for compensation of design and sensor limitations, as well as for non-steady-state O₂ concentrations in the water column. The unfiltered flow and O₂ data records were reduced from 64 to 8 Hz by averaging, which reduced noise but maintained sufficient resolution to describe the entire frequency spectrum carrying the flux signal. For each 8 Hz time point, the average signal of the three O₂ sensors of the 3OEC was calculated to determine an estimate of the O₂ concentration in the centre of the ADV measuring volume. Similarly, the signals of the two sensors of the 2OEC were averaged to produce mean concentrations. O₂ fluxes then were calculated based on these averages, as well as based on the signals of each individual sensor using the software EddyFlux 3.2 (Peter Berg, unpublished). The software determines mean O₂ base concentrations for 15 min time segments through Reynolds decomposition (Lorrai et al., 2010; Berg et al., 2009; Lee et al., 2004). Within each 15 min interval, the mean O₂ concentration $\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{2}}}$ (defined as a least-square linear fit to the data) then is subtracted from each 8 Hz O₂ data point to arrive at the instantaneous O₂ fluctuation ${\mathrm{O}}_{\mathrm{2}}^{\prime }$ for that time point. The instantaneous vertical velocity $V_z^{\prime }$ is determined using the same procedure. The flux at each 8 Hz time point is calculated by multiplying the instantaneous vertical velocity and associated instantaneous O₂ concentration. The changes in fluxes were added over time to produce cumulative flux curves. For three consecutive time intervals (42 to 95 min in length) with undisturbed flux during day and night, slopes of these curves then were calculated to determine light and dark fluxes, respectively. In the following text, fluxes based on the averaged signal of three (3OEC) or two (2OEC) O₂ sensors were termed “3S-flux” and “2S-flux”, respectively. Single sensor fluxes were termed “1S-flux” and uncorrected fluxes “raw” fluxes.

Standard corrections, abbreviated in this text by single letters, were applied to the flux data to reduce potential errors caused by instrument tilt (R), wave effects (W), time shift (T) – caused by spatial separation of sensors (2OEC) and sensor response time – and changes in water column O₂ storage (S) (Berg et al., 2015; McGinnis et al., 2008; Lorke et al., 2013; Huettel et al., 2020). Influence of potential instrument tilt (R) on flux was tested and corrected when necessary through the rotation of the velocity data so that the mean transverse and vertical velocity were nullified (Lee et al., 2004; Lorke et al., 2013; Lorrai et al., 2010). Similarly, wave rotation (W) was rectified by rotating the flow velocity field so that SD(V_y) and SD(V_z) reached a minimum (SD represents 1 standard deviation) (Berg et al., 2015; Berg et al., 2013). Time shifts (T) were rectified through applying time shift corrections to the O₂ data that produced the maximum absolute fluxes (McGinnis et al., 2008; Berg et al., 2015; Berg et al., 2003; Reimers et al., 2016; Fan et al., 1990). Effects of large-scale variations in the average water column O₂ concentration (S) were compensated for through applying an O₂ storage term (J_St = ${\int }_{\mathrm{0}}^h\mathrm{d}C/\mathrm{d}th$, with $\mathrm{d}C/\mathrm{d}t$ being the change in the average O₂ concentration over time, calculated through linear detrending of the measured O₂ data over 15 min intervals, and h = height of the measuring volume) (Holtappels et al., 2013; Rheuban et al., 2014). Furthermore, acceleration or deceleration of current flows can alter the O₂ concentration profile and thereby temporarily modulate vertical flux (Holtappels et al., 2013). Our data analysis indicated that the temporal flux variations caused by transient velocity changes largely cancelled out over time, and a correction for transient velocity changes was not applied.

All recordings in this study are referenced to eastern daylight time (EDT), which is 4 h behind coordinated universal time (UTC−4). All times in the text, graphs, and legends are thus presented in EDT. Daytime was defined as the period between sunrise and sunset. To determine the significance of differences in fluxes measured by the two instruments, the paired t test was utilised. Error margins are reported as ± 1 standard deviation unless stated otherwise.

2.4 Instrument deployments

The 3OEC and 2OEC were deployed at 9 ± 1 m water depth on an exposed backreef carbonate platform in the Florida Keys (24^∘43.523^′ N, 80^∘49.855^′ W; Fig. 1a) on 11, 13, 15, and 16 July 2017. Prior to the deployments, the instruments were synchronised in time. Scuba divers placed the two instruments 10 m apart along a transect perpendicular to the main southwest–northeast flow direction. The tripods were rotated such that the x axis of the ADVs was aligned with the main current direction, and the measuring volumes of the ADVs were adjusted to 35 cm above the average sediment surface level. The seafloor here is covered by highly permeable medium carbonate sand (median grain size: 440 µm; permeability: 3.2 × 10⁻¹¹ ± 1.2 × 10⁻¹² m²) with relatively low carbon content (0.23 % ± 0.05 $\mathrm{\%}\mathrm{sed}.\mathrm{dw}.$) and colonised by microphytobenthos (Chlorophyll a: 4.9 ± 0.1 $\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\mathrm{sed}.\mathrm{dw}.$). During the deployment week, water temperatures averaged 29.9 ± 0.3 ^∘C and salinity 35.0 ± 0.5. Bottom current velocities ranged from 5 to 14 cm s⁻¹. Waves increased from 11 to 15 July, when maximum wave heights of 90 cm were reached, and then dropped again on 16 July. The weather was mostly sunny with some scattered clouds resulting in relatively high light intensities at the seafloor reaching 392 $\mathrm{\mu }\mathrm{mol}\mathrm{photons}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. On each measuring day, the instruments were deployed during daylight time to include the effect of benthic photosynthesis and were retrieved the following day for data download.

Figure 2 Fluxes and cumulative flux trajectories recorded during the deployment week. (a) Daytime and (b) nighttime O₂ fluxes measured by the 3OEC (red circles and line) and 2OEC (black circles and line), as well as significant wave heights (light blue triangles and line) and current flows (green diamonds and line) recorded in July 2017. Note reversed y axis scale for (b) night fluxes. Grey circle in (b) indicates data point compromised by sensor deterioration. Error bars depict standard deviation except for the 11–16 July averages (single circles on right side of panels), for which error bars present standard error. (c) Comparison of cumulative fluxes measured by the 3OEC (left column) and 2OEC (right column) on 11–12, 13–14, 15–16, and 16–17 July (top to bottom) with the respective light (PAR) intensities (orange lines) at the seafloor. Sensor 1: red; sensor 2: blue; sensor 3: green; average sensor signals: thick black lines.

3 Results

3.1 Benthic fluxes

The O₂ fluxes measured by the 3OEC were lower and less variable than those recorded by the 2OEC (Fig. 2a and b). Daytime 3OEC O₂ fluxes averaged 5.2 ± 0.6(SE) $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ and nighttime fluxes −2.8 ± 0.6(SE) $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$, characterising the permeable carbonate sand bed as a site of high carbon turnover and net autotrophy in July 2017. Average 3OEC daytime fluxes were 7 % lower and nighttime fluxes 38 % lower than the respective 2OEC fluxes (day 5.6 ± 0.8(SE), night −3.9 ± 0.5(SE) $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). The difference in the nighttime fluxes between the two instruments was statistically significant (p = 0.04685, p(x≤T) = 0.02342, T = −3.268, DF = 3), while the difference in daytime fluxes was not (p = 0.08077, p(x≤T) = 0.9596, T = 2.5944, DF = 3). The trajectories of the cumulative fluxes were similar in both instruments, but in the 3OEC, the additional sensor and elimination of errors associated with time shift corrections reduced fluctuations of the averaged signal trajectories (Fig. 2c).

Figure 3 The 4 min interval of nighttime data recorded on 13 July (04:53:20–04:57:20) comparing the simultaneous elevation, oxygen, and pressure readings recorded by the 3OEC. Elevation expresses the instantaneous relative elevation of a water parcel that is moved up and down in the water column (see text). (a) Elevation (blue line) and (b–d) the associated O₂ concentrations recorded by the three optodes (red, blue, and green lines) and their (e) average (black line). Brown line in (e) depicts pressure P (scale right y axis) at the height of the ADV. The data were smoothed by a 2.5 s running average. In the absence of a time shift, a minimum in O₂ change occurs when water displacement is zero, and points with zero O₂ change therefore are connected in this graph to points with no elevation change. Since this is a nighttime recording, O₂ minima cross-correlate with elevation maxima. Vertical red lines connect elevation maxima and associated O₂ minima as identified by the OriginLab 2017 software peak-finding algorithm (analysis of 2nd derivative). Vertical grey lines do the same, but in these cases one of the minima or maxima could not be identified by the peak-finding algorithm (manual fit). The instance of the largest temporary time shift (−1.79 s) between elevation maximum and the corresponding 3S O₂ minimum observed within this 4 min interval is indicated by the vertical purple line. The data listed below the graph reveal how the averaging of the O₂ signals reduces the variance of the O₂ signal in the flow measuring volume relative to the individual O₂ signals.

3.2 Time shift

In the 3OEC, elimination of the time shift caused by spatial separation of O₂ and flow measurements does not completely remove time shift errors from the raw fluxes. The remaining time shift errors are caused by the response time of the O₂ sensors (0.2–0.3 s) and temporary distortions of the O₂ concentration field (Fig. 3).

Within the oxygen gradient near the seafloor, vertical water movement associated with wave orbital motion causes O₂ oscillations at a fixed point above the sediment, i.e. at the ADV's flow measuring point. During nighttime, O₂ signal minima occur at the maxima of water parcel elevation as water originating near the O₂-consuming sediment surface is moved up within the water column. Elevation z expresses the instantaneous relative elevation of a water parcel that is moved up and down at the velocity V_z and can be estimated as $z=\int V_z\mathrm{d}t$ (Berg et al., 2015). In the 4 min recording example shown in Fig. 3, nearly parallel vertical connecting lines between the O₂ concentration minima recorded by the three optodes during time intervals with reduced wave activity (e.g. 40–60, 90–120, 210–220 s) confirm that the sensor response times were similar and consistent. Bending in the connecting lines during periods with increased wave activity (reflected by larger pressure oscillations; Fig. 3e) reveals distortions of the O₂ concentration field at the scale of the oxygen sensor spacing (i.e. within the triangle in Fig. 1b). Since the observed time shifts between elevation maxima and the associated O₂ signal minima are positive as well as negative, the bending cannot be attributed to optode response characteristics, implying that it is caused by distortions in the O₂ concentration field. Such distortions temporarily move the 3S signal slightly off centre in the flow measuring volume, producing varying time shifts between 3S signal and velocity data. In this example, a maximum time shift of 1.79 s briefly was reached at t = 183.49 s, lasting less than 6 s. The total time shift caused by sensor response time plus transitory shifts produced by distortion during these 4 min averaged −0.36, −1.52, and −0.20 s for the three sensors, respectively, and −0.59 s for the 3S signal. To compensate for sensor response time, a time shift correction was included in the corrections (STW) used when calculating all 3S-fluxes. The temporary time shifts caused by transitory concentration field distortions largely average out over time as reflected in the 3S variances that were 1.9 to 3.4 times smaller than the 1S variances, and a correction for concentration field distortion was not applied.

Figure 4 Comparison of fluxes based on single sensor signals and average signals, as well as the effects of flux corrections. Error bars represent standard error. (a) Comparison of the daytime (yellow/orange) and nighttime (grey/dark grey) fluxes averaged over all days based on individual sensors and sensor averages (day: red/brown, night: blue/dark blue). (b) Normalised 1S-, 2S-, and 3S-fluxes, colour code as in (a). The 3S-flux (red/blue) was set to 100 % (all data in (a) and (b) are STW-corrected). (c) Effect of corrections on 3S and 2S daytime, nighttime, and 24 h fluxes (light blue) averaged over all days (error bars SE). Raw: not corrected; S: storage-corrected; T: time-shift-corrected; W: wave rotation-corrected; R: rotation-corrected. SW, STW, and STWR are combinations of the above corrections. (d) Normalised differences between the 3OEC STW-corrected average fluxes (set to 100 %) and fluxes with no correction or different corrections recorded with the 3OEC and the 2OEC. Column colour coding in (c) and (d) as listed for (a). Dotted lines allow for the comparison of the fluxes with the 3S-fluxes. Corresponding graphs based on the individual sensor readings are included in the Supplement.

3.3 Uncertainties and flux corrections

3OEC fluxes based on the 3S signal had lower standard errors than fluxes based on individual 1S signals, i.e. standard error was reduced by 27 % ± 35 % (1 SD) during daytime and 114 % ± 158 % (1 SD) during nighttime (Fig. 4a). The normalised daytime 1S-fluxes deviated 3 % (3OEC) and 2 % (2OEC) and nighttime fluxes 18 % (3OEC) and 26 % (2OEC) from the respective normalised 3S-flux and 2S-flux (Fig. 4b). Corrections for storage (S) and time shift (T; in 3OEC to correct for response time) increased raw flux, while corrections for wave rotation (W) and instrument tilt (R) reduced it (Fig. 4c). Applying a combination of storage, time shift, and wave rotation corrections (STW) led to the best agreement between 1S-fluxes, as well as to the strongest enhancement of the raw flux (Fig. 4c and d), as previously found in 2OEC deployments conducted at the same study site (Huettel et al., 2020).

Figure 5 Effect of flow and light on flux. (a) Effect of flow velocity on day- and nighttime fluxes measured with the 3OEC (red circles and line) and 2OEC (black circles and line). The data points indicate the average fluxes calculated for light daytime and dark nighttime periods, separated at 20:00, plotted against the average flow velocity for the respective time periods. Compromised data point from the 16 July 2OEC deployment was excluded from regression (grey circle). Error bars represent standard error. (b) Increase in daytime fluxes with increasing light intensity at the seafloor. Black curves depict photosynthesis-irradiance curves (red circles, thick line: 3OEC; black dots, thin line: 2OEC) calculated using Michaelis–Menten kinetics (P = $P_{max}\left[I\right]/(K_I$ + [I]), with P being the photosynthetic rate at a given light intensity, P_max the maximum potential photosynthetic rate, [I] the light intensity, and K_I the half-saturation constant, i.e. the light intensity at which the photosynthetic rate proceeds at $\mathrm{1}/\mathrm{2}$ P_max). The horizontal black line indicates approximate level of daytime respiration. (c) Light utilisation efficiency of the benthic community based on data shown in (b) (red circles: 3OEC; black circles: 2OEC).

3.4 Effects of waves, unidirectional currents, and light

At our study site, waves were relatively high for this shallow environment (wave height up to 10 % of water depth), and wave orbital motion influenced water movement and pressure near the seafloor during the entire study (e.g. Fig. 3e). Yet, fluxes scaled with the average unidirectional bottom current velocity, which slowed during the deployment week (∼ 30 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ flux increase or decrease per metre per second flow decrease; Fig. 5a), and not with significant wave height (R² < 0.04; Fig. 2a and b) that increased during the study except the last day. On that last day, significant wave heights were nearly identical to those recorded 3 d earlier (Fig. 2a and b; day 0.31 and 0.31 m and night 0.27 and 0.24 m for 13–14 and 16–17 July, respectively), but daytime fluxes decreased by 47 % and sand nighttime fluxes by 58 % between these deployments. Light was ruled out as a cause for the decreases in the fluxes over time because light conditions were similar between deployment days (858 ± 165 mmol photon m⁻² surface PAR for the overlapping time period 17:00–20:00). Improved precision and the generally lower fluxes in the 3OEC were reflected in the community photosynthesis-irradiance (PI) curves (Bernardi et al., 2015). The 3OEC predicted a slightly lower maximum gross benthic primary production (GPP) of 9.9 $\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ (R²: 0.999) than the 2OEC (10.7 $\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$, R²: 0.998; Fig. 5b), as well as a lower light utilisation efficiency (LUE; ratio between GPP and PAR, 3OEC LUE 12.3 % lower than 2OEC LUE at 10 $\mathrm{\mu }\mathrm{mol}\mathrm{photon}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ and 7.4 % lower at 350 $\mathrm{\mu }\mathrm{mol}\mathrm{photon}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$; Fig. 5c). Due to the scatter in the data, these differences in GPP maxima and LUE were statistically not significant.

Table 1 Comparison of 3OEC and 2OEC key characteristics. Approximate costs as of August 2021. Not included are costs for support frame (custom made) and additional sensors (e.g. light, temperature, same for both instruments).

4 Discussion

4.1 3OEC, 2OEC, and advection chamber fluxes

The 3OEC improves benthic flux measurements through the addition of the third concentration sensor, which eliminates errors that can be produced by time shifts between concentration and flow measurements. The averaging of the three instantaneous concentration signals also reduces signal variance and uncertainties that can arise from the disagreement in two of three concentration sensor readings (Table 1). The measuring approach of the aquatic eddy covariance technique – determining fluxes at a distance from their origin – inherently produces data with relatively large variance. This can raise questions regarding their reliability. The fluxes recorded with the 3OEC are validated by the general agreement of the magnitudes and trends of the 3OEC, 2OEC, and benthic advection-chamber-based fluxes (Huettel et al., 2020) all measured at the same study site, as well as benthic fluxes reported by Long (2021). Long's study site with carbonate sands at 6 m water depth off Key Largo (Florida) was close (68.5 km distance) to ours, and the fluxes he recorded in June 2018 reached 5 $\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ during daytime and −3 $\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ during nighttime, similar to the fluxes we measured (5.2, −2.8 $\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$).

4.2 3OEC measurements in dynamic environments

The 3OEC improves benthic flux measurements in dynamic shelf environments. Here waves can produce artefacts in eddy covariance flux measurements (Berg et al., 2015), which in our 3OEC measurements were reduced by the elimination of errors that can be caused by the spatial separation of concentration and flow measurements. Long (2021) proposed positioning the eddy covariance measurement point higher in the water column to decrease flux bias caused by waves. Our 35 cm measuring height was identical to that Long (2021) used for his measurements over Florida Keys' sands and may have contributed to further reducing potential wave artefacts in our measurements. Our findings reveal that horizontal bottom currents dominated benthic flux modulation at our site despite the significant wave action (Figs. 2a and b and 5a) in agreement with results of earlier studies that found an enhancing effect of current on flux in shallow shelf environments with permeable sediment (Berg et al., 2013; Chipman et al., 2016; McGinnis et al., 2014). Continuous flow may be more effective than oscillating flow in driving advective pore water exchange in permeable sediments. In contrast to the steady pressure gradients that drive pore water exchange under continuous unidirectional bottom currents, wave orbital motion produces oscillating gradients which enhance turbulence in the pore space of the sand (Cardenas, 2008; Horton and Pokrajac, 2009; Jouybari et al., 2020). This turbulence and inertial losses associated with the acceleration and deceleration of the pore flows may lessen the effectiveness of the oscillating pressure gradients for driving pore flows and interfacial water exchange.

Figure 6 Change in the cumulative co-spectra for the 3S O₂ flux during the deployment week. Cumulative co-spectra were calculated for hourly intervals from 16:00 to 24:00 for the four deployment periods (0.12 $\mathrm{nmol}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ corresponds to 4.3 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). Colours indicate the time periods for which the co-spectra were calculated, with red, orange, and yellow (warm colours) depicting light periods before sunset (green). Blue and purple colours depict dark periods after sunset.

4.3 The 3OEC facilitates detailed analyses

Reduced uncertainties and higher precision achieved with the 3OEC facilitates more detailed analyses at higher temporal resolution (Fig. 6). This can produce new insights in the processes controlling fluxes at the seafloor. Co-spectra time series, plotted for hourly intervals from 16:00 to 24:00 for our four deployments indicate that during the transition from light to dark (Fig. 6, 16:00–19:00, warm colours), turbulence with a frequency < 0.1 Hz (larger eddies) still contained an upward-directed positive flux signal, while the higher-frequency turbulence (smaller eddies) already carried a downward-directed negative flux signal. As the microphytobenthos photosynthetic O₂ production declined with the decreasing light intensity at the seafloor, flux switched from benthic O₂ release to O₂ uptake. The co-spectra suggest that the ensuing negative benthic flux signal initially was transported by the faster smaller eddies generated at the rough and O₂-consuming sediment–water interface, while the slower large eddies higher in the water column still carried the positive flux signal. The co-spectra document the gradual mixing of the smaller eddies with negative flux signal into the large eddies with positive flux signal, i.e. the negative flux dip in the daytime co-spectra broadened with decreasing light conditions, expanding from the higher to lower frequencies. This eddy memory effect decreased with the general decrease in bottom current velocity during our field campaign as less high-frequency, small eddy turbulence is created at the sediment–water interface at lower flow velocities (Lee and Cheung, 1999; Sleath, 1974). Consequently, the negative dip in the daytime co-spectra disappeared, and the co-spectra appeared almost undisturbed in the last deployment (16–17 July) when bottom currents were low.

4.4 O₂ fluxes in permeable carbonate reef sands

The O₂ fluxes recorded by the 3OEC characterised the coarse carbonate reef sands as sites of intense benthic production and coupled respiration. The nighttime O₂ consumption rates of the coral sand rival respiration rates measured in shallow shelf sediments with much higher organic carbon content (Glud, 2008; Middelburg et al., 2005; Hopkinson and Smith, 2005; Laursen and Seitzinger, 2002) and are within the range reported from other coral reef sands (Cyronak et al., 2013; Eyre et al., 2013; Grenz et al., 2003; Rasheed et al., 2004; Wild et al., 2005, 2004). Since the coral sands at our site are low in organic carbon (< 0.3 % dw) and occur in an oligotrophic subtropical reef environment with low water column chlorophyll and dissolved organic carbon content (NO₃ + NO₂ < 0.2 µmol L⁻¹, NH₄ < 0.5 µmol L⁻¹, PO₄ < 0.05 µmol L⁻¹, Chl a < 0.2 µg L⁻¹, DOC < 200 µmol L⁻¹; Markus Huettel unpublished), a substantial sedimentary source of reduced compounds is required to maintain the observed high respiration rates. Our measurements point to benthic primary production as this source. The compensation light intensity (intensity at which O₂ production exceeds respiration), reached at ∼ 12 $\mathrm{\mu }\mathrm{mol}\mathrm{photons}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, and the high light utilisation efficiency of 0.09–0.10 O₂ per photon near the theoretical limit (0.12 O₂ per photon; Brodersen et al., 2014; Attard and Glud, 2020) indicated that the microphytobenthos could maintain excess production under cloudy conditions, identifying the sedimentary microalgae as the source for the intense organic matter production and export. The estimated maximum production of ∼ 10 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ in the Florida carbonate sands (Fig. 5b) is in line with rates reported for reef lagoon sediments in Moorea (P_max 6.8 ± 0.5 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; Boucher et al., 1998), New Caledonia (P_max ∼ 10 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; Clavier and Garrigue, 1999), and the Great Barrier Reef (P_max ∼ 11 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; Eyre et al., 2013). To put these rates into perspective, eddy covariance flux measurements over dense Mediterranean Posidonia seagrass meadows (13 m depth, PAR 300–400 $\mathrm{\mu }\mathrm{mol}\mathrm{photons}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) revealed daytime O₂ fluxes of 6.8 ± 0.7 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ and nighttime fluxes of −3.6 ± 0.4 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ (Koopmans et al., 2020), i.e. rates of the same magnitude as measured in the microphytobenthos communities. This suggests that the benthic metabolic activity in these shallow oligotrophic environments is largely controlled by light. The trends of nighttime respiration that mirrored those of daytime production (Figs. 2a and b and 5a) indicate that at our site microphytobenthos drove the high O₂ consumption rates through its respiration and by producing highly degradable organic matter that was promptly recycled by the benthic heterotrophic community. Factors contributing to the high microbial activity in the carbonate sands include the high specific surface area of the biogenic grains, their permeability to water and gases, the organic content of the grains, their chemical buffering capacity, and their light guiding characteristics (Marcelino et al., 2013; Huettel et al., 2014; Wild et al., 2006, 2005).

5 Conclusions

The deployments of the 3OEC demonstrate that the new instrument can improve the precision and reliability of benthic flux measurements. 3OEC fluxes in general were smaller, less variable, and had smaller error margins than those produced by the conventional 2OEC eddy covariance instrument that was deployed next to the 3OEC. The advantages of the 3OEC may be most valuable in shallow energetic environments as reflected in the nighttime fluxes recorded by the 3OEC that differed significantly from those measured by the 2OEC. We believe that especially in dynamic settings, the improvements in flux determinations clearly outweigh the downsides associated with the slightly higher complexity of the 3OEC relative to conventional eddy covariance instruments with one or two solute sensors. As summarised in Table 1, the increases in set-up time and costs are modest and may be justified by the improvement of quality and reliability of the flux data that can be achieved with the new instrument (Table 1).

O₂ flux is a key indicator for changes in benthic metabolism and ecosystem health (Glud, 2008), emphasising the need for reliable flux estimates. The aquatic eddy covariance technique arguably is the best available method for measuring flux at the seafloor as it does not alter activities of benthic fauna and flora but integrates effects of patchiness and accounts for the effects flow, light, temperature, and the supply of electron donors and acceptors that affect the fluxes. The increased precision and reliability of the 3OEC data allow for improved modelling and ecological interpretation. In many environmental measuring tasks, a basic, inexpensive instrument can produce data that are relatively close to the “true” values, and the effort and cost for improving the quality of these data typically increase exponentially with gain in data accuracy and precision. The recent developments in affordable optode technology allow a three-optode aquatic eddy covariance instrument to be set up at modest extra cost relative to a conventional two-optode instrument. The 3OEC is an improved eddy covariance instrument that requires less data post-processing and produces flux data of higher quality and reliability. It presents a hardware solution that permits flux measurements also in dynamic shallow shelf environments and is an unmatched instrument to study and improve flux extraction methodologies. We consider the O₂ fluxes produced by this instrument for the permeable reef sands as some of the most realistic flux estimates achievable with present-day technology.

Appendix A

Table A1 Specifications of the PyroScience™ OXR430-UHS retractable oxygen mini sensors used in this study.

Table A2 Specifications of the PyroScience™ FireStingO₂-Mini oxygen meter used in this study.

Table A3 Specifications of the NORTEK Vector acoustic Doppler velocimeter used in this study.