Closing the oxygen mass balance in shallow coastal ecosystems

The oxygen concentration in marine ecosystems is influenced by production and consumption in the water column and fluxes across both the atmosphere–water and benthic–water boundaries. Each of these fluxes has the potential to be significant in shallow ecosystems due to high fluxes and low water volumes. This study evaluated the contributions of these three fluxes to the oxygen budget in two contrasting ecosystems, a Zostera marina (eelgrass) meadow in Virginia, U.S.A., and a coral reef in Bermuda. Benthic oxygen fluxes were evaluated by eddy covariance. Water column oxygen production and consumption were measured using an automated water incubation system. Atmosphere–water oxygen fluxes were estimated by parameterizations based on wind speed or turbulent kinetic energy dissipation rates. We observed significant contributions of both benthic fluxes and water column processes to the oxygen mass balance, despite the often‐assumed dominance of the benthic communities. Water column rates accounted for 45% and 58% of the total oxygen rate, and benthic fluxes accounted for 23% and 39% of the total oxygen rate in the shallow (~ 1.5 m) eelgrass meadow and deeper (~ 7.5 m) reef site, respectively. Atmosphere–water fluxes were a minor component at the deeper reef site (3%) but a major component at the shallow eelgrass meadow (32%), driven by diel changes in the sign and strength of atmosphere–water gradient. When summed, the measured benthic, atmosphere–water, and water column rates predicted, with 85–90% confidence, the observed time rate of change of oxygen in the water column and provided an accurate, high temporal resolution closure of the oxygen mass balance.

Shallow coastal ecosystems lie at the interface between terrestrial and open ocean environments but also contain important internal boundaries at the benthic-water and atmospherewater interfaces. These ecosystems are important sites of biogeochemical cycling and can drive large water column oxygen and carbon concentration changes over short temporal and spatial scales (Bauer et al. 2013). The resulting water concentration changes have direct relevance to critical coastal processes such as primary production, carbon storage, hypoxia, and acidification. However, the relative contributions of fluxes across the benthicwater interface, fluxes across the atmosphere-water interface, and water column metabolic processes are often poorly known.
Oxygen (O 2 ) is often used as a tracer of biogeochemical processes due to its direct involvement in photosynthesis and oxic respiration as well as secondary oxidation of reduced metabolites (Glud 2008). The resulting changes in O 2 concentration can be large due to high metabolic rates relative to the dissolved O 2 concentration in seawater. This makes O 2 an advantageous tracer to examine exchange across coastal ecosystem interfaces compared with other dissolved gases, such as dissolved inorganic carbon, that have much higher total concentrations and a complex reaction chemistry with seawater. With the advent of optical oxygen sensors (Klimant et al. 1995;Tengberg et al. 2006), the measurement of oxygen has become increasingly precise but challenges still exist in accurately quantifying the sources and sinks of O 2 .
The time rate of change of dissolved solutes in the coastal zone can be described by the mass balance of fluxes across coastal ecosystem interfaces and water column metabolism: where C is concentration of a dissolved solute, h is the height of the water column, and the overbars indicate depth-average quantities. The left-hand term represents the overall depth-*Correspondence: mlong@whoi.edu This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Additional Supporting Information may be found in the online version of this article.
integrated change of dissolved solutes through time. The F aw and F bw represent fluxes across the atmosphere-water and benthic-water interfaces, respectively. R wc represents the integrated average rate of water column metabolism. The last term on the right, A Horz , represents the concentration change in a dissolved solute driven solely by horizontal advection and not interface fluxes or water column metabolism. Biogeochemical cycling of O 2 is estimated using a variety of methods that span a wide range of spatiotemporal resolutions and integrate different processes (e.g., F bw , F aw , and R wc ). Onedimensional formulations based on O 2 time-series measurements (i.e., the left-hand side of Eq. 1) were pioneered by Odum (1956), and variations of this method are still commonly applied today (Caffrey et al. 2014;Howarth et al. 2014;Qin and Shen 2019;Tassone and Bukaveckas 2019). These formulations are simple to implement and are easily scaled due to the prevalence of O 2 and depth measurements but provide limited information regarding the different components of the system responsible for the observed flux (Kemp and Boynton 1980;Ziegler and Benner 1998;Gazeau et al. 2005;Demars et al. 2015) and can be biased by physical processes such as advection and mixing (Swaney et al. 1999;Beck et al. 2015).
Advances in techniques and technology have increased the ability to measure different flux components (i.e., the right-hand side of Eq. 1) to better resolve the drivers and spatiotemporal variability of coastal O 2 dynamics. Sediment core incubations, benthic chambers, pore-water profiling, and boundary layer exchange techniques (e.g., Tengberg et al. 1995;Berg et al. 2003;McGillis et al. 2011;Long et al. 2015a) can provide estimates of F bw . A variety of formulations exist for the estimation of F aw from dissolved oxygen concentrations and wind, current, and turbulence parameters (Borges et al. 2004;Zappa et al. 2007;Wanninkhof 2014), and direct estimates have been conducted using floating chambers and eddy covariance techniques (Marino and Howarth 1993;Edson et al. 1998;Long and Nicholson 2018). Methods developed for determining R wc include bottle incubations, isotopic tracers, and automated incubators (Laws et al. 2000;Luz and Barkan 2000;Collins et al. 2018). Each of these different methods has specific biases and ideal conditions for their use, which has been discussed elsewhere (e.g., Tengberg et al. 1995;Berg et al. 2003;Collins et al. 2018;Long and Nicholson 2018).
Here, we quantify the contributions of atmosphere-water (F aw ) and benthic-water (F bw ) oxygen fluxes and water column (R wc ) oxygen production and consumption rates to the O 2 mass balance (Eq. 1) in two shallow coastal ecosystems-an oligotrophic coral reef in Bermuda and a mesotrophic eelgrass meadow in Virginia, U.S.A., both of which have limited nutrient loading and freshwater inputs (Bates et al. 2010;McGlathery et al. 2012). F bw was determined by aquatic eddy covariance, R wc was estimated using a mid-water-depthautomated water incubator, and F aw was parameterized from wind or turbulent kinetic energy dissipation rates. The combination of these techniques yielded closure of the oxygen mass balance in Eq. 1 across an order of magnitude concentration variability and revealed: (1) high spatiotemporal variability in F bw , F aw , and R wc , (2) that R wc can be a large component even in oligotrophic, benthic-dominated systems, and (3) that a combination of high-resolution concentration and turbulence measurements can yield substantially more information than that derived from the time rate of change of concentration alone.

