Simulation of Marine Debris Path Using Mathematical Model in the Gulf of Thailand

: Marine debris is an important environmental problem that affects aquatic animals, ecosys-tems, economy, society, and humans. This research aims to simulate the path of marine debris in the Gulf of Thailand using a mathematical model that includes two models: the Oceanic Model (OCM), which is based on the Shallow Water Equations (SWE), and the Lagrangian Particle Tracking (LPT) model. The OCM is the partial derivative equation system solved by the ﬁnite difference method to satisfy the Arakawa C-grid and the splitting method. The LPT model includes the current velocity, wind velocity at 10 m above sea level, random walk term, and the buoyancy ratio of marine debris with six cases, which are 100:1, 10:1, 1:1, 0:1, 1:10, and 1:100. The current velocity from OCM is applied to the LPT model. This research uses a garbage boat that capsized near Koh Samui on 1 August 2020 as a case study. The simulated current velocity of OCM is compared with Ocean Surface Current Analyses Real-time (OSCAR) data. The Root Mean Square Error (RMSE) of u-velocity is 0.070 m/s, and that of v-velocity is 0.058 m/s. The simulation of the marine debris’s path from the LPT model demonstrates the movement to Koh Samui, Koh Taen, Koh Wang Nai, Koh Wang Nok, Koh Rap, the east coast of Nakorn Si Thammarat province, Phu Quoc Island of Vietnam and the middle of the Gulf of Thailand with the different buoyancy ratios and time durations. by the ﬁnite difference method with Arakawa C-grid and splitting methods. The mean current velocity from the OCM is compared with the mean surface current of OSCAR; the result of the OCM reveals the eddy and direction of current velocity similar to the mean surface current of OSCAR. We compare the u-velocity and v-velocity between the OCM and OSCAR using RMSE. The RMSE of u-velocity is 0.070 m/s and that of the v-velocity is 0.058 m/s. From the LPT model, the controllable and uncontrollable factor is considered by current velocity from the OCM, wind velocity, random walk, and buoyancy ratio. The buoyancy ratio with six cases (100:1, 10:1, 1:1, 0:1, 1:10, and 1:100) considers the garbage boat that capsized on 1 August 2020 at 99.84 ◦ E and 9.41 ◦ N as a case study. The simulation shows that marine debris from the starting point moves to Koh Samui, Koh Taen, Koh Wang Nai, Koh Wang Nok, Koh Rap, the east coast of Nakorn Si Thammarat province, Phu Quoc Island of Vietnam and the middle of the GoT with the various buoyancy ratios and time duration. The closest movement of marine debris from the starting point was on the west coast of Koh Samui with a buoyancy ratio of 100:1, 10:1, and 1:1, and the furthest movement of marine debris from the starting point was on the Phu Quoc island in Vietnam with a buoyancy ratio of 1:100.


Introduction
Marine debris is the most important environmental pollution such as marine debris generated by the mismanagement of waste. Marine debris also affects aquatic animals, causing injury or death that affects the ecosystem. Finally, humans and the economy are affected. For example, through aquatic animals, humans will ingest microplastics, and increased marine debris affects tourism. The Thailand Development Research Institute (TDRI) reports that Thailand has the 10th largest amount of marine debris (2020). It is classified as the amount of plastic waste that is not properly managed, a total of 1.03 million tons/year. Plastic waste leaks into the sea at a rate of about 0.41 million tons/year. Most of the waste is plastic, accounting for 12%, while foam boxes account for 10%, food packaging 8%, plastic bags 8%, glass bottles 7%, plastic bottles 7%, and straws 5% [1] After macroplastics travel to the sea, they will decompose into microplastics and move with the current from the source of microplastics to the nearshore region [2]. Abolfathi and Pearson [3,4] study the particle hydrodynamics and consider the dispersion and diffusion of micro-substances or solute dispersion using waves and currents in the nearshore area.
The Gulf of Thailand (GoT) is the shallow inlet of the South China Sea (SCS) and is located from 6 to 13.5 • N latitude and 99 to 105 • E longitude (Figure 1a). It is a semienclosed bay by the coast of the Malay Peninsula in the west and the ground of Southeast Asia in the north and east. The southeastern part of the GoT is an opening between the Laem Yuan in Vietnam and Malaysia's Kota Bharu that connects with the SCS. The average depth is about 44 m, and the deepest part of the bay is about 86 m deep.
The Oceanic Model (OCM) is the model for studying oceanic variables such as current velocity and sea surface elevation. It has been applied to study in many ways, i.e., using the OCM to detect the storm surge or study current velocity and sea surface elevation from the tsunami. Ren et al. [5,6] study the shallow water equations (SWE) to consider wave amplitude and currents in the SCS with tsunami events that occur from earthquakes in the Manila subduction zone by using the OCM based on the SWE and the second-order Godunov method to discretize the SWE. The region with a high hazard risk is considered, and tsunami information is applied to develop a tsunami warning in 5 min.
In this paper, we simulate the path of marine debris in the GoT from a starting point (Figure 1b) using two models: the OCM and the Lagrangian Particle Tracking (LPT) model. The surface current velocity of the OCM is compared with the surface current of the Ocean Surface Current Analyses Real-time (OSCAR) model and then the LPT model used the surface current velocity as initial data. let where u, v are the current velocity on the x and y-axis, respectively, ζ is the sea surface elevation, h is the total water depth (the sum of sea surface elevation and bottom topography (H)), p a is air pressure, f is the Coriolis force, ρ is the density of water, ρ a is the density of air, g is the gravity constant, W x , W y are the wind velocity on the x and y-axis, r b is the coefficient of the bottom friction, C D is the coefficient of the wind stress, |W| is the wind speed, |V| is the current speed, τ sx , τ sy are the surface stress by the wind velocity on the x and y-axis, τ bx , τ by are the bottom stress by the current velocity, F sx , F sy are the friction force on the surface, F bx , F by are the friction force on the bottom. From Guo et al. [7], the defined transformation variables are given by Multiplying Equation (1) by g and Equations (2) and (3) by Φ, we obtain the transformation of the SWE as shown in Equations (4)- (6): ∂U ∂t ∂V ∂t Equations (4)-(6) can be written in the vector form as follows: where V = (U, V), v = (u, v) and F s = F sx , F sy . We define the domain S and boundary C of S in the OCM; the initial condition is set at t = 0, and the boundary conditions are set as shown below. Initial condition: Boundary condition: where n is the unit outward normal to C. Applying the initial, boundary condition and integrating on the domain S to Equations (7) and (8), the following equation of mass conservation and energy conservation are obtained

