The Role of Stokes Drift in the Dispersal of North Atlantic Surface Marine Debris

Understanding the physical mechanisms behind the transport and accumulation of floating objects in the ocean is crucial to efficiently tackle the issue of marine pollution. The main sinks of marine plastic are the coast and the bottom sediment. This study focuses on the former, investigating the timescales of dispersal from the ocean surface and onto coastal accumulation areas through a process called “beaching.” Previous studies found that the Stokes drift can reach the same magnitude as the Eulerian current speed and that it has a long-term effect on the trajectories of floating objects. Two particle tracking models (PTMs) are carried out and then compared, one with and one without Stokes drift, named PTM-SD and PTM-REF, respectively. Eulerian velocity and Stokes drift data from global reanalysis datasets are used for particle advection. Particles in the PTM-SD model are found to beach at a yearly rate that is double the rate observed in PTM-REF. The main coastal attractors are consistent with the direction of large-scale atmospheric circulation (Westerlies and Trade Winds). After 12 years (at the end of the run), the amount of beached particles is 20% larger in PTM-SD than in PTM-REF. Long-term predictions carried out with the aid of adjacency matrices found that after 100 years all particles have beached in PTM-SD, while 8% of the all seeded particles are still floating in PTM-REF. The results confirm the need to accurately represent the Stokes drift in particle models attempting to predict the behaviour of marine debris, in order to avoid overestimation of its residence time in the ocean and effectively guide policies toward prevention and removal.


INTRODUCTION
Marine pollution is globally recognised as one of the great environmental issues our society is currently facing, with world-leading NGOs warning us that "plastic waste is flooding our oceans (Hancock, 2019)" and launching campaigns to fight "the ocean plastic crisis" (Weyler, 2017). Ocean waters have been contaminated in the last few decades by a range of man-made pollutants, from oil spills to rubber ducks (Ebbesmeyer and Ingraham, 1994). Among all these, plastic has received significant attention in the media. Most of it is nonbiodegradable and can survive in the ocean for decades, dramatically affecting marine fauna through entanglement and ingestion (Onink et al., 2019). Plastic debris comes in a huge range of shapes, sizes, and materials, making it difficult to describe and predict its behaviour in the ocean in a reliable way. To prevent inconsistencies, plastic objects are often grouped by size: micro-(< 5mm), meso-(5 − 25mm), and macro-plastics (> 25mm) (Hinata et al., 2017). Plastic is generally highly buoyant and its maximum concentration in the water column is at the surface, where it is subject to the effect of winds, currents, and waves (Reisser et al., 2015). One particular wave effect that is prevalent close to the free surface is the Stokes drift. It is defined in Oceanography as the net drift that a particle moving in a fluid experiences in the direction of propagation of the fluid's wave field (Van den Bremer and Breivik, 2017). The Stokes drift velocity, u St , is also commonly described as the difference between the average Lagrangian velocity of the particle at a mean depth, u Lag , and the Eulerian velocity of the fluid at the same mean depth, u Eu (Röhrs et al., 2012;Breivik et al., 2016): (1) The Stokes current field generally has larger spatial scales than the Eulerian velocity field and it has been shown to affect the transport of floating objects such as oil droplets (Drivdal et al., 2014), pelagic eggs and larvae (Röhrs et al., 2014), and plastic (Onink et al., 2019). The marine debris problem has been defined as a source, pathway and sink issue (Hardesty et al., 2017). Two major physical sinks of plastic pollution in the ocean are beaching at the coast and sinking to the bottom sediment (Kaandorp et al., 2020). This study aims at investigating the timescales of the beaching process of floating particles and the main accumulation zones in the presence of Stokes drift, building directly onto existing work (see e.g., Iwasaki et al., 2017;Onink et al., 2019). The hypothesis is that particles at the ocean surface will disperse and gather at the coast at a faster rate if Stokes drift is included in the model. Coastal attractors are also expected to arise at different locations as Stokes drift has been found to change the mean direction of transport, e.g., steering particles toward the poles (Onink et al., 2019). The hypothesis is tested by comparing two Lagrangian particle tracking models (PTMs): one where particles are advected by the Eulerian flow alone and one where Stokes drift is included. Reanalysis data from the HYCOM ocean model (Eulerian current) and the ERA-interim wave model (Stokes drift) are used to advect particles in the North Atlantic. The PTMs are two-dimensional (no vertical axis). To justify this choice, a theoretical background on the vertical scale of Stokes drift and particle concentration is given in section 2.1. The model runs and the methods used for data processing are described in sections 2.2 and 2.3. Section 3.1 presents an overview of the particle simulation output. The beaching process is analysed in terms of its timescales (section 3.2) and spatial patterns (section 3.3). The mean dispersal distance is calculated (section 3.4) and further forecasting of particle distribution over a longer time period is carried out (section 3.5). Finally, an interpretation of results and their implications on marine pollution is given in section 4.

