Open Ocean Particle Flux Variability From Surface to Seafloor

The sinking of carbon fixed via net primary production (NPP) into the ocean interior is an important part of marine biogeochemical cycles. NPP measurements follow a log‐normal probability distribution, meaning NPP variations can be simply described by two parameters despite NPP's complexity. By analyzing a global database of open ocean particle fluxes, we show that this log‐normal probability distribution propagates into the variations of near‐seafloor fluxes of particulate organic carbon (POC), calcium carbonate, and opal. Deep‐sea particle fluxes at subtropical and temperate time‐series sites follow the same log‐normal probability distribution, strongly suggesting the log‐normal description is robust and applies on multiple scales. This log‐normality implies that 29% of the highest measurements are responsible for 71% of the total near‐seafloor POC flux. We discuss possible causes for the dampening of variability from NPP to deep‐sea POC flux, and present an updated relationship predicting POC flux from mineral flux and depth.

expensive (and commensurately sparse), but also because these measurements are highly variable (Buesseler & Boyd, 2009;Buesseler et al., 2007;Mouw et al., 2016). Estimates of the global export of organic carbon out of the upper ocean vary by more than a factor of three Boyd & Trull, 2007;Henson et al., 2011), with fluxes to the deep ocean being less constrained.
One long-standing question is whether the majority of mass and energy is delivered as a steady background level or from rare deposition events (Rokop, 1974;Smith et al., 2014). Addressing this question is not straightforward. Although NPP and deep-sea flux are undoubtedly linked, predicting flux from NPP is challenging on a point-by-point basis because such observations decouple over space and time and because flux-generating processes are not directly a function of primary production (Buesseler, 1998;Plattner et al., 2005; Footnote 1 in the supporting information). Yet, quantitative understanding of the vertical fluxes of carbon and other elements requires characterizing their variability. Capturing the probability distribution of deep-sea flux allows one to quantify how episodic it is, and can provide a simple yet comprehensive way to describe the flux of carbon, its variability, and its relationship with other fluxes.  showed that NPP is robustly well-characterized by a log-normal probability distribution, meaning its extensive variability can be described by just two parameters (the log-mean μ and the log-standard deviation σ; Footnote 2 in the supporting information). Similarly as for many other natural phenomena (Limpert et al., 2001), this distribution is thought to emerge because NPP is the product of both structured (e.g., average nutrient supply to the surface) and stochastic (e.g., turbulent mixing) processes. POC export out of the upper ocean also follows a log-normal distribution, though the statistical power behind this statement is limited by sample size .
Here, we show that the log-normal character of NPP propagates down to deep-sea measurements of POC flux, further demonstrating that the log-normal is a robust feature of variations in ocean carbon fluxes. We show that log-normality holds not only for a collection of global open ocean measurements, but also for measurements from individual time-series, indicating its consistency across scales. Furthermore, log-normality holds not only for POC but also for calcium carbonate and opal, which are thought to play an important role in the vertical transport of POC (Armstrong et al., 2001;Wiedmann et al., 2020). So-called "ballast" minerals are known to be strong predictors of POC flux, although it is still disputed whether this effect is mechanistic or coincidental (Footnote 3 in the supporting information; Passow & De La Rocha, 2006). Our results suggest (1) an updated global relationship between mineral and POC fluxes, (2) that POC deposition to the seafloor, and hence supplied to benthic communities, is episodic, but not enough so to meet the 80/20 Pareto principle benchmark for a heavily imbalanced phenomenon in other contexts (Sanders, 1987), and (3) that marine heterotrophs may dampen variability in NPP by consuming a higher fraction of NPP when NPP is larger (Footnote 4 in the supporting information).