Solving Method
The SWE are solved by the finite difference method and the splitting method [7] with the Arakawa C-grid [8].

The Finite Difference Method (FDM)
The derivative terms in the SWE can be replaced by finite difference. We use forwardfinite difference approximation for the time derivative term and space derivative term as follows where f is an unknown function, ∆t is the time interval, and ∆x, ∆y are the distance intervals on the x and y-axis, respectively.

The Splitting Method
According to Guo et al. [7], the SWE contain three feature terms: a linear term, a nonlinear term, and a physical term. Assume F ≡ U, V, or ϕ, and we consider the splitting method as follows where A 1 F is the linear term, A 2 F is the nonlinear term, and A 3 F is the physical term. The algorithm of the splitting method is as follows and applying the algorithm of the splitting method in Equations (4)- (6) with the Arakawa C-grid, we rearrange the sets of equation into three parts. The first part is called the adjustment process for the linear term as follows: The second part is called the development process for the nonlinear term as follows: The third part is called the dissipation process for the physical term as follows: where ∆t 3 = ∆t 2 = M∆t 1 , and the value of M is between 5 and 10. The data used in the OCM are the wind velocity at 10 m above sea level and surface pressure, which are the hourly reanalysis products derived from the Climate Data Store (CDS), with a horizontal resolution of 0.25 × 0.25 degrees. The topography dataset was obtained from the General Bathymetry Chart of the Ocean (GEBCO).

Spin-Up Method
The initial data for simulation of the OCM can be prepared by using the spin-up process, starting with the meteorological data at t = 0 to obtain U, V, and ϕ in the Equations (4)-(6) by using the splitting method for each time step with the grid size of 0.25 × 0.25 degrees and repeat until the total energy of OCM reaches an equilibrium state [9]. We define the total energy in discrete form as where n is the time step number, ndx, ndy are the number of points in longitude and latitude, respectively, and the condition for the spin-up process is defined as E n − E n−1 < 10 −3 .

Root Mean Square Error
Root Mean Square Error (RMSE) is used to measure the differences between prediction and observation. In this case, the RMSE is calculated between the surface current velocity of the OCM and the OSCAR. The formulation of RMSE is as follows where q is the surface current velocity from the OCM, q oscar is the surface current velocity from the OSCAR, and n is the number of data.

