Quantifying pelagic primary production and respiration via an automated in situ incubation system

Pelagic photosynthesis and respiration serve critical roles in controlling the dissolved oxygen (DO) concentration in seawater. The consumption and production via pelagic primary production are of particular importance in the surface ocean and in freshwater ecosystems where photosynthetically active radiation is abundant. However, the dynamic nature and large degree of heterogeneity in these ecosystems pose substantial challenges for providing accurate estimates of marine primary production and metabolic state. The resulting lack of higher‐resolution data in these systems hinders efforts in scaling and including primary production in predictive models. To bridge the gap, we developed and validated a novel automated water incubator that measures in situ rates of photosynthesis and respiration. The automated water incubation system uses commercially available optodes and microcontrollers to record continuous measurements of DO within a closed chamber at desired intervals. With fast response optodes, the incubation system produced measurements of photosynthesis and respiration with an hourly resolution, resolving diel signals in the water column. The high temporal resolution of the time series also enabled the development of Monte Carlo simulation as a new data analysis technique to calculate DO fluxes, with improved performance in noisy time series. Deployment of the incubator was conducted near Ucantena Island, Massachusetts, U.S.A. The data captured diel fluctuations in metabolic fluxes with an hourly resolution, allowed for a more accurate correlation between oxygen cycling and environmental conditions, and provided improved characterization of the pelagic metabolic state.

