Modeling submerged biofouled microplastics and their vertical trajectories

. The fate of (micro)plastic particles in the open ocean is controlled by physical and biological processes. Here, we model the effects of biofouling on the subsurface vertical distribution of spherical, virtual plastic particles with radii of 0.01–1 mm. For the physics, four vertical velocity terms are included: advection, wind-driven mixing, tidally induced mixing, and the sinking velocity of the biofouled particle. For the biology, we simulate the attachment, growth and loss of algae on particles. We track 10,000 particles for one year in three different regions with distinct biological and physical properties: 5 the low productivity region of the North Paciﬁc Subtropical Gyre, the high productivity region of the Equatorial Paciﬁc and the high mixing region of the Southern Ocean. The growth of bioﬁlm mass in the euphotic zone and loss of mass below the euphotic zone result in the oscillatory behaviour of particles, where the larger (0.1–1.0 mm) particles have much shorter average oscillation lengths ( < 10 days; 90th percentile) than the smaller (0.01–0.1 mm) particles (up to 130 days; 90th percentile). A subsurface maximum concentration occurs just below the mixed layer depth (around 30 m) in the Equatorial Paciﬁc, which is 10 most pronounced for larger particles (0.1–1.0 mm). This occurs since particles become neutrally buoyant when the processes affecting the settling velocity of the particle and the motion of the ocean are in equilibrium. Seasonal effects in the subtropical gyre result in particles sinking below the mixed layer depth only during spring blooms, but otherwise remaining within the mixed layer. The strong winds and deepest average mixed layer depth in the Southern Ocean (400 m) result in the deepest redistribution of particles ( > 5000 m). Our results show that the vertical movement of particles is mainly affected by physical 15 (wind-induced mixing) processes within the mixed layer and biological (bioﬁlm) dynamics below the mixed layer. Furthermore, positively buoyant particles with radii of 0.01–1.0 mm can sink far below the euphotic zone and mixed layer in regions with high near-surface mixing or high biological activity. This work can easily be coupled to other models to simulate open-ocean biofouling dynamics, in order to reach a better understanding of where ocean (micro)plastic ends up.


20
Observations have shown that plastic ends up everywhere in the ocean, from the Arctic (Cózar et al., 2017), to the Mariana Trench (Chiba et al., 2018;Peng et al., 2018). Two key questions that many studies address are: how much plastic is found in each ocean reservoir and how does it get there? To tackle the first question, studies using the 2010 global Jambeck et al. (2015) estimate of plastic entering the ocean from the coasts (4.8-12.7 million tons), suggest that approximately 1% is at the sea surface (Eriksen et al., 2014;Van Sebille et al., 2015) and 67-77% ends up on beaches or in coastal waters, up to 10 km 25 offshore (Lebreton and Andrady, 2019;Onink et al., 2021). Following these approximations, 20-30% of ocean plastic debris is unaccounted for and could either be in the water column or on the seafloor. The focus of this study is therefore to explore the processes affecting the vertical distribution of ocean plastic.
Logistical constraints result in limited (and sometimes only shallow and coastal) subsurface observations of plastic concentration (e.g. Reisser et al., 2015;Kanhai et al., 2017;Dai et al., 2018;Pabortsava and Lampitt, 2020;Kukulka et al., 2012). nutrient levels, diatom concentrations and primary productivity are very low in these oligotrophic waters. Lastly, the average vertical mixing profile in the SO reaches almost an order of magnitude higher than the other two regions (0.12 m 2 s −1 ) and extends to 500 m deep (Fig. 1c).
In each of the three regions, we initiate a total of 100 release locations on a 1 • x 1 • grid at z = 0.6 m (the sea surface 90 depth in MEDUSA). Within the Parcels Lagrangian framework, a C-grid interpolation scheme is used for temporal and spatial interpolation of the fields (Delandmeter and van Sebille, 2019). The three-dimensional fourth-order Runge-Kutta method is used with an integration time step of 60 seconds, and the 3-D position and biofouling state of each particle is stored every 12 hours.
We simulate spherical plastic with a radius between 0.01 and 1 mm. The upper limit is the same as in Lobelle et al. (2021), 95 to comply with the Stokes motion equations (Kooi et al., 2017), and the lower plastic radius limit is to comply with the lowest order of magnitude on which diatoms have been observed to attach (Amaral-Zettler et al., 2021a, b). We use 25 size bins within this range, releasing 4 identical particles per bin at each release location. Since we have 100 release locations in each region of study, we simulate 400 particles per bin size per region, and 10,000 particles in total per region. Table 1. The comparison between the modeled physical and biofilm dynamic specifications, progressing from the Kooi et al. (2017) model to the Lobelle et al. (2021) model and to the model in this study.