Particle Concentration
The present experiment focuses on the upper meter at the ocean surface, where the concentration of buoyant particles C is at its maximum. For a steady state, C has a vertical distribution (Drivdal et al., 2014): where C 0 is the surface concentration, ν is a constant eddy diffusivity and w r is the particles' rising speed, i.e. the speed at which particles rise due to buoyancy forces. From Equation (2), the vertical decay scale of particle concentration δ p is inversely proportional to their rising speed: A typical value for eddy diffusivity in the upper ocean is ν = 10 −2 m 2 s −1 (Talley, 2011), while w r varies significantly depending on the chemical and physical structure of the considered item. For instance, oil droplets have a typical rising speed of 400 m day −1 (Drivdal et al., 2014), yielding a depth scale of about 2 m, while Arctic cod eggs with w r = 100 m day −1 would lead to δ p = 8.5 m. Marine litter items, such as small wood pieces or plastic fragments with length scales of 1−10 cm, have an estimated upward terminal velocity of 0.1 to 1 ms −1 (Hinata et al., 2017). The particle depth scale δ p would then range between 10 and 1 cm. The Stokes drift on a rotating plane is counteracted by a Eulerian return flow which dominates at depth, while the Stokes current is strongest at the surface (Van den Bremer and Breivik, 2017). For a single wave frequency, the Stokes drift decreases exponentially with depth. Its vertical scale for a full wave field, δ St , has been estimated to be around 1.8 m for basins larger than the gulf of Mexico (Clarke and Van Gorder, 2018). This depth is at least one order of magnitude smaller than the estimated Ekman depth, δ Ek , in a large basin. Given the magnitude of the three depth scales, the following approximation for the ocean surface can be made: The Lagrangian velocity experienced by particles under these conditions is equal to the sum of the Ekman and Stokes currents (Clarke and Van Gorder, 2018) (see Appendix A for a more in depth derivation). PTM-REF then corresponds to the case where δ p > δ St and, conversely, PTM-SD corresponds to δ p ≤ δ St .

The Particle Tracking Model
To investigate how the inclusion of Stokes drift affects particle transport in a Lagrangian simulation, numerical experiments are carried out with and without Stokes drift. For simplicity, the two experiments will be hereon referred to as PTM-REF and PTM-SD according to Table 1. Daily means of eastward and northward water velocity from ocean model HYCOM are used to advect particles in PTM-REF, called u Eu and v Eu , respectively. The HYCOM model (for HYbrid Coordinate Ocean Model) is a data-assimilative hybrid isopycnal-sigma-pressure (generalised) coordinate ocean model (Chassignet et al., 2007). In PTM-SD, daily means of eastward and northward Stokes drift velocities, u St and v St , from wave model ERA-Interim (Dee et al., 2011) are added to the HYCOM velocity components. Both velocity datasets range from January 1997 to December 2012 and cover the global ocean from 80 • S to 80 • N. The domain is reduced to the North Atlantic ocean between 100 • W and 60 • E and from the Equator up to the Arctic at 80 • N. Conclusions regarding the North Atlantic can be applied to the other ocean basins, where an equivalent structure in the accumulation of marine debris is observed (Maximenko et al., 2012). The initial condition consists of 2·10 4 particles distributed randomly over the spatial domain. This initial condition is repeated five times, by releasing a batch of 2 · 10 4 particles on January 1st of years 1997-2001, achieving a total of 10 5 particles. The five batches are then synced so that all particles begin drifting at the same time (t = 0) and a potential bias due to interannual variability is avoided. Seeding particles for the first 5 years is a trade-off between releasing particles for as long as possible while obtaining a drift period of at least 10 years for all particles. The OceanParcels v2.0 framework with its 4th order Runge-Kutta advection scheme is used to advect particles (Delandmeter and Van Sebille, 2019). The particle simulations produce daily outputs which are then processed in the following ways. Firstly, the North Atlantic domain is gridded in 2 • × 2 • boxes. Southern and northern boundaries are posed at 80 • N and at the Equator, respectively. Finally, particles whose velocity drops to 10 −7 ms −1 or less are considered beached and removed. This is based on the numerical assumption that such particle velocity values would only appear on land cells. The time-averaged Eulerian and Stokes drift velocities are shown in Figure 1 to gauge the main features of the two. The maximum time-averaged Eulerian flow is one order of magnitude faster than the time-averaged Stokes velocity, with peaks along the western boundaries given by the North Equatorial current and the Gulf stream, as well as the Caribbean and Guiana currents ( Figure 1A). The North Atlantic subtropical and subpolar gyres are easily identifiable: the former revolving anticyclonically between 10 and 50 • N, and the latter rotating cyclonically in the Northern sector of the domain. The direction of the Stokes drift velocity ( Figure 1B) around the subtropical gyre is consistent with large-scale atmospheric circulation (Trade Winds and Westerlies). The effect of the Polar Easterlies is also visible between 60 and 80 • N along the coast of Greenland.

