Subsurface dynamics of buoyant microplastics subject to algal biofouling

Mathematical modeling offers an avenue to predict distributions of microplastics in the ocean. With this goal in mind this study presents a simple analytic model describing the deterministic dynamics of long‐time vertical motion of buoyant microplastics that have succumbed to biofouling in a quiescent ocean. The model couples equations governing vertical particle position and algal population dynamics and is formulated in a nondimensional frame work to determine the minimum number of system parameters. We present results which give a comprehensive understanding of the system's key dynamics when varying four parameters. Previous studies have shown that biofouling causes microplastic particles to oscillate below the ocean's surface; however, for the first time we are able to identify the fundamental processes that govern the microplastic trajectories. The model demonstrates that plastic particle properties are the biggest factor in determining the period and characteristics of the oscillation profile, while the algal population dynamics determine the maximum depth reached. As the fouled particles cross the pycnocline their long‐term behavior is to become trapped within a subsurface layer, bounded by the depth of the euphotic zone and the seasonal pycnocline. The smallest particles are extremely sensitive to algal cell attachment and growth suggesting they are always submerged at depths surrounding the base of the euphotic zone or could become trapped in large algal colonies. In general, the results suggest that a higher concentration of biofouled microplastic is expected to be found subsurface, trapped close to the euphotic zone depth rather than at the ocean's surface.

The abundance of plastic pollution in marine environments is well-documented and it is predicted that by 2025 the weight of plastic entering our oceans annually will be between 100 to 250 million metric tons (Jambeck et al. 2015). Subtropical gyres have been identified as "hot spots" for macro-sized buoyant plastic, confirmed by numerical modeling where floating drifters converge to gyre centers due to surface Ekman drift (Lebreton et al. 2012;Maximenko et al. 2012). However, plastic debris has been observed throughout the marine environment in the water column, sediments, and sea ice (Law 2017). Under marine conditions plastic fragments into smaller particles due to photodegradation, chemical degradation and biodegradation (Andrady 2011), and as such microplastics, defined as particles smaller than 5 mm, are ubiquitous throughout the ocean (C ozar et al. 2014;van Sebille et al. 2015;Isobe et al. 2017). Exposure to microplastics adversely affects marine fauna, reducing feeding, growth, and fecundity (Galloway et al. 2017), and evidence suggests that plastic pollution impacts almost all marine ecosystems which service humans, resulting in loss of food security, livelihoods, income and good health (Beaumont et al. 2019). Therefore, it is vital to gain a better insight into how microplastics are distributed throughout our oceans, not only to increase understanding of their impact but to aid in implementing effective plans for their partial removal.
Common plastics encountered in marine environments are low-density polyethylene, high density polyethylene, and polypropylene, which account for 62% of global plastic production (Andrady 2011). These plastics typically have a density lower than seawater and float on the ocean surface. Due to fragmentation, it is expected that there is an exponential increase in particle number with decreasing size (C ozar et al. 2014). However, empirical studies of plastic particles at the oceans' surface reveal that there is a significantly lower abundance of buoyant microplastic particulates than would be expected (C ozar et al. 2014;Eriksen et al. 2014;van Sebille et al. 2015). These observations suggest microplastics are lost to deeper regions of the vertical water column but the mechanisms contributing to the submergence of buoyant microplastics remain unclear (Law et al. 2010;C ozar et al. 2014).
* Correspondence: hannah.kreczak@newcastle.ac.uk This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Additional Supporting Information may be found in the online version of this article. Suggested mechanisms include ingestion by marine biota, for example phytoplankton and zooplankton, and vertical mixing (C ozar et al. 2014;Reisser et al. 2015). Vertical mixing in the ocean arises from wind-driven shear and loss of buoyancy due to surface cooling, the former having been shown to predict submergence of buoyant plastics with exponential decrease in the concentration with depth, inversely proportional to windspeed (Kukulka et al. 2012).
The process of "biofouling," where organisms grow on subsurface objects, has been proposed as another explanation since aggregates of buoyant plastic and marine organisms can become negatively buoyant. Field and laboratory experimental studies have found that biofouling can reach sufficient levels to sink microplastics in both marine (Ye and Andrady 1991;Fazey and Ryan 2016;Kaiser et al. 2017) and freshwater environments (Lagarde et al. 2016;Chen et al. 2019). The time it takes for the aggregates to first become negatively buoyant can be as short as 2 weeks (Fazey and Ryan 2016); however, there are large variations in this time due to geographical location (Kaiser et al. 2017) and seasonal effects (Chen et al. 2019). Once sufficiently fouled to initiate sinking, the fate of microplastics below the ocean's surface remains unclear since it is difficult to observe and measure experimentally. Film plastics which have succumbed to biofouling have been found at depths of up to 400 m (Holmström 1975). Defouling of biofilm colonies can occur, either through light limitation or animal grazing (Ye and Andrady 1991;C ozar et al. 2014), causing the plastic and algal aggregate to become positively buoyant and resurface. Subsurface surveys have found high densities of buoyant microplastics at sampling depths of around 30 m for relatively quiescent ocean conditions (Lattin et al. 2004;Song et al. 2018), contrasting the exponential decrease with increased depth predicted by wind mixing. These observations lead us to speculate that biological mechanisms result in microplastics becoming trapped subsurface, potentially accumulating at intermediate depths in the upper ocean column.
Development of mathematical models of microplastic biofouling in the oceans offer an avenue to predict distributions of microplastics by coupling to 3D general ocean circulation models (OGCMs). The fates of microplastics with fixed density (Jal on-Rojas et al. 2019; Mountford and Morales Maqueda 2019) and linearly increasing density with time (Jal on-Rojas et al. 2019) have been investigated in 3D ocean circulation models. Jal on-Rojas et al. (2019) concluded that density and biofilm thickness have the biggest influence on microplastic debris paths and fates. Kooi et al. (2017) devised the first deterministic model coupling microplastic hydrodynamics and biofilm population dynamics, both of which varied with vertical position and time. They found that biofouling can result in particles oscillating vertically in the upper ocean column. This study was a huge leap forward but had limitations in that the complex model restricted the ability to perform an extensive parameter study and they were unable to comment in detail on the underlying physical and biological processes that drive the observed behaviors. Recently, Lobelle et al. (2021) have incoporated the Kooi model of biofouling into an OGCM, identifying initial sink times and depth reached on first sink for fouled microplastics dependent of geographical location, season, and particle size.
Here we build on the Kooi et al. (2017) study and address the need for a canonical model, deriving a nondimensional dynamical system capturing the key dynamics of vertical particle motion, due to algal biofouling. We shall demonstrate this leads to a system whose dynamics can be fully understood through a detailed parameter study. The model is easily generalized to different ocean, freshwater, and algal properties under latitude and seasonal variations. The results are organized to investigate biofouled particle paths in both well-mixed and stratified oceans. Analytic solutions of the dynamical system highlight the main features of, and contributing mechanisms to, the fouled particles long-time motion going far beyond what was possible in the Kooi et al. (2017) model.
For a well-mixed ocean and when ignoring background algal attachment as a limiting case in the biofilm population dynamics, particles oscillate within the upper ocean vertical column and are observed to briefly resurface in their motion. Introducing either a vertical variation in ocean density, or allowing algal cell attachment from background algal concentrations, results in trajectory switching and particles do not return to the surface. Instead, long-time particle motion is restricted to oscillations within a mid-depth region in the upper ocean, centered about the euphotic zone depth, which we call the "subsurface plastic trapping layer." An important caveat to this conclusion is that ocean dynamics are neglected in this study, which will be addressed elsewhere.