Models
Kooi et al. 2017 Lobelle et al. 2021  We also assign an initial density to the virtual particles. Following results in Lobelle et al. (2021) showing that sinking 100 characteristics of biofouled plastic particle with initial densities of 30 and 920 kg m −3 show minor differences, we focus on one density here, 920 kg m −3 (representative of low-density polyethylene; one of the most commonly produced plastic polymers). We have also run a sensitivity analysis for particles with a density of 30 kg m −3 (representative of expandable polystyrene) and found that the only difference to 920 kg m −3 results is that in the oligotrophic NPSG, particles of the largest size class remain at the surface instead of sinking to the base of the mixed layer (Fig. C1b). The particles in our simulations 105 are spherical and we assume that they do not fragment, change shape or change density throughout the simulations.

Physical dynamics
There are three components that make up the physical dynamics included in this study: vertical advection (from NEMO-MEDUSA), computed wind-induced vertical mixing and computed tidally induced vertical mixing. Any horizontal motion is 115 not included in order to keep particles in their original release coordinates, which have specific biophysical profiles we are interested in analysing.
Ocean General Circulation Models (OGCMs) generally have very low-resolution representations of turbulence. However, high-resolution wind-driven turbulence can play an important role in the vertical concentration profiles of buoyant particles (Kukulka et al., 2012) and is therefore its inclusion is one of the novelties of this study (Table 1) amount of turbulence in the surface mixed layer is computed using the K-profile parametrisation (KPP) (Large et al., 1994;Boufadel et al., 2020), K z , given by: where κ = 0.4 is the Von Kármán constant, u * w is the frictional velocity of the seawater's surface [m s −1 ], φ = 0.9 is the "stability function" of the Monin-Obukov boundary layer theory, θ = 1 is a Langmuir circulation enhancement factor (which is 1 when Langmuir cells are negligible, as is the assumption here), z 0 is the roughness scale of turbulence  (2019) is: where u 10 is the 10m wind speed [m s −1 ] and g = 9.81 m s −2 is the acceleration due to gravity. We use a local form of the KPP profile, meaning that our profiles only depend on the surface wind stress at that location and do not consider any non-local effects, such as horizontal advective transport or non-locally generated swell waves. Our boundary condition at the surface is 135 set such that if a particle is about to cross the surface, we set its depth to the NEMO-MEDUSA surface depth (0.6 m). Onink et al. (subm.) show that for positively buoyant particles, 1D vertical profiles estimated with this Markov-0 approach match reasonably well with observations, where increased wind stress results in increased mixing of particles and reduced particle rise velocities.
The vertical stochastic velocity perturbation due to turbulence, w , according to Gräwe et al. (2012) and solved using an 140 Euler-Maruyama scheme (Maruyama, 1955), is given by: where ∂ z K z = ∂K z /∂z, ∆t is the integration timestep, dW is the Wiener increment with a mean of 0 and standard deviation, σ = √ ∆t and V s is the sinking velocity of the particle defined in Eq. (5).

145
Since KPP only estimates turbulent mixing above the MLD, we also include a background full-depth vertical mixing induced by internal tides (see de Lavergne et al. (2020) for the detailed methodology). The global tidal mixing maps they provide which we use to estimate tidally formed K z are available from: https://www.seanoe.org/data/00619/73082/.

Particle settling velocity
The Lobelle et al. (2021)  two factors: 1) the density difference between the biofouled plastic particle and the surrounding seawater and 2) the size and density of the particle: where ρ tot is the total density of the particle plus attached algae [kg m −3 ], ρ sw is the ambient seawater density derived from 155 NEMO-MEDUSA's temperature and salinity [kg m −3 ] using the TEOS-10 standard equation of state (see Roquet et al., 2015;McDougall et al., 2003), ω * is the dimensionless settling velocity and υ sw is the kinematic viscosity of the seawater [m 2 s −1 ].
The settling velocity, V s [m s −1 ], can therefore be computed as a function of the three spatial directions and time (x, y, z and t). The supplementary material of Lobelle et al. (2021) describes the equations behind each term in Eq.( 5). Also, here the kinematic viscosity has been computed dynamically in 3D space and time as opposed to using the same profile as defined 160 in Kooi et al. (2017), though the spatiotemporal variations are so small that this modification has a minor impact on the results (Chen et al., 1973).