Connectivity and Dispersal Distance
Adjacency matrices are a powerful tool to extract statistics on the connectivity of different locations within the domain of a given problem (Jonsson et al., 2020). Here, weighted adjacency matrices are calculated for each of the five batches of particles. Each entry p ij of such matrices represents the probability of a particle seeded in a bin i to be found in another bin j after a given time interval t. This is calculated as: where n i→j is the number of particles seeded in bin i and found in bin j at time t 1 = t and n i is the number of particles in bin i at time t 0 = 0. Only active bins are considered, i.e., bins where at least one particle was seeded at time t 0 . Each active bin is identified by a unique number ranging from 0 to m−1, where m is the total number of active bins. The resulting five matrices are then averaged to produce the stochastic matrix P ij , also known as a dispersal matrix (Botsford et al., 2009). Seasonal variability is thus captured in the dispersal matrix by mixing statistics from different years. This computation is carried out for both PTM-REF and PTM-SD with t = 1 year, to compute mean dispersal distance (or "MDD"), and with t = 10 years to make longterm predictions of particle dispersal. The MDD is calculated as (Jonsson et al., 2020): where Q ij is a matrix of geographical distance between each bin i and all other bins j. The MDD computed as in Equation (6) indicates how far from their initial position particles are likely to be found after the time interval t. Long-term predictions of particle distribution are also computed using matrix multiplication (Botsford et al., 2009). In particular, the same initial condition used for PTM-REF is reshaped into a vector of size m, called here D (0) j . This is the distribution of particles at time t 0 . Note that m is the number of 2 • × 2 • bins that make up the domain, as in section 2.2. As powers of a stochastic matrix are also stochastic matrices, it is possible to take (P ij ) n to obtain predictions for the time interval n t. The particle distribution at time t 1 = n t, called D (n) i , is then given by: Predictions obtained with the dispersal matrix are validated by comparison with the model output. Figure 2 shows the resulting anomaly. The matrix predictions overestimate concentrations in the subtropical accumulation zone compared to the model, especially in PTM-REF (dark red regions in Figures 2C-F). The maximum discrepancy for PTM-REF in the subtropical region reaches 0.4% of particles per grid-cell, 0.2% more than in PTM-SD. Underestimation occurs, on the other hand, in almost all coastal bins, indicated by the blue lining along the land boundary.
In particular, the matrix prediction inaccurately represents the concentration of particles along the coast of Central and South America, especially in PTM-SD, as well as the western coast of Europe and the southern coast of west Africa (Figures 2E,F).
The magnitude of the discrepancies increases over time as higher powers of the matrix are taken (low anomaly in Figures 2A,B). Nonetheless, the matrix and model outputs are overall in good

