Uncertainty and sensitivity in optode-based shelf-sea net community production estimates

Coastal seas represent one of the most valuable and vulnerable habitats on Earth. Understanding biological productivity in these dynamic regions is vital to understanding how they may influence and be affected by climate change. A key metric to this end is net community production (NCP), the net effect of autotrophy and heterotrophy; however accurate estimation of NCP has proved to be a difficult task. Presented here is a thorough exploration and sensitivity analysis of an oxygen mass-balance-based NCP estimation technique applied to the Warp Anchorage monitoring station, which is a permanently well-mixed shallow area within the River Thames plume. We have developed an opensource software package for calculating NCP estimates and air–sea gas flux. Our study site is identified as a region of net heterotrophy with strong seasonal variability. The annual cumulative net community oxygen production is calculated as (−5± 2.5) molm a. Short-term daily variability in oxygen is demonstrated to make accurate individual daily estimates challenging. The effects of bubble-induced supersaturation is shown to have a large influence on cumulative annual estimates and is the source of much uncertainty.


Introduction
Marine areas play a fundamental role in the cycling of carbon (Keeling and Shertz, 1992).Photo-autotrophic marine organisms fix CO 2 into organic matter.This organic matter is exported from surface waters by the biological and solubility carbon pumps (Stanley et al., 2010).
Understanding the mechanisms driving these processes is vital for predicting how marine waters will respond to and influence climate change (Guo et al., 2012;Palevsky et al., 2013).Coastal regions in particular have high value to society but are also vulnerable to anthropogenic activities (Jickells, 1998).These regions, which are typically more dynamic than the open ocean and have extensive natural variability, remain a challenge for numerical models (Polton et al., 2013).The accurate detection and prediction of long-term trends, and any response in coastal ecosystems to changing environmental conditions, require the accurate capture of this variability (Blauw et al., 2012).Effective ecosystem-based management of these vital regions requires adequate monitoring, which drives the high demand for good-quality, cost-effective observations of environmental status indicators (Platt and Sathyendranath, 2008).
The balance between dissolved inorganic carbon (DIC) fixation (i.e.autotrophy) and production of DIC through heterotrophy over a specified period is known as net community production (NCP; Williams, 1993).Net autotrophic systems occur when gross primary production is greater than respiration, and net heterotrophic systems occur when respiration is greater than primary production (Ostle et al., 2014).
NCP is a key metric for quantifying the cycling of biological carbon (Stanley et al., 2010).While interpretation of results is challenging and controversial (Williams et al., 2013;Duarte et al., 2013), the direct measurement of CO 2 in the ocean is difficult (Riser and Johnson, 2008).However, as O 2 and C are linked by a stoichiometric ratio (Anderson and Sarmiento, 1994), using in situ measurements of O 2 can offer several advantages over measuring CO 2 directly: dissolved O 2 is chemically neutral, while CO 2 reacts with water to form carbonic acid, which further reacts with other compounds such as carbonates.This buffering makes directly ob-Published by Copernicus Publications on behalf of the European Geosciences Union.

