Modeling and Parametric Simulation of Microplastic Transport in Groundwater Environments

: Efforts to reduce the toxic effects of microplastics (MPs) on the environment have increased globally in recent years. However, the existing models used for the simulation of contaminant transport in groundwater are meant for dissolved substances, which is not suitable for studying MPs. Therefore, in this study, the transport of MPs in a saturated porous medium was modeled by establishing governing equations. Simulations were performed using the ﬁnite element method to examine the effects of the parameters of the governing equations on the transport of MPs. The results suggest that it is necessary to reduce the diffusivity of MPs and increase the water ﬂow velocity, porosity, and ﬁrst-order attachment coefﬁcient to effectively contain this environmental hazard. From the simulation results, it can be derived that a combination of low diffusivity, fast water ﬂow velocity, and high soil porosity may reduce the amount of MPs that are leaked into groundwater environments. The modeling and simulations performed in this study provide a clear understanding of the transport phenomena of MPs with applications in combating water pollution.


Introduction
Microplastics (MPs) include primary MPs that are intentionally made with sizes of 5 mm or smaller, as well as secondary MPs that have been weathered by ultraviolet radiation, biodegradation, physical wear, or chemical oxidation [1]. It is estimated that 1.15 to 2.41 million tons of plastics flow into the ocean worldwide through streams and rivers because of low recycling rates and poor waste management [2,3]. Among these plastics, MPs can absorb toxic pollutants because of their high specific surface area and reactivity while releasing phthalate, bisphenol A, and brominated flame retardants, causing adverse effects on organisms that accidentally feed on these plastics [4]. Although the interest of the international community is increasing concerning this problem as MPs have emerged as a serious global environmental threat [5,6], public anxiety and concerns about MP pollution in terrestrial ecosystems are on the rise in Korea, as most pollution surveys have focused on marine and coastal environments [7][8][9]. In particular, the current conditions of MP pollution in groundwater have not been actively investigated in Korea. The transport of MPs in subsurface environments needs to be evaluated to address public concerns and anxiety over MP pollution, as well as to ensure consistency considering the quality of groundwater resources [10].
The European Union (EU) is preparing institutional strategies to monitor MPs in drinking water, including spring water, and they are devising preventive measures in circular economy packages, such as developing plastic recycling technologies and biodegradable plastics, classifying toxic substances in plastics, and preventing marine waste. The Swedish Environmental Protection Agency (SEPA) has aimed to reduce plastic usage, achieve efficient plastic recycling, and prevent plastic contamination from land-based sources by supporting projects concerning the sustainable use of plastics and adequate management 2 of 12 of plastic waste by non-profit organizations, public organizations, universities, and research institutes from 2020 to 2021 [11]. In the Netherlands' TRAMP research project, which extended from 2015 to 2019, relevant ministries, research institutes, and universities participated in the development of nanoplastic and MP detection technologies for freshwater environments; developed plastic transport, risk, and impact evaluation technologies; and predicted the current and future risks of plastics in freshwater environments [12]. In Germany, various studies on the transport of MPs in ecosystems, the biological impact of MPs, and natural attenuation of MPs have been underway for several years, and associated research findings have been published in academia [13][14][15][16]. A recent study has also focused on the existing water quality index and its modification using a multi-criteria decision-making method [11,[17][18][19].
The U.S. National Oceanic and Atmospheric Administration (NOAA) supports a variety of studies through collaborations with schools and research institutes related to MPs, and a recent study on MP discovery in a karst aquifer by a team at the University of Illinois was also supported by the U.S. National Science Foundation (NSF). The 5 Gyres Institute and the San Francisco Estuary Institute, supported by the Gordon and Betty Moore Foundation, investigated the current status, influx routes, and transport of MPs in the San Francisco Bay area for three years through the San Francisco Bay Microplastics Project (www.sfei.org; accessed on 12 July 2021). Adventure Scientists, a non-profit organization in the United States, surveyed the level of MP pollution in marine and freshwater systems worldwide from 2013 to 2017, and they established a database for the current status of pollution [20][21][22].
The analysis of the samples obtained from wells (<65 m) and springs in the karst terrains of the Salem Plateau and Drift area in Illinois, U.S., indicated the presence of MPs in 16 out of 17 samples, which were tentatively assumed to have flowed from septic tank effluents [23]. Foreign studies have been actively conducted on MPs in marine ecosystems and surface water (lakes, rivers, etc.), and soil-related research is becoming increasingly popular [24,25]. However, both domestic and foreign studies on MPs in groundwater are limited to the drinking water of karst terrains and groundwater origins [23,26,27]. In several countries, including the U.S., Germany, Canada, and Denmark, MPs were detected in drinking water of groundwater origin, as well as in spring water products (bottled water) [2,7,13].
Up to 30,000 pore-scale MPs (20 to 50 µm) were detected in hyporheic zones where surface water and groundwater actively mixed, and the possibility of MPs entering groundwater through these mixing sections has been reported [13]. Furthermore, the possible influx of MPs into groundwater through soil due to soil characteristics (pores and fractures), biological activity, human activity, and precipitation has been suggested in studies on the transport of MPs in soil [28,29]. Most studies on the subsurface transport of MPs are based on the transport of colloids or nanoparticles, using the traditional advection-dispersion equation for dissolved species, mainly involving indoor experiments [30,31]. Despite the growing worldwide interest in MPs, the existing research on the transport of MPs is confined to marine and terrestrial environments that only include coastal areas and surface water [5,9], with very few studies being related to actual soil [13] and groundwater environments [13,23,27]. The extensive physical and chemical properties and size of MPs that are smaller (or larger) than colloids but larger than nanoparticles are not reflected in such studies, and the application of existing models is limited [29].
As mentioned above, in the conventional evaluation of the transport of MPs, advectiondispersion equation models for evaluating the transport of colloids or nanoparticles in unsaturated and saturated media have been used. However, the sizes of MPs in real-life environments vary widely from sizes smaller to larger than colloids, and the conventional advection-dispersion equation for dissolved species is insufficient for modeling the transport of MPs. In addition, existing models such as MT3D [32] and RT3D [33] in MODFLOW [34], which can simulate the transport of contaminants in groundwater, are specifically used to evaluate the diffusion of dissolved species and may not be suitable for evaluating the transport of MPs. Therefore, in this study, the existing transport models are analyzed, and the details of their main parameters are identified.

