2. Qualitative collision mechanisms

The nature of encounters between marine snow and suspended objects ranges from purely diffusive to purely ballistic, because marine snow particles are highly heterogeneous (Trudnowska et al. Reference Trudnowska, Lacour, Ardyna, Rogge, Irisson, Waite, Babin and Stemmann2021), and cover a wide range of sizes and sinking speeds. The size of marine snow particles varies in the range of $1\;{\unicode{x03BC}}$ m to several millimetres (McDonnell & Buesseler Reference McDonnell and Buesseler2010; Bochdansky, Clouse & Herndl Reference Bochdansky, Clouse and Herndl2016; Williams & Giering Reference Williams and Giering2022), as visible in figure 1(a), with sedimentation speeds from zero to several hundreds of metres per day (Williams & Giering Reference Williams and Giering2022). For most particles, the Reynolds number $\textit {Re}$ is less than unity (Kiørboe et al. Reference Kiørboe2001; Alldredge Reference Alldredge1998). Thus, the Stokes approximation provides a suitable starting point, as confirmed by recent measurements (Chajwa et al. Reference Chajwa, Flaum, Bidle, Van Mooy and Prakash2024) of flow fields around sinking marine snow, shown in figure 1(b).

Figure 1.

The landscape of collision types in marine snow indicates the different physical encounter mechanisms at play. (a) Examples of different sizes and shapes of marine aggregates, imaged in situ off the coast of East Greenland. Image courtesy of E. Trudnowska, Polish Academy of Sciences. (b) Sample image of a marine snow particle collected at 80 m below sea level, with the flow field visualised by plastic microbeads. Image by R. Chajwa et al., CC BY 4.0 (Chajwa et al. Reference Chajwa, Flaum, Bidle, Mooy and Prakash2023). (c) Archetypal collision types between different objects (symbols) in the parameter space of the Péclet number $\textit {Pe}$ and relative size $\beta$ , based on experimental data in table 1. Possible collision types cover the whole space, ranging qualitatively from purely diffusive encounters, through advective-diffusive encounters, to direct (ballistic) interception. Existing collision models account for the limiting cases only.

Table 1.

Selected representative examples of particulate matter involved in marine snow encounters, with their typical size and sinking speed as reported in experimental observations and calculated the diffusion coefficient using the Stokes–Einstein relationship (assuming viscosity of seawater to be $\mu =1.6\times 10^{-3}\;\textrm {Pa}\,\textrm {s}$ ). When calculating $\beta$ and $\textit {Pe}$ for a pair, we assumed that the interaction range $b$ is the radius of a smaller particle and that the diffusion coefficient is a sum of individual diffusion coefficients.

We qualitatively describe a collision of two particles moving with relative velocity $U$ and diffusion constant $D$ considering two dimensionless numbers. With $a$ denoting the effective radius of the larger particle and $b$ the interaction range between the smaller and the larger particle (such that whenever the centres of the particles are at most $a+b$ apart they collide), we define the Péclet number $\textit {Pe}$ and $\beta$ describing the size ratio of the colliders as

(2.1) \begin{equation} {\textit {Pe}} = \frac {U (a+b)}{D}, \qquad \beta = \frac {b}{a+b}. \end{equation}

To estimate the possible values of the two parameters, we consider four illustrative actors: a small marine snow particle (smallest measured particles (McCave Reference McCave1984)); a medium-sized, mucus-laden particle (Chajwa et al. Reference Chajwa, Flaum, Bidle, Van Mooy and Prakash2024); a large particle (around the 10th particle mass quantile (Iversen et al. Reference Iversen, Nowald, Ploug, Jackson and Fischer2010)); and a non-motile bacterium (Kiørboe et al. Reference Kiørboe, Grossart, Ploug and Tang2002). This set covers a broad range of possible collision parameters, outlined in table 1. An estimate of $\textit {Pe}$ and $\beta$ for collisions between them is shown in figure 1(c), confirming that both $\beta$ and $\textit {Pe}$ span several orders of magnitude.

For very small values of $\textit {Pe}$ , diffusion dominates the encounter rate. In this scenario, collisions can be viewed as stochastic, where the large, slowly sedimenting particle ‘bumps into’ smaller particles. We refer to this collision mode as purely diffusive. When the colliders’ size ratio, $\beta$ , is vanishingly small but $\textit {Pe}$ becomes significant (e.g. $10^6$ ), advection can transport new particles into the depleted region. In this way, the sedimenting particle ‘bumps into’ more particles, although the encounter mechanism remains diffusive as there is no slip at the particle surface. We refer to this collision mode as advection-diffusion. On the right-hand side of figure 1(c), for very large values of $\textit {Pe}$ (e.g. $10^9$ ) and a non-negligible size ratio, $\beta$ (e.g. $10^{-2}$ ), previous studies have assumed that the diffusion of the smaller particle becomes negligible. In such cases, collisions can be conceptualised as a large, rapidly sedimenting particle ‘sweeping away’ stationary smaller particles. We refer to this collision mode as direct interception.