An Overview of the PTMs
Snapshots of the two particle tracking model runs are displayed in Figure 3 to gain a qualitative view of the particle's behaviour over time and highlight the differences that arise when Stokes drift is included. The results suggest that there are areas of particle convergence and areas of particle divergence. After 6 months ( Figures 3A,B), all particles seeded North of the Equator have drifted away from this region. A similar divergence is also apparent at the Arctic boundary east of Greenland, as well as in the Gulf of Mexico and Caribbean Sea. After 2 years of run time (Figures 3C,D), the formation of a subtropical accumulation zone around 30 • N becomes visible, as observed in previous studies (Maximenko et al., 2012). The area between this subtropical patch and the Equator has been entirely cleared of particles. A similar process appears after 6 years ( Figures 3E,F run ( Figures 3G,H), all remaining particles (i.e., particles that have not beached) are in the accumulation zone. Overall, particle convergence/divergence occurs in the same areas in PTM-REF and PTM-SD. However, the timescales of dispersal from the subtropical accumulation zone differ. In particular, concentrations decrease over time at a faster rate in PTM-SD in this region. Just over 14% of particles are left in the domain by the end of the PTM-SD run, compared to almost 34% in PTM-REF.

Beaching Timescales
Not enough clarity has yet been made on the residence time of plastic pollution in the ocean, with estimations ranging from a few years to millennia before it sinks to the bottom or it is washed ashore (Morét-Ferguson et al., 2010). The present experiment provides some insight on the timescales of particle dispersal from the oceanic accumulation zone. The beaching test described in section 2.2 is applied to the PTM-REF and PTM-SD cases. The resulting amount of beached particles is shown, in percentage, in Figure 4A, while Figure 4B depicts the concentration in the subtropical accumulation zone. The coastal accumulation increases rapidly within the first 2 years and it continues growing, albeit with a smaller gradient, until the end of the run. After 12 years, 86% of particles seeded in PTM-SD have beached; 20% more than in PTM-REF. New particles are seeded annually over the entire domain for the first 5 years. Some of these particles are located, at time 0, in those areas of divergence seen in Section. 3.1. The extreme increase in the number of beached particles in the Figure 4A reflects the quick removal of particles from divergence areas. As these regions empty out, the beaching process slows down. The average annual beaching rate (excluding the 5 years of seeding) is 1.6% in PTM-SD. This rate is halved in PTM-REF (0.8%). When particles are seeded, about 22% of particles are located in the offshore accumulation zone. The concentration in this subtropical patch increases within the first 4-5 years, partly due to the constant addition of particles during this time, reaching a maximum of 37% of particles in PTM-REF and 25% in PTM-SD. After the first 3-5 years, it then decreases steadily with a stronger gradient in PTM-SD. The anomaly between the two experiments increases over time, indicating that the subtropical gyre is a stronger attractor without the influence of Stokes drift (PTM-REF case). In particular, the subtropical box contains 32% of particles in PTM-REF at the end of the run, 20% more than the amount found in PTM-SD. The oscillation in Figure 4B appears as a result of seasonal variability.

Spatial Differences
Discrepancies in the beaching process emerge not only in terms of timescales, but also as far as the main coastal areas attracting particles. It is thus worth looking at their distribution at the end of the run on each side of the domain, east and west (see Figure 5). Figures 5A,B show the meridional distribution of beached particles on each side of the domain as a sum in the longitudinal direction of the amount of particles per square degree. Figure 5C depicts the same distribution for all beached particles. Percentages are calculated from the total amount of particles. From the top panels, it is clear that PTM-REF (blue dashed line) and PTM-SD (red line) have a similar distribution in the meridional direction. In particular, they present maxima and minima at similar latitudes. No particles are found anywhere on the coast of the US between 33 and 45 • N ( Figure 5A) and on the coast of west Africa between 10 and 27 • N ( Figure 5B). Particles are located mainly North of 30 • N in the east, while they appear split in half between north and south in the west. To gain a better understanding of the north-south distribution, the cumulative distribution function (CDF) is shown in the bottom panels (Figures 5D-F). It is calculated by integrating the meridional