The increasing demand for pelagic primary production and respiration data has prompted recent methodological developments, aimed at providing a more accurate and costeffective way to quantify metabolic processes (Ducklow and Doney 2013;Collins et al. 2018;Long et al. 2019). Estimates of aerobic respiration and primary production in marine ecosystems are crucial for research and modeling efforts across a wide range of disciplines in marine and environmental sciences (del Giorgio and Williams 2005;Staehr et al. 2012). For example, measurements of respiration and primary production determine the metabolic state of marine ecosystems (i.e., autotrophic or heterotrophic), help identify metabolic processes in the meso-and bathypelagic zone, and validate biogeochemical fluxes in global climate models (Reinthaler et al. 2010;Duarte et al. 2013;Williams et al. 2013). The measurements of these two metabolic parameters can also provide crucial data for management practices such as wastewater treatment and coastal/estuarine restoration (Spanjers et al. 1994;Benway et al. 2016).
Despite the growing demand and importance of metabolic rate measurements, the development of reliable new methods has not kept pace (del Giorgio and Williams 2005). Recent research has reported a large degree of uncertainty in metabolic measurements as traditional methods suffer from various biases such as bottle effects, mismatch of spatiotemporal scales, and the heterogeneity and complexity of the open ocean (Ducklow and Doney 2013;Najjar et al. 2018). Most methods of measuring metabolic rates fall into two categories: (1) in vitro incubation techniques track the rate of change in dissolved oxygen (DO) or dissolved inorganic carbon (DIC) in discrete water samples, while (2) in situ techniques utilize geochemical tracers such as continuous DO time series or isotopic compositions of DO or DIC within water masses. Early measurements of pelagic respiration rates were often generated with in vitro methods. This class of in vitro incubation approach includes two-point incubations utilizing Winkler titrations (Gaarder and Gran 1927), 14 C-spiked incubations (Nielsen 1952), or electron transfer activities (Kenner and Ahmed 1975). In vitro incubation methods are sensitive to short-term perturbations, such as diel variability, vertical plankton migration, upwelling, and weather events. Moreover, the two-point measurements provide flexibility to sampling procedures and experimental designs, especially for shipboard operations. These advantages have led to the continued dominance of in vitro approaches in metabolic process measurements, with few methodological advancements over time.
By contrast, the application of in situ methods has increased considerably due to advancements in instrumentation and sensing technologies. Recent technological advancements in optical sensors and mass spectrometry have enabled wider spatiotemporal coverage (Moore et al. 2009;Goldman et al. 2015). With continuous improvements in optical sensor technology and analysis techniques, low-power optical sensors and instruments are increasingly adapted by underwater autonomous vehicles and mooring platforms. The developments in in situ geochemical tracer techniques and the integration of sensors are essential to long-term autonomous deployments and monitoring. However, the in situ techniques also come with technical difficulties in resolving short-term perturbations inherent in the natural environment.
In situ and in vitro methods have respective advantages and faults, which have been debated for their unexplained discrepancy between the rates estimated by each method. For example, the in situ O 2 /Ar technique tends to yield positive net community production while the in vitro incubation techniques tend to indicate net heterotrophy at the same study site (Williams et al. 2004;Quay et al. 2010). The discrepancy can be organized into two categories of potential bias: (1) the time scale and processes targeted by the given methods and (2) the uncertainty inherent in the method chosen to quantify rates. First, the in vitro incubation approach tends to provide a snapshot of the metabolic state of the ocean (e.g., the instance at which a water sample is taken) while the in situ geochemical tracer approach focuses on steady-state and equilibrium on a longer timescale (e.g., integration of biological activity over days to weeks, depending on the residence time of tracer used). In this comparison, short-term natural perturbations of the water column can drastically alter in vitro incubation results. Comparatively, the in situ geochemical tracer methods tend to smooth out short-term variability by assuming a steadystate of the tracer budget (e.g., O 2 , 13 C). However, the in situ approach requires integrating across the water column, where complex and unresolved physical transport processes such as vertical mixing, lateral advective, and air-sea gas exchange increase uncertainty in the resulting rates (Williams et al. 2013). The depth of integration can also bias, attenuate, or result in a lack of resolution in rate estimates (e.g., mixed layer depth, photic depth) (Palevsky and Doney 2018). Both in vitro and in situ methods also lack sufficient temporal resolution to describe the diel cycling of metabolic processes. The in situ method can only resolve equilibrium on the longer timescale (e.g., weeks), while in vitro incubations often determine metabolic rates with two-point measurements (i.e., beginning and end of incubations, often on the scale of days).
In addition to the first category of bias in time scale and processes targeted, the uncertainty inherent in the method chosen to quantify rates also contributes to overall errors. For example, the in situ O 2 /Ar method via mass spectrometry is estimated to have 0.05% uncertainty for the O 2 /Ar ratio, better than the $ 1% uncertainty in O 2 concentrations produced by Winkler titrations (Hamme and Emerson 2006). In vitro incubations have also been plagued by contamination and disruption during sampling and preparation processes, biasing the results by hindering metabolic processes or changing the microbial community structure in unexpected ways (Suter et al. 2017). The discrete sampling of in vitro methods also suffers from "bottle effects," where unsystematic errors and heterogeneity in the ocean are amplified by sampling and analysis procedures. In vitro incubations also often suffer from unrepresentative incubation conditions where laboratory incubations fail to reproduce the complex variation in the natural environment such as temperature, pressure, and light level fluctuations. Failing to emulate key environmental conditions affects processes dependent on those conditions (i.e., photosynthesis is primarily driven by light, and respiration rates are tightly linked to temperature, pressure, and availability of nutrients).
Recently, technological advances in sensors and engineering controls have led to the development of automated water incubators that can be constructed at lower costs while having higher sampling frequencies. This new class of in situ instrumentation has sought to build upon and integrate the most advantageous components of previous in vitro and in situ methods. Pioneering examples include benthic flux chambers and recent developments of small bottle incubators (Lee et al. 2015;Collins et al. 2018;Long et al. 2019). These incubators combine commercially available sensors and embedded microcontrollers to automate measurements. The increased temporal resolution of data collected from fast response sensors also enables more rigorous statistical analyses that increase metabolic rate precision.
Here, we describe an automated light and dark chamber incubation system for measurements of respiration and primary production with a much higher temporal resolution. Our automated incubator system can resolve hourly metabolic fluxes in a self-contained unit. The automated incubation system minimizes the impacts of biofouling using ultraviolet (UV) LEDs and UV-transparent materials, which ensures a consistent and environmentally relevant starting condition for each incubation experiment. The UV biofouling control also enables longer deployment compared to previously developed incubation systems by eliminating the need for maintenance or human intervention during the deployments. The data analysis in this study uses Monte Carlo simulation techniques to determine high-frequency and robust estimation of metabolic rates. Finally, we provide recommendations for incubation lengths based on the targeted metabolic rates and signal-to-noise ratio (SNR) of the sensor used. Overall, this new class of automated, in situ incubation systems is expected to be highly useful in determining high-resolution pelagic metabolic rates that can be implemented on a range of ocean moorings, vehicles, and floats to provide vital rates for global ocean modeling.