Lagrangian Particle Tracking Model
Based on the Lagrangian Particle Tracking (LPT) model by Carlson et al. [10], Seo et al. [11,12], and Yoon et al. [13], the authors develop the LPT model by considering controllable and uncontrollable factors. The controllable factors are the current and wind velocity, and the uncontrollable factor is the random walk term. The development of the LPT model is shown in Equation (10), where x = (x, y) is the position of marine debris, u = (u, v) is the current velocity along the x and y-axis, respectively, w 10m = w x , w y is the wind velocity at 10 m along the x and y-axis, respectively, R = R x , R y is the random number between −1 and 1 along the x and y-axis, k is the coefficient of wind pressure, K h is the horizontal diffusion (here, K h = 2.09 m 2 /s for the GoT [14], ∆t is the interval time, and W, A are the cross-section of floating marine debris to the wind direction below and above the sea surface, respectively.

Oceanic Model
The resolution of OCM is 0.25 × 0.25 degrees, and the temporal scale (time step size) is 3600 s. Before simulating the OCM, we did the spin-up process and found that the total energy satisfied the condition and was stable after 1995 h ( Figure 5). Hence, the current velocity and sea surface elevation after the total energy stable are used as initial data.
There are two monsoons that affect the GoT: the southwest monsoon from May to September and the northeast monsoon from November to February. Figures 6-9 show the direction of mean current velocity from OCM in January, April, July, and October and compare it with the mean surface current of the OSCAR. The OCM (Figure 6a) simulates the clockwise eddy at the center of the GoT, and the current direction from the SCS has a westward movement and moves northward of the GoT, which is similar to the OSCAR (Figure 6b). The OCM (Figure 7a) simulates the clockwise eddy from 9 • N to 10 • N, and the current direction from the SCS has a southwestward movement similar to the OSCAR data ( Figure 7b). The OCM (Figure 8a) simulates the clockwise eddy at the center of the GoT, similar to the OSCAR (Figure 8b). Nevertheless, the current direction from the SCS is the opposite direction. The OCM (Figure 9a) simulates the counterclockwise eddy at the connected area to the SCS, southwest, and the clockwise eddy east of the GoT, which is similar to the OSCAR (Figure 9b), but the current direction is opposite between 10 • N and 13 • N. The RMSE is used to verify the current velocity between the OCM and OSCAR during January 2020 and December 2020; we found that the RMSE of the u-velocity is 0.070 m/s and that of the v-velocity is 0.058 m/s.

Marine Debris's Path
The movement of marine debris from the LPT model is calculated by using the current velocity from the OCM, wind velocity at 10 m, and random walk term with the buoyancy ratios in six cases, which are 100:1, 10:1, 1:1, 0:1, 1:10, and 1:100 ( Figure 4). We simulate the path of marine debris with a case study; the garbage boat capsized at 99.84 • E and 9.41 • N (starting point) from 1 August 2020 to 1 August 2021 ( Figure 10).
The marine debris with a buoyancy ratio of 100:1 (refers to the plastic bottle [13] floating in the sea) moves northeastward from the starting point to the west coast of Koh Samui for 1-2 days (Figure 10a). The marine debris with a buoyancy ratio of 10:1 (refers to the plastic container [13] floating in the sea) moves northeastward from the starting point to the west coast of Koh Samui for 1-2 days (Figure 10b). The marine debris with a buoyancy ratio of 1:1 (refers to the marine debris floating on the sea and the cross-sections to the wind direction below and above the sea surface half each) moves northeastward from the starting point to the southwest coast and then moves eastward to the south coast of Koh Samui for 1-5 days (Figure 10c). The marine debris with a buoyancy ratio of 0:1 (refers to the marine debris floating under the sea surface, and it has no wind that affects the surface of marine debris) moves southward from the starting point to the Koh Wang Nai, Koh Wang Nok, and Koh Rap for 15-30 days and then moves to the east coast of Nakhon Si Thammarat province for 30-55 days (Figure 10d). The marine debris with a buoyancy ratio of 1:10 (refers to the marine debris floating on the sea surface, the marine debris has sea surface currents and wind that affect the surface of marine debris) moves eastward from the starting point to the Koh Taen and south of Koh Samui for 3-5 days (Figure 10e).The marine debris with a buoyancy ratio of 1:100 (refers to the marine debris floating on the sea surface, the marine debris has sea surface currents and wind that affect the surface of marine debris) moves southward from the starting point to the Koh Taen for 3-5 days, Ko Rap for 15-20 days and then moves northeastward to the southeast of Koh Samui for 25-30 days. Moreover, marine debris moves eastward to the Phu Quoc island in Vietnam and the middle and east of the GoT for one year (Figure 10f).

Conclusions
The marine debris path is simulated by using the OCM and LPT model; the OCM based on the SWE is used to calculate current velocity and sea surface elevation using meteorological data (wind velocity at 10 m and surface pressure). The SWE are solved by the finite difference method with Arakawa C-grid and splitting methods. The mean current velocity from the OCM is compared with the mean surface current of OSCAR; the result of the OCM reveals the eddy and direction of current velocity similar to the mean surface current of OSCAR. We compare the u-velocity and v-velocity between the OCM and OSCAR using RMSE. The RMSE of u-velocity is 0.070 m/s and that of the v-velocity is 0.058 m/s. From the LPT model, the controllable and uncontrollable factor is considered by current velocity from the OCM, wind velocity, random walk, and buoyancy ratio. The buoyancy ratio with six cases (100:1, 10:1, 1:1, 0:1, 1:10, and 1:100) considers the garbage boat that capsized on 1 August 2020 at 99.84 • E and 9.41 • N as a case study. The simulation shows that marine debris from the starting point moves to Koh Samui, Koh Taen, Koh Wang Nai, Koh Wang Nok, Koh Rap, the east coast of Nakorn Si Thammarat province, Phu Quoc Island of Vietnam and the middle of the GoT with the various buoyancy ratios and time duration. The closest movement of marine debris from the starting point was on the west coast of Koh Samui with a buoyancy ratio of 100:1, 10:1, and 1:1, and the furthest movement of marine debris from the starting point was on the Phu Quoc island in Vietnam with a buoyancy ratio of 1:100.