Evaluation of new and net community production estimates by multiple ship-based and autonomous observations in the Northeast Pacific Ocean

The net flux of carbon dioxide (CO₂) between the ocean and the atmosphere plays an important role in regulating the concentration of this greenhouse gas in the atmosphere and influencing the Earth’s climate (Raven and Falkowski, 1999; Sarmiento and Gruber, 2006). The air–sea transfer of carbon is dominated by two mechanisms: the solubility pump, which draws cold, CO₂-rich waters to depth, and the biological pump, which transports organic carbon from the surface ocean to depth (Boyd et al., 2019). As part of the biological carbon pump (BCP), phytoplankton convert dissolved CO₂ in the euphotic zone into organic forms through photosynthetic carbon fixation. A fraction of this organic carbon is exported to depth through various pathways, including gravitational sinking of particles and active transport mechanisms such as subduction, mixed layer detrainment, and migration of organisms to and from depth (Ducklow et al., 2001; Boyd et al., 2019; Lacour et al., 2019). While the biological pump is central to our global climate, the mechanisms controlling its strength and efficiency remain poorly quantified.

Organic carbon export to depth is a function of the net community production (NCP) or new production (NP) in the surface ocean. NCP measures the balance between gross primary production (GPP) and community respiration (CR) in an ecosystem, representing the net production of organic matter that can be exported from the surface ocean to depth (Emerson, 1987; Li and Cassar, 2017). NP measures production that is supported by external inputs of biologically available nitrogen, primarily from below the mixed layer (ML; Eppley and Peterson, 1979). When the carbon content within the ML is at steady state, NCP and NP measurements are equivalent and represent carbon export potential (Li and Cassar, 2017). Though variable, stoichiometric relationships presented in Equation 1 allow us to estimate NCP and NP by evaluating the mixed layer budgets of oxygen, nitrogen, phosphorus, and/or carbon:

where a, b, m, n, x, y and z are coefficients that can vary based on stoichiometry. NCP and NP can be tracked either by the consumption of inorganic carbon and/or macronutrients (i.e., PO₄³⁻ and NO₃⁻) or by the generation of O₂ or organic carbon. Various constituents in Equation 1 can be measured in situ via shipboard or autonomous measurements, or in vitro using incubation experiments. Net biological oxygen production or consumption can be estimated by comparing ratios of O₂ and argon (Ar) in the seawater to those in the air. The two gases are similarly influenced by physical solubility processes, but the O₂ concentration is affected by photosynthesis and respiration, while Ar is biologically inert (Craig and Hayward, 1987). Additionally, the biological contribution to a measured change in a biogeochemical tracer (DIC, NO₃⁻ and O₂) can be estimated by constraining the abiotic, physical processes (Bushinsky and Emerson, 2015; Fassbender et al., 2016; Haskell et al., 2020). Other commonly used approaches to measuring NCP include the creation of a carbon budget in the ML using measurements of particulate organic carbon (POC) and dissolved organic carbon (DOC; Bates et al., 1998; Sweeney et al., 2000).

Historically, NCP has been measured by shipboard observations from incubations (Bender et al., 1987; Duarte et al., 2013; Emerson, 2014; Quay, 2021) or discrete sampling (Reuer et al., 2007; Palevsky and Quay, 2017). However, recent technological advancements in shipboard underway measurements (Cassar et al., 2009; Izett and Tortell, 2020), autonomous platforms (Riser and Johnson, 2008; Emerson and Stump, 2010; Alkire et al., 2012), and remote sensing (Tilstone et al., 2015; Li and Cassar, 2016) have expanded the availability of NCP data on regional and global scales. While previous studies have relied on discrete O₂/Ar sampling (Emerson et al., 1995), the development of ship-board Membrane Inlet Mass Spectrometry (MIMS; Kaiser et al., 2005; Tortell, 2005) and Equilibrator Inlet Mass Spectrometry (EIMS; Cassar et al., 2009) have allowed for the collection of continuous oceanic O₂/Ar datasets. These underway methods offer high resolution data but are confined to the spatial coverage and length of research cruises (Seguro et al., 2019). Advancements in underwater biogeochemical sensors have also allowed for increased oceanic measurements from autonomous platforms, underway measurements, and shipboard measurements. Autonomous platforms, including Lagrangian floats (D’Asaro, 2003), gliders, wirewalkers, and profiling floats outfitted with biogeochemical (BGC) sensors (e.g., BGC-float), including O₂, NO₃⁻, backscatter, and pH, enable NCP measurements over large spatial areas (Alkire et al., 2012; Plant et al., 2016; Yang et al., 2017).

In general, in situ methods allow for higher resolution data and longer-term datasets than in vitro methods; however, these methods are highly dependent on sensor accuracy and the degree to which abiotic processes in the upper ocean can be parameterized (i.e., air–sea flux, vertical mixing, advection; Huang et al., 2022). Additionally, there has been significant debate regarding the comparison between incubation-based NCP measurements and in situ biogeochemical measurements (Marra, 2009; Duarte et al., 2013; Ducklow and Doney, 2013). While environmental parameters can be controlled using in vitro methods, incubations can lead to biases caused by factors including artificial light regimes (Godoy et al., 2012), underestimates of grazing by zooplankton (Robinson and Williams, 2005), and incorrect assumptions about respiration rates throughout a 24-hour period (Bender et al., 1987). Each method of capturing NCP and NP has a unique set of benefits and limitations (Figure 1), and only a few studies have compared simultaneous measurements (Alkire et al., 2012; Manning et al., 2017; Timmerman and Hamme, 2021).

In this study, we evaluated multiple methods of estimating NCP and NP in the ML using measurements taken in the vicinity of Ocean Station Papa (Station P; 50°N, 145°W) from August 14 to September 7, 2018, as part of the NASA project EXport Processes in the Ocean from RemoTe Sensing (EXPORTS; Siegel et al., 2016; Siegel et al., 2021). Four in situ, two in vitro, and two satellite-based methods, deployed from multiple platforms, resulting in thirteen independent estimates of NCP and NP, were evaluated (Table 1). While there is no benchmark measurement for NP and NCP, we leveraged these concurrent measurements made by different instruments, platforms, and tracers to assess the factors driving differences in measured NCP and NP, thus constraining the carbon export potential of the system. To our knowledge, this comparison of NCP and NP approaches is the most comprehensive to date.

Measurements were collected as part of the EXPORTS 2018 field campaign (Siegel et al., 2021). Underway measurements were conducted aboard the R/V Roger Revelle and the R/V Sally Ride between 49.5–51.0°N and 144–146°W from August 14 to September 8, 2018. The R/V Revelle (hereafter referred to as the “process ship”) followed an instrumented Lagrangian float (LF; D’Asaro, 2003) deployed in an isopycnal around 100 m on August 14, 2018. The R/V Ride (hereafter “survey ship”) ran small- and large-scale surveys around the process ship to provide spatial context. Data from the EXPORTS campaign seaglider and two Slocum Gliders at the Ocean Observatories Initiative (OOI) Global Station Papa Array were utilized to maximize the spatiotemporal coverage of observations. The three gliders had slightly different morphologies but were functionally the same. OOI Glider 276 (G276) and the seaglider were deployed on July 23 and 27, respectively, during an NSF OOI cruise on the R/V Sally Ride near Station P. OOI Glider 469 (G469) was deployed on August 21 from the subsequent EXPORTS cruise. A BGC-float (Sea-Bird Scientific Navis, WMO number: 5905988) was deployed from the survey ship on August 16, 2018, near 50.22°N, 144.64°W. The Wirewalker was deployed twice from the process ship between August 23 and September 7. Incubations were performed on the process ship at least every two days from August 16 to September 7. Figure 2 shows the paths of all platforms throughout the cruise period. To intercalibrate O₂ sensors across platforms, Winkler titrations (Winkler, 1888; Carpenter, 1965) were performed on the survey ship with 16 CTD casts co-located with autonomous assets and the process ship. A full description of the field campaign can be found in Siegel et al. (2021). For consistency, mixed layer depth (MLD) was defined on all platforms by a density change of 0.1 kg m⁻³ with respect to 5 m (de Boyer Montégut et al., 2004; Suga et al., 2004). Methods relevant to this study are summarized in the following sections, with more detailed methods and calculations in Section S1.

