Lagrangian Modeling of Marine Microplastics Fate and Transport: The State of the Science

Microplastics pollution has led to irreversible environmental consequences and has triggered global concerns. It has been shown that water resources and marine food consumers are adversely affected by microplastics due to their physico-chemical characteristics. This study attempts to comprehensively review the structure of four well-known Lagrangian particle-tracking models, i.e., Delft3D—Water Quality Particle tracking module (D-WAQ PART), Ichthyoplankton (Ichthyop), Track Marine Plastic Debris (TrackMPD), and Canadian Microplastic Simulation (CaMPSim-3D) in simulating the fate and transport of microplastics. Accordingly, the structure of each model is investigated with respect to addressing the involved physical transport processes (including advection, diffusion, windage, beaching, and washing-off) and transformation processes (particularly biofouling and degradation) that play key roles in microplastics’ behavior in the marine environment. In addition, the effects of the physical properties (mainly size, diameter, and shape) of microplastics on their fate and trajectories are reviewed. The models’ capabilities and shortcomings in the simulation of microplastics are also discussed. The present review sheds light on some aspects of microplastics’ behavior in water that were not properly addressed in particle-tracking models, such as homoand hetero-aggregation, agglomeration, photodegradation, and chemical and biological degradation as well as additional advection due to wave-induced drift. This study can be regarded as a reliable steppingstone for the future modification of the reviewed models.


