Modeling the dynamical sinking of biogenic particles in oceanic flow

We study the problem of sinking particles in a realistic oceanic flow, with major energetic structures in the mesoscale, focussing in the range of particle sizes and densities appropriate for marine biogenic particles. Our aim is to unify the theoretical investigations with its applications in the oceanographic context and considering a mesoscale simulation of the oceanic velocity field. By using the equation of motion of small particles in a fluid flow, we assess the influence of physical processes such as the Coriolis force and the inertia of the particles, and we conclude that they represent negligible corrections to the most important terms, which are passive motion with the velocity of the flow, and a constant added vertical velocity due to gravity. Even if within this approximation three-dimensional clustering of particles can not occur, two-dimensional cuts or projections of the evolving three-dimensional density can display inhomogeneities similar to the ones observed in sinking ocean particles.


Introduction
The sinking of small particles suspended in fluids is a topic of both fundamental importance and of practical implications in diverse fields ranging from rain nucleation to industrial processes (Michaelides, 1997;Falkovich and Fouxon, 2002).
In the oceans, photosynthesis by phytoplankton in surface waters uses sunlight, inorganic nutrients and carbon dioxide to produce organic matter which is then exported downward and isolated from the atmosphere (Henson et al., 2012), a process which forms the so-called biological carbon pump. The downward flux of carbon-rich biogenic particles from the marine surface due to gravitational settling, one of the key process of the biological carbon pump, is responsible (together with the solubility and the physical carbon pumps) of much of the oceans' role in the Earth carbon cycle (Sabine et al., 2004). Although most of the organic matter is metabolized and remineralized in surface waters, a significant portion sinks into deeper horizons. It can be sequestered on various time scales spanning a few years to decades in central and intermediate waters, several centuries in deep waters and up to millions of years locked up in bottom sediments (DeVries et al., 2012). Suitable modeling of the sinking process of particulate matter is thus required to properly assess the amount of carbon sequestered in the ocean and in general to better understand global biogeochemical cycling and its influence on the Earth climate. This is a challenging task that involves the downward transport of particles of many different sizes and densities by turbulent ocean flows which contain an enormous range of interacting scales. In the oceanographic community, numerous studies approached this problem by considering biogenic particles transported in oceanic flow as passive particles with an added constant velocity in the vertical to account for the sinking dynamics (Siegel and Deuser, 1997;Siegel et al., 2008;Qiu et al., 2014;Roullier et al., 2014;van Sebille et al., 2015). They suggest that the sinking of particles may not be strictly vertical but oblique, meaning that the locations where the particles are formed at the surface may be distant from the location of their deposition in the seafloor sediment. Then Siegel et al. (2008) presented the concept of statistical funnels which describe and quantify the source region of a sediment trap (subsurface collecting device of sinkingparticles used to get estimate of vertical fluxes). The validity of this approximation and the influence of different physical processes is however poorly discussed in these analyses.
In the physical community, the framework to model sinking particles is based on the Maxey-Riley-Gatignol equation for a small spherical particle moving in an ambient flow (Maxey and Riley, 1983;Gatignol, 1983;Michaelides, 1997;Provenzale, 1999;Cartwright et al., 2010), which highlights the importance of mechanisms beyond passive transport and constant sinking velocity, such as the role of finite size, inertia and history dependence. A major outcome of these studies is that inhomogeneities and particle clustering can arise spontaneously even if the fluid velocity field is incompressible and particles do not interact (Squires and Eaton, 1991). Particle clustering and patchiness is indeed observed in the surface and subsurface of the ocean (Logan and Wilkinson, 1990;Buesseler et al., 2007;Mitchell et al., 2008) Here we consider the theory of small but finite-size particles driven by geophysical flows, which is, as mentioned above, conveniently based on the Maxey-Riley-Gatignol equation. In Sect. 2 we review the main characteristics of marine particles which are relevant for their sinking dynamics. In Sect. 3 we present the equations of motion describing this process, together with the approximations required to obtain them and the type of particles for which they are valid. In particular, we discuss its validity and the relevance of the different physical processes involved in the range of sizes and densities of marine biogenic particles. In Sect. 4 we use these equations to study the settling dynamics in a modelled oceanic velocity field produced by a realistic high-resolution regional simulation of the Benguela upwelling system (southwest Africa). We estimate the relevance of physical processes such as the Coriolis force and the inertia of the particles with respect to the settling velocity. We also observe the spatial distribution of particles falling onto a plane of constant depth above the seabed and we identify clustering of particles that is interpreted with simple geometrical arguments which do not require physical phenomena beyond passive transport and constant terminal velocity. Our main results are finally summarized in a Conclusion section.

