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

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 Holmes, 2011), biofouling (Ye and Andrady, 1991; 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 Dixon, 2018).

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; Andrady, 2017; Erni-Cassola et al., 2019; Kane and Clare, 2019). 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 Andrady, 1991; 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 Maqueda, 2019; 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 Maqueda, 2019). 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; Carr, 2017; Gago et al., 2018; Wright et al., 2020), it is not obvious to explain their abundance on abyssal oceanic plains (Kane and Clare, 2019). 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 (Carr, 2017; 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 Isachenko, 2017; 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 Riley, 1983) 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 Clare, 2019; 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 (Andrady, 2015; Corcoran, 2015), 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 Clare, 2019; 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 Clare, 2019) 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; Andrady, 2015; Filella, 2015; 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 (Andrady, 2015). 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 (Hardy, 1982; 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; Andrady, 2017; Erni-Cassola et al., 2019) in Fig. 1: PP between 850–920 kg/m³, PE of 890–980 kg/m³, PS of 1040 kg/m³ (excluding its foamed version), PEST in the range 1100–1400 kg/m³, PA within 1120–1150 kg/m³ and acrylic of 1180 kg/m³. There is also some less abundant plastic like polytetrafluoroethylene (PTFE) which has higher densities, in the range 2100–2300 kg/m³.

Thus, the full range of microplastic particle densities in the ocean, denoted here as ρ_p, is 850–2300 kg/m³, and most of them have densities within the interval 850–1400 kg/m³. This has to be compared with the seawater density, which close to the surface has a conventional mean value of ρ_f=1025 km/m³ (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/m${}^{\mathrm{3}}<{\mathit{\rho }}_{\mathrm{f}}<\mathrm{1400}$ kg/m³.

Figure 1 Polymer 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; Andrady, 2015; 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 (Filella, 2015).

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 Deuser, 1997; 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 v_s as

$$\begin{array}{}\text{(1)}& \mathit{v}=\mathit{u}+{\mathit{v}}_{\mathrm{s}},\end{array}$$

with

$$\begin{array}{}\text{(2)}& {\mathit{v}}_{\mathrm{s}}=(\mathrm{1}-\mathit{\beta })\mathit{g}{\mathit{\tau }}_{\mathrm{p}},\mathit{\beta }={\displaystyle \frac{\mathrm{3}{\mathit{\rho }}_{\mathrm{f}}}{\mathrm{2}{\mathit{\rho }}_{\mathrm{p}}+{\mathit{\rho }}_{\mathrm{f}}}}\text{ and }{\mathit{\tau }}_{\mathrm{p}}={\displaystyle \frac{a^{\mathrm{2}}}{\mathrm{3}\mathit{\beta }\mathit{\nu }}}.\end{array}$$

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 $\mathrm{1}/\mathrm{12}$^∘ and 75 s levels in the vertical with data updated every 5 d (Madec, 2008; Madec and Imbard, 1996). Salinity and temperature are also extracted from that model. The Parcels Lagrangian framework (Delandmeter and van Sebille, 2019) 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 (Michaelides, 2003; Balkovsky et al., 2001; Cartwright et al., 2010; Haller and Sapsis, 2008), we start by choosing the simplified standard form of the more fundamental Maxey–Riley–Gatignol (MRG) equation (Maxey and Riley, 1983) 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

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}}=\mathit{\beta }{\displaystyle \frac{\mathrm{D}\mathit{u}}{\mathrm{D}t}}+{\displaystyle \frac{\mathit{u}-\mathit{v}+{\mathit{v}}_{\mathrm{s}}}{{\mathit{\tau }}_{\mathrm{p}}}}.\end{array}$$