Mean Dispersal Distance
Mean dispersal distance (MDD) calculated as in section 2.3 is shown in Figure 6. Low MDD (close to 0 km) is found within the subtropical convergence zone at 30 • N and along the coastal Frontiers in Marine Science | www.frontiersin.org 9 August 2021 | Volume 8 | Article 697430 regions of West Africa, Northern Europe, and the Caribbean. A significant amount of particles was found to beach in these regions (see section 3.3). Particles seeded in any of these locations are likely to be in the same position after a year. The sector corresponding to the centre of the subtropical accumulation zone showing low MDD is larger in PTM-REF ( Figure 6A) than in PTM-SD (Figure 6B), suggesting that particles are more likely to leave the gyre within a year if Stokes drift is included in the PTM.
The highest values of MDD surround this area, following the anticyclonic circulation in the North Atlantic. From the anomaly plot (Figure 6C), it can be noted that MDD in the PTM-SD case exceeds PTM-REF almost everywhere, with the maximum difference found by the Canary islands (star marker). However, the mean dispersal distance in PTM-REF is larger along the strong Guiana and Caribbean currents. Stokes drift is weak here (see Figure 1) and pushes directly toward the coast, explaining why particles initialised here in PTM-SD are not likely to travel far. An area of negative anomaly is also found within the subpolar gyre.

Matrix Predictions
The calculation of an adjacency matrix allows for predictions that go beyond the 12 years of the model run and require minimal computational time. With this goal, matrix P ij is computed for PTM-REF and PTM-SD with t = 10 years. This allows to investigate particle distribution on a timescale of decades. It is worth pointing out that this calculation is valid for any initial condition. The calculation in Equation (7) is carried out using (P ij ) n with n = 2, 5, 10 using the same initial condition as PTM-REF.
The resulting distributions are shown in Figure 7. Previously made observations on the location of divergence and convergence zones still apply over a longer timescale. In particular, non-zero concentrations are observed only within the subtropical accumulation zone. The garbage patch at 30 • N is present in both cases and in all snapshots, albeit with varying concentrations. After 20 years (Figures 7A,B), 28.3% of particles are still floating at the ocean surface in PTM-REF.
In PTM-SD, on the other hand, only 8.4% of all particles remain in the ocean after 20 years: around 20% less than in the case without Stokes drift. This anomaly agrees with the findings in sections 3.2 and 3.3. Fifty years after initialisation (Figures 7C,D), the amount of active particles has decreased by over 10% in PTM-REF and by 7% in PTM-SD, reaching only 1.1%. All particles have beached after 100 years in PTM-SD, while almost 8% of the initial amount remain active in PTM-REF (Figures 7E,F).

Interpretation of Results
The motivation behind the present study was the hypothesis that the inclusion of Stokes drift in models investigating trajectories of floating objects is of particular importance, as implied by previous analyses (Drivdal et al., 2014;Onink et al., 2019). This hypothesis is hereby confirmed as Stokes drift is seen to have an impact both on the timescales of particle dispersal and on the concentration in the main convergence zones. An accumulation zone centred at 30 • N is observed within the first 2 years of runtime in both PTMs, covering the subtropical gyre. An equivalent aggregation of particles in the subtropics has been found in a variety of modelling experiments (Lebreton et al., 2012), as well as from trajectories of satellite-tracked Lagrangian drifters (Maximenko et al., 2012) and in-situ observations (Morét-Ferguson et al., 2010). Similar build-ups appear in the subtropics in all five ocean basins as a result of Ekman convergence (Onink et al., 2019). The distribution of beached particles at the end of the PTM runs reveals that the most affected coastal areas are northern Europe on the eastern side of the North Atlantic and Central/South America on the western side. This is in line with the direction of large-scale atmospheric circulation: particles are pushed southwest-ward by the Trade Winds and northeast-ward by the Westerlies. From the mean dispersal distance calculation (section 3.4), it can be inferred that particles released near the European coast are likely to beach locally. This is in agreement with previous studies where the densely populated north-east Atlantic region turned out to be highly impacted by beaching and coastal retention (Chenillat et al., 2021). On the other hand, small islands in Central America may be particularly affected by nonnative marine litter. These islands are located at the rim of the subtropical accumulation zone, from which they receive a constant leakage of particles which may be originating from the other side of the Atlantic and travelling across the ocean in less than a year (Lachmann et al., 2017). The increased mean dispersal distance in PTM-SD may occur as a result of the method chosen to include Stokes drift in the model. As the Stokes drift here is linearly added to the Eulerian flow, the velocity advecting particles in PTM-SD is always larger in magnitude than in PTM-REF, causing particles to travel faster. However, the consistency in direction of the Stokes drift and its predominantly zonal pattern suggest that such an increment in MDD is plausible. The adjacency matrix predictions compare well with the model results (Figure 2). Despite the coarse resolution of P ij , it succeeds in adequately representing the formation of the North Atlantic garbage patch and its decay over time. The long-term prediction (Figure 7) confirms that Stokes drift increases the rate of particle dispersal in the subtropical accumulation zone over time. The Stokes velocity measured by Lagrangian surface drifters has been shown to reach a comparable magnitude to the Eulerian current measured by ADCPs during strong wind events. However, under normal conditions, its magnitude decreases to only 20% of the Eulerian current on average (Röhrs et al., 2012). Despite this seemingly small contribution to the total current, the Stokes drift has a non-negligible effect on the rate of dispersal of floating particles. We find that the number of beached particles per year

