Bursts Characterize Coagulation of Rods in a Quiescent Fluid

Under favorable conditions, microscopic phytoplankton cells dwelling in the oceans can divide rapidly and reach high concentrations, forming blooms that span kilometers and last for weeks. When blooms collapse, dead cells settle and aggregate into “marine snow” particles, resulting in a large and climatically important vertical flux of carbon from the ocean surface to its depth, a process known as the “biological pump.” To date, the formation of marine snow has been modeled as coagulation between spherical particles driven by gravitational settling and turbulent mixing, characterized by coagulation dynamics that converge onto time-independent concentrations of aggregates. However, many phytoplankton species are elongated and how their rodlike shape affects the aggregation process has remained unknown. Here, we study marine snow formation in a quiescent fluid assuming the constituent particles are elongated and form bundles upon encounter. We derive the collision kernel between dissimilar rods settling under gravity and discover that the most frequent collisions occur between the thinnest and thickest bundles, rather than between bundles of similar size. As a consequence, in the full coagulation model that combines exponential growth with settling, the thin-thick coupling can lead to statistically stationary states where the concentrations of aggregates of different size oscillate in time, exhibiting periodic bursts. The bursts are predicted to occur on the scale of a week and eventually lead to broadening of aggregate size spectra and may thus be highly relevant for plankton dynamics and the carbon cycle in the ocean.

