Articles | Volume 17, issue 2
Research article
09 Mar 2021
Research article |  | 09 Mar 2021

Sinking microplastics in the water column: simulations in the Mediterranean Sea

Rebeca de la Fuente, Gábor Drótos, Emilio Hernández-García, Cristóbal López, and Erik van Sebille

We study the vertical dispersion and distribution of negatively buoyant rigid microplastics within a realistic circulation model of the Mediterranean sea. We first propose an equation describing their idealized dynamics. In that framework, we evaluate the importance of some relevant physical effects (inertia, Coriolis force, small-scale turbulence and variable seawater density), and we bound the relative error of simplifying the dynamics to a constant sinking velocity added to a large-scale velocity field. We then calculate the amount and vertical distribution of microplastic particles on the water column of the open ocean if their release from the sea surface is continuous at rates compatible with observations in the Mediterranean. The vertical distribution is found to be almost uniform with depth for the majority of our parameter range. Transient distributions from flash releases reveal a non-Gaussian character of the dispersion and various diffusion laws, both normal and anomalous. The origin of these behaviors is explored in terms of horizontal and vertical flow organization.

1 Introduction

Approximately 8 million tonnes of plastics end up in the oceans every year (Jambeck et al.2015). Nevertheless, only a very small percentage, around 1 %, remains on the surface (van Sebille et al.2015; Choy et al.2019). The rest leaves the surface of the ocean (Ballent et al.2013; van Sebille et al.2020) through beaching (Turner and Holmes2011), biofouling (Ye and Andrady1991; Chubarenko et al.2016; Kooi et al.2017) or sinking (Erni-Cassola et al.2019), but also wind-driven mixing presumably leads to an underestimation for the amount of particles remaining close to the sea surface (Kukulka et al.2012; Enders et al.2015; Suaria et al.2016; Poulain et al.2018). The distribution of plastic pollution in the sea is poorly understood at present but would be crucial to properly evaluate the exposure of marine biota to this material and formulate strategies for cleaning the oceans (Horton and Dixon2018).

Floating plastics and those that have beached or sedimented on the seafloor are relatively well studied through field campaigns (although an explanation is missing for many findings; Andrady2017; Erni-Cassola et al.2019; Kane and Clare2019). In contrast, the presence of plastics within the water column has received less attention, and many surveys in this realm are restricted to so-called underway samples, a few meters below the surface (e.g., Enders et al.2015). However, for example, Choy et al. (2019) reported that below the mixed layer and down to 1000 m depth in Monterey Bay, concentrations of plastics are larger than at the surface (Thompson et al.2004; Hidalgo-Ruz et al.2012). Egger et al. (2020) found more plastic between 5 and 2000 m below the North Pacific Garbage Patch than at the surface. These findings turn out to mostly concern plastic pieces that, according to their nominal material density, would be classified as positively buoyant (Egger et al.2020).

In this paper, we focus on a certain class of plastic particles, negatively buoyant rigid microplastics, excluding very small sizes, and we estimate their vertical distribution through the water column and their amount in the Mediterranean Sea. Microplastic particles are among the most important contributors to marine plastic pollution (Arthur et al.2009). Closely following the work of Monroy et al. (2017) for sinking biogenic particles but choosing particle properties to correspond to those of negatively buoyant microplastics, we first justify the use of a simplified equation of motion, in which the plastic particle velocity is the sum of the ambient flow velocity and a sinking velocity depending on particle and water characteristics. In particular, we estimate the impact of some corrections to this simple dynamics and evaluate in detail the influence of the spatial variation of the seawater density on the plastic dispersion and sinking characteristics. For our Mediterranean case study, the impact of the varying seawater density on particle trajectories can be comparable to the estimated effect of the neglected small scales below the hydrodynamical model's resolution.

We then estimate the amount of microplastic particles in the water column of the open Mediterranean. Our estimates rely on a uniform vertical distribution, which is confirmed by our numerical simulations to be a good approximation for fast-sinking particles. This can be explained by a simple model in which released particles sink with a constant velocity. Detailed consideration of the transient dynamics identifies small non-Gaussian vertical dispersion around this simple sinking behavior, with transitions between anomalous and normal effective diffusion.

2 Types of microplastics in the water column

The dynamics and the fate of microplastics in the ocean are largely determined by their material density (Erni-Cassola et al.2019). However, shape, size and rigidity are also relevant properties, characteristic transport pathways to the water column being different for different particle types.

Typically, positively buoyant plastic types will remain floating at the sea surface or close to it and then will not contribute to the microplastic content in the water column, the topic we are interested in this paper. However, it has been documented experimentally that biofouling may increase sinking rates of particles up to 81 % and enhances sedimentation (Kaiser et al.2017). So, a class of high abundance and mass may be represented by nearly neutrally buoyant microplastic particles that are generated by biofouling (Ye and Andrady1991; Chubarenko et al.2016) from positively buoyant plastic types or by other mechanisms of aggregation with organic matter, especially for small particle sizes (Kooi et al.2017).

In fact, the fallout from the North Pacific Garbage Patch almost entirely consists of plastic types nominally less dense than water (Egger et al.2020). Although some of these immersed particles finally reach the sea bottom, their proportion in sedimented plastic is minor except for the immediate vicinity of coasts where water is shallow. Most of these particles remain in the photic zone (Mountford and Morales Maqueda2019; Wichmann et al.2019; Soto-Navarro et al.2020). This suggests that reverse processes could also take place after biofouling and that the dynamics of such particles is complicated (Kooi et al.2017; Erni-Cassola et al.2019).

Particles denser than seawater dominantly accumulate at the sea bottom (Mountford and Morales Maqueda2019). A mechanism by which microplastics denser than water can also be present within the water column is the finite time taken by them to reach the bottom. Under continuous release at the surface and sedimentation at the bottom, the transient falling would lead to a steady distribution for the amount of plastic in the water column at any given time, and this distribution has never been estimated. Note that the Eulerian methodology of Mountford and Morales Maqueda (2019), treating sedimentation (i.e., deposition on the seafloor) by parametrization and thus leaving particles in the water column indefinitely long, is not suitable for this estimation. One aim of this paper is to explore this distribution by means of Lagrangian simulations.

There are different classes of microplastic particles denser than seawater. For example, dense synthetic microfibers have been found to strongly dominate in sediment samples far from the coast (Woodall et al.2014; Fischer et al.2015; Bergmann et al.2017; Martin et al.2017; Peng et al.2018) and have been detected in large proportions in deep-water samples and sediment traps in the open sea as well (Bagaev et al.2017; Kanhai et al.2018; Peng et al.2018; Reineccius et al.2020). With these mostly originating from land-based sources (Dris et al.2016; Carr2017; Gago et al.2018; Wright et al.2020), it is not obvious to explain their abundance on abyssal oceanic plains (Kane and Clare2019). Maritime-activity sources (Gago et al.2018) can contribute to that. Another reason could be that their special and deformable shape results in a strongly reduced settling velocity (Bagaev et al.2017) that allows long-distance horizontal transport (Nooteboom et al.2020). In any case, it is difficult to estimate the amount of microfibers in the oceans due to sampling issues and to their absence from statistics of mismanaged plastic waste (Carr2017; Barrows et al.2018), and we will not consider them further in this paper. We also disregard films, which are only sporadically encountered in the open ocean (Bagaev et al.2017) and thus have moderate importance.

We concentrate in the following on dense rigid microplastic particles. The most abundant particles of this class are fragments (e.g., Martin et al.2017; Peng et al.2018), which have an irregular shape, but their extension is usually comparable in the three dimensions. Experimental estimates for the settling velocities of irregular fragments or other nonspherical particles have suggested considerable deviations from values predicted by the Stokes law (Kowalski et al.2016; Khatmullina and Isachenko2017; Kaiser et al.2019), so that it is unclear how a precise full equation of motion should be constructed. For a qualitative exploration of particle transport through the water column, we will argue in Sect. 4.1 and Appendix A that the Maxey–Riley–Gatignol (MRG) equation (Maxey and Riley1983) should be appropriate for a reasonably wide range of such particles.

Whatever their precise equation of motion is, these sinking particles (directly detected by Bagaev et al.2017, and Peng et al.2018) are thought to reach the seafloor relatively fast (Chubarenko et al.2016; Kane and Clare2019; Soto-Navarro et al.2020), landing within horizontal distances of tens of kilometers from their surface location of release (see Sect. 4.2 and Appendix B). One consequence of their fast sinking is the absence of almost any fragmentation after they leave the sea surface (Andrady2015; Corcoran2015), and the influence of biological processes on the particles' properties should also be moderate, leaving their size and shape intact during sinking. Note that, in contrast to the case of floating plastics (Kooi et al.2017; Kvale et al.2020), interaction of sinking plastics with particulate matter of biological origin appears to be moderate. This is according to the absence of a need to disassemble microplastic pieces from biological aggregates during sample processing as described by Bagaev et al. (2017). Note, however, that experimental results by Michels et al. (2018) indicate that aggregation with organic material might occur within a sufficiently short time at surface layers, which would likely lead to increased sinking velocities (Long et al.2015). Transport by bottom currents (Kane and Clare2019; Kane et al.2020) is important for explaining their distribution in sediments after coastal release. However, the statements above imply that the dense rigid microplastic content of samples from deep-sea trenches and abyssal plains (van Cauwenberghe et al.2013; Fischer et al.2015; Peng et al.2018; Kane and Clare2019) must originate from sources at the surface of the open sea rather than from coastal inputs.

While methodological issues make the quantification of abundance difficult (Song et al.2014; Andrady2015; Filella2015; Lindeque et al.2020), negatively buoyant microplastic fragments have indeed been found in surface and near-surface samples of the open waters of the Mediterranean Sea (Suaria et al.2016) and the Atlantic Ocean (Enders et al.2015), respectively, from which they can contribute to microplastic content of the water column and deep-sea sediments (Fischer et al.2015; Bagaev et al.2017). Horizontal transport of these particles can be carried out by marine organisms, and spontaneous attachment to pieces of positive buoyancy is a further possibility but is not yet discussed in the literature. Composite pieces of debris or those that contain trapped air (including foams in some cases) may also represent a source of microplastic ending at the water column (Andrady2015). However, most of such particles are presumably released by local maritime activity. An example of this are flakes of paint and structural material from boats and ships, which contain negatively buoyant alkyds and poly(acrylate/styrene). Despite the particle's high density, large amounts of them may be found in the sea surface microlayer where surface tension keeps them floating (Song et al.2014). The range of horizontal transport of these particles at the sea surface is unclear but expected to be restricted to short distances because sinking from the sea surface microlayer is considerable, especially in waters disturbed by waves (Hardy1982; Stolle et al.2010).

While the idea of Kooi and Koelmans (2019) to treat all plastic particles together by means of continuous distributions is appealing, the above considerations strongly favor the separate treatment of positively buoyant pieces, negatively buoyant microfibers and negatively buoyant rigid particles of sufficiently big size, since these classes have very different dynamics and sources. In the following we concentrate on the properties, amount and dynamics of particles of the last class.

3 Considerations for modeling negatively buoyant rigid microplastics

3.1 Physical properties

