Numerical and Experimental Approach to Evaluate Microplastic Transport in Saturated Porous Media

Abstract

Under varying flow rate conditions, the transport and retention of polydisperse microplastics (MPls), with an average particle size of 16 ± 6 µm, were investigated in saturated porous media. First-order reversible and irreversible kinetic sorption models were used to describe the sorption kinetics. Sensitivity analyses provided insight into the effects of each sorption parameter. Both numerical modeling and experimental measurements were utilized to evaluate the retention rates of sand filters. The influence of flow rate on sorption was reflected in variations in the distribution coefficient (K_d), the mass transfer coefficient (β), and the irreversible sorption rate (K₁). Lower flow rates were associated with higher K_d and β values, indicating increased sorption and reduced mass transfer rates. An increase in K_d resulted in a more gradual sorption process, with a decrease in peak concentration, whereas changes in β had a comparatively smaller impact on sorption rate and peak concentration. Lower K₁ values were linked to higher peak concentrations and decreased retention efficiency. Numerical modeling revealed retention rates of 28 ± 1% at a flow rate of 31 mL min⁻¹ and 17 ± 1% at 65 mL min⁻¹. The introduction of MPls into saturated sand environments modifies the transport dynamics within the medium. Consequently, these alterations affect the hydrological characteristics of porous media, impacting groundwater quality and agricultural output. The mean absolute error (MAE) of 6% between the modeled and observed retention rates indicated a high level of accuracy. This study underscores the importance of examining retention efficiency and the accuracy of numerical models in understanding MPl transport in porous media.

1. Introduction

Plastics find extensive use across diverse sectors; however, their limited recyclability and susceptibility to degradation and weathering processes often result in fragmentation into smaller pieces. Subsequently, these plastic fragments are dispersed into natural environments through diverse transport mechanisms, including wind, runoff, and infiltration. Consequently, these plastic particles, classified as emerging pollutants, infiltrate natural ecosystems, including freshwater resources, potentially endangering organisms and disrupting natural ecological processes.

Solid plastic particles that are insoluble in water, ranging in size from 1 µm to 1000 µm, are categorized as microplastics (MPls), and those with dimensions ranging from 1 mm to 5 mm are classified as large MPls. The size classification of MPls is divided into six classes (Tabel S1) [1].

MPls have been identified in marine, terrestrial, and freshwater environments [2,3]. They have also been found in karst groundwater systems [4], alluvial sedimentary aquifers [5], and groundwater around landfill sites [6]. The filtration properties of groundwater environments play a crucial role in regulating the quality of groundwater by removing or reducing the concentration of particles and contaminants as water moves through the subsurface. The filtration mechanism of porous media is orchestrated through three key steps. The first is a transport step, characterized by a physical-hydraulic process. Following this, an attachment step ensues, driven by a physical-chemical process. The final step is detachment, wherein the retained particles are released as long as new particles or blank water continue to be supplied. These steps collectively define the dynamic processes governing the intricate filtration mechanism of porous media [7].

The physical mechanisms mainly responsible for particle transport in filter pores include Brownian diffusion, deposition, direct interception, and hydrodynamics. Diffusion assumes a critical role, especially for particles smaller than 1 micron [8]. The deposition mechanism, governed by gravity, is linked to the sedimentation rate of particles. This process guides particles through streamlines, ultimately directing them toward the grains of porous media, known as collectors. Gravity-induced deposition holds particular significance for particles larger than 1 micron, where the density of particles becomes pivotal. Interception occurs as particles traverse streamlines closely to the collector, resulting in attachment. The hydrodynamic process, influenced by the rotation and movement of particles across streamlines, is intricately linked to the shape of the particles and their interaction with the fluid [8]. For particles exceeding 10 µm in size, hydrodynamic, gravitational, and inertial forces play more significant roles, in contrast to smaller particles that are primarily influenced by physicochemical forces like electrical double-layer force, Van der Waals attraction forces, and Brownian diffusion [9].

The adaptation of approaches used for solute transport to colloidal solution transport is a common practice [10]. In the literature, there are models that integrate the classical advection-dispersion equation with various types of kinetic attachment models. The majority of these transport models presume that the attachment process adheres to either one-site kinetics, involving both forward and reverse terms [11], or two-site kinetics, encompassing reversible and irreversible attachment mechanisms [12]. In a kinetic model, the solute transport equation describes the rate at which the solute is adsorbed onto the solid surface and desorbed from it [13].