A natural question arises concerning the range of applicability of these models. This question escapes a simplistic answer, because only the limiting cases have been analysed so far. Without a quantitative model of the encounter rate, valid across a broad range of $\textit {Pe}$ and $\beta$ , assuming one limiting case over another may severely underestimate the encounter rates. We develop tools to resolve this tension in the remainder of this article.

3. Governing equations

Consider a sphere of radius $a$ sedimenting with a velocity $U$ in a quiescent fluid. In the Stokesian regime and for an incompressible fluid, the velocity field in cylindrical coordinates $(\rho ,\theta ,z)$ , with $\boldsymbol{U}$ aligned with the $z$ axis, is given by $\boldsymbol{u} = (u_\rho ,u_z)$ . The velocity components and the corresponding streamfunction $\psi$ expressed in the sphere frame of reference, are given by (Landau & Lifshitz Reference Landau and Lifshitz1987)

(3.1) \begin{align} u_\rho &= \frac {3a \rho z U}{4 R} \left ( \left ( \frac {a}{R^2} \right )^2 - \frac {1}{R^2} \right ), \end{align}

(3.2) \begin{align} u_z &= U + \frac {3aU}{4 R} \left ( \frac {2 a^2 + 3 \rho ^2}{3 R^2} -\left (\frac {\rho a}{R^2}\right )^2 - 2 \right ), \end{align}

(3.3) \begin{align} \psi &= \frac {1}{2}\, U\, \rho ^2\, \left ( 1 - \frac {3}{2} \frac {a}{R} + \frac {1}{2} \left ( \frac {a}{R} \right )^3\; \right ), \end{align}

(3.4) \begin{align} R^2 &= \rho ^2 + z^2. \end{align}

Figure 2.

Sedimenting Stokesian sphere colliding with Brownian objects with a non-zero interaction range. (a) Geometry of the collisions. Stokes flow streamlines around a particle of radius $a$ representing sedimenting marine snow, with an interaction range $b$ marked in green. The interaction radius accounts for the finite size of the suspended objects. Given the axial symmetry, we introduce sideways ( $\rho$ ) and downstream ( $z$ ) coordinates to parametrise the system. (b) Numerical solution of the advection-diffusion equation obtained using the finite element method for the steady concentration field around a sedimenting sphere (dashed area) with a size ratio $\beta = 0.2$ (dashed line) at ${\textit {Pe}} = 500$ . (c) Stochastic trajectories of objects obeying the advection-diffusion equation at ${\textit {Pe}} = 500$ . Here, $N = 3\times 10^{3}$ trajectories of objects were initially distributed uniformly on a large disk upstream of the sphere (thus more numerous further from the axis of symmetry). Trajectories that collided with the enlarged sphere (dashed) are terminated, leading to the formation of a characteristic wake free of objects behind the sphere. Note: the decreased number of trajectories upstream, at the centre, is caused by accumulating simulations in three dimensions to $\rho$ and $z$ . Thus a spatially uniform distribution transforms to a linearly increasing one, $2\pi \rho \,{\mathrm{d}}{} \rho$ .

The streamlines of the flow field around the particle are visualised in figure 2(a). Suppose now that the suspending fluid contains small objects of effective radius $b$ which undergo diffusion relative to the large particle, with a diffusion coefficient $D$ . Whenever the distance between the centre of a small object and the centre of the sphere is smaller than $a+b$ , the objects stick to the sphere and are removed from the surrounding liquid. This process can model direct contact of small spherical objects with a larger particle, but different capture mechanics can also be modelled in this fashion; one possibility is an electrostatic attraction between the particle and a small object with an effective interaction range $b$ (which can be derived by comparing the interaction potential with the typical energy of thermal fluctuations, $k_B T$ ) (Zaccone et al. Reference Zaccone, Gentili, Wu and Morbidelli2010). Alternatively, in the case of slip or mixed boundary conditions on the sedimenting sphere, an effective Stokes radius $a$ can be introduced and the difference between the true size and the Stokes size would give the effective interaction size. In these cases, $\beta = 0.2$ would cover, for example, the case where the large collider has the Stokes radius $a = 100 \;{\unicode{x03BC}} \textrm {m}$ , and the effective collision radius of $115 \;{\unicode{x03BC}} \textrm {m}$ (e.g. due to its porous nature), and it collides with a solid sphere of size $5 \;{\unicode{x03BC}} \textrm {m}$ giving $b=20 \;{\unicode{x03BC}} \textrm {m}$ , while at the same time introducing only a small perturbation to the flow field. Regardless of the specific origin of such interaction range, we focus here on the consequences of a finite interaction range described by $b$ . We further assume that the captured particles are much smaller than the sedimenting sphere, so that their influence on the flow field around the large sphere can be neglected. In this case, the steady state concentration profile of small particles $\varphi$ is governed by the advection-diffusion equation (Kiørboe & Thygesen Reference Kiørboe2001)

