The estimation of gross oxygen production and community respiration from autonomous time‐series measurements in the oligotrophic ocean

Diel variations in oxygen concentration have been extensively used to estimate rates of photosynthesis and respiration in productive freshwater and marine ecosystems. Recent improvements in optical oxygen sensors now enable us to use the same approach to estimate metabolic rates in the oligotrophic waters that cover most of the global ocean and for measurements collected by autonomous underwater vehicles. By building on previous methods, we propose a procedure to estimate photosynthesis and respiration from vertically resolved diel measurements of oxygen concentration. This procedure involves isolating the oxygen variation due to biological processes from the variation due to physical processes, and calculating metabolic rates from biogenic oxygen changes using linear least squares analysis. We tested our method on underwater glider observations from the surface layer of the North Pacific Subtropical Gyre where we estimated rates of gross oxygen production and community respiration both averaging 1.0 mmol O2 m−3 d−1, consistent with previous estimates from the same environment. Method uncertainty was computed as the standard deviation of the fitted parameters and averaged 0.6 and 0.5 mmol O2 m−3 d−1 for oxygen production and respiration, respectively. The variability of metabolic rates was larger than this uncertainty and we were able to discern covariation in the biological production and consumption of oxygen. The proposed method resolved variability on time scales of approximately 1 week. This resolution can be improved in several ways including by measuring turbulent mixing, increasing the number of measurements in the surface ocean, and adopting a Lagrangian approach during data collection.