Characteristics of marine biogenic particles
In theory, the sinking velocities of biogenic particles depend on various intrinsic factors (such as their sizes, shapes, densities, porosities) which can be modified along their fall by complex bio-physical processes (e.g. aggregation, ballasting, trimming by remineralisation) as well as by the threedimensional flow field (Stemmann and Boss, 2012). However reasonable estimates of the effective sinking velocities of marine particles can be obtained by taking into account only its size and density (McDonnell and Buesseler, 2010). In our Lagrangian setting we thus consider that the two key properties of marine particles controlling their sinking dynamics are their size and density. Here we present the standard classification of marine particles according to the typical range of size and density by compiling different bibliographical sources.  Figure 1. Size and classification of marine particles (adapted from Simon et al. (2002)).

Size
Because of the diversity of the shapes, the size of a particle refers to the diameter of a sphere of equivalent volume (Equivalent Spherical Diameter) (Guidi et al., 2008). The size of marine particles ranges from 1 nm (almost-dissolved colloids) to aggregates larger than 1 cm (Stemmann and Boss, 2012). Originally, the size classification of particles was based on the minimal pore size of the nets used for their collection, which is about 0.45 − 1.0 µm. Any material larger than 0.2 µm (thus isolated by the filtration of seawater) is regarded as particulate organic matter, while the fraction that percolates through the filter is labelled as dissolved matter. This includes colloidal and truly dissolved materials (see Fig.  1). Although this discrimination of the size-continuum observed in the real ocean is somehow arbitrary, it is usefuland we will follow it-because particles smaller than 1.0 µm are not prone to sinking (Hedges, 2002).
In the following, our focus is thus on particulate matter larger than 1.0 µm (Fig. 1). Organic matter is produced in the sun-lit layer of the ocean by the primary production through photosynthesis of autotrophic microbes (mainly bacteria and phytoplankton). During their lifetime growth they exude colloidal and small particles to finally form larger particles when they die. Dead phytoplankton are within the range of 1 µm (picoplankton, e.g. cyanobacteria) and a few hundred of micrometers (microphytoplankton, e.g. diatoms).
Thereafter zooplankton consumes alive phytoplankton and inert particles and produce fecal pellets and dead bodies. Most fecal materials have enough size to sink rapidly by their own (De La Rocha and Passow, 2007). Typical sizes of such particles are 10 µm for a pellet of copepod of 200 µm length (Jackson, 2001), krill fecal pellets are between 160 µm − 460 µm (McDonnell and Buesseler, 2010) and euphausiid fecal pellets span 300 µm − 3 mm (Komar et al., 1981), providing the total range of 10µm to 3 mm. Concerning the zooplankton dead bodies, they are divided in micro-, meso-and macro-, with sizes in the range 20µm − 1cm. A detailed summary is given in Table 1.
Finally, there are the so-called organic aggregates which occur in the size range of 1µm to 10cm. They are typically formed in situ by physical aggregation or biological coagulation and are usually composed of numerous planktonic individuals and fecal pellets sticked together within a colloidal matrice. They are often distinguished in three size classes (Simon et al., 2002): macroscopic aggregates or macro-aggregates > 5mm usually called marine snow; microscopic, from 1 to 500µm, also known as microaggregates; and submicron particles < 1µm (which do not sink).

Density
The density of marine particles depends on their composition which can be divided into a mineral and a organic fraction (Maggi and Tang, 2015). The mineral or inorganic matter consists of biogenic minerals: Particulate Inorganic Carbon (PIC), e.g. calcium carbonate produced by coccoliths with density 2700 kg/m 3 and Biogenic Silica (BSi), produced by diatoms, significantly less denser than PIC, 1950 kg/m 3 (Balch et al., 2010). The density of Particulate Organic Matter (POC) ranges widely depending on its origin. For instance, the density of cytoplasm spans from 1030 to 1100 kg/m 3 , while the one of fecal pellets ranges between 1174 kg/m 3 and 1230 kg/m 3 (Komar et al., 1981). Despite this variability, it is possible to assign a range to the density of organic matter, from 1050 to 1500 kg/m 3 .
Considering all these estimates together, the density of marine particle ranges approximately between 1050 to 2700 kg/m 3 (Maggi, 2013). This should be compared to standard values for sea water density in the interior ocean which spans roughly 1020-1030 kg/m 3 . Thus most of the particle types described previously will sink. Assuming constant size and density for each particle along its downward course, we deduce that most of the particles types described previously will sink. This holds without considering biogeochemical and (dis)aggregation processes that may occur in nature, thus lowering the particle density and resulting in clustering and trapping of particles at particular isopycnals (Sozza et al., 2016). Note that we do not consider here living organisms which show vertical movements by active swimming or by controlling their buoyancy (Moore and Villareal, 1996;Azetsu-Scott and Passow, 2004).
3 Equations of motion for small spherical rigid particles