Study site
The incubation system was deployed on the south shore of Uncatena Island, a small island located on the southern end of Cape Cod in Massachusetts, U.S.A. (Fig. 1). The instrument was attached to an anchor $ 70 m offshore from a sandy beach at $ 1.5 m deep (41 31 0 3 00 N, 70 42 0 3 00 W) over four test deployments between August and October 2021. The benthic surface at the deployment site was dominated by seagrass (Zostera marina), underlaid by mixed gravel, and mud sediments.

Automated incubation system
An automated incubation system was developed to conduct hourly incubations and produce oxygen consumption/ production rates. This design improved upon recently developed systems for both open ocean mooring and shallow coastal applications (Collins et al. 2018;Long et al. 2019). The automated incubation system consisted of three incubation chambers, all capable of producing incubation measurements independently (Fig. 1b). The main body of the system was machined out of Delrin plastic, hosting incubation chambers, pumps, and UV LEDs. The top half of the main body was three 1-L cylindrical chambers for seawater incubations. The top side of the incubation chambers was enclosed by clear UV-transparent plastic (Arkema Oroglas UVT), quartz, or black Delrin plates to produce light and dark chambers, respectively. A low-pressure check valve was mounted near the highest level of each chamber for flushing, while another check valve acted as the chamber inlet and was located at the base of the incubation chamber. A sampling port for manual seawater extraction enabled water sampling for additional lab analyses. The bottom of each chamber was enclosed with a UV-transparent quartz plate to enable biofouling control within the chambers (see below). A PyroScience OXROB3 optode was inserted along with a thermistor to monitor DO concentration and temperature in each chamber. The bottom half of the main body was housing for UV LEDs, heatsinks, and pumps. Each chamber had a dedicated EWP-DC30A1230 pump (1.6 L min À1 ) with a 500 μm stainless steel mesh inlet filter. For UV irradiation, an IRTRONIX CA3535 chip was used for each chamber, which hosted a four-by-four array of 275 nm UV-C LEDs with 290 mW of radiant power output per array. The UV LED chips were mounted in the center of their housing at a distance of 5 cm from the chamber to ensure the viewing angle of the LED chip overlapped with the quartz window on the bottom of each incubation chamber-this ensures every surface of the incubation chamber was within the field of view of the LED chip. An external LED chip was mounted in a waterproof housing to irradiate the outside of the acrylic top plate, to prevent external biofouling that could reduce light transmission into the chamber (or a quartz top plate was used). A separate pressure housing contained a microcontroller (Teensy 3.6), electronic control units, PyroScience DO measurement modules, analog-digital converters (ADC), and batteries ( Fig. 2). All UV LEDs and pumps were controlled by latching relays. Each LED chip was powered by a 500 mA constant current driver (Recom Power RCD24-0.50). Each PyroScience module was connected to a 16-bit Adafruit ADC board (ADS1115) for differential ADC. A custom spectral photosynthetically active radiation (PAR) sensor (SparkFun AS7265X) was deployed concurrently to provide the context in ambient light. The raw spectral measurements were integrated between 400 and 700 nm to calculate PAR.
In this study, the incubation system was programmed to conduct 1-h incubations. At the beginning of each hour, the sampling pump was triggered for a 5-min flushing period, flushing each chamber approximately seven times (e.g., 0:00 to 0:05 min). After flushing, individual incubation chambers were sealed by the check valves for 50 min (e.g., 0:05 to 0:55 min), during which the integral optode and thermistor monitored DO concentrations continuously. At the end of the incubation period (e.g., minute 0:55), the UV LEDs were activated to irradiate the chamber for 5 min (e.g., 0:55 to 0:00) to minimize biofouling.