Global Open Ocean Fluxes
The POC flux through the water column at a particular location and time follows a flux profile ( ) z  (mg C m −2 d −1 ) that attenuates with depth (Suess, 1980). This implies that an initial flux from some reference depth z o , is attenuated through the water column by an "attenuation factor" a so that the flux making it to the seafloor We may further define an export efficiency ef (Eppley & Peterson, 1979) as the ratio between POC flux and NPP such that Since the log-normal emerges from the product of multiple random variables, one might expect that if POC export ( ) o z  or NPP are log-normally distributed , we should expect s  to be log-normally distributed as well, as one is just multiplying by more terms. Log-normality would propagate through depth as long as a (or a × ef) scales with ( ) o z  (or NPP) in some way, is roughly constant, or randomly varying according to some well-behaved probability distribution (i.e., a probability distribution with finite variance). Log-normality would not extend to s  if, for instance, a × ef had a more complex dependence on NPP such as if a × ef sharply increases above a threshold NPP value, or a × ef ≈ 0 for large NPP values (SI). Given that attenuation and export efficiency are both quite variable (Cael & Follows, 2016;Henson et al., 2012), we hypothesize that this log-normality of NPP should propagate all the way to the seafloor flux s  (Footnote 5 in the supporting information).
This hypothesis is readily testable on the global scale thanks to the comprehensive database of particle flux measurements assembled by Mouw et al. (2016). We extracted from this database all of the measurements which were made at or below 1 km depth, and within 1 km of the seafloor, totaling 3,337 sediment trap measurements of bathypelagic POC flux ( OC  , [mg C m −2 d −1 ]), which we consider sufficiently near the seafloor to make inferences about the statistics of seafloor fluxes (Footnote 6 in the supporting information). As OC  is generally tightly correlated with mineral fluxes (Armstrong et al., 2001;Francois et al., 2002;Henson et al., 2012;Klaas & Archer, 2002), we apply the same selection criteria for fluxes of calcium carbonate (calcite and aragonite, IC  ) and opal ( Si  ), resulting in 2,968 IC  measurements and 2,821 Si  measurements (we also corroborate these correlations here). While these data are spatially biased, for example, toward under-representing polar or shelf regions ( §Episodicity) we take them to be a representative random sample of near-seafloor bathypelagic particle fluxes for the global open ocean. We test the log-normality of these data using the standard Kolmogorov-Smirnov cumulative distribution function (CDF) metric   1 2 max ( ) ( ) D P P    , where P 1 and P 2 are the CDFs of  ; here we are comparing the empirical CDF of the data itself with a hypothesized log-normal distribution's CDF. As in  we fit a log-normal distribution to each of these three measurement sets ( OC  , IC  , and Si  ), estimating the log-moments by minimizing this fitting statistic. We use D because it is the simplest such statistic (Stephens, 1974); other choices yielded the same results.
We find that all three of these measurements are well-described as log-normal, with D values well below the thresholds to reject this hypothesis at the 5% significance level (Stephens, 1974). Figure 1 shows the CAEL ET AL.   Mouw et al. (2016), with fitted log-normal distributions overlaid (colored lines). D is the Kolmogorov-Smirnov statistic; (μ, σ) are the log-mean and log-standard deviation of the log-normal distributions corresponding to the colored lines. Bottom: Percentiles of the empirical (x-axis) versus log-normal (y-axis) distributions from the top panels. Solid black line is the 1:1 line and the dotted lines represent a standard 30% measurement error (Buesseler et al., 2000(Buesseler et al., , 2007Stanley et al., 2004). hypothesized and empirical CDFs as well as the percentile-percentile plots which indicate that the log-normal approximation predicts the percentiles of all three measurement fluxes' distributions within a standard measurement error of 30% (Buesseler et al., 2000(Buesseler et al., , 2007Stanley et al., 2004;Footnote 7 in the supporting information) except for the smallest percentiles of the mineral fluxes where the absolute deviations are still small (<0.3 mg C m −2 d −1 ). Altogether there is clear statistical evidence that globally all of these fluxes are log-normally distributed, and errors introduced by this approximation are either within standard measurement error and/or very small.

Time-Series Fluxes
Is this observed log-normality just a coincident feature of where and when these measurements have historically been taken, or of large-scale biogeochemical cycles? Or is it a robust feature of ocean particle fluxes? The argument presented in the previous section is, in principle, equally applicable on smaller scales. A stricter test of our hypothesis is whether it holds for individual locations, that is, data from time-series stations; encouragingly, these fluxes are known to be episodic (Karl et al., 2012;Smith et al., 2018), consistent with log-normality.
We apply the log-normal fitting procedure from the previous section to data from six deep-sea sediment trap time-series locations: • A Long-term Oligotrophic Habitat Assessment (ALOHA) in the North Pacific Subtropical Gyre, a system which despite being perennially oligotrophic exhibits surprising variability in community structure and export Karl et al., 2012) and even interannual variations in nutrient limitation  • The Carbon Retention in a Colored Ocean (CARIACO) time-series in the Caribbean Sea, a near-coastal site also affected by terrestrial processes (Muller-Karger et al., 2019) • Station K2 in the subarctic Northwest Pacific, a mesotrophic system with important dust-borne inputs of the limiting micronutrient iron (Honda, 2020;Lam & Bishop, 2008) • Station M in the Northeast Pacific, a site within the California Current, an upwelling system that is productive and highly variable not only seasonally but also on short timescales as well as interannually (Smith et al., 2018) • The Oceanic Flux Program (OFP) time-series in the Sargasso Sea (Conte & Weber, 2014;Conte et al., 2001Conte et al., , 2003Conte et al., , 2019, an oligotrophic system with significant blooms and interannual variability (Lomas et al., 2013;Steinberg et al., 2001) • The Porcupine Abyssal Plain (PAP) time-series in the North Atlantic Ocean, a temperate system exhibiting spring blooms, but also strong influences by turbulent ocean features as well as climatic drivers (Lampitt et al., 2010) Further details about each time-series can be found in the above references. Collectively, the near-seafloor sediment traps from these time-series span measurements from a diversity of ocean settings; for example, eutrophic/oligotrophic environments, nutrient and oxygen status, or seafloor depth (Lampitt et al., 2010 . Some of these sample sizes are somewhat small to rigorously test for log-normality for individual data sets, but collectively they can demonstrate the robustness (or lack thereof) of the log-normal distribution on different spatial scales and at different locations; our conclusions are not sensitive to the inclusion/exclusion of any particular time-series. Figure 2 shows the results of fitting the individual time-series' data. As for the global data sets, the deviations between the empirical and theoretical distributions are not statistically significant, and the percentiles tend to stay within standard measurement error and/or represent small absolute deviations in the low tails, one exception being the very highest percentiles of the PAP IC  distribution, which is likely attributable to this having the smallest sample size (N = 166), and therefore these percentiles being very sensitive to the exact values of the largest two to three measurements. The other exception is the highest few percentiles of the ALOHA Si  distribution, which are more intriguing and suggest that the processes governing the diatom-diazotroph-association-dominated summer export pulse Karl et al., 2012) in the North Pacific Subtropical Gyre may be capable of "breaking" this distributional description of particle flux and generate deep-sea fluxes substantially larger than would be otherwise predicted. Outside of these anomalously large opal fluxes, however, these time-series data corroborate the robustness of the log-normal description of deep-sea particle fluxes.

Mineral Scaling
As mineral fluxes and POC fluxes are both log-normally distributed, this implies that a scaling relationship is a more appropriate way to relate the two fluxes than multi-linear regression (Armstrong et al., 2001;Francois et al., 2002;Henson et al., 2011;Klaas & Archer, 2002). The global / / OC IC Si  measurements present an opportunity to update the classic multi-linear regression equations for mineral flux, though it is important that variation in depth is also accounted for, as depth-normalization can introduce substantial uncertainty Olli, 2015). For these regressions, we restrict the global measurements of Mouw et al. (2016) shown in Figure 1 to the 2,798 coincident measurements of all three fluxes and fit the function: where κ and γ are the parameters of a standard scaling relationship, that is, y = κx γ , b is the standard "Martin curve" exponent, and β then represents the ratio of organic carbon association per unit mass for opal versus calcium carbonate, which some may interpret as a ballast effect. In effect, this equation is asking how OC  , normalized (Footnote 2 in the supporting information) to a uniform depth of 3,500 m by a power law (Martin et al., 1987), scales with mineral flux. We also assume multiplicative (i.e., % rather than absolute) errors, consistent with these variables scaling with one another and being log-normally distributed. As both OC  and   IC Si   have errors, a model II regression is necessary; we use major axis regression, but our results are insensitive to this choice. We identify the (β, b) pair that yields the best fit regression between ( / 3500) b OC z  and   IC Si   , and repeat this procedure 1,000 times with bootstrap resampling to estimate uncertainty in all of the parameters (κ, β, γ, b). Figure 3 shows the results of fitting this model; we see that a scaling relationship is an effective description of this relationship (r 2 = 0.66 for the dashed line in Figure 3; Footnote 9 in the supporting information), with depth-normalized POC flux scaling with mineral flux to the γ = 1.06 ± 0.01. This exponent is slightly -but statistically significantly -larger than 1, corresponding to a slightly systematically increasing POC/mineral ratio with increasing mineral flux. Furthermore, we find β = 0.26 ± 0.01 suggesting that roughly four times as much organic carbon is associated with calcium carbonate per unit mass relative to opal. Interestingly, we find a slightly high (Footnote 10 in the supporting information) b value of 1.16 ± 0.03. CAEL ET AL.
10.1029/2021GL092895 5 of 10  Figure 1, but for the individual time-series' data. Color indicates station and shape indicates measurement. Solid black line is the 1:1 line and the dotted lines represent a standard 30% measurement error (Buesseler et al., 2000(Buesseler et al., , 2007Stanley et al., 2004).
Analogous scaling relationships are found for the four time-series that have measurements of both minerals alongside POC flux -ALOHA, CARIACO, K2, and OFP -though the depth-normalization parameter b does not need to be estimated as each time-series' measurements are made at a single depth ( Figure S4). CARIACO, K2, and OFP have similarly good fits (r 2 = 0.81, 0.88, and 0.84, respectively) and slightly sub-linear scaling relationships (i.e., non-linear with a scaling exponent lower than unity: γ = 0.93 ± 0.03, 0.85 ± 0.02, and 0.92 ± 0.02, respectively), corresponding to a slightly systematically decreasing POC/mineral ratio with increasing mineral flux. Interestingly, CARIACO, K2, and OFP also have much larger β values (β = 1.45 ± 0.33, 1.84 ± 0.40, and 1.63 ± 0.46, respectively), indicating that per unit mass, opal has more organic carbon associated with it per unit mass in these locations. This underscores the important role that both minerals play, and that their relative importance can vary. The ALOHA time-series yields a much higher γ = 1.40 ± 0.11 and β = 15.6 ± 10.9, and a poorer fit (r 2 = 0.47). The large β value further underscores the importance of opal flux depending on the time and place -it is known to be critical for deep export at Station ALOHA (Karl et al., 2012) -but also is highly unconstrained, so we caution against its overinterpretation. The spread in scaling exponents γ between the time-series may be indicative of spatial variation in the mineral-POC flux dynamics around the global, approximately linear relationship; the global γ = 1.06 is close to the average (1.03) of the time-series' γ values.
As the global β we find is 0.26 ± 0.01, these findings suggest that perhaps aggregates containing, for example, diatom cells are ingested and remineralized more than, for example, aggregates containing coccolithophore cells and plates (liths). We cannot say from this work if these biominerals are truly ballasting POC because in some situations, these refractory biominerals may be all that remains after the labile POC is remineralized into the water column. Nonetheless, Figures 3 and S4 demonstrate that depth-normalized POC flux does scale approximately linearly with a weighted sum of these minerals; the scaling relationship is also quite strong given the variability in location and depth considered, and also appears to hold similarly for individual locations, though the relative weights of calcium carbonate versus opal differ between the global case and the local cases.

Dampening
We find an interesting consistency in the relationship between variability in NPP and variability in deep-sea flux, globally and across time-series, which is also consistent with the NPP-shallow flux relationship found in . We estimate the log-standard deviation σ of NPP at each site in the same manner as above and as in , using 14 C measurements taken as a part of the Hawai'i Ocean Time-series (HOT; Karl & Lukas, 1996), CARIACO time-series, and Bermuda Atlantic Time-series Study (BATS; Lomas et al., 2013;Steinberg et al., 2001) respectively for those locations, and all of the measurements in the data set used in Buitenhuis et al. (2013) and  within 5° of the PAP and M sites for those locations; global σ for NPP (σ NPP ) is taken from  (Footnote 11 in the supporting information). As in , we take the statistics of 14 C measurements (Footnote 12 in the supporting information) to be reflective of the statistics of NPP, noting that all of the measurements used herein are subject to caveats, ambiguities, and uncertainties. We then compare the σ for deep-sea POC flux (   ) with the corresponding σ NPP in Figure 4. In  we showed that if shallow flux scales with NPP and both are log-normal, the scaling exponent should correspond to a constant ratio of     / NPP  , and we found α ∼ 0.65 ± 0.14. The same is true here for deep flux; we regress the deep flux data in Figure 4 using the model    NPP  and find α = 0.62 ± 0.04 (r 2 = 0.89), in excellent agreement with the previous value. This implies that deep flux scales similarly with NPP to shallow flux, which in turn implies deep flux scales roughly linearly with shallow flux. As α < 1, this also implies that all of the processes that collectively transform NPP eventually into deep-sea flux tend to dampen the variability in NPP. Ecologically, this CAEL ET AL.
10.1029/2021GL092895 6 of 10 dampening effect could likely be caused by enhanced grazing, or heterotrophy in general, when NPP is high; if heterotrophs tend to consume a higher fraction of NPP when and where NPP is high, this would result in a negative correlation between a × ef (see §Global open ocean fluxes) and NPP, leading to    NPP  (Campbell, 1995). Such ecological dampening by the heterotrophic community would likely occur at shallower depths, that is, in the epipelagic and upper mesopelagic, where flux attenuation and organism abundances are highest (Iversen et al., 2010;Jackson & Checkley, 2011;van der Jagt et al., 2020). This is corroborated by the comparatively small differences between upper ocean export σ values globally and for HOT, CARIACO, and BATS/OFP and their corresponding   values Figure 4), and the closeness of the α values derived from shallow and deep fluxes versus NPP. In Figure 4, σ for shallow flux (flux measurements at ∼150 m, a.k.a. export) is always between that of deep flux and NPP, but closer to that of deep flux, supporting the argument that most of the dampening of variability occurs at shallow depths. It has been argued that zooplankton act as gatekeepers for export flux at the base of the mixed layer (Jackson & Checkley, 2011); there is experimental evidence that zooplankton-phytoplankton dynamics follow such dampening patterns and that this dampens carbon export variability (van der Jagt et al., 2020) but this dampening has not been demonstrated on the global scale. Globally representative measurements similar to those made in van der Jagt et al. (2020) would be valuable for testing whether similar mechanisms underlie this variability dampening at the global scale. Because NPP and flux observations are taken over different space and timescales, however, it is not clear how much of these σ differences can be attributed to ecological processes. The difference in measurement scales should reduce deep-sea   relative to shallow   , and shallow   relative to σ NPP ; regardless of this measurement-produced offset, however, the coherence in the differences between σ values in Figure 4 does imply a consistent relationship between the variability in NPP and the variability in deep-sea POC flux across scales and locations.