From a meta-analysis of 39 previous studies, Erni-Cassola et al. (2019) established the proportion of the most abundant polymer types discharged into water bodies: PE (polyethylene, 23 %), PP (polypropylene, 13 %), PS (polystyrene, 4 %) and PP&A (group of polymer types formed by polyesters, PEST; polyamide, PA; and acrylics, 13 %). Note that these proportions do not distinguish between different regions (e.g., coastal region or open water; even inland water bodies of urban environments are included in the analysis) and between the particle types (size range and shape) concerned in the different studies. We organize these polymer types according to their density (Chubarenko et al.2016; Andrady2017; Erni-Cassola et al.2019) in Fig. 1: PP between 850–920 kg/m3, PE of 890–980 kg/m3, PS of 1040 kg/m3 (excluding its foamed version), PEST in the range 1100–1400 kg/m3, PA within 1120–1150 kg/m3 and acrylic of 1180 kg/m3. There is also some less abundant plastic like polytetrafluoroethylene (PTFE) which has higher densities, in the range 2100–2300 kg/m3.

Thus, the full range of microplastic particle densities in the ocean, denoted here as ρp, is 850–2300 kg/m3, and most of them have densities within the interval 850–1400 kg/m3. This has to be compared with the seawater density, which close to the surface has a conventional mean value of ρf=1025 km/m3 (red line in Fig. 1) and changes around 1 % from the surface to the sea bottom. Since we are interested in sinking material, and for the sake of maximal practicality, we restrict our study to microplastics of densities 1025 kg/m3<ρf<1400 kg/m3.

Figure 1Polymer densities for the most abundant microplastics identified in water bodies (Erni-Cassola et al.2019).


Another relevant property of plastic particles is their size. By a widely accepted definition, microplastics are particles with a diameter less than 5 mm without any lower limit (Arthur et al.2009). Some observations at the ocean surface show that the most common diameter is around 1 mm (Cózar et al.2014, 2015), with an exponential decay with increasing diameter up to 100 mm. However, the absence of this peak in other studies that show an increasing abundance with further decreasing size (Enders et al.2015; Suaria et al.2016; Erni-Cassola et al.2017) suggests (Song et al.2014; Andrady2015; Erni-Cassola et al.2017; Bond et al.2018; Lindeque et al.2020) the need for new technologies in sampling methods (which usually use trawl nets with a mesh size around 0.3 mm) and especially for the adaptation of careful and standardized analysis procedures to avoid artifacts (Filella2015).

Field data about distributions of size and quantifiers of shape for negatively buoyant rigid particles in the water column or deep-sea sediments are not available to date to the best of our knowledge, except in the Arctic from Bergmann et al. (2017). However, their results may not apply to the majority of the oceans because of the very special dynamics provided by melting and freezing of sea ice (Bergmann et al.2017). Data from Bergmann et al. (2017) and Song et al. (2014) about unspecified sedimented fragments and paint particles, respectively, exhibit an increasing abundance with decreasing size, most particles being smaller than 0.05 mm. Laboratory findings about surface degradation of individual particles also indicate such a tendency (Song et al.2017). Thus, these findings seem to indicate the prominent presence of small pieces of plastic. Nevertheless the observations of Bagaev et al. (2017), Kanhai et al. (2018) and Peng et al. (2018) do not indicate this abundance of small particles.

For these reasons, we will disregard particles of extremely small size. To keep our qualitative study sufficiently simple, we will consider all our modeled particles to have a radius a=0.05 mm (a diameter of 0.1 mm). This is a rather small size, but still within the commonly measured ranges. As we will discuss in Sect. 4.1, this radius is well within the validity range of the MRG equation.

3.2 Source estimation

In this subsection we indicate the total amount of dense microplastics entering the water column in open waters of the Mediterranean. Despite the correlation of plastic source with coastal population density, the rapid fragmentation of small particles along the shoreline (Pedrotti et al.2016) and the seasonal variability of spatial distribution of floating particles (Macias et al.2019), we focus on local maritime activity and exclude direct release from surface accumulation areas or the coast, either from urban areas or from rivers. The estimations are based on the results of Kaandorp et al. (2020). They provide a total amount of yearly plastic release into the Mediterranean in the range 2200–4000 t, from which around 37 % corresponds to negatively buoyant plastic, and 6 % are due to maritime activity. This 37 % agrees well with previous global estimations (Lebreton et al.2018).

We will take these numbers, 4000 t per year, 37 % of sinking particles and the proportion of direct release by maritime activity (6 %), to obtain in Sect. 4.3 an estimate for the basin-wide yearly release of negatively buoyant sphere-like microplastics in the open Mediterranean. Note that we choose the upper bound, 4000 t per year, in order to account for the considerable amount of unregistered particles.

3.3 Dynamics

A standard modeling approach (Siegel and Deuser1997; Monroy et al.2017; Liu et al.2018; Monroy et al.2019) for the transport of noninteracting sinking particles is to consider the time-dependent particle velocity v as the combination of the ambient fluid flow u and a settling velocity vs as

(1) v = u + v s ,


(2) v s = ( 1 - β ) g τ p , β = 3 ρ f 2 ρ p + ρ f     and     τ p = a 2 3 β ν .

g denotes the gravitational acceleration vector, pointing downwards; β is a parameter depending on the particle and the fluid densities, ρp and ρf, respectively. Particles heavier than water have β<1, and β=1 for neutrally buoyant particles. The expression given for β assumes spherical particles. τp is the Stokes time, i.e., the characteristic response time of the particle to changes in the flow, where a is the radius of the particle and ν the kinematic viscosity of the fluid. Although Eq. (1) is commonly used, we are not aware of a systematic justification of it in the microplastics context. This will be done in Sect. 4.1.

3.4 Numerical procedures

For the flow velocity u we use a 3D velocity field from NEMO (Nucleus for European Modelling of the Ocean), which implements a horizontal resolution of 1/12 and 75 s levels in the vertical with data updated every 5 d (Madec2008; Madec and Imbard1996). Salinity and temperature are also extracted from that model. The Parcels Lagrangian framework (Delandmeter and van Sebille2019) is used to integrate the particle trajectories from Eq. (1) or more complex ones to be considered in Sect. 4.1. Typical numerical experiments to obtain the results presented below consist of distributing a large number N of particles in a horizontal layer over the whole Mediterranean on the nodes of a sinusoidal-projection grid (Seong et al.2002), so that their release is with uniform horizontal density. We locate this input source at 1 m depth to avoid surface boundary conditions. After particles are released at some initial date, in a so-called flash release, they evolve under equations of motion such as Eq. (1), and the statistics of the resulting particle cloud are analyzed.

4 Results

4.1 Range of validity of Eq. (1)

We next show, closely following the treatment of Monroy et al. (2017) for biogenic particles, that possible inertial effects that would correct Eq. (1) are negligible for the sizes and densities of typical dense microplastics. To this end, similarly to many other studies (Michaelides2003; Balkovsky et al.2001; Cartwright et al.2010; Haller and Sapsis2008), we start by choosing the simplified standard form of the more fundamental Maxey–Riley–Gatignol (MRG) equation (Maxey and Riley1983) and analyze under which conditions it is valid for microplastic transport. After finding the MRG equation to be valid for an important range of microplastic particles, we will explore its relationship with Eq. (1).

The simplified MRG equation gives the velocity v(t) of a very small spherical particle in the presence of an external flow u(t) as

(3) d v d t = β D u D t + u - v + v s τ p .

Beyond sphericity, two conditions are needed for the validity of Eq. (3) (Monroy et al.2017; Maxey and Riley1983): (a) the particle radius, a, has to be much smaller than the Kolmogorov length scale η of the flow, which has values in the range 0.3 mm <η<2 mm for wind-driven turbulence in the upper ocean (Jiménez1997); (b) the particle Reynolds number Rep=a|v-u|νavsν should satisfy Rep≪1. Note that this last condition imposes restrictions on the values of the particles' density and size, partially via the settling velocity vs=|vs|. For the most abundant sinking microplastics, i.e., with densities ρp=1025–1400 kg/m3, we now determine the range of validity of Eq. (3) assuming ν=1.15×10-6 m2/s and ρf=1025 kg/m3 to be fixed. This gives β in the range 0.8–1. The possibility of small changes in the seawater density as the particle sinks, which translates to variations in vs, will also be analyzed in Sect. 4.2.

In Fig. 2 we show a diagram with the settling velocities and particle sizes for which Eq. (3) is valid. We plot the minimal value of the Kolmogorov scale η=0.3 mm with the red line (Jiménez1997), and Rep=1 with a black line, which bound the area of validity (shaded in the plot). We also indicate vs as a function of a for β=0.8 with the blue curve, corresponding to the upper bound to vs for typical microplastic densities. In total, the zone with light shading in Fig. 2 represents a parameter region where Eq. (3) applies for particles with β<0.8 (i.e., particles falling faster than the typical ones), whereas the area of our interest, corresponding to β≥0.8, is represented by dark shading, denoting the typical plastic sizes and corresponding settling velocities for which the equation is valid. As a rule of thumb, in a typical situation, validity of Eq. (3) requires vs<0.01 m/s and a<0.3 mm. As discussed in Sect. 3.1, information about particles in the validity range is particularly sparse for surface waters because of the usual sampling techniques, but sediment data indicate the prevalence of sufficiently small particles. Furthermore, in sufficiently calm waters, the Kolmogorov scale is larger (of the order of millimeters; Jiménez1997), so that a can be increased to this size without compromising the equation validity. In any case, these estimates of the Kolmogorov scale assume wind-driven turbulence and are thus restricted to the mixed layer (Jiménez1997), below which η is undoubtedly larger. Deviations from a spherical shape may lead to a more complicated motion than that described by the MRG equation, especially through particle rotation (Voth and Soldati2017). In Appendix A, we present quantitative arguments for the applicability of the MRG equation to rigid microplastic particles of common shapes in the parameter ranges of our interest.

The simplified MRG equation, Eq. (3) thus represents an appropriate basis for qualitative estimations of the transport properties of negatively buoyant rigid microplastics in the water column. Note that rigidity of the particles is an essential condition which is why the advection of microfibers is out of the scope of this paper.

Figure 2Settling velocities and particle sizes for which Eq. (3) holds. Kolmogorov scale is represented by the red line and Rep=1 with a black line, which bound the area of validity. Blue curve corresponds to vs=vs(β,a) for β=0.8, the upper bound to vs for typical microplastic densities. Dark shading denotes the plastic particle sizes and corresponding settling velocities for which application of Eq. (3) is valid.


The connection between Eq. (3) and its approximation Eq. (1) is made by noticing that τη≈1 s in the ocean (Monroy et al.2017; Jiménez1997), so that the Stokes number St=τp/τη, which measures the importance of particle inertia in a turbulent flow, is very small (of the order of 10−3–10−2). Thus an expansion of the MRG equation for small St (smallness of the Froude number, i.e., smallness of fluid accelerations with respect to gravity, is also required) can be performed. The expansion in its simplest form leads to (Balachandar and Eaton2010; Monroy et al.2017; Drótos et al.2019)

(4) v u + v s + τ p ( β - 1 ) D u D t .