When the sun shines on the Earth's surface, photosynthesis produces the organic matter that fuels the metabolism of most organisms, from bacteria to human beings. This organic matter is then utilized mostly through aerobic respiration, a key metabolic process shared by most living forms. In different ways, the rates of photosynthesis and respiration define the rates of energy processing of an ecosystem, which is its pace. For this reason, estimates of gross photosynthesis and respiration are particularly valuable, but they require a considerable investment of time and effort. This is particularly true in open ocean environments where photosynthesis and respiration are traditionally estimated by incubating natural waters for several hours and by measuring either the organic carbon synthesized from inorganic carbon (Steemann Nielsen 1951;Fitzwater et al. 1982) or the oxygen, O 2 , produced and consumed by organisms (Gaarder and Gran 1927;Bender et al. 1987;Ferrón et al. 2016). Both these approaches are time consuming, and subject to potential biases due to the enclosure of natural communities in a bottle. Furthermore, questions remain unanswered about the temporal and spatial variability of photosynthesis and respiration on scales of days to tens of days, and kilometers to hundreds of kilometers, due to the scarcity of incubation-based measurements. As a consequence, it is desirable to find alternative methods to obtain more estimates of photosynthesis and respiration, especially at fine temporal and spatial resolution. Here we explore the promising approach of using autonomous underwater vehicles to measure diel changes in O 2 concentration.
Already in 1930, Butcher et al. reported that the photosynthetic production of O 2 during daytime and its continuous consumption through respiration resulted in diel oscillations in O 2 concentration in three British rivers (Butcher et al. 1930). Years later, similar observations were collected in a productive coral reef ecosystem, where diel O 2 oscillations were used for the first time to obtain quantitative estimates of gross O 2 production (GOP) and community respiration (CR) (Sargent and Austin 1949). The procedure to calculate metabolic rates from diel O 2 changes was explicitly formulated by Odum (1956) who also applied this method to coastal marine environments (Odum and Hoskin 1958). Since then, measuring the in situ variation of O 2 concentration has become a widespread method to estimate metabolic rates in lakes (reviewed in Staehr et al. 2010), coral reef ecosystems (Barnes and Devereux 1984;Gattuso et al. 1993), and estuaries (Kemp and Boynton 1980;D'Avanzo et al. 1996;Nidzieko et al. 2014). The application of this method has been more challenging in oligotrophic ocean gyres that cover most of the planet because of the low photosynthetic and respiratory rates characteristic of these environments. For example, a typical rate of GOP in the North Pacific Subtropical Gyre (NPSG) on the order of 1 mmol O 2 m −3 d −1 is associated with a diel peak amplitude in O 2 concentration of only 0.5 mmol O 2 m −3 (assuming that the system is near steady state and a 12 h day length) which is two orders of magnitude lower than the baseline O 2 concentration of~200 mmol m −3 (Williams et al. 1983). Despite the low signal in these ecosystems, diel O 2 oscillations in oligotrophic environments have now been observed using a variety of methods including Winkler titrations (Tijssen 1979;Johnson et al. 1983;Williams and Purdie 1991), O 2 to Ar ratios , and optical O 2 sensors . The use of optical sensors (i.e., optodes) is of particular interest here as they are compact enough to be mounted on autonomous underwater vehicles and their precision has now been shown to be sufficient to resolve diel oscillations smaller than 1 mmol O 2 m −3 . For example, the O 2 concentration that we measured using an optode mounted on an array drifting in the NPSG showed clear diel oscillations with peak amplitude of about 0.5 mmol O 2 m −3 (Fig. 1). In the near future, the combined use of autonomous underwater vehicles and optical O 2 sensors could dramatically increase the number of observations of O 2 diel cycles in the global ocean. Specifically, underwater gliders are commonly used in different regions of the ocean (Rudnick 2016), where they collect multiple O 2 measurements per day. These autonomous observations have lower temporal resolution (~10 per day) than measurements collected with sensors placed at a fixed depth (such as those usually employed in freshwater and coastal studies, or in the example of Fig. 1) so they require a different method of data analysis to estimate GOP and CR. In particular, Odum and Hoskin (1958) proposed a method to calculate GOP from the integral of the O 2 rate of change during daytime without making assumptions on the diel variability of this rate. However, when diel O 2 changes are low and close to the variability due to instrument noise, this approach becomes very sensitive to O 2 measurements collected near sunrise and sunset. For this reason, the method that we propose herein is based on the assumption of a fixed diel variation in the rate of O 2 production, as done by Nicholson et al. (2015). Their study further assumed that oxygen production was balanced by respiration on a daily scale; however, changes in photosynthesis can be independent from changes in respiration. Therefore, in the new method that we propose, we obtain separate estimates for GOP and CR. Furthermore, in contrast to previous efforts, we propose to use the vertical information from glider measurements to exclude temporal O 2 variations due to entrainment at the base of the mixed layer. As done by Nicholson et al. (2015), we only analyzed O 2 variability in the surface layer because we were not able to measure marked diel O 2 cycles in deeper waters (below~40 m, on average), likely due to the lower rates of photosynthesis, compounded by higher physically driven variability. In what follows, we describe an idealized procedure to separate physical and biological contributions from the variability of O 2 in the surface layer of the ocean, and we formulate a method to estimate GOP and CR once the biological variability has been isolated. In many cases, the observations collected with autonomous vehicles will differ from the ideal observations required in our conceptual model. To account for this discrepancy, we describe some simplifications and procedures to facilitate the accurate retrieval of metabolic rates from realistic sampling scenarios. We then show results from a case study in which we applied our method to the observations from four underwater glider missions in the NPSG. Last, we discuss the validity of the proposed method and analyze how it can be improved for use on future autonomous O 2 measurements.

Materials and procedures
Physical and biological variability of O 2 in the surface layer and the ideal estimation of GOP and CR from O 2 time-series The mean concentration of O 2 dissolved in a column of near-surface seawater that is followed in a quasi-Lagrangian fashion is modified by at least four processes: (1) phytoplankton photosynthesis; (2) respiration by autotrophic and heterotrophic organisms; (3) the flux through the base of the column via entrainment or diapycnal mixing; and (4) the exchange of O 2 across the air-sea interface (Fig. 2). Consequently, in order to estimate rates of GOP and CR, the first step is to separate the variation of O 2 due to physical processes (entrainment, mixing, and air-sea flux) from the variation due to biological processes.
The O 2 variability due to entrainment can be eliminated by suitably defining the surface layer in which we analyze the diel time-series of dissolved O 2 . Specifically, we can average the concentration of O 2 above an isopycnal surface that is always deeper than the base of the surface mixed layer, Z ML . When choosing this isopycnal, one should consider that GOP in oligotrophic environments tends to decrease with depth and so does the peak amplitude of diel biological O 2 oscillations. As a consequence, in order to be able to measure diel oscillations, we suggest isolating O 2 variability in a surface layer that is as thin as possible. In our analysis of glider observations, we defined the surface layer as the water above the maximum potential density measured at Z ML on a day-to-day basis. The depth of the isopycnal surface defining the surface layer is here referred to as Z SL , and the average O 2 concentration in the surface layer, O SL , is defined as: Fluxes of O 2 between the ocean surface and the atmosphere include a diffusive component and a component due to bubble dynamics. The total flux, F atm (mmol O 2 m −2 d −1 ), can be parameterized using the approach proposed by Liang et al. (2013) when knowing the concentration of O 2 dissolved in surface waters, wind speed, atmospheric pressure at the sea surface, temperature, and salinity. Here we use the sign convention that F atm is positive for downward fluxes (from air to water).
Diapycnal O 2 fluxes at the base of the surface layer, F diff , are the product of the vertical concentration gradient at Z SL , dO 2 /dz (mmol O 2 m −4 ), and vertical eddy diffusivity, K z (m 2 s −1 ). When both these quantities are estimated, we can calculate the diffusive flux and remove its contribution from the O 2 variability in the surface layer. We define positive diapycnal fluxes as those increasing the concentration of O 2 dissolved in the surface layer.
After estimating the O 2 fluxes at the sea surface and at the base of the surface layer, the O 2 time-series can be corrected for their contribution scaled by the daily average depth of the surface layer, Z SL : where O bio is the time-series of O 2 after correcting for physical processes, t is time, and t 0 is the time of the first point of the time-series. Once the variability of O 2 concentration due to biological processes has been isolated, we can separate the contributions of photosynthesis and respiration by assuming a diel shape for the rate of each process. Here we assume that photosynthesis is linearly proportional to light intensity, E, and that respiration is constant throughout the day (the rationale for both choices is discussed in a following section). As a result, O bio (t) becomes: We can then find the three parameters (O 0 , GOP, and CR) that provide the best approximation to the solution of this system using linear least squares analysis. The uncertainty of the solution can be estimated by calculating the variancecovariance matrix of GOP and CR by bootstrapping the residuals (in our examples below we used 200 iterations) (Draper and Smith 2014). The code in MATLAB ® to solve Eq. 4 and to estimate parameter uncertainty through the variance-covariance matrix is provided in the Supporting Information. One numerical consideration to take into account when solving Eq. 4 is that the last two columns of the coefficient matrix are negatively correlated. The stronger this correlation, the larger the covariance in the estimates of GOP and CR, meaning that the two estimates are more dependent on each other; for example, measurement noise in an oxygen sample will tend to cause both GOP and CR to be overestimated, or both underestimated. In a numerical simulation of the fit outcome, we found that the choice of the starting time t 0 during the day affects the quality of the fits (as detailed in the Supporting Information).

Glider and ancillary observations
Observations were collected using three underwater gliders (Kongsberg/University of Washington Seagliders s/n 146, 148, and 512) during four separate missions in the open ocean north of the islands of Oahu and Maui. The first two missions, henceforth A and B, sampled a mesoscale anticyclonic eddy in summer 2015 while missions C and D sampled a cycloneanticyclone eddy dipole in spring 2016. Periods and coordinates of glider sampling are reported in Table 1 along with the number of profiles collected during each mission. The duty pattern varied among and within missions, with maximum depths of 500, 700, or 900 m. The average dive time was 4.1 AE 1.1 h (average AE standard deviation), and the average horizontal distance traveled was 2.8 AE 1.4 km. Gliders were equipped with temperature and conductivity sensors (Sea-Bird CT Sail) and O 2 optodes (Aanderaa 3830 or 4330). Optode sampling intervals were set to 5 s in the upper 30 m, but they were increased to 10 s in the 30-200 m range in order to decrease power requirement and increase the length of glider missions. Gliders descended at an average vertical speed of 0.11 m s −1 with an average pitch angle of 19 and ascended at an average speed of 0.21 m s −1 with an average pitch angle of 20 . As a result of the different vertical velocities during ascent and descent, the average number of O 2 measurements (from a depth of 2 m to Z SL ) used to calculate a single O SL data point were 27 AE 19 during descent and 56 AE 24 during ascent. The original hydrographic and optode measurements were binned on a 2-m spaced vertical grid before subsequent analyses.
Mixed layer depth, Z ML , was calculated from hydrographic measurements as the first depth where the potential density was at least 0.03 kg m −3 larger than the value at 10 m (de Boyer Montégut et al. 2004). During each day, we took the largest potential density value at Z ML and used it as the value defining the base of the surface layer, Z SL .
The optodes used in this study have a response time of about 30 s that causes a mismatch between O 2 profiles measured during descent and ascent (Supporting Information Fig. S2). We corrected for this offset by using an inverse filtering algorithm specifically developed for optode sensors (Bittig et al. 2014;Bittig and Körtzinger 2017;Supporting Information). Inverse filtering proved very effective in reducing the differences between concentrations measured during ascent and descent (Supporting Information Fig. S2b) although we report a residual offset in the upper water column (Supporting Information Fig. S2c). To overcome artifacts in the estimates of GOP and CR due to this residual mismatch, we modified the ascent data by subtracting its median daily value and adding the median descent daily value.
While optode performance is stable during glider deployments, optode storage before deployment has been documented to affect the accuracy of O 2 measurements (Bittig We used satellite estimates of wind speed at 10 m above the sea surface from microwave backscatter. Observations were obtained from the Blended Sea Winds data product (Zhang et al. 2006) that merges the observations from several satellites and is produced by the National Oceanic and Atmospheric Administration's National Climatic Data Center (NOAA-NCDC). Sea-level pressure was retrieved from the meteorological reanalysis by the National Centers for Environmental Prediction and the National Center for Atmospheric Research (Kistler et al. 2001). Both wind speed and sea-level pressure data have a temporal resolution of 6 h and a horizontal resolution of 2.5 both in latitude and in longitude (277 and 255 km, respectively, for our study area). Wind speed and sea level pressure were used together with glider measurements to compute the O 2 flux at the sea surface. Downward photosynthetically available irradiance (PAR) at the sea surface was obtained as the data product distributed by the National Aeronautics and Space Administration (NASA) Ocean Biology Processing Group (OBPG). PAR was distributed with a 4 km resolution as daily values calculated from instantaneous satellite irradiance measurements (Frouin et al. 2012). We averaged the observations from three satellites: MODIS Aqua, MODIS Terra, and VIIRS.
The magnitude of the flux of O 2 due to turbulent mixing at Z SL was estimated by computing the vertical gradient of O 2 using central finite differences, and by assuming a K z in the range 10 −5 to 10 −4 m 2 s −1 (Hamme and Emerson 2006).

Estimates of GOP and CR from glider measurements
Daily GOP and CR were calculated for each glider mission by selecting all O 2 measurement inside a time window of 1.2 d, chosen so that most fits span at least 24 h. As a result, the average period between the first and the last measurement used for a fit was 1.1 d, with less than 10% of the values <1 d. The time of the first measurement used for the fit was on average close to 9:00 am local time consistent with the optimal time-window estimated for the region and period of our measurements (Supporting Information). For each day, the diel O 2 anomaly was defined as O SL -O SL, where O SL is the average O SL in the daily time window selected for the fit. Some daily time-series of O SL included episodic unrealistic O 2 concentrations (whose causes were not identified), which were discarded by removing observations that exceed a difference of three median absolute deviations from the median daily O SL . Table 2. Average AE standard deviation of daily observations from the four glider missions reported in separate columns. O sat is oxygen saturation; T SL is temperature in the surface layer; u10 is wind speed at a height of 10 m above sea level. F atm and F diff are the areal fluxes of O 2 in the surface layer from air-sea exchanges and diapycnal mixing, respectively.  As gliders did not measure E in situ, we assumed that it varied as in cloud-free conditions. This means that E(t) was modeled as a function of solar elevation, θ(t), and of the irradiance value when the sun was directly overhead, E 0 (mol photons m −2 s −1 ): where E(t) was set to 0 for negative values of solar elevation. Gliders did not sample in a Lagrangian fashion so we used the day-to-day change in O 2 concentration in the surface layer as an indicator of the degree of horizontal variability present in the O 2 time-series. During the time-window used for the fit, gliders traveled a horizontal distance of 16.2 AE 7.3 km while collecting 13.2 AE 3.7 vertical profiles (henceforth data are reported as mean AE standard deviation unless stated otherwise), half during ascent and half during descent. To minimize the impact of horizontal O 2 changes on rate estimates, we discarded estimates for days when we observed day-to-day O 2 changes above 0.75 mmol m −3 d −1 , which were unlikely to be caused by ecological dynamics. During these days, the absolute value of the day-to-day O 2 change was significantly correlated with the absolute value of the day-to-day salinity change (Spearman ρ = 0.6). The correlation with seawater hydrography seems to confirm that day-to-day changes in O 2 concentration larger than 0.75 mmol m −3 d −1 were linked to spatial variations rather than to ecological dynamics.
The variability due to diapycnal flux was not subtracted from the O 2 time-series measured by gliders because K z is unknown. The diapycnal flux was estimated crudely by assuming the range of K z for the base of the surface layer proposed by Hamme and Emerson (2006) of 10 −5 to 10 −4 m 2 s −1 . We then compared the magnitude of this flux to our estimates of GOP and CR to understand if the lack of correction for diapycnal fluxes caused important biases in our estimates. We found that the average diapycnal flux caused an O 2 change in the surface layer equal to 2-19% of the average metabolic rate, for K z equal to 10 −5 and 10 −4 m 2 s −1 , respectively. The largest value for these estimates accounted for less than half of the average rate uncertainty so it is unlikely that diapycnal fluxes caused strong biases in the estimates of metabolic rates. However, strong vertical O 2 gradients associated with strong diffusivities likely affected fit quality.
Rates with different uncertainty were averaged using weights, w i , defined as the inverse of the sum of the variances along the major and minor axes of the uncertainty ellipse of GOP and CR, which were then normalized so that their sum equaled 1. The compounded variances were obtained using the formula: where X is either GOP or CR, and X is its weighted average. Similarly, the covariance between GOP and CR, was obtained using the formula: where X and Y indicate GOP and CR. The standard error of weighted average rates was computed by bootstrapping the weighted average over 1000 iterations and by calculating the standard deviation of these averages.
Fit quality was assessed using: (1) the uncertainty of GOP and CR from bootstrapping the residuals; (2) the p value of the linear correlation between the fitted O 2 time-series and O bio ; and (3) the p value of the Durbin-Watson test (Durbin and Watson 1950) to estimate the degree of autocorrelation in the residuals.
As an example of the procedure for deriving estimates of GOP and CR from glider measurements, we report the different steps for 29 March 2016, during glider mission C (Fig. 3). As described in a previous section, we started from vertically resolved O 2 measurements (Fig. 3A) and ended with a corrected surface layer O 2 time-series and the curve that approximates it (Fig. 3B, Supporting Information Table S1).

Assessment and discussion
Dissolved O 2 variability and metabolic rates from glider measurements The concentration of O 2 from aggregate observations of all glider missions during a diel cycle had a minimum near sunrise and a maximum near sunset (Fig. 4A), as would be expected due to daytime photosynthesis and nighttime respiration. O 2 increased most rapidly in the late morning hours (Fig. 4B). O 2 concentration and temperature both increased toward the end of the day leading to higher saturation and O 2 diffusion from the ocean to the atmosphere around sunset. The average air-sea flux was positive (O 2 influx) approximately between 2:30 and 9:30 am, when O 2 concentration was lowest. The air-sea flux was negative (O 2 efflux) during the rest of the day, but its average contribution to the rate of O 2 change in the surface layer was small (Fig. 4B, Supporting Information). The average diapycnal mixing flux at Z SL was positive (O 2 influx) throughout the day (Fig. 4B) and its contribution to O 2 changes in the surface layer averaged 0.02 and 0.19 mmol m −3 d −1 , if considering a K z of 10 −5 and 10 −4 m 2 s −1 , respectively. In terms of areal fluxes, the air sea exchange averaged −0.6 AE 7.4 mmol m −2 d −1 , whereas the diapycnal flux was 0.8 AE 0.9 mmol m −2 d −1 if considering a K z of 10 −5 m 2 s −1 , and 10 times higher if considering a K z of 10 −4 m 2 s −1 . Considering the importance of the variability of these areal fluxes, we analyzed the sensitivity of the air-sea flux to changes in O 2 saturation and wind speed, and the sensitivity of the diapycnal flux to changes in K z and O 2 gradients (Supporting Information).
We obtained estimates of GOP and CR from a total of 211 d of glider observations (sensor malfunctioning affected O 2 measurements during 10 d), but discarded 9% of the estimates when the day-to-day change in surface layer O 2 concentration exceeded 0.75 mmol m −3 . We also discarded results from 1 d that represented an outlier in the distribution of CR value, with an estimated rate of −9 mmol m −3 d −1 . Among the remaining 191 fits, residual autocorrelation was present in 21% of the cases (p < 0.05). The weighted averages of these fits resulted in the same value for GOP and CR of 1.0 AE 0.7 mmol m −3 d −1 (weighted mean AE weighted standard deviation). The uncertainty of the single estimate (computed as the rate standard deviation obtained by bootstrapping the residuals) averaged 0.6 and 0.5 mmol m −3 d −1 for GOP and CR, respectively. The average fit uncertainty accounted for 72% and 50% of the GOP and CR weighted variance measured among daily estimates. Rates of GOP and CR were significantly correlated (R 2 = 0.6) and quantitatively similar indicating coupling between the production and consumption of O 2 in the surface layer (Fig. 5A). The covariation of GOP and CR was not an artifact due to the covariance of the two parameters obtained from the least squares approach. When subtracting the average variance-covariance due to fit uncertainty from the weighted variance-covariance matrix from the daily estimates we obtained a new estimate of rate variability whose major axis had a slope of 1.6, indicating lower variability in GOP than in CR (Fig. 5B).   Table 2). The concentration of O 2 dissolved in the surface layer was larger in spring than in summer, but O 2 saturation with respect to the equilibrium with the atmospheric partial pressure was similar among different missions (Table 2). This indicates that changes in surface O 2 concentrations among glider missions were due to changes in O 2 solubility driven by temperature (Table 2). In summer 2015, there was a strong subsurface O 2 maximum and a steep vertical O 2 gradient at the base of the surface layer, whereas in spring 2016 the O 2 gradient was very small (Table 2). This excess O 2 beneath the surface mixed layer was a result of the seasonal accumulation due to net O 2 production (Shulenberger and Reid 1981;Riser and Johnson 2008) that led to larger estimates of the diffusive flux in summer than in spring ( Table 2). The variability of GOP and CR was also larger in summer than in spring and this was at least partly due to a larger fit uncertainty (Table 2).
To assess the importance of different factors on the quality of our fits, we calculated the Spearman correlation coefficient between the uncertainty of the estimate of GOP from the fit and: (1) daily surface PAR on the glider position from satellite observations; (2) the absolute value of the daily average of the volumetric air sea flux, 1 ZSL j F atm j; 3) diapycnal mixing as represented by the absolute value of the O 2 gradient at the base of the surface layer divided by the depth of the surface layer, 1 Z SL j F diff j; and (4) horizontal O 2 changes as represented by the absolute value of the day-to-day change in O 2 concentration in the surface layer. The only significant factors (p < 0.05) affecting fit uncertainty were diapycnal mixing (ρ = 0.4) and air-sea flux (ρ = 0.2). The same correlations were obtained if the uncertainty of CR (instead of GOP) was used as a proxy for fit quality, indicating the robustness of these results. The positive correlation between rate uncertainty and diapycnal mixing indicates that the lack of a correction for mixing in surface O 2 time-series increased fit uncertainty. This finding is consistent with the observation of higher fit uncertainty in summer than in spring, due to changes in vertical O 2 gradients (Table 2). We notice that if diapycnal mixing accounted for a constant O 2 change throughout the day, it would not affect fit quality, but it would bias the estimates of CR. Consequently, O 2 change due to diapycnal mixing must have varied during the day, possibly due to the diel cycle of turbulent mixing in the surface layer (Brainerd and Gregg 1995).