Biofilm dynamics: gain and loss terms
In this study we assume that the biofilm solely consists of diatoms. While the original model described in Lobelle et al. (2021) characterised biofilms as being comprised of diatoms and non-diatoms, Amaral-Zettler et al. (2020) suggest the use of diatom-165 dominated biofilms versus non-diatom dominated biofilms. The attached biofilm, dA dt affects the total volume and density of the particle + biofilm (Eq. (5)-(10) in Kooi et al., 2017), which determines the dimensionless settling velocity ω * in Eq. (5) above. In this study, we use two gain terms and three loss terms to define the algal biofilm dynamics: The two gain terms are identical to the Lobelle et al. (2021) model (Table 1). As such, a particle's collision with and colonisation by algae; β A is the collision rate [m 3 s−1], A A is the planktonic algal concentration [no. m −3 ] and θ pl is the surface area of the spherical plastic [m 2 ] (see Eq. (S15)-(S17) in Lobelle et al., 2021). Note that the term 'planktonic' algae is used in this study to refer to the algae present in seawater and 'attached' to refer to the algae present in the biofilm. The growth of the attached algal cells is . Since TPP3 is the total primary 175 productivity of both diatoms and non-diatoms, and the productivity of diatoms alone is not available, we assume that the rate is the same for both phytoplankton. We divide TPP3 by the total planktonic diatom + non-diatom concentration, and use that rate to multiply by the total number of attached diatoms. Therefore, after converting TPP3 to an algal growth rate by multiplying by the atomic weight of nitrogen (14.007 g) and then using the median nitrogen to algal cell ratio ( (Yool et al., 2013) to compute spatially and temporally varying grazing of diatoms in biofilms as well as adding another term, nonlinear losses, which include viral lysis (Table 1). The grazing of diatoms by mesozooplankton is available from NEMO-MEDUSA, however, only as a depth-integrated variable (in mmol N m −2 s−1).
We therefore recompute 3D depth-dependent grazing (in mmol N m −3 s−1) dynamically using Eq. (54) from Yool et al. (2013): Here, g m is the maximum zooplankton grazing rate (0.5 d −1 ), p mP d is the dimensionless mesozooplankton preference for N:algal cell conversion as described above is used. We then divide by the number of attached algal cells in Eq. (6) to get a grazing rate. This means that we assume that the mesozooplankton grazing rate is the same for the planktonic algae as for the attached algae on microplastic at a specific point in time and space.

200
The loss rate via respiration remains as in the Kooi et al. (2017) model and is therefore: where R 20 A = 0.1 d −1 with a Q 10 coefficient = 2, which represents how much the respiration rate increases by every 10 • C increase in temperature; where T is temperature ( • C) from MEDUSA.
The final loss term represents processes that depend on the abundance of diatoms, such as diseases (including viral lysis).

205
This term is represented using a saturating hyperbolic function defined in Eq. (72) of Yool et al. (2013): where λ is a nonlinear maximum loss rate of 0.1 d −1 , and k P d is the loss half-saturation constant (0.5 mmol N m −3 ). As with the grazing above, these nonlinear losses are determined relative to the abundance of diatoms in seawater. This loss rate is then applied to the number of algae attached to the particle.

210
One of the key assumptions we make in this study is that the biofilm only consists of diatoms, since this is the data available in NEMO-MEDUSA. Understanding the composition of the biofilm community (or plastisphere) is important to accurately model the effects of biofilm dynamics on the vertical motion of particles. Recent work in the Mediterranean and the North Sea's coastal waters (Amaral-Zettler et al., 2021a, b) has observed that for a spherical polyethylene particle with a radius of 30µm or smaller, small single-celled organisms (including bacteria) need to be considered. We therefore focus on simulating