Episodicity
Is the benthos supplied by episodic flux events or a steady rain of particles? This question can be assessed quantitatively in terms of how heavy the tail of the probability distribution of fluxes is. Log-normals are one of a class of "heavy-tailed" distributions for which values much larger than the mean are not uncommon, but are not as heavy-tailed as other probability distributions commonly found in nature, such as a power-law distribution. The degree of heavy-tailedness or episodicity can be quantified simply in terms of the joint ratio (Footnote 13 in the supporting information): The (100 − X)/X such that X% of the largest values contribute (100 − X)% of the sum of all the values. For the global OC  data, both the observations and the log-normal fit yield a joint ratio of 71/29, corresponding to an imbalanced distribution, but a less imbalanced one than the classical 80/20 "Pareto principle" benchmark of a heavily imbalanced distribution in social, economic, and other natural systems (Sanders, 1987). This also holds for minerals and for time-series; all of the global distributions' joint ratios were between 70/30 and 75/25, and all of the time-series' distributions' joint ratios were between 60/40 and 75/25 (Footnote 14 in the supporting information). This implies an intermediate answer to the episodicity question; the deep-sea flux distributions are heavy-tailed enough that high-flux pulses do play an important role in supplying the benthos, but not heavy-tailed enough that such pulses dominate total supply. This measure of episodicity is however only in terms of the (∼weekly or longer) time intervals over which these deep flux measurements are taken. There will be variations in flux during these time intervals and so the true flux will be more episodic than the calculated episodicity from these time-averaged values. Furthermore, other relevant supplies to the benthos, such as energy (Grabowski et al., 2019) or genes linked to particular substrates (Boeuf et al., 2019), appear to be more episodic than organic carbon; the composition, nature and value of the supplied material may therefore be more episodic as well. Finally, these global and time-series data are heavily representative of temperate and subtropi-CAEL ET AL.  cal systems. Polar or shelf systems are known to be extremely episodic, sometimes with total export/supply dominated by a single regular seasonal event (e.g., Amiel and Cochran (2008), Ducklow et al. (2015), and references therein). We thus expect different degrees of episodicity in these systems, and potentially non-log-normally distributed fluxes.