UV biofouling control
Biofouling occurs on all oceanographic instrumentation, often affecting the precision and accuracy of measurements that are exacerbated by extended deployment times (Ward 2021). For measurements of pelagic productivity, the incubation system requires biofouling control to ensure the measurements do not exhibit a baseline shift in magnitude due to enhanced growth on surfaces or reductions in natural light supplied to each chamber. To assess the effectiveness of UV LED as biofouling control, experiments were designed based on a recent study reporting that biofilm growth on plastic films in the ocean was strongly linked to reductions in light transmittance (Nelson et al. 2021). The incubation system was submerged in a freshwater tank, supplemented with a natural microbial community. The incubation system was programmed to function as a normal field deployment. The UV LEDs irradiated each chamber hourly for different periods of time before the pumps flushed and renewed water samples within each incubation chamber. Triplicate low-density polyethylene (LDPE) films were placed in chambers programmed with different UV dosages for 27 d. Biofouling on the LDPE films was subsequently measured using the light transmission port of a UV-visible spectrophotometer (PerkinElmer Lambda 650s; Nelson et al. 2021). The UV dosages were varied by controlling irradiation time: 0, 5, and 10 min/h-long incubation. The UV dosage is calculated from the UV irradiation level measured by a spectral radiometer (StellarNet Solar-Rad). In addition to UV treatments, LDPE samples were placed in a saturated mercury chloride solution as an abiotic control.

Metabolic rate calculation with Monte Carlo simulation
The automated incubation system offered improved temporal resolutions within each 50-min incubation period with $ 1 Hz sampling frequency compared to existing methods (e.g., two-point linear regression). The data density also provided more statistical context for understanding metabolic processes. For data analysis, the first 5 min of each incubation period was excluded to account for sensor response time and mixing within the incubation volume. Monte Carlo simulation was utilized to extract the predominant oxygen consumption/production rates from the DO concentration data. For rate extraction via the Monte Carlo simulation, a sample set of two DO concentrations was randomly sampled within the defined incubation time series, and the slope (R s ) between the given two points was calculated by Eq. 1: where R s is the calculated slope for the given sample set, [O 2 ] is the DO concentration and t represents the elapsed timestamp for each point in the random sample set. The f subscript marks the endpoint while the i subscript marks the beginning point (earlier time) of the sample set. The technique then repeated random sampling and slope calculations over 2000 designated repetitions. Histograms were generated from all 2000 R s values for each incubation period. All R s were then binned by predefined binning resolutions (i.e., 500 bins for 2000 repetitions of R s in this study). The mean value of peak bins with a combined weight of over 80% sample size was then calculated as the representative oxygen production/ consumption rate (R c ) of the given incubation period by Eq. 2: where E high and E low are the higher and lower boundary of R set selected for R c calculation. In addition to calculating R c , this study also investigated the optimal incubation length with the Monte Carlo simulation technique. The goal was to quantify the shortest possible incubation length given the expected DO flux and SNR of the sensor used. The shortest possible incubation length was calculated by applying Monte Carlo simulation on an artificially generated dataset (A) with a given rate of change (a) and SNR (i.e., a straight line with a known slope and randomly generated noise according to a given SNR from the sensor used). The Monte Carlo simulation was applied on increasing portions of the dataset until the simulation results correctly predicted the rate of change (i.e., the model repeated the simulation with increasing lengths of the dataset until the 5% error criteria were achieved). The shortest incubation length required (t optimal ) was then calculated based on the sampling frequency of the sensor used (Eq. 3). Assessment

Efficiency of UV biofouling control
The 27-d testing of different UV treatments (0, 5, or 10 min/h) yielded substantial differences in biofouling, as inferred through the decrease of light transmission through the plastic films. The total UV dose (kJ) at the top of the incubator was calculated by multiplying the radiant power output by exposure time; the 5-min UV irradiation ($ 5 W m À1 ) was equivalent to $ 15 kJ, and the 10-min UV irradiation was equivalent to 30 kJ. The extent of biofouling in each treatment was quantified at two major chlorophyll absorption peaks, 430-450 and 665-685 nm (Fig. 3a).
In the 430-450 nm region, the integrated transmittance exhibited significant differences between the UV-treated samples and the non-UV-treated samples (Fig. 3b). On the third day of this time series, the transmittance of non-UV-treated samples decreased to $ 84% while the transmittance of UV-treated samples stayed at $ 88% for the first week of the trial. A general trend of decreasing transmittance was observed in the second half of the 27-d trial. On day 4, control, UV treated, and nontreated showed significant differences (p < 0.05), while the UVtreated groups (5 and 10 min irradiation) did not show significant differences via ANOVA test. The UV-treated groups showed efficacy for controlling biofouling until after day 24. On day 27, the 5-min UV irradiation treatment group was not significantly different from the nontreated group.

