Methods

Unless specified otherwise, measurements were taken on board the R/VSikuliaq , cruise number SKQ201617S, from 07 January 2017 through 13 January 2017 at a single station 16.5°N 106.9°W, which was located in an oligotrophic region of the Eastern Tropical North Pacific Oxygen Deficient Zone (ETNP Station P2; Figure 1A). Data are compared against measurements taken at 16.5°N 152.0°W on 08 May 2015, collected on the GO-SHIP CLIVAR/CARBON P16N Leg 1 Cruise (CCHDO Hydrographic Cruise 33RO20150410). This station was at the same latitude as ETNP Station P2, west of the ODZ, but was not anoxic (P16 Transect Station 100; Figure S1).

Water property measurements

We measured water properties of temperature, salinity, fluorescence, oxygen concentration and turbidity using the shipboard SeaBird 911 CTD. Auxiliary sensors included a WetLabs C-Star (beam attenuation and transmission) and a Seapoint fluorometer. Data were processed with Seabird Software, (programs–data conversion, align, thermal mass, derive, bin average and bottle summary) using factory supplied calibrations. Data was analyzed and visualized in R (Team 2011). Processed data are available under NCEI Accession number 1064968 (Rocap et al., 2017).

Water mass analysis

Evans et al. (2020) previously employed optimum multiparameter analysis to map the percent identity of the water observed at each depth to three water masses: the 13 Degree Celsius Water (13CW), North Equatorial Pacific Intermediate Water (NEPIW), and Antarctic Intermediate Water (AAIW). We subset and examine only the portion of these data that correspond to our site.

Acoustic Measurements

Backscattering signals from the ship-board EK-60 were collected and archived by UNOLS as raw data files. We used Echopype software (Lee et al., 2021) to convert these raw files to netcdf files, which were down-sampled to five minute time-step resolution, saved as a text file, and later visualized in R. The acoustic data appeared to be off by one hour and so one hour was subtracted from all time measurements. This correction resulted in zooplankton vertical migrations being synchronized with the diel light cycle as was recorded on board the ship by JAC.

Particle size measurements