Areas of Improvement for the PTM
The present simulation succeeds in capturing the effects of large-scale phenomena acting on particle transport under the chosen conditions, in particular the formation of convergence and divergence zones. However, the model does not account for (1) the sinking of particles to the ocean bottom and (2) the possibility for particles that have washed ashore to drift back into the ocean, thus overestimating the amount of beached particles. The residence time at the surface before sinking due to biofouling is estimated in Chubarenko et al. (2016) for different types of micro-plastic particles. The study finds that polyethylene fibres spend 6-8 months within the upper layers, while spherical particles and plastic pieces remain at the surface for up to 10-15 years. In the latter case, the inclusion of sinking in the present model would not significantly affect the result, as the timescale of sinking is larger than the total run time and estimated residence time. However, in the case of polyethylene fibres or-more generally-of objects whose sinking rate is comparable to the beaching rate or faster, the inclusion of a sinking flux in the model would be of great importance. In future work, a more realistic budget of particle density should be employed with the inclusion of a flux transporting particles from the surface to the ocean bottom. In the framework of adjacency matrices, this can be represented by the probability p srf →btm of transitioning from the surface to the ocean bottom during a time interval t (Liubartseva et al., 2018): where s is the rate of sedimentation, τ age is the particle's age and T surf is its residence time at the ocean surface. Note that the probability of sinking increases exponentially with the particle's age, giving a realistic representation of the weathering of marine debris at the surface and the consequent buoyancy loss. In terms of beaching, a strict boundary condition at the coast was chosen: particles reaching a land-cell and a certain velocity (< 10 −7 ms −1 ) are considered beached and erased. Beaching is thus an irreversible process in the PTMs. This is however not necessarily realistic. Items of marine litter that have washed ashore may in fact drift back into the ocean as a consequence of e.g., a change in wind direction or strong wave action at the coast. With that in mind, a less rigid boundary condition was tested in previous runs of the PTMs, where beaching was reversible. In this configuration, particles whose velocity became smaller than 10 −7 ms −1 were counted as "beached" but not erased, allowing them to drift back into the flow. The annual rate of dispersal from the subtropical accumulation zone and onto the coast remained doubled in PTM-SD compared to PTM-REF.
A seasonal pattern was discovered in PTM-SD with a higher number of beached particles in summer and a decrease in winter. The same oscillation was not observed in PTM-REF and is likely due to a change in wind direction affecting the Stokes drift. However, there is no indication that the way the model represented this flux in and out of coastal bins is realistic and was thus neglected here. As a future improvement to the PTM, a horizontal diffusivity term such as the one used in e.g., Iwasaki et al. (2017) would have to be included to obtain a more realistic image of the beaching process.