(3.5) \begin{equation} 0 = D {\nabla} ^2 \varphi - \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\varphi, \end{equation}

with the boundary condition of constant concentration, $\varphi = \varphi _0$ , upstream of the ball, and $\varphi = 0$ on a sphere with an effective radius of $a+b$ . To make (3.5) dimensionless, we choose the time scale to be $(a+b) / U$ , use $a+b$ as the length scale and $\varphi _0$ as the concentration scale and arrive at the dimensionless form

(3.6) \begin{equation} 0 = {\nabla} ^2 \varphi - {\textit {Pe}}\,(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\varphi ), \end{equation}

with the boundary conditions for concentration

(3.7) \begin{align} &\varphi (r = 1) = 0,\notag\\ &\varphi (r \to \infty ) = 1. \end{align}

The dimensionless velocity field, (3.1)–(3.2), can be written in terms of the rescaled large particle radius $\alpha = a/(a+b)$ as

(3.8) \begin{align} u_\rho &= \frac {3 \alpha \rho z}{4 R} \left ( \left ( \frac {\alpha }{R^2} \right )^2 - \frac {1}{R^2} \right ), \notag\\ u_z &= 1 + \frac {3\alpha }{4 R} \left ( \frac {2 \alpha ^2 + 3 \rho ^2}{3 R^2} -\left (\frac {\rho \alpha }{R^2}\right )^2 - 2 \right ). \end{align}

Crucially, the boundary conditions of no slip and no concentration are imposed here on different surfaces, with $\boldsymbol{u}(r=\alpha )=0$ , rather than the classical case, where $\boldsymbol{u}(r=1)=0$ .

We now calculate the mass of small particles intercepted by the large particle per unit time. Expressing (3.6) as $0 = \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{J}$ , we can write the particle current $\boldsymbol{J}$ as

(3.9) \begin{equation} \boldsymbol{J} = {\textit {Pe}} \, \boldsymbol{u} \varphi - \boldsymbol{\nabla } \varphi , \end{equation}

where we used the fact that the flow is incompressible, $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0$ . The total, dimensionless flux of particles onto the sphere, $\varPhi ^*$ , is given by

(3.10) \begin{equation} \varPhi ^* = \oint _{\textit{sphere}} \boldsymbol{J} \boldsymbol{\cdot }{\mathrm{d}} \boldsymbol{S}. \end{equation}

Using the Stokes theorem, and given the lack of source terms in (3.6), the integration surface in (3.10) can be changed to any other surface enclosing the sphere. When the integration surface is taken as a large cylinder that extends far from the sphere, we have $\partial _z \varphi \ll \boldsymbol{u} \varphi$ , and we can avoid calculating the gradient in the numerical evaluation of $\boldsymbol{J}$ . The dimensional flux is calculated as $\varPhi = \varphi _0 D (a+b)\varPhi ^*$ . With precise definitions, we are ready to discuss different approximations used to compute $\varPhi$ for given values of $\textit {Pe}$ and $\beta$ .

4. Divergent asymptotics of the limiting cases

Before solving (3.6) in the most general case ( ${\textit {Pe}}\geqslant 0$ , $\beta \geqslant 0$ ), we first briefly discuss its important limiting cases, namely the widely used direct interception limit and the advection-diffusion description with a zero interaction range. We also highlight important scenarios under which these limiting cases break down and yield asymptotically divergent predictions.

In the absence of flow ( ${\textit {Pe}} \to 0$ , purely diffusive encounters), (3.6) becomes a spherically symmetric Laplace’s equation, giving the classical expression for the diffusive flux $\varPhi _{{D}}$ onto a sphere (Karp-Boss et al. Reference Karp-Boss, Boss and Jumars1996)

(4.1) \begin{equation} \varPhi _{{D}} = 4 \pi D (a+b)\varphi _0. \end{equation}

When flow is present ( ${\textit {Pe}}\gt 0$ , advective-diffusive encounters) it is customary to normalise $\varPhi$ by this limiting case to obtain the dimensionless quantity proportional to the encounter rate. This ratio defines the Sherwood number (Karp-Boss et al. Reference Karp-Boss, Boss and Jumars1996), given by

(4.2) \begin{equation} {\textit {Sh}} =\frac {\varPhi }{\varPhi _{{D}}}. \end{equation}

Note that there is no consensus on the definition of $\textit {Sh}$ and $\textit {Pe}$ in the literature and the definitions might differ by a factor of two.

