Numerical Prediction of the Short-Term Trajectory of Microplastic Particles in Laizhou Bay

Microplastic particles are easily captured by microorganisms and enter the food chain, which poses a threat to ecological health. These particles are abundant in coastal areas because of the influence of anthropic activities and the interaction between the sea and land. Although much research on microplastics has been done, predicting the transportation of microplastic particles in coastal zones is still a challenge. In this paper, the trajectories of microplastic particles released from four river mouths around Laizhou Bay are investigated using the lattice Boltzmann method coupled with the Lagrangian particle-tracking method, involving inter-particle and particle-wall collisions. The trajectories of particles released from four river mouths are recorded within 30 days.


Introduction
Plastic debris is frequently found in rivers, lakes and oceans. The transportation of plastic debris needs the entrainment of water flows. Although its degradation process is extremely sluggish, plastic can be broken down into small particles during its transportation. Plastic particles with a diameter less than 5 mm are recognized as microplastics [1]. These particles can remain in the soil and sediment for hundreds of years and even longer [2]. In the meantime, microplastics may be captured by marine life and then transferred into the food chain. The microplastics themselves are harmful substances the are poisonous to living beings [3]. Hence, as a persistent pollutant, the study of microplastics has provoked much attention.
Currently, many microplastics are accumulated in coastal zones, released by land sources (tourist beaches, factory sewage and rivers) and marine sources (ship transport, fishing and marine cultivation). Browne et al. [4] tested the impact of the wind and the deposition mode on the spatial distribution of microplastics in the Tamar Estuary, the U.K. Law et al. [5] studied the concentration of the floating microplastics in the surface water of the North Atlantic. Dubaish and Liebezeit [6] studied the distribution of microplastics and black carbon particles in the surface water of Jade System, Germany. Nor and Obbard [7] investigated and categorized the microplastics extracted from sediment in several mangrove habitats in Singapore. Còzar [2] provided the worldwide distribution of microplastics and predicted the accumulation sites. In the following year, Còzar [8] studied the suspended plastics in the Mediterranean Sea and evaluated the pollution level. The existing research is focused on the sampling process, categorizing, identifying the source and assessing the ecological effect of the microplastics. However, the trajectories of microplastic particles need more attention. The area with more particle trajectories enhances the chance captured by the marine creature and should not be selected as a Considering that the vertical scale is far smaller than the horizontal scale in Laizhou Bay, this two-dimensional hydrodynamic field can be simulated based on shallow water equations [27], stated as ∂h ∂t where h is the water depth; t is time; u is the flow velocity; x is the Cartesian coordinate; g is the gravitational acceleration; and F i is the force term in the i direction, determined as where z b is the bed elevation; ρ is the water density; τ wi is the wind shear stress and is governed by τ wi = ρ a C w u wi √ u wj u wj , in which ρ a is the density of air; C w is the resistance coefficient; u wi is the component of the wind velocity in i direction; τ bi is the bed shear stress and is given by where C b is the coefficient of bed friction and is expressed as C b = gn 2 b /h 1/3 and n b is the Manning's coefficient at the bed; E i is the Coriolis term, defined as in which f c = 2ωsinφ; ω is the rotation of earth, given by ω ≈ 7.3 × 10 −5 rad/s; φ is the latitude of the study area.

Microplastic Dynamics
Considering that the microplastic particles are very small and hydrophobic, these particles are treated as Lagrangian particles. They are assumed floating on the surface constantly, and the position of particles can be predicted by in which x mj is the position of microplastic particle j; Subscript m donates the microplastic particles; u mj is the velocity of microplastic particle j; u j is the flow velocity; d m is the diameter of microplastic particles and S is the ratio of microplastic particle density to fluid density. Re m is the Reynolds number of microplastic particles and C D is the drag coefficient, defined as below.