Model development
The process of biofouling from marine algae can be described as follows. A buoyant microplastic particle which is part submerged at the ocean's surface is exposed to algae which attach and grow on the particles' surface as a biofilm. The size of this biofilm is determined by both the rate of attachment and the relative growth and death rates of the biofilm population. Algal density typically exceeds that of seawater, so once a biofilm is sufficiently large, combined plastic particle and algal aggregates will become negatively buoyant and begin to sink. The settling velocity is determined by the relative buoyancy and the size of the biofouled particle. Once the particle is submerged below the euphotic zone, algal growth from photosynthesis is not possible due to limited light exposure. Algal population in the biofilm reduces due to respiration and fauna grazing, implying that the biofouled particle can become positively buoyant once more, rising toward the ocean's surface. Rising back into the euphotic zone and experiencing an increase in light intensity means the algae can photosynthesize and grow once more. Since algal growth rates are dependent on the particle's position in the vertical oceanic water column, and settling velocity is dependent on the fouled particle size and relative buoyancy which changes with depth, this process can be captured by the coupled system where z is the vertical position of the biofouled particle, N is the population number of algal cells in the biofilm, and t is time. We choose the origin z = 0 to represent the mean sea surface with positive z axis directed vertically upwards. The choice of functions H and G, and the nondimensionalization, are now briefly described with a more detailed discussion left to the Supporting Information S1. There has been much effort dedicated to determining the settling velocities of microplastics of varying shapes and sizes (Dioguardi and Mele 2015;Khatmullina and Isachenko 2017;Waldschläger and Schüttrumpf 2019) including those which have succumbed to biofouling (Van Melkebeke et al. 2020). Since the drag coefficient varies with both Reynolds number and particle shape, approximating settling velocities for the full range of microplastics is complex and no universal relation has been found (Khatmullina and Isachenko 2017). Therefore, it is impossible for us to currently derive a universal model of biofouled particle transport dynamics for all microplastic types. In the model, we approximate the fouled particles settling velocities by assuming they are small, spherical particles subject to Stokes flow in a stratified ocean. The vertical displacement is then given by where g is the gravitational acceleration (m s À2 ), ν is the kinematic viscosity of seawater (m 2 s À1 ), and ρ f ðÞ is the ocean density profile (kg m À3 ). The factor D = 86,400 (s d À1 ) is included such that the velocity timescale aligns with the algal population net growth rate, commonly measured in the literature with units d À1 (Butterwick et al. 2005;Sandnes et al. 2005). The settling velocity changes due to the relative mass difference between the fluid displaced by the fouled particle and the total mass of the clean plastic particle with a biofilm consisting of N algal particles. The radius a B is the total radius of the fouled particle (m), which we approximate as constant for simplicity, the choice of which given in the Supporting Information S1. The background fluid density increases with depth and can be represented by ρ f (z) = ρ 0 + Δf(z), where ρ 0 is the density at the ocean surface (z = 0), Δ = ρ b À ρ 0 and ρ b is the density at the ocean floor (z = Àz b ). The function f(z) captures the normalized profile of ocean density variation, with f(0) = 0 and f (À z b ) = 1. Here we take the simplest approximation for an ocean density profile, the piecewise constant function capturing the seasonal pycnocline at a depth h. To avoid the discontinuity and reduce complexity in solving the ordinary differential equations (ODE) numerically, the function (3) is smoothed, details of which are given in the Supporting Information S2.
The algal population equation G in (1) is given by a simple population model (Kooi et al. 2017) where λ A is the algal reattachment rate, λ G the growth rate, and λ D the death rate, all with units d À1 .
The growth rate λ G (z) is dependent on algal photosynthesis, the rate of which is governed by light intensity. Light intensity in the ocean is depth dependent and varies according to the well-known Beer-Lambert relation I z ð Þ ¼ I S exp ϵz ð Þ (Béchet et al. 2013). The euphotic zone depth z E reflects the depth where light intensity is 1% of its surface value I S (W m À2 ) and is given by is dependent on the conditions of the ocean. The vertical relation for the growth rate we represent as where λ max G (d À1 ) is the maximum growth rate and α (d À1 /μ moles m À2 s À1 ) governs how the growth rate changes with light intensity. Details on the choice of (5) for the modeling can be found in the Supporting Information S1.
Algal growth rate is also temperature and nutrient dependent. For simplicity, we shall assume that there are sufficient nutrients for algal growth at all depths since fouling is able to occur in the first place. Algal growth is optimum at a specific temperature (which is species dependent), but the growth rate decreases with temperatures greater and lower than this optimum (Dauta et al. 1990;Sandnes et al. 2005). Models of light and temperature dependence on algal growth use only a function fit to approximate temperature dependence (Bernard and Rémond 2012;Béchet et al. 2013). We assume ocean temperature effects on algal growth are either negligible within the upper ocean of interest, or since temperature decreases with ocean depth, it could be assumed that the function g z ð Þ captures the negative impact of increased depth on algal growth rates for both light and temperature dependence, assuming the ocean surface temperature is equal to or less than the optimum temperature.
For the algal attachment rate, we adopt where β A is a collision kernel rate (m 3 s À1 ) of a particle encountering a single algal cell, c(z) is a normalized function capturing the vertical profile of the background algal cell concentration, and A max is the maximum value of algal concentration. The choice of c (z) is a Gaussian centered about the depth z max , where the maximum concentration occurs, with standard deviation δz.
A decrease in the biofilm algal population can occur from respiration or fauna grazing, captured in the negative death rate À λ D . For simplification the rate λ D is assumed to be independent of depth and held constant.
Equations (2) and (4) describe a coupled system for a simple dimensional model of microplastic transport from algal biofouling. Even as a simplified model there are several governing parameters for microplastic properties, algal population dynamics, and stratification, many of which change temporally with latitude and seasonal variations. We therefore proceed to nondimensionalize (2) and (4) to determine the minimum number of dimensionless parameters that govern the system. The nondimensionalized process is detailed in the Supporting Information S1, but summarized here for brevity. We nondimensionalize the system with z ¼z=ϵ , are the algal population number at which neutral buoyancy occurs at the ocean surface (no.), and the rise velocity of a clean plastic particle (m d À1 ), respectively. This rescaling yields the fully nondimensionl system In (8a), the nondimensional parameters are Clearly, R 1 ( 0), is the oceanic stratification parameter. The dimensionless parameter R 2 ( 0) is the quotient of the pycnocline density jump and the density difference between clean microplastic and the surface ocean density. Similarly, R 3 ( 0) is the quotient of the pycnocline density jump and the density difference between algal and the surface ocean density. In (8b), the nondimensional parameters are These parameters capture the ratio between timescales of the system, for example the timescale of background algal cell attachment is β A A max , while the particle motion timescale is W * ϵ, that is the time it takes for a clean particle to traverse the distance 1/ϵ. Therefore, the parameter Ω is the quotient between algal attachment and particle motion timescales. Φ is the quotient between the algal population timescale, parameterized here by the death rate λ D , and the particle motion timescale, while Λ is the quotient between algal growth and death rates.
In addition to the six dimensionless parameters given in (9) and (10), there are four parameters governing the normalized, vertical-profile functions gz ð Þ, fz ð Þ, and cz ð Þ. These are Ψ, as defined in (5), H = ϵh, which is the dimensionless pycnocline depth, and dimensionless parameters which control the Gaussian representing the background algal concentration, centered on Z max = ϵz max with standard deviation Δz = ϵδz. The realistic ranges for each of the 10 dimensionless parameters are given in Table S1 in the Supporting Information, interpreted from the oceanographic variations of their individual components.
In an effort to reduce the parameter space and improve intuition in the model behavior, we set This sets the condition that above Àz E the algae can photosynthesize and grow, while below Àz E the light intensity is weak and algal population dynamics are dominated by mortality resulting in biofilm mass reduction. Uitz et al. (2006) report high concentrations of phytoplankton about the euphotic zone depth with a wide distribution, therefore we choose to fix Z max ¼ Àz E and choose Δz = 1 to be sufficiently wide to match observations. The functions fz ð Þ, gz ð Þ, and cz ð Þ are plotted in Fig. 1 for typical oceanic parameters. The parameters R 1 -R 3 are chosen to be fixed as the median values of those listed in Table S1 of the Supporting Information, R 1 = 0.0015, R 2 = 0.018, and R 3 = 0.024. Understanding the microplastic trajectories as a function of the four remaining parameters, Φ, Λ, H and Ω, is the focus of this study.
The system is initialized withz 0 ¼ 0 andÑ 0 ¼ 1, when the particle initially starts to sink. Numerical simulations are performed using a Runge-Kutta-Fehlberg scheme via Python Scipy, with an adaptive timestep and a minimum time step of dt = 10 À3 .