Discussion
It is important to recognize that while the propagation of log-normality from NPP, to export out of the upper ocean, to deep-sea flux is intuitive when cast in terms of an export efficiency and an attenuation factor; in reality particle fluxes are very complicated and that this over-simplification appears to work so well is rather surprising. Deep-sea traps integrate over large spatial areas (Siegel & Deuser, 1997) that can be strongly influenced by lateral advection and stirring. Productivity and even export can vary greatly on much smaller scales (Estapa et al., 2015;Mahadevan, 2016). Particles sinking with a wide range of settling velocities contribute appreciably to total flux (Trull et al., 2008) meaning a deep-sea flux measurement is also a convolution of productivity at different times; deep-sea particle fluxes are also measured on timescales of weeks, unlike those of NPP or upper ocean export, which are measured on timescales of a day or days, respectively. NPP occurs throughout the euphotic layer and its depth-dependence varies (Behrenfeld & Falkowski, 1997). All of these factors complicate any one-to-one relationship between the value of a NPP measurement and the value of a corresponding deep-sea flux measurement. Nonetheless, our results indicate that these complications do not prevent log-normality from propagating into deepsea fluxes. We also note that the log-normal description is agnostic as to the driver(s) of these fluxes and their variations.
To conclude, we have shown that the log-normality of NPP measurements propagates from upper-ocean fluxes through the water column to near-seafloor fluxes, both globally and for individual time-series, both for organic carbon and for opal and calcium carbonate. This provides a simple way to characterize the variability of these fluxes, the organic carbon associations of these minerals, and the link between NPP and deepsea flux on large scales. The issue of spatial and temporal measurement scale is an important one to address in future work, and will make it possible to capitalize on these distributional relationships, for example, to link satellite models of NPP with upper-ocean and even deep-sea carbon fluxes.

Data Availability Statement
OFP data are available from the NSF Biological and Chemical Oceanography Data Management Office (BCO-DMO). Data and code are available at https://doi.org/10.5281/zenodo.4675075. For review, the global data can be found from the  and Mouw et al. (2016) sources cited in the text, and the time-series data can be found at https://github.com/bbcael/grl-deep-flux-preliminary-data.