The Maxey-Riley-Gatignol equation
To describe the sedimentation of biogenic particles, we need to study the motion of single particles driven by fluid flow. A milestone to analyze the dynamics of a small spherical rigid particle of radius a subject to gravity acceleration g in an unsteady fluid flow u(r, t) is given by the Maxey-Riley-Gatignol (Maxey and Riley, 1983;Gatignol, 1983;Michaelides, 1997;Cartwright et al., 2010) equation (MRG in the following): (1) The velocity of the particle is denoted by v = v(t). The particle and fluid densities are ρ p and ρ f , respectively, and ν denotes the fluid kinematic viscosity. The time derivative operators d dt = ∂ ∂t + v · ∇ and D Dt = ∂ ∂t + u · ∇ denote the time rate of change following the particle itself and the time rate of change following a fluid element in the undisturbed flow field u(r, t) respectively. This equation of motion gives the balance between the different forces acting on the particle, which correspond to the right-hand-side terms: the pressure force (the force exerted on the particle by the undisturbed flow), the buoyancy force, the drag force (Stokes drag), the added mass force resulting from the part of the fluid moving with the particle, and the history force. As will be discussed below the validity of this equation requires several conditions, being the main one the small size of the particles. The terms with a 2 ∇ 2 u are the Faxén corrections (Faxén, 1922).
The full MRG is very complicated to manage. A further simplification is usually performed based on the single assumption of very small particles (what this exactly means will be discussed later on). With this, the Faxén corrections and, as commented below, also the history term (since a/ √ ν << 1) can be neglected (Maxey and Riley, 1983;Michaelides, 1997;Haller and Sapsis, 2008). Note however that the history term can be relevant under some conditions, as for example larger particle size (Daitche and Tél, 2011;Guseva et al., 2013Guseva et al., , 2016Olivieri et al., 2014). Thus we obtain the standard form of the MRG equations (Maxey and Riley, 1983): where β = 3ρ f 2ρp+ρ f , the Stokes time is τ p = a 2 3βν , and v s = (1 − β)gτ p is the settling velocity in quiescent fluid. Equation (2) is the starting point for most inertial particle studies (Michaelides, 1997;Balkovsky et al., 2001;Cartwright et al., 2010).
-Submicron: size < 1 µm. Table 1. Simplified categorization of marine biogenic particles, and their associated sizes. and densities of marine organisms. We do so in the context of open-ocean flows, which are typically most energetic at the mesoscale (scales of about 100 km), and where there is a strong stratification, with vertical velocities three or four orders of magnitude smaller than horizontal ones. The motion becomes more three-dimensional, and then the concepts of three-dimensional turbulence more relevant, below scales l of some hundred of meters, with typical velocities decreasing as l 1/3 for decreasing scale and velocity gradients increasing as l −2/3 until the Kolmogorov scale l = η below which flow becomes smooth. Because of its direct exposure to wind, turbulence intensity is typically larger at the ocean surface, with values of turbulent energy dissipation in the range 1·10 −6 m 2 /s 3 < < 3·10 −5 m 2 /s 3 (Jimenez, 1997), than at depth. The first condition for the validity of the MRG equation that was originally discussed by Maxey and Riley (1983) is that the particles have to be much smaller than the typical length scale of variation of the flow. This means that for multiscale (turbulent) flows the radius of the particle a has to be much smaller than the Kolmogorov scale η, which according to the previous values of , is typically 0.3mm < η < 2mm in the ocean surface (Okubo, 1971;Jimenez, 1997). Note that we only have to consider worst-case situations for assessing the validity of the different approximations. Another condition to be fulfilled is that the shear Reynolds number must be small Re ∇ = a 2 U/νL << 1, where U and L are typical velocity and length scales. For a turbulent ocean with multiple scales and velocities, the most restrictive condition arises when they take the values of the Kolmogorov velocity v η and length η, respectively, since then the velocity gradients are maxima. In this case the condition becomes Re ∇ = a 2 /η 2 << 1, which again is satisfied for small particles. We note that Guseva et al. (2013) found that the relative importance of the history term in Eq. (1) with respect to the drag force is of the order of a parameter which in our notation is (Re ∇ ) 1/2 . This justifies neglecting the history term for small particles, although its importance increases for increasing size (Daitche and Tél, 2011;Guseva et al., 2013). Another condition to be satisfied for the validity of the MRG equation is that the so-called Reynolds particle number, Re p = a|v−u| ν should fulfill Re p << 1. Considering that gravity force dominates over other forces one has |v − u| |v s | ≡ v s , where v s is, as introduced before, the settling velocity of particles in a quiescent fluid due to Stokes drag. The Reynolds particle number is then Re p = avs ν . Note that the settling velocity depends only on the densities of particles via the parameter β. Assuming a mean density of sea water in the upper ocean as ρ f = 1025kg/m 3 the parameter β has values within the range [0.5, 0.99] for the typical values of the density of marine particles previously discussed. Fig. 2 shows v s for different sizes and the regions where Re p > 1 (and other parameter regions where MRG is not a good ap- Sinking velocity versus particle radius for different β, which is determined by densities. The blue zone determines the values of the settling velocities at a given radius, as determined by the typical marine particle densities. The green area is determined by the condition Rep > 1 for which the MRG equation is not valid. Use of the MRG equation is also unjustified for particles larger than the Kolmogorov length scale also plotted in the figure. We also show the region τp > τη ≈ 1s where the additional approximation leading to Eq. (6) becomes invalid.
proximation) as a function of particle radius and for the limiting values of β. It reveals that Eq. (1) can not describe ocean particles larger than 300µm of any density, and for a limited range of densities when the particle radius exceeds approximately 100µm. In fact, the range of application of MRG to marine particles is plotted in the blue area, which at the same time gives an estimate of the typical sinking velocities for a given particle size. Summarizing, both the MRG and its approximation Eq. (2) are valid for marine particles with size within the range 1µm and 200µm. That is, it is valid for all particulate organic matter in Fig. 1 except the largest of the micro-aggregates and meso-and macro-bodies of zooplankton. The sinking velocities range from 1mm/day-1km/day.