Dissolved O₂/Ar was measured from the underway system (approximate 5 m depth) on both ships using Equilibrator Inlet Mass Spectrometry (EIMS; Cassar et al., 2009). Biological oxygen saturation anomaly was calculated following:

where O₂/Ar_meas is the ratio of O₂ and Ar measured in the seawater and O₂/Ar_sat is the ratio of O₂ and Ar in the air measured every 3–4 hours. MLDs were obtained from CTD casts on both ships and interpolated linearly between casts to obtain a 2-minute averaged MLD for each ship over the course of the cruise. NCP was calculated for the EIMS systems on the process (P) and survey (S) ships (NCP_EIMS(P), NCP_EIMS(S), respectively) following:

where $O2sat$ represents the oxygen concentration at saturation (μmol kg⁻¹) (Garcia and Gordon, 1992), ρ is the density of the seawater (kg m⁻³), and k is the weighted gas transfer velocity (m d⁻¹). Hourly wind speed data (10 m) were downloaded from the NOAA PMEL air–sea interaction mooring (https://www.pmel.noaa.gov/ocs/Papa). Gas transfer velocity was calculated following Wanninkhof (2014). Weighted piston velocities were calculated using these gas transfer velocities following methods in Reuer et al. (2007), with the modification presented in Teeter et al. (2018). Temperature and salinity data (required to compute these parameters) were collected from the ship’s underway instruments. O₂/Ar data from the first three weeks of the cruise on the process ship were not saved properly, due to an error with the instrument. Data for this ship were digitized in Matlab using screenshots of the instrument-collected data. The digitized data align well with properly saved data (Figure S1).

The O₂/Ar method is sensitive to uncertainties in the gas transfer velocity (Bender et al., 2011; Wanninkhof, 2014), the potential for vertical mixing of oxygen-depleted water from below the mixed layer (Cassar et al., 2014), and the loss of biological oxygen due to microbial respiration in the underway lines of the ship (Juranek et al., 2010). To account for respiration in the underway lines of the ship, 11 pairs (process ship) and 13 pairs (survey ship) of discrete samples for O₂/Ar were collected simultaneously from the ship’s underway system and from Niskin bottles collected at 5 m to obtain a correction factor, following methods outlined in Timmerman and Hamme (2021; Figures S2 and S3). The EIMS measurements also assume that argon is at saturation in the surface ocean (Li and Cassar, 2016). Additional EIMS methods are outlined in Section S1.1.

Data from three autonomous gliders were used to make estimates of NCP in the ML during the cruise period (August 15 to September 7, 2018). All gliders sampled every 4–6 hours down to 1000 m, with G276 and G469 sampling only on dives. All three gliders carried a CTD, an Aanderaa dissolved oxygen optode (4831), and an optical triplet puck for backscatter, chlorophyll fluorescence, and colored dissolved organic matter (SBE/Wetlabs ECO Puck). The seaglider also carried a multispectral radiometer (SBE OCR-504-ICSW). The gliders were not outfitted for oxygen optode air calibration (e.g., Nicholson and Feen, 2017) but were calibrated using a combination of Winkler titrations performed on the survey ship and intercalibration with the air-calibrated LF.

NCP can be estimated in the mixed layer using a mass balance model of dissolved oxygen. Oxygen concentration is predominantly controlled by air–sea gas exchange, entrainment of water during mixed layer deepening, vertical fluxes across the base of the mixed layer, and biological processes (Nicholson et al., 2008; Emerson and Stump, 2010; Alkire et al., 2012). The mixed layer consistently deepened over the cruise period, so the entrainment term was assumed to be the dominant flux term across the base of the mixed layer while diapycnal exchange was assumed to be negligible in comparison. The speed of atmospheric exchange is rapid compared to other physical processes, causing any horizontal gradients of oxygen to have minimal impact (Emerson et al., 1995). There is weak horizontal advection in the mixed layer near OSP (but more important at depth), shown by glider measurements by Pelland et al. (2018). Consequently, horizontal advection and diffusion fluxes were assumed to be small in comparison to other processes and were omitted. Additionally, for the seaglider, which followed the LF, the Lagrangian framework for the experiment was designed to minimize horizontal advective fluxes. The change in oxygen content over time is due to the sum of physical and biological fluxes:

where $MLD×d([O2])dt$ is the change in mixed layer oxygen inventory through time, [O₂] is mixed layer oxygen concentration (mmol O₂ m⁻³), F_gas (mmol O₂ m⁻² d⁻¹) is the sum of gas exchange terms, F_phys is the sum of vertical mixing fluxes, and $∫0MLD(P−R)dz$ is photosynthesis (P) and respiration (R) in mmol O₂ m⁻³ d⁻¹ integrated from the base of the ML to the surface. At steady-state, this difference is equivalent to NCP (Emerson et al., 2008). Rearranging, NCP can be represented by:

The first term on the right-hand side is measured directly while the other two fluxes are modeled as outlined below (Section 2.2.8).

Surface values of oxygen, potential temperature, and salinity were calculated as an average of the values from 3–10 m (standard deviations of 1.94 × 10⁻⁴ mol m⁻³ for oxygen, 6.03 × 10⁻⁴°C for potential temperature, and 5.00 × 10⁻⁴ for salinity; n = 281), and the final NCP was shown to be insensitive to the surface depth range selected. NCP was evaluated for each 24-hour period down to the deepest daily mixed layer.

The gas flux is the combination of diffusive air–sea gas exchange plus the contributions of both partially and completely collapsing bubbles (Emerson and Bushinksy, 2016). Meteorological data used with the gas exchange parameterization come from the NOAA PMEL mooring at Station P, the closest observations to the OOI gliders. The mooring data are also the best match to the wind speed observations collected by the survey ship, which was co-located with the Seaglider during the cruise period. For consistency, a weighted gas exchange velocity was calculated following methods outlined in Section 2.2.1. The individual fluxes were calculated for each run of a Monte Carlo simulation by randomly choosing between the parameterizations of Nicholson et al. (2016) or Liang et al. (2013), with an adjustment to the bubble flux from Emerson et al. (2019). The differences between the two are modest (1.49 mmol O₂ m⁻² d⁻¹, or 4.79% on average).

Wirewalker (WW) beam attenuation computed from drift-corrected transmissivity (Wetlabs C-Star) was used as a linear proxy for POC based on direct POC observations collected in bottle samples on the survey ship (Cetinić et al., 2012). Spatiotemporal matches were defined as WW profiles within 1 hour and 5 km of the CTD cast, and WW beam attenuation data were averaged for each depth at which bottle POC samples were collected.

For each local day from August 23 to September 7, diel gross carbon productivity (GCP) and community respiration (CR) were computed from the Wirewalker POC estimate within the mixed layer using the linear diel method from Barone et al. (2019). The method described in Barone et al. (2019) models a phytoplankton diel cycle (growth and losses), according to an idealized solar radiation curve (based on latitude, longitude and year-day), and optimizes the fit between this model and the Wirewalker data from the ML. This method assumes that there is negligible input of POC from below the ML and that, over short (daily) timescales, variations due to advection of gradients in POC are small. The diel variations in POC in the ML are therefore attributed to biological activity. The method therefore allows us to parse the rates of photosynthesis and respiration by modeling photosynthesis as a linear response to light intensity, E, and assuming that respiration is constant throughout the night and day. The optimized fit in the model provides an estimate of both GCP and CR, following Equation 3 in Barone et al. (2019):

where t is time in days (t₀ = time zero), POC₀ is the POC concentration at t₀, and E is the average irradiance during the day. Using this equation to isolate GCP and CR, we calculated NCP following:

Similar approaches for calculating rates from the diel cycles method were implemented by Siegel et al. (1989), Nicholson et al. (2015), White et al. (2017), Briggs et al. (2018), Freitas et al. (2020), and Rosengard et al. (2020). Negative GCP or CR values and days when the coefficient of determination (R²) for the diel POC trend was less than 0.5 were excluded from NCP calculations. To convert from volumetric to integrated rates, daily rates were multiplied by the average daily mixed layer value. NCP was computed using the difference of daily GCP and CR. One assumption in NCP_WW is that daily export is much less than daily NCP. We also assumed the daily DOC contribution to be negligibly small.

Confidence intervals were computed within the diel method code from Barone et al. (2019) by calculating the variance-covariance matrix. For an additional uncertainty metric, we bootstrapped the WW POC proxy to the profiling frequency of the Seaglider (4 hours, every 6th WW profile) prior to taking all the mixed layer observations into the diel fitting software. This approach resulted in 6 different estimates of daily NCP, for which the mean and standard deviation were computed.

The BGC-float was equipped with an SBE41 CTD, SBE63 optical dissolved oxygen sensor, Deep SUNA nitrate sensor (Johnson and Coletti, 2002), Deep-Sea DuraFET pH sensor (Johnson et al., 2017), and SBE ECO Triplet bio-optical sensor that measured chlorophyll fluorescence, optical backscatter (approximately 700 nm), and fluorescent dissolved organic matter. The float was programmed to execute a near daily mission cycle during the cruise period with 2 m vertical resolution for temperature, salinity, oxygen, and pH and 5 m resolution for nitrate above 150 m. The calibration and quality control for the pH and nitrate sensors followed the standard procedures proposed by the international Biogeochemical Argo program (Johnson et al., 2017). Because this float is not capable of making air oxygen measurements, the oxygen data derived from the Argo Global Data Assembly Center were quality-controlled using the climatology of oxygen from World Ocean Atlas 2018. This method usually gives rise to an accuracy of approximately 3% (Takeshita et al., 2013), which is not adequate for determination of air–sea oxygen gas flux. We recalibrated the oxygen sensor data using high-precision Winkler samples collected during the cruise period to achieve an accuracy of approximately 0.3%. Total alkalinity (TA) values were estimated using the CANYON-B global biogeochemical algorithm (Bittig et al., 2018), using float temperature, salinity, oxygen, and pressure measurements as input variables, along with location and time information. DIC values were computed from float-measured pH (total scale) and estimated TA using CO2Sys (van Heuven et al., 2011).

Lateral DIC transport was found to be small in prior studies from this region that took place over a longer time period (Emerson et al., 2011; Fassbender et al., 2016). By assuming the negligible role of horizontal advection in this region, due to the short time period of this study compared to the residence time of each tracer in the ML (Emerson and Stump, 2010; Bushinsky and Emerson, 2015), observed changes in the three tracers over time within the mixed layer can be written following Equation 8:

where the subscripts, Gas, Phys, and EP represent the change in tracer concentration induced by air–sea gas exchange, physical transport, and evaporation/precipitation, respectively. The NCP term was obtained by subtracting the abiotic terms from the observed change in tracer concentration in Equation 8. A detailed description of each abiotic term estimate is provided in Huang et al. (2022). Briefly, the air–sea gas exchange (Gas) was parameterized using the gas model of Wanninkhof (2014) for CO₂ and the model developed by Liang et al. (2013) for O₂. The evaporation and precipitation (EP) term was estimated indirectly by linking the chemical tracers with the rate of change of salinity within the mixed layer depth. The physical transport term (Phys), a sum of tracer change induced by entrainment, diapycnal diffusion, and wind-induced upwelling, were quantified by vertical tracer gradient from BGC-float measurements, combined with the meteorological and physical parameters from simultaneous satellite observations (i.e., wind stress) or prior studies (i.e., climatological diapycnal diffusivity reported by Cronin et al., 2015).

Triplicate 1 L seawater samples were collected in 1 M HCl-cleaned and Milli-Q water-rinsed clear polycarbonate bottles using a trace metal clean (TMC) CTD rosette aboard the process ship at six sample depths throughout the euphotic zone corresponding to 65%, 40%, 20%, 10%, and 1% of the incident irradiance (Io). To capture a peak in phytoplankton cell abundance as measured using flow cytometry during the cruise, sample depths were modified to 40%, 20%, 10%, 5%, and 1% Io from July 28 to September 1. Seawater samples from the TMC Go-Flo bottles were dispensed into the polycarbonate bottles following three seawater rinses within a TMC positive-pressure van.

Ambient NO₃⁻ concentrations were estimated by an in situ ultraviolet spectrophotometer (ISUS, Satlantic Inc.) attached to the underway system. Samples were inoculated with Na¹⁵NO₃ concentrations corresponding to approximately 10% of the ambient NO₃⁻ concentrations within a laminar flow hood kept in a positive-pressure TMC plastic bubble set up on the ship. Inoculated samples were then immediately placed in temperature- and light-controlled on-deck flow-through plexiglass incubators consistent with in situ surface water conditions by way of neutral density screening. Following 24 hours of incubation, samples were removed from the incubators, wrapped in dark plastic bags and stored at 4ºC until filtered. Samples were analyzed via mass spectrometry at the UC Davis Stable Isotope Facility. Atom % ¹⁵N and particulate nitrogen for all samples, blanks, and standards were measured, from which we calculated NO₃⁻ uptake following the methods described by Dauchez et al. (1995) and Slawyk et al. (1997). Samples for dissolved NO₃⁻ concentrations used for the ¹⁵NO₃⁻ uptake calculations were analyzed at the University of California Santa Barbara Marine Science Institute Analytical Lab.

NO₃⁻ NP was calculated as the total depth integration of NO₃⁻ uptake rates. Uptake rates were extrapolated to the surface (0 m) and exact ML depth based on trapezoidal integration:

where ρ is density (kg m⁻³), and $∫0ZMLN​15− ρNO3−$ is the NO₃⁻ uptake integrated over the euphotic zone and converted from nitrogen units to carbon units according to the Redfield ratio (Redfield, 1958; for more information, see Meyer et al., 2022). For this study, only ML rates are compared.

The triplicate samples exhibited low variability, with the standard deviations ranging from 0% to 26% with an average of 8%. See Meyer et al. (2022) for more information about sampling uncertainty.

A bottle incubation-based measurement of the relative rates of primary production and herbivorous grazing was used to estimate net community production (McNair et al., 2021). Phytoplankton growth rate and microzooplankton grazing rate from discrete depths throughout the mixed layer were measured concurrently in two-point dilution experiments (Landry and Hassett, 1982; Morison and Menden-Deuer, 2017) 17 times during the cruise. Water for the incubation experiments was collected from 1 to 3 depths within the mixed layer, close to midnight local time, using a standard CTD Niskin rosette (i.e., not TMC). Preparation of the diluted and whole seawater treatments followed prior protocols (Morison et al., 2019). Dilution experiment bottles (1 L) were incubated for 24 hours in the same incubators used for ¹⁵N incubations aboard the process ship.

Following the incubation, phytoplankton growth rate and microzooplankton grazing rate were determined using the net change of extracted Chl a concentration corrected for photoacclimation. The difference between the phytoplankton growth rate in the diluted and undiluted treatments provides an estimate of microzooplankton grazing. The growth rate of the phytoplankton in the absence of grazing is calculated by subtracting the grazing rate from the apparent phytoplankton growth rate in the undiluted treatment. Additional information about experimental design and calculations is detailed in McNair et al. (2021). Net community production estimates from the dilution experiments reflect the balance between phytoplankton growth and the consumption of phytoplankton by microzooplankton. This estimate of NCP does not include the carbon assimilated by the microzooplankton nor the consumption and respiration of meso- and meta-zooplankton, which are likely absent from the incubations as larger animals would be poorly represented in 1 L bottles.

NCP rates were estimated by first calculating Chl a production rates for each experiment in the mixed layer using Equation 10, which accounts for changing phytoplankton biomass during the incubation:

where k_z(d⁻¹) is the difference between phytoplankton growth rate and microzooplankton grazing rate, POC_{Chl a} is the concentration of phytoplankton carbon (mmol C m⁻³) at each depth, z, and t (d) is the length of the incubation. NCP was then integrated from the surface to the MLD using trapezoidal integration, which yields the estimated net accumulation of phytoplankton carbon within the mixed layer. Chl a concentration was converted to carbon using the POC:Chl a ratio measured for each day of the experiments.

Finally, in an effort to scale-up the in situ and in vitro measurements, and corroborate with global estimates, we also calculated NCP from satellite and modeling projects. Li and Cassar (2016) used genetic programming to model NCP using satellite products and built an algorithm to model NCP using net primary production (NPP) and sea surface temperature (SST) data from remote sensing. Eight-day averaged SST data, with 0.083° latitude and longitude resolution, were downloaded from the NASA Ocean Color website (https://oceancolor.gsfc.nasa.gov/). Monthly NPP data, with 0.167° latitude and longitude resolution, were estimated using both the Vertically Generalized Production Model (VGPM) (Behrenfeld and Falkowski, 1997) and the Carbon-based Production Model (CbPM; Westberry et al., 2008) retrieved from Ocean Productivity (http://sites.science.oregonstate.edu/ocean.productivity/standard.product.php). SST and NPP data were used to calculate NCP following Equation 7 in Li and Cassar (2016):

which is trained based on the global compilation of mixed layer NCP integration measured by O₂/Ar ratio anomaly. This equation was taken directly from Li and Cassar (2016).

ML estimates of in situ NCP and NP based on ML tracer budgets can be influenced by the mixing of tracers into the ML from below. Based on the data available, vertical mixing was calculated differently for different platforms. We assessed how the differences between each method impacted the results. Oxygen-based measurements (NCP_EIMS(P), NCP_EIMS(S), NCP_SG, NCP_OOIG276, NCP_OOIG469, and NCP_BGC(O2)) were impacted most by vertical mixing.

For the glider-based datasets, the physical fluxes across the base of the ML were split into an entrainment flux and a diapycnal diffusion flux. The entrainment flux was calculated as the change in the MLD over time (⁠$ΔzΔt$⁠; m d⁻¹), multiplied by the O₂ concentration difference between the entrained water mass and the ML value from the previous day (Emerson et al., 2008):

where ΔMLD = MLD₁ − MLD₀, the change in MLD for the time range Δt = t₁ − t₀, and $[O2]ML,t=0$ is the concentration of O₂ in the ML in the previous time step. [O₂]_ent was determined for each time step by examining the previous profile and taking the average O₂ concentration in the depth range of ΔMLD. When the ML deepened, the sign of the entrainment flux was dependent on the oxygen gradient at the base of the ML. When the ML shoaled, there was no O₂ entrainment flux. The ML consistently deepened over the course of the cruise, resulting in a positive flux from the sub-ML O₂ maximum into the ML (Figure S4). The diapycnal diffusion flux was calculated as the product of the oxygen gradient across the ML and the diapycnal coefficient, K_z. K_z at the base of the ML was inferred from the monthly climatology by Cronin et al. (2015). For the BGC-float datasets, the vertical mixing term at the base of the ML is composed of wind-induced Ekman pumping, entrainment due to the change of ML, and diapycnal diffusion:

where w is the Ekman pumping velocity derived from wind stress fields following Signorini et al. (2001) and O_2base is the oxygen concentration at the base of the ML, $O2¯$ is the averaged oxygen concentration within the ML, K_z is the diapycnal coefficient at the base of the ML (adopting the monthly climatology estimated by Cronin et al. (2015) using an ML heat budget), and $∂O2∂z$ is the oxygen gradient below the ML. The glider datasets incorporate the entrainment and dyapycnal diffusion terms, but historically have not included an Ekman pumping term in the vertical mixing calculations. In this case, the Ekman term was low and did not cause a large discrepancy between the two datasets.

For the Wirewalker dataset, other studies incorporated mixing terms to estimate the transfer of POC into or out of the mixed layer (Rosengard et al., 2020); however, with the strong stratification at the base of the mixed layer this term is assumed to be negligible. For incubation-based methods, vertical mixing was not applicable because samples were collected in the ML.

The thirteen NCP and NP datasets used in this study are outlined in Table 1. Oxygen measurements were converted to carbon units using a ratio of 1.4:1 (Laws, 1991), while nitrate measurements were converted to carbon using 6.6 for the ¹⁵N incubations (Redfield, 1958) and 7.8 for the BGC-float measurements (Haskell et al., 2020). To make the datasets comparable and test for statistically significant differences, each NCP and NP estimate was daily-averaged or interpolated to a 24-hour time step. Because the ¹⁵N incubations (NP_inc) were performed from roughly 6 AM to 6 AM local time (UTC = 9), all other measurements except the dilution incubation dataset were interpolated or averaged to this time period as closely as possible. The dilution incubations were performed from midnight to midnight local time. Five datasets were defined as “Lagrangian” assets: (1) NCP_EIMS(P), (2) NCP_SG, (3) NCP_WW, (4) NP_inc, and (5) NCP_inc. Deviations from the “Lagrangian” water mass by the process ship were removed from NCP_EIMS(P) data prior to taking the daily average using timestamps from the Ship’s Log. The other assets were not considered as part of the Lagrangian experiment because they did not specifically follow the LF that was deployed to mark the water mass of interest at the beginning of the cruise.

First, a Kruskal-Wallis test was performed to determine if the datasets differed significantly from each other. The Kruskal-Wallis test was chosen for this comparison because it is a rank-based non-parametric test that does not make prior distributional assumptions. While the area covered by this study was relatively small, the platforms in this study were not always co-located, which could lead to differences in the NCP or NP estimates. Additionally, because these datasets contain time-series data with repeated measures, the estimates may be biased due to the presence of spatial autocorrelation, temporal autocorrelation, or both. To facilitate comparison between the datasets, a generalized additive mixed model (GAMM) was fit using the daily-averaged values from each dataset to model NCP or NP as a function of the time and location of each measurement. A GAMM was used instead of a generalized linear mixed model (GLMM) because the inclusion of non-parametric smooth functions in a GAMM relaxes the assumption of a linear relationship between the response variables (NCP and NP) and the predictors (time, location). This feature in turn provides a more flexible means by which our model can identify and correct for the non-constant variance and possible temporal and spatial autocorrelation structures within each dataset.

First, we built a base model to determine how each measurement method affects the NCP. Each dataset was given a unique categorical variable “type” (NCP_EIMS(P), NCP_EIMS(S), NCP_SG, NCP_OOIG276, NCP_OOIG469, NCP_WW, NP_inc, NCP_inc, NP_BGC(NO3−), NCP_BGC(O2), or NCP_BGC(DIC)) such that each dataset was modeled separately. This base model took the form of:

where NCP and NP are the modeled values as a function of type_i with a linear link function, plus an error term, ε_i. The term α_i represents the intercept of the model, and β_i is the coefficient modifying type_i. The subscript i denotes the individual type of each method. In this base model, each dataset is modeled alone based on their individual factor “type”; therefore, we tested each dataset against the mean NCP of all the datasets to determine if they are significantly different.

Subsequently, we built three models from the base model to test the effects of space and time on NCP and NP, separately and together. The first model tested for the effect of the day the data were collected:

where f(day_i) is a spline term for the continuous variable in the model for day number on which the data were collected. The second model, built off the base model, tested the effects of spatial heterogeneity on the modeled NCP or NP. This model included the daily-averaged distance of the platform from the LF (which represented the “center” of the Lagrangian experiment):

where f(distance_i) is the spline term for the average daily distance of the asset from the LF. Modeled separately, the two continuous variables day_i and distance_i allowed us to assess if the day the measurements were collected or the location of the platform during the data collection was driving any of differences between the different NCP and NP datasets. Finally, we built a third model, where day and distance were modeled jointly:

In each model, the effects of spatial and temporal autocorrelation on modeled NCP or NP were accounted for by incorporating both smooth functions in Equations 14b and 14c. We tested the goodness of fit for each model against the base model using ANOVA. Finally, we picked the model that best homogenized the data, and a Kruskal-Wallis non-parametric test was performed on the residuals of the final model to determine if the datasets were significantly different after accounting for spatial and temporal drivers of differences between datasets. Satellite measurements were not included in either GAMM because the spatial resolution of the measurements was too coarse.

The goal of the EXPORTS 2018 field deployment was to better constrain the BCP and understand the mechanisms driving carbon export in the high-nutrient, low-chlorophyll (HNLC) region at Station P (Martin and Fitzwater, 1988; Harrison, 2002). In this regard, the EXPORTS field deployment provided a unique opportunity to compare methods that differ in approach, assumptions, and integration scales to estimate NCP and NP, and therefore carbon export potential. The NCP measurements ranged from −42.1 to 66.4 mmol C m⁻² d⁻¹ with a mean of 5.8 mmol C m⁻² d⁻¹ (n = 143). The NP measurements ranged from −1.9 to 21.8 mmol C m⁻² d⁻¹ with a mean of 8.16 mmol C m⁻² d⁻¹ (n = 17). Due to the variability seen across each method over time (Figure 3), comparing the net carbon export potential averaged from each method over the course of the cruise is important. The range, median, and mean of NP and NCP in each dataset are presented in Table 2. In general, the different methods indicated that carbon export potential from the system was low around Station P, with slight net autotrophy throughout the course of the cruise. For reference, highly productive regions of the ocean, such as coastal waters, can have NCP values over 100 mmol C m⁻² d⁻¹ (Wang et al., 2020). Below, we examine potential causes for differences between methods, but first contextualize our study by comparison to historical observations of NCP and NP in the region. Where needed, more specific results from each method are summarized in Section S2.

The measurements in this study were collected near Station P, where biogeochemical observations, including those of NCP and NP, have been conducted for many years. Station P is a long-term mooring site with over 60 years of oceanographic data (Freeland, 2007) and is of special interest when investigating the ocean carbon cycle due to the annual net flux of CO₂ into the ocean in the subarctic North Pacific (Ayers and Lozier, 2012). Previous studies in the area have shown that Station P is net autotrophic in the summer and heterotrophic in the winter, with an NCP maximum in August (Plant et al., 2016), resulting in an average mixed layer annual NCP of 2 mol C m⁻² year⁻¹ (Giesbrecht et al., 2012; Fassbender et al., 2016; Plant et al., 2016). The historic summertime NCP measurements from the Northeast Pacific have ranged from 6.3 to 18.2 mmol C m⁻² d⁻¹ (Table 3). Chl a measurements from the EXPORTS field deployment in August and September 2018 (average of 0.25 µg L⁻¹; Howard et al., 2010; Siegel et al., 2021) were lower than reported in the area previously (approximately 0.4 µg L⁻¹) for the late summer (Harrison, 2002).

Consistent with past studies in the area, the results of our study show net autotrophic conditions over the course of the cruise from most platforms. Corresponding with the lower Chl a measurements during this field deployment, average NP and NCP estimates from our study are also slightly lower than previously reported values in the literature.

Overall, all datasets showed that NCP and NP around Station P in the summer of 2018 were relatively close to zero. The medians of all datasets were positive, except for NCP_inc, though NCP_WW and NCP_BGC(DIC) also had negative values (Figure 3). Some datasets had smaller variance (NCP_EIMS(P), NCP_EIMS(S), NCP_SG, NCP_OOIG276, NCP_OOIG469, NP_inc, NCP_BGC(NO3−)), while other datasets had a much larger spread (NCP_inc, NCP_WW, NCP_BGC(DIC), NCP_BGC(O2); Figure 4a). We performed a Kruskal-Wallis test which showed that the datasets were significantly different from each other (chi-squared = 78.4, df = 10, p-value = 1.01 × 10⁻¹²).

We built three models in total, in addition to the base model outlined in Equation 14a. Of these three, only the model that accounted for temporal autocorrelation and controlled for the effects of spatial heterogeneity on the NCP or NP estimates led to a result that was significantly different from the base model (Table 4). Accounting for temporal effects in the datasets alone, or the spatial and temporal effects together, did not significantly influence the relationship between each method (ANOVA p = 0.41, p = 0.21, respectively). Similarly, correcting for only spatial or temporal autocorrelation structures in the individual datasets did not lead to a significantly better fit than the base model (further information in Supplemental Material Section S2.7). However, the model that accounted for the spatial effects on NCP or NP and controlled for temporal autocorrelation in the individual datasets performed better than the base model (ANOVA, p = 0.048; Table 4). We performed a Kruskal-Wallis test on the residuals of this model and found that the datasets were not significantly different (p = 0.35). This finding means that, with the caveats listed below, the 11 in situ and in vitro NCP and NP measurements from the EXPORTS North Pacific campaign were statistically the same, after accounting for the effects of spatial heterogeneity and removing temporal autocorrelation structures in the individual datasets (Figure 4b).

Because we averaged each dataset over the course of a day and compared to the daily-average location of the LF, our analysis includes some spatial integration that could reduce the effect of spatial heterogeneity across assets. Moreover, our method measures the distance of each asset to the LF and not the actual distance between each measurement location. Additionally, because of the vast difference in temporal resolution of each method, we could not compare across methods without taking a daily average of each dataset. This averaging drastically reduced the size of our final dataset and overall reduced the efficacy of our modeling efforts. While this approach was useful to compare across different methods, it could be especially useful in comparing much larger datasets, from the global BGC-Argo array, for instance.

While satellite images indicate significant spatial variability in surface Chl a concentration in our sampling region (Figure 2b), this oceanic system is not highly physically dynamic. In a more physically dynamic area, the presence of hydrographic frontal regions could introduce significant gradients in measured surface properties among different datasets. Conversely, spatial gradients are expected to be smaller in oligotrophic gyres, such that a greater distance between assets may have less influence on the resulting measurements. In future studies, especially in more physically dynamic areas, such as the North Atlantic or the Southern Ocean, characterizing spatial gradients before comparing measurements will be critical.

While the results from the datasets agree well and the results of the GAMM indicate that spatial heterogeneity of the measurements explains the largest discrepancies between methods, other factors that may lead to differing results from each method. These factors are considered in the following sections.

Some of the differences seen across measurement methods can be explained by the inherent differences in temporal and spatial scales of integration. NCP_EIMS(P) and NCP_EIMS(S) integrate over the residence time of O₂ in the mixed layer (Timmerman and Hamme, 2021), on the order of several days to weeks, though these measurements are heavily weighted to the recent past (Teeter et al., 2018). Over the EXPORTS cruise period, the average integration time was 9 days for the O₂/Ar datasets (based on an average MLD of 30.5, and gas transfer velocities of 3.4 m d⁻¹). Due to the Lagrangian nature of the process ship, NCP can be calculated from the O₂/Ar measurements that accounts for long-term trends in O₂ in the ML, and therefore reflects a “real-time” NCP value that is not integrated over 9 days (Hamme et al., 2012). However, to be comparable with the dataset from the survey ship, which was not Lagrangian, NCP_EIMS(P) was not calculated to find the real-time values. NCP from the seagliders was measured based on daily changes in oxygen. However, due to noise in each budget calculation, the NCP measurements were calculated using fluxes smoothed with a 30-day moving average. As such, signals due to shorter time-scale changes in oxygen will be dampened, and the measurements are more representative of NCP over weekly to monthly temporal scales. This longer temporal integration of the O₂-based measurements decreases the variability in the resulting NCP estimates.

Conversely, the residence time for the ML tracers NO₃⁻ and DIC is on the order of seasonal or even annual timescales, as they are dependent on the scale of the replenishment mechanism in the tracer pools (i.e., the vertical/horizontal mixing for both DIC and nitrate and air–sea exchange for DIC; Sarmiento and Gruber, 2006). However, the tracer budget method used for the BGC-float calculated a change in the tracer (O₂, DIC or NO₃⁻) over time. In this case, the BGC-float measurements were interpolated to 24 hours, so each daily measurement represents a 24-hour NP or NCP. The measurement of nitrate incorporation into particulate organic matter (POM) collected through the NP_inc will not be affected by the long residence time of nitrate in the ML and will reflect NP rates integrated over 24 hours. While NP_inc measures the uptake of nitrate in the experiment over time, the ML must be in steady state for NP to be a measurement of carbon export potential. The incubation rates do not spatially integrate; they reflect production and respiration rates of the planktonic community from the sampling location. Wirewalker measurements of NCP were calculated from GPP and CR rates measured from diel cycles of POC in the mixed layer. As such, NCP_WW measurements reflect daily integrated values as well. These differences in the spatial and temporal integration scales could explain some of the differences between NCP and NP data collected by these various methods. As expected, methods with 24-hour resolution showed more temporal variability than those with longer temporal integration (Figures 3 and 4). The technical uncertainty of rate estimates for NCP_inc is well below the variation observed here (Morison et al., 2017). The variability in NCP_WW and NCP_inc likely reflects true oscillations between production and loss in the ML, which are muted when tracers integrate over longer time scales. Thus, NCP_inc and NCP_WW are not inherently less precise; they are capturing different variability than other methods.

An additional source of error between datasets is in situ deviation from Redfield stoichiometry. To convert from O₂ units to C units, we used a ratio of 1.4 to 1 (Laws, 1991). If we assume that these ratios could fluctuate in situ, we can perform a rudimentary analysis to see how sensitive the measurements are to these deviations. Laws (1991) identified a possible range of 1.1 to 1.8 for the photosynthetic quotient (PQ) used to convert O₂-based measurements to C units, though the author hypothesized that 1.8 is likely too high. If we apply a PQ of 1.1 to our O₂-based measurements, it increases our estimates by approximately 27%; if we use a PQ of 1.8, it decreases our estimates by approximately 22%. Even if we assume that 1.1 and 1.8 are extreme upper and lower bounds on the PQ and use 1.25 and 1.6, our measurements would change by approximately 12%. More commonly, the PQ used in publications is 1.4 (as used in this study) or 1.45 (Fassbender et al., 2016). If we convert the measurements in this study using a ratio of 1.45, our NCP estimates decrease by approximately 3.4%.

The BGC-float nitrate uptake measurements reflect both dissolved and particulate forms of organic matter production, while the ¹⁵N incubations only include measurements of POM production. Dissolved organic matter (DOM) production is thought to be characterized by a higher C:N ratio (approximately 14) compared with the nutrient ratio of POM produced (approximately 6.6; Redfield, 1958). In our study region, Haskell et al. (2020) observed the spring and summer drawdown of DIC and nitrate at a stoichiometric ratio of 7.8 and inferred that organic carbon production was composed of 70% POM and 30% DOM. As such, we used 6.6:1 (C:N) to convert the measurements from the incubations to C units, and a ratio of 7.8:1 for the BGC-float measurements (Haskell et al., 2020). Laws (1991) presents a lower bound of 5.7:1 and an upper bound of 9:1 for the ratio of C:N. Following the same calculations as above, the lower bound would decrease our measurements by as much as 26% and the upper bound would increase our measurements by up to 36%. Previous studies in the region have shown that NP derived from ¹⁵N incubations results in lower estimates of export production than estimates from shipboard O₂/Ar measurements, likely due to the fact that the incubation method does not measure production of dissolved organic nitrogen (Giesbrecht et al., 2012; Timmerman and Hamme, 2021). Following Haskell et al. (2020) and using a C:N conversion of 7.8, NP_inc would increase by approximately 18%. Biological processes including nitrogen fixation and denitrification, and physical processes including advection, can influence the in situ C:N ratio (Tyrrell and Lucas, 2002; Carter et al., 2021). However, there is no evidence that these processes influenced the C:N ratio in this study.

The average ML POC:Chl a ratio during the cruise, 19 µmol C:µg Chl a, was used to convert the NCP_inc rate to units of carbon, which is slightly greater than the global range, 12–15.5 (Legendre and Michaud, 1999). Using a POC:Chl a conversion less than 19 decreases the NCP_inc, bringing it closer to zero net production. Thus, the stoichiometric ratio used to convert to C units can create up to 36% of the differences between measurement methods.

There is a long-standing debate in the scientific community about incubation-based measurements versus in situ observations of NCP or NP, especially in oligotrophic areas where rates are low and close to the detection limit of each method (Marra, 2009; Duarte et al., 2013; Ducklow and Doney, 2013; Williams et al., 2013). The agreement between the in situ and ¹⁵NO₃ labelled incubation measurements, and NCP_inc indicate that uncertainties and biases associated with in situ observations and in vitro measurements (e.g., bottle artifacts) did not significantly influence the measurements in our study. While environmental parameters can be controlled using in vitro methods, incubations can lead to biases caused by factors that include unnatural light regimes (Godoy et al., 2012), underestimates of grazing by zooplankton (Robinson and Williams, 2005), and incorrect assumptions about respiration rates throughout a 24-hour period (Bender et al., 1987). The NCP_inc varied from day to day, but overall the cruise average indicated that phytoplankton growth and grazing were well balanced, leading to an NCP estimate close to zero. The NCP estimated from the dilution incubations is similar, but slightly lower than the other measurements in this study. All else being equal, NCP_inc would tend to be lower than other NCP measurements because it reflects the net change in phytoplankton and does not account for any secondary production (i.e, microzooplankton biomass) or fecal pellet production. Additionally, because NP is based on the incorporation of isotopes, the production rate will always be >0 (Bender et al., 1999), which would result in greater NP_inc estimates than NCP_inc. NP estimates, however, were well within the range of the other estimates of NCP. NP, like the dilution incubation, measures the production of biomass fueled by NO₃⁻ to equate NP to NCP, assuming that the system is in steady state. The different variables being measured in the two incubation experiments (i.e., NO₃⁻ uptake and change in Chl a concentration for NP and NCP, respectively) lead to different methodological caveats. Specifically, NO₃⁻ uptake relies on accurate estimation of ambient nitrate concentrations in situ and accurate conversion from carbon to nitrogen units, whereas the dilution experiments assume that phytoplankton growth is unaltered by dilution and that grazing rate is a function of encounter rate between predator and prey.

The most difficult physical flux to constrain in this study was the vertical mixing flux at the base of the ML. The vertical mixing flux broadly includes the contributions of the tracer into or out of the ML due to MLD changes, diapycnal diffusion and wind-induced Ekman pumping. A slight supersaturation of oxygen in the ML (3%) was observed during the cruise period. Directly below the ML, supersaturation increased as a result of heating in the upper thermocline and biological production (noted before in Shulenberger and Reid, 1981; Spitzer and Jenkins, 1989; Emerson et al., 1995; Nicholson et al., 2015). Mixing into the ML would therefore increase the concentration of O₂ and might lead to higher estimates of NCP from O₂-based methods (NCP_SG, NCP_OOIG276, NCP_OOIG469, NCP_BGC(O2)) and O₂/Ar-based methods (NCP_EIMS(P), NCP_EIMS(S)) if the O₂ saturation anomaly below the mixed layer is biological. In this discussion, we focus on the effect of vertical mixing on the O₂-based measurements, because the fluxes could be large enough to impact NCP estimates due to the oxygen maximum at the base of the ML. For the O₂ budget, the flux of O₂ into the ML by vertical mixing is roughly 40% of our estimated production by NCP_SG, NCP_OOIG276, NCP_OOIG469 and NCP_BGC(O2). This vertical mixing term is not incorporated into the NCP estimates from the EIMS, indicating that these measurements are likely an overestimate, though the magnitude of the overestimate is difficult to constrain because we do not know the proportion of biological O₂ coming from depth. There was no associated POC or nitrate maxima below the ML, and our analysis shows negligible impacts of vertical mixing on measurements based on tracers other than O₂. Indeed, the concentration of POC generally decreases with depth, indicating that NCP_WW would underestimate true NCP due to vertical mixing.

The capacity for continuous measurements of vertical profiles by the gliders and BGC-float provides the opportunity to parameterize the vertical mixing flux (Equations 10 and 11). Nevertheless, the parameterization of vertical mixing between glider and BGC-float is slightly different (see Sections 2.2.1 and 2.2.3), so discrepancies across methods can result in substantially different estimates of NCP. For the two gliders that made observations the entirety of the cruise period, the means (n = 24) for total physical flux were 5.87 ± 0.46 (Seaglider) and 7.20 ± 0.36 mmol C m⁻² d⁻¹ (OOI G276,). Of these totals, the means (n = 24) for entrainment fluxes were 4.23 ± 0.59 (Seaglider) and 5.25 ± 0.57 mmol C m⁻² d⁻¹ (OOI G276), and the diapycnal diffusion accounted for 1.65 ± 0.23 (Seaglider) and 1.95 ± 0.52 mmol C m⁻² d⁻¹ (OOI G276). OOI 469 was deployed later in the cruise period and did not capture the most rapid deepening of the mixed layer. The slight difference in entrainment fluxes likely arises as a result of spatial heterogeneities, as OOI G276 and G469 ranged from 28 to 85 km away from the Seaglider during the cruise period. The MLD was roughly 5 m shallower for the OOI assets and had a stronger O₂ maximum below it (Figure S4). Over the cruise period, the consistent deepening of the mixed layer into this O₂ maximum suggested that the entrainment flux incorporated a portion of the diapycnal diffusive flux. Consequently, the glider-based measurements may be an underestimate when considering both fluxes (Equation 9).

In the BGC-float parameterization, more explicit formulations were applied, with the vertical term partitioned into the entrainment, diapycnal diffusion and wind-induced Ekman pumping. The sum of averaged vertical mixing estimated by BGC-float was 4.78 ± 2.1 mmol C m⁻² d⁻¹ (2 mmol C m⁻² d⁻¹ for the entrainment, 2.7 mmol C m⁻² d⁻¹ for the diapycnal diffusion and 0.08 mmol C m⁻² d⁻¹ for the wind-induced Ekman pumping; n = 6), aligning with the estimate of 4.20 ± 0.59 mmol C m⁻² d⁻¹ by the gliders (n = 24). The O₂/Ar-based measurements of NCP from the EIMS cannot account for vertical mixing of water into the ML. The inverse relationship between N₂O and O₂ at depth can be leveraged to account for vertical mixing of O₂-depleted waters into the surface ocean (as proposed by Cassar et al., 2014; Izett et al., 2018). Briefly, without a flux of N₂O from depth, we would expect the surface ocean N₂O concentration to be in equilibrium with the atmosphere. Because N₂O concentration generally increases with depth and O₂/Ar generally decreases with depth, we can use depth profile measurements of each tracer to calculate a linear relationship between N₂O and O₂/Ar at depth. Finally, if N₂O is supersaturated in the surface water, we can use the previously derived slope to estimate what the effect would be on the O₂/Ar ratio in the ML. During this field campaign, there were no drastic changes in MLD, indicating no significant upwelling. While depth profiles of N₂O and O₂/Ar were collected from the process ship, the sampling resolution around the base of the ML was not high enough to properly quantify the biological supersaturation of O₂ directly below the ML. As such, NCP_EIMS(P) and NCP_EIMS(S) have not been corrected for vertical mixing fluxes and thus are generally higher than the O₂-based NCP estimates. In future research campaigns, we recommend collecting more frequent discrete N₂O and O₂/Ar measurements with a finer depth resolution to correct EIMS measurements for vertical mixing.

Finally, as a thought exercise, we can put an upper bound on the vertical mixing correction by comparing our NP and NCP estimates from the ML to export at depth. Various measurements can be used to estimate POC flux at depth, but here we have used estimates from the EXPORTS 2018 field deployment reported in Buesseler et al. (2020). The flux at 50 m (F₅₀) is a combination of $FPOCML′$⁠, the POC flux at 50 m that can be attributed to the ML after accounting for the attenuation that occurs between the base of the ML and 50 m, and $FPOCMLto 50m$⁠, the POC flux that is generated between the base of the ML and 50 m, following Equation 15a:

Rearranged, the proportion of POC flux at 50 m that can be attributed to the ML follows Equation 15b:

For simplicity, we assumed that there is no attenuation acting on $FPOCMLto 50m$ (which will result in a slight underestimate of $FPOCML′$⁠). We used the Martin Curve (Equation 1; Martin et al., 1987) to back-calculate the POC flux at the base of the ML before POC attenuation, following Equation 16a:

where $FPOCML​$ is the POC flux at the base of the ML, z_ML is the average MLD, and b is the coefficient taken from table 1 in Martin et al. (1987). The net POC production that occurred between the base of the mixed layer and 50 m (⁠$FPOCML to 50m$⁠) was estimated by multiplying the measured NPP by the e-ratio (⁠$NCPNPP)$⁠, following Equation 17:

where R is the e-ratio. Combining Equations 16a and 17, the flux of POC at the base of the ML can be calculated following Equation 16b:

Over the course of the cruise, NCP will be equal to the net flux not only of POC but also DOC from the ML to depth, plus any changes in POC and DOC standing stock in the ML (⁠$d(POCML),d(DOCML)$⁠, respectively), following Equation 18a:

where $FDOCML$ represents the net flux of DOC out of the ML. Conceptually, if the standing stock of POC in the ML declines over the experimental period, the NCP estimated from POC export flux will be an overestimate, because it is reflecting a flux of organic carbon to depth that was produced prior to sampling. Conversely, if the standing stock of POC increases over the experimental period, NCP estimated from POC export flux will be an underestimate. Because this argument is based on POC flux at depth, the portion of DOC production that stays in the ML and the portion that is exported to depth do not matter for these calculations, so we can combine $FDOCML+d(DOCML)$ into one term, $FDOCNCP$⁠, following Equation 18b:

The proportion of DOC produced from NCP can be estimated experimentally and represented mathematically as Equation 19:

where $(d(DOC)d(POC+DOC))NCP$is the percentage of ML NCP that results in DOC production. Adding this percentage to Equation 18b gives Equation 18c:

Following algebraic rearranging gives Equation 18d:

For the sake of argument, we applied Equation 15a to data collected throughout the EXPORTS 2018 field deployment to estimate a lower bound on C flux from the ML. Buesseler et al. (2020) estimated a mean sinking particle flux around Station P during the EXPORTS 2018 field deployment of 5.5 mmol C m⁻² d⁻¹ at 50 m (F₅₀). For the purposes of these estimates, we used average MLD of 27 m. NPP was measured on the cruise through incubation experiments. Briefly, NPP estimates were measured during incubation experiments where NPP was calculated as DI¹³C uptake over 24 hours following the addition of 180 µM NaH¹³CO₃, reflecting 10% ambient DIC concentrations. The cruise average NPP below the ML was 0.25 mmol C m⁻³ d⁻¹ (Meyer et al., 2022). Integrated between 27 m and 50 m results in an NPP of 5.75 mmol C m⁻² d⁻¹. Assuming an e-ratio of 0.25, then $FPOCMLto 50m$ was approximately 1.4 mmol C m⁻² d⁻¹ (Equation 17). Using these values in Equation 16a, $FPOCML′$ was approximately 4.1 mmol C m⁻² d⁻¹. In Equation 16b, applying b = −0.863 (table 1 in Martin et al., 1987), we estimated that $FPOCML​$ was approximately 7.0 mmol C m⁻² d⁻¹. Moving to Equation 18a, $d(POCML)$ was derived experimentally during the 2018 EXPORTS field deployment as approximately 0.8 mmol C m⁻² d⁻¹, and $(d(DOC)d(POC+DOC))NCP$ was taken from the literature. Bif and Hansell (2019) estimated that around Station P, approximately 26% of NCP is released as DOC into the ML (i.e., $(d(DOC)d(POC+DOC))NCP=0.26$⁠). Using these values in Equation 18d, we estimated that, based on POC flux estimates at depth, NCP in the ML was approximately 10.5 mmol C m⁻² d⁻¹. This rough estimate must be a lower bound on the carbon flux at 27 m, because $FPOCML to 50m$ will also be affected, necessarily, by flux attenuation. Based on these rough calculations of export flux at the base of the ML, we have increased confidence in our ML estimates of NCP. This exercise indicates that, overall, our NP and NCP estimates from all platforms are in a similar range as the C export from Station P.

To conclude our evaluation of NP and NCP measurement methods, we also assessed how the field-based measurements compare to remote sensing and modeling techniques. Changes in the BCP caused by climate change could lead to impacts on the natural and anthropogenic carbon cycle. Therefore, establishing the baseline capabilities of the BCP to monitor for future change is important. Due to the limitations of in situ measurements, satellite and model-based estimates of NCP are powerful tools for understanding the global oceanographic carbon cycle and the role that the BCP plays in removing CO₂ from the atmosphere. Though the temporal and spatial resolution of satellite- and model-based measurements is coarser than the in situ and shipboard measurements obtained during the EXPORTS 2018 field deployment, there is good agreement between all measurement methods over the months of August and September, 2018. Our mean NCP estimates following Li and Cassar (2016) are 9.6 ± 1.3 and 10.6 ± 2.0 mmol C m⁻² d⁻¹ (n = 24) for the VGPM- and CbPM-based models respectively, with the caveat that the Li and Cassar (2016) model was developed for use with the VGPB-based model, and has not been tested with the CbPM model NPP output.

In this study, we leveraged multiple concurrent measurement methods to constrain carbon export potential around Station P during the EXPORTS 2018 field campaign. Many methods used the same tracer but from different platforms (NCP_EIMS(P) and NCP_EIMS(S); NCP_SG, NCP_G276, NCP_G469, and NCP_BGC(O2); NP_inc and NCP_BGC(NO3−)), or employed similar methods using different tracers (NCP_SG, NCP_G276, NCP_G469, NCP_BGC(O2) and NCP_BGC(DIC), and NCP_BGC(NO3−); NP_inc and NCP_inc). Many measurements were based on estimating biological O₂ production in the ML, either using an Ar tracer to parse out the physical and biological components (NCP_EIMS(P) and NCP_EIMS(S)), or using a budget approach (NCP_SG, NCP_G276, NCP_G469, and NCP_BGC(O2)). Some methods (NCP_EIMS(P) and NCP_EIMS(S)) integrate over a relatively long timescale, while others (NCP_SG, NCP_G276, NCP_G469, NCP_BGC(O2), NCP_BGC(DIC), and NCP_BGC(NO3−)) measure NP and NCP as the change in the in situ tracer over time. Overall, all datasets indicated that carbon export potential from the ML throughout the EXPORTS 2018 field campaign was relatively low, with slight net autotrophy dominating the region throughout the month of the experiment. GAMM results indicate that the primary driver of differences between the datasets is spatial variability around Station P.

In addition to comparing the datasets outright, we characterized the ways in which the results from each method might differ from each other, including in situ deviations from the stoichiometric relationship used to convert to C units, differences in methods of calculating vertical mixing, and differences in temporal scale of integration for each tracer and method. Methods that integrate over larger temporal and spatial scales generally had a smaller variance, while methods that integrate over 24 hours showed more variability from day to day, possibly reflecting biological variation around the mean rates.

Despite the considerable spatial variation in surface ocean Chl a concentration, Station P is generally a low productivity, low variability system. It contrasts with more dynamic areas like the North Atlantic or the Southern Ocean. In these more variable areas, we expect that the multi-method approach may be particularly important to understanding spatial variation of NCP and characterizing carbon export potential in the region.

Data for this project can be found in SeaBASS (10.5067/SeaBASS/EXPORTS/DATA001). Data from the BGC-float (WMO ID: 5905988) can be found at: https://www.mbari.org/science/upper-ocean-systems/chemical-sensor-group/floatviz/.

The supplemental files for this article can be found as follows:

Supplemental Material. PDF

The authors thank Abhishek S. Bhatia, Ed Iversen and the Fall 2020 Introduction to Statistical Consulting class for assistance with statistical analysis. We thank Robert Izett, Ross McCulloch, Cara Manning, Erinn Raferty and Jai Cunningham for assistance with dissolved gas analysis in the Tortell Lab at the University of British Columbia and the Hamme Lab at the University of Victoria. We thank Hans Gabathuler and Yajuan Lin for technical field assistance. We acknowledge David Siegel, Ivona Cetinic, Andrew Thompson and others in the EXPORTS team for science leadership, data processing and collaboration. We thank the captain and crews of the R/V Roger Revelle and R/V Sally Ride. We thank two anonymous reviewers for their feedback and we thank the editors of the journal.

AKN, NC, AM, MGM, WT, and WG were funded through NASA Grant 80NSSC17K0552. RH was funded by Natural Sciences and Engineering Research Council of Canada (NSERC) grant number 328290-2017. PT was funded by NSERC. HMM and SMD were funded through NASA grant 80NSSC17K0716. DN, ST, MMO, and MF are funded on NASA grant 80NSSC17K0663. MMO and MF are also funded on NASA grant 80NSSC17K0662. AJF and YH were supported by NOAA’s Global Ocean Monitoring and Observing Program. This is PMEL Contribution No. 5291.

The authors declare no competing interests.

• Contributed to conception and design: AKN, ST, YH, MF, MGM, HMM, DN, AJF, MMO, AM, SMD, NC.

• Contributed to acquisition of data: AKN, ST, YH, MF, MGM, HMM, DN, AJF, MMO, AM, SMD, WT, WG, RH, NC.

• Contributed to analysis and interpretation of data: AKN, ST, YH, MF, MGM, HMM, DN, AJF, MMO, AM, SMD, WG, RH, NC.

• Drafted and/or revised the article: AKN, ST, YH, MF, MGM, HMM, DN, AJF, MMO, AM, SMD, WT, WG, PT, RH, NC.

• Approved the submitted version for publication: AKN, ST, YH, MF, MGM, HMM, DN, AJF, MMO, AM, SMD, WT, WG, PT, RH, NC.