For an arbitrary value of $\textit {Pe}$ and in the limit of a zero interaction range ( $\beta \to 0$ ), (3.6) also describes heat transfer to or from a sphere in a slowly flowing fluid. In the limit of small or high Péclet number, analytical solutions of this problem were obtained using perturbative methods (Acrivos & Taylor Reference Acrivos and Taylor1962; Acrivos & Goddard Reference Acrivos and Goddard1965; Rimmer Reference Rimmer1968; Gupalo & Ryazantsev Reference Gupalo and Ryazantsev1972; Brunn Reference Brunn1982; Bell et al. Reference Bell, Byrne, Whiteley and Waters2013). Intermediate cases require numerical simulations, pursued since the early work of Friedlander (Reference Friedlander1957). To date, numerical solutions of heat transfer between a sphere immersed in a colder, flowing liquid are well established (Westerberg & Finlayson Reference Westerberg and Finlayson1990; Feng & Michaelides Reference Feng and Michaelides2000; Clift et al. Reference Clift, Grace and Weber2013), and a simple, closed-form approximation of numerical results provided by Clift et al. (Reference Clift, Grace and Weber2013) gives the Sherwood number ${\textit {Sh}}_{{Cl}}$ in terms of $\textit {Pe}$ with great accuracy, up to 2 %, for the entire range $0\leqslant {\textit {Pe}}\lt \infty$

(4.3) \begin{equation} {\textit {Sh}}_{{Cl}} = \frac {\varPhi _{{Cl}}}{\varPhi _{{D}}} = \frac {1}{2}\left (1+(1+2{\textit {Pe}})^{1/3}\right ){= \frac {1}{2}+\left ( \frac {1}{8}+\frac {{\textit {Pe}}}{4}\right )^{1/3}}. \end{equation}

These are in agreement with experimental results considering a heated sphere in slowly flowing fluid (Kramers Reference Kramers1946) and absorption of ions by a conducting sphere (Kutateladze, Nakoryakov & Iskakov Reference Kutateladze, Nakoryakov and Iskakov1982). However, we stress that all of these approaches assume a zero interaction range $\beta \to 0$ .

When advection dominates over diffusion (i.e. diffusion can be neglected) and the interaction range is finite ( ${\textit {Pe}} \to \infty$ , $\beta \gt 0$ , direct interception), a different family of approximations has been derived. In this regime, the Laplacian term in (3.6) can be neglected, and the stationary solution takes only two values: $0$ inside a critical streamline and $1$ outside of it. To determine $\varPhi$ , one simply calculates the cross-section of the critical stream tube in the far-field flow in (3.4). This gives an expression for direct interception flux $\varPhi _{{A}}$ as (Friedlander (Reference Friedlander1957), Supplementary Materials Sec. S.1)

(4.4) \begin{equation} \varPhi _{{A}} = U \pi b^2 \frac {3-\beta }{2}\varphi _0, \end{equation}

leading to the direct interception Sherwood number

(4.5) \begin{equation} {\textit {Sh}}_{{A}}=\frac {\varPhi _{{A}}}{\varPhi _{{D}}}={\textit {Pe}}\frac {\beta ^2(3-\beta )}{8}. \end{equation}

It is clear that (4.3) has different asymptotics than (4.5), as the former gives ${\textit {Sh}}_{{Cl}} \approx ( {\textit {Pe}}/4 )^{1/3}$ when ${\textit {Pe}}\to \infty$ . Therefore, predicting encounters for fast sinking particles based on the two models depends strongly on how the limit of high $\textit {Pe}$ is approached.

To better highlight the disagreement between the two models, we consider the number of collisions $\varLambda$ that occurred with a particle sedimenting over a fixed travel distance $\Delta Z$ (for example depth of the ocean). For a given diffusion coefficient $D$ , from (4.2) we get that

(4.6) \begin{equation} \varLambda =\varPhi _{{Cl}}\frac {\Delta Z}{U}=4\pi a^2 \varphi _0 \, \frac {\Delta Z}{{\textit {Pe}}}{\textit {Sh}}_{{Cl}}, \end{equation}

which is vanishingly small for large $\textit {Pe}$ because, from (4.3), ${\textit {Sh}}_{{Cl}} \approx ( {\textit {Pe}}/4 )^{1/3}$ as ${\textit {Pe}}\to \infty$ . In other words, the advection-diffusion model in (4.3) predicts that marine snow stops intercepting new particles when it becomes fast enough ( ${\textit {Pe}}\to \infty$ , by taking $U\to \infty$ with $D$ and $a$ constant, which can be achieved by increasing the density of the larger particle). In contrast, the number of collisions predicted by the direct interception model, (4.4), is always finite and given by $\varLambda =\varphi _0\pi b^2 (3-\beta ) \Delta Z / 2$ . This divergence of the predicted number of collisions between the two models arises because the advection-diffusion model assumes an infinitely small interaction range ( $\beta = 0$ ).

Thus, for any interaction range, however small, the interaction range cannot be neglected at very large values of $\textit {Pe}$ . The true asymptotics of ${\textit {Sh}}({\textit {Pe}})$ are determined by $\beta$ and, eventually, to maintain a constant, non-zero $\varLambda$ from (4.6), ${\textit {Sh}}({\textit {Pe}})$ should scale as ${\textit {Sh}}({\textit {Pe}}) \sim {\textit {Pe}}$ . To explore this discrepancy, we next calculate $\textit {Sh}$ numerically in the general case ( ${\textit {Pe}}\geqslant 0$ and $\beta \geqslant 0$ ).