The MRG equation in a rotating frame and further simplications
We are interested in applying Eq.
(2) in oceanic flows, where the particle v and flow u velocities are expressed in a frame rotating with the Earth angular velocity Ω (Elperin et al., 2002;Biferale et al., 2016;Tanga et al., 1996;Provenzale, 1999;Sapsis and Haller, 2009). Both time derivatives d dt and D Dt have to be corrected following the rule Where Ω = |Ω| and r is the particle position vector whose origin is in the rotation axis. So that Eq.
(2) is now Two apparent forces arise in the equation, the Coriolis force 2Ω × (v − βu) and the centrifugal force, which is included in a modified sinking velocity v s = (1 − β)(g − Ω × (Ω × r))τ p . The effect of the centrifugal force is very small (of order 10 −3 compared to gravity) and can be absorbed in a redefinition of g. Thus, in the following we take v s = v s with the properly chosen g.
The ratio between the particle response time and the Kolmogorov time scale is the Stokes number St = τ p /τ η , which measures the importance of particle's inertia because of its size and density. According to the range of in the ocean mentioned before, we get 0.1 s < τ η < 5 s, and for our range of particle sizes 10 −6 s < τ p < 10 −2 s so we can assume that St << 1 (see Fig. 2). This motivates us to make a second (standard) approximation (Balkovsky et al., 2001;Haller and Sapsis, 2008) of the MRG equation expanding in powers of τ p (note that it would be more natural to make the expansion in powers of the non-dimensional St but we prefer to do it in τ p to control on the time scales of the problem). Assuming first the solution to Eq. (2): v = u + u 1 τ p + u 2 τ 2 p + . . . , and using dv dt = Du Dt +O(τ p ), we get that the particle velocity at first order in τ p is It is worth recalling that τ p (1 − β) = v s /g, so that all dependencies on particle size and density appear in Eq. (6) through the combination of parameters defining v s . Different combinations of size and density, taken within the ranges reported in Sect. 2, follow the same dynamics if they have the same undisturbed settling velocity v s .
A further discussion of Eq. (6) follows. At this order only three physical processes correct the particle velocity with respect to the fluid velocity: the Stokes friction determining the settling velocity v s , the inertial term given by τ p (β − 1) Du Dt whose major effect is to introduce a centrifugal force pulling particles away from vortex cores (Maxey, 1987;Michaelides, 1997), and the influence of the Coriolis force 2τ p (β − 1)Ω × u. Concerning sinking dynamics, the v = u + v s is the most relevant approximation, and many other studies consider it, mainly in oceanographic contexts (e.g. Siegel and Deuser, 1997). Note that we can use the right-hand-side of Eq. (6) with u = u(r, t) to define the particle velocity v as a velocity field in three-dimensional space v = v(r, t). If one uses the lowest-order approximation v ≈ u we have ∇·v = ∇·u = 0 when the fluid velocity field u is incompressible (which is the case for ocean flows). This means that when considering this term alone, one cannot obtain a compressible particle velocity whereas this was the main reason invoked to explain the clustering of finite-size particles (Squires and Eaton, 1991;Bec, 2003). For this reason, numerous studies (Tanga et al., 1996;Michaelides, 1997;Bec et al., 2007Bec et al., , 2014Cartwright et al., 2010;Guseva et al., 2013;Gustavsson et al., 2014;Beron-Vera et al., 2015) consider the role of the additional terms. With them ∇ · v = τ p (β − 1)∇ · ( Du Dt + 2Ω × u) = 0, and inertia-induced clustering may occur. In the following sections we address two main questions: a) how relevant for the sinking dynamics are the Coriolis and centrifugal terms?; and b) are they essentical ingredients for the clustering of biogenic particles? We will study the relevance of the different terms in Eq. (6) in a realistic oceanic setting.

