Reconciling the Size‐Dependence of Marine Particle Sinking Speed

Sinking particles are critical to the ocean's “biological pump,” sequestering carbon from the atmosphere. Particles' sinking speeds are a primary factor determining fluxes and subsequent ecological and climatic impacts. While size is a key determinant of particles' sinking speeds, observations suggest a variable size‐sinking relationship, affected by other particle properties, resulting in substantial spread in parameterizations of particle sinking and fluxes. We compile particle size‐sinking observations and apply hierarchical Bayesian statistical models to resolve the size‐sinking relationship while accounting for other factors. We find an overall scaling close to the general Navier‐Stokes drag equation, and differences between particle types, open ocean versus coastal/laboratory particles, and in situ versus ex situ methods. These results can help harmonize how Earth system models parameterize particle fluxes and support a weaker size‐dependence than often assumed, with implications for the flux contribution of small particles and the predicted future shrinking of marine particle populations.

where g is the gravitational constant and μ is the dynamic viscosity of seawater (kg/ms). This simple equation suggests that sinking speeds scale positively and nonlinearly with size and is practically appealing because of its simplicity and basic theoretical underpinnings. However, it has been noted that observations of sinking particles typically deviate widely from the basic relationship (Diercks & Asper, 1997;Laurenceau-Cornec et al., 2015 Figure 1). Individual studies of different assemblages of particles over the past few decades find size-sinking relationships ranging from a Stokes'-Law-like relationship (Azetsu-Scott & Johnson, 1992) to an overall decrease of sinking velocity with size (McDonnell & Buesseler, 2010). The diversity of assumed relationship demonstrates the importance of various factors that obscure the overall size-dependency of sinking speeds such as inclusion of ballast minerals (Klaas & Archer, 2002) or "sticky" transparent exopolymers (Engel & Schartau, 1999), ambient temperature (Iversen & Ploug, 2013), or particles' geometry (Johnson et al., 1996). Several variations to Equation 1 have been proposed for nonspherical and porous particles, generally with scaling exponents <2 (Tang & Raper, 2002), but it is not clear whether any of these is appropriate for the large diversity of natural marine particles.
As a result, Earth system models use a variety of different size-sinking relationships to parameterize particulate organic matter fluxes (e.g., Dunne et al., 2010;Dutkiewicz et al., 2015;Gehlen et al., 2006;Kriest & Oschlies, 2008;Yool et al., 2013; Figure 1). These size-sinking relationships yield order-of-magnitude differences in sinking speeds for particles of the same size ( Figure 1) and therefore predict widely different fluxes and introduce uncertainty into models of ocean biogeochemistry (Gehlen et al., 2006;Niemeyer et al., 2019). Furthermore, size-sinking relationships have additional implications for climate change projections, as plankton models suggest a shrinking of cell size with warming, which would positively feedback onto atmospheric CO 2 due to reduced particle sinking speeds predicted from size-sinking relationships (Bopp et al., 2005;Finkel et al., 2010). Though it has become increasingly clear that environmental context, particle type, and other factors can limit our ability to extrapolate observed empirical relationships to estimate large scale fluxes (Laurenceau-Cornec et al., 2015), these uncertainties have not yet been treated in a statistical fashion, which could help harmonize biological pump parameterizations in ocean models. Determining particles' size-sinking relationships is not only of interest for global modeling but also for the increasingly prevalent technique of estimating particle fluxes from optically derived ambient particle size distributions (Guidi et al., 2008) using, for example, the Underwater Vision Profiler (Picheral et al., 2010).  Table S2 and from models. Individual studies are randomly colored, with ordinary least squares linear regressions of log-transformed data overlaid. The solid orange line is the median result of the hierarchical regression to the overall data set (log 10 (α) = 2.11, β = 0.63) and the solid black line is the result of an ordinary least squares regression (log 10 (α) = 2.11, β = 0.45). Dashed, dotted, and dash-dotted lines are common model parametrizations of size-sinking relationships (Dunne et al., 2010;Dutkiewicz et al., 2015;Gehlen et al., 2006;Niemeyer et al., 2019;Yool et al., 2013). Thin black line is a Stokes'-Law-like β = 2. Right: 100-m normalized fluxes (flux at a given depth divided by flux at 100 m) for particle populations of different sizes with a remineralization rate of k = 0.24 days −1 (Collins et al., 2015) and sinking speeds from relationships in the left panel (β = 2 not shown).
These methods require knowledge of how particle sinking speed and mass scale with size, and how these might vary for different particles or environments.
Here, we attempt to resolve the size-sinking relationship in marine particles by assembling the largest data set of particle size and sinking speed to date. Our objective is to quantify size-sinking scaling relationships from particle observations, provide confidence bounds on this scaling, and its dependence on other factors. We develop a Bayesian hierarchical model framework, motivated from a generalized Navier-Stokes drag equation, to estimate a universal size-scaling relationship while taking into account context-dependent factors. For nonspherical, heterogeneous particles like "marine snow" aggregates or zooplankton fecal pellets, a more general Navier-Stokes drag equation can be used to describe sinking speed (Laurenceau-Cornec et al., 2019;Stokes et al., 1851): ( 2) where C D is the dimensionless drag coefficient. Here both Δρ and C D are size-dependent. The former dependence is typically expressed in terms of the fractal dimension f as Δρ ∝ d f−3 , whereas C D scales inversely with Reynolds number for solid spheres. (f can be thought of in less technical terms as quantitatively expressing how much less volume-filling larger particles are than smaller particles.) In general, f and C D are not well-characterized for complex natural particle shapes and structures (Loth, 2008;Tang & Raper, 2002 where β is the dimensionless scaling exponent and α is the prefactor coefficient (m/d), or the speed of a 1 mm particle. Viewed this way, α and β are unknown parameters with variability that depends on particle type and environmental factors. We formulate Equation 3 as a hierarchical Bayesian model and allow deviations from an overall size-scaling relationship according to: (1) aggregate particles versus. fecal pellets, coastal versus open ocean versus laboratory measurements, (3) in situ versus. indirect in situ versus flow chamber versus other measurement methodologies, (4) "warm" (≥14°C) versus "cold" (<14°C) environments, (5) diatom versus nondiatom for aggregate particles, and (6) ballast versus nonballast for aggregate particles (see Methods for more detail). Hierarchical Bayesian regression has enjoyed success in other scientific fields (Gelman & Hill, 2006), but has not previously been applied to marine particle datasets. We fit the model to the compiled data set, representing sinking speeds from a wide variety of contexts, and report posterior probability distributions for the overall scaling and context-specific deviations. We discuss the implications of our estimated size-sinking relationships for biogeochemical models and for the global carbon cycle.

Overall Scaling
The hierarchical model suggests a posterior median overall scaling exponent β = 0.63 (with a standard deviation ± 0.06) and a median coefficient α = 129 ± 26 m/d. The overall and study-specific posterior distributions for α and β are plotted in Figure 3. Despite large study-to-study variability in both parameters, the hierarchical model finds an overall scaling exponent much smaller than Equation 1 (β = 2) and roughly in line with Equation 2 (β = 1/2). The difference between the 1/2-scaling in Equation 2 and the β = 0.63 ± 0.06 scaling in Figure 3 may be due to environmental and/or methodological discrepancies (see below) but in general suggests that the ratio Δρ/C D changes only weakly with particle size. We note that the overall parameter distributions are virtually unchanged across all the different model networks we used, indicating that these overall α and β parameter values are robust (top level in Figure 2). The hierarchical model also finds a higher β than the β = 0.45 ± 0.01 found by an ordinary least squares (OLS) regression of the total aggregated data set (a posterior probability of p > 99%), but almost exactly the same α, suggesting faster sinking velocities than OLS for particles >1 mm and slower sinking velocities for particles <1 mm ( Figure 1).
The total flux F (g m −2 d −1 ) of an assemblage of particles can be expressed (Guidi et al., 2008) as is the average size-sinking speed relationship, and m(d) is the average size-mass relationship. Thus substituting Equation 3 into Equation 4 shows that total particle flux F is proportional to the coefficient α, and that the relative contribution of small and large particles to F will be determined by β and how both N and m scale with d. The utility of the hierarchical approach is primarily in constraining β, as seen by the fact that the overall α = 129 ± 26 m/d matches both the overall α found by a simple OLS regression and a common "rule of thumb" estimate that particles sink on the order of 100 m/d. An overall scaling exponent of β = 0.63 ± 0.06 suggests a much more substantial contribution of small particles to total flux than if particles sank like solid spheres or individual cells (Laws, 1975;Niemeyer et al., 2019), as has been increasingly recognized (Alonso- González et al., 2010;Durkin et al., 2015;Riley et al., 2012). This is further supported by the fact that a low β is consistent with a low fractal dimension f, that is, a low mass-size scaling exponent, even if the quantitative link between β and f is uncertain (Johnson et al., 1996;Tang & Raper, 2002). Furthermore, this suggests particle fragmentation is less effective at attenuating total flux; if a particle becomes fragmented, the change in average sinking speed of the resulting (smaller) particle fragments will be less than if β were larger. A β = 0.63 ± 0.06 also suggests that future decreases in particle size distributions (Bopp et al., 2005;Finkel et al., 2010) would not change the total flux F nearly as much as if β were more akin to the larger values for individual cells (Laws, 1975) or solid objects (Clift et al., 2005;Stokes et al., 1851).
The differences in relative contribution to total flux of small and large particles between our overall result and example model parameterizations can be seen in Figure 1. For a parameterization using β = 1.17, effectively all of the flux below the upper few hundred meters of the water column is due to very large (∼3 mm) particles, as smaller particles are attenuated rapidly. In contrast, a step-function parameterization vastly overpredicts flux due to ∼300 μm particles, and vastly underpredicts flux due to ∼30 μm particles, compared to the hierarchical model's estimated relationship. Our hierarchical β of 0.63 ± 0.06 is both globally plausible and biogeochemically relevant as supported by a recent study (Niemeyer et al., 2019) which showed that a β = 0.62 appreciably outperforms a β = 1.17 in a global biogeochemical model's ability to simulate tracer distributions. The value of 0.62 in that study was derived from one of the equations that have been proposed to link β to particles' fractal dimension f (Kriest, 2002).

Group-Level Scalings
Our hierarchical models' group-level results largely confirm general expectations about the deviations from the overall size-sinking relationship ( Figure 4; Table S1). First, both α and β are larger (p > 99% and p ≈ 86%, respectively) for fecal pellets than for aggregates. The higher α is expected as fecal pellets are generally thought to be denser (meaning larger Δρ) and less rough than aggregates (meaning smaller C D ; Ebersbach et al., 2011;Yoon et al., 2001). The higher β is also partly explained by lower porosity, or higher fractal dimension, of the fecal pellets, as the effective density of fecal pellets should decrease comparatively less with CAEL ET AL.
10.1029/2020GL091771 4 of 11 Individual studies' scaling coefficients are drawn from groups (i.e., dividing studies by type of particle, type of environment, type of method, or temperature -or no grouping is imposed in the three-level network case) of studies' distributions, which in turn are drawn from the overall scaling coefficients' distribution. Here p refers to a probability density function, i is an index indicating different groups, j is an index indicating different studies within groups, and k is an index indicating different individual particle observations within studies, that is, w ijk is the kth sinking speed observation within study j within group i. The likelihood of each individual observation given the parameters' values is evaluated on the "data" level. See Materials and Methods for details.
increasing particle size, resulting in a higher β. In any case these results provide supporting evidence that aggregates' and fecal pellets' sinking should be parameterized differently in biogeochemical models (Table S1).
We find larger α in warm environments (p ≈ 76%), expected due to the lower viscosity of seawater at higher temperatures (Taucher et al., 2014) resulting in faster sinking speeds for all particles. α is also higher for ballasted (p ≈ 94%) and diatom-enriched (p ≈ 75%) aggregates, as expected from the higher densities of these materials (Engel et al., 2009). Diatom-enriched and ballasted aggregates are not only expected to be composed of more dense material but also more compact (Laurenceau-Cornec et al., 2019), and therefore expected to have a higher β, also borne out by these analyses (p ≈ 75% and p ≈ 69% respectively). It is unclear whether β would be expected to be higher, lower, or the same for particles from warmer environments; β is found to be larger for the warm group (p ≈ 72%), further indicating faster sinking speeds for these particles.
We find larger α (p ≈ 91%) and β (p ≈ 77%) for coastal particles compared to open ocean particles. Particles from coastal environments are expected to sink faster because coastal environments are subject to loading of minerals and other dense material from rivers, sediment resuspension, dust deposition, and other sources, and also because coastal ecosystems tend to be more dominated by diatoms (Armbrust, 2009). Similarly, as laboratory-generated particles are typically generated in volumes with artificially high collision rates (Jackson, 2015), from diatom cultures or from phytoplankton collected in close proximity to the laboratory, they should also have a larger α and β for the same reasons as above, which we also find (p ≈ 85% and p ≈ 83%, respectively). Though these differences are unsurprising, they are important because much of our understanding and parameterization of open ocean particles comes from coastal and laboratory studies, which may influence not only the magnitude of sinking speeds (α) but also their size dependence (β). We find a β of 0.56 ± 0.09 for the open ocean group, that is, within one standard deviation of the 1/2 scaling from Equation 2. When modeling sinking speeds of open ocean particles, we therefore recommend the β = 0.56 ± 0.09 exponent. As above, we interpret the lower open ocean β to result from higher porosity or lower fractal dimension of open-ocean particles.
Finally, surprising and potentially illuminating differences are revealed by the hierarchical model which groups by methodology. The in situ direct methods (those which measure sinking speeds, e.g., Carder et al., 1982) and indirect methods (those which infer speeds from the ratio of flux to concentration, Mc-Donnell & Buesseler, 2010) are in very good agreement in terms of α and β. The indirect group's parameters have larger uncertainties, likely due to smaller number of studies using these methods, but their central estimates are indistinguishable. In contrast, the "flow chamber" (Ploug & Jørgensen, 1999) and "other ex situ methods" (e.g., roller tank (Shanks & Edmondson, 1989) or settling column (Smayda, 1969)) groups have substantially different posterior distributions for both parameters. Both have larger α values than the "in situ" group (p ≈ 98% and p > 99%, respectively), but also larger β values (p ≈ 93% and p ≈ 75%, respectively).
As with the open ocean group above, the in situ direct and indirect methods groups are consistent with the 1/2 scaling from Equation 2. While the differences between methods distributions may be in part explained by the types of particles considered in the studies employing these methods, these ex situ methods have been employed widely and in a variety of contexts, and the evidence for these differences is strong. Fragile particles that would be disturbed by collection and therefore not easily captured by ex situ techniques (Giering et al., 2020) should be expected to sink slower and be less compact (i.e., have lower fractal dimensions), accounting for these parameter differences, but these results do suggest that such particles are a non-negligible component of in situ sinking particle populations.  Figure 2) with medians and standard deviations given in the legends; colored lines are the posteriors for the individual studies ("study" level in Figure 2).