The vertical distribution of particles
Our work supports the findings from previous studies that floating particles can sink in the open ocean due to biofouling. The sinking behaviour (and vertical distribution) of particles varies for different locations and particle sizes (Kaiser et al., 2017;220 Choy et al., 2019;Lobelle et al., 2021). We compare vertical one-year trajectories in three regions (Equatorial Pacific; EqPac, North Pacific Subtropical Gyre; NPSG, and Southern Ocean; SO) and for two particle radius size classes (0.01-0.1 mm and 0.1-1.0 mm; Fig. 2). We reiterate that our simulations only include vertical motion (advection and mixing) in order to isolate localised biological and physical effects on vertical particle displacement.
There are a few main results to highlight. Firstly, particles can sink far below the euphotic zone depth (EZD) and mixed Andrady (1991) had theorised that eventually particles can be fouled so heavily that they would permanently sink. We discuss potential reasons why so few of our particles reach the seafloor in Sect. 3.4. One of the key novelties of our study, introducing vertical wind-driven mixing, is to verify whether particles still oscillate as in Kooi et al. (2017) and Kreczak et al. (2021). We demonstrate that particles initiate their oscillations upon sinking below the mixed layer and in Sect. 3.3.1 we characterise these oscillations relative to particle size. In all three regions and for both size classes the average MLD seems to affect the vertical 235 distribution of particles more than the average EZD. This seems contrasting to the findings in Kreczak et al. (2021) that the EZD defines the vertical displacement of particles. Since they do not include advection or mixing in their model, this supports how important wind-driven mixing is for vertical displacement of particles <1 mm in the ocean.

High productivity region
In the EqPac, high biological activity means a biofilm can quickly form and increase the density of a particle. The most 240 distinct feature is a subsurface maximum particle concentration throughout the year which generally appears just below the average MLD (around 30 m; Fig. 2a and d). The subsurface maximum is seen more prominently for the largest particles (0.1 to 1.0 mm) than the smallest ones (0.01 to 0.1 mm) and can reach a relative annual average concentration of 3 times the surface concentrations ( Fig. 3a; up to 5 times considering the 95% percentile). Such an accumulation of particles at a certain depth can occur when the processes affecting the particles' upwards and downwards movement are in balance. In Kreczak et al. (2021),   was much deeper, at 200 m). The maximum depth reached by the smallest particles is around 1700 m, and the EqPac is the 250 clearest example of larger particles oscillating with a higher frequency than the smaller ones. In a sensitivity analysis where we allow algal cell walls to remain attached to the particle after the biofilm dies, smaller particles no longer show a subsurface maximum and instead can reach much deeper depths (to over 4800 m), since even a small change in density can affect its vertical transport and would cause it to sink deeper (Fig. D1d).

Low productivity region 255
In the oligotrophic NPSG, most of the particles of both size classes remain above the MLD throughout the year (Fig. 2b and e).
Most of the largest particles remain at the sea surface, and the rest are distributed throughout the mixed layer (Fig. 2b). This is clearly seen in the annual average distribution profile, where below 10 m the particles' concentration is almost 0 relative to the surface (Fig. 3b). The smallest particles in the NPSG are more evenly distributed from the surface to the base of the mixed layer than the larger particles ( Fig. 2e; with an average annual MLD above 50 m; Fig. 3b). These results are comparable to NPSG 260 observations showing a power law decline in microplastic concentrations with depth, where most particles are found in the upper tens of metres (Egger et al., 2020). Seasonality plays a role in the NPSG, where the boreal spring bloom (February and May) causes enough biofouling to occur for particles to sink below the MLD. A feature to highlight are a couple of horizontal lines representing subsurface maximum particle concentrations between 200 and 300 m in the larger particles (Fig. 2b). One possible explanation for this is a local equilibrium of biological and physical processes causing vertical displacement. During 265 a sensitivity test using an initial particle density of 30 kg m −3 (representing expanded polystyrene) instead of 920 kg m −3 , the NPSG is the only region where vertical distribution changes drastically (Fig. C1b). Apart from during the spring bloom, all larger 30 kg m −3 particles remain at the surface. Following results from Lobelle et al. (2021), this suggests that even under surface-mixing conditions, plastic with a radius of 0.1-1 mm and with a very low density might very rarely sink in oceanic regions with low algal concentrations. In the sensitivity analysis including the dead cell attachment as explained above, the 270 smaller particles also have much longer oscillation times below the MLD, suggesting that the slower loss of biofilm mass affects the smaller particles much more than larger particles ( Fig. D1e and b).