Beyond sphericity, two conditions are needed for the validity of Eq. (3) (Monroy et al., 2017; Maxey and Riley, 1983): (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 $<\mathit{\eta }<\mathrm{2}$ mm for wind-driven turbulence in the upper ocean (Jiménez, 1997); (b) the particle Reynolds number $R{e}_{\mathrm{p}}=\frac{a|\mathit{v}-\mathit{u}|}{\mathit{\nu }}\approx \frac{a{v}_{\mathrm{s}}}{\mathit{\nu }}$ should satisfy Re_p≪1. Note that this last condition imposes restrictions on the values of the particles' density and size, partially via the settling velocity $v_{\mathrm{s}}=\left|{\mathit{v}}_{\mathrm{s}}\right|$. For the most abundant sinking microplastics, i.e., with densities ρ_p=1025–1400 kg/m³, we now determine the range of validity of Eq. (3) assuming $\mathit{\nu }=\mathrm{1.15}\times {\mathrm{10}}^{-\mathrm{6}}$ m²/s and ρ_f=1025 kg/m³ 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 v_s, 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énez, 1997), and Re_p=1 with a black line, which bound the area of validity (shaded in the plot). We also indicate v_s as a function of a for β=0.8 with the blue curve, corresponding to the upper bound to v_s 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 v_s<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énez, 1997), 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énez, 1997), 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 Soldati, 2017). 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 2 Settling velocities and particle sizes for which Eq. (3) holds. Kolmogorov scale is represented by the red line and Re_p=1 with a black line, which bound the area of validity. Blue curve corresponds to $v_{\mathrm{s}}=v_{\mathrm{s}}(\mathit{\beta },a)$ for β=0.8, the upper bound to v_s 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énez, 1997), so that the Stokes number $St={\mathit{\tau }}_{\mathrm{p}}/{\mathit{\tau }}_{\mathit{\eta }}$, which measures the importance of particle inertia in a turbulent flow, is very small (of the order of 10⁻³–10⁻²). 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 Eaton, 2010; Monroy et al., 2017; Drótos et al., 2019)

$$\begin{array}{}\text{(4)}& \mathit{v}\approx \mathit{u}+{\mathit{v}}_{\mathrm{s}}+{\mathit{\tau }}_{\mathrm{p}}(\mathit{\beta }-\mathrm{1}){\displaystyle \frac{\mathrm{D}\mathit{u}}{\mathrm{D}t}}.\end{array}$$

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 (Pawlowicz, 2010). 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/m³. This implies that for a nominal water density of ρ_f=1025 kg/m³ the value of β would be β=0.99, giving a sinking velocity v_s=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/m³. 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/m³, so that a slightly slower sinking is expected to occur.

Figure 3 Normalized 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/m³, and fluid density is obtained from the TEOS-10 equations. The vertical line indicates a conventional seawater density of 1025 kg/m³.

Figure 4 The distance, as a function of time, between trajectories obtained with constant nominal fluid density of ρ_f=1025 kg/m³ 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/m³, 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 $|{\mathbf{\text{x}}}^{\left(\mathrm{0}\right)}-{\mathbf{\text{x}}}^{\left(\mathrm{1}\right)}|$ and vertical $|z^{\left(\mathrm{0}\right)}-z^{\left(\mathrm{1}\right)}|$ distances and its average among particles. Trajectories x⁽⁰⁾(t) and x⁽¹⁾(t) are obtained with constant nominal fluid density (1025 kg/m³) 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/m³. 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/m³.

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/m³, 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/m³) 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 v_s.

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 $\times \mathrm{0.37}\times \mathrm{0.06}=\mathrm{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 Sharman, 2010; GEBCO Compilation Group, 2020) and estimate a residence time τ as the time of sinking to that mean depth. The residence time depends on the sinking velocity, $\mathit{\tau }=h/v_{\mathrm{s}}$, and thus on the physical properties of the microplastic particles. Assuming a seawater density ρ_f=1025 kg/m³, 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 $\stackrel{\mathrm{‾}}{\mathit{\tau }}\approx \mathrm{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\stackrel{\mathrm{‾}}{\mathit{\tau }}$: 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×10⁶ km³ (Eakins and Sharman, 2010) the estimated density would be ${\mathit{\rho }}_{\mathrm{V}}\approx \mathrm{7.7}\times {\mathrm{10}}^{-\mathrm{11}}$ kg/m³ (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 $\mathit{\rho }\left(z\right)/A\left(z\right)$, 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 v_s 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 5 The 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 $\mathit{\rho }\left(z\right)/A\left(z\right)$, 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 v_s=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 Deuser, 1997). 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 $\stackrel{\mathrm{̃}}{z}=\frac{z-t{v}_{\mathrm{s}}}{{\mathit{\sigma }}_z}$ where ${\mathit{\sigma }}_z^{\mathrm{2}}\equiv \langle (z_i-\langle z_i\rangle {)}^{\mathrm{2}}\rangle$ 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 6 The probability density function, estimated from a histogram of bin size 0.1, of all particles released in the Mediterranean in the rescaled variable $\stackrel{\mathrm{̃}}{z}=\frac{z-t{v}_{\mathrm{s}}}{{\mathit{\sigma }}_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ía, 2009), we now analyze this aspect in detail by considering how the variance of the vertical particle distribution, ${\mathit{\sigma }}_z^{\mathrm{2}}\left(t\right)$, 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 ${\mathit{\sigma }}_z^{\mathrm{2}}\sim t^{\mathit{\nu }}$ in different time intervals.

Figure 7 Variance 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 McWilliams, 2002; Reynolds, 2002). 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