Particle size data were collected by an Underwater Vision Profiler 5 (UVP) that was mounted below the CTD-rosette and deployed for all CTD casts shallower than 2500 m. A UVP is a combination camera and light source that quantifies the abundance and size of particles from 100 μm to several centimeters in size (Picheral et al., 2010). UVP data were processed using the Zooprocess software, which prepares the data for upload to the Ecotaxa database (Picheral et al., 2017); data from all UVP instruments are located on this online database for ease of access. Detailed descriptions for installation of the Zooprocess software can be found on the PIQv website (https://sites.google.com/view/piqv/zooprocess-uvpapp). Zooprocess uses the first and last image number selected by the user in metadata to isolate the downcast and process this subset for both particle size distribution and image data. The processed files and metadata are then uploaded to a shared FTP database where it is available for upload to Ecotaxa. This project required the extra step of filtering out images due to the discovery of an issue with the lighting system, where only one of the two LEDs would illuminate, resulting in an incomplete sample. The filtering procedure is documented in a link available at the same location as the Zooprocess download. Images where only a single light illuminated were removed from the dataset before it was uploaded on to Ecotaxa. Once uploaded to Ecotaxa, data were downloaded from EcoPart (the particle section of the database) in detailed TSV format, and analyzed in R. The UVP provided estimates of abundances of particles in different size-bins, as well as information about the volumes over which those particle numbers had been calculated. Particles were categorized into bins starting at 102-128 μm in size, with the width of each particle size bin 1.26 times larger than the previous bin. We observed particles in 26 distinct size bins, with largest, mostly empty, bin covering particles from 26-32 mm.
The instrument is capable of observing smaller particles (down to 60 μm), but these tend to be underestimated and so we only consider particles ≥102 μm in this analysis. The instrument can in principle also measure larger particles (up to the field of view of the camera), though these tend to be scarce enough to be not detected. In this paper, we do not have an upper size cut-off for our analysis and rather implement statistics that are robust to non-detection of scarce large particles (section 5.5.1). Visual inspection of images larger than 1 mm suggests that these large particles are primarily “marine snow” but about 5% are zooplankton. We did not quantify the size distribution of these images.

Flux measurements

Free floating, surface tethered particle traps were used to quantify carbon fluxes from sinking particles. Arrays, consisting of a surface float and two traps, were deployed and allowed to float freely during which time they collected particles. Trap deployments began on 07 January, concurrently with the beginning of the UVP sampling, and continued through 12 January. Trap recovery began on 08 January and continued through 13 January. Trap depths spanned the photic zone and mesopelagic, with the shallowest at 60 m and the deepest at 965 m. Trap deployments lasted between 21 and 91 hours, with deeper traps left out for longer, to collect more biomass. Two types of traps were deployed. One set of traps, generally deployed in shallower water, had a solid cone opening with area 0.46 m2. The second set had larger conical net with opening of 1.24 m2 area made of 53 μm nylon mesh similar to the description in Peterson et al. (2005). All equipment were combination incubators and particle traps, but in this study we only use trap data. No poisons were used, and both living and dead zooplankton, or ‘swimmers’, were manually removed prior to POC analysis.
Sediment trap material was filtered immediately upon trap recovery onto pre-combusted GF-75 45 mm filters (nominal pore size of 0.3 µm) and preserved until further analysis at -80°C. These filters were split into several fractions for other analyses not discussed here. Total carbon content of particles in each trap were measured by isotope ratio mass spectrometry. Elemental analyses for particulate carbon and nitrogen quantities as well as 13C and 15N isotopic compositions were conducted at the U.C. Davis Stable Isotope Facility (http://stableisotopefacility.ucdavis.edu) on acidified freeze-dried trap samples to capture organic elemental contributions. Carbon was below mass spectrometry detection limits in four traps – these traps were excluded from further analysis. Traps at similar depths did detect carbon, lending confidence to the idea that these non-detections were technical in nature, due to splitting of samples for multiple analyses, rather than reflecting environmental conditions.

Analysis

Particles were binned by depth with 20 m resolution between the surface and 100 m, 25 m resolution between 100 m and 200 m depths and 50 m resolution below 200 m. This increasing coarseness of the depth bins helped account for more scarce particles deeper in the water column, while maintaining higher depth resolution near the surface. To perform this binning, particle numbers, and volumes of water sampled of all observations within each depth bin were summed prior to other analyses. Most analyses focused on the mesopelagic, defined here as the region between the base of the secondary chlorophyll maximum layer (160 m — hereafter the base of the photic zone), which is within the ODZ, and 1000 m.
Two normalized values of particle numbers were calculated. In the first, particle numbers were divided by volume sampled, to generate values inparticles/m3 . In the second, particles were divided by both volume sampled and the width of the particle size-bins to generate values in particles/m3/mm .

Particle size distribution

We determined the slope and intercept of the particle size distribution spectrum by fitting a power law to the data, which is a common function for fitting particle size distributions (Buonassissi & Dierssen, 2010). Because large particles were infrequently detected, we used a general linear model that assumed residuals of the data followed a negative-binomial (rather than normal) distribution. We fit the equation
\(\ln\left(\frac{E\left(\text{Total}\,\text{Particles}\right)}{\text{Volume}*\text{Binsize}}\right)=b_{0}+b_{1}\,ln\left(\text{Size}\right)\)(Eqn 1).
to solve for the Intercept (\(b_{0}\)) and particle size distribution slope (\(\text{PSD}=b_{1}\)). On the left-hand side of Eqn 1.E(Total Particles) refers to the expected number of particles in a given depth and particle size bin assuming a negative binomial distribution of residuals (Date, 2020; Ooi, 2013). Volumeindicates the volume of water sampled by the UVP, or in the case of depth-binned data, the sum of the volumes of all UVP images in that depth interval. Binsize indicates the width of the particle-size bin captured by the UVP. Thus, if particles between 0.1 and 0.12 mm are in a particle size bin, the Binsize is 0.02 mm. On the right-hand side of Eqn 1, Size corresponds to the lower bound of the particle size-bin. We use the lower bound of a particle size-bin, rather than its midpoint, because, due to the power-law particle size distribution slopes, the average size of particles in each size-bin is closer to the size-bin’s lower bound.

Estimating particle flux

We estimated particle flux throughout the water column, by fitting particle data to trap measurements. We assumed that particle flux in each size bin (j) followed the equation
\(Flux=\sum_{j}{[\frac{\text{Total}\,\text{Particle}s_{j}}{\text{Volume}*\text{Binsiz}e_{i}}*C_{f}*({\text{Size}_{j})}^{A}]}\)(Eqn. 2)
Such that flux at a given depth is the sum of all size-bin specific values.
We used the optimize() function in R ’s stats package to identify values for the Cf and Acoefficients in Eqn 2. that yielded closest fits of the UVP estimated flux to each particle trap.
We also estimated the exponent of the particle size to biomass exponent\(\alpha\) and size to sinking speed exponent \(\gamma\) per the equations \(\text{Biomas}s_{j}\sim Size_{j}^{\alpha}\,\) and\(\text{Spee}d_{j}\sim{S\text{iz}e}_{j}^{\gamma}\). This is done by assuming a spherical drag profile, in which case \(A=\alpha+\gamma\)and \(\gamma=\alpha-1\) (Guidi et al., 2008); with “A” referring to the exponent in Eqn 2.

Size specific information

We separately analyzed total particle numbers, particle size distribution, and particle flux for particles larger than or equal to 500 μm, and those smaller than 500 μm, to determine the relative contributions of these two particle classes to particle properties. 500 μm was chosen as it has been previously defined as the cutoff point between microscopic “microaggregates” and macroscopic “marine snow” (Simon et al., 2002).

Variability

To explore the timescales of temporal variability in the POC flux, we determined how well the flux at each depth horizon can be described by the sum of daily and hourly temporal modes. This was achieved by fitting the general additive model of form
\(\text{Flu}x^{1/5}\sim s\left(\text{Depth}\right)+s\left(\text{Day}\right)+s\left(\text{Hour}\right)\)(Eqn. 3)
This model explored whether estimated flux levels appeared to vary by decimal day and decimal hour, holding the effects of depth constant, in the 250 m to 500 m region. The smooth terms s for Depthand Day were thin plate splines, while the s term forHour was a cyclic spline of 24-hour period.

Smoothing for Comparison to Model Results

Normalized particle abundance data, from the only UVP cast that traversed the top 2000 m of the water column, taken on January 13 at 06:13, was smoothed with respect to depth, time, and particle size using a general additive model of the form
\(\ln\left(\frac{E\left(\text{Total}\,\text{Particles}\right)}{\text{Volume}*\text{Binsize}}\right)\,\sim\,s\left(Depth,\,ln\left(\text{Size}\right)\right)\)(Eqn. 4)
In this case, there is a single, two-dimensional, smooth term, rather than additive one-dimensional terms as in Eqn. 3 so that the smooth term can consider interactions between the two parameters, rather than assuming that the terms are additive. The predicted particle numbers at each particle size and depth, as well as particle size distribution spectra, and estimated particle masses of all particles smaller than 500 μm and all particles larger than or equal to 500 μm were then compared to each of Weber and Bianchi’s (2020) models, corresponding to our H1-H3 .

Modeling remineralization and sinking

To quantify disaggregation, our goal was to compare the particle size-abundance spectrum at each depth to a prediction of the null hypothesis, that it is simply governed by the effects of sinking and remineralization reshaping the spectrum observed shallower in the water column. This prediction is generated using the particle remineralization and sinking model (PRiSM), modified from DeVries et al. (2014), which we applied to the shallower spectrum as an initial condition. The difference between the null hypotheses prediction and observation indicates the role of processes not accounted for in PRiSM, such as disaggregation, aggregation, and active or advective transport of particles with a different size spectrum than the ones seen at the deeper depth.
In practice we expanded the previous numerical implementation of PRiSM to allow for particle size distribution spectra with particle-size bins that match those obtained by the UVP, and to return estimates of the number of particles in those same size bins (Text S1). The model accepts inputs of particle size distributions at each depth, and changes in particle flux between each depth-bin and the next, deeper, depth-bin. The model optimizes a particle remineralization rate that would result in that observed flux loss. It finally returns a “predicted” particle size distribution spectrum that has total flux equal to the flux of the observed deeper spectrum that would be expected if the shallower spectrum only sank and remineralized. In cases where flux increased with depth, particles are assumed to put on mass rather than lose mass following a negative remineralization rate. Here, “negative remineralization” stands in for chemoautotrophy, active transport, and other processes that result in flux increases with depth. While these processes likely have more complex effect on the particle size distribution than is accounted for in our model, we note that flux increases with depth are very rare, and that allowing for negative remineralization allows our null model to be robust in those rare cases.