5. Numerical investigation

5.1. Simulation methods

Equation (3.6) can be solved directly using finite element method (FEM) for moderate values of Péclet number such as $1 \lt {\textit {Pe}} \lt 10^{6}$ . We used the Python package scikit-fem (Gustafsson & McBain Reference Gustafsson and McBain2020), working in cylindrical symmetry. We expressed (3.6) in a weak form and obtained solutions for the concentration profiles around the sedimenting sphere, as shown in figure 2(b). We provide more details on this method of solution in figure S1 in the Supplementary Materials (SM) and present convergence tests in figure S2 and figure S3 in the SM. The concentration profiles were post-processed to obtain $\textit {Sh}$ . We did it in two ways to further validate the numerical solution. In the first approach, we calculated $\varPhi$ by integrating ${\textit {Pe}}\, \varphi (\rho ) \boldsymbol{u}(\rho )$ far downstream of the sphere. In the second approach, we calculated $\varPhi$ by integrating $({\textit {Pe}}\, \varphi \boldsymbol{u} - \boldsymbol{\nabla }\varphi ) \boldsymbol{\cdot }{\mathrm{d}} \boldsymbol{S}$ on the surface of the capturing sphere. Whenever the simulation box was sufficiently large, these methods were in agreement (see figure S4 in the SM). For smaller values of $\textit {Pe}$ (e.g. ${\textit {Pe}} \sim 0.1$ ) the FEM becomes impractical because it requires a very large simulation box to set the boundary conditions properly. However, this regime is well described by ${\textit {Sh}}_{{Cl}}$ , defined in (4.3).

We validated the FEM approach with solutions and experiments available in the literature for $\beta =0$ . The comparison of our numerical results for ${\textit {Sh}}({\textit {Pe}})$ is shown in figure 3. After harmonising the definitions of $\textit {Sh}$ and $\textit {Pe}$ between different sources and this publication, we see that our simulations are in good agreement with previous numerical and experimental works, and agree with the approximate relationship of (4.3).

Figure 3.

Validation of our numerical results in the case of zero interaction range. Comparison of our simulations with earlier experimental (Kramers Reference Kramers1946; Kutateladze et al. Reference Kutateladze, Nakoryakov and Iskakov1982) and numerical (Friedlander Reference Friedlander1957; Westerberg & Finlayson Reference Westerberg and Finlayson1990; Feng & Michaelides Reference Feng and Michaelides2000; Clift et al. Reference Clift, Grace and Weber2013) solutions of the advection-diffusion problem around a sphere in Stokes flow (thus with a zero interaction range, $\beta =0$ ). Closed-form approximation of Clift et al. (Reference Clift, Grace and Weber2013) is shown as dashed line. Results from our numerical model show excellent agreement with earlier numerical works (excluding the work of Friedlander (Reference Friedlander1957) which is an outlier). Deviations from experimental data are likely driven by finite Reynolds number effects (Kramers (Reference Kramers1946) measured $\textit {Sh}$ with $0.42\lt {\textit {Re}}\lt 26.5$ ). Different definitions of Péclet and Sherwood numbers were harmonised.

For even higher $\textit {Pe}$ , FEM calculations proved ineffective due to a numerical stability problem: for very large values of $\textit {Pe}$ , ringing artefacts appear near sharp gradients which follow the streamlines. Refining the mesh can extend the stability of the FEM method to slightly larger values of $\textit {Pe}$ ; similarly, more sophisticated FEM algorithms can somewhat alleviate the problems of sharp gradients. These improvements, however, are insufficient to cover the limit of very large $\textit {Pe}$ . Since the calculations for ${\textit {Pe}} \gt 10^7$ are key to understanding marine snow collisions, as seen in figure 1(c), we employed another simulation method.

Given that (3.5) is of Fokker–Planck type, it is possible to obtain integrals of $\varphi$ by simulating the corresponding Itô stochastic differential equation (SDE) with appropriate initial and boundary conditions, as described in Van Kampen (Reference Van and Godfried1992) and Öttinger (Reference Karp-Boss, Boss and Jumars1996). The simulated trajectories are used to calculate the hitting probability for a given initial position, which quantifies the number of small objects absorbed by the large particle. The total flux is then obtained by integrating the hitting probability over a disk placed upstream of the particle.

Specifically, the positions of suspended objects, $\boldsymbol{q}$ , are described by the equation

(5.1) \begin{equation} {\mathrm{d}} \boldsymbol{q} = \boldsymbol{u} {\mathrm{d}} t + \sqrt {2 D}\, {\mathrm{d}} \boldsymbol{W}_t, \end{equation}