Intense mixing region
On average in the SO, particles mix to much deeper depths; the SO is the only region where the annual average particle depth distribution is not close to 0 at 450 m, relative to the surface (Fig 3c). The MLD varies greatly with the seasonal cycle (Fig 2c   275 and f), gradually deepening from 50 m in the austral summer (January to March) to 400 m by the austral spring (October). The maximum depth reached by the particles in the SO (around 5000 m for both particle size classes) is in October for the largest size class. For the smallest particles, the maximum depth is delayed by a month or two due to the particles' longer oscillation lengths (further explored in Sect. 3.3.1).

The influence of vertical advection and mixing 280
To visualise the processes that determine the vertical displacement of particles we display the ratio between the particles' absolute settling velocity, V s in Eq. (5), and the ambient water's vertical movement (Fig. 4). As described in the Methods section, the settling velocity is dependent on the initial size and density we assign to the particle and the biofouling dynamics, and the ambient water's movement includes wind-driven mixing, tidally induced vertical mixing and vertical advection. This means that particles present below the MLD are subject to vertical advection and vertical tidally induced background mixing only 285 (which are generally two orders of magnitude lower than wind-driven mixing in all regions; see Fig. A1-A3). Throughout all the simulations, two distinct horizontal layers are formed, one above the MLD, where the ambient vertical velocities dominate (and particles move passively with the flow) and one below the MLD, where the particles' settling velocity dominates (and particles move actively, relative to the flow); Fig. 4. Below the MLD is where a particle can oscillate due to sinking when the particle+biofilm's density exceeds surrounding seawater density and subsequent resurfacing once the biofilm's loss causes 290 positive buoyancy to be restored (Kooi et al., 2017). At the maximum depth of the oscillation, a short moment of neutral density (the green patches in Fig. 4a, b and c compared to d, e and f). Since smaller particles are more easily mixed above the MLD, the ratio of ambient vertical velocities to the settling velocity is larger (the red patches in Fig. 4d, e and f compared to a, b 295 and c). Areas where the ratio is one show depths at which both the ambient velocity and the relative velocity are important to determine the vertical motion of the particles.

The dominant depth-dependent processes in the biofouling dynamics
We also analyse the results by identifying which of the five biofilm gain and loss terms from Eq. (6) dominate for different depths, particle sizes and regions of our study (Fig. 5). Biofilm growth (G grow ) mostly dominates from the surface down to although not as prominently (Fig 5b). During the spring bloom, G grow dominates down the the base of the EZD, as well as surprisingly between 150 and 250 m (Fig 5b). A potential reason for this is that in subtropical gyres, where less plankton grow at the surface due to nutrient limitation, a deeper primary productivity maximum can be formed where nutrient concentrations are higher.
We can also sum the mass flux of each of the five terms over the entire trajectories for each region, to compare the overall 315 dominance of each term (Fig 6). Since most of the particles stay within the top 50 m in all three regions, and G growth dominates at these depths (as seen in Fig 5), 50% of the one-year biofilm mass flux is controlled by growth. The next largest overall term is L resp , for 30% of the datapoints in the EqPac and SO, and almost 50% in the NPSG. Finally the G coll , L graze and L nonlin are all below 15% in all three regions (with collisions almost at 0).

Characterising oscillations 320
As in Kooi et al. (2017), particles in our model oscillate in the water column. As described above, as soon as a particle sinks below the euphotic zone, algal loss via grazing and respiration dominate the biofilm dynamics. This results in the particle+biofilm eventually reaching a depth where it no longer has a density exceeding that of the surrounding seawater and rising. A point that Kooi et al. (2017) makes regarding different sizes of particles is also supported in our results; the smaller the particle, the lower the frequency of its oscillation. Larger particles have a higher sinking velocity than smaller particles (V s in Eq. (5)) and 325 hence sink and rise faster than smaller particles. This is represented clearly in the EqPac and SO (Fig. 7a and c), where oscilla- tion lengths can reach about 220 days (99th percentile) in the EqPac for particles with a radius of 0.01 mm and about 280 days (99th percentile) in the SO. In the NPSG, since very few particles sink below the mean EZD and therefore do not oscillate, this pattern of smaller particles having longer oscillation lengths than larger particles is not seen (Fig. 7b). The oscillation length for NPSG particles is just under 100 days (99th percentile), probably occurring during the 3 months of the spring bloom when 330 particles do sink below the EZD ( Fig. 5b and e). In all three regions, for particles of 1 mm, the mean oscillation length is less than 10 days. One of the key differences to the Kooi et al. (2017) results is that the oscillations reach much deeper depths. The main reason for this is that the Kooi et al. (2017) study uses a very shallow EZD (around 20 m) and with their limitation of using a constant grazing term at all depths, the biofilm dies and resurfaces very rapidly below the EZD.

