Modeling Microplastic Transport in the Marine Environment: Testing Empirical Models of Particle Terminal Sinking Velocity for Irregularly Shaped Particles

Microplastic (mP) pollution has been indicated as an area of concern in the marine environment. However, there is no consensus on their potential to cause significant ecological harm, and a comprehensive risk assessment of mP pollution is unattainable due to gaps in our understanding of their transport, uptake, and exchange processes. This research considers drag models that have been proposed to calculate the terminal settling velocity of regularly and irregularly shaped particles to assess their applicability in a mP modeling context. The evaluation indicates three models that predict the settling velocity of mPs to a high precision and suggests that an explicit model is the most appropriate for implementation in a mP transport model. This research demonstrates that the mP settling velocity does not vary significantly over time and depth relevant to the scale of an ocean model and that the terminal settling velocity is independent of the initial particle velocity. These findings contribute toward efforts to simulate the vertical transport of mPs in the ocean, which will improve our understanding of the residence time of mPs in the water column and subsequently their availability for uptake into the marine ecosystem.


SI 1. Description of empirical models evaluated.
The existing empirical models for particle terminal sinking velocity that were evaluated during this research are presented and discussed in relation to the limitations and method applicability to modelling the fate and impacts of mPs within the marine environment. The empirical models are classed as either implicit or explicit depending on their mode of calculating the terminal settling velocity.

Implicit Models
The empirical models that give an expression for the drag coefficient CD are classed as implicit models. In this study, an iterative method is then applied to calculate the terminal settling velocity in which the drag coefficient CD is used to calculate the drag force and the particle acceleration is then calculated using the net vertical force. The process is repeated until the acceleration becomes negligible, at which point it is assumed the terminal settling velocity has been attained.
Stokes' law 1 states that the viscous force acting on a spherical body falling through a fluid depends directly on the particle radius, r, and velocity, w, and the viscosity of the fluid, µ. Using dimensional analysis and empirical results, Stokes proved that for spherical particles the drag force FD can be expressed as = 6 µ and used this knowledge to define the following expression for the drag coefficient: = , where Re is the particle Reynolds number (further information on the derivation is included in SI 2). Stokes' law is used as the reference law for spheres during this evaluation but is expected to overestimate the terminal settling velocity of non-spherical particles which will experience a higher drag force than the equivalent spherical particles. The expression is also only valid in the Stokes Regime where Re<0.1 and starts to fail when Re>1 2, 3 . Therefore, Stokes law may be used to estimate the terminal settling velocity of spherical mP particles but an alternative model may be more appropriate for irregular particles.
Bagheri and Bonadonna 4 developed an empirical model to predict the drag coefficient of nonspherical irregularly and regularly shaped solid particles falling in a fluid based on the particle Reynolds number, shape and orientation and the particle-to-fluid density ratio. The expression proposed for the drag coefficient is as follows: Where Re is the particle Reynolds number, KS is the Stokes drag correction and KN is the Newtons drag correction. The drag corrections were parameterised experimentally, and the steps required for their calculation are outlined in SI 3. As this model and all the correlations within it are derived empirically, it is not recommended that it is applied beyond the limits of the experimental data. The data was collected using both a settling column and a wind tunnel and so the model is valid in a range of fluid regimes (Re<<1 to Re<3x10 5 ). The particles used during the settling column experiments had a large size range, from 155µm to 1.8mm in the settling column experiments and 10.9mm to 61.2mm in the wind tunnel experiments, and covered a range of morphologies, including irregular volcanic particles, cylinders, ellipsoids, discs, and parallelepipeds.
The implicit model by Dioguardi et al. 5 presents an empirical expression for the drag coefficient that is dependent on the particle shape. The resulting expression for the drag coefficient is: Where Re is the particle Reynolds number and Ψ is the Dellino Shape Factor, which is the ratio of the particle sphericity to the particle circularity. The steps required for the calculation of these parameters are outlined in SI 4. The expression converges to the drag coefficient of a perfect sphere given by Haider and Levenspiel 6 when Ψ=1. As with the model by Bagheri and Bonadonna 4 , the model by Dioguardi et al. 5 should not be applied beyond the tested conditions. As the database used to derive the model includes Reynolds numbers ranging from 0.03 to 1x10 5 and particles with Dellino shape factors from 0.335 to 0.943, it is valid over a range of fluid regimes and applicable to irregular particle shapes.
The newest implicit model considered was developed by Zhang and Choi 7 and proposes that the Aschenbrenner shape factor is used to calculate the drag coefficient to more readily distinguish the behaviour of fibrous particles. Their expression for the drag coefficient was derived using linear regression and is as follows: = .
. . where Re is the particle Reynolds number and ASF is the Aschenbrenner shape factor. The process required to implement this model to calculate the terminal settling velocity, including the steps required to estimate the required particle properties, is outlined in SI 5. This model was derived using 70% of the fibre data in the dataset by Van Melkebeke et al. 8 , for which the Reynolds number varies from 6.9 to 3630. It has been shown that, as well as fibres, the model reasonably predicts the settling velocity of fragment and film particles, but it is not suited to predicting the settling velocity of weathered particles.

