Physics kernels#
The parcels Lagrangian framework is a tool for advecting virtual particles that are assumed to be spherical in shape. It works by numerically integrating the velocity fields from a hydrodynamic model while including any additional \textit{behaviour} of the particle. Mathematically, particle trajectories are computed by solving the following equation:
where \(\mathbf{x}(t)\) describes the particle position at time \(t\), \(\mathbf{v} = (u,v,w)\) is the hydrodynamic model velocity field, and \(\mathbf{B}(\mathbf{x}(t),t)\) describes any displacements to the particle position caused by additional behaviour the particle exhibits or experiences. When performing a plastic dispersal simulation with plasticparcels, users have the explicit option of choosing which additional behaviour to include. Examples of these additional behaviours are described below.
Numerically, we solve the above equation using a time-stepping approach, where we compute the displacements in the particle position as
and updating the particle position at each timestep. For simplicity, by default we use the fourth-order Runge-Kutta scheme of parcels to solve the advection of the particle from the hydrodynamic model velocity field \(\mathbf{v}\), and an Euler-forward scheme for all other additional behaviours realised in \(\mathbf{B}\).
Stokes drift#
An important process that affects plastic particle dispersal in the upper ocean is the Stokes drift, whereby a particle subjected to a surface wave will experience a net displacement in the direction of wave propagation. We include a kernel to parameterise the effect of Stokes drift on a particle, based on the Phillips spectrum approximation developed in @Breivik2016. Specifically, we model this additional behaviour as \(\mathbf{B}_{\text{Stokes}}\), where the change in the particle position is described by
Here, \(z\) is the depth of the particle, \(\mathbf{v}_{\text{Stokes}}(\mathbf{x}_{z=0}(t),t)\) is the surface Stokes drift velocity, \(\beta=1\) (as we assume a Phillips spectrum), and erfc is the complementary error function. The peak wave number \(k_p\) is computed as \(k_p = \omega_{p}^2/9.81\), where \(\omega_p\) is the peak wave frequency \(\omega_p = 2 \pi / T_p\), using the peak wave period \(T_p = T_p(\mathbf{x}_{z=0}(t),t)\).
Our particular implementation of the Stokes drift kernel requires a surface Stokes velocity field \(\mathbf{v}_{\text{Stokes}}(\mathbf{x}_{z=0}(t),t)\), as well as a peak wave period field \(T_p(\mathbf{x}_{z=0}(t),t)\). Earlier versions of this kernel have been used in @Onink2022.
Wind-induced drift / Leeway#
Plastic particles at the ocean surface that are not completely submersed will experience a force from the relative wind due to a wind drag, leading to a wind-induced drift. This wind-induced drift of the particle is called leeway @Allen1999, which can be decomposed into a downwind component (in the direction of the wind), and a crosswind component (which is typically non-zero for asymmetric objects). As we assume that each plastic particle is spherical, we can ignore the crosswind component of leeway, and only consider the downwind component of leeway. The downwind component follows an almost linear relationship with the relative 10m wind speed @Allen2005, so we model the leeway as
where \(\mathbf{v}_{\text{Wind}}\) is the wind velocity 10m above sea level, and \(c\) is the leeway rate (a percentage of wind speed, which we refer to in the code as the windage coefficient). Ignoring all additional behaviour of the particle, then \(\mathbf{v}_{\text{Wind}} - \mathbf{v}\) is the relative wind acting on the particle.
A version of this kernel has been used in @Manral2024.
Biofouling#
Plastic particles in the ocean can be a hotbed for the accumulation and growth of organisms, known as biofouling. The formation of a biofilm on the surface of a plastic particle can result in a density change, affecting the buoyancy of the particle. An initially buoyant particle may become negatively buoyant, and sink or settle, depending on the surrounding seawater density.
We model the biofouling of a plastic particle following the approach of @Kooi2017, where the settling velocity of a particle is computed from the relative density difference of the plastic particle and the surrounding seawater. Here, we assume the biofilm growth (and decay) is primarily microbial algae, and is distributed homogeneously over the particle surface. The density of the biofouled plastic particle depends on the radius and density of the particle, and the thickness and density of the algal biofilm. The primary component of the biofouling kernel is modelling the change in the number of attached algae (denoted by \(A\)) on the surface of the plastic particle. As in @Kooi2017, we model the attached algal growth as
The first term models growth of algae due to collisions of the particle with algae in the surrounding seawater, where \(A_A\) is the ambient algal amount, \(\beta_A\) is the encounter kernel rate, \(\theta_{\text{Plastic}}\) is the surface area of the plastic particle. The second term models the growth of the biofilm, where the growth term \(\mu_A\) is computed from the total productivity provided by model output. The third and fourth terms model the (grazing) mortality and respiration of the biofilm respectively. As in @Kooi2017, we use constant mortality \(m_A\) and respiration \(R_{20}\) rates, with a temperature dependent term \(\big(Q_{10}^{(T-20)/10}\big)\) included in the respiration component (see @Kooi2017 for more details).
As described above, the modelled attached algal growth drives a change in the settling velocity of the biofouled particle, \(\mathbf{v}_{\text{Biofouling}}\). Hence, we model the additional behaviour of the particle due to biofouling as
This kernel has been used in various forms in @Lobelle2021, @Fischer2022, and @Kaandorp2023.
Vertical mixing#
An important process that is unresolved in even high-resolution ocean models is wind-driven turbulent mixing, which occurs at scales far smaller than a typical model ocean grid cell. In the vertical direction, this turbulent mixing can distribute even positively buoyant plastic particles throughout the mixed layer. To model this process, we take the approach of @Onink2022, by employing a Markov-0 styled stochastic parameterisation.
Denote by \(K_z = K_z(\mathbf{x}(t))\) the vertical diffusion coefficient profile based on a \(K\)-profile parameterisation (KPP) model @Large1994. Then the displacement of a particle (in the vertical direction) with a settling velocity \(w\) can be modelled as an SDE @Grawe2012,
where \(\text{d}W(t)\) is a Wiener noise increment with zero mean and a variance of \(\text{d}t\). In our case, the displacement due to the settling velocity of a particle is already accounted for in the biofouling kernel, hence we only model the stochastic term (by setting \(w=0\)). To numerically solve this equation, we use the stochastic generalisation of the Euler-forward scheme, called the Euler-Maruyama scheme @Maruyama1955.
A version of this kernel has been used in @Onink2022.
Sea-ice capture#
A sea-ice capture kernel is currently under development and will be released soon.