We can now take the results of Monroy et al. (2017) for biogenic particles of sizes and densities similar to the microplastics considered here to show that the inertial corrections (the term proportional to τp) in Eq. (4) are negligible, so that the simpler Eq. (1) correctly describes sinking of microplastics in the considered parameter range. For completeness, we report in Appendix B the explicit numerical calculations showing this (in which the influence of the Coriolis force is also taken into account, since it is known to be of the same order as or larger than the inertial term when a large-scale flow is used for u). In particular we find, from experiments of releasing particles with β in the range 0.8–1 from 1 m below the surface in the whole Mediterranean, that differences in horizontal particle positions between the dynamics of Eq. (4) and Eq. (1) after 10 d of integration are just 0.26 % of the horizontal displacements. For the vertical motions the difference is of about 0.05 %. Thus, Eq. (1) provides a proper description of the dynamics.

Even if an equation of motion is accurate, the accuracy of its solution is limited by that of the input data. In particular, small-scale flow features are absent from oceanic velocity fields u simulated on large-scale domains, which is an important limitation of the respective solutions of Eq. (1). The NEMO velocity field of our choice is not an exception, but a rigorous evaluation of the corresponding errors of particle trajectories is not possible without knowing the actual small-scale flow. Nevertheless, one can roughly estimate the effect of these small scales by adding a stochastic term to Eq. (1) with statistical properties similar to the expected ones for a small-scale flow (Monroy et al.2017; Kaandorp et al.2020). Results similar to those by Monroy et al. (2017), summarized in Appendix B, indicate that after 10 d of integration the relative difference between particle positions given by Eq. (4) with and without this “noise” term modeling small scales (using β=0.99) is around 8 % for the horizontal displacements and 5 % for the vertical ones. The figures become 12 % and 5 %, respectively, when evolving the particles for 20 d. These errors are moderate, although they may be of importance under some circumstances (Nooteboom et al.2020). We consider these figures as a baseline to evaluate corrections to the simple Eq. (1): adding more complex particle-dynamics terms to it will not improve plastic-sedimentation modeling unless the effect of these corrections is significantly larger than the above estimations for the effect of the unknown small-scale flow. In the following we consider the simple Eq. (1), but we estimate the implications of assuming or not assuming a constant value of the water density.

4.2 Effect of variable seawater density

In this section we analyze the role of a variable seawater density on the particle settling dynamics. Fluid density is calculated from the TEOS-10 equations, which is a thermodynamically consistent description of seawater properties derived from a Gibbs function, for which absolute salinity is used to describe salinity of seawater and conservative temperature replaces potential temperature (Pawlowicz2010). In the simulations described in this section, as particles move in the ocean they encounter different temperatures and salinities, as given by the NEMO model described in Sect. 3.4, and then they experience different values of the ambient-fluid density.

We consider particles of a fixed density ρp=1041.5 kg/m3. This implies that for a nominal water density of ρf=1025 kg/m3 the value of β would be β=0.99, giving a sinking velocity vs=6.2 m/d for our particles of radius a=0.05 mm, but this sinking velocity will be increased or decreased in places where water density is lower or higher, respectively, so that we have a spatially and temporally dependent velocity in Eq. (1). The particle density and size have been chosen to be representative of the slowly sinking microplastic particles, for which we expect the seawater density variations to have the largest impact. In this way we find some upper bound for the importance of variability in seawater density for particle trajectories.

We release N=78 803 particles over the whole Mediterranean Sea and monitor their trajectories under Eq. (1). Figure 3a shows the histogram of water densities encountered by the particles when the release is performed on 8 July 2000. On this summer date, the Mediterranean is well stratified, at least in its upper layers. Initially the particles are in surface waters with a range of salinities that average approximately to the nominal ρf=1025 kg/m3. But as they sink in time they reach layers with higher densities (and more homogeneous across the Mediterranean). When the release is done in winter (Fig. 3b) the water column is more mixed, so that the range of water densities experienced by the particles released at different points is always narrow. But the mean water density turns out to be always larger than the conventional surface density of ρf=1025 kg/m3, so that a slightly slower sinking is expected to occur.

Figure 3Normalized histogram of the seawater density ρf at the positions of N=78 803 particles after t=0, 10, and 20 d of being released over the whole Mediterranean. (a) Summer release (release date 8 July 2000). (b) Winter release (release date 8 January 2000). The particles' density is fixed at ρp=1041.5 kg/m3, and fluid density is obtained from the TEOS-10 equations. The vertical line indicates a conventional seawater density of 1025 kg/m3.


Figure 4The distance, as a function of time, between trajectories obtained with constant nominal fluid density of ρf=1025 kg/m3 and the actual variable fluid density, both starting at the same initial location. The range of the values among all particles released in different points of the Mediterranean is indicated by the shaded area, while the solid line indicates the average over the particles. Particles have ρp=1041.5 kg/m3, and all parameters are the same as for the summer release in Fig. 3. (a) Horizontal distances. (b) Vertical distances.


We illustrate the impact of this variable density on particle trajectories for the summer release in Fig. 4. Here we compute, as a function of time, the range of horizontal |x(0)-x(1)| and vertical |z(0)-z(1)| distances and its average among particles. Trajectories x(0)(t) and x(1)(t) are obtained with constant nominal fluid density (1025 kg/m3) and position-dependent fluid density, respectively, using the same release location and date, 8 July 2000, in both cases. Particle density is fixed at ρp=1041.5 kg/m3. The difference between the two calculations (and thus the error of considering that constant value for the density) should be compared to average horizontal and vertical displacements of 95 km and 124 m, respectively, at t=20 d. At that time, we thus find that the influence of variable fluid density on the dynamics is about 3 % for the horizontal movement and 6 % for the vertical displacement on average.

Table 1 Relative effect on horizontal and vertical particle positions after 10 and 20 d of integration, averaged over 78 803 particles released over the whole Mediterranean at 1 m depth, of replacing the actual seawater density by a nominal value ρf=1025 kg/m3.

Download Print Version | Download XLSX

A summary of the average relative differences on horizontal and vertical particle positions between using the location-dependent seawater density and a nominal constant value ρf=1025 kg/m3, both in winter and summer periods, is displayed in Table 1. The relative error produced by assuming a constant density is larger in the vertical direction. It is also larger for the release in winter, but this is a consequence of taking a value for the reference density that is not representative of winter waters but is strongly biased (see Fig. 3, right). If using a reference value more appropriate for winter waters (say ρf≈1027 kg/m3) the relative error remains quite small, due to the weaker stratification of the sea during this season. In fact, the reference value is also biased in the summer unless the investigation is restricted to the surface.

In brief, we see that the effect of location-dependent density may be a relevant effect to evaluate microplastic transport. At least, the traditional value of seawater density may be biased, which may be reflected in the particle trajectories. We recall, however, that we used parameters for the particle properties for which they are slowly falling. The impact of variable density on particles that sink faster will be smaller. Also, the effects reported in Table 1 remain of the order of the estimations of the effects of unresolved small scales of the flow (Sect. 4.1). As a consequence, in the following we will not consider variable seawater density, but restrict our modeling to Eq. (1) with a constant nominal value of the sinking velocity vs.

4.3 Total mass and vertical distribution of microplastics

We will first estimate the total mass of negatively buoyant rigid microplastics in the water column of the open Mediterranean Sea by assuming a uniform vertical distribution; then we will justify this assumption by running numerical simulations according to the conclusion of Sect. 4.1 regarding the equation of motion.

For estimating the total mass, we take the quantities of Sect. 3.2 (4000 t/year of plastic release, with 37 % being negatively buoyant and of which 6 % originates from maritime activities) to compute the rate r at which microplastic particles of our interest enter the water column in the open sea: r=4000 t/year ×0.37×0.06=88.8 t/year, or r=0.24 t/d.

The next step is to estimate the time during which these microplastic particles remain in the water column before reaching the sea bottom. We take the mean depth for the Mediterranean to be h=1480 m (Eakins and Sharman2010; GEBCO Compilation Group2020) and estimate a residence time τ as the time of sinking to that mean depth. The residence time depends on the sinking velocity, τ=h/vs, and thus on the physical properties of the microplastic particles. Assuming a seawater density ρf=1025 kg/m3, and the range of plastic densities and their proportions described in Sect. 3.1, we see from Eq. (2) that for microplastic particles of radius a≈0.05 mm the range of sinking velocities is 6.20–509.23 m/d, giving a residence time in the range 3.1–255 d. Averaging these times weighted by the proportion of each type of plastic we get τ14 d. Combining the input rate r with this mean residence time we get an estimate for the total amount present in the water column at any given time as Q=rτ: the result is Q≈3.36 t of dense rigid microplastics if all of them would be in the form of particles of size a=0.05 mm. This is below but close to 1 % of the estimated upper bound of 470 t of floating plastic in the Mediterranean (according to the corresponding estimation of Kaandorp et al.2020).

We emphasize the many uncertainties affecting this result (Sect. 3.1 and 3.2), and we highlight the one related to particle size: because of the quadratic dependence of the sinking velocity on the particle radius a, Eq. (2), choosing the particle size to be half of the one used here will lead to an estimate 4 times larger for the mass if the same release rate is assumed. This enhanced retention of smaller particles in the water column may imply, depending on the actual size distribution, a dominance of very small particles on the plastic mass content of the water column. However, our estimates of plastic input into the ocean (we use mainly Kaandorp et al.2020) rely on observations that do not catch extremely small particles. These considerations further justify our choice of a radius a=0.05 mm, small but still easily detectable, as convenient to provide reasonable estimations of negatively buoyant rigid microplastic mass in the water column within commonly quoted size ranges. We can not exclude larger plastic content at smaller sizes. Another source of bias may be that in this study the impact of small-scale turbulence and convective mixing events is not considered. While small-scale turbulence might cause an increase in lifetimes of particles in the water column, dense water formation and rapid convection, a process reported in areas such as the Gulf of Lion, would likely reduce particle retention time. These events take place in winter and were shown to transfer particles from the ocean surface to mid-waters (1000 m) and deep ocean (>2000 m) in a very short time (1–2 d) and lead to the formation of bottom nepheloid layers (de Madron et al.1999; Vidal et al.2009; Heussner et al.2006; Stabholz et al.2013).

The result for the total mass is independent of the horizontal distribution of particle release, which is quite inhomogeneous (Fig. 1 of Liubartseva et al.2018). However, for a rough estimate of the density of these microplastics in the water column, we assume a uniform particle distribution over the whole Mediterranean both in the horizontal and in the vertical. Since the volume of the Mediterranean is about 4.39×106 km3 (Eakins and Sharman2010) the estimated density would be ρV7.7×10-11 kg/m3 (with the above-discussed scaling issues with a). We remind the reader that this is a value for the open sea, and our study does not address coastal areas, where the density would likely be higher.

The above estimates are rather rough as a result of the uncertainties mentioned. The assumption of a uniform distribution in the vertical direction has not yet been justified either, but we will show it to be appropriate by means of our simulations of particle release starting at 1 m depth over the whole Mediterranean. Instead of performing a continuous release of new particles at each time step, and computing statistics over this growing number of sinking particles, we approximate this by the statistics of all positions at all time steps of a set of particles deployed in a single release event. This assumes a time-independent fluid flow, but this approximation is appropriate, since the dispersion of an ensemble of particles released in a single event follows rather well-defined statistical laws (see Sect. 4.4) and is thus independent of the time-varying details of the flow. Particles are removed when touching the bottom. For our estimate, we use β=0.8, i.e., assuming the fastest sinking velocity of typical plastic particles, for which particles reach vertical depths deeper if compared with the slower sinking velocities used in our study.