Modeling of Microplastics Transport in Porous Media
Modeling the movement, attachment, and detachment of particles within saturated porous media is a subject that has been addressed in many areas of research [35][36][37][38][39]. Previous studies are based on the concept of mass balance, and the modeling is described using continuity equations. The same concept applies to the modeling of MPs, and some researchers have adopted colloid transport models [40][41][42][43]. Figure 1 shows the transport of MPs in a saturated porous medium. When MPs flow through a pore, with a volume of θ, within saturated soil with a volume of 1, the MPs floating in the pore are transported downstream by diffusion and advection. Some MPs inside the pores are attached to the surface of the soil particles, while some of the attached MPs become detached from the surface of the soil particles.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 3 of 12 [33] in MODFLOW [34], which can simulate the transport of contaminants in groundwater, are specifically used to evaluate the diffusion of dissolved species and may not be suitable for evaluating the transport of MPs. Therefore, in this study, the existing transport models are analyzed, and the details of their main parameters are identified.

Modeling of Microplastics Transport in Porous Media
Modeling the movement, attachment, and detachment of particles within saturated porous media is a subject that has been addressed in many areas of research [35][36][37][38][39]. Previous studies are based on the concept of mass balance, and the modeling is described using continuity equations. The same concept applies to the modeling of MPs, and some researchers have adopted colloid transport models [40][41][42][43]. Figure 1 shows the transport of MPs in a saturated porous medium. When MPs flow through a pore, with a volume of θ, within saturated soil with a volume of 1, the MPs floating in the pore are transported downstream by diffusion and advection. Some MPs inside the pores are attached to the surface of the soil particles, while some of the attached MPs become detached from the surface of the soil particles.

Microplastics Balance Equations in the Aqueous Phase and for the Solid Phase
The balance equations for MPs are as follows: Here, Ct = total microplastics (N/m 3 ) f = total microplastics flux (N/m 2 s) C = microplastics in the aqueous phase (N/m 3 ) S = microplastics of the solid phase (N/kg) θ = volumetric water content (-) ρ = soil bulk density (kg/m 3 ), D = diffusivity (m 2 /s) ν = average pore water velocity (m/s) By substituting Equations (2) and (3) into Equation (1) for one-dimensional flow, we obtain:

Microplastics Balance Equations in the Aqueous Phase and for the Solid Phase
The balance equations for MPs are as follows: Here, Ct = total microplastics (N/m 3 ) f = total microplastics flux (N/m 2 s) C = microplastics in the aqueous phase (N/m 3 ) S = microplastics of the solid phase (N/kg) θ = volumetric water content (-) ρ = soil bulk density (kg/m 3 ), D = diffusivity (m 2 /s) ν = average pore water velocity (m/s) By substituting Equations (2) and (3) into Equation (1) for one-dimensional flow, we obtain: MPs in the solid phase can be expressed using the following attachment-detachment model for colloids [38]: Here, k a = the first-order attachment coefficient (1/s) 2d 50 ηαv η = collector efficiency (-) α = microplastic sticking efficiency (-) d 50 = the median porous medium grain diameter (m) k d = the first-order detachment coefficient (1/s) The k d term is usually ignored because it is very small compared to k a [10], and Equation (5) can be simplified as follows: In this study, the transport of MPs was modeled by writing and executing finite element method (FEM) codes to solve the governing equations shown in Equations (4) and (6). Finite element analysis procedures are widely used to determine the solutions of balance equations. Galerkin's residual method is employed to determine the governing partial differential equations of Equations (4) and (6). The detailed development of this procedure is described in [44]. Figure 2 shows an example of an analysis that examines the effects of the main parameters manifested in the aforementioned governing equations for this system. MPs are supplied between 2 and 10 min upstream of a 2 m long water channel, which is assumed to exhibit a one-dimensional flow, to maintain a concentration of C = 10,000 particles/m 3 . A concentration of C = 0 particles/m 3 was maintained at the downstream boundary. The concentration condition at the downstream boundary was appropriate because the distance of 2 m was long enough to prevent MPs from reaching the downstream area during the analysis period. Table 1 shows the parameter values used in this simulation, and Table 2 lists the analysis conditions for each parameter. parameters manifested in the aforementioned governing equations for this system. MPs are supplied between 2 and 10 min upstream of a 2 m long water channel, which is assumed to exhibit a one-dimensional flow, to maintain a concentration of C = 10,000 particles/m 3 . A concentration of C = 0 particles/m 3 was maintained at the downstream boundary. The concentration condition at the downstream boundary was appropriate because the distance of 2 m was long enough to prevent MPs from reaching the downstream area during the analysis period.

S0
Reference case See Table 1. Figure 3 shows the simulation results for the transport of MPs in the reference case. MPs are supplied upstream at the 2 min mark and transported downstream over time, as shown in Figure 3a. As the concentration C increases at each location, so does the concentration S (Figure 3c). After the upstream supply of MPs stops at the 10 min mark, the concentration C decreases at each location. Here, Equation (6), which concerns only the occurrence of attachment, is used instead of Equation (5), which concerns both attachment and detachment. The MPs are therefore not detached from the surface of the soil once they become attached (Figure 3d). Figure 4 shows the simulation results for the transport of MPs with increased MP diffusivity. As the transport rate of MPs increased, the downstream concentration increased more rapidly (Figure 4a,b). The attachment rate and amount of attached MPs also increased (Figure 4c,d).