Stability of Monte Carlo simulation technique
The Monte Carlo simulation technique relied on randomized sampling to generate histograms of the estimated rates for each incubation-each simulation may generate slightly different estimations as a result of randomized sampling (Fig. 4). It is important to ensure the stability of the simulation outcome by using appropriate simulation parameters. Simulation parameter testing was conducted on a set of artificially generated samples with known statistical properties such as linear slope (R s = 1) and noise-to-signal ratio (AE 20% noise). The Monte Carlo technique was applied with various parameters (e.g., number of repetitions, the width of bins, mean value averaging criteria). The mean rate simulation results indicated the general simulation outcome reached stability as soon as the number of repetitions exceeded 2000, stability was indicated by the mean rate prediction (Fig. 5a) and error of simulation outcome. The error of the simulation outcome was suppressed by increasing the number of repetitions, declining to lower than 2% uncertainty once the number of repetitions exceeded 4000 (Fig. 5b). This test conceptually supported the use of the Monte Carlo simulation technique in this application, since the noise magnitude from the optodes was much lower than the artificial dataset, and each incubation period provides an ample number of data points ($ 3000) for randomized sampling.

Field deployments and data products
The automated incubator was deployed for approximately 72 h four times throughout the summer of 2021. The incubator was programmed to generate DO measurements at $ 1 Hz; each incubation period consisted of roughly 3000 measurements of DO concentration. The Monte Carlo simulation was then applied to the time series of each incubation period respectively, generating a DO consumption/production rate (μmol L À1 h À1 ) for each hourly incubation. Field data shows hourly DO production rates from the light chamber (with quartz top) and the DO consumption rate from the dark chamber (Fig. 6).
With hourly O 2 fluxes calculated, more in-depth analyses were conducted to resolve the magnitude of photosynthesis and respiration. In respiration analysis, the goals were to determine (1) the mean respiration rate and (2) the required incubation length to achieve stable DO flux estimates. Instead of applying Monte Carlo simulation to the full length of DO time series in experiments, this analysis estimated rates with only part of the incubation length. Here, the Monte Carlo simulation was applied to the first 5 min of the time series in each dark incubation experiment (n = 73), and then 10 min, and so on until the full length of the incubation experiment is reached (Fig. 7a). A mean hourly respiration rate was calculated to be À1.62 AE 1.17 μmol L À1 h À1 O 2 at the deployment site.
For photosynthesis rate estimation, a trend between PAR and the oxygen production rate in the light chamber was observed in Fig. 6. The oxygen production rate showed sensitivity to short-term fluctuation in hourly averaged PAR (Fig. 7b). The relationship between PAR and photosynthetic DO fluxes can be characterized by comparing all in situ hourly DO fluxes from the light chamber and the hourly integrated PAR measurements. In this analysis, DO flux from gross primary production (GPP) during the daytime (i.e., when underwater PAR measurements are nonzero) was calculated by adding the previously calculated mean respiration to net primary production (Fig. 7b). A logarithmic curve was fitted to indicate the general correlation between PAR and photosynthesis. The fitted curve indicated a light compensation point of 82.2 μmol s À1 m À2 and decreased photosynthetic efficiency as PAR increased. However, light saturation was not apparent and GPP continued to increase across all measured light levels.
Furthermore, the DO flux rate could be compiled into an averaged diel DO flux estimate (Fig. 8). The diel DO flux roughly follows the trend of averaged PAR and suggested that hourly photosynthetic fluxes were higher in the morning when compared to those in the afternoon with similar PAR levels (also see Fig. 7b). Notably, the DO fluxes during dawn and dusk hours had higher uncertainties, reflecting the variability of DO flux directions during crepuscular periods. Finally, the daily net ecosystem metabolism (NEM) was calculated directly from the summation of the DO fluxes over 24 h (Eq. 4) and then the GPP was calculated by Eq. 5:  low bars) and dark chambers (blue bars). The red line represents hourly integrated PAR measurements from the spectral PAR sensor. Data suggests that the automated incubation system is sensitive to short-term environmental conditions (e.g., diel cycles, PAR fluctuations).