Figure 5 shows ρ(z), the density of plastic particles per unit of depth z in the whole Mediterranean, and also A(z), the amount of area that the Mediterranean has at each depth z. Both functions have been normalized such that the value of their integrals with respect to z is one; in this way the functions can be displayed in the same plot. We see that both curves are nearly identical (in fact, just proportional, because of the normalization), indicating that the variation of the number of particles with depth is essentially due to the decrease in sea area with depth. A clearer way to see that is to plot ρ(z)/A(z), which is proportional to the mean plastic concentration per unit volume at each depth z, ρV(z). We see that this quantity is nearly constant, at least in the first 3000 m. At larger depths a weak increase seems to occur, but this is made unclear by the poor statistics arising from the small area and number of particles present at these depths. Thus, the hypothesis of a uniform distribution of plastic in the water column seems to be a reasonable description of the simulation of the fastest-sinking particles.

A uniform distribution of plastics in z is what is expected if the vertical velocity of the particles is exactly constant (since each particle will spend exactly the same amount of time at each depth interval). The equation of motion used (Eq. 1) corrects this constant sinking velocity vs with a contribution u from the ambient flow. Thus, the close-to-constant character of the plastic concentration may imply that the flow correction u is negligible, at least when considering its effect over the whole Mediterranean. Another possibility is that the fluctuating flow component u in Eq. (1) results in a vertical dispersion compatible with a constant concentration. Although the former explanation predicts an alteration from a constant if the settling velocity is sufficiently small to allow u to induce a stronger vertical dispersion, we will see in the next section that a nearly constant concentration may be assumed for the majority of our parameter range.

Figure 5The continuous line is the area that the Mediterranean has at each depth z. The dashed line is microplastic density per unit of depth ρ(z) under continuous release of particles with β=0.8 at 1 m depth. Both curves have been normalized to have unit area, so that they can be compared on the same scale. The binning size is 100 m. The inset shows the ratio ρ(z)/A(z), proportional to the mass density of microplastic per unit of volume ρV(z).


4.4 Transient evolution

We now analyze in detail the transient evolution of particle clouds initialized by flash releases at a fixed depth. Numerically we proceed by releasing N=78 803 particles uniformly distributed over the entire Mediterranean surface at 1 m depth in the winter season, as already described. They evolve according to Eq. (1) using a constant water density. We take three examples for the particle density, which correspond to β=0.8, 0.9 and 0.99, or vs=153.48, 68.21 and 6.20 m/d for our particles of radius a=0.05 mm, respectively; in what follows, these setups shall be denominated as v153, v68 and v6. The horizontal displacements during the particle sinking times are much larger (of the order of 60 km) than the sea depth, so that in fact the particles are sinking sideways (Siegel and Deuser1997). However, the horizontal displacements still remain very small compared to the basin size, and we concentrate on the vertical motion. Even though the vertical steady distribution has been found to be close to uniform in Sect. 4.3, the reason for this is not evident, and we will give support here for the pertinence of this finding to most of the relevant parameter range.

Figure 6 shows the vertical particle distribution at different times (upper plot is for v153, middle for v68 and bottom for v6). The plot is given in terms of a rescaled variable z̃=z-tvsσz where σz2(zi-zi)2 is the variance of the particles' z coordinate. Here the subindex i refers to the particle and 〈…〉 denotes averaging over different particles. Thus, we plot in the figure the rescaled distribution of the particles around the average depth of the particles at any given time. For comparison, the normal distribution is plotted with dashed lines. This figure shows deviations from Gaussianity for early times. The deviation from normal distribution decreases for later instants but remains considerable, especially for the tails, which may also be indicative of anomalous diffusive behavior. For reference, particles reach the mean Mediterranean depth, h=1480 m, at times τ=9.64, 21.8 and 246.7 d for v153, v68 and v6, respectively.

Figure 6The probability density function, estimated from a histogram of bin size 0.1, of all particles released in the Mediterranean in the rescaled variable z̃=z-tvsσz for the different setups (v153, v68 and v6) and times (in days) as indicated. For comparison, normal distributions of zero mean and unit variance are shown with a dashed line.


Since a non-Gaussian distribution is usually linked to anomalous dispersion (Neufeld and Hernández-García2009), we now analyze this aspect in detail by considering how the variance of the vertical particle distribution, σz2(t), evolves. Although there is a continual loss of particles because of reaching the seafloor with a varied topography, we illustrate in Appendix C that our conclusions are likely unaffected by this effect.

According to Fig. 7, dispersion appears to be governed by different laws in different regimes, which we shall distinguish by the approximate effective exponents ν, defined through approximate behaviors σz2tν in different time intervals.

Figure 7Variance of depth reached by the particles as a function of time. Straight lines represent power laws for reference, with exponents 1 (in green, corresponding to standard diffusion) and 2 (in purple, corresponding to ballistic dispersion).


We start our analysis with the fastest-sinking particles (v153; Fig. 7a). At the very beginning, superdiffusion takes place with ν>2, which may be related to autocorrelation in the flow, but we will iterate on this question when comparing different settling velocities. Around t=1 d, the evolution seems to become consistent with normal diffusion (ν=1), usual after initial transients in oceanic turbulence (Berloff and McWilliams2002; Reynolds2002). However, around t=4.5 d, we can observe a crossover to ballistic dispersion (ν=2).

We explain this last crossover as resulting from a different mean sinking velocity in diverse regions of the Mediterranean, associated with up- and down-welling. This can be modeled in an effective way by writing the vertical position of particle i as

(5) z i = z i + ω ¯ i t + W i ,

where 〈…〉 denotes, as before, an averaging over different particles. Here we are assuming that zi−〈zi evolves according to the sum of a constant average velocity contribution ω¯i for sufficiently long times (a characteristic of the flow region traversed by particle i), and of Wi, a Wiener process representing fluctuations with zero mean and defining a diffusion coefficient Di for each trajectory by Wi2=Dit. The overbar refers to temporal averaging for asymptotically long times along the trajectory of a given particle (but assuming that the particle remains in a region with a well-defined ω¯i0), and Di characterizes the strength of the fluctuations. Assuming ω¯iWi=0,

(6) σ z 2 ( z i - z i ) 2 = ω ¯ i 2 t 2 + D i t ,

that is, the variance is a sum of a ballistic and a normal diffusive term, associated with regional differences in the mean velocity and with fluctuations, and dominating for long and short times, respectively. Writing D=〈Di, the crossover between the two regimes is obtained by equating the two terms as

(7) t * = D ω ¯ i 2 .

To evaluate Eq. (7), we first estimate ω¯i for each particle from the “asymptotically” long time of t=10.83 d, which is the latest time after the crossover still in the ballistic regime in Fig. 7a, when the contribution of fluctuations should already have become negligible. The horizontal pattern of the estimated ω¯i is presented in Fig. 8, which confirms its patchiness throughout the Mediterranean, associated with mesoscale features. Computing ω¯i2 and fitting a line to σz2(t) between t=1.4 and 4 d to estimate D, we obtain t*4.5 d from Eq. (7), which remarkably agrees with Fig. 7a. After approximately t=12 d there is hardly any dispersion, since most of the particles are close to the sea bottom (cf. Fig. 9) where the vertical fluid velocity is nearly zero. Note also a small drop in σz2 at the very end of the time series, where the results may actually be subject to artifacts; see Appendix C. However, this is of minor importance, since the distribution of particles so close to the bottom should nevertheless be strongly influenced by resuspension and remixing by bottom currents (Kane et al.2020).

Figure 8ω¯i as estimated from t=10.83 d plotted at the initial position of each particle i in the v153 simulation. The black rectangle in the Western Mediterranean is the area of large depth considered in Appendix C.


The different regimes are not as clear in the v68 case as for v153; see Fig. 7b. One evident novel feature is a subdiffusive regime during the transient from the initial superdiffusion (as in the case of horizontal tracer dispersion in the ocean studied by Berloff and McWilliams (2002) and Reynolds (2002)). Approximate normal diffusion is then observed until t=10 d, when a crossover to a faster dispersion does seem to take place; see the inset. A fit of normal diffusion from t=4 to 8 d and the velocity variance at t=12.5 d give an estimate t*11.7. However, the long-term ballistic regime is not clear. In fact, a long-term return from such a ballistic regime to normal diffusion is expected as a result of increasing horizontal mixing, which renders ω¯i time dependent and makes it approach zero. According to a careful visual inspection of the inset in Fig. 7b, this may take place already around t=14 d.

For v6 (Fig. 7c), the transition from the initial superdiffusive regime to that of normal diffusion appears to not involve subdiffusion. This is already informative: the fluid velocity field is the same for the three simulations with different settling velocity, so the differences must originate from the different rates of sampling of the different fluid layers by particles while they sink. In particular, the decay of autocorrelation is obviously faster for faster-sinking particles, since it is determined by the spatial structure of the velocity field.

While this is one possible explanation for the earlier timing of the initial transition from anomalous to normal diffusion for higher settling velocity, one cannot exclude that a depth-dependent organization of the flow is more in play; note that ν>2 at the beginning, which might not be explained by simple autocorrelation but might be characteristic of properties of the velocity field at those depths. The governing role of the spatial structure is supported by Fig. 9: the transition in question takes place at the same depth (≈100 m) in the different simulations, which seems to point to mixed-layer processes. Depth dependence might also govern the suppression of ballistic dispersion for long times, but it is very unclear.

Figure 9Variance of depth reached by the particles as a function of their mean depth.


Note in Fig. 9 that the vertical variance is not expected to grow much larger for v68 than for v153 even if the simulation were longer. Therefore, even if the constancy of the steady vertical distribution relies on the weak vertical dispersion for v153 (see Sect. 4.3), constancy is expected to hold in most of our parameter range. A considerably stronger dispersion and a possible corresponding deviation from constancy may arise only for extremely low settling velocities, like for v6 in Fig. 9.

5 Conclusions

We have discussed the different types of plastics occurring in the water column, pointing out gaps in our knowledge about the sources, transport pathways and properties of such particles. It would be highly beneficial to have distributions of size, polymer type and quantifiers of shape recorded separately for the dynamically different classes of microplastics.

We have focused our attention on rigid microplastic particles with negative buoyancy. We have argued that the simplified MRG equation approximates the dynamics of such particles sufficiently well for qualitative estimations.

We have then analyzed the importance of different effects in this equation and concluded that the Coriolis and the inertial terms are negligible. When a velocity field of large-scale nature is input to the equation (such that small-scale turbulence is not resolved), or when the variability in seawater density is neglected, moderate but possibly non-negligible errors emerge (Nooteboom et al.2020). However, our conclusions about the vertical distribution and dispersion of microplastics rely on robust features of the large-scale flow and must remain unaffected by moderate errors. We also note that the traditional value of seawater density, ρf=1025 kg/m3, is representative only for near-surface layers in the summer, and correcting for the bias could reduce the error of simulations with a constant seawater density. A suitable equation of motion for the particles considered is constructed by adding to the external velocity field a constant settling term, as also found by Monroy et al. (2017) for marine biogenic particles.

When the velocity field of the Mediterranean Sea is approximated by realistic simulation, this equation of motion results in a nearly uniform steady distribution along the water column, except perhaps at extremely low settling velocities. The corresponding total amount of plastic present in the water column is relatively small, close to 1 % of the floating plastic mass, but it may be an important contribution to the microplastic pollution in deep layers of the ocean and is subject to several uncertainties.