$$\begin{array}{}\text{(5)}& z_i=\langle z_i\rangle +{\overline{\mathit{\omega }}}_{i}t+W_i,\end{array}$$

where 〈…〉 denotes, as before, an averaging over different particles. Here we are assuming that z_i−〈z_i〉 evolves according to the sum of a constant average velocity contribution ${\overline{\mathit{\omega }}}_i$ for sufficiently long times (a characteristic of the flow region traversed by particle i), and of W_i, a Wiener process representing fluctuations with zero mean and defining a diffusion coefficient D_i for each trajectory by $\stackrel{\mathrm{‾}}{W_i^{\mathrm{2}}}=D_{i}t$. 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 ${\overline{\mathit{\omega }}}_i\ne \mathrm{0}$), and D_i characterizes the strength of the fluctuations. Assuming $\langle {\overline{\mathit{\omega }}}_i{W}_i\rangle =\mathrm{0}$,

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_z^{\mathrm{2}}\equiv \langle (z_i-\langle z_i\rangle {)}^{\mathrm{2}}\rangle =\langle {\overline{\mathit{\omega }}}_i^{\mathrm{2}}\rangle t^{\mathrm{2}}+\langle D_i\rangle t,\end{array}$$

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=〈D_i〉, the crossover between the two regimes is obtained by equating the two terms as

$$\begin{array}{}\text{(7)}& t^{*}={\displaystyle \frac{D}{\langle {\overline{\mathit{\omega }}}_i^{\mathrm{2}}\rangle }}.\end{array}$$

To evaluate Eq. (7), we first estimate ${\overline{\mathit{\omega }}}_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 ${\overline{\mathit{\omega }}}_i$ is presented in Fig. 8, which confirms its patchiness throughout the Mediterranean, associated with mesoscale features. Computing $\langle {\overline{\mathit{\omega }}}_i^{\mathrm{2}}\rangle$ and fitting a line to ${\mathit{\sigma }}_z^{\mathrm{2}}\left(t\right)$ between t=1.4 and 4 d to estimate D, we obtain $t^{*}\approx \mathrm{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 ${\mathit{\sigma }}_z^{\mathrm{2}}$ 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 ${\overline{\mathit{\omega }}}_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^{*}\approx \mathrm{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 ${\overline{\mathit{\omega }}}_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 9 Variance 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/m³, 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 Maqueda, 2019; 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 Clare, 2019). 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 v_s. 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 v_s, 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 v_s can be obtained by balancing the drag force F_drag(v−u) (a function of the difference of the particle and the fluid velocities, v and u, respectively) with the resultant of gravitational and buoyancy forces:

$$\begin{array}{}\text{(A1)}& \mathrm{0}={\mathit{F}}_{\text{drag}}(\mathit{v}-\mathit{u})+V\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{f}}\right)\mathit{g},\end{array}$$

with $\mathit{v}-\mathit{u}={\mathit{v}}_{\mathrm{s}}$, 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

$$\begin{array}{}\text{(A2)}& {\mathit{F}}_{\text{drag}}^{\text{(sph)}}(\mathit{v}-\mathit{u})=-\mathrm{6}\mathit{\pi }\mathit{\mu }(\mathit{v}-\mathit{u})a,\end{array}$$

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

