the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A comparison of Eulerian and Lagrangian methods for vertical particle transport in the water column
Ruben Kristiansen
Raymond Nepstad
Erik van Sebille
Andy M. Booth
A common task in oceanography is to model the vertical movement of particles such as microplastics, nanoparticles, mineral particles, gas bubbles, oil droplets, fish eggs, plankton, or algae. In some cases, the distribution of the vertical rise or settling velocities of the particles in question can span a wide range, covering several orders of magnitude, often due to a broad particle size distribution or differences in density. This requires numerical methods that are able to adequately resolve a wide and possibly multimodal velocity distribution.
Lagrangian particle methods are commonly used for these applications. A strength of such methods is that each particle can have its own rise or settling speed, which makes it easy to achieve a good representation of a continuous distribution of speeds. An alternative approach is to use Eulerian methods, where the partial differential equations describing the transport problem are solved directly with numerical methods. In Eulerian methods, different rise or settling speeds must be represented as discrete classes, and in practice, only a limited number of classes can be included.
Here, we consider three different examples of applications for a water column model: positively buoyant fish eggs, a mixture of positively and negatively buoyant microplastics, and positively buoyant oil droplets being entrained by waves. For each of the three cases, we formulate a model for the vertical transport based on the advection–diffusion equation with suitable boundary conditions and, in one case, a reaction term. We give a detailed description of an Eulerian and a Lagrangian implementation of these models, and we demonstrate that they give equivalent results for selected example cases. We also pay special attention to the convergence of the model results with an increasing number of classes in the Eulerian scheme and with the number of particles in the Lagrangian scheme. For the Lagrangian scheme, we see the $\mathrm{1}/\sqrt{{N}_{p}}$ convergence, as expected for a Monte Carlo method, while for the Eulerian implementation, we see a secondorder ($\mathrm{1}/{N}_{k}^{\mathrm{2}}$) convergence with the number of classes.
 Article
(3677 KB)  Fulltext XML
 BibTeX
 EndNote