Note that only those microplastic particles are considered here that have not yet sedimented on the bottom, and the plastic amount sedimented on the seafloor is large (Fischer et al.2015; Liubartseva et al.2018; Peng et al.2018; Mountford and Morales Maqueda2019; Soto-Navarro et al.2020). The suitability of our equation of motion to describe the sinking of a class of microplastic particles implies that advection by the flow may contribute to large-scale horizontal inhomogeneity of deep-sea plastic sediments by means of recently described noninertial mechanisms (Drótos et al.2019; Monroy et al.2019; Sozza et al.2020). This may be especially so in regions where redistribution by bottom flows is restricted to small distances, like abyssal plains (Kane and Clare2019). Resuspension and redistribution may be dominant in forming sedimented patterns (Kane et al.2020), and a future investigation should take all processes into account to identify zones of high plastic concentration on the sea bottom.

As for the vertical distribution profile, its approximate uniformity may be linked to the weak vertical dispersion of particles that is found in our simulations, started with a flash release over the whole surface of the Mediterranean sea. The shape of the emerging transient vertical distribution exhibits deviations from a Gaussian, which are related to anomalous diffusive laws that dominate the vertical dispersion process in some phases.

The different diffusive laws are related to the properties of the decay in the Lagrangian velocity autocorrelation defined along the trajectories of the sinking particles. An important example is the transition from initial superdiffusion to a longer phase of normal diffusion, occurring around 100 m depth, which indicates that the particles enter into a different flow regime. Another characteristic of the velocity field is a horizontal patchiness, which results in a long-term ballistic dispersion as long as the particles' horizontal displacements remain small. The vertical diffusion returns to the normal type when horizontal mixing becomes more developed. These results suggest regional differences in the sinking process, so that regional modeling might be more appropriate than a whole-basin approach. Future studies will include different areas of the oceans and analyze the role of Lagrangian coherent structures on the different vertical dispersion regimes.

Appendix A: Deviations from a spherical particle shape

We quantitatively assess the impact of deviations from a spherical shape through a correction to the settling velocity vs. The simplified MRG equation, Eq. (3), or its first-order approximations in the Stokes number, Eqs. (4) and (B1), are affected by particle geometry through the drag force and the added mass term; however, accelerations are irrelevant for vs, so that the added mass term does not appear in its formulation or in the simple approximation of Eq. (1). We will compare values of the settling velocity describing nonspherical and spherical particles with the same density, then finally we will comment on the results' relevance for Eqs. (3), (4) and (B1).

Most generally, the settling velocity vector vs can be obtained by balancing the drag force Fdrag(vu) (a function of the difference of the particle and the fluid velocities, v and u, respectively) with the resultant of gravitational and buoyancy forces:

(A1) 0 = F drag ( v - u ) + V ρ p - ρ f g ,

with v-u=vs, where V is the particle's volume, ρp and ρf are the densities of the particle and the fluid, respectively, and g is the gravitational acceleration vector. For a spherical particle with radius a, the Stokes drag force reads as

(A2) F drag (sph) ( v - u ) = - 6 π μ ( v - u ) a ,

where μ is the dynamical viscosity of the fluid. According to Leith (1987) and Ganser (1993), an appropriate approximation for small nonspherical particles is

(A3) F drag (non) ( v - u ) = - 6 π μ ( v - u ) 1 3 a n + 2 3 a s ,

where an is the radius of the sphere with equivalent area projected on the plane perpendicular to the relative velocity vu, and as is the radius of the sphere with equivalent total surface. From either of the last two equations, the settling velocity is obtained by substituting v-u=vs and solving Eq. (A1) for vs. We denote the magnitudes of the settling velocities obtained from Eq. (A2) and Eq. (A3) by vs(sph) and vs(non), respectively.

To characterize the correction in the settling velocity for a given nonspherical particle (with a given density ρp) with respect to assuming a spherical shape with a radius a, we will consider

(A4) q v s (non) v s (sph) = 3 4 π V (non) a 2 1 3 a n + 2 3 a s ,

where V(non) is the real volume of the given particle. In order to evaluate Eq. (A4), one has to specify the shape and the size of the particle, its orientation with respect to its relative velocity, and also how a is derived from its real size.

Note that it is always possible to define an a for which q=1, i.e., for which there is no correction arising from the deviation from a spherical shape. In this sense, any choice of a representing a spherical shape, including ours in this paper, describes the settling velocity of certain nonspherical particles; the question is just their shape and size, which will mutually depend on each other for a given a. We will nevertheless proceed by choosing a shape class and defining a along independent considerations, because we intend to link a given a to a single particle size identified during the processing of field observations.

The shape of rigid microplastic particles is not usually described in the literature, but we can see photographs of some examples in, for example, Song et al. (2014), Fischer et al. (2015), and Bagaev et al. (2017). For an explorative computation, a reasonable choice seems to be a rectangular cuboid with edges A, B=B^A<A and C=C^A<B<A, where one or both of B^ and C^ are less than 1 but greater than, say, 0.1.

Under this assumption, the particle size will correspond to the longest edge, A, of the cuboid if the size is identified through microscopy as the largest extension (“length”; e.g., Cózar et al.2015), and it may be related more to the middle edge, B, if one thinks of a sieving technique (e.g., Suaria et al.2016). The naive choice will be a=A/2 and a=B/2 in these two cases.

We can substitute either of these choices of a in Eq. (A4), as well as the appropriate formulae describing the actual cuboid. V=ABC is unique, and so is as,

(A5) a s = A B + A C + B C 2 π 1 2 .

However, an depends on the particle's orientation with respect to the relative velocity. Implications will be discussed when interpreting the results, and we take here all three directions parallel to edges A, B and C to represent different possibilities. The corresponding expressions for an read as

(A6) a n , X = A B C π X 1 2

for X=A, B and C.

After substituting all these expressions in Eq. (A4), we obtain


where X=X^A has been introduced.

Figure A1qA(A/2) as a function of B^ and C^ for B^,C^[0.1,1] with B^>C^. Dotted lines represent qA(A/2)=0.1.


We plot qA(A/2) in Fig. A1 as a function of B^ and C^ for B^,C^[0.1,1] with B^>C^. Its range extends from 0.07 to 1.5 on this domain, but it drops below 0.1 only for C^ very close to 0.1 and B^ below 0.2, i.e., for extremely thin rod-like particles, which do not appear to be common based on photographs. The range of qA(B/2) (not shown) on the same domain is between 0.2 and 7, and values above 4 are again restricted to very small C^ and to B^<0.3. The results are very similar for other choices of X; deviations beyond 20 % with respect to X=A are only found for small C^ and do not reach beyond 40 % even there.

We have left the question of which orientation is relevant open so far. In small-scale isotropic turbulence, which is certainly present in the ocean, nonspherical particles have a preferential alignment with certain characteristics of the flow but undergo rotation (Voth and Soldati2017). This is why we have chosen to simply cover three perpendicular orientations in our analysis and have found that differences that may arise from changes in the orientation are minor in most of the domain that describes shapes. The only relevant exception is small C^ with B^1. This regime may characterize paint flakes (Song et al.2014; Bagaev et al.2017), but the relative difference remains below 40 % even there.

Even though the real advection of the particles will become more complicated as a result of the ever-changing orientation and may thus be beyond the scope of the MRG equation (cf. the discussion in the main text about the settling velocity of irregular particles), we have found that changing orientation introduces minor variations in the value of the settling velocity. Together with the absence of order-of-magnitude corrections that may arise from a nonspherical shape (but comparing shapes under the assumption of the same particle density), this gives quantitative support for the applicability of a spherical shape in Eq. (1) of the main text.

Finally, we briefly comment on the more general Eqs. (4) and (B1), in which effects from nonsphericity arise in the inertial term through added mass. Since corrections in added mass with respect to a sphere are of order 1 for all common shapes (see Kaneko et al.2014, for an overview), we believe that the finding of Appendix B on the negligible importance of inertial effects in these equations is not affected by a deviation from sphericity. The Stokes number, which is proportional to the settling velocity and also depends on the coefficient of added mass (Stτp with τp given by Eq. 2), is estimated to be 10−3–10−2 for spherical particles in Sect. 4.1, so that it will not increase to 1 due to a nonspherical shape either, hence leaving the approximation (4) of Eq. (3) valid.

Appendix B: Importance of different physical effects in the dynamics

We present here the detailed numerical analysis of the relevance of a finite time of response (Stokes time, τp) of the particle to the fluid forces, the Coriolis force and scales unresolved by the NEMO velocity field.

We incorporate the first two effects to a single equation,

(B1) v = u + v s + τ p ( β - 1 ) D u D t + 2 Ω × u ,

which is identical to Eq. (4) except for the addition of the Coriolis force, 2Ω×u. Ω is Earth's angular velocity vector. We include the Coriolis force because it can be more important than the other inertial term, given by the fluid acceleration Du/Dt, in large-scale ocean flows u (Haller and Sapsis2008; Monroy et al.2017).

The effect of unresolved scales will also be estimated by keeping the original NEMO velocity field u but modifying the equation of motion, Eq. (1), by adding a stochastic noise term:

(B2) v = u + v s + W .

W(t)=(2Dhξx(t),2Dhξy(t),2Dvξz(t)), where ξ(t) is a vector Gaussian white noise process (independent for each particle) of zero mean and with correlations given by ξi(t1)ξj(t2)=δijδ(t1-t2),i,j=x,y,z. Thus, the horizontal and vertical intensities of this term are given by Dh by Dv, respectively.

Figure B1Solid line indicates the average horizontal distance dhI (a) and the average vertical distance dvI (b) of particles released from the same initial location but integrated with Eqs. (1) and (B1) in the NEMO velocity field. Shaded region indicates the range of the distances among the individual pairs of particles.


Figure B2Same as Fig. B1 for the comparison of Eqs. (1) and (B2).


The statistical properties are chosen to be similar to the ones expected for oceanic motions below the scales resolved by the numerical model (Monroy et al.2017; Kaandorp et al.2020). To do so, we use for Dh Okubo's empirical formulation (Okubo1971) that parameterizes the effective horizontal eddy-diffusion below a spatial scale as Dh()=2.055×10-41.55 m2/s. Taking for the horizontal resolution of our numerical model (1/12) we obtain Dh=7.25 m2/s. Since the Okubo formula is an empirical fit to surface motions, and effective horizontal diffusivity should be weaker below the thermocline, our results provide an upper bound for the error associated with unresolved scales of fluid motion. For the vertical diffusivity we take Dv=10-5 m2/s.

Table B1Average horizontal and vertical pairwise particle distances d and single-particle displacements r after an integration time of 10 d. See text. d/r is indicated by percentages in parentheses.

Download Print Version | Download XLSX

Figure B3Average horizontal distance dhI (a) and average vertical distance dvI (b) of particles released from the same initial location but integrated with Eqs. (1) and (B1) in the NEMO velocity field. The different lines are obtained from particles from different ranges of densities characterized by the indicated ranges in β.


Figure B4Same as Fig. B3 for the comparison of Eqs. (1) and (B2).