335
As with all models, our results depend on our parametrisations, assumptions and model design. The main assumption is that our modelled biofilm only consists of photosynthetic algae (diatoms), while observed biofilm community structures are shown to vary both spatially (e.g. between the Atlantic and Pacific; Amaral-Zettler et al., 2015) and temporally (Amaral-Zettler et al., 2020). In the latter study, diatoms dominate within the first week of colonisation of microplastic, after which other groups, including Rhodobacteriaceae, become more prevalent. Other observations show that diatoms are still one of the most frequently 340 found species after 14 weeks, though the total biofilm species richness is high (Bravo et al., 2011). Furthermore, Amaral-Zettler et al. (2021b) mention that the plastisphere consists of animals and heterotrophic protists that would not necessarily be impacted by the dark and as diatoms die, this could actually stimulate the growth of heterotophs, affecting the particle's buoyancy in Nevertheless, one possibility for future work is to apply a distribution function to biofilm density, growth and death parameters and use a Monte Carlo sampling approach to model these parameters probabilistically. Though this could be computationally demanding, a sensitivity analysis for a simplified scenario could be the starting point. Kreczak et al. (2021) perform a sensitivity analysis using the quotient between algal growth and death rates, for example, and this type of study could be expanded to 350 include the full community diversity of a biofilm.
Some of the biofouling dynamics could also be further developed in future work. For example, we currently assume that the attached algal mesozooplankton grazing rate and growth rate are the same as planktonic algae. We have made this assumption because little is known about the dynamics of plastic biofilms from laboratory settings or in situ observations. Our sensitivity test regarding algal frustules remaining attached to the particle once the biofilm dies (Appendix D) could be improved by using 355 data from biofouling experiments in the dark. One could even test the effects of diatoms entering a dormant phase without nutrients or light. Another aspect to investigate is whether the colonisation or growing of cells could result in a boundary layer effect in which nutrient supply is reduced for biofilm cells below the newly attached cells. Further, the fact that the nonlinear diatom losses in NEMO-MEDUSA include grazing by unmodeled higher trophic levels suggests that the entire plastic plus biofilm could be ingested (instead of individual algal cells within the biofilm, as we assume here). This could be addressed 360 by coupling our biofouling model with other modeled biological interactions such as ingestion (e.g. Cole et al., 2013;Kvale et al., 2021) and egestion of plastic or the merging of particles to model plastic trapped in marine snow. Again, a probabilistic analysis, as described above, could be one approach.
We only use spherical plastic particles since the equations to determine their settling velocity only apply to spheres (Kooi et al., 2017) and a universal model that fits all (micro)plastic types is currently unfeasible to derive (Kreczak et al., 2021). A 365 recent study shows that fibers make up a significant part of plastic particles in aquatic environments and their high surface area to volume ratio leads to a larger area for contact with biology (Kooi et al., 2021). Biofouled ellipsoid-shaped particles could become denser faster than spheres (which have the smallest surface area for any given volume), meaning they could possibly reach the seafloor before the biofilm dies. We hypothesise that upon including different shapes and biological interactions (described above) we might see some vertical distributions such as in Peng et al. (2018), where concentrations of plastic 370 <5 mm are several times higher in the deep-sea (Mariana Trench) than near-surface waters.
Being a process study, we have chosen not to include horizontal advection in our simulations so that the particles do not move away from regions with defined biophysical profiles that we are interested in. Future work with different aims, such as estimating where particles (that are subject to biofouling) end up when releasing them from source locations, would benefit from having 3D advection incorporated.