Introduction
Aquatic environments, including oceans, coastal regions, estuaries, seas, and rivers, as major planetary water resources, have been affected by marine plastic debris accumulation and distribution, which has triggered global environmental concern and irreversible pollution consequences [1][2][3]. Macro (> 1 cm) [4], meso (5 mm-2 cm) [4], micro (< 5 mm) [5], and nano-size (1-100 nm) [6] plastics pollution has become a critical environmental problem due to their persistence, toxicological properties, and destructive effects on not only Earth's hydrosphere [7] but also aquatic organisms [8][9][10][11], wildlife, and human health [12][13][14][15]. The tendency of microplastics to adsorb and react with wastewater agents is relatively high due to their physical properties, such as hydrophobicity [16], as well as the large ratio of their surface area to volume, allowing them to adsorb toxic pollutants (e.g., heavy metals) and invasive species [17]. Microplastics are classified into two categories: primary and secondary. Primary microplastics are mainly detected with an original size of less than 5 mm in textiles and some personal care products [2,18]. They can be transported through discharges of water treatment plants, rivers, surface run-off, and wind into aquatic environments [19]. Large plastic debris under fragmentation may produce secondary microplastics through some processes including physical, chemical, and biological interactions and photo-degradation [20,21]. The secondary microplastics originate mainly from industrial resin pellets, fishing nets, and other thrown-away plastic litters [22], and they can also be formed from the fragmentation of any larger plastic particle in the ocean, regardless of its source. Microplastics in the environment majorly consist of secondary microplastics [23]. The possibility of microplastics breaking down into nanoplastics may increase during their exposure to the environment, leading to higher environmental risks given the nano-size nature of these materials [24]. In addition, plastic pollution impacts have been observed at economic and social scales [25]. The tourism industry accounts for a considerable percentage of the gross domestic product (GDP) of many countries, which indicates the market value of all services provided within a country's territory during a specific temporal period; thus, governments schedule various plans to attract more tourists and introduce different sightseeing places in order to boost tourism income. Nevertheless, beaches confront some challenges with respect to their sustainability in terms of preservation and operation from the touristic and economic perspectives. The uncontrolled release of microplastics causes the pollution of the seawater and the shore.
Plastic debris enters the ocean directly or through rivers by uncontrolled dumping from land-based or maritime sources, such as fishing or shipping sectors, and recreational activities [26][27][28]. Plastics will be degraded and broken down into microplastics (i.e., particles less than 5 mm in size) after entering the ocean if not treated at early stages. Microplastics are found throughout the water column in oceans and also on the seabed [29,30]. Numerous studies have been conducted to quantify and characterize microplastics in aquatic environments [4,[21][22][23] including lakes [23,31], rivers [32][33][34][35][36], and shorelines [37] by employing different field surveys and using a set of lab experiments. Recent advances in computer science and technology have also promoted the application of Machine Learning (ML) approaches in the automated quantification and characterization of microplastics [38][39][40][41][42][43].
Microplastics in water undergo some transformational processes such as biofouling [44], degradation and fragmentation [16,18,45], and hetero-and homo-aggregation [45,46], known also as particle-particle interactions [47], which influence their fate in their receiving environment. Accordingly, their physical (mainly size, density, and shape), chemical, and biological properties may change as a result of their ambient environment properties, including sea surface temperature, sunlight, salinity, etc. The growth of ultrafine bacteria or microorganisms and the attachment of organic material, invertebrates, and algae on the exterior surfaces of these particles cause biofouling that results in not only density and size increases in the particles but consequently decreases buoyancy [48]. Thus, it can be concluded that the sinking process can be accelerated by biofouling [49] and may result in the sedimentation of fouled particles. In addition, microplastics have been widely pervasive in marine habitats [26,50] due to physical transport processes including horizontal and vertical advection and diffusion. Previous studies have shown that microplastics have reached very remote regions in the Arctic [51] and Antarctic [52,53]. They can be found on the sea surface [54], the seabed [55], and also in deep-sea sediments [56,57]. Predicting such transport processes, given their dependencies on different hydrological and atmospheric variables spanning over a wide temporal and spatial scale, is important for the accurate simulation of microplastics' fate and transport in water bodies. The transport of microplastics in water bodies is controlled by advection and dispersion driven by ambient currents, Stokes drift, sinking, resuspension, beaching and washing-off, deposition, aggregation, degradation, and fouling, which need to be simulated by numerical models.
Four sources for transport processes can be considered: (1) ambient currents (e.g., tidal and wind-driven flows), leading to transport by advection; (2) wind-induced drift, i.e., windage, which is an additional advection term because of wind drag and has an undeniable effect on plastic particles' displacement [58]-in other words, the wind-driven drag force exerted on the water's surface causes the drifting of floating particles exposed to wind; (3) wave action; and (4) buoyancy, which influences the vertical transport. Indeed, the influence of buoyancy on vertical transport and associated vertical mixing can be significant [59].
Physical processes involved in the transport and fate of microplastics can occur at different spatio-temporal scales, which makes the simulation challenging. In this respect, short-term processes (e.g., particle-particle interactions at local and regional scales) may last for several hours or a few days [60]; however, medium-term processes including biofouling, sedimentation, and advection from freshwater bodies to marine environments may occur within several days or weeks [44]. Long-term processes need months or even years to be completed [61]. For instance, some processes, such as degradation (or fragmentation), sinking, and transport (due to advection and diffusion) from coastal water bodies to oceans are long-term processes that increase computational costs. From the spatial scale perspective, the simulation of physical processes in large study areas with global scales is of orders of thousands of square kilometers, while the computational domains for regional and local areas are of orders of hundreds and tens of square kilometers, respectively.
Numerical simulation has been widely employed in the past decades to predict marine microplastics' behavior [62][63][64][65][66][67][68]. There are also some studies that have predicted the main sources, transport, and fate of microplastics under different environmental conditions based on numerical simulations [69][70][71][72][73]. Some of these studies have employed particle-tracking [61] and hydrodynamic models [74], individually or coupled with each other [75][76][77]. Numerical modeling is inexpensive by far in comparison with experimental measurements and in situ observations and can predict potential accumulation zones in remote areas where water quality monitoring is not feasible. Furthermore, numerical models can simulate some key physical processes using fine grid resolution, given the everincreasing computational capacity and the ongoing development of numerical techniques. Moreover, these models are capable of simulating different vertical and horizontal flow patterns encompassing microplastics with good accuracy.
Numerical simulations of microplastics' fate and transport in water are generally carried out using Eulerian and Lagrangian models. Unlike Lagrangian trajectory models that use a movable coordinate system in space, Eulerian grid models typically employ a non-movable frame of reference [58]. In other words, a Lagrangian framework, known as the particle-tracking method, tracks only individual particles, while an Eulerian approach addresses advection and diffusion for given locations (computational nodes) [78]. Numerical Eulerian hydrodynamic models can be coupled with Lagrangian models. Accordingly, hydrodynamic modeling provides the ambient hydrological parameters as inputs for particle-tracking models (PTM), which are known as promising, cost-effective tools for interpreting the fate and transport of microplastics and providing better insights into microplastics' behavior in dynamic flow systems [79]. Although some studies have used an Eulerian approach for the mixing and dispersion processes in dynamic water systems [80], some other studies have used Lagrangian models for tracking particles to investigate plastic particles' fate and transport [58]. Lagrangian models solve for a lower number of equations per node so that they do not need to solve the equations for each location in the mesh. However, given the pros and cons of both Eulerian and Lagrangian approaches, they can be employed to complement each other. The particle-tracking model coupled with the hydrodynamic model can provide invaluable details about the potential trajectories of plastic litter from local [81][82][83] to global scales [70] that lead to the assessment of microplastics' transportation, distribution, residence time, and fate under various conditions [58,67,70,[84][85][86][87].
Particle-tracking methods simulate the trajectory of particles and their transport in water bodies. They are applicable to a wide range of plastic particles' transport in their surrounding environment [78,88]. PTM has been employed in many studies as a reliable alternative for field surveys. This is due to the logistic challenges involved in field surveys [70,72,[89][90][91][92]. Most Lagrangian particle-tracking models do not simulate the flow hydrodynamics. They can be considered as an important tool to simulate physical processes and transport mechanisms that microplastics undergo.
Unlike the structured mesh, the unstructured mesh enables numerical models not only to mesh complex geometries more easily but also to address arbitrary positions, even though it should be noted that greater memory is needed for executing the simulation for unstructured mesh. Hydrodynamic values in each element of the Eulerian mesh are interpolated in the same location of the particle (as the host element) in PTM.
On the other hand, the spatio-temporal interpolation of the flow characteristics (mainly including the velocity, diffusivities, salinity, sea surface height, temperature, etc.) in both horizontal and vertical directions may be necessary, especially when the duration between successive hydrodynamic outputs is different in comparison with the time step of particles' motion. Thus, the hydrodynamic information provided on an Eulerian grid, as an input to the PTM model, should be interpolated for the spatio-temporal location of the particles before being read by the PTM [61]. For instance, the particle's velocity at each location is found using the spatial and temporal interpolation of the hydrodynamic model's velocity.
Most available numerical models have limitations in considering all effective factors relevant to microplastics' fate and transport. For example, some models have ignored beaching, sinking, and fragmentation processes as well as biofouling, deposition, resuspension, washing-off from the beach, or wind drift [67,85,93].
Among the numerical models that are employed for the simulation of marine microplastics' behavior, we studied four PTMs, i.e., D-WAQ PART, Ichthyop, TrackMPD, and CaMPSim-3D. D-WAQ PART is the particle-tracking module of Delft-3D, which is an open-source simulating software equipped with different modules for predicting hydrodynamics in an aquatic environment. Ichthyop, as a free Java-based numerical model, has been employed for tracking plastic particles, although it has been originally used for studying the dynamics of ichthyoplankton. TrackMPD has been developed in recent years to simulate microplastic debris fates and transport. The latest developed numerical model that can be employed as a particle-tracking model is Canadian microplastic simulation (CaMPSim). It utilizes a coupling of an Eulerian and Lagrangian numerical model in the simulation of marine microplastics' behavior.
The current paper systematically reviews recent progress in the numerical modeling of the transport and distribution of microplastics in marine environments with the aims of (1) reviewing the structure of Lagrangian particle-tracking models, (2) determining the substantial improvements that are still needed, and (3) discussing the consideration of physical processes that microplastics undergo and their physical properties in the simulation of microplastics' fate and transport. In this respect, different 2D and 3D numerical models are categorized, as shown in Table 1.

Materials and Methods
As stated above, the concept and framework of different numerical models for modeling microplastics' fate and transport are reviewed in this study in detail, considering various ranges of microplastics' behavior characteristics and key physical processes affecting their variations. Lagrangian particle-tracking models are investigated in terms of considering physical properties (and hydrodynamic aspects), physical transformational processes, physical transport processes, and numerical schemes ( Figure 1). In other words, the reviewed models have been compared based on capability in considering the physical properties of microplastics (mainly size, shape, and density), and in simulating the physical transformational and transport processes that microplastics undergo. Moreover, the numerical methods that have been employed by each of them in terms of considering the type of mesh they can read, their computational schemes, and the methods they employ for implementing interpolation have been addressed. Finally, limitations and gaps in the reviewed models and needed areas of improvement are determined and discussed in full.   . The structure of Lagrangian particle-tracking models in the simulation of microplastics' physical properties, physical transformation processes, physical transport processes, and numerical methods.

Physical Properties
The behavior and transport of microplastics mainly depend on major physical properties including the size, density, and shape [44,66,96,99] of the particles. Microplastics are identified in varied sizes and shapes, and denser particles normally settle faster even with similar shapes and sizes [44]. Moreover, larger particles in terms of size but with the same density experience an increase in the ratio of the gravitational force applied on the particle to the fluid viscous resistance, which results in higher settling velocities and faster settling [61]. In addition, the shape of microplastics can affect their movements and settling velocity as well [100]. It has been observed that elongated particles such as films and fibers, known as larger particles, settle faster than smaller, circular particles [96]. Microplastics come in various regular (e.g., spherules, cylinders, and beads) and irregular shapes with different geometries (pellets and fibers). Regular shapes usually belong to primary microplastics while irregular shapes often can be seen among the secondary particles as a result of degradation [99]. Moreover, the physical properties of microplastics undergo changes due to the impacts of different processes and ambient factors such as algae invasion [101], salinity [102], UV index [103], temperature [103], and Stokes drift generated by wave action [104]. For instance, unlike freshly emitted particles, the edges of microplastics become smooth during aging as a result of mechanical degradation since they are under polishing imposed by other fragments [105]. Different equations used in the reviewed models for modeling hydrodynamic aspects (including microplastics' physical properties) pertinent to the transport of microplastics have been summarized in this study in Table 2.
is the settling velocity of each particle at time t (hours), indicates the local concentration of particles kgm −3 ) is the exponent for adjusting concentration-dependent settling velocities, A 0 [106] is the non-cyclic component of the settling velocity ms −1 ), A 1 represents the amplitude of a periodic sinusoidal variation in time ms −1 ), is the period of the sinusoidal variation (hours) and ∅ is the phase lag for the sinusoidal variation (hours ).
Settling velocity of spherical particles ϑ is the kinematic viscosity of water m 2 s −1 ), d * represents the dimensionless diameter of the particle, and R is the radius of the particle (mm) .
TrackMPD [99] 3 Dimensionless diameter of the particle g is the gravitational acceleration ms −2 and, ρ p and ρ w indicate the density of the particle and the suspending medium, respectively, with the same units gm −3 ) . [99] Settling velocity of cylindrical particles L indicates the length of cylinders (mm) . [99] and ρ w gm −3 ) indicate the density of the particle and the surrounding medium, respectively, g is the gravitational acceleration ms −2 and, Drag coefficient Re p is particle Reynolds number and Ψ is the particle shape factor [108]. Equation (6) has been reported as the best fit for modeling drag when interpreting microplastics' sinking [109]. [98] * The Ichthyop model uses the Stokes law in calculating the settling velocity.