In recent studies, the analysis of MPl transport in porous media has predominantly involved the utilization of Class 1 (1 μm to 5 μm) MPIs, as outlined in Table S1 [14,15,16]. Hou et al. [17] focused on Class 3 (10 μm to 50 μm) MPls. They revealed that increased particle size in the porous medium, along with lower ionic strength and the presence of fulvic acid, can potentially augment the movement of MPIs (40–48 µm) within porous media. Conversely, smaller-sized porous media and lower velocities demonstrate superior filtration performance. The highest mass recovery of MPls was determined to be 18%. Another study involving MPIs (200–500 µm), which correspond to Class 5 (100 µm to 500 µm), found that velocity is a significant parameter influencing filtration [18]. Two models, namely the attachment model and the attachment–detachment model, were formulated. It was found that the attachment–detachment model, incorporating a time-dependent detachment coefficient, better predicted the transport of MPls (10 µm and 20 µm) within the range of Class 3 (10 µm to 50 µm) compared to the traditional attachment model [19]. This observation underscores the importance of MPl desorption in the transport process. The calculated total mass recovery under various conditions ranged from 0.08% to 35.9% [19].

Significantly, there exists a notable research gap concerning the filtration efficacy of porous media for non-reference MPls and their accompanying transport dynamics. These MPl particles, distinct from reference materials, are presumed to exhibit more shape and size distributions. The size and shape of the MPls used were not commonly used in existing studies. The main features that distinguish this study from other studies are the properties of the microplastics used and the modeling approach. The primary objective of this study is to investigate the filtration capacity of porous media during the transport of polyethylene-type, polydisperse, Class 3-sized MPl particles within a saturated quartz sand-type porous medium. The investigation utilized numerical models derived from flow and transport equations. This encompassed the computation of sorption parameters using both reversible and irreversible kinetic sorption models, determination of retention efficiencies, and examination of the impact of flow rate on the transport process.

2. Materials and Methods

2.1. Materials

The polyethylene (PE)-type MPls, in the form of fine powder with a density of 0.91 g cm⁻³, were purchased from Goonvean Fibres (Devon, UK). The PE particles are polydisperse, and their size measurements were performed under optical microscopy, with an average particle size of 16 ± 6 µm (Figure 1 and

Table 1. Statistics of MPl size distribution (All values are given in µm).

Size of Population	Max	Min	Median	Mean	Standard Deviation	Quantile 1	Quantile 3
400	77	2	16	16	6	13	19

). The graph in Figure 1a shows that the mean value remains within the range dominated by the overall size distribution of MPls. The shape of MPls can be described as fragmental (Figure 1b). A commercially available dish detergent containing anionic surfactant was employed to generate a uniform suspension of MPls. Injection concentrations were prepared with MPls and surfactant in a ratio of 1:1, with a concentration of 2 g L⁻¹. Prior to reaching the 1:1 ratio, lower ratios were attempted, but successful dispersion could not be achieved. Additionally, the 1:1 ratio was favored in the literature [20].

A stock concentration of 25 g L⁻¹ of rhodamine WT (Cole-Parmer) fluorescent dye was utilized in tracer experiments to assess the hydrodynamic characteristics of the porous medium. The primary concentration was diluted to 400 µg L⁻¹ to perform dye tracer experiments right before starting the experiment to avoid the direct influence of sunlight.

Quartz sand was kindly provided by Esan Eczacıbaşı Industrial Raw Materials Co. (Istanbul, Türkiye) with a density of 2.60 g mL⁻¹ and bulk density of 1.33 g mL⁻¹. The determination of particle size was performed through sieve analysis, resulting in a D₅₀ value of 2.9 mm (Figure 2a). A notable proportion, specifically 84%, of the sand exhibited retention on a 2 mm mesh-sized sieve (Figure 2b). The hydraulic conductivity of the sand was determined with a constant head permeameter and found to be 11.0 ± 0.3 cm min⁻¹. The effective porosity (θ) was calculated as 0.49 ± 0.009 by fluorescent dye tracer experiments and numerical modeling. Quartz sand was preferred because it is a natural material and chemically inert.

