The Dynamics of Buoyant Microplastic in the Ocean Forced by Unsteady Insolation

Abstract

Building on the deterministic mathematical models of Kooi et al. (Environ. Sci. Technol. 51, 2017) and Kreczak et al. (Limnol. Oceanogr. 66, 2021), this study investigates the trajectories of biofouled microplastic particles forced by unsteady insolation. A new, non-dimensional system of governing equations is derived to predict the particle trajectory in a stratified, quiescent ocean subject to unsteady insolation. In the absence of stratification, unsteady insolation drives surface-to-depth oscillations with amplitude varying seasonally, attaining a maximum/minimum in the summer/winter, respectively. At high latitudes, a particle spends an increasing length of time floating on the sea surface in the winter when biofilm production is minimal or absent altogether. We demonstrate that, at 70N, the oscillations are modulated; in summer they are briefly subsurface, while in spring/fall they reach the sea surface and exhibit the largest amplitude throughout the year. In contrast, forcing the particle motion with constant, annually averaged insolation, at any given latitude, always generates persistent surface-to-depth periodic oscillations. In a stratified ocean, the previously reported persistence of subsurface particle oscillations forced by constant insolation is no longer exhibited for unsteady solar forcing. At lower latitudes, surface-to-depth oscillations with seasonally varying amplitude occur in a stratified ocean. In polar latitudes, the particle dynamics displays three regimes: (i) floating at the sea surface in winter, (ii) surface-to-depth oscillations in spring/fall with time-varying amplitude, (iii) subsurface oscillations around the compensation depth, where biofilm production and mortality rates balance. Decreasing the particle size leads to longer oscillation periods, and at high latitudes the particle either floats or performs subsurface oscillations with seasonally varying amplitude about the compensation depth.

1. Introduction

Marine microplastic pollution is of major societal concern because of its impact on the functioning and health of marine ecosystems [1,2]. Macro-sized plastics entering the ocean break up into microplastics, defined as particles between 1 and 5000 $\mathsf{\mu }$m in size, which have been found throughout the ocean [3]. This fragmentation occurs when surface plastic is exposed to solar radiation, causing photodegradation, as well as via chemical degradation, biodegradation and mechanical fracture due to the action of waves and turbulence [4,5]. Microplastics can fragment further into nanoplastics (<1 $\mathsf{\mu }$m), which are potentially the most hazardous form of marine litter [6].

The increasing levels of microplastics in the ocean are having a significant impact on marine life. Microplastics can be ingested by organisms such as zooplankton, which provide an essential trophic link in the marine ecosystem [7]. Once microplastics have accumulated in the body tissue of marine organisms, a reduction in reproduction and feeding has been observed, affecting growth rates [8]. Microplastics can also attach to water pollutants such as dyes, heavy metals and other chemicals, causing further harm when consumed [9]. Contamination of these aquatic organisms provides a pathway for microplastics to enter the human body via seafood consumption. We have yet to grasp the full impact that microplastics can have on human health [10]. One positive effect of plastic in the oceans is that they have been found to create new habitats for some marine species [11]. Given the negative impacts of microplastics in the marine environment, what realistic technologies are available for their removal? A starting point to address this question is to identify microplastic accumulation “hotspots”, either in the water column [12,13] or on the seabed [14,15,16], with the aim of targeting these sites first.

It is estimated that 62% of plastic produced is less dense than sea water [17]; however, this is not reflected in the estimated amount of microplastic floating on the ocean surface. Instead, the microplastic is distributed in the water column and on the seabed. All material entering the ocean is subject to biofouling, where micro-organisms consisting of bacteria, algae and crustaceans (collectively referred to as a “biofilm”) grow on the surface. Microplastics are no exception, and when subject to biofilm growth, clean microplastic with density less than that of seawater can become negatively buoyant and sink [18]. In this study, we will focus our attention on the long-term fate of buoyant microplastics in their clean state, subject to biofouling, which will eventually trigger sinking. Will the biofouled microplastics ever return to the sea surface?

It is extremely difficult to observe the trajectories of biofouled microplastics in the marine environment in order to assess their long-term distribution. The paucity of observations of biofouled microplastics in the ocean interior has encouraged the development of mathematical models predicting their trajectories. These models are currently the best tools we have for predicting the location of “high concentration hotspots”, either within the water column or on the seabed, and how they might evolve in space and time.

The mathematical models fall into two broad categories: deterministic and stochastic. Examples of the former type [19,20] describe the vertical motion of biofouled microplastics in the water column in the absence of ocean dynamics. Deterministic microplastic trajectory models, in which the concentration of microplastic per unit volume of seawater is treated as a continuous, differentiable field, have been coupled to ocean general circulation models [21,22]. Deterministic Lagrangian (i.e., particle tracking) models have been utilised to predict the motion of debris remaining at the sea surface [23] and for determining the three-dimensional motion of plastic in the water column [24,25]. In a stochastic modelling approach [26], the one-dimensional model [19] is driven by vertical velocities derived from a global NEMO ocean circulation experiment, with the addition of turbulent stochastic transport in the surface mixed layer parameterised using a Markov-0 random walk model.