Validity of the estimates of GOP and CR and statistical considerations
Gross photosynthesis can be quantified using the rate of several steps in the photosynthetic process including photon absorption, electron flow, O 2 production, and organic carbon synthesis. Each approach provides a different answer to the quantification of photosynthesis as not all photons are used for photochemistry, not all electrons come from water splitting, and not all reductants are used for the reduction of inorganic carbon. For this reason, we validated our method by comparing our estimates of GOP and CR to previous estimates exclusively based on O 2 production or consumption. Even so, O 2 -based techniques use different approaches to measure rates of photosynthesis and respiration. A thorough assessment of these different methods is beyond the scope of the present study, but previous method comparisons have reported differences between GOP values obtained using incubations and values obtained without incubations (Williams and Purdie 1991;Juranek and Quay 2005;Quay et al. 2010). Nonincubation methods generally measure greater GOP than incubation-based methods and this could be due to factors including: (1) the enclosure of natural communities in finite volumes with solid boundaries that affect their photosynthetic potential in vitro (Gieskes et al. 1979); (2) the difficulty of reproducing in vitro the in situ environmental conditions (e.g., Marra 1978); (3) the requirement for a correction to nonincubation estimates to account for physical O 2 fluxes such as entrainment, mixing, and exchanges with the atmosphere (Hamme and Emerson 2006;Nicholson et al. 2012); and (4) the variable integration time of different techniques (particularly in reference to the longer integration time characteristic of the triple oxygen isotopes method) (Juranek and Quay 2005). In the region north of Hawaii, the average GOP rate measured in the surface layer during different periods and using different methods varies about twofold from 0.8 to 1.9 mmol O 2 m −3 d −1 (Table 3). Our estimate of 1.0 mmol m −3 d −1 for the average GOP is well within this range, surprisingly closer to the average of incubation-based estimates (1.1 mmol m −3 d −1 ) than to the average of other non-incubation-based estimates (1.6 mmol m −3 d −1 ) ( Table 3). Our average CR of 1.0 mmol m −3 d −1 is also similar to all previous estimates except from one study that reported considerably larger CR in the 2.4-4.6 mmol m −3 d −1 range based on the nighttime decline of O 2 (Wilson et al. 2014) (Table 3).
The average rates of GOP and CR that we calculated for the region north of Hawaii are sensitive to the statistical approach used to select and average daily estimates. In our study, we excluded estimates calculated in the presence of strong horizontal O 2 gradients, but we included both estimates resulting in negative metabolic rates and estimates from poor fits, here defined as those with a nonsignificant correlation (p > 0.05) between modeled and measured O 2 time-series. Considering that negative metabolic rates have no physical meaning and that poor fits yield unreliable rates, it would be tempting to exclude them from the calculation of the average rates. However, this approach would produce a bias and should be avoided. Specifically, excluding rates from poor fits (based on the p-value of the correlation between modeled and measured O 2 time-series) would lead to an overestimation of the average rates because it would disproportionately exclude rates from days when the amplitude of the diel O 2 change is small. In these cases, the modeled O 2 time-series is not very representative of the measured O 2 time-series because noise is responsible for a larger fraction of the variation. Similarly, the exclusion of negative rates would lead to an overestimation of the average rates of GOP and CR because different sources of noise (instrument precision, uncorrected physical variability, Table 3. Estimates of GOP and CR from the region north of Hawaii based on measurements of O 2 production and/or consumption. Values are average AE standard deviation, average AE (standard error), if the latter value is in parenthesis, or range when two values are separated by an en dash. ML stands for mixed layer. imperfect model assumptions) should symmetrically bias estimates toward larger and smaller values, but excluding negative rates would only remove estimates biased toward smaller values. We quantified the biases due to the exclusion of different data by calculating average GOP and CR rates on subsets of the entire data set (Table 4). Our calculation shows that the exclusion of poor fits and of fits producing negative rates would each approximately result in a 10% overestimation of GOP and CR (Table 4). Conversely, we verified that the criterion used to exclude rates based on horizontal O 2 gradients did not impact the average GOP and CR estimates (Table 4). A separate consideration concerns the way to calculate rate averages and variances starting from daily estimates with different uncertainty. The use of weighted estimators allowed us to decrease the sensitivity of our means and variances to highly uncertain daily estimates. In this study, both average rates and rate variances were lower when calculated with weighted estimators than when calculated with unweighted estimators (data not shown). This result suggests a nonuniform distribution of fit uncertainties that can bias the calculation of GOP and CR, if unweighted estimators are adopted.