with $\boldsymbol{q}(t=0) = \boldsymbol{q_0}$ , and where $\boldsymbol{W}_t$ is the Wiener process. Scaling the lengths by $(a+b)$ and the time by $(a+b)/U$ leads to the dimensionless expression

(5.2) \begin{equation} {\mathrm{d}} \boldsymbol{q} = \boldsymbol{u} {\mathrm{d}} t + \sqrt {\frac {2}{{\textit {Pe}}}}\, {\mathrm{d}} \boldsymbol{W}_t. \end{equation}

We simulated $N = 10^4$ particles for each of 50 starting points, $\boldsymbol{q_0} = (x_0,0,-h)$ , with $h = 5$ , using the Python package pychastic (Waszkiewicz et al. Reference Waszkiewicz, Bartczak, Kolasa and Lisicki2023a , Reference Waszkiewicz, Bartczak, Kolasa and Lisickib ). The typical resulting trajectories are shown in figure 2(c). The trajectories were then used to calculate the hitting probability $p_{\textit{hit}}(x_0)$ as a function of $x_0$ for an ensemble of initial locations $x_0$ . Specifically, $p_{\textit{hit}}(x_0)$ is the fraction of trajectories that approach the particle at a distance smaller than 1. The flux $\varPhi$ was calculated as

(5.3) \begin{equation} \varPhi =\varphi _0 U(a+b)^2 \int _0^\infty 2 \pi \rho \, p_{\textit{hit}}(\rho ) \, u_z(\rho ) \, {\mathrm{d}} \rho . \end{equation}

To efficiently calculate the integral, we sampled $p_{\textit{hit}}$ on a non-uniform grid (figure 4; for details, see discussion of figure S5 in the SM) and used linear interpolation of $p_{\textit{hit}}$ . We note that (5.3) includes the contribution of the advective flux ( $\boldsymbol{u}\varphi$ ) in (3.10) and ignores the diffusive component ( $D\boldsymbol{\nabla }\varphi$ ); This approximation utilises the fact that ${\partial }_z\varphi$ is negligible far downstream of the sphere.

In summary, FEM is an effective tool for simulating the dynamics of small objects for ${\textit {Pe}}\lt 10^{7}$ , while in the high-Pe regime, ${\textit {Pe}}\gt 10^5$ , SDE trajectories are increasingly practical due to faster convergence. We used the intermediate regime, where the two methods overlap, as a test ground to compare their results. In figure 4 we present an estimation of the hitting probability profile for ${\textit {Pe}}=7\times 10^5$ and $\beta =0.09$ using FEM and SDE. In this case, both methods converge. The details of the compatibility tests performed are included in the SM, in figure S6 and figure S7.

Figure 4.

Hitting probability of diffusive objects with a non-zero interaction range, calculated using SDE and FEM simulations. Numerical calculations for a representative, intermediate value of $\textit {Pe}$ and $\beta$ show good agreement between the FEM (scikit-fem) and SDE (pychastic) methods. In the SDE approach, each point was computed using $N=10^4$ trajectories. To reduce computational time, a set of 50 initial conditions was chosen. In a region near the critical streamline, sampling density was increased (more details in the SM, in figure S5 and figure S6). The profile of hitting probability was used to calculate the total flux onto the particle using (5.3).

5.2. Flux calculation – results

Figure 5.

Sherwood number for varying $\textit {Pe}$ and $\beta$ , obtained from numerical simulations. Filled dots correspond to FEM using the scikit-fem package, empty dots correspond to solutions of the SDE using the pychastic package. Colour denotes the value of $\beta$ . (a) Sherwood number, $\textit {Sh}$ , as a function of $\textit {Pe}$ . The solutions range from ignoring the radius of objects (aligning with the solution of Clift et al. (Reference Clift, Grace and Weber2013)) to ignoring diffusion of objects (parallel straight lines in the high- $\textit {Pe}$ limit). (b) Results of the same calculation presented in terms of the modified Sherwood number, $\widetilde {{\textit {Sh}}}$ , defined in (5.4) . With this parametrisation, all solutions approach $1$ as ${\textit {Pe}} \to \infty$ . (c) Relative error between numerical results for $\textit {Sh}$ and the analytical approximation ${\textit {Sh}}_{{f}}$ , proposed in (5.5).

The calculated flux $\varPhi$ reveals a different asymptotic scaling with large $\textit {Pe}$ for objects with a non-zero interaction range ( $\beta \gt 0$ ) compared with objects with zero interaction range ( $\beta =0$ ). We calculated $\varPhi$ for different values of $\textit {Pe}$ and $\beta$ and compared our calculations with the solution from Clift et al. (Reference Clift, Grace and Weber2013) ( $\beta =0$ ), given by (4.3) (figure 5 a). We find that, for small $\textit {Pe}$ , all solutions with different $\beta$ converge to the relation ${\textit {Sh}}_{{Cl}}$ . However, for higher $\textit {Pe}$ , the solutions diverge significantly from the curve with $\beta =0$ after a transition region, eventually reaching an asymptote different from ${\textit {Sh}}_{{Cl}} \approx ( {\textit {Pe}}/4 )^{1/3}$ . Instead, we observe that the solutions for $\beta \neq 0$ approach the asymptotics of ${\textit {Sh}} \sim {\textit {Pe}}$ in the high- $\textit {Pe}$ limit.