Explicit Models
Empirical models which give an expression to directly calculate the terminal settling velocity of particles falling in a fluid are known as explicit models. These models are computationally more efficient than implicit models since they do not require an iterative calculation.
Dietrich 9 developed an empirical equation to explicitly calculate the settling velocity of natural particles using the results of settling velocity experiments. The model expresses separately and quantitatively the effects of the particle density, size, shape, and roundness on the particle settling velocity but has no dependence on the particle Reynolds number Re. The final equation expresses the settling velocity in its dimensionless form and is stated to be reasonably accurate in calculating the settling velocity of natural particles: * = 10 Where W* is the dimensionless settling velocity, R1 is a fitted equation for particle size, R2 is a fitted equation for particle shape and R3 is a fitted equation for particle roundness. The definition of these equations for the implementation of the model is outlined in the SI 6. The model uses the Corey Shape Factor (CSF) to quantify particle shape. However, due to the definition of the fitted equation for particle shape R2, the model fails when CSF<0.15 and therefore is not recommended for use for particles with CSF<0.2.
A method to explicitly predict the settling velocity of plastic and natural particles of different shapes, ranging from one dimensional (fibres) to three dimensional (pellets, spheres), across a range of flow regimes was developed empirically by Francalanci et al. 10 and validated using an independent dataset. The terminal settling velocity is expressed in its dimensionless form as: * = * + 0.75 * Where W* is the dimensionless settling velocity, D* is the dimensionless particle size, C1 and C2 are coefficients, and n is an exponent which was obtained during calibration of the equation. The full procedure required to implement this model is outlined in SI 7. The experiment included particles with 1D, 2D and 3D morphologies with sizes of 1.68 to 5.44mm and so the model is applicable to a wide range of particle shapes and sizes. However, the model is limited in that it only considers quiescent fluids and neglects the effect of a bulk of particles on the settling velocity. The model is independent of the particle Reynolds number.
An alternative drag law was developed by Yu et al. 11 to explicitly calculate the settling velocity using the Corey Shape Factor (CSF) and the particle sphericity Φ. This method provides an expression for the drag coefficient Cd which can be used to calculate the terminal settling velocity with the following equation: Where ws is the particle settling velocity, is the fluid viscosity, g is the gravitational acceleration, is the particle density, is the fluid density, d* is the dimensionless particle diameter and Cd is the drag coefficient proposed by Yu et al. 11 . The process to implement this model, including the calculation of the required parameters, is outlined in SI 8. As the model was derived using data that lies predominantly in the transitional regime (1<Re<1000), its applicability to the Laminar regime (Re<1) and turbulent regime (Re>1000) is uncertain. Since the dataset used to derive the model contained 699 mPs ranging in size from 0.0025mm to 0.0356mm and included pellets, cylinders, spheres, fragments, fibres, fishing lines and angular and nodular particles, the model is applicable to a wide range of particle properties.

Derivation:
Stokes' law states that the viscous force acting on a spherical body falling through a fluid depends directly on the radius r (m), and velocity w (m/s) of the sphere and the viscosity µ (kgm -1 s -1 ) of the liquid: ∝ Using dimensional analysis, it can be shown that:

= µ
Where k is a coefficient of proportionality, which has been shown empirically to equal 6π. Therefore, the viscous drag force is: Equating this to the expression for drag force containing the non-dimensional drag coefficient CD gives: Where S is the effective area of the particle (m 2 ), ρw is the fluid density (kg/m 3 ) and w is the particle settling velocity (m/s). For a spherical particle, defining the effective area as the projected area 3 gives = : Which results in the expression:

= 24
Where Re is the dimensionless Reynold's number = . Implementation:

Description
The expression for the drag coefficient of non-spherical irregularly and regularly shaped solid particles falling in a fluid (gas or liquid) is:

+ 5330
Where Re is the particle Reynolds number = .
The shape descriptors particle flatness, f, and elongation, e, are used in the model to quantify particle shape since they are straightforward to measure and calculate: = = Where S, I and L are the shortest, intermediate, and longest length of the particle. These descriptors are then used to define the two new shape descriptors: Stokes form factor FS and Newtons form factor FN: The experimental results demonstrated that the shape dependent Stokes drag correction ks correlates with the Stokes form factor FS as below: Non-linear regression was used to derive the correlation of the shape dependent Newtons drag correction kN with Newtons form factor FN and the particle to fluid density ratio ρ' = :