Short-time dynamics for an ocean of constant density
The parameter Ω, which sets the importance of algal cell attachment, can be written in terms of Φ as Empirical studies of biofouling on plastics have found that it usually takes over 2 weeks before plastics are sufficiently fouled to start sinking (Fazey and Ryan 2016), suggesting N * ) β A A max , where β A A max represents the maximum number of algae attached daily. Estimate values for λ D are between 0.05 and 0.75 (O'Brien 1974), therefore it is a reasonable to assume that Ω Φ. We first proceed assuming Ω ( Φ andÑ remains greater than or O 1 ð Þ, thus setting Ω = 0 gives an approximation of the systems short-time dynamics. The final results section addresses the validity of this assumption with a computational study.
When Ω = 0, an analytic implicit solution for the system trajectories can be found for constant density profiles. Since the system is autonomous, the solutions can be described by For both fz ð Þ ¼ 0 and fz ð Þ ¼ 1 , (12) is separable and has solutions respectively, where c 0 and c 1 are constants. The remainder of this results subsection will focus on an ocean of constant density such that fz ð Þ ¼ 0 for allz, which is an appropriate approximation of the density profile in well-mixed shelf sea regions. The trajectory initialized byz 0 ¼ 0 andÑ 0 ¼ 1 is found from (13a) as  Figure 2b plots both the particle vertical position and population growth in time giving a better interpretation of the physical processes in the system. The periodic orbit in the phase space implies that the fouled plastic particle oscillates vertically in the ocean. The fouled particle begins to sink att ¼ 0 withÑ ¼ 1 and while it remains in the euphotic zone the algal population continues to grow such thatÑ > 1. As the particle cross the euphotic zone depth, the biofilm can no longer support photosynthesis and algal mortality dominates the population dynamics. While submerged below the euphotic zone depth the excess algal population dies off until the density of the biofouled microplastic matches that of the surrounding oceanÑ ¼ 1 À Á and begins to rise again. While rising the algae continues to die offÑ < 1 À Á before the particle crosses back into the euphotic zone and the biofilm can begin to grow once again. The particle briefly resurfaces before the motion begins again. Now that the basic oscillatory motion of fouled particles is understood, we consider the effects of the parameters Φ and Λ with Ω = 0 and fz ð Þ ¼ 0 for allz.
For six values of Φ, Fig. 3 shows three oscillatory periods of the vertical displacement for fouled particles with fixed Λ = 4. Fig. 2. (a) Solutions of the dynamical system (8) are shown in phase space (colored) when Φ = 0.1 and Λ = 4, and the trajectory (black) for the initial conditionÑ,z À Á ¼ 1,0 ð Þ. In (b), the black trajectory in (a) is shown as the time evolution of particle vertical displacementz (blue) and algal populationÑ (orange). The line Àz E ≈ À 4:605 is plotted with thez axis to show when the particle crosses the euphotic zone boundary and the algal growth rate changes from positive to negative.

A B
The dimensionless parameter Φ varies over orders of magnitude for typical microplastic sizes (see Table 1; Supporting Information S1). Small Φ corresponds to the largest microplastic particles, since it is expected they have fast rise velocities. Small particles have a large Φ, since their rise velocities are small due to increased viscous forces, and it takes them a long time to cover the distance 1/ϵ. For small Φ the oscillatory behavior appears to be regular with not much qualitative difference between the vertical displacement behavior across values of Φ. Note the change in timescales, where for the smallest value of Φ the time taken to resurface is extremely large in dimensionless timet.
As Φ is increased the oscillatory period gets shorter and appears to plateau for the largest values of Φ. The sink and rise behaviors become highly asymmetric, where the particle appears to plummet from the ocean's surface after initial fouling, but once it begins to rise again the motion is extremely slow. The dashed red lines show the relationz /t which indicates that the particle is traveling at its clean rise speed. Further insight into this can be gained by looking at the algal population dynamics. Figure 4 shows the changes in the dimensionless algal population in time for the same parameter values presented in Fig. 3. For small Φ there is little variation in the relative population aboutÑ ¼ 1. As Φ increases the changes in algal population become more significant over the oscillatory period, achieving large maximum values over the initial value when the biofouled particle begins to sink. Most notably the population appears to almost die off withÑ ! 0, explaining why the particle achieves its clean rise speed on ascent. For the largest value of Φ, the algal population is around 80 N * implying that the microplastic would be engulfed by algae on its descent from the ocean's surface.
Recall that these observations arise from an idealized model so here we make clear how the results may vary for real microplastic and algal aggregates. First of all, the model is continuous in algal population numberÑ , when in fact such a process is discrete in nature. When Φ is large, corresponding to the smallest microplastic particles, only a small N * would be needed to initiate sinking. Therefore,Ñ ! 0 would likely imply that the algal population dies off completely and the plastic can only be refouled by algal attachment, which has been neglected in our modeling framework so far. The regrowth of biofilm as the particle enters the euphotic zone in the simulations is due to the dimensionless algal population achieving only numerically small values ofÑ rather than exactly 0. It is also important to note that for large Φ, such large values of N * also suggest that the assumption the aggregate radius a B remains roughly constant no longer applies and as a consequence the velocities in the simulation will be over approximated. The real settling velocities are expected to be slower, resulting in shallower particle submergence than simulated.
One final assumption in the model formulation is that the model assumes that once the algal cells die, they are removed from the biofilm and no longer effect the dynamics. While Φ is small (larger microplastics), when it is expected that N * is large, having some lingering deceased algal cells is unlikely to affect the dynamics significantly. For the smaller particles with large Φ, where the microplastics can become completely engulfed in algae during its slow descent into the ocean depths, remaining dead algal cells may trap the microplastic. Dead algal clusters fall to the ocean floor as marine snow (Sanders et al. 2014), therefore the smallest microplastics may  descend as marine snow in aggregates dominated by algae. This mechanism may result in buoyant microplastics reaching the deep ocean floor where sedimentation has been hypothesized to be a major sink for marine microplastics (van Cauwenberghe et al. 2015).
While it was seen that changing Φ by several orders of magnitude significantly changed the asymmetry and time period of particle vertical oscillations, in Fig. 3 we observe that changing Φ has no effect on the maximum depth reached by the biofouled particle. This is only achieved by changing the dimensionless parameter Λ, which captures the ratio between the maximum growth rate and death rate of the algae. Figure 5a shows how the vertical position in time changes with different values of Λ. As Λ is increased the maximum depth the particle reaches, which we denote z D , becomes deeper. The explanation for why larger Λ achieves deeper particle submergence can be understood by considering fixing λ D and increasing λ max G . While in the euphotic zone, a larger growth rate λ max G implies that the algal biofilm population grows larger, as seen in Fig. 5b, resulting in a larger negative buoyancy for the fouled particle. The particle can travel faster and reach larger depths, but also a longer submergence time below the euphotic zone is needed for neutral buoyancy to be achieved, after which the fouled particle will rise again. There are slight variations in the oscillatory period when changing Λ, which can be observed by looking at the time it takes for the particle to resurface after one oscillation, which we denote by T. Figure 6a plots the dimensionless time-period of one oscillation T when changing Φ where each of the lines corresponds to a different value of Λ. It is seen that the most significant change in the time period arises from variations in Φ, rather than Λ, which for Φ ≤ 10 À1 varies as T / Φ À1/2 represented by the black dashed line. The relation corresponds to values of Φ where the oscillatory behavior is qualitatively similar in that the motion of the particle appears symmetric in the sink and rise phases. In dimensional time this equates to a time period relation t T / W Ã ϵ ð Þ 3=2 implying the smallest particles take significantly long times to resurface due to their slow rise speeds. A mathematical reasoning for this scaling is presented in the Supporting Information S3. Figure 6b shows how the maximum depth z D reached by the biofouled particle varies linearly with Λ. A line of best fit is plotted in dashed black with relation z D = 2.412Λ + 4.415 to 4 significant places. The implicit solution (13a) can be used to analytically predict the linear relationship, details of which are given in the Supporting Information S4.
We briefly discuss the dimensional time period t T ¼ Tt and the maximum depth reached z ¼ z Dz by an oscillatory biofouled particle in an ocean of constant density. For example, taking z E = 20 m, which gives ϵ ≈ 0.23, and median values of the oceanic parameters listed in Table S1, the dimensionless parameter Φ can be approximated from The dimensional time period can then be approximated from t T = (W * ϵ) À1 T = (Φ/λ D )T. Table 1 lists a handful of the time periods T from Fig. 3 and gives an approximate average dimensional timescale t T (d) for a single oscillation of a particle with radius a B . These values are on the order of those observed in the model of Kooi et al. (2017) for particles of the same size. The depth z D = 14.144 reached in Fig. 3 implies a maximum depth reached in the ocean is around 60 m when the euphotic zone depth is z E = 20 m. This approximation of the euphotic zone depth would be for an ocean with poor light penetration, and a euphotic zone depth of up to 80 m would imply fouled particles could reach depths of 245 m. At such depths, in shallow seas plastic particles and algal aggregates are likely to reach the sea floor and become trapped in sediment as a mechanism contributing to their removal from the water column. For example, average depth of the central North Sea is only 74 m (GEBCO Compilation Group 2020) and most biofouled microplastics would reach the sea floor.

Subsurface oscillations in a stratified ocean
In the last section, density variations in the ocean column were neglected by assuming fz ð Þ ¼ 0 for allz . In this section we introduce density variations but continue to neglect algal attachment (Ω = 0).
For the presented results R 1 , R 2 , and R 3 remain fixed with values assigned as stated in the Model Development. We set Λ = 4 and choose a couple of values of Φ, the only remaining parameter governing the dynamics (8) is H, which represents the dimensionless pycnocline depth. Figure 7 shows how the fouled particles' vertical position changes in time for the two pycnocline depths (a,b) H = 6 and (d,e) H = 2, and two values of Φ. The pycnocline depth is marked by the solid orange line while the euphotic zone depth is marked by the dashed black line. In Fig. 7a,b, when H >z E , for both Φ = 1 and Φ = 0.01 the particle continues to oscillate about the euphotic zone depth; however in the long-time limit, the amplitude of the oscillatory motion is reduced until the motion is bounded below by the pycnocline depth H. It takes many more oscillations for Φ = 1 to reach this long-time behavior than for Φ = 0.01, with the latter settling into the new regime after just four oscillations. In Fig. 7d,e, H <z E , and the oscillatory amplitude is also reduced, except now the pycnocline depth is an upper bound on the long-time oscillatory motion.
The long-time reduction in the amplitude of the oscillatory motion from interaction with the pycnocline can be understood dynamically by looking at the phase space for N,z À Á plotted in Fig. 7c,f. Recall that two analytic implicit functions are found for both fz ð Þ ¼ 0 and fz ð Þ ¼ 1 , Eqs. 13a and 13b, respectively. While the particles vertical position satisfiesz > À H , the trajectories are determined to follow solutions to (13a), but whenz < À H , the trajectories follow solutions of (13b). At the boundaryz ¼ ÀH there is a switching between regimes which causes the perceived reduction in amplitude for the vertical position. Eventually, the particle settles onto a trajectory which no longer crosses the boundary, and whether this is above or below the pycnocline depends on the relative positions of the euphotic zone and pycnocline depths. Figure 7c shows the trajectory (black) for the system when Φ = 0.01 and H >z E . The colored periodic orbits represent solutions to the analytic expressions above and below the pycnocline and the trajectory follows the respective solutions depending on whetherz > H orz < H, until it settles onto a single orbit in the long-time limit bounded below by the pycnocline. The long-time orbit is the largest orbit which circles the euphotic zone fixed point 1,Àz E ð Þ without crossing the pycnocline. The implicit, analytic curve which defines this orbit is given by (13a) which passes through the pointÑ,z À Á ¼ 1,ÀH ð Þ. In Fig. 7f, H <z E , and the long-time orbit is bounded above by the pycnocline since it encircles the euphotic zone depth without crossing the pycnocline. The implicit, analytic solution for this long-time orbit is given by (13b) which passes through the pointÑ,z À Á ¼Ñ Ã 1 ,ÀH where the dimensionless value at which dz=dt ¼ 0 when fz ð Þ ¼ 1. Fig. 6. In (a), the time period of one oscillation presents a power-law relation T $ Φ À1/2 , while (b) shows a linear relationship between algal growth and death rate ration Λ and maximum depth reached z D .

A B
The physical interpretation of these results is as follows. Once the biofouled particle crosses the pycnocline its velocity slows down due to the decrease in relative buoyancy. The algal population number needed to achieve neutral buoyancy is now larger than above the pycnocline. When the particle is submerged below the euphotic zone and the algal population starts to decrease, the particle will begin to rise sooner due to this new condition. Once it has crossed back into the euphotic zone the fouled particle weighs more than when it initially left the euphotic zone. When the biofilm begins to grow again less time is needed before neutral buoyancy is reached to initiate sinking and the particle does not resurface. This process is repeated until eventually the particle resides either above the pycnocline for all time, or below it for all time. These results predict that when biofouled microplastics interact with large vertical density gradients in the ocean they do not resurface and become trapped at subsurface depths in the upper ocean column. Therefore, we hypothesise the existence of a "subsurface plastic trapping layer" where the seasonal pycnocline depth is shallow, such as in equatorial and tropical regions. For a given Λ it is required that H < z D for the fouled particle to interact with the pycnocline, such that they no longer resurface. If H < z D does occur, changes in parameters Λ, R 1 , R 2 , and R 3 vary the oscillatory time period of the long-time bounded oscillations.

Subsurface oscillatory kick from algal attachment
In this section we present example results with Ω > 0 to show how algal attachment can impact the short and longtime dynamics. For ease of discussion the ocean is once again modeled as well-mixed, setting fz ð Þ ¼ 0 for allz. Figure 8 shows the vertical displacement of fouled particles with both Ω = 0 and Ω = 10 À3 for (a) Φ = 10, (b) Φ = 1, and (c) Φ = 0.1, and the growth and death rate ratio fixed as Λ = 4. The normalized concentration function cz ð Þ parameters are as given in Fig. 1. In Fig. 8a, when Φ = 10 the particle initially sinks to the full depth before ascending at the clean particle rise speed. Once it crosses into the euphotic zone the particle does not ascend to the surface and instead oscillates with small amplitude close to the euphotic zone depth with a shortened oscillatory time period when Ω = 10 À3 . For Φ = 1 in Fig. 8b similar dampening behavior occurs and the particle does not resurface; however after an initial significant loss in oscillatory amplitude, the amplitude gradually decreases with further time. In Fig. 8c, when Φ = 0.1 there is no drastic reduction in the oscillation amplitude after the first descent although on closer inspection there is a slight, gradual reduction in the oscillation amplitude with increasing time.
The dynamic change in behavior is due to the fact that when Ω > 0, the fixed point of the system, which isÑ,z Þ , is a stable spiral assuming Ω ( Φ (details in the Supporting Information S5). Trajectories in the long-time limit are continually dampened by algae attachment until they reside at the fixed point indefinitely, close to the point 1, Àz E ð Þ, thez value of which is plotted as the dashed black line in Fig. 8. Comparing the phase space trajectories of the system (8) with Ω = 0 and Ω > 0 reveals the dynamic mechanism contributing to the initial and long-time dampening from algal attachment. Figure 9a plots the trajectories (black) in phase space for the fouled particle trajectory with Φ = 1, Ω = 10 À3 for times t 0, 100 ½ . Solution curves (blue-green) to the implicit relation (13a) when Ω = 0 are also plotted and show that for N ) 0 the trajectories follow the solution curves. Figure 9b shows a portion of the phase space for smallÑ . During the first oscillation,Ñ ! 0 while the particle is below the euphotic zone. As the particle begins to rise and encounter the peak in algal concentration there is an increase inÑ and the trajectory is kicked from the solution curve onto another, which it follows for the remainder of the oscillation period resulting in a reduction of oscillation amplitude. When the trajectory approaches small values ofÑ during the particle rise phase another kick occurs but the reduction in the amplitude is smaller than the first. This continues for all oscillation periods with a decreasing effect on amplitude reduction. From Figs. 8 and 9, it appears that the system is more sensitive to algal attachment when Φ is large, corresponding to smaller microplastic particles. The results suggest algal attachment is most significant whenÑ ≈ 0 , evidenced by the decreasing impact of algal attachment on reducing the amplitude of the oscillations after the initial kick, most notable in Fig. 8a for Φ = 10, Ω = 10 À3 . In Fig. 8c when Φ = 0.1, Ω = 10 À3 , there is some reduction in the oscillation amplitude, but it is minimal. Comparing with Fig. 4, we see that indeed algal population in this case does not approach 0 when Ω = 0 and therefore less sensitive to algal attachment.
Although the mechanism of algal attachment differs from the impact of ocean density variations, it also results in particles becoming trapped subsurface in the particle trapping layer, particularly for the smallest particles. When both mechanisms are included in the simulations the pycnocline still bounds the subsurface long-time motion for all particles, with algal attachment producing additional dampening for the smallest. Note that the model presented is idealized and as such stochastic variations in algal population dynamics is expected, likely counteracting the minimal dampening from algal attachment experienced for the largest particles.

Discussion
We have studied the behavior of a nondimensional dynamical system capturing the dynamics of a buoyant microplastic which have been biofouled by algae in an ambient ocean. In a well-mixed ocean, biofouled particles were shown to oscillate vertically, the depth of which is determined by algal growth properties and their sensitivity to the penetration of light intensity. The time-period of oscillation is mainly determined by microplastic particle size which can vary over orders of magnitude. Smaller particles experience extremely slow descents and as such can become engulfed by an algal colony, suggesting that the smallest microplastics may become trapped in marine snow and descend to the ocean floor. The  Fig. 9. The particle trajectory for initial condition (1, 0) (black) is shown in phase spaceÑ,z À Á fort 0,100 ½ , with Φ = 1, Λ = 4, and algal concentration parameters Ω = 10 À3 , Z max ¼z E and Δz = 1.0. Solutions of the dynamical system (8) with Ω = 0 are also plotted in phase space (colored). In (a), the full phase space shows the trajectory is no longer a closed curve when Ω > 0. In (b), a zoomed in section nearÑ ¼ 0 reveals that as the trajectory encounters the background algal concentration, it is kicked to another solution curve of the system when Ω = 0.

A B
periods of oscillation were seen to range from several hours to tens of days.
For a stratified ocean, interaction with ocean background density changes resulted in fouled particles not returning to the ocean's surface, instead their long-time behavior was to oscillate within a subsurface plastic trapping layer. While jumps in density bounded the oscillatory behavior in a stratified ocean, algal reattachment in a well-mixed ocean resulted in continued dampening and an attraction of trajectories to the fixed point that lies close to the euphotic zone depth, particularly for the smallest microplastics. Kooi et al. (2017) reported similar subsurface oscillations but were unable to comment on the contributing mechanisms. These results suggest that microplastics are distributed throughout the upper-ocean, subsurface vertical column with increased density likely below the surface in a region which we have referred to as the subsurface plastic trapping layer, positioned about the euphotic zone depth. This aligns with the observations of microplastic distributions from oceans with relatively quiescent conditions where larger densities of microplastics were found at depths of around 30 m rather than the oceans' surface (Lattin et al. 2004;Song et al. 2018). The maximum depth fouled particles reached in the model suggests that in coastal and shelf seas regions microplastics and algal aggregates have the potential to reach the ocean floor and be trapped through sedimentation if the euphotic zone depth is deep enough.
The euphotic zone depth can vary from 5 to 80 m, changing with geographical location and season (Lee et al. 2007), as does the depth of the pycnocline. The results of the model herein suggest that the relative relationship between the euphotic zone depth and pycnocline determine the amplitude of microplastic and algal aggregate long-time oscillatory behavior in the trapping layer. In equatorial regions, the depth of the pycnocline is shallow and therefore may bound from above microplastic debris. At midlatitudes the depth of the pycnocline is much deeper (Talley et al. 2011), and more likely to bound the trajectories of fouled particles from below.
Effects of ocean temperature, which impacts both stratification and biological processes, are neglected in our model. Geographical variation and seasonal and diurnal temporal variation in temperature influence the depth of the pycnocline. Algal growth rates are also temperature dependent (Sandnes et al. 2005), leading Kooi et al. (2017) to investigate the impact of geographical location on the onset of sinking in their deterministic biofouling model for a single algal species. They reported that particles did not settle in the North Atlantic for their chosen parameters as a result of low temperatures. Similarly, Chen et al. (2019) observed that biofouled microplastics in freshwater did not become sufficiently fouled to initiate sinking during winter months. As a biofouled particle oscillates in the ocean vertical column, algal growth rate may also be sensitive to vertical changes in temperature. The impact of geographical and temporal temperature changes on stratification and algal growth rates would be an interesting avenue to pursue in the future.
Also neglected in our model are variations in light intensity associated with the diurnal solar cycle. The parameters in the model relating to light intensive growth can be taken as the average values when considering time variations in the surface light intensity where t T is the dimensional time period of one fouled particle oscillation. Since I S t ð Þ is periodic, a larger t T implies I S is better approximated by (16). Since t T /Φ 3/2 , it is expected for smaller particles, which have large Φ, that the average light intensity (16) is an accurate approximation when calculating the particle trajectories. For larger particles, for which Φ is small, we may expect to see deviation from the orbits predicted by the model and as such modulation in the vertical displacement. This is indeed the case in the deterministic model of Kooi et al. (2017) who include diurnal variation in light intensity and observe aperiodicity in the vertical oscillations of the larger particle sizes, although they assigned no contributing mechanism.
Beyond the scope of this study was investigating the impact of fluid motion on the settling velocities of the particles. Fluid motion in the ocean mixed layer (OML) is complex, with stirring arising from winds, wave breaking and buoyancy fluxes at the ocean surface, while large turbulent eddies drive entrainment at the bottom (Kantha and Clayson 2014). Wind stress at the surface is the primary mixing agent and has already been attributed to the sinking of buoyant microplastics (Kukulka et al. 2012). In addition, Langmuir cells are organized rotating cells in the surface layer and are important to upper-ocean mixing and transport. Vorticial motions have the potential to counteract the inertia of falling, dense particles implying that continued suspension of biofouled microplastics near the ocean surface could occur (Maxey 1987). If biofouled particles are suspended in regions of high light intensity, larger algal colonies can grow. Eventually the density will become large enough for the aggregates to escape the wind mixed vorticial motions and sink as predicted by the model here but reaching larger depths. Assessing the impact, the range of organized and turbulent motions of the OML has on the fate of fouled microplastics would allow for a more complete picture of microplastic vertical distributions in the ocean. This will be addressed in a future study.