Biofouling and Degradation (Physical Transformation Processes)
According to previous studies, microplastic particles are subjected to some natural transformation processes including fouling, degradation, and mechanical breakdown during their residence time in an aquatic ecosystem that may lead to changes in their physical properties [58]. Biofouling is controlled by different ambient environmental factors such as temperature, salinity, and the availability of sunlight, which is necessary for biofilm (including fine and ultrafine organisms and plankton) growth [76]. On the other hand, particles become more brittle and break apart into smaller pieces as a result of mechanical degradations by weathering, wave actions, interaction with other particles, the seabed, and the shore [110]. Degradation is typically a slow process that takes over 50 years to be fully completed for plastics since they are highly resistant to weathering [7,110]. However, in the swash zone [111] and some energetic coastal areas (e.g., salt marshes), natural abrasives including sediment and rock are present, and the degradation may be completed sooner, even within 8 weeks based on [112], thereby affecting seasonal microplastics' transport. The main equations employed in the reviewed models in this study for considering biofouling and degradation processes are provided in Table 3. It should be noted that none of the reviewed models are capable of simulating different types of aggregation.
Biofouling of spherical particles ρ p gm −3 ) is the density of a fouled spherical particle, ρ 0 gm −3 ) is the density of a plastic particle, R (mm) indicates microplastic radius, and BT and ρ D represent the thickness ( µm) and density gm −3 ) of a biofilm layer, respectively.
TrackMPD-CaMPSim-3D [96] 8 Biofouling of cylindrical particles is the density of a fouled cylindrical particle, ρ 0 gm −3 ) is the density of a plastic particle, R (mm) indicates microplastic radius, and BT (µm) and ρ D gm −3 ) represent the thickness and density of a biofilm layer, respectively.
TrackMPD-CaMPSim-3D [61] 10 Size(D or L) = Size 0 (1 − DR·T/100) Degradation D (mm) and L (mm) are the particle's diameter and length, respectively. Size 0 indicates the particle's initial diameter and length. DR as a constant rate shows the percentage of decrease in size of the particle per day [112] and T represents the duration of degradation from the beginning to the present time step.