375
Finally, as mentioned in the Introduction, there have been very few subsurface plastic concentration observations to date.
Overcoming some of the logistical challenges to measure and monitor the 3D movement of plastic in the open ocean would help to validate that: (1) biofouled particles oscillate and (2) their distribution is similar to our modeled results (down to 5000 m) in regions with similar biological and physical properties.

380
We have explored the vertical distribution of ocean plastic spherical particles between 0.01 and 1 mm that are initially buoyant and have been submerged due to biofouling. We present one-year trajectories with more realistic physical and biological dynamics than in the Kooi et al. (2017) and Lobelle et al. (2021) models. The three regions in this study (Equatorial Pacific, North Pacific Subtropical Gyre and Southern Ocean) represent areas in the ocean with different biological activity and windinduced mixing which have varying impacts on the average vertical distribution of particles. In upwelling regions with high 385 productivity and algal concentrations (the EqPac in our study), a subsurface maximum particle concentration is present just below the climatological MLD (which can reach up to five times the concentration of the surface). In regions with very low productivity and low wind-induced mixing (the NPSG in our study), particles remain in the upper few tens of metres and can only sink below the mixed layer if a spring bloom occurs. In areas with very high wind-induced mixing and hence a deep mixed layer (down to 400 m in the Southern Ocean in our case), particles can reach depths of thousands of metres (about 5000 m).

390
This model has been developed in order to gain further understanding of the mechanisms driving the vertical distribution of marine plastic particles. Our model can be incorporated into more sophisticated model studies mapping the total global budget of marine (micro)plastic debris. It can therefore be coupled to models that wish to, for example, add other species to the biofilm community, add sinking due to marine snow and fecal pellet aggregates, as well as how biofouling can affect the weathering Figure A1. Equatorial Pacific's annually averaged vertical velocities [m s −1 ]; stochastic Kz (red) which comprises of tidally induced mixing and wind-driven mixing, potential settling velocity of the particle (green) which is the Vs in Eq. (5) and NEMO-MEDUSA vertical advection (blue). The top row is for particles of the smallest size class (0.01-0.1 mm) and the bottom row is the largest size class (0.1-1 mm). The left column is for particles within the wind-mixing region (above the mixed layer) and the right column is for particles below the mixed layer.
The vertical black line represents a vertical velocity of 0; to the right of that line is for upward velocities and to the left is for downward velocities, each bar indicating an increase by an order of magnitude. of particles and sorption or release of persistent organic pollutants (Rummel et al., 2017). Finally, our model can provide 395 estimations for subsurface concentrations if surface samples have been collected and the full-depth biophysical properties are known.
attached in the dark. Miklasz and Denny (2010) report that the frustule can be denser than the cytoplasm (with a frustule's median density of 1800 g m −3 vs 1065 g m −3 for the cytoplasm). This alternative model means that a particle could keep sinking once the biofilm is dead due to the attached frustules. We simulate that the living algal loss terms (L graze in Eq. (7), L resp in Eq. (8) and L resp in Eq. (9)) is proportional to the gain of mass of the frustule, hence: where v cy is the volume of the cytoplasm [m 3 ] and r A is the radius of an algal cell [m]. The volume of the frustule (v f r ) is 415 therefore v A -v cy and the volume of the whole dead biofilm (v bf dead ) is v f r * A f r * θ pl , where θ pl is the surface area of the plastic particle [m 2 ]. The total volume of the plastic particle plus the biofilm (Eq. (7) in Kooi et al., 2017) is then: v tot = v pl +v bf dead + v bf A , where v pl is the volume of the plastic particle [m 3 ] and v bf A is the volume of the living biofilm (v bf A = (v A * A) * θ pl ).
The equations following these adaptations are then as described in the supplementary information in Lobelle et al. (2021).
As mentioned in the main text, the effect of this alternative model on larger size classes is almost negligible, probably due 420 to this last component, that the radius of the frustule is so much smaller than that of the cytoplasm. The smallest size class, however, is affected due to the fact that the frustules could reach the sizes of the plastic particles and hence affect the density difference between the particle plus frustule and surrounding seawater. It is important to note, however, that since little is known about what happens to biofilms in the dark, we have used some basic assumptions in constructing this model addition and these results should be considered with caution.
like to thank Clément Vic from Ifremer for suggesting the background tidal-induced vertical mixing and providing the script to interpolate Kz