Smoluchowski equations of coagulation theory describe aggregate formation between colliding particles [1,2]. Developed a century ago to quantify coagulation of colloids, they have since been applied to study the formation of rain droplets [3,4], gelation of polymers [5], large-scale structures in the Universe [6], and marine snow particle formation in the ocean [7]. Collision kernels are the central element of coagulation theory since they quantify the frequency of encounters. While many physical mechanisms can lead to encounters, including diffusion [1,2], gravitational settling [7], active swimming [8], and turbulence [3,4,7], the corresponding collision kernels have been typically computed only for spherical particles, neglecting the effect of different particle shapes.
Here, we study formation of marine snow by submillimetric elongated phytoplankton, which is key to understanding the ocean's carbon cycle [9]. Under favorable conditions, phytoplankton cells dwelling in the ocean can divide rapidly, forming blooms that can span tens of kilometers and last for weeks. When blooms collapse, dead cells settle and aggregate into particles termed marine snow, whose settling establishes a vertical flux of carbon from the ocean surface to its depth-the ocean's "biological pump" [10]. Models to date have considered the constituent particles of marine snow to be spherical and have identified gravitational settling and turbulence as the main drivers of the aggregation process [7]. However, many phytoplankton species are elongated, often with large aspect ratios [11]. For example, filaments of Skeletonema costatum or Trichodesmium are typically several micrometers wide and hundreds of micrometers long. Their settling speed reaches up to a hundred micrometers per second, keeping the Reynolds number small (< 0.1) and the Péclet number large (> 10 4 ). Additionally, upon encounters in the water column, filaments of Trichodesmium can form bundles of filaments arranged in parallel, called "tufts" [12], essentially realizing a collision model in which thin rods form thicker ones. During blooms, the concentration of individual cells reaches thousands of cells per milliliter [13], but the relatively small volume fractions (< 10 −5 ) justify coagulation models based on binary collisions. Yet, the effect of the elongation of individual filaments and their tendency to form bundles on the magnitude and dynamics of marine snow formation has remained unknown.
We have recently derived the collision kernel for identical, randomly oriented rods settling under gravity [14]. Here, we first extend these results to the case of dissimilar rods, which can differ between each other by density, length, or aspect ratio. Motivated by marine snow formation between elongated phytoplantkon filaments, we then apply the kernel to study coagulation of elongated particles in a quiescent fluid. Assuming filaments form bundles upon encounter, we discover a strong coupling between the thinnest and thickest bundles, rather than between bundles of similar size. By contrast, models based on spherical particles predict strong coupling between aggregates of similar size. In the full Smoluchowski model for coagulation that combines exponential growth of cells with settling, this thin-thick coupling leads to statistically stationary states where the concentrations of aggregates of different size oscillate in time [15,16], characterized by periodic bursts occurring on the scale of a week. The presence of bursts in the coagulation of rods is in stark contrast to time-independent coagulation predicted by models based on spherical particles.
In a quiescent fluid, the settling velocity of an elongated ellipsoid ("rod" or "filament") pointing along the unit vector p is [17] uðpÞ ¼ 1 8πμl where I is the identity matrix, l is the rod length, μ is the fluid's dynamic viscosity, and the gravity force is f ¼ ΔρV rod g, where V rod is the volume of the rod, g is the gravitational acceleration, and Δρ is the density difference between the rod and the fluid. The dimensionless quantities β 0 and β 1 are functions of the rod's aspect ratio A, given by (A > 1) [17], where ξ AE ¼ ln½A AE ðA 2 − 1Þ 1=2 =ðA 2 − 1Þ 3=2 . Equation (1) predicts that a rod's settling velocity varies with its orientation. For an ensemble of identical, randomly distributed and randomly oriented rods each settling according to Eq. (1), the collision kernel reads [14] Γ identical rods ¼ lβ 1 ΔρV rod g=ð16AμÞ: Below, we generalize Eq. (3) to the case of two ensembles of dissimilar rods, which we refer to as type I and II, characterized by the two triplets of parameters ðΔρ 1 ; l 1 ; A 1 Þ and ðΔρ 2 ; l 2 ; A 2 Þ, respectively. We expand the kernel as where Γ 0 rods is the "perfect-rod" limit of thin rods, Γ 1 rods is the finite-width contribution of the order OðA −1 Þ, etc. Neglecting terms of order OðA −2 Þ, we next find a closed-form formula for Γ 0 rods and provide an integral formula for Γ 1 rods that can readily be computed numerically. To find Γ rods ≈ Γ 0 rods þ Γ 1 rods , we first compute the collision kernel Γ p 1 p 2 between two subensembles: subensemble I of rods of type I oriented parallel to arbitrary direction p 1 and subensemble II of rods of type II oriented parallel to arbitrary direction p 2 [ Fig. 1(a)]. Averaging Γ p 1 p 2 over all orientations p 1 and p 2 yields Γ rods . We focus on the subensemble I and a single "target rod" from the Coagulation model We find the collision kernel Γ rods between two ensembles of randomly distributed and randomly oriented dissimilar rods of type I and II settling in a quiescent fluid. The two types can differ by density offset, length, or aspect ratio. (a) We first compute the collision kernel Γ p 1 p 2 between two subensembles of rods of the two types oriented parallel to two arbitrary directions p 1 and p 2 ; averaging Γ p 1 p 2 over p 1 and p 2 yields Γ rods . Shown is the subensemble I (blue) oriented along p 1 and a single target rod (red) from the subensemble II oriented along p 2 . The relative velocity between the type-I subensemble and the target is Δu 12 and their collision cross section is determined by the projections of the rods onto the collision plane (purple plane) perpendicular to Δu 12 . The green plane is spanned by Δu 12 and which is the area swept by the parallelogram (blue broken line) formed by the projections of the rods onto the collision plane. (c) For rods of finite width, approximated as cylinders, the collision cross section has an octagonal shape (blue broken line). (d) We use Γ rods to study coagulation dynamics of rods in a quiescent fluid and contrast it with coagulation of equal-volume spheres under differential settling and turbulence. For rods, aggregate i ¼ 1 is a filament of length l and aspect ratio A. Aggregate i > 1 has the same length l i ¼ l but smaller aspect ratio, representing filaments that form bundles upon collision. For spheres, aggregate i ¼ 1 is a sphere with radius r, aggregate i > 1 is a sphere with radius r i ¼ ri 1=3 .
subensemble II [blue and red rods in Fig. 1(a), respectively]. The relative velocity between the type-I subensemble and the target rod is Δu 12 and Γ p 1 p 2 is then the area swept out by the corresponding collision cross section. For perfect rods, the collision cross section is the parallelogram formed by the projections of p 1 and p 2 onto the plane perpendicular to Δu 12 [ Fig. 1(b)]. Thus, Γ 0 p 1 p 2 is given by the triple product Averaging over p 1 and p 2 gives [18] where Γ identical rods is computed using Eq. (3) with parameters ðΔρ 1 ; l 1 ; A 1 Þ and the density and shape mismatch factors For rods of finite width, the collision cross section has an octagonal shape [ Fig. 1(c)] and Γ 1 rods must be evaluated numerically [18].
Motivated by marine snow formation by elongated phytoplankton, we insert the collision kernel Γ rods in the Smoluchowski coagulation equation to study aggregation of filaments in a quiescent fluid. We take an individual phytoplankton cell to be a filament of length l and aspect ratio A [ Fig. 1(d)]. Upon collisions, filaments form bundles and we denote by c i the concentration of i bundles, i.e., bundles containing i filaments. We assume that all filaments are parallel to each other in a bundle, the length of an i bundle is l i ¼ l, the same as the length of an individual filament, and the aspect ratio is A i ¼ Ai −1=2 . This choice represents a maximum-overlap collision model in which, upon encounter, filaments maximize overlap with each other-in particular, a collision between an i bundle and a j bundle results in a thick (i þ j) bundle of length l and aspect ratio Aði þ jÞ −1=2 . All bundles are taken to have the same density offset Δρ and are assumed to be randomly oriented. The bundle concentrations evolve according to the Smoluchowski equation, where Γ ij is the kernel Γ rods evaluated using Eq. (4) with the parameter triplets ðΔρ; l; A i Þ, ðΔρ; l; A j Þ, and volume (7) represents the mean-field evolution of bundle concentrations in a water column of depth z. Terms proportional to Γ ij represent the nonlinear coagulation process mediated by collisions between bundles. Additionally, individual filaments (c 1 ) grow at rate α and bundles are removed from the system at rate w i =z due to settling, where w i ¼ ðΔρV i gÞ=ð8πμlÞ½β 0 ðA i Þ þ β 1 ðA i Þ=3 is the z component of Eq. (1) averaged over p. We focus on model parameters representative of elongated phytoplankton, which have typical aspect ratios around A ¼ 5, with A ¼ 50 being not uncommon [11]. Some cyanobacteria, such as the bloom-forming Trichodesmium, live as individual filaments, up to 1 mm long and a few micrometers wide [19], with aspect ratios reaching A ¼ 100. Additionally, Trichodesmium forms bundles of filaments called tufts, essentially realizing the maximum-overlap collision rule introduced above. In the remainder, we set the water column depth at z ¼ 100 m, single filament length at l ¼ 200 μm, aspect ratio at A ¼ 50, and filament density offset at Δρ ¼ 0.02ρ water ¼ 20 kg m −3 , where ρ water ¼ 1000 kg m −3 .
The collision kernel for rods favors interactions between individual filaments and the thickest bundles in the system. This thin-thick coupling is strongest in the perfect-rod limit, as can bee seen from the dependence of Γ 0 rods [Eq. (6)] on Comparing the collision kernels for dissimilar rods settling in a quiescent fluid (a)-(c) and spherical particles colliding under differential settling and turbulence (d)-(f) reveals that rods interact most strongly when their width mismatch is largest, in contrast to spheres. (a)-(c) The collision kernel for rods Γ rods can be approximated by the sum of the contributions due to the perfect-rod limit (Γ 0 rods ) and to the rods' finite width (Γ 1 rods ); the sum (Γ 0 rods þ Γ 1 rods ) is our estimate of Γ rods . (d)-(f) Kernels for spherical particles encountering each other under differential settling (Γ ds ), turbulence (Γ turb ), and both of these mechanisms together (Γ ds þ Γ turb ). The aggregate size is represented by the number of individual cells (filaments or spheres) it contains; the volumes of elongated and spherical aggregates of size i are equal. The density offset Δρ is identical for all aggregates. Parameters are Δρ ¼ 20 the two bundle sizes ði; jÞ [ Fig. 2(a)]. In particular, Γ 0 rods is largest for i ≈ 1 and j ≈ N, which is a consequence of its dependence on the square of the width mismatch between the bundles [Eq. (6)]. The contribution due to the finitebundle width Γ 1 rods [ Fig. 2(b)] was computed numerically [18]. Γ 1 rods also favors thin-thick coupling, although to a lesser degree than Γ 0 rods [ Fig. S.2(a) of Supplemental Material [18] ]. The two contributions add up to the total kernel Γ rods [ Fig. 2(c)]. The largest bundles we consider contain N ¼ 2000 filaments and have effective aspect ratio A 2000 ¼ 1.12.
Before studying the full coagulation dynamics described by Eq. (7), we compare Γ rods with collision kernels for spherical particles, used classically to model marine snow formation by either differential settling or turbulence [7]. For differential settling, the collision kernel between spheres of radii r 1 and r 2 reads where we take the spheres to settle according to Stokes's law, u sphere ðrÞ ¼ 2Δρgr 2 =9μ. In turbulent flow of intensity ϵ, the collision kernel reads [3,20] where ν ¼ μ=ρ water is the kinematic viscosity of water. To compare (Γ rods ) with (Γ ds þ Γ turb ), we take an individual spherical phytoplankton cell to have the same volume and density offset as an individual filament described above, and consider a turbulence intensity ϵ ¼ 10 −6 W kg −1 , characteristic of the ocean's surface layer [21]. A spherical aggregate containing i cells has radius r i ¼ ri 1=3 ; when an i aggregate and a j aggregate collide, they form an (i þ j) aggregate of radius r iþj ¼ rði þ jÞ 1=3 .  7) for the two kernels, Γ rods and Γ ds þ Γ turb , reveals fundamental differences in marine snow formation between elongated and spherical phytoplankton. The most striking difference is the emergence of statistically stationary periodic bursts in coagulation by elongated cells [ Fig. 3(a)], compared to the time-independent concentration profiles for spherical cells [ Fig. 3(b)]. The statistically stationary stage is characterized by the balance between the growth of individual cells (α ¼ 1 day −1 ) and the coagulation and settling of all aggregates. The periodic bursts in coagulation of elongated cells are a direct consequence of the thin-thick pairing characterizing Γ rods : when sufficiently many thick bundles develop in the system (c i rises for large i), they momentarily decrease the concentration of individual filaments (c 1 ) since thick bundles are most efficient at collecting individual filaments; because all aggregates settle, depletion of c 1 eventually leads to depletion of thick bundles, creating oscillatory dynamics. This run-and-chase pattern is akin to limit cycle solutions in predator-prey dynamics, with thick aggregates playing the role of predator and thin aggregates that of prey [22]. In contrast, the coagulation dynamics of spherical aggregates tend to pair aggregates of similar size, leading to a local coagulation cascade, whereby similarly sized aggregates create larger ones without significant contribution from the interactions between largest and smallest aggregates. For the parameters we considered, rods, A=50 rods, A=100 spheres  7) with the collision kernels Γ rods (a) and Γ ds þ Γ turb (b), respectively, for the same parameters as in Fig. 2. The individual cells i ¼ 1 grow exponentially at rate α ¼ 1 day −1 ; the initial conditions are c 1 ð0Þ ¼ 1 cm −3 and c i>1 ð0Þ ¼ 0. (c) Increasing the aspect ratio A with all other parameters fixed shows that the amplitude of bursts in the stationary stage grows with A, whereas the oscillation period remains nearly constant (inset). (d) The particle size spectra are broader for rods, a consequence of the thin-thick coupling between aggregates. The inset represents the sedimentation flux αhc 1 i as a function of A normalized by the sedimentation flux for spherical particles.
the bursts occur on the scale of a week and may thus be relevant for plankton dynamics and the carbon cycle in the ocean.
Making the rods thinner increases the amplitude of bursts [ Fig. 3(c)], while the corresponding oscillation period remains approximately constant (inset). This rise in the oscillations' amplitude is a direct consequence of the increase of the relative importance of the perfect-rod kernel Γ 0 rods over the finite-width contribution Γ 1 rods , which leads to stronger thin-thick coupling. More broadly, oscillations in Smoluchowski-type equations have only recently been discovered [15,16]. Our results establish that elongated particles that form bundles are a concrete realization of such periodic states in the case when the new aggregates are provided by exponential growth of individual filaments, as opposed to a constant injection rate [15] or closed systems in which aggregates undergo disintegration [16]. Additionally, these oscillations lead to the broadening of particle size spectra [ Fig. 3(d)], resulting in states with relatively many large bundles, which may be a relevant factor contributing to the observed particle size spectra in the ocean [23,24]. Finally, since in the steady state the mean total sedimentation rate is equal to the phytoplankton growth rate αhc 1 i, we observe that the sedimentation flux for rods is larger than for spheres [inset in Fig. 3(d)], by about 30% for A ¼ 50 and about fivefold for the thinnest rods.
We specifically investigated how the fundamental characteristics of phytoplankton cells affect the oscillatory burst dynamics [18]. This analysis revealed that increasing the growth rate α and decreasing the density offset Δρ increases both the oscillation's amplitude and frequency (Fig. S.3 [18]). For example, doubling α and halving Δρ corresponding to ranges of variation easily expected to occur in nature increases the oscillation amplitude and frequency approximately by a factor of 2. Thus, rapidly dividing and nearly neutrally buoyant elongated phytoplankton cells are expected to exhibit stronger, more rapid burst dynamics.
Our coagulation model builds on several assumptions. First, we assume all bundles in the volume considered are randomly oriented. While well-separated rods at low Reynolds number maintain their orientation as they settle, thus preserving the initial random orientation, hydrodynamic interactions might preferentially orient bundles [25,26]. However, owing to the low volume fraction of cells in typical blooms conditions (< 10 −5 ), these ordering effects are quenched by the mixing induced by weak turbulence [18]. Second, we assume perfect collision efficiency, reflecting the fact that phytoplankton tend to become sticky upon dying, often as a result of the excretion of polymeric substances [27]. Nevertheless, additional effects, such as hydrodynamic interactions or elasticity [28], might affect the collision efficiency [4]. If p ij is the collision efficiency between i, j bundles, then the collision kernel in Eq. (7) must be replaced with p ij Γ ij . In the simplest case, p ij ¼ p, the concentrations are rescaled by a factor p −1 , without otherwise altering the dynamics. Third, since we focus on the dynamics of rods, the largest bundle size is chosen such that the largest bundle is still elongated, rather than oblate. Should the largest bundle size be even larger, bundle configuration with filaments in parallel to each other may no longer be optimal from the overlapmaximization standpoint, and the largest aggregates could assume a more isotropic arrangement, such as the "puff" colonies formed by Trichodesmium [12]. Extension of the model in this case would require calculating collision kernels between rods and spheres.
To conclude, we have derived the collision kernel for dissimilar rods settling in a quiescent fluid. We found that collisions between thin and thick rods are favored, which leads to markedly different coagulation dynamics than coagulation of spherical particles. For parameters relevant to marine snow formation by phytoplankton, we predict that concentrations of elongated phytoplankton exhibit periodic bursts on the scale of a week. Such oscillatory dynamics of elongated phytoplankton is in stark contrast to the dynamics of spherical cells, which converge on time-independent concentrations. Our findings indicate that elongation can drive marine snow formation at rates comparable to or higher than those predicted for spherical particles, and that it does so even under quiescent conditions. Finally, the bursts in coagulation of rods are predicted to occur on the scale of a week and eventually lead to broadening of aggregate size spectra and may thus be highly relevant for plankton dynamics and the carbon cycle in the ocean.