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, Erik van Sebille
aa r X i v : . [ phy s i c s . f l u - dyn ] F e b Sinking microplastics in the water column: simulations in theMediterranean Sea
Rebeca de la Fuente , Gábor Drótos , Emilio Hernández-García , Cristóbal López , and Erik vanSebille IFISC (CSIC-UIB), Campus Universitat de les Illes Balears, Palma de Mallorca, Spain MTA–ELTE Theoretical Physics Research Group, Budapest, Hungary Institute for Marine and Atmospheric Research, Utrecht University, Utrecht, the Netherlands Centre for Complex Systems Studies, Utrecht University, Utrecht, the Netherlands
Correspondence:
Emilio Hernández-García (emilio@ifisc.uib-csic.es)
Abstract.
We study the vertical dispersion and distribution of negatively buoyant rigid microplastics within a realistic circula-tion model of the Mediterranean sea. We first propose an equation describing their idealized dynamics. In that framework, weevaluate the importance of some relevant physical effects: inertia, Coriolis force, small-scale turbulence and variable seawaterdensity, and bound the relative error of simplifying the dynamics to a constant sinking velocity added to a large-scale velocityfield. We then calculate the amount and vertical distribution of microplastic particles on the water column of the open oceanif their release from the sea surface is continuous at rates compatible with observations in the Mediterranean. The verticaldistribution is found to be almost uniform with depth for the majority of our parameter range. Transient distributions from flashreleases reveal a non-Gaussian character of the dispersion and various diffusion laws, both normal and anomalous. The originof these behaviors is explored in terms of horizontal and vertical flow organization.
Copyright statement.
To be inserted by journal staff
Approximately 8 million tonnes of plastics end up in the oceans every year (Jambeck et al., 2015). Nevertheless, only avery small percentage, around , remains on the surface (van Sebille et al., 2015; Choy et al., 2019). The rest leaves thesurface 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-drivenmixing presumably leads to an underestimation for the amount of particles remaining close to 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 understoodat present but would be crucial to properly evaluate the exposure of marine biota to this material, and formulate strategies forcleaning the oceans (Horton and Dixon, 2018). loating plastics and those that have beached or sedimented on the seafloor are relatively well studied through field cam-paigns (although 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 arerestricted to so-called underway samples, a few meters below the surface (e.g., Enders et al., 2015). However, e.g., Choy et al.(2019) reported that below the mixed layer and down to 1000 m depth in Monterey Bay, concentrations of plastics are largerthan at the surface (Thompson et al., 2004; Hidalgo-Ruz et al., 2012). Egger et al. (2020) found more plastic between 5 m and2000 m below the North Pacific Garbage Patch than at the surface. These findings turn out to mostly concern plastic piecesthat, 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 size,and we estimate their vertical distribution through the water column and their amount in the Mediterranean Sea. Microplasticparticles are among the most important contributors to marine plastic pollution (Arthur et al., 2009). Closely following thework of Monroy et al. (2017) for sinking biogenic particles but choosing particle properties to correspond to those of negativelybuoyant microplastics, we first justify the use of a simplified equation of motion, in which the plastic particle velocity is thesum of the ambient flow velocity and a sinking velocity depending on particle and water characteristics. In particular, weestimate the impact of some corrections to this simple dynamics and evaluate in detail the influence of the spatial variation ofthe seawater density on the plastic dispersion and sinking characteristics. For our Mediterranean case study, the impact of thevarying seawater density on particle trajectories can be comparable to the estimated effect of the neglected small scales belowthe hydrodynamical model’s resolution.We then estimate the amount of microplastic particles in the water column of the open Mediterranean. Our estimates relyon 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. Detailedconsideration of the transient dynamics identifies small non-Gaussian vertical dispersion around this simple sinking behavior,with transitions between anomalous and normal effective diffusion. 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 columnbeing 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 contributeto the microplastic content in the water column, the topic we are interested in this paper. However, it has been documentedexperimentally 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 aregenerated by biofouling (Ye and Andrady, 1991; Chubarenko et al., 2016) from positively buoyant plastic types or by othermechanisms of aggregation with organic matter, especially for small particle sizes (Kooi et al., 2017). n fact, the fallout from the North Pacific Garbage Patch almost entirely consists of plastic types nominally less densethan water (Egger et al., 2020). Although some of these immersed particles finally reach the sea bottom, their proportion insedimented plastic is minor except for the immediate vicinity of coasts where water is shallow. Most of these particles remainin the photic zone (Mountford and Morales Maqueda, 2019; Wichmann et al., 2019; Soto-Navarro et al., 2020). This suggeststhat 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 mech-anism by which microplastics denser than water can also be present within the water column is the finite time taken by themto reach the bottom. Under continuous release at the surface and sedimentation at the bottom, the transient falling would leadto a steady distribution for the amount of plastic in the water column at any given time, and this distribution has never beenestimated. Note that the Eulerian methodology of Mountford and Morales Maqueda (2019), treating sedimentation (i.e., depo-sition on the seafloor) by parametrization and thus leaving particles in the water column indefinitely long, is not suitable forthis 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 beenfound 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 sedi-ment traps in the open sea as well (Bagaev et al., 2017; Kanhai et al., 2018; Peng et al., 2018; Reineccius et al., 2020). Mostlyoriginating from land-based sources (Dris et al., 2016; Carr, 2017; Gago et al., 2018; Wright et al., 2020), it is not obviousto 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 settlingvelocity (Bagaev et al., 2017) that allows long distance horizontal transport (Nooteboom et al., 2020). In any case, it is difficultto estimate the amount of microfibers in the oceans due to sampling issues and to their absence from statistics of mismanagedplastic 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 thethree dimensions. Experimental estimates for the settling velocities of irregular fragments or other nonspherical particles havesuggested 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 qualitativeexploration of particle transport through the water column, we will argue in Sect. 4.1 and App. 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 App.B). One consequence of their fast sinking is the absence of almost any fragmentation after they leave the sea surface (Andrady, From a meta-analysis of 39 previous studies, Erni-Cassola et al. (2019) established the proportion of the most abundant polymertypes 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 notdistinguish between different regions (e.g., coastal region or open water; even inland water bodies of urban environments are
00 1000 1200 1400 1600 1800 2000 2200ρ p [kg/m³]OthersAcrylicPAPESTPSPEPP S e a w a t e r d e n s i t y Most abundant polymer types located at oceans
Polymer Densities
Figure 1.
Polymer densities for the most abundant microplastics identified in water bodies (Erni-Cassola et al., 2019). included in the analysis) and between the particle types (size range and shape) concerned in the different studies. We organizethese polymer types according to their density (Chubarenko et al., 2016; Andrady, 2017; Erni-Cassola et al., 2019) in Fig. 1:PP between − kg/m , PE − kg/m , PS with kg/m (excluding its foamed version), PEST in the range − kg/m , PA within − kg/m and acrylic with kg/m . There is also some less abundant plasticlike polytetrafluoroethylene (PTFE) which has higher densities, in the range − kg/m .Thus, the full range of microplastic particle densities in the ocean, denoted here as ρ p , is − kg/m , and most ofthem have densities within the interval − kg/m . This has to be compared with the seawater density, which close tothe surface has a conventional mean value of ρ f = 1025 km/m (red line in Fig. 1) and changes around from the surfaceto the sea bottom. Since we are interested in sinking material, and for the sake of maximal practicality, we restrict our study tomicroplastics of densities kg/m < ρ f < kg/m .Another relevant property of plastic particles is their size. By a widely accepted definition, microplastics are particles witha diameter less than mm without any lower limit (Arthur et al., 2009). Some observations at the ocean surface show thatthe most common diameter is around mm (Cózar et al., 2014, 2015), with an exponential decay with increasing diame-ter up to mm . However, the absence of this peak in other studies that show an increasing abundance with further de-creasing 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 (whichusually use trawl nets with a mesh size around . mm ) and especially for the adaptation of careful and standardized analysisprocedures to avoid artifacts (Filella, 2015).Field data about distributions of size and quantifiers of shape for negatively buoyant rigid particles in the water columnor deep-sea sediments are not available to date to the best of our knowledge, except in the Artic for Bergmann et al. (2017).However, their results may not apply to the majority of the oceans because of the very special dynamics provided by meltingand freezing of sea ice (Bergmann et al., 2017). Data from Bergmann et al. (2017) and Song et al. (2014) about unspecifiedsedimented fragments and paint particles, respectively, exhibit an increasing abundance with decreasing size, most particles eing smaller than . 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 theobservations of Bagaev et al (2017), Kanhai et al. (2018) and Peng et al. (2018) do not indicate this abundance of smallparticles.For these reasons, we will disregard particles of extremely small size. To keep our qualitative study sufficiently simple, wewill consider all our modelled particles to have a radius a = 0 . mm (a diameter of . mm ). This is a rather small size, butstill within the commonly measured ranges. As we will discuss in Section 4.1, this radius is well within the validity range ofthe MRG equation. In this subsection we indicate the total amount of dense microplastics entering the water column in open waters of the Mediter-ranean. Despite the correlation of plastic source with coastal population density, the rapid fragmentation of small particles alongthe 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 urbanareas or from rivers. The estimations are based on the results of Kaandorp et al. (2020). They provide a total amount of yearlyplastic release into the Mediterranean in the range − tonnes, from which around 37% corresponds to negativelybuoyant 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 tonnes per year, 37% of sinking particles, and the proportion of direct release by maritimeactivity (6%) to obtain in Sect. 4.3 an estimate for the basin-wide yearly release of negatively buoyant sphere-like microplasticsin the open Mediterranean. Note that we choose the upper bound, 4000 tonnes per year, in order to account for the considerableamount of unregistered particles. A standard modeling approach (Siegel and Deuser, 1997; Monroy et al., 2017; Liu et al., 2018; Monroy et al., 2019) for thetransport of noninteracting sinking particles is to consider the time-dependent particle velocity v as the combination of theambient fluid flow u and a settling velocity v s as: v = u + v s , (1)with v s = (1 − β ) g τ p , β = 3 ρ f ρ p + ρ f , and τ p = a βν . (2) g denotes the gravitational acceleration vector, pointing downwards; β is a parameter depending on the particle and the fluiddensities, ρ p and ρ f , respectively. Particles heavier than water have β < , and β = 1 for neutrally buoyant particles. Theexpression given for β assumes spherical particles. τ p is the Stokes time, i.e., the characteristic response time of the particle o changes in the flow, where a is the radius of the particle and ν the kinematic viscosity of the fluid. Although Eq. (1) iscommonly used, we are not aware of a systematic justification of it in the microplastics context. This will be done in Section4.1. For the flow velocity u we use a 3D velocity field from NEMO (Nucleus for European Modelling of the Ocean), which im-plements a horizontal resolution of 1/12 degrees and 75 s-levels in the vertical with updates data every 5 days (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 con-sidered 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 bound-ary conditions. After particles are released at some initial date, in a so-called flash release, they evolve under equations ofmotion such as Eq. (1), and the statistics of the resulting particle cloud are analyzed. We next show, closely following the treatment of Monroy et al. (2017) for biogenic particles, that possible inertial effectsthat would correct Eq. (1) are negligible for the sizes and densities of typical dense microplastics. To this end, similarly tomany other studies (Michaelides, 2003; Balkovsky et al., 2001; Cartwright et al., 2010; Haller and Sapsis, 2008), we start bychoosing 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 validfor 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 d v dt = β D u Dt + u − v + v s τ p . (3)Beyond sphericity, two conditions are needed for the validity of Eq. (3) (Monroy et al., 2017; Maxey and Riley, 1983): a) theparticle radius, a , has to be much smaller than the Kolmogorov length scale η of the flow, which has values in the range . mm < η < mm for wind-driven turbulence in the upper ocean (Jiménez, 1997); b) the particle Reynolds number Re p = a | v − u | ν ≈ av s ν should satisfy Re p ≪ . Note that this last condition imposes restrictions on the values of the particles’ densityand size, partially via the settling velocity v s = | v s | . For the most abundant sinking microplastics, i.e., with densities ρ p =1025 − kg/m , we now determine the range of validity of Eq. (3) assuming ν = 1 . × − m /s and ρ f = 1025 kg/m o be fixed. This gives β in the range . − . 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 Section 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 minimalvalue of the Kolmogorov scale η = 0 . mm with the red line (Jiménez, 1997), and Re p = 1 with a black line, which boundthe area of validity (shaded in the plot). We also indicate v s as a function of a for β = 0 . with the blue curve, correspondingto the upper bound to v s for typical microplastic densities. In total, the zone with soft shading in Fig. 2 represents a parameterregion where Eq. (3) applies for particles with β < . (i.e., particles falling faster than the typical ones), whereas the area of ourinterest, corresponding to β ≥ . , is represented by a dark shading, denoting the typical plastic sizes and corresponding settlingvelocities for which the equation is valid. As a rule of thumb, in a typical situation, validity of Eq. (3) requires v s < . m/s and a < . mm . As discussed in Section 3.1, information about particles in the validity range is particularly sparse for surfacewaters because of the usual sampling techniques, but sediment data indicates 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. These estimates of the Kolmogorov scale anywayassume wind-driven turbulence and are thus restricted to the mixed layer (Jiménez, 1997), below which η is undoubtedlylarger. 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 App. A, we present quantitative arguments for the applicabilityof 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 prop-erties of negatively buoyant rigid microplastics in the water column. Note that rigidity of the particles is an essential conditionwhich is why the advection of microfibers is out of the scope of this paper.The connection between Eq. (3) and its approximation Eq. (1) is made by noticing that τ η ≈ s in the ocean (Monroy et al.,2017; Jiménez, 1997), so that the Stokes number St = τ p /τ η , which measures the importance of particle inertia in a turbulentflow, is very small (of the order of − − − ). Thus an expansion of the MRG equation for small St (smallness of theFroude number, i.e. smallness of fluid accelerations with respect to gravity, is also required) can be performed. The expansionin its simplest form leads to (Balachandar and Eaton, 2010; Monroy et al., 2017; Drótos et al., 2019): v ≈ u + v s + τ p ( β − D u Dt . (4)We can now take the results of Monroy et al. (2017) for biogenic particles of sizes and densities similar to the microplasticsconsidered here to show that the inertial corrections (the term proportional to τ p ) in Eq. (4) are negligible, so that the simplerEq. (1) correctly describes sinking of microplastics in the considered parameter range. For completeness, we report in App. Bthe explicit numerical calculations showing this (in which the influence of the Coriolis force is also taken into account, sinceit is known to be of the same order or larger than the inertial term when a large-scale flow is used for u ). In particular we findfrom release experiments from m below the surface of a large number of particles with β in the range . − in the wholeMediterranean that the difference between horizontal particle positions after 10 days of integration calculated from Eq. (4) and v s [ m / s ] a [mm] K o l m ogo r o v sc a l e β = . R e p = Figure 2.
Settling velocities and particle sizes for which Eq. (3) holds. Kolmogorov scale is represented by the red line and Re p = 1 witha black line, which bound the area of validity. Blue curve corresponds to v s = v s ( β, a ) for β = 0 . , the upper bound to v s for typicalmicroplastic densities. Dark shading denotes the plastic particle sizes and corresponding settling velocities for which application of Eq. (3)is valid. the simpler Eq. (1) is just a 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 limitationof the respective solutions of Eq. (1). The NEMO velocity field of our choice is not an exception, but a rigorous evaluation ofthe corresponding errors of particle trajectories is not possible without knowing the actual small-scale flow. Nevertheless, onecan roughly estimate the effect of these small scales by adding a stochastic term to Eq. (1) with statistical properties similar tothe 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 App. B, indicate that after 10 days of integration the relative difference between particle positionsgiven by Eq. (4) with and without this ‘noise’ term modeling small scales (using β = 0 . ) is around 8% for the horizontaldisplacements and 5% for the vertical ones. The figures become 12% and 5%, respectively, when evolving the particles for 20days. These errors are moderate, although they may be of importance under some circumstances (Nooteboom et al., 2020). Weconsider these figures as a baseline to evaluate corrections to the simple Eq. (1): adding more complex particle-dynamics termsto it will not improve plastic-sedimentation modeling unless the effect of these corrections is significantly larger than the aboveestimations for the effect of the unknown small-scale flow. In the following we consider the simple Eq. (1), but we estimate theimplications of assuming or not a constant value of the water density. In this section we analyze the role of a variable seawater density on the particle settling dynamics. Fluid density is calculatedfrom the TEOS-10 equations, which is a thermodynamically consistent description of seawater properties derived from a Gibbsfunction, for which absolute salinity is used to describe salinity of seawater and conservative temperature replaces potentialtemperature (Pawlowicz, 2010). In the simulations described in this section, as particles move in the ocean they encounterdifferent temperatures and salinities, as given by the NEMO model described in Sect. 3.4, and then they experience differentvalues of the ambient-fluid density.We consider particles of a fixed density ρ p = 1041 . kg/m . This implies that for a nominal water density of ρ f = 1025 kg/m the value of β would be β = 0 . , giving a sinking velocity v s = 6 . m/day for our particles of radius a = 0 . mm , but thissinking velocity will be increased or decreased in places where water density is lower or higher, respectively, so that we havea spatially- and temporally-dependent velocity in Eq. (1). The particle density and size have been chosen to be representativeof the slowly-sinking microplastic particles, for which we expect the seawater density variations to have the largest impact. Inthis way we find some upper bound for the importance of variability in seawater density for particle trajectories.We release N = 78 , particles over the whole Mediterranean Sea, and monitor their trajectories under Eq. (1). The leftpanel of Fig. 3 shows the histogram of water densities encountered by the particles when the release is performed on July 8th2000. On this summer date, the Mediterranean is well stratified, at least in its upper layers. Initially the particles are in surfacewaters with a range of salinities that average approximately to the nominal ρ f = 1025 kg/m . But as they sink in time theyreach layers with higher densities (and more homogeneous across the Mediterranean). When the release is done in winter (rightpanel of Fig. 3) the water column is more mixed, so that the range of water densities experienced by the particles released atdifferent points is always narrow. But the mean water density turns out to be always larger than the conventional surface densityof ρ f = 1025 kg/m , so that a slightly slower sinking is expected to occur.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 ( kg/m ) and position-dependent fluiddensity, respectively, using the same release location and date 8 July 2000 in both cases. Particle density is fixed at ρ p =1041 . kg/m . The difference between the two calculations (and thus the error of considering that constant value for thedensity) should be compared to average horizontal and vertical displacements of km and m , respectively, at t = 20 days. At that time, we thus find that the influence of variable fluid density on the dynamics is about for the horizontalmovement and for the vertical displacement on average.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 displayedin Table 1. The relative error produced by assuming a constant density is larger in the vertical direction. It is also larger for therelease in winter, but this is a consequence of taking a value for the reference density that is not representative of winter watersbut is strongly biased (see Fig. 3, right). If using a reference value more appropriate for winter waters (say ρ f ≈ kg/m )
023 1024 1025 1026 1027 1028 1029ρ f [kg/m³]0.00.51.01.52.02.5 (a) t = 0 dayst = 10 dayst = 20 days 1023 1024 1025 1026 1027 1028 1029ρ f [kg/m³]0.00.51.01.52.02.5 (b) t = 0 dayst = 10 dayst = 20 days Figure 3.
Normalized histogram of the seawater density ρ f at the positions of N = 78803 particles after t = 0 , , and days of beingreleased 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 . kg/m , and fluid density is obtained from the TEOS-10 equations. The vertical line indicatesa conventional seawater density of kg/m . t [days] ⟨ | x ⟨ ) ⟨ t ) − x ⟨ ) ⟨ t ) | ⟩ [ m ] ⟨a) t [days] ⟨ | z ⟨ ) ⟨ t ) − z ⟨ ) ⟨ t ) | ⟩ [ m ] ⟨b) Figure 4.
The distance, as a function of time, between trajectories obtained with constant nominal fluid density of ρ f = 1025 kg/m andthe actual variable fluid density, both starting at the same initial location. The range of the values among all particles released in differentpoints of the Mediterranean is indicated by the shaded area, while the solid line indicates the average over the particles. Particles have ρ p = 1041 . kg/m , and all parameters are the same as for the summer release in Fig. 3.(a): horizontal distances; (b): vertical distances. the relative error remains quite small, due to the weaker stratification of the sea during this season. In fact, the reference valueis also biased in the summer unless the investigation is restricted to the surface. able 1. Relative effect on horizontal and vertical particle positions after 10 and 20 days of integration, averaged over 78803 particles releasedover the whole Mediterranean at m depth, of replacing the actual seawater density by a nominal value ρ f = 1025 kg/m .10 days 20 daysSummer release Horizontal: 1.12 % 2.75 %Vertical: 3.19 % 6.25 %Winter release Horizontal: 1.88 % 5.62 %Vertical: 8.14 % 9.32 % In brief, we see that the effect of location-dependent density may be a relevant effect to evaluate microplastic transport. Atleast, 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 densityon particles that sink faster will be smaller. Also, the effects reported in Table 1 remain of the order of the estimations of theeffects of unresolved small scales of the flow (Sect. 4.1). As a consequence, in the following we will not consider variableseawater density, but restrict our modeling to Eq. (1) with a constant nominal value of the sinking velocity v s . We will first estimate the total mass of negatively buoyant rigid microplastics in the water column of the open MediterraneanSea by assuming a uniform vertical distribution, then we will justify this assumption by running numerical simulations accord-ing to the conclusion of Section 4.1 about the equation of motion.For estimating the total mass, we take the quantities of Section 3.2 ( tonnes/year of plastic release, with 37% beingnegatively buoyant of which 6% originates from maritime activities) to compute the rate r at which microplastic particlesof our interest enter the water column in the open sea: r = 4000 tonnes/year × . × .
06 = 88 . tonnes/year , or r =0 . tonnes/day .The next step is to estimate the time during which these microplastic particles remain in the water column before reaching thesea 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 sink-ing velocity, τ = h/v 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 mi-croplastic particles of radius a ≈ . mm the range of sinking velocities is . − . m/day , giving a residence timein the range . − days . Averaging these times weighted by the proportion of each type of plastic we get τ ≈ days .Combining the input rate r with this mean residence time we get an estimate for the total amount present in the water columnat any given time as Q = rτ : the result is Q ≈ . tonnes of dense rigid microplastics if all of them would be in the form of articles of size a = 0 . mm . This is below but close to 1% of the estimated upper bound of tonnes of floating plastic inthe Mediterranean (according to the corresponding estimation of Kaandorp et al. (2020)).We emphasize the many uncertainties affecting this result (Sections 3.1 and 3.2), and we highlight the one related to particlesize: because of the quadratic dependence of the sinking velocity on the particle radius a , Eq. (2), choosing the particle size to behalf of the one used here will lead to a four times larger estimate for the mass if the same release rate is assumed. This enhancedretention of smaller particles in the water column may imply, depending on the actual size distribution, a dominance of verysmall particles on the plastic mass content of the water column. However, our estimates of plastic input into the ocean (we usemainly Kaandorp et al. (2020)) rely on observations that do not catch extremely small particles. These considerations furtherjustify our choice of a radius a = 0 . mm , small but still easily detectable, as convenient to provide reasonable estimationsof negatively buoyant rigid microplastic mass in the water column within commonly quoted size ranges. We can not excludelarger plastic content at smaller sizes. Another source of bias may be not considering in this study the impact of small-scaleturbulence and convective mixing events. While small-scale turbulence might cause an increase of lifetimes of particles in thewater column, dense water formation and rapid convection, a process reported in areas such as the Gulf of Lions, might likelyreduce particle retention time. These events take place in winter and were shown to transfer particles from the ocean surface tomid-waters (1000 meters) and deep ocean (>2000 meters) in a very short time (1-2 days) and lead to the formation of bottomnepheloid 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, weassume a uniform particle distribution over the whole Mediterranean both in horizontal and in vertical. Since the volume of theMediterranean is about . × km (Eakins and Sharman, 2010) the estimated density would be ρ V ≈ . × − 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 doesnot address coastal areas, where the density would likely be higher.The above estimates are rather rough as a result of the mentioned uncertainties. The assumption of a uniform distributionin the vertical direction has not yet been justified either, but we will show it to be appropriate by means of our simulations ofparticle release starting at m depth over the whole Mediterranean. Instead of performing a continuous release of new particlesat each time step, and computing statistics over this growing number of sinking particles, we approximate this by the statisticsof all positions at all time steps of a set of particles deployed in a single release event. This assumes a time-independent fluidflow, but this approximation is appropriate, since the dispersion of an ensemble of particles released in a single event followsrather well-defined statistical laws, see Section 4.4, and is thus independent of the time-varying details of the flow. Particlesare removed when touching the bottom. For our estimate, we use β = 0 . , i.e., assuming the fastest sinking velocity of typicalplastic particles, for which particles reach vertical depths deeper if compared with the slower sinking velocities used in ourstudy.Figure 5 shows ρ ( z ) , the density of plastic particles per unit of depth z in the whole Mediterranean, and also A ( z ) , theamount of area that the Mediterranean has at each depth z . Both functions have been normalized such that the value of theirintegrals with respect to z is one; in this way the functions can be displayed in the same plot. We see that both curves are ρ ( z ) , A ( z ) z [m] ρ ( z ) / A ( z ) z [m] Figure 5.
The continuous line is the area that the Mediterranean has at each depth z . The dashed line is microplastic density per unit ofdepth ρ ( z ) under continuous release of particles with β = 0 . at m depth. Both curves have been normalized to have unit area, so that theycan be compared on the same scale. The binning size is m . The inset shows the ratio ρ ( z ) /A ( z ) , proportional to the mass density ofmicroplastic per unit of volume ρ V ( z ) . nearly identical (in fact, just proportional, because of the normalization), indicating that the variation of the number of particleswith depth is essentially due to the decrease of sea area with depth. A clearer way to see that is to plot ρ ( z ) /A ( z ) , whichis proportional to the mean plastic concentration per unit volume at each depth z , ρ V ( z ) . We see that this quantity is nearlyconstant, at least in the first 3000 m. At larger depths a weak increase seems to occur, but made unclear by the poor statisticsarising from the small area and number of particles present at these depths. Thus, the hypothesis of a uniform distribution ofplastic 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 a constant(since each particle will spend exactly the same time at each depth interval). The equation of motion used, Eq. (1) correctsthis constant sinking velocity v s with a contribution u from the ambient flow. Thus, the close-to-constant character of theplastic concentration may imply that the flow correction u is negligible, at least when considering its effect over the wholeMediterranean. Another possibility is that the fluctuating flow component u in Eq. (1) results in a vertical dispersion compatiblewith a constant concentration. Although the former explanation predicts an alteration from a constant if the settling velocityis sufficiently small to allow u to induce a stronger vertical dispersion, we will see in the next section that a nearly constantconcentration may be assumed for the majority of our parameter range. We now analyze in detail the transient evolution of particle clouds initialized by flash releases at a fixed depth. Numericallywe proceed by releasing N = 78803 particles uniformly distributed over the entire Mediterranean surface at m depth in thewinter season, as already described. They evolve according to Eq. (1) using a constant water density. We take three examples for he particle density, which correspond to β = 0 . , . , . , or v s = 153 . , . , . m / day for our particles of radius a =0 . mm , respectively; in what follows, these setups shall be denominated as v153, v68 and v6. The horizontal displacementsduring the particle sinking times are much larger (of the order of km ) than the sea depth, so that in fact the particles aresinking sideways (Siegel and Deuser, 1997). However, the horizontal displacements still remain very small compared to thebasin size, and we concentrate on the vertical motion. Even though the vertical steady distribution has been found to be closeto uniform in Section 4.3, the reason for this is not evident, and we will give support here for the pertinence of this finding tomost 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 forv6). The plot is given in terms of a rescaled variable ˜ z = z − tv s σ z where σ z ≡ h ( z i − h z i i ) i is the variance of the particles’ z coordinate. Here the subindex i refers to the particle and h . . . i denotes averaging over different particles. Thus, we plot in thefigure the rescaled distribution of the particles around the average depth of the particles at any given time. For comparison, thenormal distribution is plotted with dashed lines. This figure shows deviations from Gaussianity for early times. The deviationfrom normal distribution decreases for later instants but remains considerable, especially for the tails, which may also beindicative of anomalous diffusive behavior. For reference, particles reach the mean Mediterranean depth, h = 1480 m at times τ = 9 . , . and . days for v153, v68 and v6, respectively.Since a non-Gaussian distribution is usually linked to anomalous dispersion (Neufeld and Hernández-García, 2009), we nowanalyze this aspect in detail by considering how the variance of the vertical particle distribution, σ z ( t ) , evolves. Although thereis a continual loss of particles because of reaching the seafloor with a varied topography, we illustrate in App. C that ourconclusions 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 bythe approximate effective exponents ν , defined through approximate behaviors σ z ∼ t ν in different time intervals.We start our analysis with the fastest-sinking particles (v153, Fig. 7a). At the very beginning, superdiffusion takes placewith ν > , which may be related to autocorrelation in the flow, but we will iterate on this question when comparing differentsettling velocities. Around t = 1 day, the evolution seems to become consistent with normal diffusion ( ν = 1 ), usual afterinitial transients in oceanic turbulence (Berloff and McWilliams, 2002; Reynolds, 2002). However, around t = 4 . days, wecan 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 z i = h z i i + ¯ ω i t + W i , (5)where h . . . i denotes, as before, an averaging over different particles. Here we are assuming that z i − h z i i evolves according tothe sum of a constant average velocity contribution ¯ ω i for sufficiently long times (a characteristic of the flow region traversedby particle i ), and of W i , a Wiener process representing fluctuations with zero mean and defining a diffusion coefficient D i foreach trajectory by W i = D i t . The overbar refers to temporal averaging for asymptotically long times along the trajectory of agiven particle (but assuming that the particle remains in a region with a well-defined ¯ ω i = 0 ), and D i characterizes the strength pd f (z-tv s ) / σ z [m] t=1.16t=4.16t=10.12 0.001 0.01 0.1 1 10-8 -6 -4 -2 0 2 4 6 8(b) v68 pd f (z-tv s ) / σ z [m] t=1.16t=4.16t=10.12t=20.16t=29.8 0.001 0.01 0.1 1 10-8 -6 -4 -2 0 2 4 6 8(c) v6 pd f (z-tv s ) / σ z [m] t=1.16t=4.16t=10.12t=20.16t=29.8 Figure 6.
The probability density function, estimated from a histogram of bin size . , of all particles released in the Mediterranean in therescaled variable ˜ z = z − tv s σ z for the different setups (v153, v68 and v6) and times (in days) as indicated. For comparison, normal distributionsof zero mean and unit variance are shown with a dashed line. of the fluctuations. Assuming h ¯ ω i W i i = 0 , σ z ≡ h ( z i − h z i i ) i = h ¯ ω i i t + h D i i t, (6) σ z [ m ] t [day](a) v153 0.01 0.1 1 10 100 1000 0.1 1 10 σ z [ m ] t [day](b) v68 0.001 0.01 0.1 1 10 100 1000 0.1 1 10 σ z [ m ] t [day](c) v6 10 100 10 t [day] Figure 7.
Variance of depth reached by the particles as a function of time. Straight lines represent power laws for reference, with exponents (in green, corresponding to standard diffusion) and (in purple, corresponding to ballistic dispersion). that is, the variance is a sum of a ballistic and a normal diffusive term, associated with regional differences in the mean velocityand with fluctuations, and dominating for long and short times, respectively. Writing D = h D i i , the crossover between the two
30 35 40 45 0 10 20 30 ϕ [ d e g r ee ] λ [degree] ω i [m/day]-1 0 1 Figure 8. ¯ ω i as estimated from t = 10 . days plotted at the initial position of each particle i in the v153 simulation. The black rectangle inthe Western Mediterranean is the area of large depth considered in App. C. regimes is obtained by equating the two terms as t ∗ = D h ¯ ω i i . (7)To evaluate Eq. (7), we first estimate ¯ ω i for each particle from the “asymptotically” long time of t = 10 . days, which is thelatest time after the crossover still in the ballistic regime in Fig. 7a, when the contribution of fluctuations should already havebecome negligible. The horizontal pattern of the estimated ¯ ω i is presented in Fig. 8, which confirms its patchiness throughoutthe Mediterranean, associated with mesoscale features. Computing h ¯ ω i i and fitting a line to σ z ( t ) between t = 1 . and daysto estimate D , we obtain t ∗ ≈ . days from Eq. (7), which remarkably agrees with Fig. 7a. After approximately t = 12 daysthere is hardly any dispersion, since most of the particles are close to the sea bottom (cf. Fig. 9) where the vertical fluid velocityis nearly zero. Note also a small drop in σ z ( at the very end of the time series, where the results may actually be subject toartifacts, see App. C. However, this is of minor importance, since the distribution of particles so close to the bottom shouldanyway be strongly influenced by resuspension and remixing by bottom currents (Kane et al., 2020).The different regimes are not as clear in the v68 case as for v153, see Fig. 7b. One evident novel feature is a subdiffusiveregime during the transient from the initial superdiffusion (as in the case of horizontal tracer dispersion in the ocean studiedby Berloff and McWilliams (2002); Reynolds (2002)). Approximate normal diffusion is then observed until t = 10 days, whena crossover to a faster dispersion does seem to take place, see the inset. A fit of normal diffusion from t = 4 to days and thevelocity variance at t = 12 . days give an estimate t ∗ ≈ . . However, the long-time ballistic regime is not clear. In fact, along-term return from such a ballistic regime to normal diffusion is expected as a result of increasing horizontal mixing, whichrenders ¯ ω i time dependent and makes it approach zero. According to a careful visual inspection of the inset in Fig. 7b, this maytake place already around t = 14 days.For v6 (Fig. 7c), the transition from the initial superdiffusive regime to that of normal diffusion appears to not involvesubdiffusion. This is already informative: the fluid velocity field is the same for the three simulations with different settlingvelocity, so the differences must originate from the different rate of sampling of the different fluid layers by particles whilethey sink. In particular, the decay of autocorrelation is obviously faster for faster-sinking particles, since it is determined by thespatial structure of the velocity field. σ z [ m ] 〈 z i 〉 [m] v153v68v6 Figure 9.
Variance of depth reached by the particles as a function of their mean depth.
While this is one possible explanation for the earlier timing of the initial transition from anomalous to normal diffusion forhigher settling velocity, one cannot exclude that a depth-dependent organization of the flow is more in play; note that ν > atthe beginning, which might not be explained by simple autocorrelation but might be characteristic of properties of the velocityfield at those depths. The governing role of the spatial structure is supported by Fig. 9: the transition in question takes placeat the same depth ( ≈ m) in the different simulations, which seems to point to mixed-layer processes. Depth-dependencemight also govern the suppression of ballistic dispersion for long times, but it is very unclear.Note in Fig. 9 that the vertical variance is not expected to grow much larger for v68 than for v153 even if the simulation werelonger. Therefore, even if the constancy of the steady vertical distribution relies on the weak vertical dispersion for v153 (seeSection 4.3), constancy is expected to hold in most of our parameter range. A considerably stronger dispersion and a possiblecorresponding deviation from constancy may arise only for extremely low settling velocities, like for v6 in Fig. 9. We have discussed the different types of plastics occurring in the water column, pointing out gaps in our knowledge about thesources, transport pathways and properties of such particles. It would be highly beneficial to have distributions of size, polymertype 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 simplifiedMRG 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 inertialterms are negligible. When a velocity field of large-scale nature is input to the equation (such that small-scale turbulence isnot 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 obust features of the large-scale flow and must remain unaffected by moderate errors. We also note that the traditional valueof seawater density, ρ f = 1025 kg/m , is representative only for near-surface layers in the summer, and correcting for thebias could reduce the error of simulations with a constant seawater density. A suitable equation of motion for the particlesconsidered is constructed by adding to the 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 anearly uniform steady distribution along the water column, perhaps except at extremely low settling velocities. The correspond-ing total amount of plastic present in the water column is relatively small, close to 1% of the floating plastic mass, but it maybe 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 plasticamount 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 microplasticparticles implies that advection by the flow may contribute to large-scale horizontal inhomogeneity of deep-sea plastic sedi-ments by means of recently described noninertial mechanisms (Drótos et al., 2019; Monroy et al., 2019; Sozza et al., 2020).This may especially be 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 seabottom.As for the vertical distribution profile, its approximate uniformity may be linked to the weak vertical dispersion of particlesthat is found in our simulations started with a flash release over the whole surface of the Mediterranean sea. The shape of theemerging transient vertical distribution exhibits deviations from a Gaussian, which are related to anomalous diffusive laws thatdominate the vertical dispersion process in some phases.The different diffusive laws are related to the properties of the decay in the Lagrangian velocity autocorrelation definedalong the trajectories of the sinking particles. An important example is the transition from initial superdiffusion to a longerphase of normal diffusion, occurring around 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 asthe particles’ horizontal displacements remain small. The vertical diffusion returns to the normal type when horizontal mixingbecomes more developed. These results suggest regional differences in the sinking process, so that regional modeling might bemore appropriate than a whole-basin approach. Future studies will include different areas of the oceans, and analyze the roleof Lagrangian coherent structures on the different vertical dispersion regimes. Data availability.
The velocity field from the NEMO simulation used in our study can be downloaded from http://opendap4gws.jasmin.ac.uk/thredds/nemo/root/nemo_scan_catalog.html.Parcels is a set of Python classes developed under the TOPIOS project and is accessible at https://oceanparcels.org/. The scripts for runningthe particle transport simulations are available upon request from Rebeca de la Fuente, rebeca@ifisc.uib-csic.es.
We quantitatively assess the impact of deviations from a spherical shape through a correction to the settling velocity v s . Thesimplified MRG equation, Eq. (3), or its first-order approximations in the Stokes number, Eqs. (4) and (B1), are affected byparticle geometry through the drag force and the added mass term; however, accelerations are irrelevant for v s , so that theadded mass term does not appear in its formulation or in the simple approximation of Eq. (1). We will compare values of thesettling velocity describing nonspherical and spherical particles with the same density, then finally 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 thedifference of the particle and the fluid velocities, v and u , respectively) with the resultant of gravitational and buoyancy forces: F drag ( v − u ) + V ( ρ p − ρ f ) g (A1)with v − u = v 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 F (sph)drag ( v − u ) = − πµ ( v − u ) a, (A2)where µ is the dynamical viscosity of the fluid. According to Leith (1987); Ganser (1993), an appropriate approximation forsmall nonspherical particles is F (non)drag ( v − u ) = − πµ ( v − u ) (cid:18) a n + 23 a s (cid:19) , (A3)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 isobtained by substituting v − u = v s , and solving Eq. (A1) for v s . We denote the magnitudes of the settling velocities obtainedfrom Eq. (A2) and Eq. (A3) by v (sph) s and v (non) s , respectively.To characterize the correction in the settling velocity for a given nonspherical particle (with a given density ρ p ) with respectto assuming a spherical shape with a radius a , we will consider q ≡ v (non) s v (sph) s = 34 π V (non) a (cid:0) a n + a s (cid:1) , (A4)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 ofthe 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 devi-ation from a spherical shape. In this sense, any choice of a representing a spherical shape, including ours in the manuscript,describes the settling velocity of certain nonspherical particles, the question is just their shape and size, which will mutuallydepend on each other for a given a . We will nevertheless proceed by choosing a shape class and defining a along indepen-dent considerations, because we intend to link a given a to a single particle size as identified during the processing of fieldobservations. he shape of rigid microplastic particles is not usually described in the literature, but we can see photographs of someexamples in, e.g., Song et al. (2014); Fischer et al. (2015); Bagaev et al. (2017). For an explorative computation, a reasonablechoice seems to be a rectangular cuboid with edges A , B = ˆ BA < A and C = ˆ CA < B < A , where one or both of ˆ B and ˆ C are less than but greater than, say, . .Under this assumption, the particle size will correspond to the longest edge, A , of the cuboid if the size is identified throughmicroscopy as the largest extension (“length”; e.g., Cózar et al., 2015); and it may be related more to the middle edge, B , ifone thinks of a sieving technique (e.g., Suaria et al., 2016). The naive choice will be a = A/ and a = B/ 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 , a s = (cid:18) AB + AC + BC π (cid:19) . (A5)However, a n depends on the particle’s orientation with respect to the relative velocity. Implications will be discussed wheninterpreting 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 a n ,X = (cid:18) ABCπX (cid:19) (A6)for X = A , B and C .After substituting all these expressions in Eq. (A4), we obtain q ( A/ X = 9 π − ˆ B ˆ C ˆ B ˆ C ˆ X ! + 2 (cid:16) ˆ B + ˆ C + ˆ B ˆ C (cid:17) − , (A7) q ( B/ X = ˆ B − q ( A/ X , (A8)where X = ˆ XA has been introduced.We plot q ( A/ A in Fig. A1 as a function of ˆ B and ˆ C for ˆ B, ˆ C ∈ [0 . , with ˆ B > ˆ C . Its range extends from . to . on thisdomain, but it drops below . only for ˆ C very close to . and ˆ B below . ; i.e., for extremely thin rod-like particles, whichdo not appear to be common based on photographs. The range of q ( B/ A (not shown) on the same domain is between . and , and values above are again restricted to very small ˆ C and to ˆ B < . . The results are very similar for other choices of X ,deviations beyond with respect to X = A are only found for small ˆ C and do not reach beyond even there.We have left the question which orientation is relevant open so far. In small-scale isotropic turbulence, which is certainlypresent in the ocean, nonspherical particles have a preferential alignment with certain characteristics of the flow but undergorotation (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 describingshapes. The only relevant exception is small ˆ C with ˆ B ≈ . This regime may characterize paint flakes (Song et al., 2014;Bagaev et al., 2017), but the relative difference remains below even there. q A ( A /2) q A ( A /2) B ˆ C ˆ q A ( A /2) Figure A1. q ( A/ A as a function of ˆ B and ˆ C for ˆ B, ˆ C ∈ [0 . , with ˆ B > ˆ C . Dotted lines represent q ( A/ A = 0 . . Even though the real advection of the particles will become more complicated as a result of the ever-changing orientationand may thus be beyond the scope of the MRG equation (cf. the discussion in the main text about the settling velocity ofirregular 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 shapesunder the assumption of the same particle density), this gives quantitative support for the applicability of a spherical shape inEq. (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 inertialterm through added mass. Since corrections in added mass with respect to a sphere are of order for all common shapes (seeKaneko et al., 2014, for an overview), we believe that the finding of App. B about the negligible importance of inertial effectsin these equations is not affected by a deviation from sphericity. The Stokes number, which is proportional to the settlingvelocity and also depends on the coefficient of added mass ( St ∼ τ p with τ p given by Eq. (2)), is estimated to be − − − for spherical particles in Sect. 4.1, so that it will not increase to due to a nonspherical shape either, hence leaving theapproximation (4) of Eq. (3) valid. We present here the detailed numerical analysis of the relevance of a finite time of response (Stokes time, τ p ) of the particle tothe fluid forces, the Coriolis force, and scales unresolved by the NEMO velocity field.We incorporate the first two effects to a single equation, v = u + v s + τ p ( β − (cid:18) D u Dt + 2 Ω × u (cid:19) , (B1)which is identical to Eq. (4) except for the addition of the Coriolis force, Ω × u . Ω is Earth’s angular velocity vector. Weinclude the Coriolis force because it can be more important than the other inertial term, given by the fluid acceleration D u /Dt ,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 theequation of motion, Eq. (1), by adding a stochastic noise term: v = u + v s + W . (B2) W ( t ) = ( √ D h ξ x ( t ) , √ D h ξ y ( t ) , √ D v ξ z ( t )) , where ξ ( t ) is a vector Gaussian white noise process (independent for eachparticle) of zero mean and with correlations given by h ξ i ( t ) ξ j ( t ) i = δ ij δ ( t − t ) , i, j = x, y, z . Thus, the horizontal andvertical intensities of this term are given by D h by D v , respectively.The statistical properties are chosen to be similar to the ones expected for oceanic motions below the scales resolved by thenumerical 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 h ( ℓ ) = 2 . × − ℓ . m /s .Taking for ℓ the horizontal resolution of our numerical model (1/12 degrees) we obtain D h = 7 . m /s . Since the Okuboformula is an empirical fit to surface motions, and effective horizontal diffusivity should be weaker below the thermocline, ourresults provide an upper bound for the error associated with unresolved scales of fluid motion. For the vertical diffusivity wetake D v = 10 − m /s .In order to compare the different equations of motion, we release a large number N = 78803 of particles on the wholeMediterranean at m depth on 8 January 2000. We associate a β parameter to each particle by selecting it from a randomuniform distribution in the range β ∈ [0 . , , and once it is selected it remains fixed at all times for the corresponding particle.High values of β (close to one) correspond to more buoyant plastic particles whereas low values of β correspond to highsettling velocities. The corresponding range of velocities is v s ∈ [1 . × − , m/s . We integrate the particle trajectoriesusing Eq. (1) and also, from the same initial conditions, using the corrected dynamics in Eq. (B1) or (B2), all with the sameNEMO velocity field u . The horizontal and vertical distances between particles released from the same point but evolved withdifferent equations are compared by calculating the following averages over particles: d Ih ( t ) = 1 N N X k =1 | r k ( t ) − r Ik ( t ) | , d Iv ( t ) = 1 N N X k =1 | z k ( t ) − z Ik ( t ) | . (B3) r k = ( x k , y k ) is the horizontal position of particle k , z k is its vertical position, and the superindices and I indicate thatthe particle trajectory has been integrated by using the reference equation, Eq. (1) or the one containing inertial corrections, .0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0t [days]020406080100120140160 ⟨ | x ⟨ ) ⟨ t ) − x ⟨ ) ⟨ t ) | ⟩ [ m ] ⟨a) ⟨ | z ⟨ ) ⟨ t ) − z ⟨ ) ⟨ t ) | ⟩ [ m ] ⟨⟩) Figure B1.
Solid line indicates the average horizontal distance d Ih (a) and the average vertical distance d Iv (b) of particles released from thesame initial location but integrated with equations (1) and (B1) in the NEMO velocity field. Shaded region indicates the range of the distancesamong the individual pairs of particles. ⟨ | x ⟨ ) ⟨ t ) − x ⟨ ) ⟨ t ) | ⟩ [ m ] ⟨a) ⟨ | z ⟨ ) ⟨ t ) − z ⟨ ) ⟨ t ) | ⟩ [ m ] ⟨⟩) Figure B2.
Same as Fig. B1 for the comparison of equations (1) and (B2).
Eq. (B1). The quantities d Wh ( t ) and d Wv ( t ) , comparing Eq. (1) with the dynamics (B2) modeling small-scale flow effects, aredefined analogously.In Fig. B1, we display the average distances d Ih ( t ) and d Iv ( t ) as a function of time, characterizing the corrections by inertialterms to the simple dynamics of Eq. (1). Analogously, Fig. B2 displays the average distances d Wh ( t ) and d Wv ( t ) as a functionof time, characterizing the estimated corrections arising from small scales unresolved by the NEMO velocity field. The effectinduced by the inertial terms is very small and clearly negligible. The impact of W , and thus of small unresolved scales islarger. able B1. Average horizontal and vertical pairwise particle distances d and single-particle displacements r after an integration time of days. See text. d/r is indicated by percentages in parentheses. β v s ( − m/s) d Ih (m) d Iv (m) d Wh (m) d Wv (m) r h (m) r v (m) [0 . , .
85) [1 . , .
60 (0.26%) 0.016 (0.001%) 2675 (11.8%) 2.1 (0.16%) 22601 1291 [0 . , .
9) [1 . , .
42 (0.14%) 0.010 (0.001%) 2831 (9.7%) 2.0 (0.24%) 29278 870 [0 . , .
95) [0 . , .
31 (0.08%) 0.009 (0.002%) 3313 (8.0%) 2.1 (0.43%) 41587 493 [0 . ,
1) [0 . ,
16 (0.03%) 0.008 (0.005%) 4668 (8.0%) 2.7 (1.79%) 58320 151 ⟨ | x ⟨ ) ⟨ t ) − x ⟨ ) ⟨ t ) | ⟩ [ m ] ⟨a) β ∈ [0.8⟩0.85)β ∈ [0.85⟩0.9)β ∈ [0.9⟩0.95)β ∈ [0.95⟩1] ⟨ | z ⟨ ) ⟨ t ) − z ⟨ ) ⟨ t ) | ⟩ ⟨b) β ∈ [0.8,0.85)β ∈ [0.85,0.9)β ∈ [0.9,0.95)β ∈ [0.95,1] Figure B3.
Average horizontal distance d Ih (a) and average vertical distance d Iv (b) of particles released from the same initial location butintegrated with equations (1) and (B1) in the NEMO velocity field. The different lines are obtained from particles from different ranges ofdensities characterized by the indicated ranges in β . To evaluate the different effects more quantitatively, we summarize in Table B1, considering particles separately in differentdensity 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 days . To fully appreciate the importance of these numbers, in thetwo 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 andan absolute sense for any realistic application, more care has to be taken with regards to the unresolved scales. Although thevertical 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 biggestfor the fastest-sinking particles and is well above for them.We can also observe the time evolution of the d entries of this table in Figs. B3 and B4. The overall effect of the unresolvedscales is confirmed to be much larger than that of the inertial terms, and the differences between particles with different densities(ranges of β ) are less noticeable (except for the smallest densities considered, i.e. β ≈ ). .0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0t [days]0200040006000800010000 ⟨ | x ⟨ ) ⟨ t ) − x ⟨ ) ⟨ t ) | ⟩ [ m ] ⟨a) β ∈ [0.8⟩0.85)β ∈ [0.85⟩0.9)β ∈ [0.9⟩0.95)β ∈ [0.95⟩1] ⟨ | z ⟨ ) ⟨ t ) − z ⟨ ) ⟨ t ) | ⟩ ⟨b) β ∈ [0.8,0.85)β ∈ [0.85,0.9)β ∈ [0.9,0.95)β ∈ [0.95,1] Figure B4.
Same as Fig. B3 for the comparison of equations (1) and (B2).
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 inSect. 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 bottomboundary (i.e., the seafloor) extends from the surface (close to the coast) to the deepest point of the basin, it is thus relevantat 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 havealready been lost, (iii) decrease in the geographical area sampled by the particles.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, and that computed over a restricted set of particles. This restricted set contains only thoseparticles at the positions of which the bathymetry Z satisfies Z > h z i i− σ z , where the average h z i 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 issufficiently far for its effect on the particle distribution to be negligible at the positions of all particles kept for the computation.Fig. C1 shows that the difference in the results between the full set and the restricted one is negligible for all three settlingvelocities 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, thedistribution of particles so close to the bottom (cf. Fig. 9) should anyway be strongly influenced by resuspension and remixingby bottom currents (Kane et al., 2020). For the rest of the time evolution, we can be confident that a possible distortion of theparticle distribution induced by the boundaries has no impact on the variance curves. Such distortion may be present, but thenumber 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 on the originallycomputed variance. This is a result of the relatively weak dispersion; artifacts might be found for broader distributions, possibly σ z [ m ] (a) v153 whole Mediterraneanadaptive restriction N t [day] 0.01 0.1 1 10 100 1000 σ z [ m ] (b) v68 0 50000 0.1 1 10 N t [day] 0.001 0.01 0.1 1 10 100 1000 σ z [ m ] (c) v6 0 50000 0.1 1 10 N t [day] Figure C1.
Variance of depth reached by all sinking particles (in black) and an adaptively restricted subset of them (in yellow) as a functionof time. See text for details. Straight lines represent power laws for reference, with exponents and . The number N of particles consideredfor the computation of the variance is also shown. ncluding a hypothetical continuation of the v6 simulation. The time evolution of N in Fig. C1 also assures us that poor-qualitystatistics, 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 (withrespect to the population variance) very low. This is so because, under standard assumptions, the ratio between the sample andthe population variance should be close to a chi-square random variable with N − degrees of freedom, divided by N − (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 numberof particles during the simulations that could result in the changes in the slope of the curves of variance versus time. To furthersupport this conclusion, we compute the time evolution of the variance over the particles initialized in a subregion of thewhole Mediterranean, see in Fig. 8. In particular, we choose the box of longitudes [5 , . degrees, and latitudes [37 . , degrees, corresponding to the Sea of Sardinia, where the bathymetry is deep enough to prevent particles from reaching theseafloor (except for the very last few time steps of the v153 configuration), so that the horizontal area sampled by the particlesapproximately remains constant (remember that horizontal displacements are small compared to geographical features).According to Fig. C2, the character of the dispersion in the Sea of Sardinia is nearly identical to that in the whole Mediter-ranean, except after the crossover from normal diffusive to ballistic dispersion in the v153 case (Fig. C2a). The smaller variancein the velocity should naturally lead to a later crossover to the long-time behavior in the Sea of Sardinia (see Section 4.4), itis nevertheless doubtful that the crossover should fall outside the simulation time. As the crossover is not observed, one mightspeculate that horizontal mixing might just get strong enough to suppress ballistic dispersion, similarly to the v68 case (seethe discussion of Fig. 7), for which the time evolution of the variance is very similar in the Sea of Sardinia and the wholeMediterranean (Fig. C2b). Fig. 8 suggests that the characteristic patch size is smaller in the western basin of the Mediter-ranean, including the Sea of Sardinia, than in the (considerably larger) eastern one, which makes inter-patch mixing easier. Weconclude that the only substantial difference between the dispersion in the Sea of Sardinia and the whole Mediterranean mayoriginate from the smaller extension and some special characteristics of the former, and the decrease in the area sampled in thewhole-Mediterranean simulation presumably has not effect on the results.Based on these analyses, we believe that the findings of Section 4.4 are unaffected by the boundary and thus have generalrelevance for oceanic dispersion. Author contributions.
All authors designed research. RF performed the simulations. RF and GD analyzed data. RF, GD, EHG and CLprepared the first draft, and all authors reviewed and edited the manuscript.
Competing interests.
The authors declare that they have no conflict of interest. σ z [ m ] t [day](a) v153 whole MediterraneanSea of Sardinia σ z [ m ] t [day](b) v68 0.001 0.01 0.1 1 10 100 1000 0.1 1 10 σ z [ m ] t [day](c) v6 Figure C2.
Variance of depth reached by all sinking particles (in black) and those initialized in the Sea of Sardinia (in blue; longitudesin [5 , . degrees, and latitudes in [37 . , degrees) as a function of time. See text for details. Straight lines represent power laws forreference, with exponents and . cknowledgements. This work is part of the “Tracking of Plastics in Our Seas” (TOPIOS) project, funded by the European Research Councilunder the European Union’s Horizon 2020 research and innovation program (grant agreement no. 715386). We acknowledge publication feesupport from the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI). R.F.,G.D., E.H-G and C.L acknowledge financial support from the Spanish State Research Agency through the María de Maeztu Program forUnits of Excellence in R&D (MDM-2017-0711). R.F. also acknowledges the fellowship no. BES-2016-078416 under the FPI program ofMINECO, Spain. G.D. also acknowledges financial support from the European Social Fund through the fellowship no. PD/020/2018 underthe postdoctoral program of CAIB, Spain, and from NKFIH, Hungary (grant agreement no. NKFI-124256). ≥ µ m in the Atlantic Ocean and their modelled vertical distribution, Marine Pollution Bulletin, 100, 70 – 81,https://doi.org/https://doi.org/10.1016/j.marpolbul.2015.09.027, 2015.Erni-Cassola, G., Gibson, M. I., Thompson, R. C., and Christie-Oleza, J. A.: Lost, but Found with Nile Red: A Novel Method for Detectingand Quantifying Small Microplastics (1 mm to 20 µ m) in Environmental Samples, Environmental Science & Technology, 51, 13 641–13 648, https://doi.org/10.1021/acs.est.7b04512, pMID: 29112813, 2017.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, Journal of Hazardous Materials, 369, 691–698, 2019.Filella, M.: Questions of size and numbers in environmental research on microplastics: methodological and conceptual aspects, Environmen-tal Chemistry, 12, 527–538, https://doi.org/10.1071/EN15012, 2015.Fischer, V., Elsner, N. O., Brenke, N., Schwabe, E., and Brandt, A.: Plastic pollution of the Kuril–Kamchatka Trench area (NW pacific),Deep Sea Research Part II: Topical Studies in Oceanography, 111, 399 – 405, https://doi.org/https://doi.org/10.1016/j.dsr2.2014.08.012,2015.Gago, J., Carretero, O., Filgueiras, A., and Viñas, L.: Synthetic microfibers in the marine environment: A review on their occurrence inseawater and sediments, Marine Pollution Bulletin, 127, 365 – 376, https://doi.org/https://doi.org/10.1016/j.marpolbul.2017.11.070, 2018.Ganser, G. H.: A rational approach to drag prediction of spherical and nonspherical particles, Powder Technology, 77, 143 – 152,https://doi.org/https://doi.org/10.1016/0032-5910(93)80051-B, 1993. ozza, A., Drotos, G., Hernández-García, E., and López, C.: Accumulated densities of sedimenting particles in turbulent flows, Physics ofFluids, 32, http://dx.doi.org/10.1063/5.0003614, 2020.Stabholz, M., Durrieu de Madron, X., Canals, M., Khripounoff, A., Taupier-Letage, I., Testor, P., Heussner, S., Kerhervé, P., Delsaut, N.,Houpert, L., et al.: 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.Stolle, C., Nagel, K., Labrenz, M., and Jürgens, K.: Succession of the sea-surface microlayer in the coastal Baltic Sea under natural and ex-perimentally induced low-wind conditions, Biogeosciences, 7, 2975–2988, https://doi.org/https://doi.org/10.5194/bg-7-2975-2010, 2010.Suaria, G., Avio, C. G., Mineo, A., Lattin, G. L., Magaldi, M. G., Belmonte, G., Moore, C. J., Regoli, F., and Aliani, S.: The MediterraneanPlastic Soup: synthetic polymers in Mediterranean surface waters, Scientific Reports, 6, 37 551, 2016.Thompson, R. C., Olsen, Y., Mitchell, R. P., Davis, A., Rowland, S. J., John, A. W., McGonigle, D., Russell, A. E., et al.: Lost at sea: whereis all the plastic?, Science(Washington), 304, 838, 2004.Turner, A. and Holmes, L.: Occurrence, distribution and characteristics of beached plastic production pellets on the island of Malta (centralMediterranean), Marine Pollution Bulletin, 62, 377–381, 2011.van Cauwenberghe, L., Vanreusel, A., Mees, J., and Janssen, C. R.: Microplastic pollution in deep-sea sediments, Environmental Pollution,182, 495 – 499, https://doi.org/https://doi.org/10.1016/j.envpol.2013.08.013, 2013.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, Environmental Research Letters, 10, 124 006, 2015.van Sebille, E., Aliani, S., Law, K. L., Maximenko, N., Alsina, J. M., Bagaev, A., Bergmann, M., Chapron, B., Chubarenko, I., Cózar, A.,et al.: The physical oceanography of the transport of floating marine debris, Environmental Research Letters, 15, 023 003, 2020.Vidal, A. S., Pasqual, C., Kerhervé, P., Heussner, S., Calafat, A., Palanques, A., Durrieu de Madron, X., Canals, M., and Puigc, P.: Acrossmargin export of organic matter by cascading events traced by stable isotopes, northwestern Mediterranean Sea, Limnology and Oceanog-raphy, 54, 1488–1500, 2009.Voth, G. A. and Soldati, A.: Anisotropic Particles in Turbulence, Annual Review of Fluid Mechanics, 49, 249–276,https://doi.org/10.1146/annurev-fluid-010816-060135, 2017.Wichmann, D., Delandmeter, P., and van Sebille, E.: Influence of near-surface currents on the global dispersal of marine microplastic, Journalof Geophysical Research: Oceans, 124, 6086–6096, 2019.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 Society Open Science, 1, 140 317,https://doi.org/10.1098/rsos.140317, 2014.Wright, S., Ulke, J., Font, A., Chan, K., and Kelly, F.: Atmospheric microplastic deposition in an urban environment and an evaluation oftransport, Environment International, 136, 105 411, https://doi.org/https://doi.org/10.1016/j.envint.2019.105411, 2020.Ye, S. and Andrady, A. L.: Fouling of floating plastic debris under Biscayne Bay exposure conditions, Marine Pollution Bulletin, 22, 608 –613, https://doi.org/https://doi.org/10.1016/0025-326X(91)90249-R, 1991.ozza, A., Drotos, G., Hernández-García, E., and López, C.: Accumulated densities of sedimenting particles in turbulent flows, Physics ofFluids, 32, http://dx.doi.org/10.1063/5.0003614, 2020.Stabholz, M., Durrieu de Madron, X., Canals, M., Khripounoff, A., Taupier-Letage, I., Testor, P., Heussner, S., Kerhervé, P., Delsaut, N.,Houpert, L., et al.: 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.Stolle, C., Nagel, K., Labrenz, M., and Jürgens, K.: Succession of the sea-surface microlayer in the coastal Baltic Sea under natural and ex-perimentally induced low-wind conditions, Biogeosciences, 7, 2975–2988, https://doi.org/https://doi.org/10.5194/bg-7-2975-2010, 2010.Suaria, G., Avio, C. G., Mineo, A., Lattin, G. L., Magaldi, M. G., Belmonte, G., Moore, C. J., Regoli, F., and Aliani, S.: The MediterraneanPlastic Soup: synthetic polymers in Mediterranean surface waters, Scientific Reports, 6, 37 551, 2016.Thompson, R. C., Olsen, Y., Mitchell, R. P., Davis, A., Rowland, S. J., John, A. W., McGonigle, D., Russell, A. E., et al.: Lost at sea: whereis all the plastic?, Science(Washington), 304, 838, 2004.Turner, A. and Holmes, L.: Occurrence, distribution and characteristics of beached plastic production pellets on the island of Malta (centralMediterranean), Marine Pollution Bulletin, 62, 377–381, 2011.van Cauwenberghe, L., Vanreusel, A., Mees, J., and Janssen, C. R.: Microplastic pollution in deep-sea sediments, Environmental Pollution,182, 495 – 499, https://doi.org/https://doi.org/10.1016/j.envpol.2013.08.013, 2013.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, Environmental Research Letters, 10, 124 006, 2015.van Sebille, E., Aliani, S., Law, K. L., Maximenko, N., Alsina, J. M., Bagaev, A., Bergmann, M., Chapron, B., Chubarenko, I., Cózar, A.,et al.: The physical oceanography of the transport of floating marine debris, Environmental Research Letters, 15, 023 003, 2020.Vidal, A. S., Pasqual, C., Kerhervé, P., Heussner, S., Calafat, A., Palanques, A., Durrieu de Madron, X., Canals, M., and Puigc, P.: Acrossmargin export of organic matter by cascading events traced by stable isotopes, northwestern Mediterranean Sea, Limnology and Oceanog-raphy, 54, 1488–1500, 2009.Voth, G. A. and Soldati, A.: Anisotropic Particles in Turbulence, Annual Review of Fluid Mechanics, 49, 249–276,https://doi.org/10.1146/annurev-fluid-010816-060135, 2017.Wichmann, D., Delandmeter, P., and van Sebille, E.: Influence of near-surface currents on the global dispersal of marine microplastic, Journalof Geophysical Research: Oceans, 124, 6086–6096, 2019.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 Society Open Science, 1, 140 317,https://doi.org/10.1098/rsos.140317, 2014.Wright, S., Ulke, J., Font, A., Chan, K., and Kelly, F.: Atmospheric microplastic deposition in an urban environment and an evaluation oftransport, Environment International, 136, 105 411, https://doi.org/https://doi.org/10.1016/j.envint.2019.105411, 2020.Ye, S. and Andrady, A. L.: Fouling of floating plastic debris under Biscayne Bay exposure conditions, Marine Pollution Bulletin, 22, 608 –613, https://doi.org/https://doi.org/10.1016/0025-326X(91)90249-R, 1991.