Discussion
These results help resolve the size-sinking relationship of marine particles, with implications for particulate carbon fluxes and our understanding of the global carbon-climate system. Our approach establishes a baseline for a universal size-sinking relationship and quantifies important context dependent deviations. The hierarchical framework allowed us to account for unbalanced sampling across different contexts and provides a rigorous uncertainty bound that can also be propagated through subsequent carbon flux analyses. With this new hierarchical analysis across datasets, we have provided evidence that size-sinking relationships hold at broad scales and under various environmental conditions, and have quantified them for use in biogeochemical modeling and optical flux estimation.
These parameterizations can inform not only size-spectrum-resolving biogeochemical models (Kriest & Oschlies, 2008) but also those with only a few (e.g., two) size classes (Yool et al., 2013), as well as those with implicit flux profiles (Cael & Bisson, 2018;Najjar et al., 2007). The parameterizations are also important for the rapidly expanding field of optically derived particle flux estimation (Giering et al., 2020). It also has important implications for particle fluid dynamics that the study-specific scaling exponents vary greatly from study to study while the central scaling exponent is close to that of the generalized Navier-Stokes' drag equation exponent of 1/2 (Equation 2). This suggests that the size-dependencies of excess densities and drag coefficients of natural aggregates vary between particle populations, but that overall they are in approximate balance.
To date, carbon flux models generally combine sinking speed and remineralization relationships to achieve a predicted distribution of biogeochemical tracers that approximates observations (Najjar et al., 2007). However, the variety of size-sinking relationships in use (Figure 1) suggests that the two processes are being combined in highly uncertain ways. Uncertainty in size-sinking and remineralization relationships have direct implications for climate change (Völker et al., 2016). Warmer oceans are predicted to enhance water column stratification and reduce nutrient fluxes, which is predicted to shift cells toward smaller sizes (Finkel et al., 2010). Shifts in average cell size would lead to decreases in sinking speed and therefore carbon flux, which feedbacks onto atmospheric CO 2 and global climate. Our results suggest a comparatively weak size-dependence of sinking speed, which suggests this climate feedback may be weaker than previously assumed. A better parameterization of the size-sinking relationships should also allow better isolation in Earth System Models of remineralization relationships that are expected to vary with warming due to faster rates of heterotrophic metabolism (Cael & Follows, 2016;Cavan & Boyd, 2018;Cavan et al., 2019), especially given the large vertical gradients in temperature in the upper ocean. Improved parameterization of particle flux thereby provides a means to better predict changes in the ocean carbon cycle with climate change.
Our results underscore the importance of not only size but also composition, structure, and/or density in determining particle sinking speed; as seen in Figure 1, for example, a small aggregate can sink at similar velocity as that of a large aggregate. While much remains unknown or uncertain about the controlling factors for sinking speeds of in situ formed aggregates, comparatively much more is known about the plankton communities from which these aggregates are formed, including where and when different types of plankton are prevalent or dominant (Barton et al., 2013). Therefore, the distinct parameterizations we provide for fecal pellets, ballasted aggregates, and nonballasted aggregates (Figure 4, Table S1) that extract information CAEL ET AL.   Figure 2). Vertical lines show the medians of each distribution. n. b. the bottom two rows do not include fecal pellet data, and legends from the left column apply to the right column. See Table S1 for medians and confidence intervals. about the size-sinking relationships for different particle/aggregate types may be particularly useful for improving particle flux estimates based on community composition. This is an especially appealing possibility as it is becoming possible to differentiate community structure information (e.g., whether diatoms vs. coccolithophores vs. flagellates are dominant) from remote sensing (Kramer & Siegel, 2019), though challenges persist for the rigorous application of such approaches (Cael et al., 2020). Given suitable data, it would be particularly valuable in this context to estimate scaling relationships for aggregates formed from plankton communities dominated by different functional groups such as diatoms versus coccolithophores.
While our results help bridge the gap between modeled fluxes and observed sinking speed variability, more work remains to be done in understanding the unexplained variance (Figure 1). The statistical framework herein is useful for propagating the uncertainty due to explained variance in the data. Some of the unexplained variance is likely due to spatial sampling variability and measurement error, but other more mechanistic relationships should also be developed to explain the variance that may provide more predictive relationships for modeling sinking particle fluxes. Note also that the relationships we find are bulk scalings that are of interest for predicting particles' impact on biogeochemical fluxes; individual particles' sinking speeds will undoubtedly vary around these scalings.
In summary, our results advance our understanding of ocean carbon fluxes through an improved, data-driven description of the gravitationally driven sinking of marine particles. Our data compilation and hierarchical modeling suggests an approximately 0.6-scaling of sinking speed with particle size, with implications for the role of small particles in the biological pump and for how carbon flux will vary with climate change. Despite these results, the data compilation also points to significant unexplained variability that requires further study to confidently predict changes in these climatically important ocean processes.