Description
The expression proposed for the drag coefficient is: Where Re is the particle Reynolds number = and Ψ is the Dellino 12 shape factor which is expressed as: where Φ is the particle sphericity (dimensionless), Χ is the particle circularity (dimensionless), Asph is the surface area of the volume equivalent sphere (m 2 ), Ap is the particle surface area (m 2 ), Pmp is the perimeter of the maximum projection area (m) and Pc is the perimeter of the circle equivalent to the maximum projection area Amp of the particle (m). = − Δ i. If acceleration>0.001m/s 2 , step forward one timestep and restart the loop at point (5a). ii. If acceleration<0.001m/s 2 , assume that terminal settling velocity has been obtained and end the calculation iterations.
SI 5: Procedure to implement the explicit model by Zhang and Choi (2021) 7

Description
The proposed drag model is as follows: . Where Re is the particle Reynolds number = and ASF is the Aschenbrenner shape factor = , a is the longest length of the particle (m), b is the intermediate length of the particle (m) and c is the shortest length of the particle (m).  Derivation: Dietrich 9 expresses the dimensionless settling velocity W* as: * = 10 Where R1 is a fitted equation for particle size: R2 is a fitted equation for particle shape: R3 is a fitted equation for particle roundness: . .
The nondimensional parameters used in these expressions are: Corey Shape Factor CSF:

= √
Where a, b and c are the longest, intermediate and shortest axis of the particle. Note that a, b and c are perpendicular. The CSF varies from 0.0 to 1.0, with a smaller CSF indicating a flatter particle.
Powers roundness index, P, which is determined by observation and varies from 0.0 (perfectly angular) to 6.0 (perfectly round).
Dimensionless particle size D*: * = ( − ) Where ρs is the particle density (kg/m 3 ), ρ is the fluid density (kg/m 3 ), Dn is the nominal spherical diameter i.e., the diameter of a sphere with the same volume as the particle (m) and ν is the kinematic viscosity of the fluid (m 2 /s).
Using the dimensionless settling velocity W*, the particle settling velocity ws (m/s) can then be calculated as follows: where Dg is a reference diameter that includes the influence of particle shape:

= √
Where CSF is the Corey shape factor, a, b, and c are the longest, intermediate, and shortest dimension of the particle respectively (m), g is the gravitational acceleration (m/s 2 ), is the kinematic viscosity (m 2 /s), ρp is the particle density (kg/m 3 ) and ρf is the fluid density (kg/m 3 ).
The dimensionless settling velocity W* can then be used to calculate the actual settling velocity using the following equation: The settling velocity of irregularly shaped particles is then defined as: Where is the fluid kinematic viscosity (m 2 /s), ρs is the particle density (kg/m 3 ) and ρf is the fluid density (kg/m 3 ).
The dimensionless particle diameter is defined as: * = − where dn is the volume equivalent spherical diameter                 Figure S23: The impact of the choice of initial velocity on the modelled settling velocity of six particles when using Zhang and Choi's model 7 and taking the particle projection area as the effective area in the calculation of the drag force that were randomly extracted from the dataset by Van Melkebeke et al. 8 . Figure S24: The impact of the choice of initial velocity on the modelled settling velocity of six particles when using Zhang and Choi's model 7 and taking the particle surface area as the effective area in the calculation of the drag force that were randomly extracted from the dataset by Van Melkebeke et al. 8 .

SI 12
: Results for each model tested during the evaluation of the variation in the terminal settling velocity over the range of density in the ocean and the impact on the models of assuming a constant terminal settling velocity. Figure S25: The influence of fluid density on the terminal settling velocity of six random particles using Stokes' model 1 .                 Figure S42: Comparison of the distance travelled in attaining the terminal settling velocity to the distance travelled if the particle sank constantly at the terminal settling velocity in the equivalent period of time when using Zhang and Choi's model 7 and taking the particle projected area as the effective area in the calculation of the drag force.The solid line indicates the ideal fit where there is no difference in the distance travelled and the dotted lines indicate the distance travelled at a constant velocity is ±30% of the distance travelled whilst attaining terminal settling velocity. The dashed line indicates the best fit line in the form y=mx that was obtained using linear regression. Figure S43: Comparison of the distance travelled in attaining the terminal settling velocity to the distance travelled if the particle sank constantly at the terminal settling velocity in the equivalent period of time when using Zhang and Choi's model 7 and taking the particle surface area as the effective area in the calculation of the drag force.The solid line indicates the ideal fit where there is no difference in the distance travelled and the dotted lines indicate the distance travelled at a constant velocity is ±30% of the distance travelled whilst attaining terminal settling velocity. The dashed line indicates the best fit line in the form y=mx that was obtained using linear regression.