Results and Discussion
These simulation results were obtained using a constant upstream concentration. Therefore, these results correspond to the conditions in which sufficient MPs are supplied upstream. Under the boundary condition in which the amount of supplied MPs is constant, the surface values at the upstream boundary (x = 0) may be lower than C = 10,000 and S = 1000. The boundary condition in this case can be modeled using the time-dependent boundary condition of C. However, to simulate situations such as the accidental bulk leakage of MPs into groundwater environments, the current boundary condition of a constant upstream concentration of C is more appropriate. Therefore, the diffusivity is an important parameter for MP transport in groundwater environments because a higher diffusivity means greater and faster leakage into groundwater environments.        Figure 5 shows the simulation results for the transport of MPs with an increased water flow velocity through a pore space. A noticeable difference compared to the reference case, S0, is that the concentration S almost doubles. This is attributed to the increased influx of MPs with a higher flow velocity, as well as a sufficient upstream supply of MPs (see Equation (3)). The rate at which S grows over time increases with an increase in ka, as shown in Equation (6). In summary, an increase in the flow velocity of water leads to an increase in the amount of attached MPs. Considering water management and leakage accidents, the pressurization of the contaminated soil surface area may be utilized such that the leaked MPs may be contained within the soil. The contaminated soil could then be removed later. Many geotechnical engineering solutions for achieving such pressurization have been developed for application [45]. Figure 5 shows the simulation results for the transport of MPs with an increased water flow velocity through a pore space. A noticeable difference compared to the reference case, S0, is that the concentration S almost doubles. This is attributed to the increased influx of MPs with a higher flow velocity, as well as a sufficient upstream supply of MPs (see Equation (3)). The rate at which S grows over time increases with an increase in ka, as shown in Equation (6). In summary, an increase in the flow velocity of water leads to an increase in the amount of attached MPs. Considering water management and leakage accidents, the pressurization of the contaminated soil surface area may be utilized such that the leaked MPs may be contained within the soil. The contaminated soil could then be removed later. Many geotechnical engineering solutions for achieving such pressurization have been developed for application [45].  Figure 6 shows the simulation results for the transport of MPs with an increased porosity. Similar to the previous case, wherein an increased water flow velocity is investigated, the concentration of the attached MPs, S, increased noticeably compared to the reference case, S0. However, the amount of attached MPs increased more rapidly in this case, and the final amount of attached MPs was larger for an increased water flow velocity. It is not possible to remove soil into which MPs have leaked while a leakage accident is occurring. However, it may be a good precaution to initially construct a base soil with a large porosity for industries affected by MPs. Figure 7 shows the simulation results for the transport of MPs with a decreased firstorder attachment coefficient, ka. The simulation results were obtained by setting d50 = 2 d50,0. The concentration of attached MPs, S, decreased as the attachment performance in Equation (6) decreased. A combination of decreased S and a constant value of C indicates that a relatively larger amount of floating MPs existed in this case, and therefore adverse effects in terms of controlling the leakage of MPs were present.   Figure 6 shows the simulation results for the transport of MPs with an increased porosity. Similar to the previous case, wherein an increased water flow velocity is investigated, the concentration of the attached MPs, S, increased noticeably compared to the reference case, S0. However, the amount of attached MPs increased more rapidly in this case, and the final amount of attached MPs was larger for an increased water flow velocity. It is not possible to remove soil into which MPs have leaked while a leakage accident is occurring. However, it may be a good precaution to initially construct a base soil with a large porosity for industries affected by MPs. Figure 7 shows the simulation results for the transport of MPs with a decreased firstorder attachment coefficient, ka. The simulation results were obtained by setting d 50 = 2 d 50,0 . The concentration of attached MPs, S, decreased as the attachment performance in Equation (6) decreased. A combination of decreased S and a constant value of C indicates that a relatively larger amount of floating MPs existed in this case, and therefore adverse effects in terms of controlling the leakage of MPs were present.      Figure 8 shows the simulation results for the transport of MPs with an increased soil bulk density. Substituting Equation (6) into Equation (4) results in Equation (7), in which the soil bulk density term is eliminated, and its influence does not manifest. Figure 8 shows the simulation results for the transport of MPs with an increased soil bulk density. Substituting Equation (6) into Equation (4) results in Equation (7), in which the soil bulk density term is eliminated, and its influence does not manifest.
The simulation results may be summarized as follows. When the diffusivity is varied (S1), the change in the concentration of attached MPs, S, is small, while the change in the concentration of MPs floating in the water, C, is more noticeable. Furthermore, changes in the soil bulk density (S5) do not affect the transport of MPs.
The other examined cases (S2, S3 and S4) showed significant differences in terms of the amount of attachment. From the above findings, it can be derived that a combination of low diffusivity, fast water flow velocity, and high soil porosity may reduce the amount of MPs being leaked into groundwater environments.

Conclusions
Governing equations were established to model the transport of MPs in a saturated porous medium, and FEM codes were written and executed to perform simulations. Based on these simulations, the effects resulting from varying the parameters of the governing equations for the transport of MPs were examined. To reduce the influx of MPs into the ocean, the simulation results indicate that it is necessary to reduce the diffusivity of MPs and increase the water flow velocity, porosity, and first-order attachment coefficient. The properties of MP materials and pore water in groundwater sources are key elements in decreasing diffusivity, but these characteristics are virtually impossible to change.
However, it was found that the application of a pressurization system can increase the flow velocity of water in areas identified as sources of MP influx. In addition, highly contaminated soil at the top of these areas can be replaced with high-porosity soil to  The simulation results may be summarized as follows. When the diffusivity is varied (S1), the change in the concentration of attached MPs, S, is small, while the change in the concentration of MPs floating in the water, C, is more noticeable. Furthermore, changes in the soil bulk density (S5) do not affect the transport of MPs.
The other examined cases (S2, S3 and S4) showed significant differences in terms of the amount of attachment. From the above findings, it can be derived that a combination of low diffusivity, fast water flow velocity, and high soil porosity may reduce the amount of MPs being leaked into groundwater environments.

Conclusions
Governing equations were established to model the transport of MPs in a saturated porous medium, and FEM codes were written and executed to perform simulations. Based on these simulations, the effects resulting from varying the parameters of the governing equations for the transport of MPs were examined. To reduce the influx of MPs into the ocean, the simulation results indicate that it is necessary to reduce the diffusivity of MPs and increase the water flow velocity, porosity, and first-order attachment coefficient. The properties of MP materials and pore water in groundwater sources are key elements in decreasing diffusivity, but these characteristics are virtually impossible to change.
However, it was found that the application of a pressurization system can increase the flow velocity of water in areas identified as sources of MP influx. In addition, highly contaminated soil at the top of these areas can be replaced with high-porosity soil to prepare for future MP influxes. Lastly, to increase the attachment coefficient of MPs, it may be necessary to apply an electrical treatment to the groundwater system of a contaminated area. To reduce the influx of various MPs into the ocean, the establishment of governing equations for the transport of MPs in groundwater environments and subsequent analysis of the effects of the main parameters of these equations that were provided by this study may aid in improving the safety of these systems.