2.2. Methods

Transport experiments and constant head permeameter experiments were conducted in an acrylic column with 100 cm height and 7 cm inner diameter (Figure S1). The column was packed with sand which was washed and overnight dried at 105 °C. During the packaging process, water was pumped from bottom to top, and the column was gently tapped to prevent possible air bubbles.

The experiments began with constant head permeameter tests to determine the hydraulic conductivity, followed by dye tracer experiments to evaluate the hydrodynamic properties of the porous media, and ended with MPls tracer experiments aimed at determining the sorption parameters and retention efficiencies of the porous media. This sequence was consistently followed for each column packing, and the porous media were replaced for a new set of tests after the completion of the series of experiments. In total, four sets of experiments were performed.

2.2.1. Permeameter Experiments

Following the column packing process, preceding each dye tracer test, constant head permeameter experiments were conducted. The hydraulic conductivity of the packed sand was determined with Dacy’s Law [21], as follows:

$$K=-{\displaystyle \frac{Q}{A(dh/ds)}}$$

(1)

where K is the hydraulic conductivity (L/T), Q is the discharge (L³/T), A is the cross-sectional area of the column (L²), and dh/ds is the hydraulic gradient. The hydraulic gradient was altered three times, and the resulting average values were used in the numerical models.

2.2.2. Transport Experiments

After the permeameter experiments, the column underwent thorough flushing with distilled water, reaching a state of saturation characterized by the effluent exhibiting negligible turbidity measurements.

The volume of one pore volume in each sand-packed column varies slightly depending on the sand height and porosity. A 1.8 pore volume of fluorescent dye tracer was continuously injected, followed by 2 pore volumes for the 65 ± 0.3 mL min⁻¹ flow rate, and 5 pore volumes for the 31 ± 0.4 mL min⁻¹ flow rate of distilled water to flush the column. This procedure was used for injection of 2 g L⁻¹ MPl suspension. Each flow rate condition was replicated. The concentrations of the collected effluents were measured with a fluorometer for the dye experiments and with a turbidimeter for the MPl experiments. These measurements were plotted as breakthrough curves (BTCs). Detailed information is presented in the Supplementary Information, Section S1.

The flow rates of 31 mL min⁻¹ and 65 mL min⁻¹ correspond to linear velocities (v) of 0.4 cm min⁻¹ and 0.8 cm min⁻¹, respectively. Groundwater moves from higher to lower elevations and from higher pressure locations to lower pressure locations. Typically, this movement is quite slow, ranging from less than 0.3 m per day to a few tens of meters per day. For karst aquifers, which are not considered porous media, it can be faster [22]. The typical groundwater velocity in a sandy or gravelly aquifer may range from approximately 0.01 to 1 cm per min [23].

There is nearly a twofold difference between the two employed flow rates, allowing for the analysis of their effects on the experiments. When the flow rate is too high, advective transport becomes dominant, overshadowing the dispersion and sorption effects. Therefore, a maximum flow rate of 65 mL min⁻¹ was chosen. For the selection of the lowest flow rate, the reliability of operational conditions was considered, particularly focusing on the duration of the experiment, which was set to six hours for a flow rate of 31 mL min⁻¹. These selected flow rates align with the general conditions of groundwater flow.

For the tracing test simulations, MODFLOW and MT3DMS packages were used jointly. MODFLOW is a block-centered finite-difference flow model. MT3DMS (where MT3D stands for the Modular 3-Dimensional Transport model, and MS denotes Multi-Species) is a transport model that allows modeling of advection, dispersion, and chemical reactions within groundwater systems. MT3DMS is commonly utilized alongside MODFLOW [24,25]. Following the development and calibration of the flow model with MODFLOW, the required data for the transport model were retrieved accordingly and used in MT3DMS. The numerical models developed in this study were created with Groundwater Vistas Version 7 software (student license) (Leesport, PA, US), which provides user interfaces for MODFLOW and MT3DMS [25]. The transport of fluorescent dye was modeled using advection-dispersion equations incorporating first-order reversible kinetic (non-equilibrium) sorption. The transport of MPls was modeled with advection-dispersion equations incorporating first-order reversible kinetic (non-equilibrium) sorption and first-order irreversible sorption terms.