Size-Sinking Data Compilation
We compiled observations of 5,655 particles' sizes and sinking speeds from 54 studies published in 32 papers (Alldredge & Gotschalk, 1988Azetsu-Scott & Johnson, 1992;Belcher et al., 2016aBelcher et al., , 2016bCarder et al., 1982;Deibel, 1990;Diercks & Asper, 1997;Engel & Schartau, 1999;Gibbs, 1985;Engel et al., 2009;Hawley, 1982;Hill et al., 1998;Iversen & Robert, 2015;Iversen et al., , 2017Jouandet et al., 2011;Kajihara, 1971;Laurenceau-Cornec et al., 2015McDonnell & Buesseler, 2010;Nowald et al., 2009;Shanks & Trent, 1980;Smayda, 1969;Syvitski et al., 1995;Van der Jagt et al., 2018); see Table S2. Size data were standardized to equivalent spherical diameter in units of millimeters. Sinking data were standardized to sinking speed in units of meters per day. Data that we were otherwise not able to access were digitized (https://automeris.io/WebPlotDigitizer/). Previous analyses of digitized data suggest that the error introduced by digitization is negligible (<1%; Cael & Bisson, 2018). Particles were classified as aggregate or fecal pellet, diatom-enriched or not, ballasted or not, and by their environment according to information from the publications. Temperatures were classified as "warm" if ≥14°C and "cold" otherwise. We chose this coarse temperature classification because many of the publications did not report environmental temperature, in which cases we supplemented by referencing the sampling locations and times to a monthly climatology (Reynolds et al., 2007); 14°C is the median of that climatology. Methods were classified into in situ direct observations, indirect in situ observations (i.e., estimation by the ratio of flux to concentration, e.g., Jouandet et al., 2011), flow chamber observations (Ploug & Jørgensen, 1999), and other ex situ methods (which measure speed by distance traveled per unit time, e.g., settling column (Smayda, 1969)). See Table S2 for the sources and classifications of the data.