Lattice Boltzmann Method
The lattice Boltzmann method is selected as the numerical method, by virtue of its efficiency and flexibility. The lattice Boltzmann method is a completely discrete mesoscopic model that treats the fluid as fractional particles in regular lattices, providing a promising approach to complex partial differential equations such as the shallow water equations [28].
To simulate the hydrodynamic field, the lattice Boltzmann method was used, and the lattice Boltzmann equation is given by in which f α is the distribution function of fictional particles in the α direction, f eq α is the local equilibrium function and τ is the relaxation time. W α is the weight coefficient is given by Equation (10), employing the D2Q9 lattice pattern as Figure 1. F α is the force component in α direction.
The lattice speed e is determined by e = ∆x/∆t. (11) Figure 1. Sketch of D2Q9 lattice pattern.
The particle velocity vector e α is shown as The local equilibrium distribution function is expressed as The macroscopic water depth and velocity can be calculated by According to Equation (5), the position of microplastic particle at time t + ∆t can be obtained as where x t+∆t mj is the position of microplastic particle j at time t + ∆t, x t mj is the position of microplastic particle j at time t and u t mj is velocity of microplastic particle j at time t, given as

Inter-Particle Collision
To investigate the influence extent of particles collision to trajectories, the inter-particle collision model is employed. The microplastic particles are assumed floating on the water surface and their trajectories could be affected by inter-particle collisions. When the concentration of particles is low, the binary collision model is adopted [12]. In this study, the particles are rigid, and the collision occurs instantaneously. This critical step of this model is obeying the law of conservation of momentum.

Collision Identification
The total number of particles is N in the computational field. In each computational step, N(N − 1) × 2 pairs of particles should be checked for the collision identification. Take particles i and j as an example to show the identifying progress; assume particle i is advanced to particle j and the relative motion of these two particles is demonstrated in Figure 2. The distance of these two particles is R 0 in time t and R ∆t in time t + ∆t. The relative position is R c = R 0 + k(R ∆t − R 0 ), in which k is a factor of the nondimensional time scale. The occurrence of the collision requires that the following equation has two real roots, k1 and k2 (k1 < k2, 0 ≤ k1 ≤ 1).
in which d m is the diameter of microplastic particles.

Post-Collision Quantities
After collision, the velocities of the collision pair can be obtained as where u mi and u mj are precollision velocities of particles i and j. J is the impulsive force, consisting of the normal J n and the tangential J t components.
in which c is the precollision relative velocity, c = u m j − u m i , n is the normal unit velocity of the relative precollision position, n = R c /|R c |, t is the tangential unit vector of the slip velocity, and c f c is the tangential component of c.
The post-collision position of particles i and j is updated as in which x m is the post-collision position, and x m is the precollision position.

Particle-Wall Collision
The particle-wall collision frequently occurs near riverbanks and seashores. To address the complex irregular boundary shape, the lattice Boltzmann method requires that the computational domain is covered with lattices, and then, the solid boundary is converted to a denoted value stored in the center of the lattices. To begin, considering the particles are in a meshless system, it is necessary to determine the lattice where the particle located to decide whether the particle will collide with wall or not. Taking the particle collision with lower boundary as an example, except for the particle on the solid boundary, there are four lattices surrounding a particle. Select the nearest lattice position and check its solid data. A solid lattice indicates that this particle collided with wall during the time (t − ∆t, t). Second, determine the particle position (x mi , t − ∆t) at time t − ∆t and the location of this collision ((x collide , t + tw − ∆t)), as shown in Figure 3. t w is the time before the collision, defined by t w = |x mi − x collide |/u mi , where u mi is the velocity before the collision. Finally, post-collision velocity is calculated by Equations (19), (21)-(23), where m j is equal to positive infinity [29]. The post-collision position at time t is