Physical Transport Processes
In addition to the transformation processes (biofouling, degradation, aggregation), numerical fate and transport models can simulate different physical processes driving microplastics' transport in water such as advection, diffusion, windage, resuspension, beaching, and washing-off. The transport of microplastics is generally governed by a generic advection-diffusion reaction, as follows [98]: where C indicates the concentration at time t, U is the velocity of a surrounding medium, K represents the dispersion coefficient, and ρ(C) shows the reaction function. Given the aforesaid concept, the main equations used in the reviewed models for addressing the main physical processes controlling microplastics' transport in aquatic systems are provided in Table 4. Table 4. The mathematical or empirical equations employed in the reviewed models for considering the physical transport processes that microplastics undergo in aquatic environments.

No.
Mathematical or Empirical Relationships Application Defining Parameters Employed in Ref.
Trajectory Calculation X is the 3D position vector with components in the horizontal plane, d is the total derivative, u is the 3D velocity vector, and X is the fluctuation of the position vector.
D-WAQ PART [58] 13 V adv = V f low + V windage Advection velocity V adv ms −1 represents the velocity field, which is responsible for advection. V f low ms −1 is the velocity of flow, and is an empirical wind drag coefficient that is related to particle characteristics and V w ms −1 represents the wind velocity at 10 above sea level. [58] 14 are calibrating coefficients. It should be noted that the upper and lower limits of K h(x,y) are equal to at and a respectively, and K h(x,y) increases with time. It has been used in Equation (11).

Turbulent vertical diffusion coefficient
D z is the vertical dispersion coefficient, C µ indicates a calibrating constant for local equilibrium shear layers that is equal to 0.09, approximately. L is the mixing length (m), represents turbulent kinetic energy, and σ C is the Prandtl-Schmidt number. It has been used in Equation (11).
D-WAQ PART [107] 16 P = 0.5 −t/T The probability of being washed-off P is the probability of a particle being washed off, t is the time step, and T is the half-life of plastic litter.
D-WAQ PART-TrackMPD [87] 17 v part = v water + 1 24 gd 2 ∆ρ ρw ϑ −1 ln 2I d + 1 2 Vertical velocity of particle v water is vertical velocity of water, g is the gravitational force ms −2 , d and I are prolate spheroid axes (indicating the area of a particle), ∆ρ = ρ part − ρ w , ρ part gm −3 ) is particle density and ρ w gm −3 ) is the density of water, and ϑ indicates the kinematic viscosity.
Ichthyop [114] 18 Horizontal turbulent diffusion coefficient ε is the turbulent dissipation rate and l is the cell size. It has been used in Equation (11).

19
x n+1 = x n + U∆t i y n+1 = y n + V∆t i z n+1 = z n + W∆t i First-order Euler method for advection ∆t i is an internal time interval and U = (U, V, W) is the velocity vector. This is a first-order discretization of Equation (11).
TrackMPD [61] 20 Calculation of horizontal displacement due to diffusion Equation (11) R indicates a random number with an average of zero and r is its standard deviation that equals 1, and K h is horizontal diffusivity m 2 s −1 ). [61] 21 Calculation of vertical diffusion term Equation (11) K v is vertical diffusivity m 2 s −1 ). [61] 22 Velocity field as a function of current and wave velocities U f is the current velocity, U w is the wind velocity, ρ air and ρ w indicate the density of air and water, S above and S below represent the cross-sectional areas of spherical and cylindrical particles in dry and wet conditions, respectively. Note that Equation (17) and (21) use the same principles and their difference is only in some coefficients. TrackMPD The ratio of the dry cross-sectional area of the particles to its wet cross-sectional area.
α is defined by Equation (24). [61] 24 α = 2arcccos 1 − h R Parameter for Equation (23) h is the Archimedean force, R is the radius of the particle. [61] Archimedean force ρ gm −3 ) and ρ w gm −3 ) are the density of the particle and water, respectively. [61] 26 Resuspension flux R jmax is the maximum resuspension constant for microplastics. τ and τ crit are shear stress and critical shear stress, respectively.
CaMPSim-3D [116] 27 Actual shear stress ρ w gm −3 ) is the density of seawater, g ms −2 is gravitational acceleration, v w is the mean velocity of the seawater, and Chezy is the Chezy coefficient m 0.5 s −1 ). [116] where v p is the vertical displacement velocity of particles, v i is particle velocity at the computational cells, ∝ i is vertex, and v s is particles settling/rising velocity (obtained from Equation (5)). CaMPSim-3D [98] 29