Bayesian Hierarchical Analysis
We consider the model that individual observations vary according to study-specific size-sinking scaling relationships. These study-specific scalings then deviate from an overall scaling relationship for marine sinking particles. We also consider models where the study-specific scalings deviate from group-specific scaling relationships, which in turn deviate from an overall scaling relationship. We formalize each of these by fitting a hierarchical model to these data, which captures group, within-group, and within-study variance in a single statistical framework. Each hierarchical model can be represented as a network of conditional probabilities (see Figure 2). At the "overall" level we assume that the overall particle size-sinking relationship can be represented as a scaling relationship (Equation 3) governed by overall parameters  10 log and  . Group level parameters (indexed by i) then deviate from these overall distributions according to the group parameter distributions We assume broad uniform priors over all α, β, and σ parameters. Broad bounds were imposed to keep the parameters from nonphysical values. No posteriors were close to the bounds and so these bounds did not impact the inference.
Models are fit via Markov Chain Monte Carlo (MCMC) sampling, drawing random samples of model parameters and propagating them through the hierarchical network. Doing so many times yields a set of samples from the posterior distribution, that is, the distribution of the model parameters after the data are taken into account. We use the RStan and brms software package to perform the MCMC sampling Bürkner et al., 2017;Carpenter et al., 2017). We run 4 chains with 5,000 iterations per chain. We discarded the first 2,500 samples from each chain as warmup, leaving a total of 10,000 samples from each posterior. All R statistics for the MCMC chains were between 1.0 and 1.1 which indicated convergent mixing.

Data Availability Statement
All code and data will be made shared on GitHub and archived on Zenodo, and the data will also be archived on Pangaea, should this manuscript be accepted for publication. The data are also available via the supporting information. Council (NE/R015953/1) and the Horizon 2020 Framework Program (820989). BBC has also received funding from the European Union's Horizon 2020 Research and Innovation Program under grant agreement No. 820989 (project COMFORT, our common future ocean in the Earth system-quantifying coupled cycles of carbon, oxygen, and nutrients for determining and achieving safe operating spaces with respect to tipping points). The work reflects only the authors' view; the European Commission and their executive agency are not responsible for any use that may be made of the information the work contains.