NEM ¼
The NEM at the deployment site was 2.56 AE 10.02 μmol O 2 L À1 d À1 , and the GPP was 49.59 AE 15.26 μmol O 2 L À1 d À1 . These estimates indicated that the water column was near the metabolic balance, driven by the high photosynthesis during the day. We observed no apparent seasonal trend across four deployments. Mean respiration rates and photosynthesis-PAR relationship were comparatively stable between late August and early October as the daily temperature range in the water was stable roughly between 24 C and 29 C.

Optimal incubation length
The optimal incubation length can be estimated by Monte Carlo simulation using artificially generated time series (sampling frequency of 1 Hz), the sensor SNR, and the expected rate of change (e.g., O 2 flux) (Fig. 9). The SNR and incubation length exhibited a relationship of exponential decay while the incubation length and expected O 2 flux share a polynomial correlation. The relationships between all three variables can be summarized in Eq. 6: t optimal ¼ 16:3 þ 3:46 Â 10 À4 SNR 2 Â R c þ 0:0713 Â SNR Â R c À 0:332R c þ 47:51e À16:5Rc Â SNR À 1:36SNR: The t optimal on the left side of the equation represents the shortest incubation time required for accurate rate prediction. With modern optodes, the SNR is generally low enough to   The diel signal was a result of the averaging of four testing deployments, each deployment spans $ 72 h. The error bar indicates the standard deviations of the respective hourly DO fluxes. The red line is the hourly average PAR over all test deployment dates. allow a very short incubation period according to Eq. 6. The incubation length can be as short as 15 min, as shown by the "flat plane" area in Fig. 9.

Discussion
This study demonstrates the ability to quantify pelagic DO fluxes and their hourly variation in a shallow coastal ecosystem using an automated incubator. The DO optodes provide a high temporal resolution time series, enabling high-frequency (hourly) measurements of pelagic primary production and respiration. The high data density of the DO concentration time series allows the use of the Monte Carlo simulation technique as opposed to the traditional calculation (e.g., two endpoints or linear regression). The Monte Carlo simulation technique can resolve the dominant DO rates without the error apparent in the linear regression models, especially when the time series does not resemble a straight line shape. The hourly fluxes were useful for modeling the correlations between GPP and available PAR, which was not possible with previous, lowerresolution techniques. The compilation of hourly flux rates also enabled the estimation of mean pelagic rates at the deployment site. Resolving these fluxes with higher temporal resolutions is essential for correctly assessing parameters of ecosystem metabolism such as NEM and GPP and determining the environmental drivers of these processes.

UV treatment
UV treatment is the best candidate for a noninvasive, nonchemical approach for biofouling control. UV disinfection has been used for wastewater treatment and drinking water disinfection since the 1950s (Kruithof et al. 1992). Some UV wavelength is absorbed by nucleic acids in microorganisms, hindering the reproductive ability of targeted microorganisms to achieve inactivation (Gates 1930). To mitigate the potential bias from biofouling in our deployments, the most effective UV dosage of 30 kJ m À2 was chosen. In this study, UVtransparent plastic or quartz plates were installed over the light chamber to ensure sufficient UV dosage for the outside surface. Biofouling control outside of the chamber is especially important for light incubations. We observed no visible biofouling in or outside of the incubation chambers for both the UV-transparent plastic-covered chamber (with supplemental, external UV light) and the quartz-covered chamber. The UVtransparent plate needed external UV light because it only transmitted $ 50% of UV irradiation while the quartz plate transmitted $ 95% of UV irradiation. The UV dosage of our system exceeds the EPA recommendations for UV disinfection by $ 100 times for eliminating 99% of pathogens in drinking water (Schmelling 2006), suggesting that a higher UV dosage may be required to maintain long-term biofouling control in marine applications. Such a high UV dosage relies on high power input, posing significant challenges for in situ instrumentation. Recent studies have shown that given sufficient UV dosage, pulsing UV irradiation is more effective than continuous UV irradiation (Olsen et al. 2016;Zou et al. 2019). Pulsing UV irradiation also provides the technical advantage of driving UV LEDs with a higher current to provide higher output, since LEDs can endure short periods of higher forward current. With pulsing UV irradiation, the irradiation duty cycles, and pulsing frequencies can be adjusted to minimize power consumption while maintaining the high efficacy of biofouling control in future applications.