Calculation of diffusion term
where K is the diffusion coefficient (obtained from Equation (14)), ∆t is the time step, R is a random number, and K p is the diffusion coefficient. [98]

Numerical Lagrangian Models
In PTM models, a velocity value is allocated to each particle. This is typically done by the interpolation of values from an Eulerian hydrodynamic model onto the position of microplastic particles. The advection velocity typically includes the settling velocity, the velocity of flow provided by the hydrodynamic model, and a diffusion-based velocity that is randomly generated, typically based on a Gaussian distribution [76,115]. Many PTM models use numerical methods such as the fourth-order Runge-Kutta method to solve equations in time.
In most PTM models, the number of required particles will be determined based on several factors such as the simulation period, the extent of particle dispersion, the cells' size in the computational domain, computational resources, etc. [107].
The majority of PTM models calculate the settling velocity considering the impacts of the aforementioned physical properties such as density, size, and shape. The settling velocity of each particle, v s , is defined based on the relationships, as shown in Table 2.

D-WAQ PART
Delft3D, as open-source software, can be applied for simulating flows, the transport of sediments, tidal currents, water quality, waves, etc. in various water settings. Delft3D has several modules as shown in Figure 2. The modules are linked to each other in order to simulate the multi-dimensional transport and hydrodynamics (in 2D or 3D) [117]. Particle transport and dynamical spatial distribution are simulated by the PART module of Delft3D using particle-tracking modeling using the inputs provided by the FLOW module (i.e., two-or three-dimensional flow data). In other words, velocity and other hydrodynamic parameters are computed by the FLOW module using an offline coupling and provided as inputs for the D-WAQ PART module.   D-WAQ PART is not capable of considering the impact of the biofouling process on the trajectories of micro-size particles via increasing the density and size of the particles since it cannot consider the thickness of the biofilm layer, and so, the density of a fouled particle. D-WAQ PART also lacks a specific tool to provide insight and appropriately predict complex processes including degradation as a long-term and aggregation as a short-term process.
Particle trajectories are computed over time in three dimensions by modeling the advection and dispersion processes [118]. As stated above, each individual particle can be advected by water currents and wind or wave drag as additional sources and also displaced randomly due to diffusion in order to explain unresolved processes.
The velocity of particles due to advection includes the local current velocity prepared by the hydrodynamic model in the Eulerian velocity field and the wind velocity at the location of a particle as shown in Equation (13) [58]. The wind-drag coefficient is related to object characteristics as an empirical parameter. In addition, D-WAQ PART employs a random-walk particle-tracking model in order to describe the movement of particles due to diffusion processes. The diffusion is described using kinetic gas theory and the theory of Brownian motion [107]. D-WAQ PART addresses beaching and washing-off using the probability of being washed off, as shown in Equation (16).
Many numerical models, including D-WAQ PART, are not capable of reading the unstructured mesh information provided by an Eulerian model. However, some of the other Delft-3D modules (e.g., D-Water Quality and GPP) can solve the advection-diffusion equation numerically on arbitrary, unstructured grids and structured, shaped grids such as triangles, rectangles, and curved elements.
On the other hand, D-WAQ PART utilizes the linear and block interpolation implemented for all time functions between time intervals [107]. In the linear interpolation, D-WAQ PART interpolates the hydrodynamic inputs linearly between the time records while it interpolates the two different subsequent inputs with similar time as the block interpolation.

Ichthyop
Ichthyop is a free Java-based particle tracking numerical model which was originally developed to study ichthyoplankton dynamics [97]. Ichthyop has also been employed for modeling the fate and transport of microplastics and toxic algae [76,119]. The model can use three-dimensional hydrodynamic information including temperature, velocity, and salinity, simulated by oceanic numerical models such as the Regional Oceanic Modelling System, called ROMS, and the Model for Applications at Regional Scale, named MARS, as input [97]. Ichthyop is capable of providing two types of microplastics simulations: (1) modeling plastic particle trajectories, and (2) ensemble simulation (i.e., executing a series of numerous particle tracking simulations) [97]. Ichthyop can also simulate the trajectory of more than one million particles as passive tracers [119].
Similar to D-WAQ PART, Ichthyop does not consider the integration of complex particle transformation processes including biofouling, different types of aggregation, degradation, and fragmentation.
Ichthyop uses the Runge-Kutta method [97] for solving the advection equation [76]. Horizontal and vertical advection and dispersion can be simulated in Ichthyop [76,97,113,119,120]. The model follows D-WAQ PART and uses Equation (11) to simulate the horizontal advection. Ichthyop utilizes a terminal velocity (v part ) in order to address the vertical movement of particles, as defined in Equation (17) [114].
Ichthyop uses Equation (18) for the consideration of horizontal diffusion. Accordingly, horizontal diffusion is defined as a function of the turbulent dissipation rate and the cell size of the computational domain [115]. Furthermore, the vertical diffusion of particles is predicted based on a proposed random-walk simulation [121].
Ichthyop is incapable of simulating some other effective physical processes on microplastics' fate and transport, such as beaching and washing-off, windage, resuspension, and aggregation, which can be regarded as an important limitation.
Hydrodynamic parameters are simulated on a rectangular structured mesh in some hydrodynamic models such as ROMS and MARS. Then, they are interpolated linearly and extracted at time steps in order to be prepared as inputs for Ichthyop.