The purpose of this study is to investigate the impact of time-varying insolation on the vertical motion of biofouled microplastic particles using the model of Kreczak et al. (2021) [20], which, in turn, was inspired by that of Kooi et al. (2017) [19]. For simplicity, Kreczak et al. (2021) [20] assume constant insolation, and they predict that the large-time behaviour of biofouled microplastic trajectories in a stratified ocean will exhibit subsurface oscillations bounded by the pycnocline and the base of the euphotic layer. Clearly, light is a key factor in the rate of growth of the biofilm, and the question arises of whether the large-time subsurface oscillations predicted by the former study will persist for unsteady insolation. We will demonstrate that rich oscillatory behaviour is exhibited by biofouled microplastics driven by unsteady insolation, especially at high latitudes where there are lengthy periods of little or no sunlight. Furthermore, the computationally efficient model developed in this study is well suited for coupling to three-dimensional ocean circulation models for a more comprehensive investigation of microplastic trajectories.

The plan of the paper is as follows. Section 2 reviews the dimensionless equations governing the behaviour of vertical oscillations of a biofouled microplastic particle in the presence of unsteady insolation. In contrast with Kreczak et al. (2021), the scalings used to derive the dimensionless governing equations are different in this study; a natural choice for the scaling of time comes from the prescribed insolation function, namely, one day. Section 3 proceeds to non-dimensionalise the equations governing the vertical motion of biofouled microplastic particles. The dimensionless equations are solved numerically, and the results are presented in Section 4, which is sub-divided for ease of presentation as follows. Section 4.1 examines the particle trajectory in a homogeneous ocean with unsteady insolation and compares the solutions with those forced by steady, annually averaged insolation. The particle dynamics are studied in a stratified ocean with a mixed layer depth deeper/shallower than the constant euphotic layer depth in Section 4.2. Section 4.3 investigates the effect of clean microplastic’s particle size on its dynamics in a stratified ocean. Finally, the study concludes with a summary and discussion of the results in Section 4.4.

2. Equations Governing the Vertical Motion of Biofouled Microplastics in the Ocean

The dimensional governing equations describing the vertical motion of biofouled microplastic particles in the absence of ocean dynamics are discussed in detail in Kreczak et al. (2021) [20]. For completeness, we briefly summarise these equations before proceeding to non-dimensionalise them in a new manner, suitable for unsteady insolation. With respect to a z-axis aligned with the local gravitational acceleration, the infinitely deep stratified ocean occupies $z<0$, with $z=0$ corresponding to the undisturbed sea surface. The settling velocity of a spherical biofouled microplastic particle is given by

$$\frac{dz}{dt}=\frac{Dg}{6\pi \nu a_B{\rho }_f\left(z\right)}[{\rho }_f\left(z\right)(V_p+N{V}_A)-V_p{\rho }_p-N{V}_A{\rho }_A],$$

(1)

where t is time, g is the gravitational acceleration, $\nu$ is the kinematic viscosity of seawater, and ${\rho }_f\left(z\right)$ is the ocean density profile, $a_B$ is the radius of the fouled microplastic particle, V and $\rho$ denote, respectively, the volume and density of the microplastic particle (denoted by subscript P) and a single algal cell (denoted by subscript A), N is the number of particles in the biofilm growing around the plastic particle, and $D=$ 86,400 (sd ${}^{-1}$) is a factor ensuring that the velocity timescale aligns with the algal population net growth rate, commonly expressed in the literature with units $d^{-1}$ [11,12]. For simplicity, the radius of the fouled particle $a_B$ (m) is assumed to be constant, and the biofilm is assumed to have grown uniformly over the particle’s surface. To date, there is no accepted universal expression for the settling velocity appropriate for microplastics of arbitrary shape. Employing (1) enables us to directly compare the fouled particle behaviour driven by unsteady insolation with that in [20], where insolation is steady.

Kreczak et al. (2021) [20] represented the ocean stratification with a step density profile in which the pycnocline thickness is infinitesimally thin. In this study, we use a more realistic open ocean density profile [27]:

$${\rho }_f\left(z\right)=\left\{\begin{array}{cc}{\rho }_0,\hfill & -h<z\le 0\hfill \\ {\rho }_0+\frac{{\rho }_0{N}^2}{\gamma g}[e^{-\gamma h}-e^{\gamma z}],\hfill & z\le -h\hfill \end{array}\right.$$

(2)

In (2), h is the uniform mixed layer depth, N is the constant buoyancy frequency, ${\gamma }^{-1}$ is the characteristic depth scale of the pycnocline, and ${\rho }_0$ is a constant surface density. As $z\to -\infty$, we assume that ${\rho }_f\to {\rho }_b$ a constant deep abyssal density. Then, we can rewrite (2) as ${\rho }_f\left(z\right)={\rho }_0+\mathsf{\Delta }f\left(z\right)$, where $\mathsf{\Delta }={\rho }_b-{\rho }_0$, and $0\le f\left(z\right)\le 1$ is a normalised function defined as

$$f\left(z\right)=\left\{\begin{array}{cc}0,\hfill & -h\le z\le 0,\hfill \\ 1-e^{\gamma h+\gamma z},\hfill & z<-h.\hfill \end{array}\right.$$

(3)

Substituting this expression for the oceanic density profile into (1) yields the following expression for the settling velocity of a spherical biofouled microplastic particle:

$$\frac{dz}{dt}=\frac{Dg({\rho }_A-[{\rho }_0+\mathsf{\Delta }f\left(z\right)])V_A}{6\pi \nu a_B[{\rho }_0+\mathsf{\Delta }f\left(z\right)]}\left(\frac{([{\rho }_0+\mathsf{\Delta }f\left(z\right)]-{\rho }_p)V_p)}{({\rho }_A-[{\rho }_0+\mathsf{\Delta }f\left(z\right)])V_A}-N\right).$$

(4)

Following Kooi et al. (2017) [19] and Kreczak et al. (2021) [20], algal biofilm is governed by a simple population model

$$\frac{dN}{dt}={\lambda }_A\left(z\right)+({\lambda }_G(z,t)-{\lambda }_D)N,$$

(5)

where ${\lambda }_A$ is the algal reattachment rate, ${\lambda }_G$ the growth rate, and ${\lambda }_D$ is the death rate, all with units $d^{-1}$. In contrast with Kreczak et al. (2021) [20], the growth rate, ${\lambda }_G(z,t)$, is a function of time as well as depth, reflecting the fact that photosynthesis is controlled by the unsteady light intensity in the water column. Light intensity is dependent on ocean depth due to light absorption and scattering by organic and inorganic particles in the water column, and it is modelled using the Beer–Lambert relation [28]:

$$I(z,t)=I_S\left(t\right)exp\left(ϵz\right).$$

(6)

In (6), $I(z,t)$ is the light intensity in the water column at a fixed latitude, $I_S\left(t\right)$ is the insolation at the sea surface, which varies on the diurnal and seasonal timescales as well as latitude, and ${ϵ}^{-1}$ is the e-folding decay scale of light intensity. We decompose the surface insolation as

$$I_S\left(t\right)=\widetilde{I_S}S\left(t\right),$$

(7)

where $\widetilde{I_S}$ is the insolation at the summer solstice, and $S\left(t\right)$ is a normalised function capturing the daily and seasonal variation of the insolation:

$$\mathsf{\Delta }\le S\left(t\right)\le 1,\mathsf{\Delta }={\tilde{I}}_W/{\tilde{I}}_S.$$

(8)

In (8), ${\tilde{I}}_W$ is the minimum insolation at the winter solstice.

Insolation as a function of latitude, $\varphi$, and time is calculated from the experimentally determined equation

$$I_S=1.353\times {0.7}^{{AM}^{0.678}},$$

(9)

where $AM$ is airmass, the length of the path which light takes through the atmosphere normalised by its shortest path length when the sun is directly overhead [29]. The airmass quantifies the reduction in the power of light as it passes through the atmosphere and is absorbed by air and dust. It is given by

$$AM=\frac{1}{cos\theta },$$

(10)

where $\theta$ is the solar zenith angle, defined such that when the sun is directly overhead, $\theta =0$. Following Honsberg and Bowden (2019) [29],

$$cos\theta =sin\varphi sin\delta +cos\varphi cos\delta cosh,$$

(11)

where $\delta$ is the declination angle of the Sun, which varies seasonally due to the tilt of the Earth on its axis of rotation and the rotation of the Earth around the Sun, and h is the hour angle. These two angles are calculated using the expressions

$$\delta =-23.45°\times cos\left(\frac{360}{365}\times (d+10)\right),$$

(12)

$$h=15°(t-12),$$

(13)

where d is the day of the year, defined such that $d=1$ on 1 January, and t is time in hours. In the Supporting Information (S1), examples of the time-dependent insolation at $55°$ N are presented, calculated using (9)–(13). Specifically, Figure S1 is a plot of the summer (21 June, day 172) and winter (21 December, day 355) solstice insolation at $55°$ N, and Figure S2 plots the daily maximum insolation over a year at this same latitude.

Light sensitive algal growth is modelled using a Monod-like system, capturing growth rate saturation with increasing light intensity, ${\lambda }_G(z,t)={\lambda }_G^{max}g(z,t)$ [15]:

$$g(z,t)=\frac{I(z,t)}{I(z,t)+{\lambda }_G^{max}/\alpha {\tilde{I}}_S}=\frac{S\left(t\right)}{S\left(t\right)+exp(-ϵz){\lambda }_G^{max}/\alpha {\tilde{I}}_S},$$

(14)

using (6) and (7). In (14), ${\lambda }_G^{max}$ is the maximum growth rate (with units $d^{-1}$), and $\alpha$ ($d^{-1}\mu moles{m}^{-2}{s}^{-1}$) governs how the growth rate changes with light intensity.

The attachment of algae onto the biofouled microplastic particle from a pre-existing distribution in the ocean’s water column is represented by the equation:

$${\lambda }_A\left(z\right)={\beta }_A{A}_{max}c\left(z\right),$$

(15)

where ${\beta }_A$ is the rate of collision between the microplastic particle and a single algal cell, $A_{max}$ is the maximum background algae concentration, and $0\le c\left(z\right)\le 1$ is a dimensionless normalised function describing how algae concentration changes with depth.

3. Derivation of the Dimensionless Governing Equations for the Vertical Motion of Biofouled Microplastic Particles in the Ocean

We now proceed to non-dimensionalise the governing equations developed in the previous section, thereby revealing the number of independent parameters that govern the behaviour of biofouled microplastic trajectories. Kreczak et al. (2021) [20] carried out the equivalent exercise for this model, albeit with constant insolation, with a timescale corresponding to the time taken for a clean microplastic particle to traverse the euphotic layer. For time-dependent insolation, it is appropriate to use a timescale associated with this field, namely, 24 h. At the ocean surface, $f\left(z\right)=0$, and using (5) we can find the number of algal cells $N^{*}$ such that the biofilm makes the particle neutrally buoyant ($dz/dt=0$):

$$N^{*}=\frac{V_p({\rho }_0-{\rho }_p)}{V_A({\rho }_A-{\rho }_0)}.$$

(16)

We now introduce the following dimensionless variables, denoted by tilde; $N=N^{*}\tilde{N}$, $t=T_s\tilde{t}$, where $T_s$ is the timescale, taken to be equal to 1 day, and $z=\tilde{z}/ϵ$. The dimensionless equation for the rise speed becomes

$$\frac{d\tilde{z}}{d\tilde{t}}=\frac{\delta }{1+R_{1}f\left(\tilde{z}\right)}[1+R_{2}f\left(\tilde{z}\right)-(1-R_{3}f\left(\tilde{z}\right))\tilde{N}],$$

(17)

where

$$\delta =T_sϵW^{*},R_1=\frac{\mathsf{\Delta }}{{\rho }_0},R_2=\frac{\mathsf{\Delta }}{{\rho }_0-{\rho }_p},R_3=\frac{\mathsf{\Delta }}{{\rho }_A-{\rho }_0},$$

(18)

and $W^{*}$ is the rise speed of a clean microplastic particle, namely,

$$W^{*}=\frac{Dg{V}_p({\rho }_0-{\rho }_p)}{6\pi \nu a_B{\rho }_0}.$$

(19)

In (17),

$$f\left(\tilde{z}\right)=\left\{\begin{array}{cc}0\hfill & -H\le \tilde{z}\le 0\hfill \\ 1-e^{\mathsf{\Gamma }H+\mathsf{\Gamma }\tilde{z}}\hfill & \tilde{z}<-H\hfill \end{array}\right.$$

(20)

where $H=ϵh$, and $\mathsf{\Gamma }=\gamma /ϵ$.

The dimensionless biofilm growth equation is given by

$$\frac{d\tilde{N}}{d\tilde{t}}=\mathsf{\Omega }c\left(\tilde{z}\right)+\mathsf{\Phi }\left[\mathsf{\Lambda }\frac{S\left(\tilde{t}\right)}{S\left(\tilde{t}\right)+\mathsf{\Psi }exp(-\tilde{z})}-1\right]\tilde{N},$$

(21)

where

$$\mathsf{\Omega }=\frac{T_s{\beta }_A{A}_{max}}{N^{*}},\mathsf{\Phi }=T_s{\lambda }_D,\mathsf{\Lambda }=\frac{{\lambda }_G^{max}}{{\lambda }_D},\mathsf{\Psi }=\frac{{\lambda }_G^{max}}{\alpha {\tilde{I}}_s}$$

(22)

are dimensionless parameters.

The base of the euphotic layer is defined to be the location where 99% of the incident solar radiation at the surface of the ocean has been extinguished. When insolation is steady, the depth of this layer, $z_E$, is given by:

$$z_E=ln\left(0.01\right)/ϵ.$$

(23)

For unsteady insolation, given by (7) and (8), this layer depth is modified as follows. The unsteady euphotic layer depth, $z_E\left(t\right)$, is the solution of the equation

$$0.01{\tilde{I}}_S=I_S\left(t\right)exp\left(ϵz_E\left(t\right)\right).$$

(24)

Using (7), we can readily show that this depth, non-dimensionalised by $1/ϵ$, is given by

$${\tilde{z}}_E=ln\left[0.01\right]-ln\left[S\right(t\left)\right].$$

(25)

At higher latitudes, insolation partly or totally vanishes in winter, corresponding to $S\left(t\right)$ vanishing. We thus modify (25) to take this into account. Let $\overline{S}\left(d\right)$ be the daily average value of $S\left(t\right)$. Then, we modify (25), the unsteady non-dimensional euphotic layer depth, as follows:

$${\tilde{z}}_E=\left\{\begin{array}{cc}\overline{S}\left(d\right),\hfill & \overline{S}\left(d\right)\le 0.01\hfill \\ ln\left[0.01\right]+ln[1/\overline{S}\left(d\right)],\hfill & \overline{S}\left(d\right)>0.01\hfill \end{array}\right.$$

(26)

4. Results

Before exploring the impact of unsteady insolation on the vertical motion of a bio-fouled microplastic particle,

Table 1. Summary of the biofouled microplastic trajectory behaviour driven by steady insolation in the model of Kreczak et al. (2021) [20].

Description of the Ocean	Large-Time Trajectory Behaviour
Unstratified ocean.	Particle oscillates from the surface to a fixed depth. When the size of the particle is decreased, the oscillation period increases.
Stratified ocean.	Particle initially sinks, then becomes trapped in a subsurface layer. Oscillations are bounded between the base of the euphotic layer and the pycnocline.
Background algal field with a fixed attachment rate onto the particle in either a unstratified/stratified ocean.	In a homogeneous ocean, the particle exhibits damped oscillates around the base of the euphotic zone layer. The damping rate increases as the clean particle size decreases. In a stratified ocean, the particle exhibits subsurface oscillations bounded between the pycnocline and the base of the euphotic layer.

briefly summarises the particle behaviour discussed in Kreczak et al. (2021) [20] when the insolation is constant.

We will contrast the above behaviour with the results discussed below for unsteady insolation.

4.1. Homogeneous Ocean with Background Algal Field Absent

We begin with a highly idealised example in which the ocean is unstratified, (i.e., $\mathsf{\Delta }=0$), corresponding to $R_1=R_2=R_3=0$, in (18). Upon dropping the tilde superfix denoting a non-dimensional quantity, the governing equations for the biofouled microplastic trajectory become:

$$\frac{dz}{dt}=\delta [1-N],$$

(27)

$$\frac{dN}{dt}=\mathsf{\Phi }(\mathsf{\Lambda }g(z,t)-1)N,$$

(28)

upon setting $c\left(z\right)\equiv 0$, corresponding to no background algal field.

In analysing the biofouled particle trajectory, a useful diagnostic is the “compensation depth” where production and mortality of the algae are in balance. At this depth, $d\tilde{N}/d\tilde{t}=0$, and (28) implies that the compensation depth, $z_C\left(t\right)$ satisfies

$$\mathsf{\Lambda }\frac{S\left(t\right)}{S\left(t\right)+\mathsf{\Psi }{e}^{-z}}-1=0.$$

(29)

From (29), the non-dimensional compensation depth, $z_C$, is given by

$$z_C=ln\left[\frac{\mathsf{\Psi }}{(\mathsf{\Lambda }-1)S\left(t\right)}\right].$$

(30)

If a background algal field is present (i.e., $c\left(z\right)$ is not identical to zero), we postulate that no modification of (30) is required because $c\left(z\right)$ is not associated with “new photosynthetic algae production”.

The governing Equations (27) and (28), and more generally (17) and (21), are solved numerically using Python. A Runge–Kutta–Fehlberg adaptive time step scheme is employed, with a minimum time step of $dt={10}^{-3}$. The following dimensionless parameter values $\mathsf{\Phi }=0.5$, $\mathsf{\Lambda }=4$, and $\delta =350$ are used, corresponding to an oscillation period of the clean particle of approximately one day. The particle trajectories, forced by time-varying insolation, are calculated at latitudes $30\mathrm{N}$, $55\mathrm{N}$ and $70\mathrm{N}$. Figure 1 shows the particle trajectory (blue curve) at $30\mathrm{N}$ over a 5-day interval, where the red line is the compensation depth and the green line is the constant euphotic zone depth. We note that, over such a short time interval, the compensation depth is approximately constant.

Figure 2 depicts the 5-year trajectory of a microplastic particle at a latitude of 30N, where the seasonal behaviour of the compensation depth is apparent. Figure 1 reveals that the particle remains at the sea surface for a short period until it acquires sufficient biofilm to become negatively buoyant. It then performs vertical oscillations in the upper water column with a daily period. Once the particle has sunk beneath the euphotic layer, defouling starts to occur, due to a lack of sunlight. In addition, the periods of darkness, when no photosynthesis takes place, also contribute to the decay of the biofilm. As the biofilm dies off, the particle trajectory starts to slow down. The particle reaches a depth where its buoyancy is neutral, and thereafter it begins to rise. Upon reaching the surface, the particle remains there until it acquires sufficient biofilm to become negatively buoyant again. In summer, the particle sinks to a deeper depth than in winter, due to higher amounts of insolation being present during the summer months, supporting greater biofilm mass than in the winter, as seen in Figure 2.

Figure 2 also shows that the microplastic particle trajectory mirrors the time-varying compensation depth, being shallower in winter and deeper in summer, which is to be expected as the latter is calculated from daily average insolation and represents the depth at which the algal growth is balanced by the mortality rate.

Figure 3 plots the 5-year trajectory of a microplastic particle at a latitude of 55N. The particle trajectory displays the same qualitative behaviour as a particle at 30N. Initially, the particle begins sinking, but to a shallower depth than the particle observed at 30N due to there being less insolation, and thus less biofilm growth, during winter at 55N. We observe that the particle exhibits oscillatory behaviour, resting at the surface for a period of days during the winter months when biofilm growth is at a minimum. As the year progresses from winter to summer, the particle begins to sink deeper, due to an increase in daily integrated insolation supporting greater biofilm mass. Once again, the particle trajectory over each year mirrors the behaviour of the compensation depth.

Figure 4 plots the 5-year trajectory of a microplastic particle at a latitude of 70N. Now, the particle displays aperiodic behaviour, unlike the trajectories in Figure 2 and Figure 3.

At this high latitude, the particle floats at the ocean surface for a prolonged period of time during the winter and over the early spring, when there is little or no insolation present and hence no biofilm production to make the particle negatively buoyant.

Once there is sufficient insolation for biofilm growth to occur, the particle sinks suddenly. The cause of this sudden drop is revealed in Figure 5, where we observe that the algal population increases until it reaches a critical level at around day 105, when the particle becomes negatively buoyant. During the summer, there is a period of time during which the microplastic does not resurface. This is because the sun does not set, and so there is always a certain amount of penetrative radiation present in the upper ocean water column. When the biofilm begins to die off below the base of the euphotic layer, the particle begins to rise, and once it has re-entered the euphotic layer there is sufficient summer insolation for photosynthesis to begin biofilm growth once more, thus preventing the particle resurfacing. Unlike Figure 2 and Figure 3, the particle trajectory does not mimic the behaviour of the compensation depth. However, it is observed that at times when the compensation depth is at a maximum, the particle begins to sink after long periods of rest at the ocean surface. At these times, the biofilm production is near its maximum, and in the water column above this depth, biofilm production exceeds mortality.

In summary, it can be seen that there are significant variations between the particle trajectories in Figure 2, Figure 3 and Figure 4. Microplastic trajectories at latitudes 30N and 55N exhibit time-periodic behaviour, whereas the trajectory of a microplastic particle at 70N displays aperiodic oscillations. At both 55N and 70N, there are prolonged periods of time in winter where the microplastic rests at the ocean surface, whereas this behaviour is not seen at 30N. It is also shown that there is no year-on-year accumulating effect of insolation on biofilm growth. Each year, the trajectory of the microplastic particle is the same, with the only exception being at a latitude of 70N for the first year, due to the initial conditions for the numerical integration of (27) and (28). The maximum depth achieved by the microplastic at these three latitudes is approximately the same, namely $\tilde{z}=14$, or around three times the euphotic zone depth.

We now contrast the particle trajectories plotted in Figure 3 and Figure 4 with those forced by constant insolation considered by Kreczak et al. (2021) [20]. At latitudes 55N and 70N, the annual average insolation is used to force the microplastic particle. Figure 6 and Figure 7 show the corresponding particle trajectory over a 30-day time interval at latitudes 55N and 70N, respectively. The 30-day time interval used in these plots clearly resolves the particle oscillation.

In Figure 6 and Figure 7, the maximum depth reached by the particle is significantly shallower compared with the summer maximum depth acquired by the particle when insolation varies with time (see Figure 3 and Figure 4). This is because the summer biofilm production in Figure 3 and Figure 4 exceeds the constant biofilm production rate in Figure 6 and Figure 7. In Figure 6 and Figure 7, the maximum depth attained by the particle is greater at a latitude of 55N compared with 70N, due to the slightly greater average annual insolation at the lower latitude. Naturally, this behaviour is also reflected in the compensation depth. Note also that the microplastic particle does not float at the ocean surface for a period of time during winter when insolation is constant, due to biofilm growth occurring throughout the year when the particle enters the euphotic zone. With unsteady insolation, there is no biofilm growth in the absence of sunlight, irrespective of whether the particle is in the euphotic layer or not.

4.2. Stratified Ocean with Background Algal Field Absent

In the previous subsection, we assumed that there was no stratification in the ocean by setting $\mathsf{\Delta }=0$. In this section, we introduce a stratified ocean with density profile given by $f\left(\tilde{z}\right)$ in (20). Once again, we neglect algal attachment ($\mathsf{\Omega }=0$) onto the microplastic from a pre-existing field in the upper ocean.

In Kreczak et al. (2021) [20], a step stratification was adopted, where the pycnocline depth was infinitesimally thin. In this study, the stratification (20) has a more realistic finite-depth pycnocline, which scales as ${\gamma }^{-1}$. Recall from (20) that in dimensionless form this becomes $\mathsf{\Gamma }=\gamma /ϵ$. We now examine the microplastic trajectories as a function of the pycnocline depth for two cases: (i) mixed layer deeper than the euphotic layer; (ii) mixed layer shallower than the euphotic layer. In the following numerical results, $R_1=3/1025$, $R_2=1/100$, and $R_3=1/120$. A euphotic layer depth of 75m is fixed throughout the results using a stratified ocean, corresponding to $ϵ=-ln\left(0.01\right)/75$. For comparison with Kreczak et al. (2021) [20], we will first use constant insolation fixed at the annual average value at 55N, followed by results for time-varying insolation at this latitude.

Figure 8 plots the particle trajectory when insolation is constant, with a pycnocline thickness of 50 m and a mixed layer depth below the euphotic zone. Unlike the behaviour observed in Figure 3, the particle no longer reaches the ocean surface and instead oscillates in a subsurface duct, bounded below by the base of the mixed layer, similar to the behaviour observed in Kreczak et al. (2021) [20]. This subsurface trapping behaviour is due to the density jump in the pycnocline, displayed in Figure 9. Once the particle enters the pycnocline, its velocity decreases due to its relative buoyancy decreasing.

Fewer algal cells need to be removed in order to attain neutral buoyancy in the pycnocline. The particle begins to rise, crossing back into the mixed layer, and now weighing more than it did on its descent. When the particle enters the euphotic zone, biofilm growth restarts, and now, less time is required for the particle to attain neutral buoyancy. Thus, the particle does not resurface. This process repeats until the particle becomes trapped above the pycnocline.

Figure 10 plots the particle trajectory for a pycnocline thickness of 500 m, for comparison with Figure 8. There are some notable differences between the two figures. For instance, the adjustment time for the particle in Figure 10 to exhibit trapped subsurface oscillations of fixed amplitude, bounded below by the base of the mixed layer, is longer compared to Figure 8. With the smaller pycnocline density gradient in Figure 11, compared with Figure 9, the particle reaches neutral buoyancy at a greater depth, thus extending the damping timescale.

We now turn to the case when the mixed layer depth is shallower than the euphotic layer and examine the particle trajectories forced by constant, annually averaged insolation at 55N, as shown in Figure 12 and Figure 13.

Unlike Figure 8 and Figure 10, the particle’s oscillation is no longer bounded below by the mixed layer depth, and occurs within the pycnocline. The particle now displays rapidly damped oscillations about the average compensation depth, whereas in Kreczak et al. (2021) [20], the particle oscillated about the base of the euphotic layer. Indeed, the large-time behaviour of the particle is one of equilibrium, as predicted by (17). From (17), this equilibrium depth within the pycnocline is given in dimensionless form as

$$\tilde{z}=-H+{\mathsf{\Gamma }}^{-1}ln\left(\frac{1+R_2+\tilde{N}(R_3-1)}{R_2+\tilde{N}{R}_3}\right).$$

(31)

Again, it can be seen that for a thicker pycnocline, Figure 13, the adjustment time to reach equilibrium is longer compared with the thinner pycnocline example shown in Figure 12.

We now introduce time-varying insolation, to assess its impact on the particle trajectory in a stratified ocean. For brevity, we set the pycnocline thickness to 500 m in all the results discussed in the remaining part of this section. Figure 14 and Figure 15 plot the particle trajectory at latitudes 55N and 70N, respectively, when the mixed layer is deeper than the euphotic layer.

Comparing Figure 10 and Figure 14, we see that the particle is no longer performing trapped subsurface oscillations when forced by unsteady insolation. However, at 70N (see Figure 15), the particle displays three distinct dynamical regimes: (i) floating at the surface throughout the winter, (ii) spring and autumn large-amplitude surface-to-pycnocline oscillations; (iii) subsurface trapped oscillations about the compensation depth. In the summer, when the insolation is quasi-steady and at maximum, there is sufficient biofilm attached to the microplastic to prevent surfacing. At its deepest point within the pynocline during the summer, biofilm remains on the particle. Thus, on its ascent into the euphotic zone, it does not reach the surface because new, vigorous biofilm has grown once the particle is above the compensation depth. If the insolation were to remain fixed at its summer value throughout the following autumn and winter, the particle would perform damped subsurface oscillations trapped around the compensation depth. Of course, this does not occur in Figure 15 because as summer gives way to autumn and biofilm production tails off, the particle reverts to oscillations that reach the sea surface.

We now consider the trajectory behaviour throughout the year at 70N when the mixed layer is shallower than the euphotic layer, more typical of that in summer. For brevity, we omit the plot for 55N using this stratification, because it is very similar to that of Figure 14.

Figure 16 plots the annual trajectory of a particle forced by time-varying insolation at 70N and with a shallow mixed layer. Comparing Figure 15 and Figure 16, we see that the shallower mixed layer leads to more strongly damped subsurface oscillations about the compensation depth. Once again, the particle in Figure 16 exhibits three regimes of trajectory behaviour throughout the year.

4.3. Varying the Clean Microplastic Particle Size

In all of the above results, the clean microplastic particle size was set so that its oscillation period was equal to one day. In this section, we now examine the particle trajectory behaviour at 70N for smaller clean microplastic particles which have slower rise velocity [4]. The following results include stratification with the mixed layer depth greater than the euphotic layer depth, a pycnocline thickness of 500 m and time-varying insolation. For simplicity, we again neglect the background algal field $\mathsf{\Omega }=0$. We keep the growth and mortality rate ratio, $\mathsf{\Lambda }=4$, and the algal population and particle motion timescale ratio, $\mathsf{\Phi }=0.5$, both fixed.

Figure 17 plots the trajectory of a particle that is half the original particle size used in Section 4.1 and Section 4.2 ($\delta =87.5$). There are three distinct timescales at play in Figure 17 that result in a highly non-linear oscillation. These are the diurnal and seasonal timescales of the insolation and the clean particle oscillation period. Qualitatively, the particle trajectories in Figure 15 and Figure 17 are identical, albeit with significant fluctuations in the amplitude seen in the former plot. During the summer, the slower rise speed for the biofouled particle in Figure 17 leads to larger amounts of biofilm when it is in the euphotic zone compared with Figure 15, which leads to a longer period of subsurface oscillations about the compensation depth.

Decreasing the size of the clean particle to one tenth of the original size, $\delta =3.5$, used in Section 4.1 and Section 4.2 drastically alters the trajectory behaviour, as shown in Figure 18 and Figure 19. With the onset of spring biofilm growth, the particle sinks and exhibits subsurface oscillations with an 8-day period (see Figure 19) about the compensation depth. Now, the length of time during which the particle exhibits subsurface oscillations is approximately 140 days, following the trend noted in Figure 17. In the autumn, the oscillations are confined to the mixed layer and reach the surface. As expected, during the winter months the oscillations cease entirely and the particle floats at the sea surface.

4.4. Discussion

In this study, the impact of time-dependent insolation on the vertical trajectory of a biofouled microplastic particle in the ocean is explored using the model of Kreczak et al. (2021) [20], which in turn is a simplified version of the model by Kooi et al. (2017) [19]. The latter study introduced the first one-dimensional deterministic model for predicting the vertical motion of a biofouled microplastic particle in the ocean. There are two notable departures from the model formulation of Kreczak et al. (2021) [20] employed in this study. First, the daily timescale is used in formulating the dimensionless governing equations for the particle trajectory, rather than the time taken for a clean (buoyant) microplastic particle to traverse the euphotic layer. Insolation, which drives the particle motion by supporting biofilm production, varies on the daily and seasonal timescales, and it is natural to choose one of these to non-dimensionalise time. Second, a more realistic density profile with a non-zero finite depth pycnocline is used in this study, whereas Kreczak et al. (2021) [20] employed a step density profile. The dimensionless governing equations for the microplastic particle motion reveal that seven dimensionless parameters govern its behaviour when neglecting a pre-existing algal field in the upper ocean. Fortunately, once the stratification and densities of the clean plastic and the biofilm are prescribed, the number of dimensionless parameters reduces to four.

A new diagnostic variable is introduced in this study, namely, the compensation depth, where biofilm growth is balanced by its mortality. Because growth is a function of irradiance, the compensation depth has the same time-dependent behaviour as the insolation. We demonstrate that, under certain circumstances, microplastic particles oscillate about the compensation depth, as discussed below.

Initially, the study examines the microplastic trajectories driven by unsteady insolation in a homogeneous ocean. In this case, the particle exhibits vertical oscillations that reach the sea surface, as in Kreczak et al. (2021) [20], although with amplitude that varies seasonally. In the summer, the amplitude is at a maximum, reflecting the fact that biofilm production is at a maximum. Indeed, the seasonal behaviour of the oscillation amplitude is reflected in the seasonal behaviour of the compensation depth. Moving poleward, the particle spends an increasing amount of time floating at the sea surface, and the oscillations become aperiodic. For example, at 30N, the particle spends several hours at the sea surface during the night when biofilm production ceases. Throughout the year, at this latitude the oscillations persist. Moving to 55N, the oscillations are absent entirely during a period centred around the solstice, because there is insufficient biofilm production to make the particle negatively buoyant. This trend continues at 70N, where during the winter, insolation is absent altogether for several weeks, and the particle floats on the sea surface. When the oscillations commence at this high latitude, they are highly modulated, exhibiting maximum amplitude in the spring and fall. In both of these seasons, the particle reaches the sea surface during each oscillation cycle. However, in a neighbourhood centred around the summer solstice, the particle oscillations briefly become subsurface in response to the large quasi-steady biofilm production.

For comparison with Kreczak et al. (2021) [20], results are also presented for a particle driven by constant, annual average insolation at 55N and 70N in a homogeneous ocean. In both cases, the particle exhibits periodic oscillations that reach the sea surface during every cycle, as noted in Kreczak et al. (2021) [20]. The depth of these oscillations decreases with increasing latitude. Also of note, for a given latitude, the maximum summer amplitude of the oscillations driven by unsteady insolation is greater than that driven by the annual average insolation.

Particle trajectories are calculated in a stratified ocean when: (i) the pycnocline thickness varies, (ii) the depth of the mixed layer is deeper/shallower than the fixed euphotic layer depth. In contrast with this study, the step-density profile used by Kreczak et al. (2021) [6] did not permit an investigation of how the pycnocline density gradient impacts on the particle motion. Forcing the particle motion with annually averaged insolation at 55N leads to large-time, constant-amplitude subsurface oscillations bounded below by the base of the mixed layer, irrespective of the pycnocline thickness. With a thick pycnocline, corresponding to a smaller vertical density gradient, the particle penetrates the pycnocline during the initial adjustment phase. When the mixed layer is shallower than the euphotic layer, significant biofilm growth can take place in the pycnocline. In this case, forcing the particle with annually averaged insolation at 55N leads to highly damped oscillations around the compensation depth, and for large times, the particle is trapped at this depth.

Particle trajectories are calculated using unsteady insolation at 55N and 70N for a fixed pycnocline thickness of 500 m and with a mixed layer deeper/shallower than the euphotic layer. At 55N, the particle exhibits surface-to-depth oscillations with amplitude varying seasonally. However, at 70N, the particle dynamics displays three regimes: (i) floating at the sea surface in winter, (ii) surface-to-depth oscillations in spring/fall with time-varying amplitude, (iii) subsurface oscillations around the compensation depth. The time interval during which subsurface oscillations persist increases/decreases if the mixed layer is shallower/deeper than the euphotic layer. Subsurface oscillation amplitude increases/decreases when the mixed layer is deeper/shallower than the euphotic layer.

Finally, the study investigates how decreasing the clean particle size impacts its trajectory when forced with time-varying insolation in a stratified fluid. As the clean particle size decreases, its rise speed decreases, leading to longer-period oscillations in the presence of biofilm attachment. This is borne out in the results. In addition, as the clean particle size decreases, the particle transitions from exhibiting the three dynamical regimes noted above to one, in which it either floats at the sea surface or performs subsurface oscillations about the compensation depth with seasonally varying amplitude.

In this study, we have not considered the impact of a pre-existing algal field on the microplastic trajectories (i.e., allowing $\mathsf{\Omega }$ to be non-zero in Equation (21)), because this warrants a study in its own right. The specification of the vertical distribution of this field, captured by $c\left(z\right)$ in (15), is highly uncertain and depends on location and season. Preliminary investigations indicate that earlier sinking of the biofouled microplastic particle to a greater depth is a consequence of a non-zero $c\left(z\right)$. We have also assumed that the stratification is steady throughout the year, again to aid understanding of the dynamics of the biofouled particle forced by unsteady insolation. Local wind and buoyancy forcing gives rise to an upper mixed layer depth that varies on the diurnal and seasonal timescales [30,31]. Coupling a one-dimensional mixed layer model to the model described in this study to assess the impact of time-dependent stratification on the particle trajectory would be worthwhile.

The clean microplastic density has not been altered in this study, with the focus, instead, on time-varying insolation as a function of latitude driving the variable buoyancy of the biofouled particle. Increasing the clean particle density from below towards the density of the ocean surface mixed layer will lead to smaller-amplitude, shorter-period oscillations. Conversely, decreasing the clean particle density will lengthen the oscillation period and increase the oscillation amplitude. Clearly, these predictions are dependent on the rates of growth/decay of the biofilm, water column stratification and the depth of the euphotic layer.

We have also neglected the impact of a background ocean velocity field on the particle dynamics. Biofouled microplastic particles are expected to be strongly coupled to three-dimensional upper ocean turbulence. Kreczak et al. (2023) [32] have addressed one aspect of this complex problem by examining the statistics of a large ensemble of biofouled microplastic particles in the presence of a Taylor–Green vortex flow. This latter flow is a prototypical representation of a Langmuir cell. The key question in such studies is whether the vortical cell(s) delay the descent of the microplastic particles, and if so, by how long. A worthwhile future study would be to couple the microplastic trajectory model presented in this paper to a numerical model for upper ocean three-dimensional turbulence in order to address the question of delayed particle fall-out.