Site description
The study sites for this research were two shallow coastal ecosystems: Hog Reef (32.46 N,64.83

Benthic fluxes
The flux across the benthic-water interface was determined with an eddy covariance system (Fig. 1a) that consisted of an acoustic Doppler velocimeter (ADV, Nortek) coupled to a Fire-stingO 2 Mini fiber-optic O 2 meter with a temperature-compensated, fast-response (< 0.3 s), 430 μm diameter optode (Pyroscience, GE), similar to previous designs (Long et al. 2015a;Long and Nicholson 2018). The mean turbulent O 2 flux was calculated over 0.25 h periods from the product of the instantaneous variations in the vertical velocity and O 2 concentration by: where the prime values indicate the turbulent fluctuating components determined from Reynolds decomposition and the overbar indicates temporal averaging. A flow-through O 2 sensor design used a microfluidic volume (0.33 cm 3 ) and KNF micropump (model NF10, 100 mL min −1 ) that resulted in a quick flush rate (5 Hz) while shielding the optode from ambient light. This design also prevented any disruption of ADV-measured flow rates (Long et al. 2015a) as it was located 2.5 cm behind the measuring volume (see Berg et al. 2015;Donis et al. 2015). The active pumping past the sensors created a constant-flow environment negating zero-crossing velocities and minimized any artifacts related to boundary layer and pressure fluctuations and associated changes in sensor response times (see Reimers et al. 2016). The system logged data from the ADV (three-dimensional velocity), O 2 optode, and inertial measurement unit (IMU) sensor (see below) at a frequency of 16-64 Hz. A pressure sensor on the ADV measured water depth. The height of the velocimeter measuring volume above the sediment surface was chosen using relationships determined by Rheuban and Berg (2013) with measuring heights of 0.8 and 0.5 m over the high surface roughness environments of the reef and eelgrass meadow, respectively Rheuban and Berg 2013). These high measuring heights and the presence of a biological canopy required the use of a storage correction to account for changes in the mean concentration (see Rheuban et al. 2014;Long et al. 2015b).
All instruments were packaged onto a rotating base (e.g., Fig. 1a), which allowed the precise correction for the separation between the sensors (Donis et al. 2015;Holtappels et al. 2015;Reimers et al. 2016) using the known sensor separation, current flow rate, and the fact that sensors were always oriented in line with the flow. The IMU (housing a triaxial accelerometer, gyroscope, and magnetometer) was integrated within the ADV and measured the exact instrument orientation, movement, and acceleration to allow for coordinate matrix transformation that accounts for platform rotation and movement (see Long and Nicholson 2018). The design and motion correction is based on similar approaches used in atmospheric eddy covariance measurements to correct for platform motion (e.g., Edson et al. 1998;Flügge et al. 2016). This new measurement approach including microfluidics (Long et al. 2015a) and IMU integration (Long and Nicholson 2018) allows the EC technique to overcome previous challenges related to platform motion, sensor bias, and sensor separation (Berg et al. 2015;Donis et al. 2015;Holtappels et al. 2015;Reimers et al. 2016). The specific configurations, data treatment, and validation can be found in Long et al. (2015a) and Long and Nicholson (2018). The eddy covariance system was deployed in South Bay, VA, and recorded continuous time series of benthic flux. Periodic gaps in the data, due to low water levels (~1 h), where instruments were exposed at low tide, were linearly interpolated for mass balance calculations (Rheuban et al. 2014). In Bermuda, the eddy covariance system was deployed continuously for 3.5 d without interruption.