Nonlinear behavior in DO time series
The Monte Carlo simulation technique provided an important opportunity to improve time series analysis workflow. Segments of short DO time series data (e.g., 1 h) often exhibited complicated behavior rather than an ideal straight line. The DO fluctuation represented not only the effects of photosynthesis and respiration but also the potential biases that contribute to the unsystematic behavior in the DO time series (Ducklow and Doney 2013;Collins et al. 2018). Monte Carlo simulation techniques were introduced to investigate the behavior of such complex systems (Harrison et al 2010) and were employed here to reduce errors associated with simpler regression techniques. For example, a sharp decrease in DO concentration is occasionally observed from the dark chambers during daytime (Fig. 10b). This is consistent with the observed "deep breath" phenomenon in a range of respiration research (Robert 2012;Guillemette et al. 2013). In the deep breath phenomenon, it is proposed that highly labile dissolved organic carbon (DOC) is rapidly respired during the initial phase of the incubation. The rate of respiration subsequently slows as the community shifts from the highly labile to semi-labile fractions of DOC. In other instances, the DO concentration showed net production for a portion of the incubation, then reverted to DO consumption after the initial production phase in the dark chamber. This pattern resulted in sporadic positive DO flux estimates in the dark chamber during the daytime. In this example, the DO flux estimate can be potentially biased by taking the sunlit phytoplankton community and instantly placing them in the dark. The abrupt light changes are not representative of natural conditions. In addition, this observation can also potentially be attributed to microbubble accumulation inside the chamber or on the optodes (Wikner et al. 2013;Collins et al. 2018).
While these examples could be ascribed to some other potential processes such as photorespiration, tidal influences, variability in temperature, variability in ambient light, mixing within the chamber, and inherent bottle effects, it is difficult to untangle the probable cause via hourly incubations (Bender et al. 1999;Jacquet et al. 2001;Robinson et al. 2002). Regardless, the Monte Carlo simulation, as a probability simulation technique, does not depend on the shape of the entire time series to make predictions. As shown in Fig. 10, Monte Carlo simulation tends to generate rate estimates for the "predominant" trend in time series data-enabling quantification of the driving processes in time series data. The Monte Carlo simulation technique is also important for developing a universal data analysis workflow since it requires minimal quality control for the raw data. This workflow is also essential in developing a protocol for autonomous incubation and operation (see Comments and recommendations).
It is worth mentioning that the probability distribution of DO consumption/production rate can change depending on the random sampling function used (Harrison et al 2010), and the number of repetitions set for simulation. In this study, a uniform probability distribution function is used to give equal probability in random sampling throughout the given incubation time series. The number of repetitions controls the density of the outcome histogram. The Monte Carlo simulation tends to produce a "sharper" peak in the probability distribution spectrum with a significant proportion of calculated rates residing within the peak bin, while an insufficient number of repetitions often leads to a "flatter" probability distribution spectrum. An insufficient number of repetitions produces results more susceptible to the stochastic nature of random sampling. The convergence of simulation results with an increasing number of repetitions also suggests the linearity of respiration and photosynthesis rates under stable environmental conditions on an hourly timescale (del Giorgio and Williams 2005).