Numerical simulations
The velocity flow u of the Benguela region was produced by a regional simulation of a hydrostatic free-surface primitive equations model called ROMS (Regional Ocean Modelling System). The configuration used here extends from 12°S to 35°S and from 4°E to 19°E (blue rectangle in Fig. 3) and was forced with climatological atmospheric data (Gutknecht et al., 2013). The simulation area extends from 12 • S to 35 • S and from 4 • E to 19 • E (blue rectangle in Fig. 3). The velocity field data set consists of 2 years of daily averages of zonal (u), meridional (v) and vertical velocity (w) components, stored in a three-dimensional grid with a horizontal resolution of 1/12 o and 32 vertical terrain-following levels using a stretched vertical coordinate where the layer thickness increases from surface/bottom to the ocean interior. In order to integrate particle trajectories from the velocity in Eq. (6) we interpolate linearly u(r, t) from the closest space-time grid points to the actual particle locations. Given the huge disparity between the model resolution and the small particle-sizes considered, it is pertinent to parameterize in some way the unresolved scales. This can be done by different approaches, from stochastic Lagrangian modeling (Brickman and Smith, 2002), to deterministic kinematic fields (Palatella et al., 2014). The first approach is adopted by adding a simple white noise to the particle velocity (Tang et al., 2012), with different intensity in the vertical and horizontal directions. Thus, we consider this noisy version of the a three-dimensional vector Gaussian white noise with zero mean and correlations W i (t)W j (t ) = δ ij δ(t − t ), i, j = x, y, z. We consider an horizontal eddy diffusivity, D h , depending on resolution length scale l according to Okubo formula (Okubo, 1971;Hernandez-Carrasco et al., 2011): D h (l) = 2.055 × 10 4 l 1.55 (m 2 /s). Thus, if taking l ∼ 8 km = 8000 m (corresponding to 1/12 • ) we obtain 10m 2 /s. In the vertical direction we use a constant value of D v = 10 −5 m 2 /s (Rossi et al., 2013).
In order to obtain quantitative assessment of the relative effects of the different physical terms in Eq. (8), we will compare trajectories obtained from the following expressions which only consider some of the terms of the full expression Eq. (8): Besides the random noise term, the first expression (9) only considers the settling velocity, equation (10) resolves the set-tling velocity plus the Coriolis effect, and equation (11) considers the settling plus the inertial term.
For the numerical experiments we will consider a set of six values of v s ranging from 5m/day to 200m/day, with different integration times to have in all the cases a sinking to about 1000−1100 m depth. The stochastic equation (7) with expressions (8)-(11) is written in spherical coordinates and numerically integrated using a second-order Heun's method with time step of 4 hours (Toral and Colet, 2014). We use R = 6371 km for the Earth radius, g = 9.81m/s 2 , and the angular velocity Ω is a vector pointing in the direction of Earth axis and modulus |Ω| = 7.2722×10 −5 s −1 . We take v s and τ p constant in each experiment because, although water density may increase with depth, this variation is at most of 10kg/m 3 in the range of depths we are considering here and then the impact on v s is below 0.1%. We use as initial starting date 17 September 2008. The numerical experiments consist in launching N = 6000 particles from initial conditions randomly chosen in a square of size 1/6 • centered at 10.0 • E 29.12 • S and −100.0m depth (red rectangle in Fig. 3), and in letting them evolve for a given time t f (stated in Table 2) following Eq. (7) with expressions (8)-(11). We use in each case identical initial conditions and the same sequence of random numbers for the noise terms. In this way we guarantee that any difference in particle trajectories arise from the inclusion or not of the inertial and Coriolis terms. We obtain the time-dependent positions of all the particles for each approximation to the dynamics: r i (t), r    Table 2 gives the mean and the standard deviation of the depths attained by the set of particles in each numerical experiment as obtained from Eqs. (7) and (8). We find that the use of the different approximations (9)-(11) gives virtually the same results. The only differences larger than 1 cm in mean or standard deviation are the ones for the smallest unperturbed settling velocity considered, v s = 5 m/day, and are also reported in Table 2. The measured differences are negligible as compared with the traveled distance or even with the model grid size. Indeed small changes in the ROMS model configuration or in the velocity interpolation procedure would have an impact larger than this. The mean displacements in the horizontal obtained with the different approximations (not shown) differ also in less 0.1%. We thus conclude that the simplest approximation Eq. (9) which only considers passive transport and an added constant sinking velocity already provides a good description of the sinking process for the type of marine particles and the range of space and time scales considered here. Note that the depth attained by the particles is always slightly shallower than z = −1100 m, which is the depth that would be reached in a still fluid. It is still debated under which conditions fluid flows enhances or reduces the settling velocity (Maxey, 1987;Wang and Maxey, 1993;Ruiz et al., 2004;Bec et al., 2014).  (7) and (8), by the set of particles released from the red rectangle in Fig. 3 at z = −100 m for the different values of vs and integration times used. The results labeled (co), (in), and (0) are obtained from the different approximations in Eqs. (9)-(11), which differ more than 1 cm from the ones obtained from Eq. (8) only in the vs = 5 m/day case.
We perform now a more stringent test going beyond the analyses of mean displacements by considering differences between individual particle trajectories. To assess the impact of the Coriolis and of the inertial effects we compare the positions r   i (t) for each time t. To do so we compute the root mean square difference in position per particle and time, which we separate in vertical and horizontal components: with x i = (x i , y i ), the horizontal position vectors, and the superindex (k) takes the values (co) or (in). Fig. 4 shows the influence of the Coriolis term in the horizontal component for each sinking velocity as a function of time. We observe an exponential growth in a wide range of times, which reveals the chaotic behavior of each of the compared trajectories. The value of the exponent 0.08days −1 is in agreement with the order of magnitude of the Lyapunov exponent calculated using the same ROMS velocity model and region (Bettencourt et al., 2012). Similar exponential growth with the same growth rate were observed for the inertial terms and the vertical components (not shown), although the absolute magnitude of these mean root square differences was much smaller.
The horizontal and vertical differences r (co) h,v at the final integration time t f (i.e. the time at which the particles reach an approximate depth of 1000 m for each value of v s ) are displayed in Fig. 5, both as a function of v s and of t f . Similarly, the values of r (co) h,v are presented in Fig. 6. The behavior can be understood as resulting from two factors: on the one hand smaller v s requires larger t f to reach the final depth, and larger integration time t f allows for accumulation of larger differences between trajectories. On the other hand the Coriolis and inertial terms in Eqs. (10)-(11) are proportional to τ p (β − 1) = v s /g so that their magnitude decreases for smaller v s . The combination of these two competing effects shapes the curves in Figs. 5 and 6, which for the vertical-difference case turn-out to be non-monotonic in v s or t f .
In all cases, the differences (both in vertical and horizontal) between the simple dynamics (9) and the corrected ones in Eqs. (10) and (11) are negligible when compared with typical particle displacements, or even with model grid sizes. For example, we imposed in our simulations a vertical displacement close to 1000 m, whereas the mean root square difference with respect to simple sinking is below 1 m for the Coriolis case (Fig. 5) and below 1 cm for the inertial case (Fig. 6). In the horizontal direction, displacements during those times are of the order of hundreds of km, whereas the corrections introduced by the Coriolis and inertial terms are in the worst cases of the order of a few kilometers or of tens of meters, respectively. In particular, the most important impact (horizontal differences of tens of kilometer) is attributed to the Coriolis term for particles sinking at 5 m/day (Fig. 5). This indicates that the inclusion of the Coriolis term would be required to properly model slowly sinking particles at high latitudes. It is worth noting that although the small value of Rossby number 0.01 for mesoscale processes might indicate a strong influence of the Coriolis force in Eq. (8), its influence on particle dynamics becomes negligible because it is multiplied by τ p or equivalently, the Stokes number, which is significantly small for biogenic particles. Nevertheless Rossby number coincides with the ratio of inertial term to Coriolis term in Eq. (8) and its value 0.01 explains the difference of two orders of magnitude among the corrections arising from the inertial force and from Coriolis. The trajectories of the full dynamics ruled by Eq. (8) are nearly identical to the ones under the approximation which keeps only the sinking term and Coriolis, so that the corresponding comparison to r (0) i gives a figure essentially identical to Fig. 5 (not shown).
In summary, for the range of sizes and densities of the marine particles considered here, the sinking dynamics is essentially given by the velocity v = u + v s , which has been the one used in some oceanographic studies (Siegel and Deuser, 1997;Siegel et al., 2008;Roullier et al., 2014). Note however that a new question arises: what is then the reason for the observed clustering of falling particles (Logan and Wilkinson, 1990;Buesseler et al., 2007;Mitchell et al., 2008)? The argument of the non-inertial dynamics of the particles does not  Table 2) computed with and without the Coriolis term (Eqs. (10) and (9), respectively). Data are presented as a function of the unperturbed sinking velocity vs used (upper horizontal scale) and of the final integration time t f (lower horizontal scale). Upper violet line, the horizontal difference r  serve since ∇·v = ∇·u = 0. A possible response is explored in the next section.