T. Hull et al.: Shelf-sea NCP uncertainty
serving changes in CO 2 difficult.By comparison O 2 can be measured accurately and at high resolution over long periods with relative ease (Wikner et al., 2013).
Estimating net community production rates in the ocean is notoriously difficult (Williams et al., 2013;Duarte et al., 2013).This is in part because the net state is finely balanced between large opposing fluxes and measurements have large uncertainties (Ducklow and Doney, 2013).Approaches have broadly fallen into three categories: in vitro incubation experiments, ocean colour remote-sensing products, and in situ geochemical mass-balance methods.Mouriño-Carballido and Anderson (2009) noted that with in vitro incubation experiments the captured biota may not exhibit the same behaviour as they would in situ.Furthermore bottle samples may be spatially disparate from the source of production.For instance, where deep chlorophyll maxima form, the organisms of interest may not be captured unless specifically targeted (Weston, 2005).Karl et al. (2003) suggested that short, intensive bursts of photosynthesis driven by short-duration changes in light climate are regularly missed with traditional sampling techniques.Kaiser et al. (2005) also concluded that bottle incubations are not suitable to correctly represent the net metabolic balance over larger temporal and spatial scales.
The remote sensing of NCP via ocean colour is in its infancy and requires calibration against reliable in situ measurements (Tilstone et al., 2015;Reuer et al., 2007).These methods are further hampered by insufficient spatial and temporal resolution or obscuring cloud cover (Thomas et al., 2002).Satellites only observe surface waters; they are thus unable to observe the deep chlorophyll maximum, which can contribute up to 60 % of the primary production (Fernand et al., 2013).
Given that production is episodic rather than continuous (Emerson et al., 2008) and the sites of increased production are patchy in nature (Alkire et al., 2012), high temporal resolution in situ sampling is needed (Blauw et al., 2012) Oxygen mass-balance techniques utilise measured changes in oxygen saturation and attempt to quantify the biological contribution to those changes in saturation.The approach to teasing apart the physical and biological drivers to these saturation changes can be subdivided into two groups: those which use a biologically inert analogue to oxygen, typically argon (Kaiser et al., 2005), and those which utilise gas solubility/transfer parametrisations to estimate air-sea exchange.The dual measurement of oxygen and an inert analogue tracer allows determination of solubility changes with fewer uncertainties than using gas solubility parametrisations; however the equipment required for this is not yet in widespread use.
The gas transfer parameterisation approach can be applied to historic data sets; given that the concentration of dissolved oxygen is the most widely measured property of seawater after temperature and salinity (McNeil and D'Asaro, 2014), oxygen-based methods offer many opportunities to reveal new insights into data collected for other purposes.
To date, the majority of oxygen-based NCP estimates have focused on oceanic waters (Alkire et al., 2012).Emerson (2014) noted that coastal NCP values can be 3 times greater than open-ocean values; however, there are too few measurements to be confident in geographical variability.Palevsky et al. (2013) also found during their Gulf of Alaska O 2 /Ar survey that the transitional coastal zone contributed 58 % of the total NCP whilst representing only 20 % of the total area surveyed.The nature of the metabolic balance is particularly important in river-dominated margins, where high carbon and nutrient inputs stimulate primary production and microbial respiration with large seasonal variations (Guo et al., 2012).
The Cefas (Centre for Environment, Fisheries and Aquaculture Science) SmartBuoy network consists of autonomous data collection moorings placed at key locations in the UK shelf seas (Mills et al., 2005;Greenwood et al., 2010).The long-term, high-temporal-resolution multi-parameter data sets produced by the programme provide unique opportunities for observing biogeochemical processes in temperate coastal and shelf seas (Neukermans et al., 2012;Blauw et al., 2012;Foden et al., 2010).
In this paper we present new estimates of NCP from a long-term SmartBuoy mooring situated in the southern North Sea.We explore the uncertainty in these estimates and their sensitivity to uncertain input parameters.Lastly we make our algorithms available as open-source tools for readers to perform their own NCP calculations.

Study site
The SmartBuoy sensor package consists of a Cefas ESM2 data logger coupled with Falmouth Scientific OEM conductivity and temperature sensors (Falmouth Scientific, USA), an Aanderaa 3835 series optode (Aanderaa Data Instruments, Norway), a chlorophyll fluorometer (Seapoint Inc., USA), and a quantum photosynthetically active radiation meter (PAR; LiCor Inc., USA).The ESM2 includes a threeaxis roll and pitch sensor with a internal pressure sensor (PDR1828 -Druck Inc).The data logger was configured to sample for a 10 min burst every half hour.Salinity, temperature, chlorophyll, and PAR are sampled at 1 Hz during the measurement period; oxygen is sampled at 0.2 Hz.
The Warp Anchorage SmartBuoy site, shown in Fig. 1, is located on a shallow bank in the mouth of the River Thames.The site is highly turbid with significant riverine inputs and experiences a 15-day spring-neap cycle with 12 h 25 min semidiurnal tides.Conductivity-temperature-depth (CTD) profiles taken over the last 15 years (Cefas data) have always shown the Warp site to be vertically well mixed.This mix-  ing, together with the shallow water depth, has important implications to the application of oxygen-based NCP methods, which will be discussed later.The main characteristics of the study site are summarised in Table 1.

Data processing
SmartBuoy data undergo rigorous automated and manual quality assurance processes.Automated processes apply a quality flag to data which fall outside realistic value bounds.Manual processes assess the instrument performance and apply flags where the data quality is compromised, e.g.due to biofouling or sensor damage.The CT sensor salinity data are corrected using in situ bottle samples analysed using a Guildline Portsal 8410A (Guildline, Canada) standardised with IAPSO standard seawater.Water depth was calculated using a global tidal model forced with European shelf area constituents (TPX08-atlas).Tidal waves have been shown to arrive almost simultaneously at both the Sheerness and the Warp SmartBuoy site (Blauw et al., 2012); thus model output was validated against the nearby Sheerness tide gauge (UK National Tide Gauge Network) and demonstrated good agreement visually.Windspeed and sea level air pressure were taken from ECMWF MACC (Monitoring Atmospheric Composition and Climate) reanalysis with a 0.125 • grid.ECMWF data were found to compare well with in situ shipborne anemometers used during mooring servicing (see Fig. 2).Details of the ECMWF and tidal model validations and their bearing on the sensitivity analysis are discussed later.
Continuity of the 10-year Warp oxygen data set is hampered primarily by biofouling of the instrumentation.To avoid extrapolation or interpolation of the data, only periods of complete data were used in the analysis.Two contrasting periods were selected, a spring-summer period of 150 days from January to June 2008 and an autumn-winter period of 95 days from September to December of the same year.The 10 min half-hourly burst data from the buoy and the tidal model output were combined with the 6-hourly ECMWF data.These burst means were further smoothed to 25 h averages to remove any structural biases in the data caused by the tidal cycle (Blauw et al., 2012).

Optodes
Aanderaa Instruments model 3830 and 3835 optodes (Aanderaa, Norway) have been fitted to the Cefas SmartBuoys since 2005.Optodes drift due to foil photobleaching in a predictable way (Tengberg et al., 2006), which is well described by a decaying exponential with a decay constant of approximately 2 years (McNeil and D'Asaro, 2014).All optodes used were fitted with the opaque black silicon protective coating.Thus drift is significantly reduced after a burning-in period, and the temperature correction is unaffected (D'Asaro and McNeil, 2013).Sensor drift was corrected with an offset calculated from frequent discrete samples measured with volumetric Winkler titrations (Hansen, 1999).Titrations were performed using an automatic photometric end-point detection system (Metrohm Dosimat 665 Autotitrator); the thiosulfate is intermittently standardised with a standard potassium iodate solution (Wiliams and Jenkinson, 1982).The classical Winkler method if executed with care by a skilled operator offers very low uncertainty (Helm et al., 2009), typically better than 0.2 % (Emerson and Stump, 2010;Ostle et al., 2014).It is however a demanding task that is affected by numerous uncertainty sources, such as contamination of the sample and reagents by atmospheric oxygen and iodine volatilisation.Photometric end-point detection is further affected in highly turbid waters, which can limit the number of successful samples.

Model implementation
NCP is calculated here using a modified version of the zerodimensional oxygen mass-balance (box) model of Emerson (1987) and Emerson et al. (2008).This describes the oxygen mass balance in the mixed layer, assuming no vertical or horizontal advection and no turbulent diffusion across any mixed-layer boundary.
Given that the Warp site is permanently mixed, there is in effect direct connection between the atmosphere and the benthos.It is thus an important distinction from prior studies that our community productivity estimate considers both the pelagic and benthic processes as one system.This method assumes that other oxygen-consuming processes in the water column such as nitrification, methanotrophy, and photooxidation are negligible relative to respiration (Reuer et al., 2007).In our discussion we explore the implications for a site, such as the Warp site, where all of these assumptions may not hold.
The model (Eq. 1) is used to predict the concentration of oxygen at a subsequent point in time given measured physical parameters.Any deviation from the predicted value is assumed to be from biological activity, with a positive value corresponding to net production.This method of NCP estimation makes no distinction between matter which is imported then locally respired and that which is fixed locally.All of these terms introduced below and their estimated uncertainties are summarised in Table 2.
where h is the mixed-layer depth, C is the oxygen concentration in the mixed layer, E is entrainment of oxygen through changes in the mixed-layer depth (Eq.2), G is the gas exchange through diffusive and bubble processes (Eq.3), and J is the net community production.
where C b is the oxygen concentration below the mixed layer.
where k w is the parametrisation of Wanninkhof (2014) (Eq.4).C * is the concentration of oxygen in equilibrium with the one atmosphere as per García and Gordon (1992) using the Benson and Krause (1984) data, B is supersaturation caused by bubble processes (Eq.5), and P slp is sea level pressure, P atm is standard atmospheric pressure (101 325 Pa).
where U is the wind speed at 10 m and Sc O 2 is the dimensionless Schmidt number for oxygen.The typically quoted Schmidt number for CO 2 at 20 • C in salt water (S = 35) is 660.Note the result of Eq. ( 4) is converted from cm h −1 to m s −1 for use in Eq. ( 3).
The square root of the squared mean was used for wind speed to fit with the quadratic k w parametrisation used.Wanninkhof et al. (2009) argue that comprehensive surface forcing models provide little to no improvement over simple wind speed algorithms, and although simple parametrisations cannot capture all the processes that control gas transfer, they appear to capture most.
The injection of bubbles into the mixed layer through wave action can supersaturate the surface waters even if net gas exchange is zero (Liang et al., 2013).Here we utilise a modern k w parametrisation with an explicit bubble equilibrium fractional supersaturation parametrisation B, which enables the influence of the two elements on the NCP estimate to be quantified independently.For B the bubble supersaturation parametrisation of Woolf and Thorpe (1991) is used: where U i is the wind speed at which the equilibrium supersaturation is 1 %.For oxygen Woolf and Thorpe (1991) report this value to be 9 m s −1 .Liang et al. (2013) argue that bubble supersaturation effects at a given temperature differ significantly among parametrisations, and their comparison between Stanley et al. (2009), Woolf andThorpe (1991), and their own parametrisation demonstrates differences in the order of 50 % for argon.The Woolf and Thorpe (1991) parametrisation does not account for any temperature or solubility dependence and is derived from calculated bubbled fields; implementation is however straightforward and the large relative uncertainties in the bubble term will be accounted for in the sensitivity analysis outlined below.
We solve Eq. (1) for NCP (J ) using the analytical solution shown in Eq. ( 6), providing mean values for each variable except oxygen concentration and assuming a constant rate of NCP over the time step, which for this study corresponds to 25 h.The numerical scheme used in this paper was implemented using R, the open-source language and environment for statistical computing (R Foundation for Statistical Computing, www.r-project.org).The analytical solution along with k w and B parametrisations are included in the "airsea" package (Hull and Johnson, 2015).The scheme was validated in silico using numerical estimation; air-sea fluxes were simulated every half second forced with a known value of NCP; the resultant change in oxygen concentration was provided to our model; and the calculated value of NCP compared to the known forced value.This was repeated over a range of input scenarios.
where C 0 is the oxygen concentration at the initial time step (t = 0), and C 1 is the concentration at t.
It should be noted that for this study the entrainment ( dh dt ) term is neglected as the Warp site is a perpetually fully mixed site; as such the entrainment term of Eqs. ( 7) and ( 8) are set to 0.

Sensitivity analysis methods
Accurately assessing the sensitivity of a model output to uncertain input variables has many uses.Primarily it is to determine the precision of the model output and the sources of output uncertainty, knowledge of which informs future research in targeting the main sources of uncertainty if robustness is to be increased (Saltelli et al., 2000).
Local sensitivity analysis methods, such as the so-called one-at-a-time techniques, are limited to providing information only in a very specific location of the parameter space.These methods rely on the selection of an applicable baseline and varying a single input parameter, which ignores the effects of covariant parameter uncertainty (Saltelli et al., 2000).
Global methods such as Latin hypercube sampling with partial rank correlation coefficients (LHS/PRCCs) and the extended Fourier amplitude sensitivity test (eFAST) are capable of assessing multiple locations across the entire parameter space; thus covariant parameter uncertainty is captured.
LHS/PRCC and eFAST have proven to be two of the most efficient and reliable methods in each of their classes, sampling-based and variance decomposition-based respectively (Marino et al., 2008).These two popular methods have differing strengths and weaknesses and measure different properties of the model which together can provide a complete uncertainty analysis.LHS/PRCC is a robust technique for non-linear but monotonic relationships, assuming little to no correlation exists between inputs (Sanchez and Blower, 1997).LHS is an improved method of Monte Carlo which generates more efficient estimates of the desired parameters with far fewer simulation runs.PRCCs are a ranked measure of monotonicity after removing the linear effects of all but one of the variables.A simple one-at-a-time analysis reveals that the variables do indeed demonstrate the monotonic relationships required for effective PRCCs.eFAST provides first-and total-order Sobol indices which indicate the variance of the conditional expectation of the output for a given variable (Saltelli et al., 2000).
LHS is performed by assigning a error probability density function (PDF) to each of the parameters.Each PDF is split into n equiprobable divisions, and each area randomly sampled once without replacement.This table of input variables is then used to calculate NCP, with a new hypercube being generated for each time step.A column-wise, pair-wise algorithm is then used to generate an optimally designed hypercube, where the mean distance between each point and all other points in the hypercube is maximised (Stocki, 2005).We utilise the "improved" LHS implementation within the "lhs" R package (Carnell, 2012) together with the PRCC routine from "epiR" (Nunes et al., 2014).The eFAST scheme is provided by the "sensitivity" package (Pujol et al., 2014).
While there is no a priori exact rule for determining sensible sample size for these methods, minimum values are known to be n = k +1 for LHS/PRCC and n = 65 for eFAST (Saltelli et al., 2000), where k is the number of parameters.
Here we took the usual approach of systematically increasing sample size and checking if the sensitivity index is consistent at least for the main effects, thus demonstrating there is no advantage to increasing sample size as the conclusions remain the same.
LHS/PRCC and eFAST analyses were run 500 times for each 25 h step of the time series, and the results were aggregated.For cumulative calculations k w , B, and C * and the bias element of each measurement parameter were applied globally for the entire time series; that is to say a single hypercube (n = 500) is used to set the bias and scaling factors for multiple runs over the entire time series, while the stochastic uncertainties are applied at each time step independently.

Uncertainty distributions
Critical to the value of any sensitivity or uncertainty analysis is the selection of adequate probability distribution functions for each input parameter (Marino et al., 2008).Table 2 summarises the probability distribution functions used for each of the NCP model input parameters.
The two oxygen terms (C 0 , C) were determined through replicate anchor station Winkler samples taken close to the mooring during maintenance surveys, combined with an estimate of Winkler method error and water bath tests of optode precision.C 0 represents the precision and accuracy of the initial (t = 0) oxygen concentration.We estimate this residual standard error in oxygen determination from the corrected optode, combined with the accuracy of the Winkler samples, to be within ±0.52 mmol m −3 .The error bounds for C, unlike the other measured parameters, are derived solely from the standard error of the difference between the oxygen con-centration at each time time step.This standard error represents both the variability within each 25 h mean and the precision of the optode.
The calculation of k w is conservatively assumed to be accurate to ±15 % (Wanninkhof, 2014).The root-mean-square error (RMSE) from regressions between ECMWF and ship anemometer, shown in Fig. 2, is used to give an estimated wind speed error.For salinity we use the RMSE between the corrected CT, as detailed above, and the bottle samples (0.1).Water bath calibrations have confirmed the SmartBuoy temperature sensors to be accurate to within ±0.1 • C. García and Gordon (1992) provide an uncertainty estimate for the measurement of their oxygen solubility parameterisation of 0.3 %.We have selected a 50 % uniform uncertainty distribution for B, the equilibrium bubble supersaturation term, based on the assessment of parametrisations by Liang et al. (2013).
At the Warp site, given the assertion that it is always fully mixed, the uncertainty in h is reduced to an estimate for the inaccuracies in the tidal model.
Regressions between the predicted height from the model and the Sheerness tide gauge results in a RMSE of approximately 0.4 %.These estimates of parameter measurement uncertainty were combined, using the square root of the sum of squares, with the standard error of each mean observed value.The uniform bias was found to be relatively small compared to the observed standard errors, and thus the overall parameter error is considered to be normally distributed.
Uncertainty distributions for k w , B, and C * were applied by multiplying the parameterised output by a scaling factor sampled from a uncertainty probability distribution.This renders the uncertainty in the parametrisation independent of the input parameters; i.e. k w uncertainty is independent of u uncertainty.

NCP
The 25 h mean chlorophyll time series for the Warp site is shown in Fig. 3a, showing the low levels of chlorophyll in winter, before a marked phytoplankton bloom in late spring.This bloom is known from prior studies to be triggered by improved light climate through increased solar radiation and reduced turbidity (Blauw et al., 2012;Weston et al., 2008).The oxygen saturation anomaly (Fig. 3b), the oxygen concentration minus the solubility (C * ), demonstrates mostly undersaturated near-equilibrium conditions before the bloom, with a large degree of supersaturation during the bloom.Figure 3b illustrates the effects of the B term on increasing the equilibrium saturation concentration and, thus, on reducing the apparent saturation anomaly.Figure 3c shows the ECMWF wind speed data for our study period demonstrating a high degree of variability between days and within our 25 h mean.All NCP values are given as oxygen equivalents unless otherwise stated.NCP is characterised by small, mostly negative fluxes for the first 3 months.This is followed by a marked phytoplankton bloom (Fig. 3a) and resulting positive net community production lasting approximately 3 weeks.Large negative NCP is seen following the bloom, indicating enhanced community respiration.The observed NCP signal is in good agreement with chlorophyll fluorescence (Fig. 3a).
The maximum rate of net community oxygen production was calculated as (485 ± 129) mmol m −2 d −1 with 2σ confidence and precedes maximum observed chlorophyll by 3 days.The mean rate during the non-productive period (January-April) is estimated as (−30 ± 9.5) mmol m −2 d −1 .
The maximum rate of O 2 influx from the atmosphere was (161 ± 47) mmol m −2 d −1 , measured on 1 February 2008, which was concomitant with 14 m s −1 winds (Fig. 3c) and a −2.5 mmol m −3 oxygen anomaly.The maximal rate of oxygen out-gassing was observed on 1 May 2008 of (380 ± 102) mmol m −2 d −1 after the initial peak of the phytoplankton bloom.
Mean gas residence time for oxygen was calculated to be 5 days.Calculating the seasonal net balance (Fig. 4c) at the end of the spring study period (January-June), the cumulative NCP is estimated as (0.5 ± 1.0) mol m −2 at 2(σ ) confidence.The net balance for the winter period (Fig. 5) between 26 September and 30 December is calculated as (−3.We estimate the cumulative NCP for the missing 4-month period of 2010 (July-October) using the mean rate for this period across other years of the 10-year Warp data set, a subset of which is shown in Fig. 6.We calculate the mean value (−18.2 ± 2.3) mmol m −2 d −1 , giving a cumulative estimate for this period of (−2.2 ± 0.4) mol m −2 .There are no significant net autotrophic periods observed between June and September in any other year.
We thus determine that the Warp site is net heterotrophic with an annual oxygen NCP of (−5 ± 2.5) mol m −2 a −1 .However the validity of this assertion is discussed further later.

Sensitivity
Figure 7a shows total-order Sobol indices for the same period computed with eFAST.Here "total" is given to mean the factors' main effects on the NCP estimate, combined with all the interacting terms involving that factor as per Saltelli et al. (2000).The Sobol indices are normalised to the total variance, giving an indication of the fractional contribution to the variance for each factor.Note that, unlike first-order indices, the sum of the total indices can exceed one; in Figs.7a  and 8 we have normalised the total-order indices to one to aid visualisation.
The squared PRCC values from spring 2008 are shown in Fig. 7b.These values are ranked measures, normalised to one, of the degree of monotonicity of each variable on NCP (Sanchez and Blower, 1997).In plainer terms, these are a measure of the independent effect of each input parameter on NCP regardless of whether any input parameter variables correlate.Using squared values makes for easier compari-son with the eFAST indices as the ranked coefficients can be both negative and positive.The relationship between each of the variables and NCP is monotonic for the parameter ranges generated for each time step and thus each PRCC calculation.However, in aggregate over the data set some of the variables can demonstrate a positive and negative (non-monotonic) relationship with NCP.
Both techniques indicate the determination of the change in oxygen concentration ( C) has the largest influence on overall uncertainty, with both the highest PRCC ranking and Sobol total-order indices.The eFAST analysis indicates that C typically accounts for 53 % of the overall uncertainty.Wind speed u is the second-largest contributor, typically comprising 26 % of the uncertainty budget.The bubble supersaturation parametrisation B accounts for 9 %.The gas transfer velocity parametrisation (k w ) and the initial oxygen concentration accuracy (C 0 ) are shown to have similar contributions of 6 %.The García and Gordon (1992) oxygen saturation parametrisation contributes 4 %.Similar results from both sensitivity analyses indicate the model is well characterised by these methods.
The large confidence limits shown for u, k w , and B in Fig. 7 illustrate the large variability in PRCC ranking and Sobol indices over the period studied.This indicates how the relative importance of these factors varies greatly over the data set.The timings for this variability is illustrated in Fig. 8.Here we observe periods (early January and most of March) where C uncertainty is of minimal importance and wind speed uncertainty dominates.The uncertainty in NCP during the onset of the bloom (mid-April to mid-May) is almost completely dictated by uncertainty in C.
LHS/PRCC is not suitable for assessing the effects of measurement and parameterisation bias on the cumulative NCP estimate.Uncertainty in some of the parameters, principally u and k w , do not demonstrate monotonic relationships with the output measure.That is to say uncertainty in u can lead to both increased or decreased cumulative NCP.Thus we present only eFAST indices for cumulative uncertainty in Fig. 9. B is shown to have the largest contribution, accounting for 40 % of the uncertainty in NCP alone, with a further 7 % from interactions primarily with u.

NCP
As the water column at the Warp site is fully mixed, processes occurring at or in the seabed are incorporated into the mixed-layer mass balance and thus the NCP estimate.This includes non-respiration oxygen-consuming processes such as nitrification and the oxidation of reduced compounds other than ammonia and nitrite.A previous study at the Warp site using incubated sediment cores provides estimated rates of sedimentary oxygen uptake of 55 in July and  26 mmol m −2 d −1 in April (Trimmer et al., 2000).Braeckman et al. ( 2014) observed maximal mean rates of nitrification reaching 6 mmol m −2 d −1 and similar for mineralisation in muddy coastal North Sea sediment.This combined with sediment respiration equated to a sediment community oxygen consumption of 15 for February and 20 mmol m −2 d −1 for April.This indicates that a large fraction (perhaps 50 %) of the observed negative NCP at the Warp site could be due to sedimentary processes.
It is important to consider that chemoautotrophic processes, such as nitrification, contribute positively to the metabolic balance but negatively to the oxygen inventory.This is true not just for benthically coupled sites like Warp but for any system where these processes occur.These processes, while assumed small relative to respiration and photoautotrophy by (Reuer et al., 2007) in the Southern Ocean, are likely more important for shelf-sea systems.
There are two events -one at the start of February, another in the second week of March -where high winds appear to coincide with increased negative NCP (Fig. 4a).This could be considered non-intuitive as one may expect increased ventilation to drive the system closer to equilibrium, but this is not the case as shown in Fig. 3b.There are several possible explanations.The optode may be underestimating, or the estimation of saturation concentration incorrect, while in truth the system is supersaturated and is being driven closer to equilibrium during the windy events.We think this unlikely given our error bounds, calibration procedures, and the results from our sensitivity analysis, which indicate the bulk of the contribution to uncertainty is from the u term (Fig. 8).The windy periods could be driving resuspension events which could induce the apparent negative NCP. Lastly, this could be an artefact of the bubble supersaturation term overestimating at high wind speed.The orange line of Fig. 3b shows the effects of the bubble term, and uncertainty, relative to the uncorrected blue line.
While its use in improving our knowledge of carbon cycling is well known, NCP also represents a potential nextgeneration indicator of ecosystem health.The short duration of the bloom and the large impact that a 2-week period has on the annual budget could indicate that annual estimates, while vital for carbon cycling studies, are a less useful indicator for ecosystem health.A carefully resolved bloom period NCP may be more useful.

NCP as carbon equivalents
The commonly used "Redfield" stoichiometric ratio for O : C of 1.45 (Anderson and Sarmiento, 1994;Hedges et al., 2002) was applied to our positive oxygen NCP estimates for easier comparisons with other studies.Literature values for NCP estimates from regions similar to the Warp site are scarce.Tijssen and Eijgenraam (1982) calculated net community oxygen production in the Southern Bight of the North Sea using shipboard 4-hourly Winkler samples.They performed two surveys of 2-3 days in March and April 1980 with 24 h net community oxygen production estimates of 26 and 304 mmol m −2 d −1 respectively.
The rates of net production seen at the Warp site when expressed in units of carbon are of comparable magnitude to other estimates, with a maximal carbon NCP rate of (346 ± 92) mmol m −2 d −1 .Guo et al. (2012) report similar magnitudes of peak NCP from other studies in large river plume regions.Bozec et al. (2006) reported an annual carbon NCP estimate for the entire Thames plume region of 3 mol m −2 a −1 .Their study integrated their four seasonal survey tracks into ICES (International Council for the Exploration of the Sea) regions, of which the Thames plume is one.Our annual carbon NCP estimate of (−3.6 ± 1.8) mol m −2 a −1 represents a much smaller area, measured at considerably higher temporal resolution, for a much longer duration.

Measurement and model uncertainty
Prior oxygen NCP studies have neglected to include the production of oxygen within the time step; that is to say they assume an instantaneous production of NCP at the end of their time step when the measured oxygen concentration and abiotically predicted concentration are compared.This results in the underestimation of the magnitude of NCP.For example, oxygen produced at the start of the time step will out-gas quicker due to the increased air-sea concentration gradient, and when the degree of supersaturation is later measured at the end of the time step the true magnitude of the supersaturation will be masked.
The effect of neglecting the within-time-step NCP is negligible when conditions are near equilibrium saturation.However, during the bloom, neglecting the within-time-step NCP would result in a 45 mmol m −2 d −1 (9 %) underestimation of peak oxygen NCP.
The results from both LHS/PRCC and eFAST techniques support the conclusion that the bulk of the uncertainty in the NCP calculation is dependent on the determination of changing oxygen in the mixed layer.This is in keeping with the observations of the Emerson et al. (2008)  The mean and median value for C standard error were 1.1 and 0.6 mmol m −3 respectively.Greater variability is seen during the bloom, with values up to 7.0 mmol m −3 .During calibration in a thermostatic bath the optodes used typically demonstrated a precision of ±0.3mmol m −3 .This is within the specification from the manufacturer of ±0.4 mmol m −3 and in agreement with the findings of Wikner et al. (2013).Thus it would appear that the largest source of uncertainty constrained here is the large degree of variability captured within the 25 h mean rather than the instrument.The range of values observed within any 25 h period differed by up to 91.2 mmol m −3 during the bloom.During the non-productive period the observations within each 25 h period varied by on average 9.2 mmol m −3 .This variability is shown with the small subsection of the raw oxygen time series presented in Fig. 10.The variability seen here represents both tidal movement of water past the buoy and diel cycling of production.
Thus we believe improvements in identifying homogeneous water masses over the tidal cycle, rather than integrating it entirely, is the best approach to reducing uncertainty with this scheme.
Shipboard transect studies (typically utilising O 2 / Ar methods in open-ocean environments) observe any disequilibrium oxygen in relation to the gas residence time; that is, they assume constant NCP in the period leading up to the measurement (Kaiser and Gist, 2006).It would thus appear that single shipboard transects will struggle to fully capture the tidally induced variability found in areas such as the Warp site.
For the investigation of cumulative uncertainty we consider only the bias in each parameter.The bubble supersaturation term (B), while small in regards to PRCC and eFAST values for an individual estimate (Fig. 7), has a large effect on the cumulative mass balance (Fig. 9).We calculate a pseudocumulative spring period NCP of (2.3 ± 0.9) mol O 2 m −2 resulting from neglecting B: 4 times our true estimate.This relatively large effect is due to the biased nature of the supersaturation term, which serves to only increase the oxygen concentration in the mixed layer.
Optodes tend to drift towards underestimating oxygen concentrations (Wikner et al., 2013), which will typically result in underestimates of NCP.We re-ran our analysis, simulating a 1 mmol m −3 per month negative linear drift, which provides a pseudo-cumulative oxygen NCP estimate for the spring period of (−0.5 ± 0.8) mmol m −2 , which contrasts with our corrected value of (0.5 ± 1.0) mmol m −2 .This reinforces the requirement for well-calibrated, drift-corrected measurements.
Future studies are likely to benefit from newer optode designs than those used here.Together with the improved multipoint calibration equation (Stern-Volmer) of McNeil and D'Asaro (2014), these can offer greater accuracy and pre-cision.The in-air calibration procedures outlined by Bushinsky and Emerson (2013) can reportedly offer frequent in situ calibrations of ±0.1 %.The in-air measurements could also be used to calculate the concentration gradient between the mixed-layer waters and the air, which eliminates the requirement for a C * parametrisation Emerson et al. (2008) noted that at the Hawaii Ocean Time-Series site small daily fluctuations in the measured oxygen concentration caused large fluxes, but these were both positive and negative and had little impact on the cumulative NCP.Fluctuations around zero are seen at the Warp site.These do not tend to cancel out and combine to form a significant negative NCP flux.Emerson (2014) observed the standard deviation of the individual mean annual values is up to ±50 %, which reflects both real inter-annual variability and measurement/model error.This study has produced NCP estimates for the spring period of up to almost 100 % due primarily to the large uncertainty centred around the bloom.Our winter period estimate demonstrates a degree of uncertainty similar to that of Emerson (2014) albeit with a net heterotrophic system.

Advection and sampling uncertainty
Previous studies in open-ocean environments have ignored horizontal advection (Emerson et al., 2008;Nicholson et al., 2008).Air-sea gas exchange is typically considered to be sufficiently rapid that horizontal gradients are too small to drive a significant flux (Alkire et al., 2014).Semi-diurnal tidal systems such as at the Warp site demonstrate horizontal displacement of water masses with a periodicity of 12 h 25 min, with maxima in current speeds every 6 h 12 min, which drive significant horizontal variability (Blauw et al., 2012).
The box model presented here relies on the assumption that the instruments are measuring the same body of water twice; i.e. the comparison of two consecutive 25 h averages represents the same mass of water evolved over time.
If we assume that conditions along the path length are homogeneous on 25 h time scales, in effect the NCP estimates presented here can be thought of as integrating over a length scale proportional to the residual flow.Historic in situ acoustic Doppler current profiler data gathered over 3 months at the Warp site (see Appendix A) show a residual mean current flow estimated at 1.9-2.2cm s −1 , bearing 120 • .This combined with the average tidal excursion of 1.7 km d −1 equates to a observational window of approximately 3.5.km for t = 25 h.
While our 25 h averages and dC error bounds most likely capture the tidal and diel variability, further uncertainty is introduced by submesocale variability such as phytoplankton patches and eddies.Given that Tijssen and Eijgenraam (1982) observed horizontal oxygen gradients of up to 3 mmol m −3 over a few hundred metres, determining to what extent our assumption of homogeneity holds over 25 h and to what extent patchiness within this timescale can influences our estimates is a further step to ensuring a robust NCP estimate.
Residual currents will also affect the NCP estimates by the addition and loss of water from outside of our observational window.Alkire et al. (2014) calculated the advective flux during their glider study and observed daily mean flow of up to 2 cm s −2 .This when combined with their measured horizontal gradient produced the mean removal of (18 ± 10) mmol m −2 d −1 oxygen through horizontal advection.
We have attempted to estimate the oxygen concentration gradient from the tidally driven oxygen variability, that is, the difference between the oxygen concentration at low and high tide.We calculate this for our January period to be approximately 2 mmol m −3 , with low-tide concentration greater than that of high tide.From that we can estimate an advective flux of 51 mmol m −2 d −1 using Eq. ( 9) (Emerson and Stump, 2010).
where v is the Ekman advection velocity.This is not an insignificant flux relative to our calculated winter heterotrophy and would indicate that our site could actually be autotrophic with the heterotrophic processes occurring upstream.It is clear that consideration of advection is required to accurately estimate the annual metabolic state at this site.

Other sources of uncertainty
There are several other known contributors to NCP uncertainty which are outside the scope of this study.Kitidis et al. (2014) argue that all O 2 -based methods underestimate NCP due to photochemical processes, and they report that their modelled photochemical oxygen demand was shown to occasionally exceed respiration, with demand ranging between 3 and 16 mmol m −3 d −1 .Oxygen photolysis was found to correlate with chromophoric dissolved organic matter (CDOM) absorbance at 300 nm.While significant concentrations of CDOM can be found at the Warp site (Foden et al., 2008), the effects are likely mitigated by the typically high turbidity, the associated rapid light attenuation, and shallow (frequently < 6 m) photic depth.Tijssen and Eijgenraam (1982) observed, in the northern end of the Southern Bight of the North Sea in April, vertical oxygen gradients of up to 0.15 mmol m −3 .These can form throughout the day during the phytoplankton bloom.The gradient was reversed during the night, indicating the redistribution of oxygen by vertical mixing over a 24 h period.Takagaki and Komori (2007) found the maximum enhancement to CO 2 gas transfer by rainfall is similar in magnitude to that of high wind speeds.This enhancement is thought mainly to be through increased turbulence and surface area at the air-water interface, and as such it is likely to be most significant where heavy rain is coincident with light winds (Beale et al., 2013).Frew (1997) found that surfactants may be responsible for coastal waters having significantly lower transfer velocities than oligotrophic areas.However Nightingale et al. (2000) found no measurable change in k w during a 30-fold increase in chlorophyll during an algal bloom.We, like Wanninkhof et al. (2009), consider that practically surfactants are always in effect and are thus incorporated into empirically derived k w parametrisations.
Similarly while sea spray may also enhance gas transfer, we believe this to also already be accounted for in the parametrisation.Further uncertainties relating to the parametrisation of k w are likely of little concern without first reducing other, more significant sources.

Conclusions
Our work identifies the Warp SmartBuoy site as an annually net heterotrophic location with strong seasonal variability and autotrophy during the growth phase of the bloom.However, this assertion is brought into question due to significant unconstrained uncertainties from horizontal advection, the determination of which is outside the scope of this study.
We have demonstrated that the largest constrained source of uncertainty in our NCP estimates comes not from the selection of gas exchange parametrisation, or the quality of remote-sensed and modelled parameters, but from the measurement of the changing oxygen concentration.For cumulative annual estimates, the strongly biasing uncertainty of bubble-induced supersaturation is the dominant source of uncertainty.
Constraining the degree of horizontal advection is vital to improving long-term NCP estimates and to determining the overall metabolic balance.Further work should also focus on understanding the nature of the short-term variability associated with changing oxygen concentration to enable better NCP estimates in dynamic areas such as the Warp site.

Figure 1 .
Figure 1.Map of Warp Anchorage study site.

Figure 3 .
Figure 3. Spring 2008 Warp Anchorage time series.(a) Chlorophyll fluorometry.(b) Oxygen saturation anomaly (oxygen concentration minus the solubility).Orange and blue lines represent oxygen saturation anomaly with and without bubble supersaturation effects respectively.(c) ECMWF MACC reanalysis 10 m wind speed.For (b) and (c) thin lines represent 2σ confidence bounds.

Figure
Figure4ashows the calculated NCP for the spring 2008 study period at the Warp site.All NCP values are given as oxygen equivalents unless otherwise stated.NCP is characterised by small, mostly negative fluxes for the first 3 months.This is followed by a marked phytoplankton bloom (Fig.3a) and resulting positive net community production lasting approximately 3 weeks.Large negative NCP is seen following the bloom, indicating enhanced community respiration.The observed NCP signal is in good agreement with chlorophyll fluorescence (Fig.3a).The maximum rate of net community oxygen production was calculated as (485 ± 129) mmol m −2 d −1 with 2σ confidence and precedes maximum observed chlorophyll by 3 days.The mean rate during the non-productive period (January-April) is estimated as (−30 ± 9.5) mmol m −2 d −1 .

Figure 4 .Figure 5 .Figure 6 .
Figure 4. Spring 2008 Warp Anchorage time series.(a) Net community production (J ); negative values correspond to net respiration.(b) Oxygen air-sea gas exchange (G); negative values correspond to movement into the sea.For (a) and (b) thin lines represent 2σ confidence bounds.(c) Cumulative net community production, mean value shown in blue, each run shown in grey, 2σ confidence bounds in red.

Figure 8 .
Figure 8. Warp eFAST total-order Sobol indices over time, indicating changing fractional contributions to uncertainty from each of the main parameters.

Figure 9 .Figure 10 .
Figure 9. Warp eFAST first-order (red) and total-order (Cyan) Sobol indices for cumulative NCP, indicating relative contributions from parameter bias uncertainty to cumulative NCP uncertainty.

Table 1 .
Study site characteristics for w Winter (November-February) and s Summer (June-September), based on multi-year seasonal means.* FTU: formazin turbidity units, ISO 7027.

Table 2 .
Parameters and their uncertainty distributions used for LHS/PRCC and eFAST at the Warp site.