Studying the vertical transport of positively, negatively, or neutrally buoyant particles is a common task in oceanography. Examples include both anthropogenic and naturally occurring particles, such as microplastics, mineral particles, nanoparticles, aggregates, gas bubbles, oil droplets, fish eggs, or even particles with active swimming behaviour such as zooplankton. These particles may show a range of different behaviours, including rising, sinking, and interacting with the ocean surface and seafloor in different ways.
Vertical transport modelling may be applied at different scales. In a simple onedimensional water column model, the goal may be to investigate the timescale of settling or surfacing for a specific type of particle. However, accurate modelling of vertical transport is also key to predicting horizontal transport at large scales, due for example to the vertical variability of horizontal ocean currents (Röhrs et al., 2018; Wichmann et al., 2019). Hence, a good description of vertical behaviour is an essential part of any threedimensional model.
Commonly used transport models may be divided into two classes: Eulerian and Lagrangian. Eulerian methods consist of solving the advection–diffusion–reaction equation directly with numerical methods for partial differential equations. Ocean circulation models, such as ROMS (Shchepetkin and McWilliams, 2005) or NEMO (Madec et al., 2022), are typical examples of Eulerian models. Related are computational fluid dynamics (CFD) models, which can be used, for example, to model waves and turbulence on smaller scales (Cui et al., 2020). Eulerian models are also used to compute the transport of suspended particles in water, for example natural sediments (Warner et al., 2008), nanoparticles in rivers (Saharia et al., 2019), or microplastics in the oceans (Mountford and Morales Maqueda, 2019).
A challenge of the Eulerian approach, when it comes to particles with a distribution of rise or settling speeds, is that one must use discrete speed classes and solve the advection–diffusion equation for each class, forming a large system of equations to solve simultaneously. Hence, it is of interest to know how many classes are needed to achieve the desired accuracy.
Lagrangian particle methods are quite a popular choice for modelling particle transport in the ocean (van Sebille et al., 2018). In these methods, numerical particles, also called Lagrangian elements, are used to represent physical particles. The numerical particles will move with the current, may exhibit some form of random displacement to model turbulent diffusion, and may have a vertical rise or settling speed (or even active swimming behaviour). Lagrangian particle methods have been applied to transport modelling of a wide range of particle types and processes, including plastics (Delandmeter and Van Sebille, 2019; de La Fuente et al., 2021; Fischer et al., 2022), the residence time of water masses (Dugstad et al., 2019), the transport and sedimentation of mineral particles in the ocean (Dissanayake et al., 2014; Nepstad et al., 2020), the surfacing and entrainment of oil (Cui et al., 2018; Nordam et al., 2019b), the transport of dissolved gases (Dissanayake et al., 2012; Wimalaratne et al., 2015), produced water (Nepstad et al., 2022), harmful algae (Rowe et al., 2016), dinocysts and foraminifera (Van Sebille et al., 2015; Nooteboom et al., 2019), fish eggs (Sundby, 1983; Röhrs et al., 2014), and even fish (Scutt Phillips et al., 2018). Stochastic particle transport methods are also used in many other fields of science, including the study of cosmic rays and highenergy particles from the sun (Strauss and Effenberger, 2017), the deposition of inhaled nanoparticles in the airways (Longest and Xi, 2007), and transport modelling of airborne virus transmission (Abuhegazy et al., 2020).
One of the advantages of a Lagrangian approach to particle transport modelling is the ability to represent a wide range of properties or behaviours. By letting each numerical particle with its own properties move independently of the others, one can model physical particle distributions where the sizes, and hence the terminal velocities, can vary by several orders of magnitude.
In practice, Eulerian and Lagrangian models may offer complementary benefits, and they are often used together. This can be as simple as forcing a Lagrangian model with precalculated fields from an Eulerian model (offline Lagrangian modelling; see e.g. Young et al., 2015), or it can involve running a Lagrangian model as part of an Eulerian model (online Lagrangian modelling). In the latter case, it is also possible for the Lagrangian particles to impact the Eulerian model. As examples, Marsh et al. (2015) and Matsumura and Ohshima (2015) model, respectively, icebergs and frazil ice as Lagrangian particles with feedback to the Eulerian model.
The purpose of this paper is to compare and discuss Eulerian and Lagrangian methods, with a focus on the numerical implementation of the models including different boundary conditions and a reaction term. We demonstrate that the two implementations give the same results, and we also address questions of efficiency and convergence. The methods are illustrated using three different onedimensional cases as a basis for the discussion: fish eggs, microplastics, and oil droplets. The cases have been chosen because they represent simplified (but realistic) cases, and they highlight how particles can interact with the boundaries in different ways.
We consider the water column without background flow and investigate the vertical transport of different types of particles. In Sect. 2, we describe the advection–diffusion equation, which is a partial differential equation (PDE) that will form the basis of our transport problems. Additionally, we present a stochastic differential equation (SDE), which yields a Lagrangian particle method that is equivalent to the advection–diffusion equation. In Sect. 3, we briefly outline how we implement an Eulerian model based on numerically solving the advection–diffusion equation with PDE methods, and similarly, we show how we implement a Lagrangian model based on numerically solving the equivalent SDE for a large ensemble of particles.
In Sect. 4, we apply our Eulerian and Lagrangian water column models to three different cases: fish eggs, microplastics, and oil droplets. All three cases are modelled both with the Eulerian approach and the Lagrangian approach. In each case, the particles considered have a distribution of rising and/or sinking speeds. Additionally, the cases serve to highlight the effects of different boundary conditions and a source term. In Sect. 5, we discuss and summarise our comparison of the Eulerian and Lagrangian approaches. Finally, in Sect. 6, we present some concluding remarks.
Some additional details are given in the Appendices: Appendix A contains detailed descriptions of our Eulerian and Lagrangian model implementations, Appendix B has details on how we obtained a distribution of terminal velocities for microplastics, and Appendix C shows additional numerical convergence results.
2.1 The advection–diffusion equation
As our starting point, we assume that the movement of a collection of particles with different rise or settling velocities in the water column may be described by the advection–diffusion equation (see e.g. Hundsdorfer and Verwer, 2003). If C_{k}(z,t) is the concentration of particles with (constant) vertical rise or settling speed w_{k}, then we have
where K(z) is the diffusivity. To represent a distribution of particle sizes, densities, and/or shapes and thus a distribution of rise speeds, one must, in principle, solve a set of advection–diffusion equations, with one equation for each particle speed w_{k}.
2.2 Equation of motion of Lagrangian particles
A mathematically equivalent formulation of the advection–diffusion problem is to consider an ensemble of numerical particles, whose positions develop in time according to the stochastic differential equation (SDE):
Here, W_{t} is a standard Wiener process (Kloeden and Platen, 1992, p. 40); K(z) is again the diffusivity as a function of depth; K^{′}(z) is its derivative with respect to z; and w is the terminal velocity due to buoyancy, which is constant for a given particle.
The connection between Eqs. (1) and (2) is that the probability distribution for the position of particles moving according to Eq. (2) will develop according to Eq. (1). A solution, z(t), to Eq. (2) is called an Itô diffusion, and Eq. (1) is the corresponding Fokker–Planck equation (Kloeden and Platen, 1992, p. 37; see also Appendix A of Nordam et al., 2019b). Each numerical particle position, determined from the SDE, represents a sample from the probability distribution of particle positions, which develops according to the PDE. Solving Eq. (2) for a sufficient number of particles allows us to approximate the probability distribution, which is proportional to the concentration in the Eulerian scheme.
2.3 Velocity distribution
In many practical applications, one will have a size distribution of particles obtained from measurements, a model, or other estimates. For the model we consider here, however, the relevant property is not the size itself but the terminal rise or settling speed of the particle, which is a function of size (and other particle parameters such as shape and density, as well as the viscosity of the ambient fluid). Note that we assume here that a particle will immediately attain its terminal velocity and that particle motion can be described as the combination of a constant terminal velocity and a series of random displacements. This is a commonly used approximation, which holds well if the timescale needed to reach terminal velocity is short compared to other relevant timescales. As an example, consider a particle with some initial speed v_{0}, moving in water under the influence of Stokes' drag (see Appendix B). Then we have that the initial speed will decay exponentially, as given by
where $\mathit{\tau}=\frac{m}{\mathrm{6}\mathit{\pi}\mathit{\mu}r}$, with μ being the dynamic viscosity of water, r being the radius of the particle, and m being the mass of the particle. For typical values for small particles in water, the timescale τ is less than a second (e.g. μ=0.0017 Pa s, r=1 mm, and $m=\mathrm{4}\times {\mathrm{10}}^{\mathrm{6}}$ kg give τ=0.14 s).
In the example cases we consider in Sect. 4, we will assume that we have access to the distribution of terminal velocities, either by directly specifying a functional form of the velocity distribution or by numerical means through mapping the particle size distribution and other properties to velocities.
2.4 Boundary conditions
Depending on the application, different boundary conditions may be used to control the fluxes across the domain boundaries in a water column model. We will deal separately with the diffusive flux,
and the advective flux
For the boundary conditions for the diffusive flux, we first observe that diffusion in our model represents the mixing that occurs due to the combination of random turbulent fluctuations and molecular diffusion (Thorpe, 2005, p. 20). Being caused by the motion of the water, the diffusive flux should not allow suspended particles to leave the water column. Hence, we wish to enforce zero diffusive flux, j_{D}=0, across the boundaries at both the surface and the seafloor. This applies in all the cases we consider in this paper.
Next, we consider the advective flux. In our model, the advection velocity w_{k} represents the terminal rising or sinking velocity of a particle due to buoyancy. Depending on the application, this velocity may or may not allow particles to leave the water column. As an example, consider positively buoyant fish eggs. When fish eggs rise to the surface, they can of course rise no further. However, they remain submerged, having a hydrophilic surface and a density only slightly less than that of the surrounding seawater. Hence, we have chosen to model fish eggs with a zeroflux boundary condition at the surface. With the surface being at z=0 (and z being negative downwards), the total flux j_{tot}(z) at the surface is then given by
As another example, consider positively buoyant oil droplets. When oil droplets rise to the surface, they go from being individual droplets surrounded by water to forming small patches of floating oil at the surface or possibly to merging with a larger surface slick. When this happens, the oil droplets are no longer subject to random motion due to turbulent diffusion, and the oil will remain at the surface until some highenergy event like a breaking wave causes the surface slick to break up and form droplets.
To express the mechanism of oil droplets surfacing and leaving the water column, we have two options: we can use a loss term, which removes oil droplets close to the surface at some rate (Tkalich and Chan, 2002), or we can use a partially absorbing boundary condition. Here, we have chosen to allow the advective flux to carry oil droplets across the boundary at the air–sea interface while still enforcing zero diffusive flux across the surface:
The reasons we describe the surfacing of oil droplets through the boundary condition instead of a loss term are that it is straightforward to express the boundary condition in both the Eulerian and the Lagrangian implementations and that we do not have to define what it means to be close to the surface for the loss term. The idea behind allowing an advective flux while forcing the diffusive flux to be zero is that buoyancy is the mechanism that leads to surfacing. Additionally, if the diffusive flux out of the system is nonzero, then increased diffusivity would lead to faster surfacing, which is contrary to observations. For an additional discussion of this point, see Nordam et al. (2019a).
We note that similar reasoning to the above also applies to the boundary at the seafloor. For example, the settling out of negatively buoyant particles (microplastics, sediments, etc.) can be modelled as an advective flux leaving the model domain through the bottom boundary (i.e. settling on the seafloor or being incorporated into sediments). Note also that, here, advective flux means the flux due to the rise or settling speed of the particle and not that due to vertical currents.
2.5 Reaction terms
An additional reaction term ${S}_{k}(z,{C}_{\mathrm{1}},{C}_{\mathrm{2}},\mathrm{\dots},{C}_{{N}_{k}})$ may be included in Eq. (1), in which case it becomes the advection–diffusion–reaction equation. It can describe the addition of mass (a source); the removal of mass (a sink); or the transformation of mass between classes, for example from C_{1} to C_{2}, which means that, in general, the reaction term for class k can depend on the concentration of all classes. In an Eulerian implementation, a reaction term can simply add or remove mass in different locations. In a Lagrangian stochastic particle implementation, Eq. (2) only describes the transport of numerical particles, and reaction terms have to be implemented as a separate step, with details depending on the nature of the reaction.
In this study, we will only consider reaction terms that add mass. The reaction term S has units equal to the time derivative of concentration, and it enters our Eulerian scheme in a straightforward manner as a rate of change in the average concentration inside a cell.
For the Lagrangian implementation, adding particles is done in a stochastic manner, designed to be consistent with the source term in the Eulerian description. For a source term that adds mass at some rate, in some region of the domain, we add a random number of particles at each time step, such that the expected value of mass added at each time step is equal to the integrated source term. The position of each added particle is drawn from a distribution that corresponds to the spatial distribution of the source term. As this applies to only one of the three cases considered (entrainment of oil droplets), the details are given in the description of that case in Sect. 4.3.
2.6 Sampling error in stochastic particle methods
When we solve an advection–diffusion–reaction problem by means of a stochastic Lagrangian particle model; the position (and possibly other properties) of each numerical particle represents a sample from an underlying distribution. When we draw random samples from a distribution and calculate, for example, the mean of the samples, there will be a random error in the sample mean relative to the true (but usually unknown) mean of the distribution. In particular, if we draw N_{p} independent samples from a distribution with finite variance and take the mean of those samples, this is called the sample mean ${\mathit{\mu}}_{{N}_{p}}$, which is an approximation of the true mean μ. By the strong law of large numbers, we have that ${\mathit{\mu}}_{{N}_{p}}\to \mathit{\mu}$ as N_{p}→∞, almost surely (Billingsley, 1979, p. 85). Furthermore, according to the Lindeberg–Lévy theorem (Billingsley, 1979, p. 308), we have
where σ^{2} is the (assumed finite) variance of the distribution from which the samples are drawn, and 𝒩(0,1) is a standard normal distribution with zero mean and unit variance. In other words, in an ensemble of simulations, the error in the sample mean $\left({\mathit{\mu}}_{{N}_{p}}\mathit{\mu}\right)$ will be a Gaussian random number, with a standard deviation of $\sqrt{{\mathit{\sigma}}^{\mathrm{2}}/{N}_{p}}$.
3.1 Eulerian model
There is a wide range of numerical methods for PDEs to choose from for solving the advection–diffusion–reaction equation (see for example Hundsdorfer and Verwer, 2003). Here, we have chosen to use a finitevolume method (FVM) for discretisation in space and the Crank–Nicolson (implicit trapezoidal) scheme for discretisation in time. An FVM was selected due to the ease of applying prescribedflux boundary conditions and due to the inherent conservation of mass. Our chosen approach leads to secondorder convergence in both space and time, and the Crank–Nicolson scheme is also known to be unconditionally stable for the diffusion equation (Versteeg and Malalasekera, 2007, p. 247).
A detailed description of the discretisation, as well as of the boundary conditions and the numerical solver, is given in Appendix A. Below, we focus on the three most relevant implementation details for the later discussion: a brief review of the boundary conditions, the discretisation of the velocity distribution into classes, and a measure of the numerical error.
3.1.1 Boundary conditions
As mentioned in Sect. 2.4, we use two different types of boundary conditions: zero flux and zero diffusive flux. Using a finitevolume method, our Eulerian implementation represents space as a grid of cells, with the changing concentration in a cell expressed in terms of the sum of the fluxes across the cell faces. The fluxes are approximated numerically from the concentrations in the neighbouring cells, as well as from the diffusivity and the advection velocity for the diffusive and advective fluxes respectively. Details are given in Appendix A1.2.
With a finitevolume method, it is trivial to implement prescribedflux boundary conditions by simply setting one or both of the advective and diffusive fluxes across the cell face at the end of the domain to the prescribed value. In our case, we implement noflux boundary conditions by setting both fluxes to zero across the boundary and nodiffusiveflux boundary conditions by setting only the diffusive flux to zero, leaving the advective flux unchanged.
3.1.2 Discrete representation of velocity distribution
In this study, we are interested in the scenario where our particles have a distribution of rising and/or settling velocities. To investigate the convergence of our Eulerian method with an increasing number of velocity classes, we need an automated way to represent the velocity distribution as a given number of discrete classes. Depending on the application, different ways of dividing the velocity distribution into intervals may be preferable.
In our selected approach, we first choose lower and upper limits: w_{min} and w_{max}. We then divide the interval between the limits into N_{k} classes, with equal spacing on either a linear or logarithmic scale, depending on the case. The representative velocity for each class is chosen to be the midpoint of the class on a linear or logarithmic scale. For cases where the relevant velocities span several orders of magnitude, a logarithmic spacing may be preferable. Where both negative (sinking) and positive (rising) velocities are represented, they are handled separately if a logarithmic spacing is chosen. In Sect. 4, we will describe in detail how the velocity distribution is represented for each of the considered cases.
To set up the initial condition, the total mass must be divided among the different velocity classes. The velocity distribution might be specified directly (see e.g. Sundby, 1983), or it could be inferred from another distribution, most commonly the size distribution combined with particle density and, in some cases, other particle properties such as shape.
To set up a case with N_{k} velocity classes and constant linear spacing, we use a constant class width Δw such that ${w}_{k}={w}_{\mathrm{1}}+(k\mathrm{1})\mathrm{\Delta}w$. All particles with velocities in the interval $[{w}_{k}\mathrm{\Delta}w/\mathrm{2},{w}_{k}+\mathrm{\Delta}w/\mathrm{2})$ belong to class k and are represented by the velocity w_{k} in the Eulerian implementation. The fraction of the total mass belonging to this class is given by
where p(w) is the velocity distribution, normalised such that
Note that ${w}_{min}={w}_{\mathrm{1}}\mathrm{\Delta}w/\mathrm{2}$, and ${w}_{max}={w}_{{N}_{k}}+\mathrm{\Delta}w/\mathrm{2}$.
With equal spacing on a logarithmic scale, the calculation of mass fractions is similar, but the limits on the integral are different. For constant logarithmic spacing δw, we have it that $\frac{{w}_{k+\mathrm{1}}}{{w}_{k}}=\mathit{\delta}w$ for $k=\mathrm{1},\mathrm{\dots},{N}_{k}\mathrm{1}$, and the particles belonging to class k are those with velocities in the interval $[{w}_{k}/\sqrt{\mathit{\delta}w},{w}_{k}\cdot \sqrt{\mathit{\delta}w})$. For completeness, we note that, in this case, ${w}_{min}={w}_{\mathrm{1}}/\sqrt{\mathit{\delta}w}$, and ${w}_{max}={w}_{{N}_{k}}\cdot \sqrt{\mathit{\delta}w}$; also, here, δw is a dimensionless scale factor.
3.1.3 Measures of the numerical error
To compare our Eulerian and Lagrangian implementations, we will present direct comparisons of the predicted concentration, while an investigation of the convergence of our solutions as functions of the different numerical parameters requires a quantitative measure of the error. We have chosen to use the first moment (i.e. the centre of mass) of the distribution. The reason for this choice is primarily due to numerical methods for SDEs usually having a welldefined weak convergence in terms of the moments of the distribution (Kloeden and Platen, 1992, p. 327), making this a natural parameter to consider for the Lagrangian method. We chose to use the same metric for the Eulerian method to facilitate a direct comparison. When used to consider timescales for the surfacing or settling of particles, one could also argue that the depth of the centre of mass of a distribution is a relevant output from a water column model.
The first moment of a distribution p(z) on the interval $z\in [H,\mathrm{0}]$ (where the depth H is a negative number) is given by
where the distribution is normalised such that ${\int}_{H}^{\mathrm{0}}p\left(z\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}z=\mathrm{1}$.
Numerically, we approximate the first moment for the Eulerian results by calculating the integral in Eq. (11) with the standard midpoint Riemann sum quadrature method (Wendroff, 1969, p. 23) and by taking the sum over all velocity classes:
where Δz is the cell size in the spatial discretisation, C_{k,j} is the concentration of particles with velocity w_{k} in cell j, and z_{j} is the midpoint of cell j (see Sect. A1.1 in the Appendix for a description of the spatial discretisation).
As we do not have analytical solutions for the case studies we consider in Sect. 4, the convergence analysis is conducted purely with numerical results. To consider the convergence with the number of classes N_{k} for the Eulerian implementation, we approximate the error in the numerically calculated first moment M_{1}(N_{k}) as follows:
where M_{1} is the true (but unknown) first moment of the distribution, and ${M}_{\mathrm{1}}\left({N}_{k}^{\text{ref}}\right)$ is a reference solution calculated with a large number of classes. The other numerical parameters (time step and spatial discretisation) are kept constant, while N_{k} is varied. As long as ${N}_{k}\ll {N}_{k}^{\text{ref}}$, we assume that the approximation in Eq. (13) is a good one. See also Appendix C for additional details on this point, including convergence tests used to select suitable values of the time step and spatial discretisation.
3.2 Lagrangian model
The starting point for our Lagrangian implementation for solving the advection–diffusion equation is the SDE given by Eq. (2). We solve this equation a large number of times, where each solution describes the trajectory of a numerical particle, and the concentration can be found from the distribution of the numerical particles.
For the numerical solution, we use the Euler–Maruyama method (Kloeden and Platen, 1992, p. 305), which gives
Here, ΔW_{i} refers to independent Gaussian distributed random numbers with zero mean, 〈ΔW_{i}〉=0, and variance $\langle \mathrm{\Delta}{W}_{i}^{\mathrm{2}}\rangle =\mathrm{\Delta}t$. The Euler–Maruyama method has an order of convergence of 1 in the weak sense, which means that, for sufficiently small Δt, we have that
where z(t_{n}) is the true (but usually unknown) solution at time t_{n}; z_{n} is the numerical solution at time t_{n}; 〈 〉 means the ensemble average (expectation value); f(z) is a 4times differentiable function of, at most, polynomial growth; and γ is some constant (Kloeden and Platen, 1992, p. 327).
If a method convergences in the weak sense, this implies that the moments of the modelled distribution converge to the moments of the true distribution as Δt→0. Hence, a weakly convergent method means that the distribution of numerical solutions will converge to the true distribution, which means that we can use the method to predict concentrations. There exist methods which have higher orders of weak convergence, but we have chosen to consider the Euler–Maruyama method as it is commonly seen in the literature. See e.g. Gräwe (2011) and Gräwe et al. (2012) for reviews of higherorder methods in the context of water column models.
3.2.1 Boundary conditions
As described in Sect. 2.4, different boundary conditions are relevant for different applications. Here, we wish to employ zeroflux boundary conditions and zerodiffusiveflux boundary conditions. We also need to implement the same boundary conditions in both the Eulerian and the Lagrangian schemes. While the implementation of boundary conditions in numerical PDE methods is well known, the issue of boundary conditions in Lagrangian models is less well established and not as well described in the literature, except in simpler cases such as those with no advection and constant diffusivity (see e.g. Gillespie and Seitaridou, 2012). Hence, we describe our approach in some detail.
We split our numerical scheme (Eq. 14) into two parts and treat the advection and diffusion separately. To implement noflux boundary conditions at the sea surface (z=0) for positively buoyant particles, we use the following sequence of steps in the particle model:

random displacement ($z\to z+{K}^{\prime}\left(z\right)\mathrm{\Delta}t+\sqrt{\mathrm{2}K\left(z\right)}\mathrm{\Delta}W$)

reflection at the surface ($z\to \leftz\right$)