Model Setup
Laizhou Bay (Figure 4) is one of the three bays in Bohai Sea, belonging to Shandong Province [30]. More than ten rivers, such as the Yu River and the Wei River, flow into this bay and carry sediment particles and nutrient substance into Laizhou Bay, building the habitat of natural animals and providing fishery resources [25]. However, with the increasing influence of anthropic activities, the sea water is polluted by different pollutants [26], one of which is microplastics. The degraded microplastic particles are transported with the sewage discharging into rivers and then flows into the sea.
The position and the bottom elevation of Laizhou Bay are demonstrated in Figure 5. This computational field is divided into 140 × 230 square lattices (∆x = 460 m). The north border of this area is set to the tidal boundary of 12 h period, as shown in Figure 6. The east water border is the open boundary, next to the shoreline, where the wet-dry boundary is used [31]. The most frequent wind is considered: SSW wind with an average speed of 2.5 m/s (with ρ α = 1.205 and C w = 0.0026). The initial water depth is 20 m and the velocity is 0 m/s. The time step is 4 s. The Manning coefficient of the bottom roughness is 0.01. Four groups of microplastic particles are released from the following four river mouths: the Yu River, the Di River, the Wei River and the Jiaolai River and numbered as Groups 1-4, respectively. They are released at the same time. Each group has 81 particles, equally distributed in a lattice (Figure 7). All particles are assumed floating on the water surface and the density of particles is 850 kg/m 3 . The initial velocity of these particles is 0 m/s.       This model is validated by the measured value of two sites in Laizhou Bay, Weifang Harbor and Laizhou Port [31]. In Figure 8, the results show a good agreement with the measured data, with less than 15% error, demonstrating the well accuracy of this numerical model. The hydrodynamic model is run for 24 h to obtain a well-developed flow field and then reset the computational time to 0. Considering that the hydrodynamic condition in this bay is varying with the tidal, the trajectories of microplastic particles are varying correspondingly. Hence, the release time of 324 microplastic particles at 4 river mouths are (1) 0 h, (2) 3 h, (3) 6 h and (4) 9 h, which are used to investigate the influence of the hydrodynamic condition on particle transportation.

Particles Input at t = 0 h
There are four groups of particles released from the four river mouths at t = 0 h in this section, respectively. The particles numbered 1 of each group are selected as representative, and their trajectory are recorded within 24 h. Figures 9-16 demonstrate the trajectories of these particles in the following four time periods: 0-6 h, 6-12 h, 12-18 h and 18-24 h, respectively. x is the distance in the lateral direction and y is the distance in the longitudinal direction. The units are in meters. The origin is located at (118.86 • E, 37.09 • N). At time t = 0 (Figure 9), these particles are in the initial position and then move with the water flow. During this period, the transport distance of the particle in Group 4 is the largest, while Group 1 is the smallest. Figure 10 shows the amplified diagram of trajectories of 4 particles and the velocity fields at t = 0, 3, 6 h. At t = 0 h, the flow moves toward the offshore direction, then it moves eastward at t = 3 h and southwestward at t = 6 h. The movement direction of particles is changing along with the flow direction.  Unlike the previous situation, the traveling distances in Figure 11 are shorter, and the general tendency is towards the shore during 6-12 h. The particles in Group   The cumulative collision number of particles released at t = 3 h is the most, which is 1.01 × 10 5 , followed by those released at t = 9 h (8.48 × 10 4 ). The numbers of collision of particles released at t = 0 and 6 h are only 4.05 × 10 4 and 3.32 × 10 4 .

Particles Input at t = 0 h within 30 Days
Compared to the trajectories of particles within 24 h, it can be seen in Figure 19 that the moving ranges of particles within 30 days expands obviously. The transport distances of particles in Groups 1-4 near the releasing points are 4732 m × 1224 m, 2375 m × 1709 m, 6393 m × 2466 m and 3587 m × 3105 m, respectively. In x direction, the transport distances of particles in Group 3 is the largest, while it is a relatively narrow range in Group 2 (2.7 times less than that in Group 3). Moreover, the particles in Group 4 have the widest transport distances in y direction, 2.5 times larger than that in Group 1. Generally, the transport distance of particles in Group 3 is the largest. The scope of 6393 m × 3105 m can be the minimal collect scope of microplastic particles releasing within 30 days.

Conclusions
In this work, the movement trajectories of microplastic particles in Laizhou Bay are simulated based on the lattice Boltzmann framework coupled with the Lagrangian particle-tracking method.
In addition to the current, inter-particle collisions and particle-wall collisions are also involved. After releasing, particles move forth and back with the changing tides. The trajectories of particles released at t = 0, 3, 6, 9 h are recorded with tide boundary of 12 h period. The results show that the moving range of the particles releasing at t = 9 h is the largest, which is unfavorable for the recycle process of microplastics. The simulated results indicate that these particles normally drift within the scope of 4593 m × 8242 m near the releasing point. Collisions indeed occur with as many as 3.32 × 10 4 within 24 h, which has little effect on changing trajectories of particles. Even within 30 days, the transportation of these particles is still near the seashore and their moving ranges expand to 6393 m × 3105 m from the releasing point.

Conflicts of Interest:
The authors declare no conflict of interest.