Water column incubations
An automated water incubation system (AWIS) was developed to measure hourly rates of oxygen production and consumption in the water column ( Fig. 1b), similar to a recently designed instrument for open ocean moorings (Collins et al. 2018). The system consisted of a clear 150 mL polyethylene terephthalate plastic reaction vessel with a 2 psi check valve mounted at its highest point. The screw cap of the vessel was epoxied in a plastic mount with an optical O 2 and temperature logger (HOBO U26, Onset) and water inlet port located inside the vessel cap. A separate pressure housing contained batteries, a KNF micropump (model NF10, 100 mL min −1 ), an Arduino microcontroller, and real-time clock. The pump inlet had a 500 μm stainless steel inlet filter, and the pump outlet was attached to the reaction vessel inlet port. The oxygen logger recorded temperature-corrected O 2 concentrations in the vessel at 1 min intervals. The microcontroller operated the pump for 5 min each hour, flushing > 3 volumes of the reaction chamber each time. The reaction chamber was changed daily by divers to limit biofouling accumulation on the vessel walls. The AWIS was located at mid-water height and moored via a bungee tether to an earth anchor screwed into the seabed. A second bungee tether attached the incubation chamber to a surface float to provide turbulence-driven mixing of the water within the reaction chamber over each 55 min incubation. Rates of O 2 production or consumption (R wc ) were

Long et al.
Closing the oxygen mass balance determined by linear fit over each 55 min incubation period. The AWIS was deployed continuously at both sites, except for short data gaps (~2 min) due to daily changing of the reaction vessel by divers.
The AWIS system was developed in response to a lack of resolution and consistency from in situ incubations using biological oxygen demand (BOD) bottles (i.e., Collins et al. 2018). Four to eight glass BOD bottles (300 mL) or intravenous (IV) bags (500 mL) were deployed in Virginia and Bermuda, respectively, to compare to the AWIS system periodically. All bottles were filled and sealed underwater. Two bottles were immediately returned to the surface for oxygen concentration determination using a portable Firesting optode. The remaining bottles were left to equilibrate for 2-4 h in situ and then returned to the surface to measure end point oxygen concentration with a Firesting optode.