Application to Marine Plastic
Plastic debris in the ocean has two main sinks: the coast and the bottom sediment (Kaandorp et al., 2020). The buoyancy of an object highly influences its fate, determining which of the two sinks it will end up in. Buoyancy depends on chemical composition and physical properties, such as size and shape. Plastic objects at the ocean surface are subject to weathering due to UV radiation as well as temperature and salinity changes (Jahnke et al., 2017). Weathering leads to fragmentation, which is the process of larger objects breaking up into smaller pieces, leading to the formation of "secondary microplastic" (Jahnke et al., 2017). The rate at which microplastics sink to the bottom sediment is inversely proportional to their characteristic lengthscale: given the same buoyancy, the smaller the object, the faster it sinks (Chubarenko et al., 2016). Sinking happens mainly as a consequence of biofouling, as the growth of biofilm around plastic objects is found to increase their density (Jahnke et al., 2017). This suggests that microplastics are more likely to settle within the water column and eventually sink, while largeror in general, more buoyant-objects are more likely to be transported by winds and currents to one of the convergence zones, whether it be the coast or the subtropical garbage patch. As seen in section 2.1, the Stokes drift has its maximum at the ocean surface. The Lagrangian velocity of highly buoyant items at the ocean surface is equal to the sum of the Eulerian velocity and the Stokes drift velocity (Equation A5, Clarke and Van Gorder, 2018). If these objects are located just below the surface, where they are not directly affected by winds, their behaviour is then closely described by the PTM-SD configuration of the model. Nevertheless, the PTM-REF case is not unrealistic. As seen in section 2.1, the Stokes drift is counteracted by an Eulerian return flow which prevails below the Stokes depth δ St . The return flow is not present at the ocean surface and was thus not accounted for in PTM-SD (Equation A5). However, if the depth scale of particle concentration, δ p , were larger than δ St , the effect of the Stokes drift on particle transport would be increasingly dampened by the Eulerian return flow as depth increases. This would lead to a scenario that is closer to PTM-REF.
A similar result was derived in Iwasaki et al. (2017), where the contribution of Stokes drift to particle motion was found to depend on particle size, with buoyant large particles moving more rapidly at the surface.

CONCLUSION
This work presents a comparison between two different particle tracking models (PTMs) to shed light on the effect of Stokes drift on the dispersal of floating objects at the ocean surface.
In the first experiment (PTM-REF) particles are advected by the Eulerian current alone, while in the second (PTM-SD) the Stokes drift is added linearly to the Eulerian current to transport particles. Global renalysis products HYCOM and ERAinterim are used for the Eulerian velocity and Stokes drift, respectively. Particles are seeded with a random geographical distribution over the North Atlantic and at different times to avert biases caused by inter-annual variability. The results reveal areas of particle convergence both within the basin and along the coast. An offshore accumulation of particles develops in the subtropics, around 30 • N. This is observed in all five ocean basins and has been referred to in the literature as a "garbage patch." The PTMs and adjacency matrix predictions show that the garbage patch disperses particles faster in PTM-SD, with an annual beaching rate that is doubled compared to the one in PTM-REF. The coastal environment is highly affected as it becomes a stronger attractor for floating particles than the open ocean. The model suggests Western Europe and Central America are the coastal regions most affected by beaching of particles. However, due to enormous difficulties in sampling plastic debris in the ocean, it is difficult to validate these results against independent information. It is shown by the mean dispersal distance calculation that particles beaching in Central America, which is located at the edge of the subtropical gyre, may originate from all the way across the ocean. Particles released off the coast of West Africa in PTM-SD, where the MDD is at its maximum, take under a year to be found almost 6,000 km from their starting point. On the other hand, particles found stuck on the European coast are likely domestic. Investigating the connectivity of different regions of an oceanic basin gives insight on the distances travelled by particles and their timescales. Such considerations can be useful in informing policy change on marine pollution prevention and clean-up. It is estimated that 60-80% of marine debris is comprised of plastic (Lebreton et al., 2012). The present results suggest that the idea that marine plastic may spend hundreds of years floating in the open ocean is likely a misconception. We highlight the need to take into account the Stokes drift velocity when modelling the transport of marine debris by showing that ignoring this wave effect may be causing an overestimation of the residence time of particles in the subtropical accumulation zones. This is of particular importance when modelling the behaviour of buoyant objects of marine litter that sit just below the surface, where the Stokes drift is at its maximum. On the other hand, for debris that tends to settle lower in the water column (such as certain items of microplastic, Hinata et al., 2017), the effect of Stokes drift on particle transport is dampened by the Eulerian return flow. Under these assumptions, both PTMs are deemed valid descriptors of two different limit cases: PTM-REF applies to the case where the depth scale of particle concentration is larger than the Stokes depth (δ p > δ St ), while PTM-SD better reflects the case where particles concentrate at the very surface of the ocean (δ p ≤ δ St ).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.