To better capture the transition between the diffusion and direct interception regimes, we introduce the modified Sherwood number $\widetilde {{\textit {Sh}}}$ as

(5.4) \begin{equation} \widetilde {{\textit {Sh}}} = \frac {\varPhi }{\varPhi _{{D}}+\varPhi _{{A}}}. \end{equation}

This extended definition ensures that all solutions for $\beta \gt 0$ tend to a finite limit as ${\textit {Pe}} \to \infty$ . As shown in figure 5(b), all numerical solutions for $\beta \gt 0$ eventually detach from ${\textit {Sh}}_{{Cl}}$ at a given $\textit {Pe}$ and tend to 1, thus transitioning from advection-diffusion with zero interaction range to pure direct interception. This definition thus allows us to assess the applicability of (4.3) and (4.5) simultaneously: advection-diffusion is a good estimate of the true encounter kernel as long as $\widetilde {{\textit {Sh}}}$ is close to (4.3) and the direct interceptions approach begins to be valid, when $\widetilde {{\textit {Sh}}}$ approaches the value of 1 for high $\textit {Pe}$ .

We provide two ways of incorporating the improved estimates of $\varPhi$ for future research. First, we developed a Python package pypesh accompanying this manuscript which, for given $\textit {Pe}$ and $\beta$ , calculates $\textit {Sh}$ as output interpolated from the calculations from figure 5(a). Second, we propose a closed-form approximation of the numerical results. We approximate $\varPhi$ as the sum of $\varPhi _{{Cl}}$ and $\varPhi _{{A}}$

(5.5) \begin{align} \begin{split} {\textit {Sh}}\approx {\textit {Sh}}_{{f}} &=(\varPhi _{{Cl}}+\varPhi _{{A}})/\varPhi _{{D}} \\ & ={\textit {Sh}}_{{Cl}} + \frac {U \pi b^2}{4 \pi D (a+b) } \frac {3-\beta }{2} \\ & = \frac {1}{2}\left (1+(1+2{\textit {Pe}})^{1/3}+ {\textit {Pe}}\:\frac {\beta ^2(3-\beta )}{4}\right ). \end{split} \end{align}

Figure 5(c) shows the relative difference between the numerically calculated value of $\textit {Sh}$ and the approximation of (5.5). We see that ${\textit {Sh}}_{\text{f}}$ underestimates the numerical value by at most 20 %, with the maximum of deviations shifting towards higher $\textit {Pe}$ for smaller values of $\beta$ . We thus conclude that, the two mechanisms generating encounters – direct interception and advection-diffusion – act approximately independently. In the solution for $\beta =0.001$ , a delicate discrepancy occurs, between the two numerical schemes, caused by numerical errors of the estimated $\textit {Sh}$ . This error, however, is smaller than 5 %. In summary, the formula in (5.5) provides a closed-form approximation of the flux as a function of $\textit {Pe}$ and $\beta$ .

6. Ecological implications for pico- and nanoplankton

Figure 6.

Advection-diffusion vs. direct interception for pico- and nanoplankton. Partial contribution to the total flux from $\varPhi _{{Cl}}$ (denoted as the advection-diffusion share) and $\varPhi _{{A}}$ (denoted as direct interception share). The white region represents a family of curves reflecting the partial contribution generated by varying $\Delta \rho$ and $a$ ( $30\; \textrm {kg}\,\textrm {m}^{-3}\lt \Delta \rho \lt 200\; \textrm {kg}\,\textrm {m}^{-3}$ and $20\;{\unicode{x03BC}} \textrm {m}\lt a\lt 10^3\;{\unicode{x03BC}} \textrm {m}$ ) in the Stokes law. An increase of $\Delta \rho$ and $a$ leads to an increase of $\varPhi _{{A}}$ in comparison with $\varPhi _{{Cl}}$ . See figure S9 in the SM for parameterisation using smaller density differences. The vertical dashed lines represent two examples of non-motile pico- and nanoplankton each: Prochlorococcus (image by A. Thompson, public domain (Thompson Reference Thompson2009)) and Pelagibacter (image by L. Steindler et al., CC BY 4.0 (Steindler et al. Reference Steindler, Schwalbach, Smith, Chan and Giovannoni2011)) as well as E. huxleyi (image by M. Iglesias-Rodriguez et al., CC BY 4.0 (Iglesias-Rodriguez et al. Reference Iglesias-Rodriguez, Jones, Blanco-Ameijeiras, Greaves, Huete-Ortega and Lebrato2017)) and Thalassiosira (image reproduced from Sumper & Brunner Reference Sumper and Brunner2008 with permission).