Comments and recommendations
Automated sampling, high-frequency DO measurements, and new data analysis approaches have enabled this automated incubation system to sample and produce hourly DO fluxes. Resolving hourly fluxes is a significant improvement in temporal resolution compared to traditional in situ and in vitro methods. The recurring hourly sampling strategy also alleviates potential uncertainties by providing many replicates through time. Data with high temporal resolution is also essential for resolving DO fluxes in dynamic environments such as coastal and estuarine ecosystems. In addition to applications in coastal and estuarine systems, the automated incubation system can also provide insight into subdaily variation in the euphotic surface ocean, where metabolic rates remain largely unresolved due to the heterogeneity and complexity of biophysical interactions (Ducklow and Doney 2013). The hourly fluxes can also be integrated over time to produce the mean pelagic metabolic state. Integrating fluxes over time also provides better estimates of uncertainty/variability associated with metabolic fluxes in the open ocean (i.e., uncertainty in traditional methods is usually associated with analytical errors instead of variability in the natural environment). Automated hourly measurements can also provide more information for other acute environmental variations that influence metabolic activity in the water column as evidence suggests that metabolic rates change significantly on daily and hourly timescales (Caffrey 2004;Staehr et al. 2012). There are numerous logistical advantages of automated incubation, including the acquisition of more data without the preparation of shipboard or laboratory incubations, and sampling in more inaccessible environments where acquiring and physically tending to samples is difficult.
The UV LEDs mitigate the previous concern of biofouling and the requirements of constant maintenance (Collins et al. 2018;Long et al. 2019). Biofouling control consumes the most power in this incubation system as an essential function. However, pulsing UV LEDs at shorter duty cycles and higher frequencies presents a promising solution to extend the deployment span and versatility of this instrument (Olsen et al. 2016;Zou et al. 2019). Future research on efficacy for integrated biofouling control with pulsing UV irradiation is needed to optimize instrument design and programming. A significantly reduced power demand also allows miniaturization in future designs while keeping the inherent bottle effects from small incubation volumes in mind. Future investigation into optimized surfaceto-volume ratios for small-volume incubations is needed to guide the miniaturization of automated incubators. In an optimized design, a smaller incubation volume implies more efficiency in many engineering aspects (e.g., reducing pump/flush time, reducing UV irradiation output, reducing battery size). In addition to the smaller incubation volume, the shape and material of the chamber should be taken into consideration for optimal light exposure. The current incubator design had chambers with wide diameters and short heights to incubate water samples under downwelling light conditions. However, future designs of similar incubators should consider constructing incubation volumes with transparent material to enable capturing diffuse light. The miniaturization of a conceptualized incubation system can unlock the potential for the data collection on various platforms (e.g., floats, AUVs, moorings).
To achieve wider applicability, future designs should also take the length of incubation into account. For example, while the physical design of incubation volumes tends to be fixed based on the surface-to-volume ratio for biofouling and bottle effect control, incubation length can be programmed adaptively based on the initial observed rates to target the ideal incubation length to resolve fluxes. As a starting point for future studies, the Monte Carlo simulation technique provides recommendations for incubation length based on targeted metabolic rates and sensor SNR, which can significantly shorten incubation length compared to previous in situ incubation studies. An adaptive sampling strategy could also be programmed as part of the algorithm in the microcontroller, adjusting incubation length during deployments based on the metabolic rates calculated from the raw DO time series. The adaptive sampling strategy will allow the optimization of the incubation experiment autonomously and in real time, approaching the highest possible temporal resolution while maintaining data integrity.
In addition to improving the temporal resolution, this incubation technique also opens many opportunities for future developments. Specifically, improved spatial resolution can be achieved by deploying multiple units throughout the water column or in an array of horizontally spaced systems; the incoming light can be conditioned by using different chamber materials to investigate different processes (e.g., photosynthesis, respiration, photo-oxidation); and fluid sampling ports can be integrated for future collaboration with other in situ sensors and lab-on-a-chip technology (Wang et al. 2015;Beaton et al. 2017;Vaughan et al. 2018). The ability to accommodate other in situ sensors opens a wide range of applicability to incorporate many biogeochemical parameters (e.g., dissolved inorganic carbon, nitrate, phosphorus).
With the technological advantages discussed above, this automated incubation technique provides many possibilities for advanced conceptual designs. A design of such conceptualized incubation systems should incorporate interdisciplinary efforts to assess the uncertainty of metabolic measurements, statistical approaches for data analysis, and engineering optimization. Consideration of design should include types of sensors used, materials used to construct the incubation chamber, shape and size of the incubation chamber, and wavelengths of UV irradiation used for efficient biofouling controloptimization of these elements would be essential to the success of future prototype development. Moreover, an autonomous incubation system based on an adaptive sampling strategy and the Monte Carlo simulation technique would greatly improve the understanding of pelagic metabolism and the longstanding discrepancies between methods to specifically address the lack of spatiotemporal resolution in pelagic metabolic rates.

Data availability statement
All study data and a MATLAB script for the Monte Carlo simulation are archived at the Woods Hole Open-Access Server (https://hdl.handle.net/1912/29487).