Effective temporal resolution and directions to improve the proposed approach
The uncertainty of daily estimates of GOP and CR from glider O 2 measurements was~50% of the average rate while day-to-day variation in primary production was estimated to be 16-20% of the average rate from incubation measurements of bicarbonate incorporation done in summer in the NPSG . Although the variability measured using bicarbonate incorporation is lower than the variability of GOP (Quay et al. 2010), the measurements by Wilson et al. (2015) suggest that our approach could not resolve day-to-day changes in GOP and CR in oligotrophic waters, where the signal-to-noise ratio of diel oscillations is low. In order to understand the time scale of the variability that was resolved by this first application of our method, we compared the timeseries of GOP and CR averaged on a different number of days (Fig. 6A,B,C). As expected, standard errors of average rates decreased with the length of the period used to calculate the averages, from 5 to 9 d (from 0.30 to 0.23 mmol m −3 d −1 for GOP, and from 0.31 to 0.27 mmol m −3 d −1 for CR, on average). Furthermore, when using a time window of 7 or 9 d, maximal rates of GOP and CR coincided with high chlorophyll a concentration and high O 2 saturation in August, indicating a more active plankton community during this period (Fig. 6B,C,D). These results suggest that our method is suited to reconstruct the temporal variability of metabolic rates on a time scale of approximately 1 week, in oligotrophic waters. The same method can likely resolve shorter time scales in Table 4. Weighted mean rates AE weighted standard deviations from different subsets of the daily glider estimates. Daily estimates were excluded based on: (1) average day-to-day absolute O 2 changes >0.75 mmol m −3 d −1 (no horizontal changes); (2) GOP or CR < 0 (no negative rates); (3) fits with a correlation between observed and modeled O 2 with p > 0.05 (only significant fits). The criteria were applied sequentially meaning that each row of the table also excludes the estimates that were excluded in the row above it. All subsets also exclude one outlier estimate with CR = −9 mmol m −3 d −1 . more productive waters, where diel O 2 oscillations have larger amplitude, and daily estimates have lower relative uncertainty. In future applications, we can try to increase the temporal resolution in oligotrophic waters in different ways including changes in data collection, addition of auxiliary measurements, and changes in the fits.
The first factor that hinders the accurate measurement of the O 2 time-series is the finite precision of optical O 2 sensors. It is evident from Fig. 1 that if we only took sporadic and isolated measurements of O 2 we would not be able to reconstruct the diel oscillation that emerges from high frequency observations. On the bright side, each glider O 2 time point used to fit our model was a vertical average of several measurements taken in the surface layer (thus decreasing the uncertainty associated with its value). During future observations, we can further reduce the uncertainty due to instrument precision by increasing the number of samples collected in the surface layer. For glider measurements, this can be achieved by: (1) increasing the sampling rate in shallow waters; (2) decreasing the target depth of glider dives to collect more observations near the surface; and (3) mounting more than one optode on each glider.
O 2 mixing at the base of the surface layer also affected fit quality. We found that vertical O 2 gradients scaled by the depth of the surface layer covaried with the uncertainty of our fits, meaning that we obtained worse fits when the uncorrected diapycnal flux was larger. This phenomenon was more important in summer and fall, when a subsurface O 2 maximum develops below the mixed layer (Riser and Johnson 2008), leading to larger vertical gradients of O 2 . In our case study, the average diapycnal O 2 flux through the base of the surface layer accounted for a volumetric concentration change in the 0.02-0.19 mmol m −3 d −1 range, for diapycnal diffusivities from 10 −5 to 10 −4 m 2 s −1 . This change was smaller than the rates of GOP and CR, but similar to the net O 2 production reported for the mixed layer from the O 2 /Ar saturation ratio not corrected for diapycnal fluxes (Quay et al. 2010;Ferrón et al. 2015). These results indicate that the surface layer can be close to metabolic balance, consistent with some previous studies from subtropical gyres in stratified conditions that suggested that the export of organic material to the mesopelagic zone is not compensated by net photosynthesis within the mixed layer (Knauer et al. 1984;Coale and Bruland 1987;Haskell II et al. 2016). To confirm this hypothesis, as well as to improve our estimate of the biogenic O 2 change, we need more precise estimates of diapycnal oxygen fluxes. Specifically, this could be achieved through direct measurements of K z using microstructure profilers, which have been recently used on underwater autonomous platforms as underwater gliders (St. Laurent and Merrifield 2017) or wave-powered drifting profilers (Lucas et al. 2016).
A third limitation of the glider measurements described in this study is due to the non-Lagrangian nature of our sampling: as the gliders did not follow a water parcel, horizontal O 2 variability contributed to the O 2 change measured in daily time-series. To minimize the influence of horizontal O 2 variations, future glider missions could either follow a surface drifter, or be piloted to follow near-surface currents estimated from other observations and models. A fundamental limitation is that there is shear within and below the mixed layer: the surface layer over which one averages does not move as a column. For this reason, it is unlikely that any sampling design could completely overcome the limitations due to non-Lagrangian measurements.
A fourth phenomenon that could affect rate estimates is the variability in surface irradiance due to clouds. The presence of clouds can affect our estimates both by reducing the amplitude of the diel O 2 oscillation, and by altering the sinusoidal shape of irradiance assumed in our model (in the case of uneven cloudiness throughout the day). While both these effects were expected to reduce the quality of the fits in conditions of low irradiance, we did not observe a significant correlation between fit uncertainty and daily-integrated surface irradiance, as measured from satellite.
As a last point, in our conceptual model of O 2 dynamics in the mixed layer, we neglected processes such as horizontal O 2 exchanges through turbulent mixing, or O 2 photolysis of the organic matter (Kitidis et al. 2014). It is likely that the contribution of these processes to O 2 changes in the surface layer is generally negligible, but we acknowledge that they could be episodically important, particularly in the presence of strong O 2 fronts or anomalously high concentrations of chromophoric dissolved organic matter.