TrackMPD
Amongst the three-dimensional numerical models, TrackMPD is an open-source model that has been released in recent years [61] based on the framework of the Particle Tracking and Analysis Toolbox (PaTATO) [122]. TrackMPD is compatible with various ocean models in order to address the fate and transport of marine microplastic debris in coastal areas. TrackMPD is capable of simulating microplastics' behavior in three dimensions considering different physical processes, including beaching and washing-off, biofouling and degradation, advection, windage, turbulent dispersion, and sinking, as well as different physical properties such as density (Equations (7) and (8)), size (Equation (10)), and shape (Equations (2)-(4), (7) and (8)). The physical properties of microplastics can be parametrized in TrackMPD, particularly in the calculation of the settling velocity. TrackMPD utilizes the empirical equations of [99] and [123] in estimating the settling velocity of microplastics with different shapes as a function of particle density and radius. Accordingly, the settling velocity (v s ) of spherical particles (ms −1 ) is obtained based on a calibrated formulation provided by [123], as shown in Equations (2) and (3), while this is different, as determined in Equation (4), for cylindrical fragments (mms −1 ) [99].
The effects of biofouling on the fate and transport of microplastics can be investigated by TrackMPD. In other words, TrackMPD uses some empirical equations in order to approximate the thickness of the biofilm layer. In this respect, the thickness of biofilm layers of spherical and cylindrical particles can be estimated based on Equations (7) and (8), respectively [96]. Accordingly, biofilm thickness can be given either a constant or variable value over time. The case of constant biofilm thickness is referred to as stationary biofouling, while the case of variable biofilm thickness is considered as non-stationary biofouling.
TrackMPD simply predicts mechanical degradation by using an empirical relation (Equation (10)). In this respect, the size of plastic particles (their diameter or length) decreases based on a constant rate and time. While TrackMPD can simulate biofouling and degradation, it suffers from incapability in simulating some other transformational processes such as hetero-and homo-aggregation.
To simulate the horizontal and vertical movements of plastic particles due to advection, TrackMPD employs the Runge-Kutta method [61]. In this respect, the velocity of the flow is calculated iteratively in three directions at the location of a particle to consider velocity values (ms −1 ) at past and future times. Then, to calculate particle displacements (m) and obtain their new location in each direction, the previous location of microplastics is added to the multiplication of the internal time interval (∆t i ) and the flow velocity's components (U = (u, v, w)), as shown in Equation (19). Moreover, the horizontal (in x or y directions) and vertical diffusion of particles is predicted in TrackMPD based on a random-walk model, as demonstrated in Equation (20) and (21), respectively [61].
As with D-WAQ PART, TrackMPD can consider the windage effect in the calculation of particle velocity. Accordingly, the horizontal wind-induced drift of microplastics is considered using Equations (22)-(25) [61].
TrackMPD estimates the probability of microplastics' washing-off based on a Monte Carlo method, as shown in Equation (16) [61]. The probability (P) of being washed off in TrackMPD is an exponential function of the time step (t), and the half-life of plastic litter (T) is equal to the time that microplastics remain on the beach [87]. If P is higher than a random number, which varies from 0 to 1, beached microplastics can be transported (i.e., washed off) seaward at extreme tides. D-WAQ PART uses a basically similar approach in the calculation of the probability of wash-off.
Currently, TrackMPD can read hydrodynamic parameters off-line from some models such as ROMS in only a rectangular structured mesh. Moreover, hydrological properties are interpolated spatially at each time step, and then the hydrodynamic model assigns the values to particles located in the corresponding grid cells of TrackMPD [61].