In order to compare the different equations of motion, we release a large number N=78 803 of particles over the whole Mediterranean at 1 m depth on 8 January 2000. We associate a β parameter with each particle by selecting it from a random uniform distribution in the range β[0.8,1), and once it is selected it remains fixed at all times for the corresponding particle. High values of β (close to 1) correspond to more buoyant plastic particles, whereas low values of β correspond to high settling velocities. The corresponding range of velocities is vs[1.776×10-3,0) m/s. We integrate the particle trajectories using Eq. (1) and also, from the same initial conditions, using the corrected dynamics in Eqs. (B1) or (B2), all with the same NEMO velocity field u. The horizontal and vertical distances between particles released from the same point but evolved with different equations are compared by calculating the following averages over particles:

(B3) d h I ( t ) = 1 N k = 1 N | r k 0 ( t ) - r k I ( t ) | , d v I ( t ) = 1 N k = 1 N | z k 0 ( t ) - z k I ( t ) | .

rk=(xk,yk) is the horizontal position of particle k, zk is its vertical position, and the superindices 0 and I indicate that the particle trajectory has been integrated by using the reference equation (Eq. 1) or the one containing inertial corrections (Eq. B1). The quantities dhW(t) and dvW(t), comparing Eq. (1) with the dynamics (B2) modeling small-scale flow effects, are defined analogously.

In Fig. B1, we display the average distances dhI(t) and dvI(t) as a function of time, characterizing the corrections by inertial terms to the simple dynamics of Eq. (1). Analogously, Fig. B2 displays the average distances dhW(t) and dvW(t) as a function of time, characterizing the estimated corrections arising from small scales unresolved by the NEMO velocity field. The effect induced by the inertial terms is very small and clearly negligible. The impact of W and thus of small unresolved scales is larger.

To evaluate the different effects more quantitatively, we summarize in Table B1, considering particles separately in different density ranges (given by the ranges in β and associated vs), the average horizontal and vertical pairwise particle distances, calculated after integrating the different dynamics for 10 d. To fully appreciate the importance of these numbers, in the two final columns we compute the horizontal and vertical average of the total distance r traveled by the particles (using Eq. 1) during the same interval of time. While the influence of the inertial terms is completely irrelevant both in a relative and an absolute sense for any realistic application, more care has to be taken with regards to the unresolved scales. Although the vertical error associated with the latter remains small, its relative importance hugely increases with decreasing settling velocity. (Indeed, it would tend to infinity for approaching neutral buoyancy.) At the same time, the relative horizontal error is the biggest for the fastest-sinking particles and is well above 10 % for them.

We can also observe the time evolution of the d entries of this table in Figs. B3 and B4.

Appendix C: The effect of bathymetry on vertical dispersion

We investigate here if the finite and spatially varying depth of the basin (the bathymetry) could affect the conclusions in Sect. 4.4 about vertical dispersion. The falling particles released from different locations reach the seafloor at different times, and then computing some statistics over these particles involves an increasingly narrow set of particles. Note that the bottom boundary (i.e., the seafloor) extends from the surface (close to the coast) to the deepest point of the basin, it is thus relevant at any time during the simulation. We intend to exclude three different effects arising from the continual removal of particles: (i) distortion of the shape of the distribution close to the boundary, (ii) poor quality of the statistics when many particles have already been lost, (iii) decrease in the geographical area sampled by the particles.

Figure C1Variance of depth reached by all sinking particles (in black) and an adaptively restricted subset of them (in yellow) as a function of time. See text for details. Straight lines represent power laws for reference, with exponents 1 and 2. The number N of particles considered for the computation of the variance is also shown.


We start with effect (i) by comparing, in Fig. C1, the variance presented in the main text (Fig. 7), computed over all sinking (but not yet sedimented) particles, to the variance computed over a restricted set of particles. This restricted set contains only those particles at the positions of which the bathymetry Z satisfies Z>zi-3σz, where the average zi and the standard deviation σz are taken with respect to the original (unrestricted) set of particles. This restriction is adaptive and ensures that the seafloor is sufficiently far for its effect on the particle distribution to be negligible at the positions of all particles kept for the computation. Figure C1 shows that the difference in the results between the full set and the restricted one is negligible for all three settling velocities considered, except perhaps for the drop at the very end of the time evolution, after the constant section, in Fig. C1a. This drop is, however, very short compared to the bulk of the sinking process and thus have minor importance. Furthermore, the distribution of particles so close to the bottom (see Fig. 9) should in any case be strongly influenced by resuspension and remixing by bottom currents (Kane et al.2020). For the rest of the time evolution, we can be confident that a possible distortion of the particle distribution induced by the boundaries has no impact on the variance curves. Such distortion may be present, but the number of particles close to the seafloor is very small at any given time instant (see the small difference between the numbers N of particles considered for the two kinds of computation in Fig. C1) so that they have a negligible contribution to the originally computed variance. This is a result of the relatively weak dispersion; artifacts might be found for broader distributions, possibly including a hypothetical continuation of the v6 simulation. The time evolution of N in Fig. C1 also assures us that poor-quality statistics, effect (ii), does not arise until the final drop discussed in the previous paragraph. Except for that very final section, the variance is estimated from a sufficiently large number of particles to keep the relative error of the sample variance (with respect to the population variance) very low. This is so because, under standard assumptions, the ratio between the sample and the population variance should be close to a chi-square random variable with N−1 degrees of freedom, divided by N−1 (Douillet2009).

The time evolution of N also suggests that effect (iii) is avoided as well: in Fig. C1, there is no sudden drop in the number of particles during the simulations that could result in the changes in the slope of the curves of variance versus time. To further support this conclusion, we compute the time evolution of the variance over the particles initialized in a subregion of the whole Mediterranean; see Fig. 8. In particular, we choose the box of longitudes [5,6.9], and latitudes [37.5,42], corresponding to the Sea of Sardinia, where the bathymetry is deep enough to prevent particles from reaching the seafloor (except for the very last few time steps of the v153 configuration), so that the horizontal area sampled by the particles approximately remains constant (remember that horizontal displacements are small compared to geographical features).

Figure C2Variance of depth reached by all sinking particles (in black) and those initialized in the Sea of Sardinia (in blue; longitudes in [5,6.9], and latitudes in [37.5,42]) as a function of time. See text for details. Straight lines represent power laws for reference, with exponents 1 and 2.


According to Fig. C2, the character of the dispersion in the Sea of Sardinia is nearly identical to that in the whole Mediterranean, except after the crossover from normal diffusive to ballistic dispersion in the v153 case (Fig. C2a). The smaller variance in the velocity should naturally lead to a later crossover to the long-term behavior in the Sea of Sardinia (see Sect. 4.4); it is nevertheless doubtful that the crossover should fall outside the simulation time. As the crossover is not observed, one might speculate that horizontal mixing might just get strong enough to suppress ballistic dispersion, similarly to the v68 case (see the discussion of Fig. 7), for which the time evolution of the variance is very similar in the Sea of Sardinia and the whole Mediterranean (Fig. C2b). Figure 8 suggests that the characteristic patch size is smaller in the western basin of the Mediterranean, including the Sea of Sardinia, than in the (considerably larger) eastern one, which makes inter-patch mixing easier. We conclude that the only substantial difference between the dispersion in the Sea of Sardinia and the whole Mediterranean may originate from the smaller extension and some special characteristics of the former, and the decrease in the area sampled in the whole-Mediterranean simulation presumably has no effect on the results.

Based on these analyses, we believe that the findings of Sect. 4.4 are unaffected by the boundary and thus have general relevance for oceanic dispersion.

Data availability

The velocity field from the NEMO simulation used in our study can be downloaded from (Yool et al.2013). Parcels is a set of Python classes developed under the TOPIOS project and is accessible at (van Sebille et al.2021). The scripts for running the particle transport simulations are available upon request from Rebeca de la Fuente,

Author contributions

All authors designed research. RdlF performed the simulations. RdlF and GD analyzed data. RdlF, GD, EHG and CL prepared the first draft, and all authors reviewed and edited the paper.

Competing interests

The authors declare that they have no conflict of interest.


We thank two anonymous referees for useful comments.

Financial support

This research has been supported by the H2020 European Research Council (TOPIOS grant no. 715386) and the Spanish State Research Agency through the María de Maeztu Program for Units of Excellence in R&D (grant no. MDM-2017-0711). Rebeca de la Fuente was also supported by MINECO, Spain, through the FPI program (fellowship no. BES-2016-078416). Gábor Drótos was also supported by CAIB and the European Social Fund through its postdoctoral program (grant no. PD/020/2018) and by NKFIH, Hungary (grant no. NKFI-124256). The publication fee was supported by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

Review statement

This paper was edited by Piers Chapman and reviewed by two anonymous referees.


Andrady, A. L.: Persistence of Plastic Litter in the Oceans, pp. 57–72, Springer International Publishing, Cham,, 2015. a, b, c, d

Andrady, A. L.: The plastic in microplastics: A review, Mar. Poll. Bull., 119, 12–22,, 2017. a, b

Arthur, C., Baker, J., and Bamford, H.: Proceedings of the international research workshop on the occurrence, effects and fate of microplastic marine debris, NOAA Technical Memorandum NOS-OR&R30, 2009. a, b

Bagaev, A., Mizyuk, A., Khatmullina, L., Isachenko, I., and Chubarenko, I.: Anthropogenic fibres in the Baltic Sea water column: Field data, laboratory and numerical testing of their motion, Sci. Total Environ., 599–600, 560–571,, 2017. a, b, c, d, e, f, g, h, i

Balachandar, S. and Eaton, J. K.: Turbulent dispersed multiphase flow, Annu. Rev. Fluid. Mech., 42, 111–133,, 2010. a

Balkovsky, E., Falkovich, G., and Fouxon, A.: Intermittent distribution of inertial particles in turbulent flows, Phys. Rev. Lett., 86, 2790,, 2001. a

Ballent, A., Pando, S., Purser, A., Juliano, M. F., and Thomsen, L.: Modelled transport of benthic marine microplastic pollution in the Nazaré Canyon, Biogeosciences, 10, 7957–7970,, 2013. a

Barrows, A., Cathey, S., and Petersen, C.: Marine environment microfiber contamination: Global patterns and the diversity of microparticle origins, Environ. Poll., 237, 275–284,, 2018. a

Bergmann, M., Wirzberger, V., Krumpen, T., Lorenz, C., Primpke, S., Tekman, M. B., and Gerdts, G.: High Quantities of Microplastic in Arctic Deep-Sea Sediments from the HAUSGARTEN Observatory, Environ. Sci. Technol., 51, 11000–11010,, 2017. a, b, c, d

Berloff, P. S. and McWilliams, J. C.: Material transport in oceanic gyres. Part II: Hierarchy of stochastic models, J. Phys. Oceanogr., 32, 797–830, 2002. a, b

Bond, T., Ferrandiz-Mas, V., Felipe-Sotelo, M., and van Sebille, E.: The occurrence and degradation of aquatic plastic litter based on polymer physicochemical properties: A review, Critical Reviews in Environmental Science and Technology, 48, 685–722,, 2018. a

Carr, S. A.: Sources and dispersive modes of micro-fibers in the environment, Int. Environ. Assess. Manag., 13, 466–469,, 2017. a, b

Cartwright, J. H. E., Feudel, U., Károlyi, G., De Moura, A., Piro, O., and Tél, T.: Dynamics of finite-size particles in chaotic fluid flows, in: Nonlinear dynamics and chaos: advances and perspectives, edited by: Thiel, M., Kurths, J., Romano, M., Károlyi, G., and Moura, A., pp. 51–87, Springer, Berlin, Heidelberg,, 2010. a