Geometric clustering of particles
Compressibility of the particle-velocity field, i.e. ∇ · v = 0, which can arise from inertial effects even when the corresponding fluid-velocity field is incompressible ∇ · u = 0, has been identified as one of the mechanisms leading to preferential clustering of particles in flows (Squires and Eaton, 1991;Balkovsky et al., 2001). This is so because ρ(t), the particle  Figure 6. Root mean square difference per particle between final positions (at times t f stated in Table 2) computed with and without the inertial term (Eqs. (11) and (9)   density at time t at the location r = r(r 0 , t) of a particle that started at r 0 at time zero, satisfies ρ(t) = ρ(0)δ −1 , where δ is a dilation factor equal to the determinant of the Jacobian | ∂r ∂r0 |, which satisfies 1 δ or, using δ(0) = 1: Thus, particles will accumulate (i.e. higher ρ(t f )) in final deep locations receiving particles whose trajectories have predominantly travelled through regions with ∇ · v < 0. We have seen however that, to a good approximation ∇ · v ≈ ∇ · u = 0 since inertial effects can be neglected for the type of marine particles we consider here, and then the threedimensional particle-velocity field is incompressible. We now reproduce numerically a typical situation in which clustering of marine particles is observed. We release particles uniformly in an horizontal layer close to the surface, we let them sink within the oceanic flow and we finally observe the distribution of the locations where they touch another horizontal deeper layer. The domain chosen is the rectangle 12 • S to 35 • S and 4 • E to 19 • E (orange rectangle in Fig. 3). We divide the domain horizontally in squares of side 1/25 • , then initialize 1000 particles at random positions in each of them in August 20, 2008 at depth z = −100 m (i.e. the bottom of the euphotic layer, starting point of our biogenic particles), and then integrate each trajectory until it reaches −1000m depth. We use expression (9) for the velocity, with v s = 50m/day. In order to avoid any small fluctuating compressibility arising from the noise term we put W = 0 but we have checked that the result in the presence of noise is virtually indistinguishable (not shown). At the bottom layer (z = −1000 m) we count how many particles arrive to each of the 1/25 • boxes and display the result in Fig.  7(a). Despite ∇ · v = 0 we see clear preferential clustering of particles in some regions related to eddies and filaments. We note that our horizontal boxes have a latitude-dependent area so that distributing particles at random in them produces a latitude-dependent initial density which could lead to some final inhomogeneities. We have checked however that for the range of displacements of the particles, this effect is everywhere smaller than 5% and thus can not be responsible for the large clustering observed in Fig. 7(a). Nevertheless, this effect will be taken into account later. We explain the observed particle clustering by considering the field displayed in Fig. 7(a) as a projection in two dimensions of a density field (the cloud of sinking particles) which evolves in three-dimensions. Even if the threedimensional divergence is zero, and then an homogeneous three-dimensional density will remain homogeneous, a twodimensional cut or projection can be strongly inhomogeneous. This mechanism has been proposed to explain clustering and inhomogeneities in the ocean surface (Huntley et al., 2015;Jacobs et al., 2016), but we show here that it is also relevant for the crossing of a horizontal layer by a set of falling particles.
A simple way to confirm that this clustering arises from the two-dimensionality of the measurement is to estimate the changes in the horizontal density of evolving particle layers as if they were produced just by the horizontal part of the velocity field. This is only correct if an initially horizontal particle layer remains always horizontal during the sinking process, which is not true. But, given the huge differences in the values of the horizontal and vertical velocities in the ocean, we expect this approximation to capture the essential physics and provide a qualitative explanation of the observed clustering. We expect the approximation to become better for increasing v s , because of the shorter sinking time during which vertical deformations could develop. Thus we compute the two-dimensional version of the dilation field, δ h (x, t f ), at each horizontal location x in the deep layer at z = −1000m: where in the second equality we have used Eq. (9) from which ∇ h · v = ∇ h · u and the third one is a consequence of ∇·u = 0. In order to get the values of δ h on a uniform grid on the −1000m depth layer at the arrival date t f of the particles in the previous simulation, we integrate backwards in time trajectories from grid points separated 1/50 • at z = −1000m until they reach −100m. The starting date (t f ) of the backwards integration was September 7, 2008, i.e. 18 days after the release date used in the previous clustering experiment. This value correspond to the average duration time of trajectories in that experiment. Then δ h was computed integrating in time the values of ∇ h · v along every trajectory using Eq. (16). Figure  7(b) displays the quantity δ(x, t f ) −1 cos(θ f )/ cos(θ 0 ), which gives the ratio between densities in the upper and lower layer, corrected with the angular factors controlling the area of the horizontal boxes so that this can be compared with the ratio between particle numbers displayed in Fig. 7(a). θ f is the latitude of point x, and θ 0 is the latitude of the corresponding trajectory in the upper z = −100 m layer. As stated before, the latitudinal corrections by the cosine terms are always smaller than a 5%. Although there is no perfect quantitative agreement, there is clear correspondence between the main clustered structures in panels (a) and (b) of Fig. 7, confirming that they originate from the horizontal dynamics in an incompressible three-dimensional velocity field. We have checked in specific cases that locations with larger differences between Figs. 7(a) and (b) correspond to places with large dispersion in the arrival times to the bottom layer, indicating deviations from the horizontality assumption.

Conclusions
We have studied the problem of sinking particles in a realistic oceanic flow, focussing in the range of sizes and densities appropriate for marine biogenic particles. Starting from a modeling approach in terms of the MRG equation (1), our conclusion is that the simplest approximation given by Eq. (9) in which particles move passively with the fluid flow with an added constant settling velocity in the vertical direction is an accurate framework to describe the sinking process in the type of flows and particles considered. A re-assessment of these assumptions may be required if more complex processes (such as aggregation/disaggregation) are included and when super-high resolution (submesoscale and below) mimicking the real ocean will become available.
Corrections arising from the Coriolis force turn out to be about 100 times larger than the ones coming from inertial effects, in agreement with the results in Sapsis and Haller (2009) or in Beron-Vera et al. (2015), but both of them are negligible when compared to the effects of passive transport by the fluid velocity plus the added gravity term, except for very slowly sinking particles in high latitudes.
If the fluid flow field u(r, t) has vanishing divergence then the same is true for the particle velocity field defined by the approximation in Eq. (9). Then, no three-dimensional clustering can occur within this approximation. Nevertheless, we have shown that two-dimensional cuts or projections of evolving three-dimensional particle clouds display horizontal clustering.