rise due to buoyancy ($z\to z+w\mathrm{\Delta}t$)

stop particles at the surface ($z\to max(z,\mathrm{0})$).
For noflux boundary conditions at the seafloor for negatively buoyant particles, the steps are analogous but with reflection and stopping at the bottom boundary instead.
To implement a nodiffusiveflux boundary condition while allowing an advective flux across the boundary at z=0, as discussed in Sect. 2.4, we instead use the following sequence of steps in the particle model:

random displacement ($z\to z+{K}^{\prime}\left(z\right)\mathrm{\Delta}t+\sqrt{\mathrm{2}K\left(z\right)}\mathrm{\Delta}W$)

reflection at the surface ($z\to \leftz\right$)

rise due to buoyancy ($z\to z+w\mathrm{\Delta}t$)

remove particles above the surface (z→Surfaced).
When particles are set to surfaced in the last step, they are removed from the simulation. To get the same boundary conditions for e.g. sinking plastics or mineral particles that settle on the seafloor, we again have to reflect at the seafloor in the second step and remove those particles that settle out when they cross the seafloor boundary due to sinking in the last step shown above. In Nordam et al. (2019a), this approach to implementing boundary conditions in a Lagrangian implementation was found to give identical results to an Eulerian implementation using prescribed fluxes in an FVM.
Depending on the application, particles that have been removed from the simulation may be reintroduced with some probability to represent, for example, breaking waves entraining material from the water surface or strong currents resuspending material from the seafloor.
3.2.2 Velocity distribution
In a Lagrangian model, a distribution of terminal rising or sinking velocities can be represented very naturally simply by allowing each particle to have its own vertical velocity drawn from the desired distribution. Hence, any distribution can be represented, and the quality of the representation will depend on the number of particles. The requirement for assigning a random velocity to each numerical particle is that we can draw samples from the velocity probability distribution. Depending on the available information, different approaches might be suitable (see e.g. Press et al., 2007, Chap. 7). For the three case studies considered in Sect. 4, we use three different approaches. A justification for the selection and details are given in the description for each case.
3.2.3 Obtaining a concentration profile from a distribution of particles
Mathematically, the link between the Lagrangian and the Eulerian methods is that each solution of the SDE in the Lagrangian method represents a sample from the distribution that evolves under the advection–diffusion equation in the Eulerian method. However, in practical applications, we are often interested in the distribution itself and not just samples. Numerous methods exist for reconstructing the distribution from samples as this is a common problem not just in applied geoscience but in statistics in general. For a review of a few different approaches in the context of applied oceanography, see e.g. Lynch et al. (2014, Chap. 8).
Here, we have chosen to use the simple approach of constructing a histogram of particle positions. This works well when the number of particles is large relative to the number of cells, as is the case for our onedimensional examples below. An additional benefit is that the only free parameter is the bin size of the histogram, and a very direct comparison between Lagrangian and Eulerian results is obtained by setting the bin size to be equal to the cell size of the Eulerian approach.
3.2.4 Measures of the numerical error
As described in Sect. 3.1.3, we will use a direct comparison of particle concentration, as well as of the first moment of the concentration profile, to assess our results. Calculating the first moment (centre of mass) of the spatial distribution of particles is trivial as this is simply the average position over all particles:
where N_{p} is the number of particles, and z_{j} is the position of particle j. Here we have assumed that all particles represent an equal mass, but if this is not the case, the average is simply weighted by the mass m_{j} of each particle.
For the Lagrangian implementation, convergence with the number of particles N_{p} is of a stochastic rather than a deterministic nature. Running a simulation once will give an approximate result for e.g. the first moment but with some random error. Running the simulation again will give another result. To investigate the convergence with the number of particles, we ran 100 repeated simulations for each value of N_{p}. For each of those 100 simulations, we calculate the first moment M_{1}(N_{p}). As described in Sect. 2.6, M_{1}(N_{p}) will be a Gaussian random variable whose mean is the true (but unknown) mean of the distribution M_{1}, with a standard deviation of $\mathit{\sigma}\left({N}_{p}\right)\sim \mathrm{1}/\sqrt{{N}_{p}}$. Hence, we consider σ(N_{p}) as a measure of the error when considering convergence with the number of particles in the Lagrangian implementation.
We present three different cases and simulate all three with both an Eulerian and a Lagrangian implementation of our model. Then, we compare the results and discuss convergence in terms of numerical parameters for both implementations. The cases are described below, and an overview of some aspects of the cases is presented in Table 1. Table 2 lists the numerical parameters we have used to investigate convergence.
For case 1, we consider positively buoyant fish eggs with a distribution of rise speeds. Fish eggs will float to the surface but do not leave the water column; hence, we model these with a noflux boundary condition at the surface (i.e. the particles cannot leave the water column).
For case 2, we consider microplastic particles with a velocity distribution that includes both rising and sinking particles, representing the diversity of densities associated with different polymer types. In this case, we again use a noflux boundary at the surface, while the boundary at the seabed is zero flux in diffusion, though an advective flux is allowed to leave the domain. This represents negatively buoyant microplastic particles that are removed from the water column by settling onto the seabed.
For case 3, we consider oil which is entrained as droplets when a surface slick is broken up by waves. The droplets are all positively buoyant with the same density but have a size distribution which varies in time and which leads to a distribution of rising speeds. In this case, the boundary condition at the surface is zero diffusive flux, though an advective flux that represents droplets that merge with the surface slick is allowed. Additionally, we include the effect of entrainment of oil from the surface slick by breaking waves in our model.
A spatially variable diffusivity profile is used in all three cases. To calculate eddy diffusivity as a function of depth, we have used the GOTM turbulence model (Umlauf et al., 2005). The model was forced with a surface stress corresponding to a wind speed of 9 m s^{−1} (Gill, 1982, p. 29), and it was run with a kω turbulence closure with a flux condition for turbulent kinetic energy (TKE) at the surface. The TKE flux at the surface accounts for wave breaking (Craig and Banner, 1994).
To support regridding and differentiation, an analytical expression has been fitted to the discrete diffusivity output from GOTM. In Fig. 1, the fitted analytical profile is shown together with the output from GOTM. The analytical expression is given by
where β=0.00636 m s^{−1}, γ=0.088 m^{−1}, δ=1.54, z_{0}=1.3 m, and z is negative downwards. The choice of using an analytical profile was mainly pragmatic: for the Lagrangian model, the diffusivity and its derivative must be evaluated at arbitrary positions; hence, some form of interpolation with continuous derivatives is needed in any case. Additionally, we made the profile nonzero at the surface, which seems physically reasonable for conditions with some wave breaking. Also, if the diffusivity is zero at the surface, positively buoyant particles would eventually get stuck at the surface, and this is not observed in surveys of e.g. fish eggs (Sundby, 1983).
Note that both implementations make use of the analytical diffusivity profile: in the Eulerian implementation, the diffusivity is given by the value of Eq. (17) evaluated at the cell boundaries, while for the Lagrangian implementation, Eq. (17) is evaluated at the position of each particle at each time step, and the derivative is evaluated by secondorder central finite difference.
For each of the three cases, our initial condition will be a submerged particle distribution. The spatial distribution will be a Gaussian distribution centred at −20 m depth, with a standard deviation of 4 m. The distribution of rise and settling speeds will be different for each case and will be explained in the description of each case below.
In all cases, we will present a direct comparison of the predicted total water column concentrations at different times from the Eulerian and the Lagrangian implementations. Total concentration in this case means that it is integrated over the velocity distribution, showing the total suspended concentration regardless of rising or settling velocity. Additionally, we show a numerical convergence analysis separately for the two schemes. For the Eulerian implementation, we consider the error to be a function of the number of velocity classes. For the Lagrangian implementation, we present the error as a function of the number of particles. As a measure of the error, we consider the first moment (the centre of mass) of the spatial concentration distribution (see Sect. 3.1.3 and 3.2.4). Additional numerical convergence tests are also presented in Appendix C.
4.1 Case 1: pelagic fish eggs
In case 1, we consider pelagic fish eggs with a distribution of rise velocities. These positively buoyant particles will rise towards the surface but do not form a surface slick. Rather, they rise to the surface in stagnant water but stay submerged and can be mixed back down by the eddy diffusivity (Sundby, 1983; Röhrs et al., 2014; Sundby and Kristiansen, 2015).
Taking an example from Sundby (1983, Table 3), we consider the eggs of the ArctoNorwegian cod (Gadhus morhua L.), which is typically found to have a Gaussian distribution of terminal rise velocities with average $\overline{w}=\mathrm{0.96}$ mm s^{−1} and a standard deviation of σ=0.38 mm s^{−1}. We furthermore assume that the distribution is truncated symmetrically about the mean such that $w\in [\overline{w}\pm \mathrm{2}\mathit{\sigma}]$. To divide this distribution into N_{k} velocity classes for use in the Eulerian implementation, we divide the interval $[\overline{w}\pm \mathrm{2}\mathit{\sigma}]$ into N_{k} subintervals with equal spacing Δw on a linear scale, and we calculate the mass fraction in each class, as described in Sect. 3.1.2. For the Lagrangian implementation, we sample the particle velocities at random from the truncated Gaussian distribution. Hence, a numerical particle may have any rise speed $w\in [\overline{w}\pm \mathrm{2}\mathit{\sigma}]$.
The results of our simulations for fish eggs are shown in Fig. 2. In panel (a), concentration profiles are shown for four different times (including t=0), with the Lagrangian results shown as continuous lines and the Eulerian results for the same times shown as dashed black lines. The concentration profiles from the Eulerian implementation were produced with time step Δt_{E}=10 s, N_{k}=128 classes, and N_{z}=8000 spatial cells, while the Lagrangian results used time step Δt_{L}=10 s and ${N}_{p}=\mathrm{3}\times {\mathrm{10}}^{\mathrm{6}}$ particles converted to a concentration by bin count for 8000 bins. We observe that the two implementations produce essentially identical predicted concentrations.
In panel (b) of Fig. 2, we show the convergence of the Eulerian results with the number of velocity classes, with the other numerical parameters kept constant at time step Δt_{E}=10 s and N_{z}=8000 spatial grid cells. The plot shows the error in the first moment, where the error is calculated relative to a reference solution from a simulation with N_{k}=256 classes. We observe that the convergence with the number of classes appears to be of order 2, with the error going down as $\mathrm{1}/{N}_{k}^{\mathrm{2}}$.
Finally, in panel (c), we show the convergence of the Lagrangian results as a function of the number of particles N_{p}, keeping a fixed time step of Δt_{L}=10 s. The figure shows the standard deviation of the modelled first moment, calculated over an ensemble of 100 simulations for each value of N_{p}. As noted above, the first moment is just the average particle position, and as expected from the Lindeberg–Lévy theorem, the standard deviation of the first moment goes down as $\mathrm{1}/\sqrt{{N}_{p}}$ (see Sect. 2.6).
4.2 Case 2: microplastics
In this case, we consider a distribution of microplastic particles. As the only parameter describing the numerical particles in our model is the terminal rise and settling velocities, we have obtained a velocity distribution based on published descriptions of microplastics.
As our starting point, we consider particles with a distribution of densities from 0.8 to 1.5 kg L^{−1}, a distribution of sizes from 20 µm to 5 mm, and a distribution of shapes described by the Corey shape factor. We assume the density of seawater to be 1.025 kg L^{−1}. The distributions of density, size, and shape are taken from Kooi and Koelmans (2019). The terminal rise and settling velocities of these particles are calculated from empirical relations obtained by Waldschläger and Schüttrumpf (2019), which give terminal velocities as a function of density, size, and shape. To obtain a distribution of velocities, we used a Monte Carlo method that consists of generating a large number of random combinations of density, size, and shape, drawn from the distributions given in Kooi and Koelmans (2019), and mapping those to velocities via the relationships given by Waldschläger and Schüttrumpf (2019).
For the Eulerian simulations, discretisations of the velocity distribution into different numbers of classes were obtained as histograms with different numbers of bins of 10^{9} randomly generated terminal velocities. For the Lagrangian simulations, the terminal velocities of the particles were obtained by directly generating random velocities as described above. A detailed description of the steps involved in obtaining the terminal velocities for microplastics is given in Appendix B. See also Isachenko (2020) for another example of a similar approach.
In this case, particles are allowed to leave the water column via the seafloor at 50 m depth, where we enforce zero diffusive flux while allowing the advective flux to carry particles across the boundary. This represents particles that leave the water column (and thus the simulation) by settling onto the seafloor, where they are no longer able to be resuspended by the eddy diffusivity. Different resuspension mechanisms can of course be included in the model, but we have chosen to ignore that here.
The results of our simulations for microplastics are shown in Fig. 3. In panel (a), concentration profiles are shown for four different times, with the Lagrangian results shown as continuous lines and the Eulerian results for the same times shown as dashed black lines. As in case 1, the concentration profiles from the Eulerian implementation were produced with time step Δt_{E}=10 s, N_{k}=128 classes, and N_{z}=8000 spatial cells, while the Lagrangian results used time step Δt_{L}=10 s and ${N}_{p}=\mathrm{3}\times {\mathrm{10}}^{\mathrm{6}}$ particles converted to concentration by bin count for 8000 bins. Again, we observe that the two implementations give very similar predictions as the curves are visually indistinguishable.
In panel (b), we show the convergence of the first moment in the Eulerian results as a function of the number of classes. The $\mathrm{1}/{N}_{k}^{\mathrm{2}}$ trend is less clear than in case 1, and the size of the error appears to be somewhat sensitive to the number of classes – for example, 32 classes give a smaller error than 36 classes for two of the times considered.
In panel (c), we show the convergence with the number of particles of the standard deviation of the first moment for the Lagrangian implementation. Again, as expected from the discussion in Sect. 2.6, this follows a $\mathrm{1}/\sqrt{{N}_{p}}$ trend.
In contrast to case 1, here we have particles with both positive and negative terminal velocities. As the particles are allowed to leave the water column, the total mass in suspension (i.e. the integral of the concentration profile) decreases with time. Hence, we also present the remaining suspended mass as a function of time. Panel (a) of Fig. 4 shows the remaining mass fraction suspended in the water column for both implementations, while panel (b) shows the difference between the two implementations for 10 different Lagrangian runs. The oscillatory nature of the difference is due to the randomness of the stochastic Lagrangian particle model. The results were produced with the same numerical parameters as stated in the previous paragraph.
4.3 Case 3: oil droplets with breakingwave entrainment
In this case, we consider a simplified onedimensional oil spill model, which includes entrainment of surface oil by breaking waves; rising of oil droplets due to buoyancy; turbulent diffusive mixing; and oil droplets rising to the surface, creating a surface slick. As the entrainment mechanism is unique to this case, it is described in some detail below.
To model the entrainment of surface oil by breaking waves, we need the entrainment rate, the entrainment depth, and the size distribution of the entrained droplets. We model entrainment as a firstorder decay process (Johansen et al., 2015):
where Q_{s} is the amount of oil at the surface, and α is the entrainment rate coefficient. The coefficient α is calculated according to Li et al. (2017). For the entrainment depth, we use the classic result by Delvigne and Sweeney (1988), where the entrained oil is evenly distributed across the interval
H_{s} is the significant wave height, here assumed to be 1.74 m. The size distribution of the entrained droplets is also taken from Johansen et al. (2015) using modified Weber number scaling. The size distribution is lognormal, with a fixedwidth parameter and a median size d_{50}, which depends on the wave height and the properties of the oil. For a detailed description of this scheme, as well as of the values of the relevant oil properties used here, see Nordam et al. (2019b).
From Eq. (18), we see that the entrainment rate depends on the amount of oil at the surface. Since all the oil in this model is either submerged or at the surface, the amount of oil in the surface slick at time t, Q_{s}(t), can be found as follows:
where Q_{tot} is the total amount of oil present.
Using uniform entrainment throughout the interval $z\in [\mathrm{1.15}{H}_{\mathrm{s}},\mathrm{1.85}{H}_{\mathrm{s}}]$ (Delvigne and Sweeney, 1988), we find the source term in the Eulerian advection–diffusion–reaction equation from the entrainment rate and the length of the entrainment interval:
In the Lagrangian implementation, we use the approximation that the analytical solution of Eq. (18) for constant α is an exponential decay with lifetime $\mathit{\tau}=\mathrm{1}/\mathit{\alpha}$. This means we submerge a fraction of the surface slick at every time step, with that fraction being given by
Thus, entrainment is implemented as a stochastic process where particles that have surfaced are reentrained with probability p at every time step. If entrained, a particle is assigned a depth drawn from a uniform random distribution at the interval $z\in [\mathrm{1.15}{H}_{\mathrm{s}},\mathrm{1.85}{H}_{\mathrm{s}}]$ and a random size drawn from the lognormal size distribution described above.
Surfacing of oil droplets is described as a boundary condition where an advective flux is allowed to pass through the boundary at the surface while enforcing zero diffusive flux, as described in Sect. 2.4. For additional details on the implementation of surfacing and entrainment of oil in a Lagrangian model, see also Nordam et al. (2019a, b).
The results of our simulations for oil droplets are shown in Fig. 5. In panel (a), concentration profiles are shown for four different times, with the Lagrangian results shown as continuous lines and the Eulerian results for the same times shown as dashed black lines. As in cases 1 and 2, the concentration profiles from the Eulerian implementation were produced with time step Δt_{E}=10 s, N_{k}=128 classes, and N_{z}=8000 spatial cells, while the Lagrangian results used time step Δt_{L}=2 s and ${N}_{p}=\mathrm{3}\times {\mathrm{10}}^{\mathrm{6}}$ particles converted to concentration by bin count for 8000 bins.
In panel (b), we show the convergence of the Eulerian results as a function of the number of classes. As in case 1, the error scales with $\mathrm{1}/{N}_{k}^{\mathrm{2}}$, though here, the absolute value of the error is about 2 orders of magnitude larger. Panel (c) shows the convergence of the Lagrangian results with the number of particles, which, as before, follows the $\mathrm{1}/\sqrt{{N}_{p}}$ trend closely.
The total amount of submerged oil changes in time as oil droplets surface and are reentrained. Over time, the submerged size distribution shifts towards smaller droplets as these take longer to reach the surface, causing the centre of mass to shift downwards over time. As before, the two schemes give very similar results. As for case 2, we also present here the remaining suspended mass as a function of time. Panel (a) of Fig. 6 shows the remaining mass fraction suspended in the water column for both implementations, while panel (b) shows the difference between the implementations for 10 different runs of the Lagrangian implementation. The oscillatory nature of the difference is due to the randomness of the stochastic Lagrangian particle model. The results were produced with the same numerical parameters as stated in the previous paragraph.
In this paper, we have conducted a comparison of an Eulerian and a Lagrangian implementation of a water column model for particles with different distributions of rising and sinking speeds. To highlight different choices of boundary conditions and a reaction term, we have chosen to use fish eggs, microplastics, and oil droplets as our example cases. The model can also easily be applied to other cases, such as mineral particles, nanoparticles, algae, and chemicals. More complex reaction terms, such as agglomeration, can also be added to the same modelling framework.
5.1 Boundary conditions
Boundary conditions are always an essential part of the problem for any model based on PDEs, and indeed, the boundary conditions must be specified before the problem can be said to be well posed (Gustafsson, 2008, Chap. 2). The implementation of different boundary conditions in Eulerian models is amply addressed in the applied numerical literature (see e.g. Hundsdorfer and Verwer, 2003; Versteeg and Malalasekera, 2007). Hence, any Eulerian models for environmental transport problems may draw on a wide body of standard reference works describing the implementation of different boundary conditions.
Transport modelling with stochastic particle methods, on the other hand, is perhaps more of a niche endeavour, and there exists less general applied literature on the topic of numerical methods for SDEs compared to PDEs. The standard reference on numerical solution of SDEs, Kloeden and Platen (1992)^{1}, does not directly address the question of boundary conditions but mentions a few references in the section on bibliographical notes (Kloeden and Platen, 1992, pp. 593–594). These are all of a rather technical nature.
In the mathematical literature, an Itô diffusion (of which Eq. 2 is an example) on a domain with a reflecting boundary is known as a Skorokhod problem, after the Ukrainian mathematician Anatoliy V. Skorokhod, who wrote two early papers on the topic (Skorokhod, 1961, 1962). In particular, Skorokhod (1961) describes how an SDE with a reflecting boundary can be formally described by adding a term which is zero everywhere except on the boundary and which causes any trajectory that touches the boundary to be immediately reflected. However, the focus of that paper was on proving the existence and uniqueness of solutions to this modified system and not on numerical solutions.
In the applied literature on atmospheric transport modelling with random flight models, the issue of boundary conditions has received some attention (see e.g. Wilson and Flesch, 1993, and Rodean, 1996, Chap. 11). Less has been written about boundary conditions in the applied literature on random walk models, which are more commonly used in applied oceanography. For reflecting boundary conditions, a pragmatic and common choice is to simply take any particles that have been randomly displaced to a point outside the boundary and to reflect them to an equal distance inside the boundary (see e.g. Israelsson et al., 2006). This approach has also been shown to be exact in the special case of no advection and constant diffusivity (Gillespie and Seitaridou, 2012, pp. 58–60 and 75–76). In the more general case of variable diffusivity, some details are discussed in Ross and Sharples (2004), including a numerical artifact that appears when the diffusivity has a nonzero derivative at the boundary. Nordam et al. (2019a) discuss the separate treatment of boundary conditions for advection and diffusion and conduct a numerical comparison with Neumann and Robin boundary conditions in an Eulerian scheme. In this paper, we have used the same approach as in Nordam et al. (2019a).
We have provided a detailed description of our implementation of two different types of boundary conditions in the Lagrangian scheme. While our implementation may not have a rigorous foundation in the theory of SDEs, we do present a comparison between our Lagrangian and Eulerian implementations, where the Eulerian method uses a standard approach for FVMs. In cases 2 and 3 (microplastics and oil droplets), we show that the two implementations give very similar predictions of the amount of mass which remains suspended in the water column as a function of time, which of course indicates that the net effect of the boundary conditions is the same in the two implementations. For an additional discussion of this topic, see also Nordam et al. (2019a).
5.2 Convergence of numerical results
We have paid special attention to the representation of particles with a distribution of terminal rising and settling velocities. In a Lagrangian particle model, a distribution of terminal velocities may be straightforwardly represented since each numerical particle can have its own velocity. In an Eulerian model, we have to introduce discrete classes with different velocities to represent the true distribution. We have investigated how the error is reduced with an increasing number of particles and classes respectively.
As a measure of the error in the Lagrangian implementation, we have used the standard deviation of the first moment of the position distribution, M_{1} (i.e. the centre of mass). As is usual for a Monte Carlo method, the standard deviation of M_{1}(N_{p}) goes down proportionally to $\mathrm{1}/\sqrt{{N}_{p}}$, where N_{p} is the number of particles. If the particles can be considered as independent samples, this behaviour follows from the Lindeberg–Lévy theorem, as discussed in Sect. 2.6. We note, however, that the usual $\mathrm{1}/\sqrt{{N}_{p}}$ scaling also seems to be followed very closely in case 3, where the particles cannot strictly be said to represent independent samples due to the inclusion of a distributiondependent reaction term.
For the Eulerian implementation, we have considered the error in the modelled first moment M_{1}(N_{k}) as a function of the number of classes N_{k} used to discretise the velocity distribution representing the particles. We found that the error scales as $\mathrm{1}/{N}_{k}^{\mathrm{2}}$, where the velocity of each class is represented by the midpoint of the interval spanned by the class. To the best of our knowledge, this result has not previously been published in this context. Models where a distribution of velocities (e.g. due to a distribution of sizes) is represented by a finite number of discrete classes are often used in both marine and atmospheric transport modelling (see e.g. Tegen and Lacis, 1996; Zender et al., 2003; Wichmann et al., 2019; Cui et al., 2020). Hence, additional understanding of how the number of velocity classes affects the accuracy of predictions should be of some interest within these fields.
We note, however, that the result showing $\mathrm{1}/{N}_{k}^{\mathrm{2}}$ scaling of the error is perhaps not so surprising. Consider as an example the mean of a velocity distribution p(v), which is given by
If the integral is evaluated numerically, there will be a numerical error, which will depend on the chosen approach. Our approach has been to divide the underlying velocity distribution into equalsized bins and to let each bin be represented by its midpoint (on either linear or logarithmic scales, depending on the case). This is equivalent to the midpoint quadrature method, which is known to have an error proportional to Δx^{2}, where Δx is the bin size (Wendroff, 1969, p. 23). This, in turn, is of course proportional to $\mathrm{1}/{N}_{k}^{\mathrm{2}}$, where N_{k} is the number of bins.
5.3 Relative merits of Eulerian and Lagrangian methods
Numerous studies have discussed different aspects and applications of Lagrangian and Eulerian methods and have compared the two approaches (see e.g. Fay et al., 1995; Benson et al., 2017; Nordam et al., 2019a; Nepstad et al., 2022). Here, we consider the question of which implementation should be preferred when modelling particles with a distribution of rising or settling speeds.
First we consider case 1, in which we modelled positively buoyant fish eggs with a Gaussian distribution of terminal rise velocities. For the Eulerian implementation, the numerical error in the first moment due to a finite number of velocity bins starts out between 0.01 and 0.001 m for N_{k}=4 size classes and reduces very consistently as $\mathrm{1}/{N}_{k}^{\mathrm{2}}$ from there (see Fig. 2). A numerical error of less than 0.01 m in the position of the centre of mass is clearly far smaller than the error due to other uncertainties and model approximations. Hence, with such a wellbehaved velocity distribution in this particular case, an Eulerian model with a small number of velocity classes will give more than adequate numerical accuracy.
For the Lagrangian implementation, on the other hand, we need to use 10^{6} particles to reach an error of about 0.01 m, a result of the slower $\mathrm{1}/\sqrt{{N}_{p}}$ convergence. For modelling fish eggs in a water column model, an Eulerian approach seems to be the most efficient choice in terms of the balance between error and computational effort. This is not necessarily the case in a full threedimensional model due to the computational effort involved in tracking the concentration fields for the individual velocity classes across the entire computational domain.
In case 2 with the microplastics, we have a more complicated velocity distribution. The positive and negative parts of the distribution are shown separately on log scales in Fig. B1, where we see that the positively buoyant part of the distribution is bimodal. Due to the very wide range of velocities present, we chose to use logspaced velocity bins (separate for positive and negative velocities) in the Eulerian implementation. Nevertheless, we see from Fig. 3 that the error in the first moment shows some oscillation with the changing number of classes. We also note that the absolute value of the error is far higher than in case 1, and we need around 64 velocity classes to reach an error of 0.01 m in the modelled first moment M_{1}(N_{k}). Compare this to case 1, where the same accuracy was achieved with only four classes.
For the Lagrangian implementation, the results of case 2 look far more similar to those of case 1, and we find that 10^{6} to 3×10^{6} particles will yield a sampling error of 0.01 m. In the onedimensional case, the Eulerian implementation is still more efficient than the Lagrangian implementation, but in three dimensions, this is probably no longer the case. The results will, of course, depend on parameters such as time step, grid resolution, and the choice of numerical solver. However, it is clear that the computational effort of solving the advection–diffusion equation for 64 individual velocity classes is considerable.
In case 3 with the oil droplets, we consider a simplified onedimensional oil spill model. Lagrangian approaches have a long history in oil spill modelling, both for horizontal transport (Tayfun and Wang, 1973) and vertical transport (Elliott et al., 1986). A critical issue in oil spill modelling is that the droplet size distribution is both broad and changing over time. This is easily modelled with a Lagrangian approach, while an Eulerian approach might require either some form of dynamic size class scheme or a very large number of size classes. From panel (b) of Fig. 5, we see that the error in the first moment starts at a little above 1 m for N_{k}=4 classes and goes down fairly consistently as $\mathrm{1}/{N}_{k}^{\mathrm{2}}$. As such, the scaling is very similar to case 1, while the prefactor is again 2–3 orders of magnitude larger, and an error of 0.01 m is achieved with N_{k}=32 classes. The Lagrangian approach, on the other hand, behaves quite similarly to cases 1 and 2 and reaches an error of 0.01 m with 10^{6} particles.
It is important to point out that the errors discussed above are of different natures for the Eulerian and the Lagrangian implementations. In the Eulerian approach, solving with a finite number of speed classes, N_{k}, leads to a systematic error that goes to zero as N_{k}→∞. In contrast, a finite number of particles, N_{p}, leads to a stochastic sampling error in the Lagrangian approach. If the particles are independent, the sampling error is not systematic, and when N_{p}→∞, the standard deviation of the sample error goes to zero as $\mathrm{1}/\sqrt{{N}_{p}}$. See also a further discussion on this point in Sect. 2.6 and Appendix C.
In this paper, we have implemented and compared two different versions of a water column model for particles that undergo diffusion and that rise or sink with different distributions of terminal velocities. Our Eulerian implementation used discrete velocity classes to represent the velocity distribution, while the Lagrangian implementation allowed each particle to have its own velocity.
We have studied the rate of convergence of the two different implementations considering the centre of mass (i.e. the first moment of the concentration profile) as a measure of the error. Our main interest has been to show that we can implement different boundary conditions in an equivalent manner in the two schemes and to demonstrate how the numerical error is reduced with an increasing number of velocity classes and particles for the Eulerian and Lagrangian implementations respectively. Convergence results for varying time steps and (in the Eulerian case) spatial resolutions are also shown in Appendix C.
Three different example cases were considered: positively buoyant fish eggs with a Gaussian velocity distribution, microplastics with a distribution of positive and negative velocities obtained by means of a Monte Carlo approach (see Appendix B), and positively buoyant oil droplets with a timevarying velocity distribution. In all three cases, we find that the stochastic sampling error in the Lagrangian method scales as $\mathrm{1}/\sqrt{{N}_{p}}$ (where N_{p} is the number of particles) with about the same prefactor for all cases. Interestingly, this also appears to hold in the case of oil droplets, where the particles cannot strictly be said to be independent samples.
For the Eulerian implementation, we find that the error appears to scale as $\mathrm{1}/{N}_{k}^{\mathrm{2}}$ (where N_{k} is the number of velocity classes), though the prefactor varies by several orders of magnitude between the cases. We also see that, in the case of microplastics, the exact choice of the velocity classes has a significant impact on the error, possibly due to the multimodal nature of the velocity distribution.
Owing to current shear and density gradients, vertical behaviour is very important for correctly modelling horizontal transport in the ocean. While it might be hard to draw strict conclusions about threedimensional simulations from onedimensional examples, we would still argue that a onedimensional model captures a large part of the complexity of modelling particle distributions: horizontal advection and diffusion affect all particles equally, and the differences in vertical behaviour are well represented in a water column model. Hence, it is clear from our results that the number of classes needed to get accurate results in an Eulerian simulation will depend strongly on the velocity distribution of the relevant particles. We observed that far fewer classes were needed for the fish egg case with a Gaussian velocity distribution than for the microplastic case, which used a wider, multimodal distribution. Similarly, the timedependent size distribution (and thus velocity distribution) seen in the oil case needed more classes to give accurate results than the fish egg case. Hence, our results would suggest that Lagrangian particle models have a particular advantage when wider, less normal velocity distributions are considered or when the size distribution changes with time.
A1 Eulerian implementation
A1.1 Discretisation
The advection–diffusion–reaction (ADR) equation describes how the concentration C_{k} of a substance affected by diffusive, advective, and reactive processes develops in time. In this study, we consider the onedimensional version of the ADR equation and consider vertical motion only. The diffusivity is assumed to be a function of position. Similarly, the advection velocity, which here represents particles rising or sinking due to buoyancy, may also be spatially dependent such that particles stop at the surface and/or bottom boundaries. Finally, the reaction term is assumed to depend on both concentration and position. Consequently, Eq. (1) with all dependencies explicitly shown becomes
where C_{k}(z,t) is the concentration of the material of class k, K(z) is the spatially dependent diffusivity, w_{k}(z) is the spatially dependent velocity of the material in class k, $S(z,{C}_{\mathrm{1}}(z,t),\mathrm{\dots},{C}_{{N}_{k}}(z,t\left)\right)$ is the source term (which may depend on the concentration of all classes), and N_{k} is the number of classes.
In finitevolume methods (FVMs), the PDE is converted to an integral form by integration over control volumes or cells. As illustrated in Fig. A1, cell centres are given integer indices, and cell faces are given halfinteger indices such that cell j extends from ${z}_{j\mathrm{1}/\mathrm{2}}$ to ${z}_{j+\mathrm{1}/\mathrm{2}}$, with the cell width being equal to Δz_{j}. Considering a spatially discretised version of Eq. (A1), we use the divergence theorem to write the volume integrals of the terms on the righthand side as surface integrals over the cell faces. By applying the Leibniz integral rule to change the order of the spatial integral and the differentiation with respect to time on the lefthand side, the average concentration of class k in cell j, ${\stackrel{\mathrm{\u203e}}{C}}_{k,j}$, may be defined. Similarly the average reaction term for class k in cell j, ${\stackrel{\mathrm{\u203e}}{S}}_{k,j}$, is found by integration.
The integral form of Eq. (A1) for the average concentration of the material of class k in cell j thus reads as
where the subscripts on the brackets indicate where the variables are evaluated, e.g. $({w}_{k}{C}_{k}{)}_{j+\mathrm{1}/\mathrm{2}}={w}_{k}({z}_{j+\mathrm{1}/\mathrm{2}}\left){C}_{k}\right({z}_{j+\mathrm{1}/\mathrm{2}})$. This is noted to be an exact equation for the rate of change of the cellaveraged concentrations, ${\overline{C}}_{k,j}$, given in terms of the advective and diffusive fluxes through the cell faces, and of the cellaverage source term. The fact that the concentration within a cell is explicitly expressed as the sum of the fluxes into and out of the cell gives the FVM its massconserving properties (in the case when S=0).
Next, we make an approximation by assuming the concentration to be constant within each control volume such that the cell averages ${\stackrel{\mathrm{\u203e}}{C}}_{k,j}$ and ${\stackrel{\mathrm{\u203e}}{S}}_{k,j}$ may be approximated by the values located at the cell centres. The concentration and source term at z_{j} are denoted by C_{k,j} and S_{k,j} respectively. Choosing a constant grid spacing Δz, the z axis is discretised into ${N}_{z}=H/\mathrm{\Delta}z$ equally sized cells, where H is the length of the domain. For another example of a similar derivation, see Versteeg and Malalasekera (2007, pp. 243–246).
A1.2 Numerical approximation of fluxes
In Eq. (A2), the first two terms on the right represent the diffusive fluxes through the faces of cell j, and the next two terms represent the advective fluxes. The diffusive fluxes at the cell faces are approximated by a secondorder central difference of the two adjacent cells (see e.g. Versteeg and Malalasekera, 2007, p. 117) – that is, in the form
where the diffusivity K_{j+½} is determined on the cell faces either explicitly by an analytic expression, if available, or by interpolation in the case of a discrete diffusivity.
For the advective fluxes, however, linear numerical schemes of secondorder accuracy and higher are known to yield numerical oscillations for advectiondominated problems (i.e. $Pe=\frac{\leftv\right\mathrm{\Delta}z}{K}\gg \mathrm{1}$) for nonsmooth solutions (Hundsdorfer and Verwer, 2003, pp. 66–67, 118–119). To be able to handle advectiondominated cases, we implemented a more flexible secondorder numerical scheme for advection. A flux limiter approach was used, where we approximate the advective fluxes with the firstorder upwind scheme with a secondorder correction that depends on a limiter function (Hundsdorfer and Verwer, 2003, pp. 216–217; Versteeg and Malalasekera, 2007, pp. 165–171). The advective fluxes are thus approximated as
Here, positive and negative velocities are handled separately by letting ${w}_{k,j+\mathrm{1}/\mathrm{2}}^{+}=max(\mathrm{0},{w}_{k}({z}_{j+\mathrm{1}/\mathrm{2}}\left)\right)$ and ${w}_{k,j+\mathrm{1}/\mathrm{2}}^{}=min(\mathrm{0},{w}_{k}({z}_{j+\mathrm{1}/\mathrm{2}}\left)\right)$. The limiter function ψ(ρ^{±}) determines the correction, with ρ^{+} and ρ^{−} given by
The flux limiting was done with the UMIST limiter function (Lien and Leschziner, 1994; Versteeg and Malalasekera, 2007, pp. 170–178), given by
We see that ${\mathit{\rho}}_{j+\mathrm{1}/\mathrm{2}}^{\pm}$ is the ratio between the concentration gradients at the upstream or downstream sides and the gradient across the cell face at ${z}_{j+\mathrm{1}/\mathrm{2}}$ (Versteeg and Malalasekera, 2007, pp. 167, 171–172; Hundsdorfer and Verwer, 2003, p. 216). If ρ is close to 1, which will be the case for reasonably smooth concentration profiles, then ψ(ρ) is also close to 1, in which case Eq. (A4) is approximately equal to a regular secondorder central finite difference. On the other hand, if there are large differences in concentration gradients between neighbouring pairs of cells, then ψ(ρ) will be close to either 0 or 2, in which case Eq. (A4) will be approximately equal to a firstorder upwind or downwind scheme to avoid numerical oscillations. It is also reduced to upwind near the boundaries. This approach was found to result in secondorder accuracy in space for the cases considered here (see Appendix C for examples of grid refinement tests).
A1.3 Boundary conditions
As described in Sect. 2.4, we consider two different types of boundary conditions, namely zero flux (Eq. 6) and zero diffusive flux (Eq. 7). In our Eulerian implementation, which uses a finitevolume discretisation scheme, it is straightforward to use different boundary conditions by simply setting the required fluxes to zero.
First, we consider the cell adjacent to the surface boundary. Note that the uppermost cell is indexed j=0 and that the cell centre is at ${z}_{\mathrm{0}}=\mathrm{\Delta}z/\mathrm{2}$, while the surface boundary is found at ${z}_{\mathrm{1}/\mathrm{2}}=\mathrm{0}$ (see also Fig. A1). The insertion of the zerodiffusiveflux surface boundary condition of Eq. (7) into Eq. (A2) yields
We note that setting $(K\frac{\partial {C}_{k}}{\partial z}{)}_{\frac{\mathrm{1}}{\mathrm{2}}}=\mathrm{0}$ implicitly means that we are assuming zero concentration gradient across the surface since K(z)>0 everywhere.
To implement a zeroflux boundary condition, we force both the diffusive and advective fluxes to be zero. This we can do by setting the diffusive flux to zero as above and additionally setting ${w}_{k}\left({z}_{\mathrm{1}/\mathrm{2}}\right)=\mathrm{0}$. We note that, in our model, the advection velocity represents the rise or settling velocity of the suspended matter; hence, this is equivalent to introducing a spatially variable velocity, saying that a positively buoyant particle has a rise velocity of zero if it is at the surface:
The same reasoning applies to the boundary conditions at the seafloor.
A1.4 Numerically solving the discretised equation
The scheme was discretised in time with a fixed time step, such that
and it was integrated by the Crank–Nicolson method (Gustafsson, 2008, p. 39), given by
where ${C}_{k,j}^{i}$ is the concentration of class k in cell j at time t_{i}, and ${F}_{k,j}^{i}$ is the righthand side of Eq. (A2) evaluated at time t_{i}. The Crank–Nicolson scheme was chosen for its secondorder accuracy in time and its favourable stability (Versteeg and Malalasekera, 2007, pp. 247–248).
With the chosen discretisation schemes in space and time, our numerical solver can be rewritten into a matrix equation in the form
where L and R are matrices to be described below. The vectors C^{i} and S^{i} contain, respectively, the concentration in each cell for each class and the source term for each cell and each class at time t_{i} (see schematic illustration in Fig. A2).
The matrices L and R can each be written as sums of two terms:
Here, L_{AD} and R_{AD} contain the coefficients that implement the finitedifference approximations of the spatial derivatives of the concentration, specifically the central finite difference for diffusion and the upwind for advection. L_{FL} and R_{FL} contain the coefficients for the flux limiter correction to the advection scheme.
In all the cases we consider below, the matrices describing the advection and diffusion terms remain constant in time. The flux limiter correction and the reaction terms, on the other hand, are themselves functions of the concentration C. Since this leads to a nonlinear system of equations, the equation system given by Eq. (A11) cannot be solved directly, as in the linear case, and thus we have adopted an iterative scheme.
Writing out Eq. (A11) with the dependence on concentration shown explicitly, we get
This equation must be solved at every time step to find the concentration C^{i+1} at time t_{i+1} from the current concentration C^{i} at time t_{i}, taking the reaction terms S^{i} and S^{i+1} into account. Since C^{i+1} is of course not yet known at time t_{i}, we start by using the current concentration as an initial guess at the next concentration, ${\stackrel{\mathrm{\u0303}}{\mathit{C}}}^{i+\mathrm{1}}={\mathit{C}}^{i}$, and then we solve the following system:
We then refine our guess by letting ${\stackrel{\mathrm{\u0303}}{\mathit{C}}}^{i+\mathrm{1}}\to {\stackrel{\mathrm{\u0303}}{\stackrel{\mathrm{\u0303}}{\mathit{C}}}}^{i+\mathrm{1}}$, we solve again, and then we repeat. At each iteration, we calculate a measure of the error,
and we terminate the iterative procedure when $\mathrm{err}<\mathit{\eta}\cdot {\mathrm{err}}_{\mathrm{0}}$, where err_{0} is the error calculated from the first guess, for some tolerance η.
For the cases we consider, L and R are tridiagonal matrices, which means that Eq. (A14) can be solved efficiently using the tridiagonal matrix algorithm (TDMA; see, for example, Press et al., 2007, pp. 56–57).
A1.5 Solving for multiple classes
Equation (A1) describes the advection and diffusion of a concentration field C_{k}(z,t), where the advection velocity w_{k} represents the rising or settling speed of particles of class k. In this paper, we wish to consider a distribution of particles with different rising or settling speeds, represented by discrete classes.
If there are no reaction terms, there is no interaction between the different classes of particles; hence, the advection–diffusion equation can be solved independently for each class, as described in the previous section. Nevertheless, for reasons of flexibility, it may be convenient to implement an approach where the advection–diffusion equation is solved simultaneously for all classes since this allows interacting reaction terms to be added.
When there are reaction terms that allow the classes to interact, i.e. if the reaction term S_{k} of class k depends on the concentration of any other class, then we have to solve the advection–diffusion–reaction equation simultaneously for all classes.
Our chosen approach is first to set up the equation system for each individual class. The concentration of class k in the different spatial cells is then represented by a vector C_{k}, and the matrices in Eq. (A11) for class k we call L_{k} and R_{k}. We then create the full matrices L and R as blockdiagonal matrices, with the blocks being made up of the matrices ${\mathbf{L}}_{\mathrm{1}},\mathrm{\dots},{\mathbf{L}}_{{N}_{k}}$ and ${\mathbf{R}}_{\mathrm{1}},\mathrm{\dots},{\mathbf{R}}_{{N}_{k}}$. We next combine the concentrations for all the classes into an overall vector $\mathit{C}=[{\mathit{C}}_{\mathrm{1}},\mathrm{\dots},{\mathit{C}}_{{N}_{k}}]$, where N_{k} is the total number of classes. This leads to a linear system of equations, schematically illustrated in Fig. A2, which can be solved to find the concentrations in all cells for all classes.
As a practical matter, we note that since the matrices L and R are tridiagonal, using a sparse array class in the implementation means that the full matrices do not require any additional memory compared to storing the matrices ${\mathbf{L}}_{\mathrm{1}},\mathrm{\dots},{\mathbf{L}}_{{N}_{k}}$ and ${\mathbf{R}}_{\mathrm{1}},\mathrm{\dots},{\mathbf{R}}_{{N}_{k}}$ separately.
To obtain a distribution of rising and settling velocities for microplastic particles, we combine the distributions for density, size, and shape presented by Kooi and Koelmans (2019) with the empirical relations for terminal velocity obtained by Waldschläger and Schüttrumpf (2019). In the following, we describe the steps leading to a full distribution of rising and settling velocities. The Python code used to obtain the distribution is available on GitHub, both as a script and as a Jupyter Notebook (https://github.com/nordam/EulerianandLagrangianmethods, last access: 11 September 2023).
We employ a Monte Carlo approach to obtain a velocity distribution, where we draw a sample from each of the distributions describing density, size, and shape, which we can then map to a terminal velocity via the relations of Waldschläger and Schüttrumpf (2019). By repeating this process a large number of times, we obtain a distribution of samples of rising and settling speeds which are consistent with the underlying distributions of density, size, and shape. In the Lagrangian implementation, we can use the sampled velocities directly, while in the Eulerian implementation, we draw a very large number (10^{9}) of samples and map them to discrete classes by creating a histogram with the required number of bins.
B1 Density distribution
Kooi and Koelmans (2019) found density to be accurately described
by a normal inverse Gaussian distribution, with the parameters μ=0.84,
δ=0.097, α=75.1, and β=71.3. Generating random samples
from this distribution is straightforward using, for example, the class
scipy.stats.norminvgauss
from the SciPy library.
B2 Size distribution
The size distribution found by Kooi and Koelmans (2019) is a truncated powerlaw distribution for the number of particles, where the probability number density for particle size s in micrometres is
with α=1.6, where s is the size in micrometres and where the range of the distribution is from 20 µm to 5 mm. P_{0} is a normalisation factor given by
To draw random samples from the distribution with the probability density function (PDF) given by Eq. (B1), we use the inversion method (see e.g. Devroye, 1986, pp. 27–28). We let
be the cumulative distribution function (CDF) of the size distribution. Then, random samples from this distribution can be generated by evaluating F^{−1}(U), where F^{−1} is the inverse of F(s) (which can be found analytically in this case), and $U\in [\mathrm{0},\mathrm{1}]$ is a uniformly distributed random number.
B3 Shape distribution
Kooi and Koelmans (2019) provide a shape distribution for microplastics in terms of the Corey shape factor (CSF), which is defined as
where L, W, and H are respectively the length, width, and height of a particle. We assume that $L\ge W\ge H$ and that the values are scaled by L such that L=1 and $W,H\in [\mathrm{0},\mathrm{1}]$. While Kooi and Koelmans (2019) obtain a single distribution for the CSF in the end, the underlying data are presented as distributions for the width and height for different shape categories, with different relative abundances: fibres (48.5 %), fragments (31 %), beads (6.5 %), films (5.5 %), and foam (3.5 %). Note that these relative abundances only add up to 95 %; hence, we have normalised them by dividing each one by 0.95, and we thus achieve the probabilities listed in Table B1.
For each of the shape categories, the distributions for width and height are
assumed to be symmetric, triangular distributions, as described in
Kooi and Koelmans (2019), and the distribution parameters for each shape are
given in Table S2 of the supplementary materials of Kooi and Koelmans (2019).
We also list the distribution parameters in Table B1 for
completeness. Generating random samples from a triangular distribution is
straightforward with, for example, scipy.stats.triang
from the SciPy
library.
Based on the above, our algorithm for the generation of a random CSF for a microplastic particle is as follows:
B4 Terminal velocities
Waldschläger and Schüttrumpf (2019) provide relations for the rising and settling speeds of microplastic particles in terms of their density, size, and shape. Specifically, they give the terminal velocity as
where d_{eq} is the equivalent diameter (i.e. the diameter of a sphere with the same volume as the particle), ρ_{s} and ρ are respectively the densities of the particle and the ambient water, g is the acceleration of gravity, and C_{D} is the drag coefficient.
Four different expressions are given for the drag coefficient, distinguishing between particles that are lighter or denser than water, as well as distinguishing fibres from other shapes (pellets and fragments). For sinking nonfibre particles, the expression is
For rising nonfibre particles, the expression is
For sinking fibres, the expression is
For rising fibres, the expression is
The drag coefficients are given in terms of the CSF, the Reynolds number, and in one case, the Powers' roundness, with the Reynolds number given by
where ν is the kinematic viscosity of the ambient fluid. Given that the drag coefficients depend on the Reynolds number, which in turn depends on the velocity, Eq. (B5) was solved numerically by means of an iterative approach to find the velocity.
The Powers' roundness P is a set of six classes from 1 (very angular) to 6 (well rounded), where particles are assigned a class by comparison with a published photograph of reference particles (Powers, 1953). Powers' roundness is not addressed by Kooi and Koelmans (2019), but Waldschläger and Schüttrumpf (2019, Fig. S6, Supplement) show a histogram of the Powers' roundness distribution. Based on this, we have chosen to use a uniformly distributed random integer between 1 and 6 to represent Powers' roundness. Note that the Powers' roundness only appears in the expression for the drag coefficient for rising pellets and fragments and is thus not used for fibres.
Even though Waldschläger and Schüttrumpf (2019) provide separate drag coefficients for rising and sinking particles, we have chosen to use the same coefficients in both cases. Our reasoning is twofold: first, from a physical point of view, it is expected that two particles with density ±Δρ relative to water should have terminal velocities with opposite directions but equal magnitude (certainly as long as $\left\mathrm{\Delta}\mathit{\rho}\right\ll \mathit{\rho}$). Second, the drag coefficient of sinking nonfibre particles (pellets and fragments) provided by Waldschläger and Schüttrumpf (2019) scales as $R{e}^{\mathrm{1}/\mathrm{3}}$, while Stokes' law is
For small particles, in the lower part of the microplastic size range^{2}, the sinking velocities predicted with Eq. (B6) can be 1–2 orders of magnitude higher than those predicted by Stokes' law (even with CSF=1, as would be the case for spherical particles). The drag coefficient provided for rising nonfibre particles, on the other hand, has the same scaling as Stokes' law for small Reynolds numbers. Hence, we have chosen to use this drag coefficient for both rising and sinking particles. In the case of fibres, Eqs. (B9) and (B8) give very similar velocities, but for reasons of symmetry, we have chosen to use Eq. (B9) for both rising and sinking fibres.
B5 Terminal velocity distribution
Based on the description above, we generate random realisations of particle properties, where each particle has the properties shown in Table B2. From these properties, we can calculate terminal velocities using the expressions found by Waldschläger and Schüttrumpf (2019), as described above.
For the Lagrangian implementation, N_{p} random velocities were generated at the start of each simulation. For the Eulerian implementation, 10^{9} random velocities were generated and mapped to N_{k} different classes. Given the wide range of velocities, we chose to use equally spaced velocity bins on a logscale, with positive and negative velocities treated separately. For positive (rising) velocities, we used bins ranging from 10^{−8} to 0.1 m s^{−1}, and for the negative velocities, we used bins ranging from −0.3 to ${\mathrm{10}}^{\mathrm{8}}$ m s^{−1}. We used the same number of classes to represent the positive and the negative velocities. By always choosing a total number of classes N_{k}, which is divisible by 2, we get a consistent and automated way to divide the velocity distribution into discrete classes. Note that the approach described here was chosen for simplicity and consistency as our goal is to investigate the numerical convergence with the number of classes.
To get the discrete velocity classes, we generated 10^{9} random velocities and mapped them onto the bins described below, with the velocity of each bin being represented by its midpoint on the logarithmic scale, as described in Sect. 3.1.2. None of the randomly generated velocities were smaller than −0.3 m s^{−1} or greater than 0.1 m s^{−1}. Approximately 0.02 % of the generated velocities fell between ${\mathrm{10}}^{\mathrm{8}}$ and 10^{−8} m s^{−1}. These were ignored. The bin counts were then converted to normalised mass fractions, such that the total mass sums to 1. Generating 10^{9} random terminal velocities in this manner took approximately 15 min on a fairly standard PC.
An example histogram showing the distribution of 10^{9} terminal velocities is shown in Fig. B1 using the same 256 bins that were used for the most accurate Eulerian velocity discretisation. Note that this is a number distribution, giving the number of particles with a particular speed. In many cases, a volume (or mass) distribution may be more useful as this would give the volume of microplastics present as particles with a particular speed. However, obtaining a joint volume probability distribution for size and shape from the number distributions presented by Kooi and Koelmans (2019) is outside the scope of this study.
Here, we present the convergence of the first moment of the distribution with respect to the time step for the Lagrangian implementation and with respect to the time step and grid cell size for the Eulerian implementation. These results were used to choose suitable numerical parameters for the investigation of convergence with the number of classes and number of particles. They also serve as a sanity check to demonstrate that the numerical schemes behave as expected.
The results are shown in Fig. C1 (case 1), Fig. C2 (case 2), and Fig. C3 (case 3). In each figure, panel (a) shows the convergence of the first moment with the time step in the Lagrangian implementation. For each time step, the first moment was calculated for each value of the time step Δt as the average particle position over 100 simulations, each with 3×10^{6} particles. The error was then found by using the average particle position across 100 simulations with 3×10^{6} particles using a shorter time step of Δt=2 s. The convergence appears to be of the first order, which is as expected for the Euler–Maruyama method as this method has an order of convergence of 1 in the weak sense.
We also note that the error in the first moment of the particle distribution for a Lagrangian stochastic method has two terms: a discretisation term due to the time step and a stochastic term due to the finite number of samples (Pavliotis, 2014, p. 151). The expected magnitude of the stochastic error term can be estimated from the Lindeberg–Lévy theorem, as described in Sect. 2.6. In panel (a) of Figs. C1, C2, and C3, the magnitude of the stochastic sample error is shown as coloured dashed lines, where we have used ${N}_{p}=\mathrm{3}\times {\mathrm{10}}^{\mathrm{8}}$ samples, and where the variance of the particle positions for each of the times considered have been assumed to be a good approximation for the (unknown) true variance.
For the Eulerian implementation, the convergence with time step is shown in panel (b) of Figs. C1, C2, and C3, where the number of classes and cells have been kept fixed at N_{k}=128 and N_{z}=8000. Finally, the convergence with spatial discretisation for the Eulerian method is shown in panel (c) of Figs. C1, C2, and C3. Here, the time step and number of classes have been kept fixed at Δt=10 s and N_{k}=128. We observe that convergence is of the second order in both space and time, as expected.
The code used to run the simulations and to create the plots shown in this paper is available under a GPL3.0 license and can be found at https://github.com/nordam/EulerianandLagrangianmethods (last access: 11 September 2023) or at https://doi.org/10.5281/zenodo.8164812 (Nordam, 2023). No additional data have been used.
TN wrote the Lagrangian implementation, ran all simulations, and wrote the first draft of the paper. RK wrote the Eulerian implementation. All the authors contributed significantly to the review and editing of the paper.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors would like to thank Svein Sundby for the helpful advice on the velocity distribution of pelagic fish eggs.
The work of Tor Nordam was supported in part by the Norwegian Research Council project INDORSE (grant no. 267793) and in part by the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 814426 – NanoInformaTIX. The work of Erik van Sebille was supported through funding from the Netherlands Organization for Scientific Research (NWO), Earth and Life Sciences (grant no. OCENW.KLEIN.085) and through funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 715386).
This paper was edited by Sylwester Arabas and reviewed by Alethea Mountford and one anonymous referee.
Abuhegazy, M., Talaat, K., Anderoglu, O., and Poroseva, S. V.: Numerical investigation of aerosol transport in a classroom with relevance to COVID19, Phys. Fluids, 32, 103311, https://doi.org/10.1063/5.0029118, 2020. a
Benson, D. A., Aquino, T., Bolster, D., Engdahl, N., Henri, C. V., and FernandezGarcia, D.: A comparison of Eulerian and Lagrangian transport and nonlinear reaction algorithms, Adv. Water Resour., 99, 15–37, https://doi.org/10.1016/j.advwatres.2016.11.003, 2017. a
Billingsley, P.: Probability and measure, John Wiley & Sons, New York Chichester Brisbane Toronto, ISBN 0471031739, 1979. a, b
Craig, P. D. and Banner, M. L.: Modeling waveenhanced turbulence in the ocean surface layer, J. Phys. Oceanogr., 24, 2546–2559, https://doi.org/10.1175/15200485(1994)024<2546:MWETIT>2.0.CO;2, 1994. a
Cui, F., Boufadel, M. C., Geng, X., Gao, F., Zhao, L., King, T., and Lee, K.: Oil Droplets Transport Under a DeepWater Plunging Breaker: Impact of Droplet Inertia, J. Geophys. Res.Oceans, 123, 9082–9100, https://doi.org/10.1029/2018JC014495, 2018. a
Cui, F., Zhao, L., Daskiran, C., King, T., Lee, K., Katz, J., and Boufadel, M. C.: Modeling oil dispersion under breaking waves. Part II: Coupling Lagrangian particle tracking with population balance model, Environ. Fluid Mech., 20, 1553–1578, https://doi.org/10.1007/s10652020097591, 2020. a, b
de la Fuente, R., Drótos, G., HernándezGarcía, E., López, C., and van Sebille, E.: Sinking microplastics in the water column: simulations in the Mediterranean Sea, Ocean Sci., 17, 431–453, https://doi.org/10.5194/os174312021, 2021. a
Delandmeter, P. and van Sebille, E.: The Parcels v2.0 Lagrangian framework: new field interpolation schemes, Geosci. Model Dev., 12, 3571–3584, https://doi.org/10.5194/gmd1235712019, 2019. a
Delvigne, G. and Sweeney, C.: Natural dispersion of oil, Oil Chem. Pollut., 4, 281–310, 1988. a, b
Devroye, L.: Nonuniform random variate generation, SpringerVerlag, New York, https://doi.org/10.1007/9781461386438, 1986. a
Dissanayake, A. L., DeGraff, J. A., Yapa, P. D., Nakata, K., Ishihara, Y., and Yabe, I.: Modeling the impact of CO_{2} releases in Kagoshima Bay, Japan, J. HydroEnviron. Res., 6, 195–208, https://doi.org/10.1016/j.jher.2012.02.001, 2012. a
Dissanayake, A. L., Yapa, P. D., and Nakata, K.: Simulation of hydrothermal vents in the Izena Cauldron, Mid Okinawa trough, Japan and other Pacific locations, J. HydroEnviron. Res., 8, 343–357, https://doi.org/10.1016/j.jher.2014.05.003, 2014. a
Dugstad, J. S., Koszalka, I. M., Isachsen, P. E., Dagestad, K.F., and Fer, I.: Vertical Structure and Seasonal Variability of the Inflow to the Lofoten Basin Inferred From HighResolution Lagrangian Simulations, J. Geophys. Res.Oceans, 124, 9384–9403, https://doi.org/10.1029/2019JC015474, 2019. a
Elliott, A., Hurford, N., and Penn, C.: Shear diffusion and the spreading of oil slicks, Mar. Pollut. Bull., 17, 308–313, 1986. a
Fay, B., Glaab, H., Jacobsen, I., and Schrodin, R.: Evaluation of Eulerian and Lagrangian atmospheric transport models at the Deutscher Wetterdienst using ANATEX surface tracer data, Atmos. Environ., 29, 2485–2497, https://doi.org/10.1016/13522310(95)00144N, 1995. a
Fischer, R., Lobelle, D., Kooi, M., Koelmans, A., Onink, V., Laufkötter, C., AmaralZettler, L., Yool, A., and van Sebille, E.: Modelling submerged biofouled microplastics and their vertical trajectories, Biogeosciences, 19, 2211–2234, https://doi.org/10.5194/bg1922112022, 2022. a
Gill, A. E.: AtmosphereOcean Dynamics, vol. 30 of International Geophysics Series, Academic Press, New York, ISBN 9780122835223, 1982. a
Gillespie, D. T. and Seitaridou, E.: Simple Brownian Diffusion, Oxford University Press, Oxford, ISBN 0199664501, 2012. a, b
Gräwe, U.: Implementation of highorder particletracking schemes in a water column model, Ocean Model., 36, 80–89, https://doi.org/10.1016/j.ocemod.2010.10.002, 2011. a
Gräwe, U., Deleersnijder, E., Shah, S. H. A. M., and Heemink, A. W.: Why the Euler scheme in particle tracking is not enough: the shallowsea pycnocline test case, Ocean Dynam., 62, 501–514, https://doi.org/10.1007/s102360120523y, 2012. a
Gustafsson, B.: High order difference methods for time dependent PDE, SpringerVerlag, Berlin Heidelberg, https://doi.org/10.1007/9783540749936, 2008. a, b
Hundsdorfer, W. and Verwer, J. G.: Numerical solution of timedependent advectiondiffusionreaction equations, SpringerVerlag, Berlin, Heidelberg, https://doi.org/10.1007/9783662090176, 2003. a, b, c, d, e, f
Isachenko, I.: Catching the variety: Obtaining the distribution of terminal velocities of microplastics particles in a stagnant fluid by a stochastic simulation, Mar. Pollut. Bull., 159, 111464, https://doi.org/10.1016/j.marpolbul.2020.111464, 2020. a
Israelsson, P. H., Do Kim, Y., and Adams, E. E.: A comparison of three Lagrangian approaches for extending near field mixing calculations, Environ. Modell. Softw., 21, 1631–1649, https://doi.org/10.1016/j.envsoft.2005.07.008, 2006. a
Johansen, Ø., Reed, M., and Bodsberg, N. R.: Natural dispersion revisited, Mar. Pollut. Bull., 93, 20–26, https://doi.org/10.1016/j.marpolbul.2015.02.026, 2015. a, b
Kloeden, P. E. and Platen, E.: Numerical Solution of Stochastic Differential Equations, SpringerVerlag Berlin Heidelberg, https://doi.org/10.1007/9783662126165, 1992. a, b, c, d, e, f, g
Kooi, M. and Koelmans, A. A.: Simplifying microplastic via continuous probability distributions for size, shape, and density, Environ. Sci. Tech. Let., 6, 551–557, https://doi.org/10.1021/acs.estlett.9b00379, 2019. a, b, c, d, e, f, g, h, i, j, k, l
Li, Z., Spaulding, M. L., and FrenchMcCay, D.: An algorithm for modeling entrainment and naturally and chemically dispersed oil droplet size distribution under surface breaking wave conditions, Mar. Pollut. Bull., 119, 145–152, 2017. a
Lien, F.S. and Leschziner, M.: Upstream monotonic interpolation for scalar transport with application to complex turbulent flows, Int. J. Numer. Meth. Fl., 19, 527–548, https://doi.org/10.1002/fld.1650190606, 1994. a
Longest, P. W. and Xi, J.: Effectiveness of direct Lagrangian tracking models for simulating nanoparticle deposition in the upper airways, Aerosol Sci. Tech., 41, 380–397, https://doi.org/10.1080/02786820701203223, 2007. a
Lynch, D. R., Greenberg, D. A., Bilgili, A., McGillicuddy Jr, D. J., Manning, J. P., and Aretxabaleta, A. L.: Particles in the Coastal Ocean: Theory and Applications, Cambridge University Press, ISBN 9781107061750, 2014. a
Madec, G., BourdalléBadie, R., Chanut, J., Clementi, E., Coward, A., Ethé, C., Iovino, D., Lea, D., Lévy, C., Lovato, T., Martin, N., Masson, S., Mocavero, S., Rousset, C., Storkey, D., Müeller, S., Nurser, G., Bell, M., Samson, G., Mathiot, P., Mele, F., and Moulin, A.: NEMO ocean engine, vol. 27 of Scientific Notes of IPSL Climate Modelling Center, Zenodo, https://doi.org/10.5281/zenodo.6334656, 2022. a
Marsh, R., Ivchenko, V. O., Skliris, N., Alderson, S., Bigg, G. R., Madec, G., Blaker, A. T., Aksenov, Y., Sinha, B., Coward, A. C., Le Sommer, J., Merino, N., and Zalesny, V. B.: NEMO–ICB (v1.0): interactive icebergs in the NEMO ocean model globally configured at eddypermitting resolution, Geosci. Model Dev., 8, 1547–1562, https://doi.org/10.5194/gmd815472015, 2015. a
Matsumura, Y. and Ohshima, K. I.: Lagrangian modelling of frazil ice in the ocean, Ann. Glaciol., 56, 373–382, https://doi.org/10.3189/2015AoG69A657, 2015. a
Mountford, A. and Morales Maqueda, M.: Eulerian Modeling of the ThreeDimensional Distribution of Seven Popular Microplastic Types in the Global Ocean, J. Geophys. Res.Oceans, 124, 8558–8573, https://doi.org/10.1029/2019JC015050, 2019. a
Nepstad, R., Liste, M., Alver, M. O., Nordam, T., Davies, E., and Glette, T.: Highresolution numerical modelling of a marine mine tailings discharge in Western Norway, Regional Studies in Marine Science, 39, 101404, https://doi.org/10.1016/j.rsma.2020.101404, 2020. a
Nepstad, R., Nordam, T., Ellingsen, I. H., Eisenhauer, L., Litzler, E., and Kotzakoulakis, K.: Impact of flow field resolution on produced water transport in Lagrangian and Eulerian models, Mar. Pollut. Bull., 182, 113928, https://doi.org/10.1016/j.marpolbul.2022.113928, 2022. a, b
Nooteboom, P. D., Bijl, P. K., van Sebille, E., von der Heydt, A. S., and Dijkstra, H. A.: Transport Bias by Ocean Currents in Sedimentary Microplankton Assemblages: Implications for Paleoceanographic Reconstructions, Paleoceanogr. Paleoclim., 34, 1178–1194, https://doi.org/10.1029/2019PA003606, 2019. a
Nordam, T.: nordam/EulerianandLagrangianmethods: 1.2, Zenodo [code], https://doi.org/10.5281/zenodo.8164812, 2023. a
Nordam, T., Kristiansen, R., Nepstad, R., and Röhrs, J.: Numerical analysis of boundary conditions in a Lagrangian particle model for vertical mixing, transport and surfacing of buoyant particles in the water column, Ocean Model., 136, 107–119, https://doi.org/10.1016/j.ocemod.2019.03.003, 2019a. a, b, c, d, e, f, g
Nordam, T., Nepstad, R., Litzler, E., and Röhrs, J.: On the use of random walk schemes in oil spill modelling, Mar. Pollut. Bull., 146, 631–638, https://doi.org/10.1016/j.marpolbul.2019.07.002, 2019b. a, b, c, d
Pavliotis, G. A.: Stochastic Processes and Applications, Springer Science+Business Media, New York, ISBN 9781493913220, 2014. a
Powers, M. C.: A new roundness scale for sedimentary particles, J. Sediment. Res., 23, 117–119, https://doi.org/10.1306/D42695672B2611D78648000102C1865D, 1953. a
Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: Numerical Recipes, Cambridge university press, New York, 3rd Edn., ISBN 9780521880688, 2007. a, b
Rodean, H. C.: Stochastic Lagrangian Models of Turbulent Diffusion, vol. 26 of Meteorological Monographs, American Meteorological Society, Boston, MA, ISBN 1878220233, 1996. a
Röhrs, J., Christensen, K. H., Vikebø, F., Sundby, S., Saetra, Ø., and Broström, G.: Waveinduced transport and vertical mixing of pelagic eggs and larvae, Limnol. Oceanogr., 59, 1213–1227, https://doi.org/10.4319/lo.2014.59.4.1213, 2014. a, b
Röhrs, J., Dagestad, K.F., Asbjørnsen, H., Nordam, T., Skancke, J., Jones, C. E., and Brekke, C.: The effect of vertical mixing on the horizontal drift of oil spills, Ocean Sci., 14, 1581–1601, https://doi.org/10.5194/os1415812018, 2018. a
Ross, O. N. and Sharples, J.: Recipe for 1D Lagrangian particle tracking models in spacevarying diffusivity, Limnol. Oceanogr.Meth., 2, 289–302, https://doi.org/10.4319/lom.2004.2.289, 2004. a
Rowe, M., Anderson, E., Wynne, T., Stumpf, R., Fanslow, D., Kijanka, K., Vanderploeg, H., Strickler, J., and Davis, T.: Vertical distribution of buoyant Microcystis blooms in a Lagrangian particle tracking model for shortterm forecasts in Lake Erie, J. Geophys. Res.Oceans, 121, 5296–5314, https://doi.org/10.1002/2016JC011720, 2016. a
Saharia, A. M., Zhu, Z., Aich, N., Baalousha, M., and Atkinson, J. F.: Modeling the transport of titanium dioxide nanomaterials from combined sewer overflows in an urban river, Sci. Total Environ., 696, 133904, https://doi.org/10.1016/j.scitotenv.2019.133904, 2019. a
Scutt Phillips, J., Gupta, A. S., Senina, I., van Sebille, E., Lange, M., Lehodey, P., Hampton, J., and Nicol, S.: An individualbased model of skipjack tuna (Katsuwonus pelamis) movement in the tropical Pacific ocean, Prog. Oceanogr., 164, 63–74, https://doi.org/10.1016/j.pocean.2018.04.007, 2018. a
Shchepetkin, A. F. and McWilliams, J. C.: The regional oceanic modeling system (ROMS): a splitexplicit, freesurface, topographyfollowingcoordinate oceanic model, Ocean Model., 9, 347–404, https://doi.org/10.1016/j.ocemod.2004.08.002, 2005. a
Skorokhod, A. V.: Stochastic equations for diffusion processes in a bounded region, Theor. Probab. Appl+., 6, 264–274, https://doi.org/10.1137/1106035, 1961. a, b
Skorokhod, A. V.: Stochastic equations for diffusion processes in a bounded region II, Theor. Probab. Appl+., 7, 3–23, https://doi.org/10.1137/1107002, 1962. a
Strauss, R. D. T. and Effenberger, F.: A hitchhiker’s guide to stochastic differential equations, Space Sci. Rev., 212, 151–192, https://doi.org/10.1007/s112140170351y, 2017. a
Sundby, S.: A onedimensional model for the vertical distribution of pelagic fish eggs in the mixed layer, DeepSea Res. Pt. I, 30, 645–661, https://doi.org/10.1016/01980149(83)900420, 1983. a, b, c, d, e
Sundby, S. and Kristiansen, T.: The principles of buoyancy in marine fish eggs and their vertical distributions across the world oceans, PLoS One, 10, e0138821, https://doi.org/10.1371/journal.pone.0138821, 2015. a
Tayfun, M. A. and Wang, H.: Monte Carlo simulation of oil slick movements, J. Waterway Div.ASCE, 99, 309–324, https://doi.org/10.1061/AWHCAR.0000197, 1973. a
Tegen, I. and Lacis, A. A.: Modeling of particle size distribution and its influence on the radiative properties of mineral dust aerosol, J. Geophys. Res.Atmos., 101, 19237–19244, https://doi.org/10.1029/95JD03610, 1996. a
Thorpe, S. A.: The turbulent ocean, Cambridge University Press, Cambridge, UK, https://doi.org/10.1017/CBO9780511819933, 2005. a
Tkalich, P. and Chan, E. S.: Vertical mixing of oil droplets by breaking waves, Mar. Pollut. Bull., 44, 1219–1229, https://doi.org/10.1016/S0025326X(02)001789, 2002. a
Umlauf, L., Burchard, H., and Bolding, K.: GOTM – Scientific Documentation: version 3.2, Marine Science Reports, LeibnizInstitute for Baltic Sea Research, Warnemuende, Germany, https://www.iowarnemuende.de/files/forschung/meereswissenschaftlicheberichte/mebe63_2005gotm.pdf (last access: 11 September 2023), 2005. a
Van Sebille, E., Scussolini, P., Durgadoo, J. V., Peeters, F. J., Biastoch, A., Weijer, W., Turney, C., Paris, C. B., and Zahn, R.: Ocean currents generate large footprints in marine palaeoclimate proxies, Nat. Commun., 6, 1–8, https://doi.org/10.1038/ncomms7521, 2015. a
van Sebille, E., Griffies, S. M., Abernathey, R., Adams, T. P., Berloff, P., Biastoch, A., Blanke, B., Chassignet, E. P., Cheng, Y., Cotter, C. J., Deleersnijder, E., Döös, K., Drake, H. F., Drijfhout, S., Gary, S. F., Heemink, A. W., Kjellsson, J., Koszalka, I. M., Lange, M., Lique, C., MacGilchrist, G. A., Marsh, R., Mayorga Adame, C. G., McAdam, R., Nencioli, F., Paris, C. B., Piggott, M. D., Polton, J. A., Rühs, S., Shah, S. H., Thomas, M. D., Wang, J., Wolfram, P. J., Zanna, L., and Zika, J. D.: Lagrangian ocean analysis: Fundamentals and practices, Ocean Model., 121, 49–75, https://doi.org/10.1016/j.ocemod.2017.11.008, 2018. a
Versteeg, H. K. and Malalasekera, W.: An introduction to computational fluid dynamics: the finite volume method, Pearson/Prentice Hall, Harlow, 2nd Edn., ISBN 9780131274983, 2007. a, b, c, d, e, f, g, h
Waldschläger, K. and Schüttrumpf, H.: Effects of particle properties on the settling and rise velocities of microplastics in freshwater under laboratory conditions, Environ. Sci. Technol., 53, 1958–1966, https://doi.org/10.1021/acs.est.8b06794, 2019. a, b, c, d, e, f, g, h, i, j
Warner, J. C., Sherwood, C. R., Signell, R. P., Harris, C. K., and Arango, H. G.: Development of a threedimensional, regional, coupled wave, current, and sedimenttransport model, Comput. Geosci., 34, 1284–1306, https://doi.org/10.1016/j.cageo.2008.02.012, 2008. a
Wendroff, B.: First Principles of Numerical Analysis: An Undergraduate Text, AddisonWesley Publishing Company, Reading, Menlo Park, London, Don Mills, ISBN 9781114499157, 1969. a, b
Wichmann, D., Delandmeter, P., and van Sebille, E.: Influence of nearsurface currents on the global dispersal of marine microplastic, J. Geophys. Res.Oceans, 124, 6086–6096, https://doi.org/10.1029/2019JC015328, 2019. a, b
Wilson, J. D. and Flesch, T. K.: Flow boundaries in randomflight dispersion models: enforcing the wellmixed condition, J. Appl. Meteorol. Clim., 32, 1695–1707, https://doi.org/10.1175/15200450(1993)032<1695:FBIRFD>2.0.CO;2, 1993. a
Wimalaratne, M. R., Yapa, P. D., Nakata, K., and Premathilake, L. T.: Transport of dissolved gas and its ecological impact after a gas release from deepwater, Mar. Pollut. Bull., 100, 279–288, https://doi.org/10.1016/j.marpolbul.2015.08.039, 2015. a
Young, E. F., Belchier, M., Hauser, L., Horsburgh, G. J., Meredith, M. P., Murphy, E. J., Pascoal, S., Rock, J., Tysklind, N., and Carvalho, G. R.: Oceanography and life history predict contrasting genetic population structure in two Antarctic fish species, Evol. Appl., 8, 486–509, https://doi.org/10.1111/eva.12259, 2015. a
Zender, C. S., Bian, H., and Newman, D.: Mineral Dust Entrainment and Deposition (DEAD) model: Description and 1990s dust climatology, J. Geophys. Res.Atmos., 108, 4416, https://doi.org/10.1029/2002JD002775, 2003. a
More than 11 000 citations according to Google Scholar (as of June 2023).
Note that the expressions given by Waldschläger and Schüttrumpf (2019) were found by fitting to experimentally measured terminal velocities, but as the smallest sinking nonfibre particle in their experiments was 1 mm, the expression is not necessarily valid for much smaller particles.
 Abstract
 Introduction
 Components of the advection–diffusion reaction model
 Implementation
 Case studies
 Discussion
 Conclusions
 Appendix A: Additional implementation details
 Appendix B: Distribution of terminal velocities for microplastics
 Appendix C: Additional convergence results
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Components of the advection–diffusion reaction model
 Implementation
 Case studies
 Discussion
 Conclusions
 Appendix A: Additional implementation details
 Appendix B: Distribution of terminal velocities for microplastics
 Appendix C: Additional convergence results
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References