Choy, C. A., Robison, B. H., Gagne, T. O., Erwin, B., Firl, E., Halden, R. U., Hamilton, J. A., Katija, K., Lisin, S. E., Rolsky, C., and S. Van Houtan, K.: The vertical distribution and biological transport of marine microplastics across the epipelagic and mesopelagic water column, Sci. Rep., 9, 1–9, 2019. a, b

Chubarenko, I., Bagaev, A., Zobkov, M., and Esiukova, E.: On some physical and dynamical properties of microplastic particles in marine environment, Mar. Poll. Bull., 108, 105–112, 2016. a, b, c, d

Corcoran, P. L.: Benthic plastic debris in marine and fresh water environments, Environ. Sci. Processes Impacts, 17, 1363–1369,, 2015. a

Cózar, A., Echevarría, F., González-Gordillo, J. I., Irigoien, X., Úbeda, B., Hernández-León, S., Palma, Á. T., Navarro, S., García-de Lomas, J., Ruiz, A., Fernández-de Puelles, M. L., and Duarte, C. M.: Plastic debris in the open ocean, P. Natl. Acad. Sci., 111, 10239–10244,, 2014. a

Cózar, A., Sanz-Martín, M., Martí, E., González-Gordillo, J. I., Ubeda, B., Gálvez, J. Á., Irigoien, X., and Duarte, C. M.: Plastic accumulation in the Mediterranean Sea, PLoS One, 10, e0121762,, 2015. a, b

de Madron, X. D., Castaing, P., Nyffeler, F., and Courp, T.: Slope transport of suspended particulate matter on the Aquitanian margin of the Bay of Biscay, Deep Sea Res. Pt. II, 46, 2003–2027, 1999. a

Delandmeter, P. and van Sebille, E.: The Parcels v2.0 Lagrangian framework: new field interpolation schemes, Geosci. Model Dev., 12, 3571–3584,, 2019. a

Douillet, P. L.: Sampling distribution of the variance, in: Proceedings of the 2009 Winter Simulation Conference (WSC), pp. 403–414, 2009. a

Dris, R., Gasperi, J., Saad, M., Mirande, C., and Tassin, B.: Synthetic fibers in atmospheric fallout: A source of microplastics in the environment?, Mar. Poll. Bull., 104, 290–293,, 2016. a

Drótos, G., Monroy, P., Hernández-García, E., and López, C.: Inhomogeneities and caustics in the sedimentation of noninertial particles in incompressible flows, Chaos: An Interdisciplinary J. Nonlin. Sci., 29, 013115,, 2019. a, b

Eakins, B. and Sharman, G.: Volumes of the World's Oceans from ETOPO1, NOAA National Geophysical Data Center, Boulder, CO, available at: (last access: 2 March 2021), 2010. a, b

Egger, M., Sulu-Gambari, F., and Lebreton, L.: First evidence of plastic fallout from the North Pacific Garbage Patch, Sci. Rep., 10, 1–10,, 2020. a, b, c

Enders, K., Lenz, R., Stedmon, C. A., and Nielsen, T. G.: Abundance, size and polymer composition of marine microplastics 10 μm in the Atlantic Ocean and their modelled vertical distribution, Mar. Poll. Bull., 100, 70–81,, 2015. a, b, c, d

Erni-Cassola, G., Gibson, M. I., Thompson, R. C., and Christie-Oleza, J. A.: Lost, but Found with Nile Red: A Novel Method for Detecting and Quantifying Small Microplastics (1 mm to 20 μm) in Environmental Samples, Environ. Sci. Technol., 51, 13641–13648,, 2017. a, b

Erni-Cassola, G., Zadjelovic, V., Gibson, M. I., and Christie-Oleza, J. A.: Distribution of plastic polymer types in the marine environment; A meta-analysis, J. Hazard. Mater., 369, 691–698, 2019. a, b, c, d, e, f, g

Filella, M.: Questions of size and numbers in environmental research on microplastics: methodological and conceptual aspects, Environ. Chem., 12, 527–538,, 2015. a, b

Fischer, V., Elsner, N. O., Brenke, N., Schwabe, E., and Brandt, A.: Plastic pollution of the Kuril–Kamchatka Trench area (NW pacific), Deep Sea Res. Pt. II, 111, 399–405,, 2015. a, b, c, d, e

Gago, J., Carretero, O., Filgueiras, A., and Viñas, L.: Synthetic microfibers in the marine environment: A review on their occurrence in seawater and sediments, Mar. Poll. Bull., 127, 365–376,, 2018. a, b

Ganser, G. H.: A rational approach to drag prediction of spherical and nonspherical particles, Powder Technol., 77, 143–152,, 1993. a

GEBCO Compilation Group: General Bathymetric Chart of the Oceans (GEBCO) grid,, 2020. a

Haller, G. and Sapsis, T.: Where do inertial particles go in fluid flows?, Physica D, 237, 573–583, 2008. a, b

Hardy, J.: The sea surface microlayer: Biology, chemistry and anthropogenic enrichment, Prog. Oceanogr., 11, 307–328,, 1982. a

Heussner, S., de Madron, X. D., Calafat, A., Canals, M., Carbonne, J., Delsaut, N., and Saragoni, G.: Spatial and temporal variability of downward particle fluxes on a continental slope: Lessons from an 8-yr experiment in the Gulf of Lions (NW Mediterranean), Mar. Geol., 234, 63–92, 2006. a

Hidalgo-Ruz, V., Gutow, L., Thompson, R. C., and Thiel, M.: Microplastics in the marine environment: a review of the methods used for identification and quantification, Environ. Sci. Technol., 46, 3060–3075, 2012. a

Horton, A. A. and Dixon, S. J.: Microplastics: An introduction to environmental transport processes, Wiley Interdisciplinary Reviews: Water, 5, e1268,, 2018. a

Jambeck, J. R., Geyer, R., Wilcox, C., Siegler, T. R., Perryman, M., Andrady, A., Narayan, R., and Law, K. L.: Plastic waste inputs from land into the ocean, Science, 347, 768–771, 2015. a

Jiménez, J.: Ocean turbulence at millimeter scales, Sci. Mar., 61, 47–56, 1997. a, b, c, d, e

Kaandorp, M. L. A., Dijkstra, H. A., and van Sebille, E.: Closing the Mediterranean Marine Floating Plastic Mass Budget: Inverse Modeling of Sources and Sinks, Environ. Sci. Technol., 54, 11980–11989,, 2020. a, b, c, d, e

Kaiser, D., Kowalski, N., and Waniek, J. J.: Effects of biofouling on the sinking behavior of microplastics, Environ. Res. Lett., 12, 124003,, 2017. a

Kaiser, D., Estelmann, A., Kowalski, N., Glockzin, M., and Waniek, J. J.: Sinking velocity of sub-millimeter microplastic, Mar. Poll. Bull., 139, 214–220,, 2019. a

Kane, I. A. and Clare, M. A.: Dispersion, Accumulation, and the Ultimate Fate of Microplastics in Deep-Marine Environments: A Review and Future Directions, Front. Earth Sci., 7, 80,, 2019. a, b, c, d, e, f

Kane, I. A., Clare, M. A., Miramontes, E., Wogelius, R., Rothwell, J. J., Garreau, P., and Pohl, F.: Seafloor microplastic hotspots controlled by deep-sea circulation, Science, 6495, 1140–1145,, 2020. a, b, c, d

Kaneko, S., Nakamura, T., Inada, F., Kato, M., Ishihara, K., Nishihara, T., and Langthjem, M. A.: Chapter 8 – Vibrations in Fluid–Structure Interaction Systems, in: Flow-induced Vibrations (Second Edition), pp. 359–401, Academic Press, Oxford, UK, second edition edn.,, 2014. a

Kanhai, L. D. K., Gårdfeldt, K., Lyashevska, O., Hassellöv, M., Thompson, R. C., and O'Connor, I.: Microplastics in sub-surface waters of the Arctic Central Basin, Mar. Poll. Bull., 130, 8–18,, 2018. a, b

Khatmullina, L. and Isachenko, I.: Settling velocity of microplastic particles of regular shapes, Mar. Poll. Bull., 114, 871–880,, 2017. a

Kooi, M. and Koelmans, A. A.: Simplifying Microplastic via Continuous Probability Distributions for Size, Shape, and Density, Environ. Sci. Technol. Lett., 6, 551–557,, 2019. a

Kooi, M., van Nes, E. H., Scheffer, M., and Albert, A.: Ups and downs in the ocean: effects of biofouling on vertical transport of microplastics, Environ. Sci. Technol., 51, 7963–7971,, 2017. a, b, c, d

Kowalski, N., Reichardt, A. M., and Waniek, J. J.: Sinking rates of microplastics and potential implications of their alteration by physical, biological, and chemical factors, Mar. Poll. Bull., 109, 310–319,, 2016. a

Kukulka, T., Proskurowski, G., Morét-Ferguson, S., Meyer, D., and Law, K.: The effect of wind mixing on the vertical distribution of buoyant plastic debris, Geophys. Res. Lett., 39, L07601,, 2012. a

Kvale, K., Prowe, A., Chien, C.-T., Landolfi, A., and Oschlies, A.: The global biological microplastic particle sink, Sci. Rep., 10, 1–12,, 2020. a

Lebreton, L., Slat, B., Ferrari, F., Sainte-Rose, B., Aitken, J., Marthouse, R., Hajbane, S., Cunsolo, S., Schwarz, A., Levivier, A., Noble, K., Debeljak, P., Maral, H., Schoeneich-Argent, R., Brambini, R., and Reisser, J.: Evidence that the Great Pacific Garbage Patch is rapidly accumulating plastic, Sci. Rep., 8, 1–15,, 2018. a

Leith, D.: Drag on Nonspherical Objects, Aerosol Sci. Technol., 6, 153–161,, 1987. a

Lindeque, P. K., Cole, M., Coppock, R. L., Lewis, C. N., Miller, R. Z., Watts, A. J., Wilson-McNeal, A., Wright, S. L., and Galloway, T. S.: Are we underestimating microplastic abundance in the marine environment? A comparison of microplastic capture with nets of different mesh-size, Environ. Pollut., 265, 114721,, 2020. a, b

Liu, G., Bracco, A., and Passow, U.: The influence of mesoscale and submesoscale circulation on sinking particles in the northern Gulf of Mexico, Elementa, Science of the Anthropocene, 6, 36,, 2018. a

Liubartseva, S., Coppini, G., Lecci, R., and Clementi, E.: Tracking plastics in the Mediterranean: 2D Lagrangian model, Mar. Poll. Bull., 129, 151–162,, 2018. a, b

Long, M., Moriceau, B., Gallinari, M., Lambert, C., Huvet, A., Raffray, J., and Soudant, P.: Interactions between microplastics and phytoplankton aggregates: impact on their respective fates, Mar. Chem., 175, 39–46,, 2015. a

Macias, D., Cózar, A., Garcia-Gorriz, E., González-Fernández, D., and Stips, A.: Surface water circulation develops seasonally changing patterns of floating litter accumulation in the Mediterranean Sea. A modelling approach, Mar. Poll. Bull., 149, 110619,, 2019. a