Assumptions on the diel shape of photosynthesis and respiration
The model used in this study works under the simple assumptions that O 2 production is proportional to light intensity and that respiration is constant throughout the day, even though both of these processes could have more complicated diel variations. For example, there is evidence of light saturation in the photosynthetic rate at Station ALOHA (Li et al. 2011), and this process was included in the model by Nicholson et al. (2015). The disadvantage of this approach is the requirement for an additional parameter for light saturation that could be highly variable depending on factors including surface irradiance, the depth of the base of the mixed layer, and phytoplankton community composition. When comparing the results of the fits using a photosynthetic rate linearly proportional to light with those obtained using light saturation with the same parameters as in Nicholson et al. (2015) we obtained similar values for GOP and CR, but worse fits in the case of light saturation (data not shown). We conclude that the simpler model with photosynthesis proportional to light is generally to be preferred.
As for the assumption of constant respiration throughout the day, this is traditionally adopted in diel-based approaches (Odum 1956;Williams and Purdie 1991;Staehr et al. 2010;Ferrón et al. 2015;Nicholson et al. 2015) due to the limited observations and the lack of agreement on the diel variability of respiration in the ocean. In coastal waters, Pringault et al. (2007) estimated that respiration in the light was on average 3.5 times greater than respiration in the dark. Conversely, Grande et al. (1989) observed slightly greater respiration in the dark than in the light using in situ incubations in the open ocean of the NPSG. These conflicting results reflect our limited knowledge of aquatic respiration, leading us to believe that a simple assumption of uniform respiration throughout the day is still the most justifiable one.
As a last point, the detection of residual autocorrelation from several of our fits may indicate that our model assumptions do not accurately reproduce the shape of diel O 2 variations. This might be related to the observation that the average rate of O 2 change is higher in the morning than in the afternoon (Fig. 3B) rather than being symmetrical around noon. We propose that at least two processes can cause this morning-enhanced O 2 increase: (1) higher photosynthetic efficiency in the morning than in the afternoon (Doty and Oguri 1957;Lorenzen 1963) or (2) higher respiration rates during the afternoon due to the progressive accumulation of organic matter from dawn to dusk (Beyers 1963). The proposed processes would need to be experimentally verified before being incorporated in a model, but they are an example of working hypotheses that could lead to improvements in the accuracy of GOP and CR estimates from diel O 2 time-series.

Conclusions
The variability of photosynthesis and respiration in the open ocean is not well characterized due to the considerable efforts required to measure these processes using traditional techniques. The situation could soon improve by using autonomous vehicles to measure diel O 2 cycles that can then be used to calculate rates of photosynthesis and respiration. This approach can be used to estimate rates with a spatial and temporal coverage that would be difficult and expensive to obtain by means of shipboard observations. In the first application of the method proposed in this study, we resolved temporal changes taking place in oligotrophic systems on time scales of approximately 1 week. This temporal resolution can improve with some modifications in the sampling strategy, which have been identified in this study. With such a careful data collection, we could use autonomous observations to study important dynamics such as those taking place at the mesoscale and the submesoscale, or those linked to episodic events such as phytoplankton blooms or the passage of storms.