We now apply our results to environmentally realistic scenarios. We quantify the relative contributions of advection-diffusion (4.3) and direct interception (4.4) to the total encounter kernel (5.5) as a function of the size of the suspended objects in figure 6. We assume the smaller objects to be neutrally buoyant and have a diffusion coefficient described by the Stokes–Einstein relationship

(6.1) \begin{equation} D=\frac {k_{{B}}T}{6\pi \mu b}, \end{equation}

where $k_{{B}}=1.38\times 10^{-23}\;\textrm {J}\,\textrm {K}^{-1}$ is the Boltzmann constant, $T=277\;\textrm {K}$ is the approximate temperature of seawater and $\mu = 1.6\times 10^{-3}\;\textrm {Pa}\,\textrm {s}$ is the dynamic viscosity of seawater. To calculate $D$ for small objects, we make a common assumption that the interaction range and the size of the object are identical (Humphries Reference Humphries2009; Kiørboe et al. Reference Kiørboe, Grossart, Ploug and Tang2002; Jackson Reference Jackson1990; Burd & Jackson Reference Burd and Jackson2009; Alcolombri et al. Reference Alcolombri2025).

We assume for simplicity that larger particles sink according to the Stokes’ law

(6.2) \begin{equation} U = \frac {2 g \Delta \rho }{9\mu }a^2, \end{equation}

where $\Delta \rho$ is the density difference between water and the particles, and $g = 9.81\;\textrm {m}\,\textrm {s}^{-2}$ is the gravitational acceleration. We consider $ a$ from 20 to $ 10^3 \; {\unicode{x03BC}} \textrm {m}$ to cover a realistic range of marine snow sizes (DeVries et al. Reference DeVries, Liang and Deutsch2014). To mimic the variability in sinking velocities, we consider two ranges of $\Delta \rho$ . First, we set $ 30 \; \textrm {kg}\,\textrm {m}^{-3} \lt \Delta \rho \lt 200 \; \textrm {kg}\,\textrm {m}^{-3}$ based on previously reported values (Chase Reference Chase1979; McDonnell & Buesseler Reference McDonnell and Buesseler2010). Second, we consider much smaller values of $\Delta \rho$ to represent slowly sinking particles.

Our formula implies that encounters between sinking marine snow and submicrometre-sized suspended objects, such as non-motile picoplankton, are driven primarily by advection-diffusion, whereas the capture mechanism for larger objects, such as non-motile nanoplankton, depends sensitively on the properties of the sinking marine snow. Figure 6 shows the relative contribution of the two terms (advective-diffusive in blue and direct interception in orange) to the total flux as a function of the size of suspended objects for the case of marine snow density differences in the range of $ 30 \; \textrm {kg}\,\textrm {m}^{-3} \lt \Delta \rho \lt 200 \; \textrm {kg}\,\textrm {m}^{-3}$ . In this regime, the share of advection-diffusion vs. direct interception contributions takes a sigmoidal shape, centred at approximately $1.5\; {\unicode{x03BC}} \textrm {m}$ , a location that is close to the traditional size threshold separating marine microorganisms into picoplankton and nanoplankton (Omori & Ikeda Reference Omori and Ikeda1992). For a given size of suspended objects, the contributions of $\varPhi _{{Cl}}$ and $\varPhi _{{A}}$ are weakly dependent on the size of the marine snow particle $a$ and vary mainly with $\Delta \rho$ (with the range of relative shares of $\varPhi _{{Cl}}$ and $\varPhi _{{A}}$ constrained by the white region; see figure S8 in the SM for more details). For example, when the smaller object has a radius of $b \approx 1 \; {\unicode{x03BC}} \textrm {m}$ , the share of advection-diffusion is between $50\,\%$ and $83\,\%$ , and the share of direct interception is between $17\,\%$ and $50\,\%$ , depending on the density of the particle. On the other hand, for $b \approx 2.5 \; {\unicode{x03BC}} \textrm {m}$ , direct interception contributes between $64\,\%$ and $91\,\%$ of the total flux. For excess density of marine snow less than $ \Delta \rho = 30\;\textrm {kg}\,\textrm {m}^{-3}$ , the white region in figure 6 shifts towards larger values of $b$ , while maintaining its sigmoidal shape. In the case of $ \Delta \rho \approx 10^{-2}\;\textrm {kg}\,\textrm {m}^{-3}$ , which corresponds to $ U = 1.2\;\textrm {m}\,\textrm {day}^{-1}$ for $ a = 10^3\;{\unicode{x03BC}} \textrm {m}$ (such slowly sinking particles have also been observed in situ (Williams & Giering Reference Williams and Giering2022)), the 50 % threshold on the right-hand side of figure 6 occurs at $ b \approx 12\;{\unicode{x03BC}} \textrm {m}$ (see the figure S9 in the SM for more details). This estimation shows that, as the particle slows down, direct interception may underestimate the encounter rate with suspended objects as large as nanoplankton.