Madec, G.: the NEMO team, Note du Pôle de modélisation, Institut Pierre-Simon Laplace (IPSL), France, NEMO ocean engine, 27, 2008. a

Madec, G. and Imbard, M.: A global ocean mesh to overcome the North Pole singularity, Clim. Dynam., 12, 381–388, 1996. a

Martin, J., Lusher, A., Thompson, R. C., and Morley, A.: The Deposition and Accumulation of Microplastics in Marine Sediments and Bottom Water from the Irish Continental Shelf, Sci. Rep., 7, 10772,, 2017. a, b

Maxey, M. R. and Riley, J. J.: Equation of motion for a small rigid sphere in a nonuniform flow, Phys. Fluids, 26, 883–889, 1983. a, b, c

Michaelides, E. E.: Hydrodynamic force and heat/mass transfer from particles, bubbles, and drops–the Freeman Scholar Lecture, J. Fluids Eng., 125, 209–238,, 2003. a

Michels, J., Stippkugel, A., Lenz, M., Wirtz, K., and Engel, A.: Rapid aggregation of biofilm-covered microplastics with marine biogenic particles, P. Roy. Soc. B, 285, 20181203,, 2018. a

Monroy, P., Hernández-García, E., Rossi, V., and López, C.: Modeling the dynamical sinking of biogenic particles in oceanic flow, Nonlin. Processes Geophys., 24, 293–305,, 2017. a, b, c, d, e, f, g, h, i, j, k, l

Monroy, P., Drótos, G., Hernández-García, E., and López, C.: Spatial Inhomogeneities in the Sedimentation of Biogenic Particles in Ocean Flows: Analysis in the Benguela Region, J. Geophys. Res., 124, 4744–4762, 2019. a, b

Mountford, A. and Morales Maqueda, M.: Eulerian Modeling of the Three-Dimensional Distribution of Seven Popular Microplastic Types in the Global Ocean, J. Geophys. Res.-Oceans, 124, 8558–8573, 2019. a, b, c, d

Neufeld, Z. and Hernández-García, E.: Chemical and biological processes in fluid flows: a dynamical systems approach, Imperial College Press, London, UK, 2009. a

Nooteboom, P. D., Delandmeter, P., van Sebille, E., Bijl, P. K., Dijkstra, H. A., and von der Heydt, A. S.: Resolution dependency of sinking Lagrangian particles in ocean general circulation models, PLOS ONE, 15, e0238650,, 2020. a, b, c

Okubo, A.: Oceanic diffusion diagrams, Deep Sea Res. Oceanogr. Abstracts, 18, 789–802,, 1971. a

Pawlowicz, R.: What every oceanographer needs to know about TEOS-10 (The TEOS-10 Primer), available at: (last access: 2 March 2021), 2010. a

Pedrotti, M. L., Petit, S., Elineau, A., Bruzaud, S., Crebassa, J.-C., Dumontet, B., Martí, E., Gorsky, G., and Cózar, A.: Changes in the floating plastic pollution of the Mediterranean Sea in relation to the distance to land, PloS one, 11, e0161581,, 2016. a

Peng, X., Chen, M., Chen, S., Dasgupta, S., Xu, H., Ta, K., Du, M., Li, J., Guo, Z., and Bai, S.: Microplastics contaminate the deepest part of the world's ocean, Geochem. Perspect. Lett., 9, 1–5,, 2018. a, b, c, d, e, f, g

Poulain, M., Mercier, M. J., Brach, L., Martignac, M., Routaboul, C., Perez, E., Desjean, M. C., and Ter Halle, A.: Small microplastics as a main contributor to plastic mass balance in the North Atlantic subtropical gyre, Environ. Sci. Technol., 53, 1157–1164, 2018. a

Reineccius, J., Appelt, J.-S., Hinrichs, T., Kaiser, D., Stern, J., Prien, R. D., and Waniek, J. J.: Abundance and characteristics of microfibers detected in sediment trap material from the deep subtropical North Atlantic Ocean, Sci. Total Environ., 738, 140354,, 2020. a

Reynolds, A. M.: On Lagrangian stochastic modelling of material transport in oceanic gyres, Physica D, 172, 124–138, 2002. a, b

van Sebille, E., Delandmeter, P., Lange, M., Kehl, C., Brouwer, R., Reijnders, D., Reint, Rath, W., pdnooteboom, Phillips, J. S., Simnator101, Gibson, A., Kronborg, J., Petersik, P., Thomas-95, Wichmann, D., pierrick-giffard, Busecke, J., Sterl, M., Edwards, R., SLYpma, Mol, W., Klink, D., hart-davis, Vandenbulcke, L., Treep, J., Walbridge, S., Novaes de Santana, A., Ham, D. A., and Vettoretti, G.: OceanParcels/parcels: Parcels v2.2.2: a Lagrangian Ocean Analysis tool for the petascale age (Version v2.2.2), Zenodo,, 2021. a

Seong, J. C., Mulcahy, K. A., and Usery, E. L.: The sinusoidal projection: A new importance in relation to global image data, The Professional Geographer, 54, 218–225, 2002. a

Siegel, D. and Deuser, W.: Trajectories of sinking particles in the Sargasso Sea: modeling of statistical funnels above deep-ocean sediment traps, Deep-Sea Res. I, 44, 1519–1541,, 1997. a, b

Song, Y. K., Hong, S. H., Jang, M., Kang, J.-H., Kwon, O. Y., Han, G. M., and Shim, W. J.: Large Accumulation of Micro-sized Synthetic Polymer Particles in the Sea Surface Microlayer, Environ. Sci. Technol., 48, 9014–9021,, 2014. a, b, c, d, e, f

Song, Y. K., Hong, S. H., Jang, M., Han, G. M., Jung, S. W., and Shim, W. J.: Combined Effects of UV Exposure Duration and Mechanical Abrasion on Microplastic Fragmentation by Polymer Type, Environ. Sci. Technol., 51, 4368–4376,, 2017. a

Soto-Navarro, J., Jordá, G., Deudero, S., Alomar, C., Amores, A., and Compa, M.: 3D hotspots of marine litter in the Mediterranean: A modeling study, Mar. Poll. Bull., 155, 111159,, 2020. a, b, c

Sozza, A., Drotos, G., Hernández-García, E., and López, C.: Accumulated densities of sedimenting particles in turbulent flows, Phys. Fluids, 32, 075104,, 2020. a

Stabholz, M., Durrieu de Madron, X., Canals, M., Khripounoff, A., Taupier-Letage, I., Testor, P., Heussner, S., Kerhervé, P., Delsaut, N., Houpert, L., Lastras, G., and Dennielou, B.: Impact of open-ocean convection on particle fluxes and sediment dynamics in the deep margin of the Gulf of Lions, Biogeosciences, 10, 1097–1116,, 2013. a

Stolle, C., Nagel, K., Labrenz, M., and Jürgens, K.: Succession of the sea-surface microlayer in the coastal Baltic Sea under natural and experimentally induced low-wind conditions, Biogeosciences, 7, 2975–2988,, 2010. a

Suaria, G., Avio, C. G., Mineo, A., Lattin, G. L., Magaldi, M. G., Belmonte, G., Moore, C. J., Regoli, F., and Aliani, S.: The Mediterranean Plastic Soup: synthetic polymers in Mediterranean surface waters, Sci. Rep., 6, 37551,, 2016. a, b, c, d

Thompson, R. C., Olsen, Y., Mitchell, R. P., Davis, A., Rowland, S. J., John, A. W., McGonigle, D., and Russell, A. E.: Lost at sea: where is all the plastic?, Science (Washington), 304, 838,, 2004. a

Turner, A. and Holmes, L.: Occurrence, distribution and characteristics of beached plastic production pellets on the island of Malta (central Mediterranean), Mar. Poll. Bull., 62, 377–381, 2011.  a

van Cauwenberghe, L., Vanreusel, A., Mees, J., and Janssen, C. R.: Microplastic pollution in deep-sea sediments, Environ. Pollut., 182, 495–499,, 2013. a

van Sebille, E., Wilcox, C., Lebreton, L., Maximenko, N., Hardesty, B. D., Van Franeker, J. A., Eriksen, M., Siegel, D., Galgani, F., and Law, K. L.: A global inventory of small floating plastic debris, Environ. Res. Lett., 10, 124006,, 2015. a

van Sebille, E., Aliani, S., Law, K. L., Maximenko, N., Alsina, J. M., Bagaev, A., Bergmann, M., Chapron, B., Chubarenko, I., Cózar, A., Delandmeter P., Egger, M., Fox-Kemper, B., Garaba, S. P., Goddijn-Murphy, L., Hardesty, B. D., Hoffman, M. J., Isobe, A., Jongedijk, C. E., Kaandorp, M. L. A., Khatmullina, L., Koelmans, A. A., Kukulka, T., Laufkötter, C., Lebreton, L., Lobelle, D., Maes, C., Martinez-Vicente, V., Morales Maqueda, M. A., Poulain-Zarcos, M., Rodríguez, E., Ryan, P. G., Shanks, A. L., Shim, W. J., Suaria, G. Thiel, M., van den Bremer, T. S., and Wichmann, D.: The physical oceanography of the transport of floating marine debris, Environ. Res. Lett., 15, 023003,, 2020. a

Vidal, A. S., Pasqual, C., Kerhervé, P., Heussner, S., Calafat, A., Palanques, A., Durrieu de Madron, X., Canals, M., and Puigc, P.: Across margin export of organic matter by cascading events traced by stable isotopes, northwestern Mediterranean Sea, Limnol. Oceanogr., 54, 1488–1500, 2009. a

Voth, G. A. and Soldati, A.: Anisotropic Particles in Turbulence, Annu. Rev. Fluid Mech., 49, 249–276,, 2017. a, b

Wichmann, D., Delandmeter, P., and van Sebille, E.: Influence of near-surface currents on the global dispersal of marine microplastic, J. Geophys. Res.-Oceans, 124, 6086–6096, 2019. a

Woodall, L. C., Sanchez-Vidal, A., Canals, M., Paterson, G. L., Coppock, R., Sleight, V., Calafat, A., Rogers, A. D., Narayanaswamy, B. E., and Thompson, R. C.: The deep sea is a major sink for microplastic debris, Royal Soc. Open Sci., 1, 140317,, 2014. a

Wright, S., Ulke, J., Font, A., Chan, K., and Kelly, F.: Atmospheric microplastic deposition in an urban environment and an evaluation of transport, Environ. Int., 136, 105411,, 2020. a

Ye, S. and Andrady, A. L.: Fouling of floating plastic debris under Biscayne Bay exposure conditions, Mar. Poll. Bull., 22, 608–613,, 1991. a, b

Yool, A., Popova, E. E., and Anderson, T. R.: MEDUSA-2.0: an intermediate complexity biogeochemical model of the marine carbon cycle for climate change and ocean acidification studies, Geosci. Model Dev., 6, 1767–1811,, 2013 (data available at:, last access: 15 December 2019). a

Short summary
Plastic pollution is a major environmental issue affecting the oceans. The number of floating and sedimented pieces has been quantified by several studies. But their abundance in the water column remains mostly unknown. To fill this gap we model the dynamics of a particular type of particle, rigid microplastics sinking rapidly in open sea in the Mediterranean. We find they represent a small but appreciable fraction of the total sea plastic and discuss characteristics of their sinking motion.