2.1 Lagrangian Model Set‐Up

The three‐dimensional transport of microplastic is simulated using the Lagrangian particle‐tracking framework, OceanParcels Version 2.1.6 (Delandmeter & van Sebille, 2019). We use the NEMO‐MEDUSA‐2.0 ORCA0083‐N06 output (Yool et al., 2013), hereafter MEDUSA, for the hydrodynamic data and biogeochemical data, which has a 1/12° global horizontal resolution and 75 depth levels. The data include the three‐dimensional ocean velocity fields, temperature, salinity, phytoplankton concentrations and primary productivity (the latter two are converted to algal concentrations and biofilm growth, explained in Section 2.2). We use five‐day averages from 2000 to 2009 (available from http://opendap4gws.jasmin.ac.uk/thredds/nemo/root/catalog.html). It is therefore important to note that we are investigating the effects of large‐scale advective features (such as upwelling, downwelling, convergence, and divergence). Previous work has shown that Lagrangian experiments using temporal means of up to nine days do not exhibit any significant degradation in these large‐scale flow characteristics (Qin et al., 2014).

Particles are released on a global 2° × 2° grid (9,620 particles) from 70°S to 80°N. The full latitudinal range is not used due to interpolation challenges when using the tripolar MEDUSA grid. The particles are released at 0.6 m (the surface depth in MEDUSA) and we simulate the particle trajectories for 90 days forward in time. The Lagrangian framework provides spatial and temporal interpolation of the fields following the C‐grid interpolation scheme (Delandmeter & van Sebille, 2019). We use the three‐dimensional fourth‐order Runge‐Kutta method with an integration time step of 30 s and the three‐dimensional position of each particle is stored every 12 h (180 time steps per simulation). Five sizes are chosen for the microplastic with a radius from 1000 to 0.1 μm (decreasing by one order of magnitude). Three different densities representing different initially buoyant plastic types are used; expanded polystyrene (EPS = 30 kgm⁻³), polypropylene (PP = 840 kgm⁻³) and low‐density polyethylene (LDPE = 920 kgm⁻³ ). We assume that all virtual microplastic particles are spherical and initially pristine, as in the Kooi model. It should be noted that our results are subject to the parametrizations and design of the model and that for example, the assumption of microplastic being present at the sea surface as clean and spherical must be taken into consideration when interpreting the results. To optimize efficiency (and reduce intensive computing time), we choose to run simulations for the “most typical year” from the 10 years mentioned above; 2004 (see Section S1 for a description on how we determine that year).

Two regional simulations are additionally generated in order to evaluate the effects of removing the biofouling process and removing advection on the sinking of the particles, where 25 particles are released on a 10° x 10° grid (with a 2° x 2° horizontal resolution). One region is in the algal‐rich zone of the equatorial Pacific (EqPac), from 4°S‐4°N and 140°W–148°W. The other region is in the algal‐scarce North Pacific subtropical gyre (NPSG), from 28°N‐36°N and 135°W–143°W. This latter region has been chosen since it is shown to have the highest concentrations of marine debris in the NPSG from Manta trawl samples (Egger, Nijhof, et al., 2020). We also choose to run this second simulation set‐up in the NPSG for 1000 days, from January 1, 2004, in order to test whether a longer simulation time allows for virtual microplastic particles to leave the accumulation zone. Lastly, these 1000 days simulations are repeated with horizontal advection only and vertical advection only, to isolate the effects of the different components of the three‐dimensional velocity fields.

2.2 Kooi Model

The equations governing the biofouling and sinking of microplastic in this study are based on the Kooi model (see Kooi et al. (2017) for the detailed method). Our code to generate the results can be found here: https://doi.org/10.5281/zenodo.4543145. The sinking velocity, V _s (ms⁻¹), is dependent on the density difference between the combined plastic plus biofilm and its surrounding seawater, as well as the buoyancy of the virtual microplastic (i.e., the particle’s size and density): V s ( x , y , z , t ) = − ρ t o t − ρ s w ρ s w g ω * υ s w 1 / 3

where ρ _tot is the total density of the particle plus attached algae (kgm⁻³), ρ _sw is the ambient seawater density (kgm⁻³), which is a function of the three spatial directions and time (x, y, z and t), g is the gravitational acceleration (ms⁻²), ω _* is the dimensionless settling velocity and υ _sw is the kinematic viscosity of the seawater (m²s⁻¹). The total density, ρ _tot (Equation S11), is calculated from the radius (hereafter referred to as ‘size’) of the particle (m), the volume of an algal cell (López‐Sandoval et al., 2014), the biofilm thickness (m), and the biofilm and plastic density (kgm⁻³ ) (Equation S6–S10). The seawater density, ρ _sw , is calculated as a function of conservative temperature and absolute salinity (see Roquet et al. (2015) and McDougall et al. (2003) for in‐depth explanations of TEOS‐10 standard equation of state). As mentioned above, we use temperature and salinity profiles from the MEDUSA model output that vary in time and three‐dimensional space. Using a more realistic representation of the global hydrodynamic estimates allows us to further develop the idealized one dimensional Kooi model that have profiles that only vary with depth. The dimensionless settling velocity, ω _* , is a function of the dimensionless particle diameter (Equation S18–S21) and we keep υ _sw profiles as defined in the Kooi model (Equation S1–S5; following Sharqawy et al. (2010)). Variations in kinematic viscosity are shown to be negligible for temperatures found in the ocean (Chen et al. (1973); see Section S2).

The attached algal growth, dA/dt (no. m⁻²s⁻¹), is dynamically estimated as: d A d t = A A β A θ p l + μ A A − m A A − Q 10 ( T − 20 ) / 10 R 20 A

where A _A is the ambient algal concentration (no. m⁻³), β _A is the encounter kernel rate (m³s⁻¹), θ _pl is the surface area of the spherical plastic (m²), μ _A is the growth of the attached algae (s⁻¹), m _A is the grazing or mortality rate (s⁻¹) and the final term accounts for respiration. The ambient algal concentration, A _A, is estimated by converting the MEDUSA output of diatom plus non‐diatom phytoplankton concentrations (mmol Nm⁻³) into number of algal cells per unit volume (no. m⁻³). This is calculated by first multiplying the total phytoplankton concentrations by the atomic weight of 1 mol of nitrogen (i.e., 14.007 g, to obtain mg Nm⁻³). Then, for the nitrogen to algal cell conversion, the median value in the literature (356.04 × 10⁹ ) is used, reported by Menden‐Deuer and Lessard (2000). The median value is also used in the Kooi model, however they use a carbon to cell conversion method since they estimate the theoretical depth‐profile of chlorophyll‐a concentration (and then compute the chl‐a/carbon and carbon/algal cell conversions). The encounter rate, β _A, represents the collision of the particle with water, and if the water contains algae (from A _A), we assume that the algae attaches to the particle. This β _A term is calculated as the sum of collision frequencies from Brownian motion, differential settling and advective shear (Equation S15–S17; which are dependent on diffusivity equations from Jackson (1990)). MEDUSA’s total primary productivity (mmol Nm⁻³ d⁻¹) is converted to attached algal growth, μ _A, by doing the same conversions as above for A _A (multiplying by 14.007 and dividing by 356.04 × 10⁹), however, the final step includes dividing by this algal cell concentration in order to isolate the growth factor (in d⁻¹, and then converting to s⁻¹ ). While in Kooi et al. (2017) only temperature and light limit algal growth, the MEDUSA output also considers nutrient limitation by nitrate, silicon, and iron (Yool et al., 2013). Another difference between the Kooi model and this study is that the hourly Kooi growth rate reflects a daily cycle of light intensity, while MEDUSA output is only available as five‐day averages, masking the daily cycle. However, we expect the impact on the average daily effects to be negligible for sinking timescales longer than 1 day. Since the Kooi et al. (2017) results are based on a light/dark ratio that is equal (12 h each), our average global MEDUSA‐based sinking timescales should be similar to their study (which we show is the case below). The mortality and respiration terms have been modeled identically to the Kooi model using constant rates of m _A = 0.39 d⁻¹ and R ₂₀  A = 0.1 d⁻¹ with a Q ₁₀ coefficient = 2, which represents how much the respiration rate increases by every 10°C increase in temperature; where T is temperature (°C) from MEDUSA. Since the aim of our study is to incrementally increase the complexity of the Kooi model with three‐dimensional data available from MEDUSA, processes such as grazing, mortality, and respiration are still parametrized following the Kooi model.

It should be mentioned at this stage that the Kooi equations produce unrealistic results for the sinking velocities of 10 mm size particles (which reach a V _s value of 20 ms⁻¹). Hence, we only focus on particles of a 1 mm particle radius (that have a maximum V _s of 0.2 ms⁻¹ ) and smaller. Furthermore, within the context of the Kooi model, any effects of turbulence are not included in this study. This has been decided due to the fact that a wind‐driven mixing term proposed by Kukulka et al. (2012) used in previous studies (e.g., Wichmann et al. (2019)) does not hold for negatively buoyant particles, which is the case for biofouled particles that sink in our study. Parametrizations of mixing (such as wind‐driven mixing or sub‐grid processes) must therefore be further developed in future work and have been left out of the current study. Furthermore, the advantage of this set‐up for our process study is to separate the effects of advection from wind‐driven mixing. We would expect that if one were to include vertical mixing to our study, particles with a smaller sinking velocity than that of the mixing term would exhibit shorter sinking timescales (for example, the smallest particles) and could sink deeper.

Our results are represented using two quantities. The first is the sinking timescale, T _s, defined as the number of days until the vertical velocity (V _s from Equation 1) plus the vertical advection from MEDUSA is downward, within the simulation time of 90 days. Since the particles are all initially positively buoyant, we force any rising particles to remain at 0.6 m (the ‘surface depth’ in MEDUSA), with a vertical velocity of 0 m s⁻¹ . This is to avoid algal attachment due to oscillations between 0 and 0.6 m (see Equation S17) that are meaningless since we treat 0.6 m as our surface. Simulations for the four seasons are first run separately, between December 2003 and November 2004 inclusive (released on December 1st, 2003 for DJF, December‐January‐February; and so on for MAM, March‐April‐May, JJA, June‐July‐August and SON, September‐October‐November). The four T _s maps produced have been averaged at each initial release location (for only those simulations where particles do sink within the 90 days) in order to produce one global map of T _s estimations. These T _s maps are therefore plotted where the particles are first released (as opposed to the location where the particles first sink, since it would be harder to interpret particles that overlap spatially). The first sinking depth, Z _s, for each individual season is defined as the depth each particle reaches until its vertical velocity is no longer downwards. We do this to highlight the initial effects of biofouling on the sinking of particles that are originally found at the surface, rather than the final position of the particles. If a submerged particle remains at a certain depth, or ascends, that part of the trajectory is not included in this study. Furthermore, if a particle continuously sinks throughout the 90 days and its vertical velocity remains downwards, its depth at 90 days is recorded as Z _s. We therefore capture each particle’s trajectory as a function of depth and time until Z _s, its first sinking depth.

3.1 Sinking Timescale

The globally mapped sinking timescale (T _s ) as a result of biofouling and three‐dimensional advection is estimated for different virtual microplastic sizes and densities (Figure 1). A particle sinks as soon as its density exceeds that of its surrounding seawater density and vertical advection is downwards. As explained in Kooi et al. (2017), the sinking timescale is a trade‐off between the particle’s size and the surface‐to‐volume ratio. The larger particles are more likely to encounter ambient algal cells, hence the collision rate is higher and the density can increase sooner. On the other hand, the smaller the object, the greater the relative surface area, meaning that the attachment of a very small number of ambient algae can cause it to sink. This concept can be visualized by the results for the smallest particles simulated here (0.1 μm in Figures 1m, 1n, and 1o) where particles start sinking almost immediately globally (median global T _s = 1 day for all three densities; 30, 840 and 920 kgm^‐3 ). The second‐smallest particles (1 μm in Figures 1j, 1k, and 1l) take slightly longer in some regions than the smallest particles, and then the three largest sizes (1, 0.1 and 0.01 mm in Figures 1a–1i) produce almost identical global T _s maps. For these three largest sized particles, the median global T _s is 40–43 days for the density particle type, LDPE (Figures 1c, 1f and 1i), 35–38 days for PP (Figures 1b, 1e and 1h) and 33–36.5 days for the least dense, EPS (Figures 1a, 1d and 1g; although these differences are not statistically significant, with a standard deviation of ±18 days for all nine global maps). Our results mirror the one dimensional Kooi et al. (2017) results in that an asymptotically shaped relationship between the log of increasing particle sizes versus the sinking timescale is reached in almost all locations.

Another result to highlight is what occurs in the five subtropical gyres. For the larger sizes in this study (1 mm–10 μm in Figures 1a–1i) the subtropical gyres are distinguished from the other regions by the fact that the particles do not sink within the simulation time of 90 days (black patches). This behavior can be attributed to the subtropical gyres being largely oligotrophic zones (Morel et al., 2010), implying that there is low biological activity (lighter green patches in Figure 3, left column). This leads to low algal collision frequencies such that the density of the particles does not exceed that of its surrounding seawater. On the other hand, since less algal attachment is required to increase the density of smaller particles (1 μm in Figure 1j–1l and 0.1 μm in Figure 1m–1o), the small particles in the subtropical gyres do sink within the 90 days (albeit later than in the other oceanic regions for 1 μm particles). Even when extending this to a regional simulation of 1000 days, only 2 of the 25 particles sink for LDPE particles of 1 mm (Figure S3). This suggests that as a result of the biofouling process, larger microplastic particles can be more present at the surface than smaller particles in subtropical gyres. This follows the size‐selective removal theory proposed by Cózar et al. (2014) and further supported by Egger, Nijhof, et al. (2020). Furthermore, for all particle sizes, the equatorial regions show the shortest sinking timescales of <10 days, which coincides with regions of higher algal concentrations (darker green patches in Figure 3, left column; regardless of the season). These subtropical and equatorial results prompt us to run regional simulations to better understand T _s in regions of extreme high and low algal concentrations and with different large‐scale physical processes, such as upwelling and downwelling (see Section 3.4).

Lastly, the qualitative differences in T _s as a function of size are more distinguishable than the differences in T _s for different densities of microplastic. The most sensitive particle size to density differences is 1 μm; in the subtropical gyres the 1 μm particles of 30 and 840 kgm⁻³ take longer to sink (up to twice as long, if they sink at all within 90 days), relative to 920 kgm⁻³ particles (Figures 1j and 1k relative to 1l). The general mapped patterns of T _s for the three densities are otherwise very similar globally for the other size classes and hence our results demonstrate that the global modeled distribution of sinking timescales of microplastic is very similar for different types of initially buoyant plastic (including EPS, PP, LDPE). This is somewhat expected due to the fact that the sizes can range over many orders of magnitude (five in this case) and the densities do not even reach two orders of magnitude (due to the definition of buoyant plastic types being below the average density of seawater: 1025 kgm⁻³).

3.2 First Sinking Depth

For the Z _s results, we use one of the seasons as an example; boreal spring (March‐April‐May; MAM). The deeper Z _s (m) and sooner Z _s timescale (days) for the particles of all sizes released in the northern hemisphere (red trajectories in Figure 2) can be explained by the deeper mixed layer depth and larger algal blooms during the spring in the northern Atlantic. To clarify, the Z _s timescale includes both T _s (days) as well as the time it takes from leaving the surface to reaching Z _s. The effect of seasonality on both T _s and Z _s is explored further in Section 3.3.

Although as stated in the previous section, T _s is almost identical globally for the three largest sizes, the time it takes to reach the first sinking depth, Z _s , is not identical for those three sizes (1–0.01 mm in Figures 2a–2i). Rather, it gradually takes longer to reach Z _s with a decrease in particle size. For example, the largest‐sized virtual LDPE particles can sink from the surface to around 100 m within a few days (Figure 2c), whereas the smallest‐sized virtual particles can take 40 days to reach a similar depth (Figure 2o), and the latter start sinking almost immediately around the globe. Once again, as with the T _s results, the difference in Z _s for the different initial densities is very small, therefore all further simulations have used the density of LDPE; 920 kgm⁻³.

3.3 Effects of Biological and Physical Seasonality on Sinking Characteristics

Sinking characteristics of particles can be influenced by biological and physical features, such as ambient algal concentrations and the mixed layer depth (Figure 3 and Figure 4). Since we release the virtual particles at 0.6 m, we show surface algal concentrations to illustrate the surrounding algae available for initial attachment and growth that influences when particles start sinking (i.e., T _s). The MLD is shown since waters maintain generally homogeneous properties from the surface to the base of the mixed layer. With a deeper mixed layer, for example, a deeper constant seawater density profile from the sea surface could result in sinking particles having a deeper Z _s since sudden changes in surrounding density are not encountered (Figure 4).

In this study, we distinguish the subtropical gyres as the regions that maintain the lowest surface algal concentrations across all seasons (which are consistent with previous work defining the convergence zones; e.g., Figure 4d in Maximenko et al. (2012)). Using one initial particle size (1 μm) to isolate effects of seasonality on the sinking timescale (Figures 3b, 3d, 3f, and 3h), the regions with higher algal concentrations generally correspond to shorter T _s (equatorial regions, Southern Ocean and northern Atlantic and Pacific), while lower concentrations correspond to a longer T _s (subtropical gyres) across the seasons (this is also seen for other sizes, e.g. 0.1 mm in Figure S1). Focusing on one season, the spring is known to influence algal concentrations due to algal blooms in MAM in the northern hemisphere and SON in the southern hemisphere (Figures 3c and 3g, respectively). During boreal spring, the North Atlantic and North Pacific show the lowest T _s across all seasons, with some particles sinking within 10 days in the subtropical gyres, while the southern hemisphere subtropical gyres show T _s values of up to 90 days or more (Figure 3d). Similarly, during austral spring (SON) the southern hemisphere subtropical gyres show comparatively lower T _s values relative to the other seasons, reducing the patches with T _s  > 90 days (Figure 3h).

Regarding the Z _s results as a function of seasonality (Figure 4), the spring months show consistency with T _s results. During MAM, 1 μm particles released in the northern hemisphere reach their Z _s sooner (Figure 4d), and during SON, particles released in the southern hemisphere reach their Z _s sooner (Figure 4h). This is also seen for all other sizes (for example, 0.1 mm particles in Figure S2). The maximum Z _s reached for particles released in the northern and southern hemisphere is also deeper for their respective spring months relative to the winter or summer months. This once again suggests that during spring, particles sink deeper with higher surface algae present (Figure 3c and 3g, respectively) and with deeper MLDs (Figure 4c and 4g, respectively).

Although T _s differences have been mentioned for different seasons, it is important to highlight that the results across the seasons are generally fairly similar. We can therefore further justify averaging the four seasons for each of the sizes in Figure 1. Furthermore, other physical and biological parameters from MEDUSA have also been analyzed, such as the euphotic layer depth and the total primary productivity (which is converted to algal growth for our model), however, they produce very minor global distribution differences across the four seasons (not shown), hence the focus on MLD and algal concentrations.

3.4 Regional Simulations Analyzing Advection and Biofouling

We now isolate the effects of advection and biofouling on the sinking characteristics by comparing regional results with and without these processes. Since the subtropical gyres are a recurring region of interest in this study, we choose to compare T _s in the oligotrophic North Pacific subtropical gyre (NPSG) to T _s in the algal‐rich equatorial Pacific (EqPac). Furthermore, the particle size that deserves a more in‐depth analysis is the 1 μm, since particles sink in the NPSG within 90 days (Figure 5a), contrary to particles larger than 1 μm (Figure 1).

In the NPSG simulation without biofouling while keeping advection (Figure 5b), T _s slightly increases relative to simulations with both biofouling and advection (Figure 5a). The simulation lacking advection while keeping biofouling (Figure 5c) results in almost all of the 25 particles staying at the surface throughout the 90 days. This suggests that three‐dimensional advection is the process with the largest impact on T _s for 1 μm in the NPSG. When repeating these simulations for 1000 days, and running one simulation with horizontal advection only and the other with vertical advection only, we see that vertical advection specifically plays a key role in transporting 1 μm particles below the surface in this downwelling zone (Figure S5d vs. Figure S4d). Since algal concentrations are low throughout all seasons (Figure 3) in the NPSG, the fact that removing biofouling has a lower impact on T _s is somewhat expected. The Z _s results further confirm this; those 1 μm particles that do sink within the 90 days of the MAM simulation, have a very shallow Z _s (<2 m in Figure 6a), and removing biofouling has a very minor impact on this first sinking depth (Figure 6b). If particles of all sizes are fixed to the northern edge of the NPSG (without horizontal advection), the particles can be exposed to more algae than those that are advected toward the center of the gyre, and hence can sink after > 700 days (Figure S5).

In the EqPac, although the presence of algae is consistently high throughout the year (Figure 3), T _s is almost not affected by the removal of biofouling or advection for 1 μm particles (T _s  < 10 days in Figures 5d–5f). When analyzing the impacts of these processes on Z _s (for MAM), however, if local advection in this upwelling region is isolated (and biofouling is removed), very shallow first sinking depths are produced (Z _s  < 2 m; Figure 6e) whereas Z _s can otherwise reach 20‐50 m when advection (mainly upwelling here) is removed (Figure 6f).