Atmosphere-water exchange
At the Bermuda reef site, atmosphere-water exchange was determined by the revised wind-parameterization method of Wanninkhof (2014). Although this approach is appropriate for the reef site at the edge of the Bermuda platform that is open to interaction with the wind, it is not appropriate for the Virginia coastal bay that is enclosed by islands, marsh, and banks that limit the influence of wind. Instead, at the Virginia site, we estimated gas exchange rates from the turbulent kinetic energy dissipation rate (ε) based on the methods of Zappa et al. (2003Zappa et al. ( , 2007, as described by Long and Nicholson (2018). The ε was determined hourly from the fit of the vertical velocity wavenumber spectrum in the inertial subrange (Fairall and Larsen 1986;Zappa et al. 2003;Long and Nicholson 2018): where S is the wavenumber spectrum of the fluctuating components of the vertical velocity, κ N is the radian wavenumber, and α is Kolmogorov's empirical constant (= 0.52; see details in Zappa et al. [2003Zappa et al. [ , 2007). The ε was scaled to 35 cm depth to be consistent with parameterizations described by Zappa et al. (2007) that were developed at 35 cm from the surface. The ε was scaled using u * = (εκz) 1/3 and law of the wall (e.g., Reidenbach et al. 2006), where u * is the friction velocity, κ is Von Karman's constant, and z is depth. The hourly gas transfer velocity (k 600 ) was determined from ε using the proportionality relationship determined by Zappa et al. (2007) of k 600 (cm h −1 ) = 0.419(Sc) −0.5 (εv) 0.25 , where Sc is the nondimensional Schmidt number and v is the kinematic viscosity of seawater. The product of the difference between the measured O 2 and the saturation O 2 concentration and the k 600 was then used to calculate F aw (Long and Nicholson 2018).
For Bermuda atmosphere-water fluxes, the 10 m wind speed data were obtained from Sta. BEPB6 on St. George Island, Bermuda, from the National Data Buoy Center. The wind-parameterized Bermuda atmosphere-water O 2 flux was calculated using the k 600 parameterization of Wanninkhof (2014) such that k 600 (cm h −1 ) = 0.251(U 10 ) 2 (Sc/600) −0.5 , where U 10 is 10-m wind speed (m s −1 ). The F aw flux, as above, was then determined by multiplying k 600 by the difference between the measured O 2 and the saturation O 2 concentration.

Supporting measurements
Over each deployment period, a data sonde (measuring salinity, temperature, depth, O 2 , and pH; SeapHOx, Satlantic) and a photosynthetically active radiation (PAR) logger (Odyssey, NZ) were attached to a simple, weighted frame, and placed on the bottom at each site. A four-channel HR-4 spectroradiometer system (HOBI Labs HydroRAD-4) was deployed at the Virginia site to measure down-welling spectral irradiance [E d (λ)] at 6 min intervals throughout the deployment. The spectral energy values (W m −2 nm −1 ) were converted to quanta (μmol photons m −2 s −1 nm −1 ) and integrated across the photosynthetically active region of the spectrum (PAR = P 700 nm λ = 400 nm E d λ ð Þ). The PAR loggers were calibrated against the spectrally integrated PAR estimates provided by the HR-4 (Long et al. 2012). Water column chlorophyll fluorescence at the Virginia site was obtained from a data sonde (EXO2, Total Algae PE Smart Sensor, YSI) immediately adjacent to our sampling site maintained by the Virginia Institute of Marine Sciences (https://stormcentral.waterlog.com/ public/vims).

Daily rates
The daily rates of community gross productivity (P g ) and respiration (R) were determined for continuous time series (24 h) of flux data with the rates weighted by the hours of light and dark as described by Hume et al. (2011). Periods with average light of > 5% of the maximum were defined as light (Hume et al. 2011). Due to the inability to separate respiration and productivity in the measured fluxes, the daytime R rates were assumed to be equal to the observed R rate during the preceding night (Falter et al. 2008;Hume et al. 2011). The presented P g and R rates are therefore conservative estimates, as daytime respiration rates are likely higher than nighttime rates (Glud 2008). The net productivity (P net ) was determined directly from the summation of fluxes over the same 24 h periods.

Uncertainty analysis
Uncertainty estimates for each of the fluxes, parameterizations, and time rate of change measurements were determined using standard error propagation techniques (Glover et al. 2011). The uncertainty for k 600 parameterizations of 20% (Wanninkhof 2014) and 31% (Zappa et al. 2007) were applied in conjunction with the uncertainty of the SeapHOx O 2 sensors (2%) and O 2 solubility parameterizations (2%; Garcia and Gordon 1992). Uncertainty estimates for R wc rates were determined from the error estimates of the linear fits, which were Long et al.
Closing the oxygen mass balance weighted by the HOBO U26 O 2 logger uncertainty (5%). An uncertainty estimate of 20% was applied for the eddy covariance technique as previous studies have determined errors of 10-20% based on both model estimates (Lorrai et al. 2010;Berg et al. 2013) and validation via complementary methods (Berg et al. 2003;Attard et al. 2015).

Results
The instruments were deployed for approximately 4 d at each site (Bermuda and Virginia), and approximately 3.5 d of reliable data were obtained at each site. For all presented fluxes, a positive flux represents an increase in dissolved oxygen in the water and a negative flux is a loss of oxygen from the water. Each of the atmosphere-water ( aw ), water column ( wc ), and benthicwater ( bw ) components is presented in the format of areal fluxes (F aw ,F wc , and F bw ; in units of mmol O 2 m −2 h −1 ) and volumetric rates (R aw , R wc , andR bw ; in units of μmol O 2 L −1 h −1 ), where the circumflex (^) indicates a value that was normalized by the time-varying water depth to allow for comparison between different components (Tables 1-2). The mean hourly flow rates were 0.146 AE 0.070 and 0.016 AE 0.008 m s −1 and the mean water depths were 1.44 AE 0.49 and 7.65 AE 0.34 m for Virginia and Bermuda, respectively (Table 1). The maximum hourly irradiance (PAR) was~3 times higher at the shallow eelgrass site (1902 μE m −2 s −1 ) than at the deeper Bermuda reef (673 μE m −2 s −1 ).

Benthic fluxes
The hourly benthic fluxes (F bw ; Figs. 2-3) show expected diel trends, consistent with photosynthetic production variability due to available irradiance. They resulted in daily gross benthic primary productivity (P g ) of 93.4 AE 16.4 and 111.2 AE 9.7 mmol O 2 m −2 d −1 , benthic respiration (R) of −91.7 AE 18.9 and −121.0 AE 15.1 mmol O 2 m −2 d −1 , and net benthic productivity (P net ) of 1.7 AE 2.5 and −9.8 AE 7.9 mmol O 2 m −2 d −1 , in Virginia and Bermuda, respectively. This yielded an average benthic P g /R ratio of 1.0 for Virginia and 0.9 for Bermuda (Table 2). Hourly benthic P net ranged from 13.0 to −10.0 in Virginia and 11.8 to −8.9 mmol O 2 m −2 h −1 in Bermuda (Table 1; Figs. 2-3).

Water column incubations
The water column rates (R wc ) were determined volumetrically (μmol O 2 L −1 d −1 ) and resulted in mean water column P g of 137.2 AE 34.9 and 14.1 AE 12.8, water column R of −126.5 AE 38.2 and −27.8 AE 13.5, and water column P net of 10.7 AE 6.4 and −13.7 AE 1.1 μmol O 2 L −1 d −1 , in Virginia and Bermuda, respectively ( Table 2). The average water column P g / R ratio was 1.1 and 0.5 for Virginia and Bermuda, respectively. Hourly water column P net ranged from 30.9 to −18.9 in Virginia and 2.9 to −2.8 μmol O 2 L −1 h −1 in Bermuda (Table 1; Figs. 4,5).
The water column rates from the AWIS were linearly correlated to traditional seawater incubations using a York fit (Glover et al. 2011) that incorporates both x and y error (R 2 = 0.54; Supporting Information Fig. S1) with a slope (1.8 AE 0.5) that was not significantly different from a 1 : 1 slope (t = 1.67, p = 0.1252). However, there was considerable variability between the bottle incubations and temporally averaged samples (AWIS).
The R wc rates initially increased with irradiance, but a large reduction in the rate was found at high irradiance coincident with high O 2 concentrations (Figs. 4, 6). When the rates were normalized by irradiance to remove light as the dominant driver, they showed a significant decrease in production with greater O 2 concentration (Fig. 6). This effect was most pronounced at the Virginia site where large diel extremes in O 2 were consistent with high irradiance and high temperatures that occurred at low tide (Figs. 2, 4). For fluxes at night, the opposite trend was observed with greater nighttime respiration during periods of high oxygen (Fig. 6).

Atmosphere-water exchange
The atmosphere-water fluxes (F aw ) were large and varied significantly in Virginia due to the more than an order of magnitude variability in ε and the large difference between the measured O 2 and the saturation O 2 concentrations (Fig. 7). The Bermuda fluxes were an order of magnitude smaller due to consistent and low wind speeds and low variability in the difference between the measured O 2 and the saturation O 2 concentrations (Fig. 8). Hourly fluxes ranged from 22.9 to Table 1. Areal flux and volumetric rate ranges of different flux components with site characteristics. Water column rate (R wc ), benthic-water flux (F bw ), areal water column flux (F wc ), and the combined areal fluxes of gross productivity (P g ), respiration (R), and net productivity (P net ) mean daily values, where the AE values are standard deviations.  18.2%; Fig. 7), and net positive fluxes in Bermuda were consistent with net undersaturation (93.5% AE 6.2%; Fig. 8).

Discussion
Biogeochemical processes in the water column often dominate the oxygen budget in the open ocean because the water volume is large relative to the area of the atmosphere-water and benthic-water interfaces (Glud 2008). In shallow-water systems, fluxes across the atmosphere-water and benthicwater interfaces are likely to have greater relative importance and, in fact, such systems are often described by their dominant benthic communities (e.g., seagrass or coral reef ecosystems). However, water column processes and (in one case) atmosphere-water exchange were important drivers of water column O 2 variability in the subtropical coral reef and a temperate eelgrass meadow studied here. Sustained highfrequency measurements revealed important temporal variability in these drivers and highlight the difficulty of attributing water column time-rate-of-change measurements to specific ecosystem processes without evaluating the full suite of contributing fluxes and rates. Closing the oxygen budget in shallow coastal systems is particularly challenging due to the relatively small volume of water and the potentially large fluxes across both the atmosphere and benthic interfaces.

Benthic fluxes
Benthic P g , R, and P net were on the same order in the Bermuda reef and Virginia eelgrass meadow and showed similar diel variability in benthic flux (~AE 10 mmol m −2 h −1 ). Daily P net indicated a net balance between autotrophy and heterotrophy (P net = 1.7 AE 2.5 mmol m −2 d −1 ) at the Virginia eelgrass site in July. This is consistent with a shift from net autotrophy to net heterotrophy during July found by Rheuban et al. (2014). Similarly, our September P net rates in Bermuda indicated slight heterotrophic conditions (−9.8 AE 7.9 mmol m −2 d −1 ), which is consistent with a shift from net autotrophy to net heterotrophy in September found by Bates et al. (2010).
The eddy covariance benthic flux estimates of Rheuban et al. (2014) are directly comparable to the eddy covariance flux estimates presented here. However, the benthic flux estimates of Bates et al. (2010) were calculated by combining onshore and offshore concentrations with water residence times, with the assumptions that F aw and R wc have a minor impact on daily rates (Bates 2002;Bates et al. 2010). The assumptions that the F aw component was small in Bermuda is consistent with our data, but the Bates et al. (2010) flux estimates are best compared with the sum of the separate benthic and water column rates determined here.

Water column fluxes
Comparison of water column oxygen production and consumption indicated net water column autotrophy at the sensors, depth (a, right), the rate of water column metabolism (R wc , b) based on the slope of line fit to each 55 min AWIS incubation (red lines in c), and the water column chlorophyll concentration (b, right) over a~4 d period. The oxygen time series inside the AWIS (black line, c) and the fit over each 55 min incubation (red, c) where the vertical jumps in the black line indicate the period when the AWIS was flushing with fresh seawater. In addition to negative fluxes at night, the red arrows highlight that negative fluxes were also observed at high PAR, coincident with low tide and strong O 2 supersaturation. eelgrass site and heterotrophy at the coral site ( Table 2). The high rates of R wc relative to F bw and the temporal variability of R wc were both surprising considering the oligotrophic nature of the Bermuda site (Bates et al. 2010) and the low water depth and chlorophyll concentrations at the mesotrophic Virginia site (~2 to 6 μg L −1 ; Zimmerman et al. 2015, McGlathery et al. 2012. Previous results comparing F bw and R wc in similar ecosystems have suggested that R wc was small compared to F bw , even in some eutrophic systems (Ziegler and Benner 1998;Howarth et al. 2014;Long et al. 2015b).
The variability in the AWIS measurements relative to the bottle incubations (Supporting Information Fig. S1) presumably reflects natural temporal changes in the rate on hourly timescales, driven by natural processes and water column parameter variability (e.g., chlorophyll, O 2 , and flow). In contrast, the bottle samples integrated across multiple hours. The AWIS-bottle comparison is further complicated by the fact that longer bottle incubations may suffer from significant "bottle effects" (Zobell 1946;Ferguson et al. 1984;Collins et al. 2018). For example, the bottle samples were typically taken at the beginning of each day and incubated through the periods when we observed a significant decrease in the AWIS rates in conjunction with increasing O 2 (Fig. 6). Likewise, whether bottle incubations were started at low tide or high tide could yield very different rates as a result of the large diel changes in chemical conditions and phytoplankton abundance.
The effect of natural environmental conditions on the AWIS rates was most apparent at the eelgrass site. There was a shift from strong autotrophic R wc in the morning to net heterotrophic R wc around mid-day when high irradiance, low tide, and low flow rates coincided to produce O 2 saturation levels of nearly 170% (Figs. 4, 6). These strong changes are counter to the expected increase in autotrophy through the day with increasing irradiance and may reflect the effect of photorespiration at high O 2 (Lex et al. 1972;Beardall et al. 2003). The nighttime fluxes showed the opposite trends, with higher O 2 consumption rates during periods of high O 2 , consistent with enhanced oxic respiration (Zimmerman et al. 1989). Similar water column patterns were observed at the Bermuda coral site-decreased daytime O 2 production coincident with high O 2 and enhanced nighttime respiration with high O 2 . Although the data do not let us distinguish between various possible causes, they illustrate that the assumption of minimal impact of R wc in systems that appear to be dominated by the benthic community is not always correct. It may be that R wc has been considered small in coastal ecosystems simply due to the lack of temporal sampling resolution and continuous measurements that are needed to accurately capture the variability and dynamics of water column metabolism.

Atmosphere-water fluxes
Atmosphere-water exchange was determined by different methods at each site due to different influences of wind and currents caused by the site characteristics (Long and Nicholson 2018). Using a parameterization based on direct measurements of turbulence is particularly important at the shallow eelgrass site, not only because of the limited influence of wind compared to tides (e.g., Borges et al. 2004) but also because the difference between the measured O 2 and the saturation O 2 concentrations can be particularly large. The large difference between the measured O 2 and the saturation O 2 concentrations in Virginia results in large F aw with substantial diel variability. The largest difference between the measured O 2 and the saturation O 2 concentration occurred at low and high tides, which were followed by maximum flow rates that enhanced turbulence and led to large temporal variability in F aw . During periods of extreme O 2 undersaturation or supersaturation (68-170%), the F aw was twice as large as the F bw . In Bermuda, the low and consistent wind speeds and small difference between the measured O 2 and the saturation O 2 concentrations resulted in F aw that were an order of magnitude lower than that of Virginia. However, the nearly constantly undersaturated water resulted in a net flux of O 2 from the atmosphere.

Mass balance
To examine the closure of the local oxygen mass balance, the time rate change of oxygen in the water column (dO 2 /dt) was compared to that calculated from the ecosystem fluxes. The fluxes F bw and F aw were converted to volumetric rates by integrating the fluxes across the time-varying water depth (R aw and R bw , respectively). The rates were then summed, and the calculated dO 2 /dt (i.e., the right-hand side of Eq. 1) was compared to that of the dO 2 /dt calculated from the data sonde timeseries measurements (lower panels of Figs. 2-3). The summation of rates accurately predicted the observed time rate change of oxygen in the water column at both sites (R 2 = 0.64 and 0.68; Fig. 9) with a normalized root-mean-squared error (NRMSE) of 10% and 15%.
The horizontal advection of O 2 (Eq. 1; A Horz ) was not quantified in this study, and some of the difference between observed and calculated dO 2 /dt (e.g., the model residuals) was likely due to A Horz . At the Virginia eelgrass site, there was a weak correlation between dDepth/dt and the model residuals when the O 2 > O 2mean (R 2 = 0.21) and when O 2 < O 2 mean (R 2 = 0.27), suggesting that advection of nearshore and offshore water, respectively, may be influencing the dO 2 /dt (Supporting Information Fig. S2). Similar comparisons for the Bermuda site produced no significant correlation, likely due to the much smaller dO 2 /dt and the associated uncertainty. In future studies, the Fig. 6. The R wc plotted as a photosynthesis to irradiance curve showing a reduction in rate at high PAR levels (a). R wc normalized by PAR to remove light as the dominant driver and plotted against the oxygen concentration of the seawater (b). Fits are second-order polynomials. At both sites, the normalized values of R wc were reduced at high oxygen levels and were very pronounced in the eelgrass meadow possibly due to photorespiration. Data were excluded from this analysis at < 10% of PAR max . The nighttime values of R wc (c) showed the opposite trend with greater rates during periods of high oxygen, consistent with enhanced oxic respiration.
presented O 2 mass balance could be further refined by directly quantifying A Horz , which can be accomplished by locating oxygen sensing systems along the dominant advection direction to quantify lateral gradients, for example, through Eulerian control volume approaches (Falter et al. 2008;Long et al. 2015b).

Fig. 7.
Atmosphere-water exchange estimated from the turbulent kinetic energy dissipation rate (ε) over an eelgrass bed in South Bay, Virginia. The ε was determined from the ADV data (a, black, left) and was used to estimate the piston velocity by the methods of Zappa et al. (2007) (a, red, right), resulting in the estimated atmosphere-water flux (b). The PAR (orange and red) from the HR4 and Odyssey (Ody) sensors and the oxygen concentration and 100% oxygen saturation (blue and black, respectively) are shown in (c). Fig. 8. Atmosphere-water exchange estimated from the wind speed on Hog Reef, Bermuda. The wind speed (a, black, left) was used to estimate the piston velocity (a, red, right) by the methods of Wanninkhof (2014) to estimate atmosphere-water flux (b). The PAR (orange) and the oxygen concentration and 100% oxygen saturation (blue and black, respectively) are shown in (c). Note different axes scales from Fig. 7.
Although the measured rates (R aw , R wc , andR bw ) combined to yield predicted changes in water column dO 2 /dt that closely matched the observed dO 2 /dt, there are a number of potential sources of bias, in addition to A Horz , that could have resulted in small offsets that were observed. One limitation is the hourly timescale that these methods resolved, which could be a limiting factor in capturing highly dynamic changes in, for example, irradiance, flow, and O 2 concentration. Water column metabolism will vary with light attenuation vertically in the water column and could be a source of bias in the presented measurements at either site due to the high turbidity in Virginia and the greater depth in Bermuda. Finally, at the Virginia site, bubble production on the surface of the eelgrass blades was observed at low tide during high irradiance and could have delayed or prevented the change in measured concentration due to storage of O 2 within the bubble. However, these bubbles were frequently transported directly to the surface and represent a loss of O 2 that is not accounted for in the dissolved O 2 measurements or fluxes (M. H. Long et al. unpubl.) and would not have contributed to the differences between the observed and calculated dO 2 /dt. The summation of rates often led to net rates that were small, for example, during crepuscular periods or during periods of high O 2 undersaturation/supersaturation whereR aw is large and opposite in magnitude to R wc andR bw . Therefore, the percent contribution of the magnitude of each rate component (jR aw j, jR wc j, jR bw j) was divided by the summation of the magnitude of the rates (jR aw j + jR wc j + jR bw j) to examine their relative contributions to the total O 2 rate through time (Fig. 10). Water column production and consumption (R wc ) was the dominant contributor to changes in the total O 2 rate at both the shallow, mesotrophic eelgrass site (45%), and Fig. 9. The observed hourly change in water oxygen concentration (black) and that calculated from the discrete rate measurements (red) through time over an eelgrass bed in South Bay, Virginia (a) and over a coral reef in Bermuda (b). Shading represents oxygen sensor (gray) and rate measurement (red) uncertainty. The PAR (orange and red) from the HR4 and Odyssey (Ody) sensors are on the right (a, b). Panels c and d show the individual atmosphere-water, water column, and benthic rates, where the shading represents rate uncertainty. The correlation coefficient (R 2 ) is derived from a linear fit of the observed and calculated data. The normalized root-mean-squared error (NRSME) is defined as (sum[calculated − observed] 2 )/n) 0.5 , where n is the number of hourly flux values.
deeper, oligotrophic reef site (58%). Benthic rates (R bw ) contributed 23% and 39% of the total O 2 rate at the Virginia and Bermuda sites, respectively. The importance of atmospherewater gas exchange (R aw ) in Virginia (32%) was an order of magnitude greater than in Bermuda (3%). The large contribution of theR bw to the total O 2 rate was expected as reefs and eelgrasses are often viewed as dominant primary producers.
TheR aw contribution scales with the difference between the measured O 2 and the saturation O 2 concentration and was thus not surprising that theR aw was large and variable in Virginia and small and positive in Bermuda based on the measured differences. The high rates of R wc in Bermuda were surprising due to our supposition that low-nutrient oligotrophic water (Bates et al. 2010) would not support high rates of planktonic metabolism. However, the rates were an order of magnitude lower in Bermuda than Virginia (~AE 3 vs. AE 30 μmol L −1 h −1 ), which reflects the higher chlorophyll concentrations typically found at the Virginia site (~2 to 6 μg L −1 ; Zimmerman et al. 2015;McGlathery et al. 2012). The similar importance of R wc to the total O 2 rate in both Virginia and Bermuda (~45% to 60%) reflects the greater water column depth at Bermuda (7.5 m vs. 1.5 m in Virginia), which reduced the relative influence ofR aw andR bw in Bermuda.
At the shallow eelgrass site, there was considerable variability through time in the relative importance of each rate component with highR aw during times of extreme O 2 saturation, highR bw during peak light, and high R wc in the morning. At a specific time, any one of the three different rates could represent 70% or more of the total O 2 rate. For example, when the AWIS flushed with fresh seawater during periods of high water column O 2 saturation (i.e., the vertical jumps in the black line of Fig. 4c), the R wc rates that followed were negligible or negative, possibly indicating photorespiration (Lex et al. 1972;Beardall et al. 2003). Therefore, the R wc did not contribute to the large increases in O 2 concentrations observed during these periods. AsR aw is negative during these periods of O 2 Fig. 10. The absolute contribution of each rate component (%) to the total rate (jR aw j + jR wc j + jR bw j) (a, b), the mean % rate contribution (c, d), and the daily net rate (e, f) of the atmosphere-water (red), water column (blue), and benthic (green) rates from the Virginia eelgrass meadow (left) and Bermuda reef (right). Differences between net daily rates were evaluated by analysis of variance with an asterisk representing significant differences determined by Tukey posttests.
supersaturation, these increased O 2 concentrations were driven largely by eelgrassR bw , which is apparent from its contributions of the total O 2 rate (Fig. 10a;at hour 192,217,and 242). These periods of high eelgrass production (e.g., Fig. 2) coincided with low tide, which allowed F bw to have a much greater effect on the O 2 concentration of a very shallow-water column. Despite this variability in the importance of the individual rates through time, the mean % of total O 2 rate revealed similar magnitudes between the three components, and the net daily rate revealed no significant differences between the three components, indicating the first-order importance of each rate component at the eelgrass site (Fig. 10). Overall, the combined net areal productivity (F bw +F wc ) indicated net autotrophy at the Virginia eelgrass site (4.9 AE 2.7 mmol m −2 d −1 ) driven by both net positive benthic fluxes and net positive water column fluxes.
In Bermuda, the rate partitioning was comparatively stable withR aw being less than~10% of the total O 2 rate and large contributions from eitherR bw or R wc that could represent up to 90% of the total O 2 rate (Fig. 10). Rates of R wc were again higher in the early morning hours, whereasR bw more closely followed the available light. The mean % of total O 2 rate of the three components was more variable in Bermuda than in Virginia, with the greatest contribution from R wc . The net daily rates in Bermuda were dominated by water column heterotrophy, with small contributions fromR bw andR aw . Overall, the combined net areal productivity (F bw +F wc ) indicated net heterotrophy at the Bermuda reef site (−116.8 AE 4.4 mmol m −2 d −1 ) driven largely by net negative water column fluxes.

Conclusions
In many shallow coastal environments, the combination of high metabolic rates, a fully photic water column, and long residence times leads to drastic diel changes in water column conditions. The presented techniques and rates provide for a high temporal resolution closure of the oxygen mass balance for coastal ecosystems. The measured rates accurately predicted the observed water column O 2 changes with 85-90% confidence at two different coastal sites that had a 10-fold difference in diel O 2 variability. These results highlight that benthic, water column, or atmosphere-water rates all have the potential to individually drive the majority of observed O 2 concentration change at short timescales and that these rates are highly dependent on the biological communities present and their biomass relative to the total water volume. Changes in mean water column concentrations reflect multiple factors, even at shallow-water sites sometimes thought of as benthic-dominated. As a result, it is difficult to partition primary production between the water column and the benthos without measuring these rates separately. Evaluating the contribution of these fluxes to changing water conditions has the potential to greatly enhance our ability to understand ecosystem dynamics and how they relate to critical processes such as primary production, carbon storage, acidification, and hypoxia.