CaMPSim-3D
Canadian microplastic simulation (CaMPSim), the newest three-dimensional Eulerian-Lagrangian framework [98], has been developed recently for simulating microplastics' fate and transport in different aquatic environments. CaMPSim-3D is based on the coupling of a three-dimensional Lagrangian particle-tracking model and TELEMAC modeling system as the Eulerian hydrodynamic numerical model. Similar to procedures followed by the aforesaid models, hydrodynamic inputs such as atmospheric pressure, the direction and speed of the wind, hydrodynamic and water quality parameters (including velocity, temperature, salinity, water level, and algae concentration), wave height and direction as well as its time period, and characteristics of suspended sediments that are all prepared by the Eulerian model are fed into the Lagrangian particle-tracking model. On the other hand, CaMPSim-3D considers the physical properties of microplastics (such as density, size, and shape) and physical processes (including the effect of positive/negative buoyancy, advection, biofouling, spatio-temporal varying diffusion, degradation, washing-off, and beaching, etc.).
In CaMPSim-3D, the physical properties of plastic particles such as density, diameter, gravity, and shape are parameterized in the calculation of the rising or settling velocity. Unlike TrackMPD, the formulation employed in CaMPSim-3D for calculating the settling velocity can be used for particles with different shapes, such as spherules, films, fragments, and fibers [98]. CaMPSim-3D calculates the settling velocity (v s ) based on force balance formulation Equations (5) and (6) provided for Newtonian fluids [98]. CaMPSim-3D can predict the formation of biofilm layers on plastic particles. In this respect, CaMPSim-3D utilizes some proposed empirical equations for estimating the biofilm thickness, e.g., Equations (7)- (9). In addition, CaMPSim-3D considers the degradation of microplastics based on empirical relationships provided in the literature, e.g., Equation (10). However, it cannot simulate both types of aggregation (i.e., hetero-and homo-aggregation) that microplastics undergo.
As with some of the aforementioned models, the Lagrangian particle-tracking model, CaMPSim-3D, can simulate not only the advection and diffusion of particles but also beaching and washing-off. It is also capable of considering the effect of buoyancy on particle movement [98]. In CaMPSim-3D, the displacement of particles due to horizontal and vertical advection is calculated by solving the first-to the fourth-order of Equations (19)- (21). In addition, the vertical displacement velocity of particles (v p ) is calculated based on the vertical component of flow velocity, and particle settling/rising velocity (v s , obtained from Equation (5)), as shown in Equation (28) [98]. Furthermore, CaMPSim-3D considers the impact of wind as a driving factor on the transport of microplastics (see Equations (20)- (22). The turbulent diffusion of plastic particles is simulated in CaMPSim-3D by using a modified relationship, as shown in Equation (29) as a function of K (diffusion coefficient, see Equation (14)), ∆t (the time step), R (a random number obtained from a normal distribution, which varies between −1 and 1), and K p (diffusion coefficient at X p as a function of t) [98].
K p is calculated using the eddy viscosity (ϑ t ) calculated by the Eulerian hydrodynamic model [98]. Moreover, CaMPSim-3D predicts resuspension and deposition as well as beaching and washing-off using the proposed relationships (see Equation (16)).
One of the advantages of CaMPSim-3D, in comparison with the aforesaid models, is its ability to employ an unstructured grid [98]. It should be mentioned that the CaMPSim-3D model reads hydrodynamic parameters from an Eulerian hydrodynamic model (TELEMAC). Although this approach has some benefits by reducing interpolation errors, it imposes more computational time for point locations in an unstructured mesh. Furthermore, CaMPSim-3D employs an inverse distance-weighted method for the interpolation (or map) of hydrodynamic inputs in space, which is implemented for each particle [98].

Results and Discussion
Despite the fact that the reviewed models have generally acceptable accuracy in simulating the physical transformational and transport processes, each has some noticeable limitations in this regard. In this section, we have tried to discuss the reviewed models' limitations in detail regarding the consideration of microplastics' physical properties and processes. It can be concluded that this section can give better insights into the current gaps in the reviewed models.
D-WAQ PART and Ichthyop assume microplastics as passive particles without inertia when they are floated at the sea surface, which is different in real ocean conditions [58,119] since the shape and geometry of microplastics have impacts on their behavior, particularly during the different physical processes that they undergo. TrackMPD considers plastic particles as both active and passive particles using different scenarios. It means TrackMPD can simulate microplastics with different shapes (spherical and cylindrical) and sizes [61]. CaMPSim-3D is more powerful in this regard in comparison with the aforesaid models. This is because it utilizes the shape factor (Ψ) which capable CaMPSim-3D in considering particles with different sphericity and circularity [98].
As stated in Section 4, D-WAQ PART and Ichthyop cannot simulate biofouling and degradation, while both TrackMPD and CaMPSim-3D are capable of simulating these processes using the same empirical relationships. It should be noted that TrackMPD and CaMPSim-3D cannot simulate biofouling and degradation together within the same simulation. On the other hand, there are some other transformational processes in the real ocean condition such as hetero-and homo-aggregation [124] and agglomeration (known as the formation of particle clusters in which weak physical interactions are dominant between particles [125]). In addition, microplastics undergo weathering, photodegradation, chemical and biological degradation, and defouling (known as the biofilm layer disintegration [44]) in real conditions. However, none of the reviewed models can simulate any of the aforementioned processes yet.
Unlike Ichthyop, D-WAQ PART, TrackMPD, and CaMPSim-3D can predict beaching and washing-off using a probability empirical theory with the same accuracy. However, actual beaching and washing-off in the real condition are highly complex processes since they depend on many other factors, such as the characteristics of winds, tides, waves, and coastline structure [58]. Moreover, microplastics may undergo stranding at the shoreline and be washed off from it repeatedly in subsequent episodes. Accordingly, the aforesaid models are currently incapable of simulating beaching and washing-off more accurately.
All of the reviewed models can predict advection and diffusion; however, their accuracy in this regard is different. All models can simulate advection and windage (except for Ichthyop). Moreover, all of them are capable of modeling horizontal and vertical diffusion, although some of them employ different methods and relationships in calculating the horizontal/vertical diffusion coefficient. However, CaMPSim-3D utilizes a modified relationship for the calculation of turbulent diffusion, which seems to be a more accurate approach in comparison with the methods employed by other models. Among the reviewed models, only CaMPSim-3D can predict resuspension. On the other hand, there are some additional important driving transport processes such as Stokes drift (as a result of wave action [58,104]), burial into the deep sediment of the waterbodies [124], and remobilization [126], which none of the reviewed models can predict.
All of the reviewed models (except for CaMPSim-3D) are incapable of reading hydrodynamic parameters from an unstructured mesh provided by an Eulerian model. It should be noted that although reading hydrodynamic data from an unstructured mesh leads to a decrease in interpolation errors, it increases the computational costs of particletracking simulations. D-WAQ PART and TrackMPD can read and then simulate spatio-temporally interpolated hydrodynamic parameters, while Ichthyop is only capable of utilizing temporally interpolated data. CaMPSim-3D can consider only spatially interpolated hydrodynamic values in the particle-tracking of plastic particles. It seems D-WAQ PART and TrackMPD utilize more accurate numerical schemes in this regard. On the other hand, the efficiency and computational speed of each of the reviewed models can be compared in systematic case studies. It should be noted that there are other numerical models such as TRACMASS and OceanParcels which have not been reviewed in this paper for the sake of briefness and to focus on models with more capabilities.
To make a general comparison among the reviewed models, the main characteristics and capabilities of the reviewed models have been summarized in Table 5. It can be seen which of the reviewed models are free for non-commercial use. Moreover, the capability of the reviewed models in terms of the particle tracking of what type and shape of particles have been compared. In addition, the main physical transformational and transport process each model can simulate are shown. * Based on the publication date of the paper which introduces the model.

Conclusions and Future Research Directions
The current study aims, as a worthwhile endeavor, to systematically review some wellknown Lagrangian numerical particle-tracking models that have been used for predicting the fate and transport of microplastics. In this regard, four Lagrangian particle-tracking models, i.e., D-WAQ PART, Ichthyop, TrackMPD, and CaMPSim-3D, have been investigated in terms of addressing physical transport processes (advection, diffusion, windage, beaching and washing-off, and resuspension) and physical transformation processes (biofouling and degradation) with consideration of the effects of the physical properties (mainly size, diameter, and shape) of microplastics on their fate and transport. Moreover, a numerical scheme of each model is described based on its spatio-temporal interpolation scheme and the structure of the adopted mesh. Accordingly, it has been shown that some of these models are more powerful than others in capturing the physical properties of particles and predicting their behavior under different physical transport and transformation processes. Some concluding remarks for the current study are summarized below: • With respect to the consideration of the physical properties of microplastics, D-WAQ PART and Ichthyop use simplified equations for calculating the settling velocity as a function of the size and density of plastic particles, while TrackMPD and CaMPSim-3D employ more accurate equations to meet the effect of particle shape as well. In other words, D-WAQ PART utilizes an equation for calculating the settling velocity of plastic particles with different shapes and types, while Ichthyop is suited for considering the settling velocity of prolate spheroids. TrackMPD can differentiate spherical and cylindrical particles in the simulation of the settling velocity. CaMPSim-3D can be applied for simulating particles with different shapes, including spherules, films, fragments, and fibers. • Among the transformation processes reviewed in this study (i.e., biofouling and degradation), although biofouling has been regarded as one of the important processes in the fate of microplastics, D-WAQ PART and Ichthyop are unable to predict the behavior of particles under the influence of this process. However, TrackMPD and CaMPSim-3D can simulate fouled particles using empirical equations. Unlike D-WAQ PART and Ichthyop, TrackMPD, and CaMPSim-3D can predict the degradation of microplastics based on the proposed relationship in the literature.

•
Four particle-tracking models employ generally universal advection-diffusion equations, while the equations for parameterizing various physical processes may be different. However, only D-WAQ PART, TrackMPD, and CaMPSim-3D consider the effect of wind-induced drift as an additional term of advection. In addition, all models are capable of meeting horizontal and vertical diffusion (dispersion). All reviewed models use spatiotemporally varying diffusion coefficients. D-WAQ PART, TrackMPD, and CaMPSim-3D can predict beaching and washing-off using some probability-based relationships. • D-WAQ PART, Ichthyop, and TrackMPD can read the extracted hydrodynamic data and simulate particle behavior and particle trajectories only on a structured mesh, while CaMPSim-3D reads hydrodynamic data from an unstructured mesh system. Moreover, the hydrodynamic data are only interpolated linearly at each time step and prepared as inputs for Ichthyop. D-WAQ PART and TrackMPD utilize the spatiotemporal interpolation of the current characteristics, while CaMPSim-3D can use only spatially interpolated data.
Given the aforesaid descriptions, there are some physical transport and transformation processes that some of the studied models are unable to simulate or exhibit less accuracy for in modeling. In fact, in order to understand the impacts of microplastics on the sustainability of water resources comprehensively, with higher levels of accuracy, we need to fully understand their behavior, which can be elucidated by using a versatile and accurate modeling system. None of the reviewed models can simulate homo-and hetero-aggregation or agglomeration. None of the four models can predict photodegradation or chemical and biological degradation as important processes in the fate of plastic particles. Wave-induced drift as an additional advection is observed only in CaMPSim-3D. It is suggested that further studies be conducted in order to expand the current knowledge of microplastics' fate and transport. There are some important processes that need to be further studied with regard to the behavior of microplastics in the real condition (e.g., hetero-and homo-aggregation, agglomeration, weathering, photodegradation, chemical and biological degradation, defouling, resuspension, Stokes drift, and remobilization). More in-depth studies, as well as more accurate numerical models in capturing the microplastics' behavior, along with significant modifications in the numerical simulation structure of the reviewed models, are recommended for the future.