The general partial differential equation governing the fate and transport of contaminants in transient one-dimensional groundwater flow systems can be expressed as follows [25]:

$${\displaystyle \frac{\partial \left(\theta C\right)}{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left(\theta D{\displaystyle \frac{\partial C}{\partial x}}\right)-{\displaystyle \frac{\partial }{\partial x}}\left(\theta vC\right)+q{C}_s+\mathrm{\Sigma }{R}_n$$

(2)

where θ is the porosity (-); C is the concentration (ML⁻³); t is time (T); x is distance (L); D is the hydrodynamic dispersion coefficient (L²T⁻¹); v is the linear pore water velocity (LT⁻¹); q is the volumetric flow rate per unit volume of aquifer, representing fluid sources (positive) and sinks (negative) (T⁻¹); Cs is the concentration of the source or sink flux (ML⁻³); and ΣRn is the chemical reaction term (ML⁻³T⁻¹). Equation (2) expresses all transport processes that can be handled by the MT3DMS package in homogeneous or heterogeneous porous media. In the present study, transport is considered in homogeneous porous media where the porosity (θ) and the dispersion coefficient (D) are constant; thus, the derivatives ∂θ/∂x and ∂D/∂x are null.

In Equation (2), the term qC_s is a sink/source term. Considering the configuration of the experimental column used for the experiments in this study, this term does not play any role.

In instances where the local equilibrium cannot be attained, it is posited that the sorption process may be effectively described through a reversible kinetic reaction of the first order, as follows [25]:

$${\mathrm{\rho }}_b{\displaystyle \frac{\partial C^{\prime }}{\partial t}}=\beta (C-{\displaystyle \frac{C^{\prime }}{K_d}})$$

(3)

where ρb is the bulk density (M/L³), β is the first-order mass transfer rate between the dissolved (C) and sorbed (C’) phases (T⁻¹), and K_d is the distribution coefficient (L³M⁻¹).

The bulk density ρb is calculated from the density ρ of sand grains (ML⁻³), which is calculated using a pycnometer. Then, the bulk density is calculated as follows [21]:

$${\mathrm{\rho }}_b=(1-\mathrm{\theta })$$

(4)

In some situations, once the solute is sorbed on the solid phase, it cannot be desorbed. The reaction is irreversible and leads to a mass loss of the dissolved phase. This process can be described through a first-order irreversible kinetic sorption model with the following equation [13]:

$${\displaystyle \frac{\partial C^{\prime }}{\partial t}}=K_{1}C$$

(5)

where K₁ (T⁻¹) is the first-order irreversible kinetic sorption rate constant.

The first-order irreversible kinetic sorption model was used by Lee et al. [26] to analyze the migration of toluene contaminant plume in 3D flow conditions. This kinetic model was also applied to simulate radioactive decay or biodegradation [27]. In such cases, the constant rate (K₁) is related to the half-life (t_1/2) of the materials, as follows:

$$K_1={\displaystyle \frac{ln2}{t_{1/2}}}$$

(6)

where t_1/2 conventionally denotes the half-life (T), which was employed to signify the retained mass rather than its conventional application to represent the temporal duration associated with the decay of a substance by half.

Rhodamine WT may not behave entirely as a conservative tracer, especially when applied to the quartz sand used [28,29]. As a result, a numerical solution incorporating sorption was selected as a more suitable option, improving the alignment of calculated BTCs with observed BTCs. Through fluorescent dye experiments, the kinematic porosity (θ), dispersivity (α), distribution coefficient for the sorbed phase (K_d), and first-order mass transfer rate (β) values that characterize the properties of porous media were estimated. The obtained θ and α values were subsequently utilized as inputs for analyzing the MPl experiments. In the MPl experiments, the parameters under investigation included K_d, β, and K₁. The identification of optimal values was performed through a trial-and-error approach. The parameters K_d, β, and K₁ were adjusted iteratively to align the calculated curve with the observed curve. Following each parameter adjustment, the residual sum of squares (RSS) error was scrutinized (Table S2).

2.2.3. Retention Efficiencies

The quantity of MPls transported through the porous medium was assessed through numerical model mass calculations and by measuring the mass of particles at the end of experiments. The numerical model evaluation involved accounting for mass released from storage as a result of a decrease in sorbed concentration and mass accumulation in storage due to an increase in sorbed concentration, incorporating considerations for first-order irreversible reactions.

Experimental measurements of retention efficiencies were conducted by collecting 7.3 ± 0.7 L of transported effluent at the end of experiment and measuring the concentration in the collected effluent. The concentration of accumulated particles was determined using a turbidimeter and was converted to concentration units. Subsequently, the transported quantity was subtracted from the total injected amount, thereby yielding insight into the retention rates of MPls in porous media. The retention efficiencies were calculated with Equation (7), as follows:

$$\mathrm{R}={\displaystyle \frac{(m_i-m_R)}{m_i}}\times 100$$

(7)

where R represents the retention efficiency (%), m_i represents the total injected mass [M] of MPls, and m_R represents the recovered mass [M].

The mean absolute error (MAE) was calculated by comparing the retention values obtained numerically with those measured experimentally, aiming to evaluate the precision of the numerical model. The MAE provides a direct measure of the average absolute deviation between predicted and actual values. It is frequently used in scenarios where all errors carry equal significance, rendering it more resilient to the influence of outliers compared to other metrics such as RMSE. Additionally, MAE shares the same unit of measurement as the original data, further enhancing its interpretability.

2.2.4. Sensitivity Analyses

The sensitivity of the model was analyzed by systematically altering the input of one parameter at a time, including the flow rate (Q), distribution coefficient (K_d), first-order mass transfer rate (β), and the first-order irreversible kinetic sorption rate K₁, to evaluate their effects on the filtration of MPls. In the sensitivity analyses, each variable was individually modified based on the predicted values of the MPl A and MPl B experiments. The influence of sorption on MPl transport was assessed based on both the MPl particle transport experiments and the dye (solute) experiments. The dye experiments served as references.

3. Results

3.1. Permeameter Experiments

In this study, a constant head permeameter and a dye experiment were conducted under consistent conditions for each MPl experiment. This allowed for the observation of the transport behavior of MPls in comparison to the solute. Hydraulic conductivity values (K, cm min⁻¹) were determined with Equation (1) (

Table 2. Hydraulic conductivity (K) values corresponding to each experimental series.

Experiment	K (cm min−1)
A	11.3 ± 0.4
A.1	10.7 ± 0.1
B	10.9 ± 0.1
B.1	11.3 ± 0.3

) for each packed sand column using constant-head permeameter tests. These K values were then used in the numerical models. Minor disparities in hydraulic conductivity values are evident, likely attributed to the utilization of freshly packed sand for each experiment. Variations in the arrangement of sand grains, discrepancies in grain sizes, and the presence of angular shapes could account for these observed differences.

3.2. Transport Experiments

The observed and computed BTCs are shown in Figure 3a, b for 31 mL min⁻¹ and 65 mL min⁻¹, respectively. The values for α, θ, β, and K_d, which were obtained through modeling and calibrated with the dye tracer experiments, are presented in

Table 3. Dye experiment results. θ: porosity, α: dispersivity, K_d: distribution coefficient, β: mass transfer rate.

Experiment	Q (mL min−1)	θ (%)	α (cm)	Kd (cm3 g−1)	β (min−1)
DYE A	31	0.49	0.1	0.08	0.0008
DYE A.1	30	0.49	0.1	0.08	0.0010
DYE B	65	0.49	0.1	0.20	0.0025
DYE B.1	64	0.47	0.1	0.05	0.0025

. The average values deduced from these calibrations are as follows: dispersivity = 0.1, porosity = 0.49 ± 0.01, K_d = 0.08 mL g⁻¹ for a flow rate of 31 mL min⁻¹, and K_d = 0.125 ± 0.075 mL g⁻¹ for a flow rate of 60 mL min⁻¹. For β values, a slight increase was noted with increasing flow rate, as follows: β = 0.0009 ± 0.0001 min⁻¹ for a flow rate of 31 mL min⁻¹ and 0.0025 min⁻¹ for a flow rate of 65 mL min⁻¹.

Subsequently, α and θ were taken as constant input parameters for modeling the corresponding MPl experiments. Dye concentrations swiftly reach a relative concentration close to unity after the arrival of the solute. When assessing the experiments conducted with MPls, there is an observable gradual increase in particle concentration. This gradual increase is more prominent at the 31 mL min⁻¹ flow rate (Figure 3).

Sorption, which is expected between MPl particles and sand grains, causes the MPls to appear later than expected at the column outlet. If the displacement of solute is held back within the column due to some chemical or physical process, the BTC will shift to the right [30]. The BTCs for MPls at 31 mL min⁻¹ exhibit a shape more indicative of transport involving sorption reactions compared to the BTCs at 65 mL min⁻¹ (Figure 3). An examination of the BTCs highlights that at high flow rate, the transport of MPls would be closer to the reference transport of dye.

Sorption is characterized through the values of K_d, β, and K₁, and the calculated values are presented in

Table 4. Computed values for distribution coefficient (K_d), mass transfer rate (β), irreversible kinetic sorption rate (K₁). Q: flow rate, C/C₀ peak: ratio of MPls concentration at the peak point.

Experiment	Q (mL min−1)	Kd (mL g−1)	β (min−1)	K1 (min−1)	C/C0 peak
MPl A	31	0.2	0.05	0.006	0.78
MPl A.1	31	0.2	0.05	0.006	0.82
MPl B	65	0.05	0.02	0.008	0.86
MPl B.1	64	0.09	0.02	0.007	0.86

. At the lower flow rate, the average K_d value is 0.2 mL g⁻¹, whereas at the higher flow rate, the average K_d value is 0.07 ± 0.02 mL g⁻¹. A reduction in flow rate corresponds to an increase in the K_d value. The β values are 0.05 min⁻¹ for the lower flow rate and 0.02 min⁻¹ for the higher flow rate. The bulk density of sand is 1.33 g mL⁻¹. This value is used by the model for the calculation of β values. Increased flow rates are associated with a decrease in β values, indicative of altered mass transfer rates. The irreversible sorption of MPls in the column is characterized by K₁. The calculated average K₁ values are 0.006 ± 0.0003 min⁻¹ for a low flow rate of 31 mL min⁻¹ and 0.007 ± 0.0004 min⁻¹ for a relatively high flow rate of 65 mL min⁻¹. An increase in the flow rate results in an increase in the K₁ value (Table 4).

3.3. Retention Efficiencies

The retained mass within the porous media was determined through both numerical simulations (R₁) and experimental measurements (R₂). Analyzing the percentage of retained mass enabled the assessment of the retention efficiency of MPls within the utilized porous medium, particularly in response to variations in flow rate. For a flow rate of 31 mL min⁻¹, the average retention rates equal 28 ± 1% for R₁ and 24 ± 2% for R₂. For a flow rate of 65 mL min⁻¹, the average retention rates are 17 ± 1% for R1 and 9 ± 4% for R₂, which are presented in

Table 5. Retention rates according to flow rate (Q). R₁: Simulated retention rates; R₂: Measured retention rates.

Experiment	Q (mL min−1)	R1 (%)	R2 (%)
MPl A	31	29	26
MPl A.1	31	27	22
MPl B	65	18	13
MPl B.1	64	16	5

. A trend may be observed indicating a decrease in retention efficiency with an increase in flow rate. The standard deviation of R₂ is notably higher for the flow rate of 65 mL min⁻¹, which may be attributed to non-homogenized effluent suspension. This could be explained by the higher flow rate leading to increased flow heterogeneity.

The mean absolute error (MAE) is 6%, indicating that, on average, the model’s predicted retention rates (R₁) differ by 6% from the experimentally measured retention rates (R₂). This can provide information on the accuracy of the numerical model simulation based on experimentally measured retention rates. It is important to acknowledge that sampling homogeneously from MPl effluent reservoirs is challenging, which may contribute to slight measurement errors that need to be accounted for while considering model accuracy.

3.4. Sensitivity Analyses

The impact of sorption parameters can be uncovered through sensitivity analyses. The sensitivity analyses were conducted for flow rates of 31 mL min⁻¹ and 65 mL min⁻¹, and their corresponding BTCs are presented in Figure 4 and Figure S4, respectively.

The variations in the K_d value were modeled while keeping β fixed at 0.05 min⁻¹ and K₁ fixed at 0.006 min⁻¹. In the context of non-equilibrium transport, the phenomenon of sorption is contingent on time and occurs gradually [31]. The increase in K_d leads to a more gradual sorption over an extended period and increased retardation.

To analyze the impact of β, K_d and K₁ were held constant at values of 0.2 mL g⁻¹ and 0.006 min⁻¹, respectively. A decrease in β leads to a reduction in peak concentration. Conversely, as the β value increases, the sorption process occurs more gradually (Figure 4b).

Finally, to examine the impact of K₁, the K_d and β values were held constant at 0.17 mL g⁻¹ and 0.06 min⁻¹, respectively. A decrease in K₁ leads to an increase in the concentration peak value, while an increase in K₁ results in a decrease in the concentration peak value (Figure 4c). Since it is irreversible kinetic sorption rate, K₁ directly influences the retained mass, with a decrease in K₁ resulting in a decrease in retention efficiency (

Table 6. Impact of sorption parameters on retention efficiency. K₁: irreversible kinetic sorption rate; K_d: distribution coefficient; β: mass transfer rate.

Q (mL min−1)	Kd (mL g−1)	β (min−1)	K1 (min−1)	Retention Efficiency (%)
31	0.2	0.05	0.006	29
31	0.2	0.05	0.03	81
31	0.2	0.05	0.001	7
65	0.25	0.02	0.01	18
65	0.05	0.02	0.04	63
65	0.05	0.02	0.002	4

).

Retention efficiencies were examined by evaluating the sensitivity to flow rate changes while holding the sorption parameters (K_d, β, and K₁) constant. The used sorption values and calculated retention efficiencies are shown in

Table 7. Impact of flow rate (Q) on retention efficiency. K₁: irreversible kinetic sorption rate; K_d: distribution coefficient; β: mass transfer rate.

Kd (cm3/g)	β (min−1)	K1 (min−1)	Q (mL min−1)	Retention Efficiency (%)
0.2	0.05	0.006	155	7
0.2	0.05	0.006	65	15
0.2	0.05	0.006	31	29
0.2	0.05	0.006	6	82

. As the flow rate decreases, an increase in the retention of MPls within the porous medium is observed. This phenomenon is attributed to the higher velocities associated with increased flow rates, leading to a stronger advection process. Consequently, when flow rates are higher, advection predominates over retention, resulting in lower retention rates (Figure S5).

4. Conclusions

The study aimed to address the knowledge gap concerning the transport of polydisperse MPls, with average particle size of 16 ± 6 µm, in saturated porous media under various flow rate conditions. The hydrophobic nature of polyethylene presents a challenge for dispersing MPls in water. Surfactants play a crucial role in facilitating their dispersion in distilled water. Moreover, surfactants can be commonly found in natural environments due to being frequently employed in cleaning products.

The study explored the impact of flow rate on sorption, examining the distribution coefficient (K_d), mass transfer coefficient (β), and irreversible sorption rate (K₁). The BTCs for dye tracer experiments showed consistent behavior, while those for MPls exhibited gradual increases in particle concentration, particularly at lower flow rates, indicating sorption effects. It was found that sorption parameters changed depending on the flow rate. The higher values of K_d and β were determined at lower flow rates. This indicated that the sorption rate increased under conditions of low flow rates.

Increasing K_d led to more gradual sorption and increased retardation. A decrease in β was followed by a decrease in peak concentration and less gradual sorption. Changes in K₁ influenced the retained mass, with lower K₁ values associated with higher concentration peak values and decreased retention efficiency when the flow rate was constant. Another significant parameter affecting the retained mass in the porous medium was the flow rate. When keeping the sorption parameters constant, an increase in flow rate led to a decrease in retention efficiency. That indicated that advection became preponderant compared to retention.

While recognizing that the amounts of surfactant and injected MPls in this study may not perfectly replicate natural conditions, it is crucial to note that the study effectively addresses gaps in understanding non-reference MPls. For future perspective, investigations could focus on modifying the water properties to simulate natural systems more accurately. This may entail introducing organic materials and various types of salt ions. Additionally, reducing the concentration of injected MPl particles and considering particle–particle and particle–collector interactions could provide further valuable insight in MPI transport characteristics and porous media filtration efficiency. These conditions can contribute to an improved understanding of freshwater source protection and aid in the development of strategies to safeguard these vital resources.