$$\begin{array}{}\text{(A3)}& {\mathit{F}}_{\text{drag}}^{\text{(non)}}(\mathit{v}-\mathit{u})=-\mathrm{6}\mathit{\pi }\mathit{\mu }(\mathit{v}-\mathit{u})\left({\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{a}_{\mathrm{n}}+{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{a}_{\mathrm{s}}\right),\end{array}$$

where a_n is the radius of the sphere with equivalent area projected on the plane perpendicular to the relative velocity v−u, and a_s is the radius of the sphere with equivalent total surface. From either of the last two equations, the settling velocity is obtained by substituting $\mathit{v}-\mathit{u}={\mathit{v}}_{\mathrm{s}}$ and solving Eq. (A1) for v_s. We denote the magnitudes of the settling velocities obtained from Eq. (A2) and Eq. (A3) by $v_{\mathrm{s}}^{\text{(sph)}}$ and $v_{\mathrm{s}}^{\text{(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

$$\begin{array}{}\text{(A4)}& q\equiv {\displaystyle \frac{v_{\mathrm{s}}^{\text{(non)}}}{v_{\mathrm{s}}^{\text{(sph)}}}}={\displaystyle \frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{V^{\text{(non)}}}{a^{\mathrm{2}}\left(\frac{\mathrm{1}}{\mathrm{3}}{a}_{\mathrm{n}}+\frac{\mathrm{2}}{\mathrm{3}}{a}_{\mathrm{s}}\right)}},\end{array}$$

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=\widehat{B}A<A$ and $C=\widehat{C}A<B<A$, where one or both of $\widehat{B}$ and $\widehat{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/\mathrm{2}$ and $a=B/\mathrm{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 a_s,

$$\begin{array}{}\text{(A5)}& a_{\mathrm{s}}={\left({\displaystyle \frac{AB+AC+BC}{\mathrm{2}\mathit{\pi }}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

However, a_n 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 a_n read as

$$\begin{array}{}\text{(A6)}& a_{\mathrm{n},X}={\left({\displaystyle \frac{ABC}{\mathit{\pi }X}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\end{array}$$

for X=A, B and C.

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

$$\begin{array}{}\text{(A7)}& {\displaystyle }{q}_X^{(A/\mathrm{2})}=\mathrm{9}{\mathit{\pi }}^{-\frac{\mathrm{1}}{\mathrm{2}}}\widehat{B}\widehat{C}{\left[{\left({\displaystyle \frac{\widehat{B}\widehat{C}}{\widehat{X}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}+{\mathrm{2}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\left(\widehat{B}+\widehat{C}+\widehat{B}\widehat{C}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\right]}^{-\mathrm{1}},\text{(A8)}& {\displaystyle }{q}_X^{(B/\mathrm{2})}={\widehat{B}}^{-\mathrm{2}}{q}_X^{(A/\mathrm{2})},\end{array}$$

where $X=\widehat{X}A$ has been introduced.

Figure A1 $q_A^{(A/\mathrm{2})}$ as a function of $\widehat{B}$ and $\widehat{C}$ for $\widehat{B},\widehat{C}\in [\mathrm{0.1},\mathrm{1}]$ with $\widehat{B}>\widehat{C}$. Dotted lines represent $q_A^{(A/\mathrm{2})}=\mathrm{0.1}$.

We plot $q_A^{(A/\mathrm{2})}$ in Fig. A1 as a function of $\widehat{B}$ and $\widehat{C}$ for $\widehat{B},\widehat{C}\in [\mathrm{0.1},\mathrm{1}]$ with $\widehat{B}>\widehat{C}$. Its range extends from 0.07 to 1.5 on this domain, but it drops below 0.1 only for $\widehat{C}$ very close to 0.1 and $\widehat{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 $q_A^{(B/\mathrm{2})}$ (not shown) on the same domain is between 0.2 and 7, and values above 4 are again restricted to very small $\widehat{C}$ and to $\widehat{B}<\mathrm{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 $\widehat{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 Soldati, 2017). 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 $\widehat{C}$ with $\widehat{B}\approx \mathrm{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⁻³–10⁻² 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,

$$\begin{array}{}\text{(B1)}& \mathit{v}=\mathit{u}+{\mathit{v}}_{\mathrm{s}}+{\mathit{\tau }}_{\mathrm{p}}(\mathit{\beta }-\mathrm{1})\left({\displaystyle \frac{\mathrm{D}\mathit{u}}{\mathrm{D}t}}+\mathrm{2}\mathbf{\Omega }\times \mathit{u}\right),\end{array}$$

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 $\mathrm{D}\mathit{u}/\mathrm{D}t$, in large-scale ocean flows u (Haller and Sapsis, 2008; 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:

$$\begin{array}{}\text{(B2)}& \mathit{v}=\mathit{u}+{\mathit{v}}_{\mathrm{s}}+\mathit{W}.\end{array}$$

$\mathit{W}\left(t\right)=\left(\sqrt{\mathrm{2}{D}_{\mathrm{h}}}{\mathit{\xi }}_x\right(t),\sqrt{\mathrm{2}{D}_{\mathrm{h}}}{\mathit{\xi }}_y(t),\sqrt{\mathrm{2}{D}_{\mathrm{v}}}{\mathit{\xi }}_z(t\left)\right)$, where ξ(t) is a vector Gaussian white noise process (independent for each particle) of zero mean and with correlations given by $\langle {\mathit{\xi }}_i\left(t_{\mathrm{1}}\right){\mathit{\xi }}_j\left(t_{\mathrm{2}}\right)\rangle ={\mathit{\delta }}_{ij}\mathit{\delta }(t_{\mathrm{1}}-t_{\mathrm{2}}),i,j=x,y,z$. Thus, the horizontal and vertical intensities of this term are given by D_h by D_v, respectively.

Figure B1 Solid line indicates the average horizontal distance $d_{\mathrm{h}}^I$ (a) and the average vertical distance $d_{\mathrm{v}}^I$ (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 B2 Same 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 D_h Okubo's empirical formulation (Okubo, 1971) that parameterizes the effective horizontal eddy-diffusion below a spatial scale ℓ as $D_{\mathrm{h}}\left(\mathrm{\ell }\right)=\mathrm{2.055}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{\ell }}^{\mathrm{1.55}}$ m²/s. Taking for ℓ the horizontal resolution of our numerical model ($\mathrm{1}/\mathrm{12}$^∘) we obtain D_h=7.25 m²/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 $D_{\mathrm{v}}={\mathrm{10}}^{-\mathrm{5}}$ m²/s.

Table B1 Average 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.

Figure B3 Average horizontal distance $d_{\mathrm{h}}^I$ (a) and average vertical distance $d_{\mathrm{v}}^I$ (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 B4 Same 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 $\mathit{\beta }\in [\mathrm{0.8},\mathrm{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 $v_{\mathrm{s}}\in [\mathrm{1.776}\times {\mathrm{10}}^{-\mathrm{3}},\mathrm{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:

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}& d_{\mathrm{h}}^I\left(t\right)={\displaystyle \frac{\mathrm{1}}{N}}\sum _{k=\mathrm{1}}^N\left|{\mathit{r}}_k^{\mathrm{0}}\right(t)-{\mathit{r}}_k^I(t\left)\right|,\\ & d_{\mathrm{v}}^I\left(t\right)={\displaystyle \frac{\mathrm{1}}{N}}\sum _{k=\mathrm{1}}^N\left|z_k^{\mathrm{0}}\right(t)-z_k^I(t\left)\right|.\end{array}\end{array}$$

${\mathit{r}}_k=(x_k,y_k)$ is the horizontal position of particle k, z_k 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 $d_{\mathrm{h}}^W\left(t\right)$ and $d_{\mathrm{v}}^W\left(t\right)$, comparing Eq. (1) with the dynamics (B2) modeling small-scale flow effects, are defined analogously.

In Fig. B1, we display the average distances $d_{\mathrm{h}}^I\left(t\right)$ and $d_{\mathrm{v}}^I\left(t\right)$ 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 $d_{\mathrm{h}}^W\left(t\right)$ and $d_{\mathrm{v}}^W\left(t\right)$ 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 v_s), 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 C1 Variance 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>\langle z_i\rangle -\mathrm{3}{\mathit{\sigma }}_z$, where the average 〈z_i〉 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 (Douillet, 2009).

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 $[\mathrm{5}{}^{\circ },\mathrm{6.9}{}^{\circ }]$, and latitudes $[\mathrm{37.5}{}^{\circ },\mathrm{42}{}^{\circ }]$, 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 C2 Variance of depth reached by all sinking particles (in black) and those initialized in the Sea of Sardinia (in blue; longitudes in $[\mathrm{5}{}^{\circ },\mathrm{6.9}{}^{\circ }]$, and latitudes in $[\mathrm{37.5}{}^{\circ },\mathrm{42}{}^{\circ }]$) 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.