CFD Analysis of Microplastic Transport over the Slopes

: Microplastics, ubiquitous in our environment, are significantly impacted by the hydrodynamic conditions around them. This study utilizes CFD to explore how various breaker types influence the dispersion and accumulation of microplastics in nearshore areas. A special focus is given to the impact of wave dynamics and particle size, particularly on buoyant microplastics in spilling breakers. It was discovered that spilling breakers, common on gently sloping seabeds, encourage broad dispersion of microplastics, notably for smaller-sized particles. Plunging breakers exhibit a similar pattern but with less dispersion and an initial forward movement of neutral and heavy particles. Surging breakers feature minimal dispersion and a distinct oscillatory motion. It has been observed that medium-sized particles with a 1 mm diameter in this work exhibit the most substantial forward movement, likely due to an optimal balance between inertia and viscosity, enabling an effective response to wave momentum. Larger particles, influenced mainly by inertia, tend to show less dispersion and advection. Meanwhile, smaller particles, more affected by viscosity, demonstrate greater dispersion, interacting extensively with wave-induced turbulence. This study reveals the significance of inertia in the behavior of microplastics over slopes, emphasizing the importance of considering inertial effects for precise modeling of microplastic movement in nearshore areas.


Introduction
Microplastics, often smaller than 5 mm, have become widespread in marine and coastal ecosystems.Their presence is influenced by a complex interplay of their physical properties and diverse oceanic processes [1][2][3].Once in the ocean, these plastics are subject to changes in size, density, and shape due to various physical, chemical, and biological factors, impacting their movement and distribution [4][5][6][7].The ubiquity and diminutive size of microplastics pose considerable risks to marine life and, by extension, human health [8][9][10].Notably, the accumulation of microplastics in nearshore zones raises concerns due to the concentration of waste from sources like river runoff and human activities, potentially leading to ingestion by marine organisms and subsequent ecosystem and human health impacts [11].
Nearshore areas act as dynamic reservoirs for microplastics [12][13][14], with wave action adding to the complexity of their behavior.A significant portion of plastic waste is subjected to wave-driven processes, influencing the dispersal and accumulation of microplastics [15][16][17][18].The dynamics of waves, coupled with particle buoyancy and coastal topography, play a crucial role in directing microplastics either onshore or offshore [19].Moreover, wave-induced mass transport, known as Stokes drift, intensifies near coastlines [20][21][22], and wave-induced turbulence from breaking waves can alter microplastic distribution beneath the surface [23,24].The interplay of wind, wave shoaling, currents, and Langmuir circulation, all of which intensify wave nonlinearity in shallow waters, further directs microplastics toward the coast [25][26][27][28].This intricate interaction between hydrodynamics and microplastic movement underlines the need for systematic research into their distribution and transport in varied marine settings.
The beach slope significantly influences how waves interact with microplastics, affecting their dispersion and accumulation [29].This interaction can be categorized into three primary wave breaker types-spilling, plunging, and surging-each associated with different levels of nonlinearity due to varying slopes [30].Spilling breakers, found on gentle slopes, evenly distribute energy, leading to the uniform dispersion of sediments and microplastics.Plunging breakers, typical of steeper slopes, release energy more suddenly, creating concentrated turbulence and localized microplastic accumulation.In contrast, surging breakers on steep slopes or against structural barriers concentrate wave energy into intense bands, driving longshore currents that laterally disperse microplastics.Despite the potential effects of slope on nearshore dynamics and corresponding microplastic behaviors, research on how varying beach slopes influence the behavior of contaminants like microplastics remains limited.
Computational Fluid Dynamics (CFD) simulations are instrumental in elucidating the interactions between waves and microplastics, providing insights into the complexities of their behavior and aiding in the development of pollution control strategies.These simulations offer a detailed understanding of nearshore hydrodynamics and the impact of beach slope variations on surf zone turbulence [31], which also affects sediment transport and deposition.It has been noted that particle size and buoyancy significantly impact dispersion, with smaller and neutrally buoyant particles experiencing greater dispersion [32].The CFD has been used to evaluate the effects of physical parameters of microplastics, such as shape, density, and size, on their settling velocity [33], and the turbulence generated by waves considerably affects this velocity [34].The behavior of microplastics within fluvial systems was also investigated using CFD to understand their infiltration and transport dynamics in the hyporheic zone [35].Although there have been advancements in predicting transport and settling velocity in various marine environments, an extensive analysis that encompasses the effects of different breaker types on microplastic dispersion and accumulation across varying beach slopes has yet to be fully explored.
In this study, we employed 2-dimensional CFD to investigate how particle size and density influence their distribution patterns in relation to different breaker types.We analyzed background Eulerian data, including mean current, Reynolds stress, and streamlines, and examined Lagrangian data using statistical tools such as calculations of moments of particle distribution and kinetic energy spectra of particle velocities.Our discussion focuses on buoyant particle behavior in spilling breakers pertinent to contaminant dispersal in gently sloping beach areas.Additionally, we demonstrate how differences in particle inertia can lead to discrepancies between Eulerian and Lagrangian velocities, thereby influencing microplastic behavior.

Numerical Wave Tank
This study aimed to investigate the behavior of microplastics in the nearshore area, where the microplastics interact with sloped beaches in wave-dominant environments, by utilizing the CFD wave flume (ANSYS FLUENT 2020 R2).The length and depth of the numerical flume were 12.0 m and 0.8 m, respectively, with a still water depth (h) of 0.4 m.The detailed dimensions of the wave flume and the reference coordinate system, with the origin set at (x = 0, y = 0), are depicted in Figure 1.The x-coordinate signifies the cross-shore direction, pointing rightward.conducted using grid sizes of 0.01, 0.02, 0.03, and 0.05 m with 99,060, 75,564, 11,745, and 4410 nodes, respectively, showed that the Root Mean Square Error (RMSE) between the mean velocity fields of the 0.01 m case and the others decreased significantly from the 0.05 m to the 0.03 m grid cases.The errors for the 0.02 m and 0.03 m grid cases were similar to each other, leading us to choose the 0.02 m grid size for this study to achieve efficient simulation while reducing computational intensity.Three wave flumes with different beach slopes were simulated to investigate the effect of breaker types on the behaviors of microplastics in nearshore areas.The Iribarren parameter, a dimensionless parameter that categorizes wave breaker types, was used to categorize both beach slope and wave steepness [40].It is calculated as the ratio of the beach slope ( ) to the square root of the deep-water wave steepness ( / ) given by = tan / / (1) The 2D mesh in our study was characterized by a maximum grid size of 0.02 m, resulting in node (grid) counts of 75,189 (24,621), 75,564 (24,744), and 76,189 (24,947) for the three numerical flumes with different slopes, respectively.The two-dimensional (2D) computational mesh was developed using a rectangular quadrilateral grid with a quadratic element order.Our focus is on observing different behaviors for various particle properties, an approach that is particularly suitable for case studies involving multiple scenarios rather than aiming for high fidelity relative to true values.In this context, 2D simulations offer a streamlined and efficient means for analyzing the impact of varying particle properties on behavior across diverse scenarios while avoiding the computational complexity and resource intensity often associated with three-dimensional (3D) modeling.Previous studies have shown that 2D numerical wave flumes effectively predict the dynamic evolution of water waves [36][37][38][39].
Simulations were conducted with a numerical time step of 0.002 s, totaling 60,000 time steps over the 120 s simulation period.Each time step involved 60 iterations, and the computation time for each simulation was approximately 12 h.A grid sensitivity test, conducted using grid sizes of 0.01, 0.02, 0.03, and 0.05 m with 99,060, 75,564, 11,745, and 4410 nodes, respectively, showed that the Root Mean Square Error (RMSE) between the mean velocity fields of the 0.01 m case and the others decreased significantly from the 0.05 m to the 0.03 m grid cases.The errors for the 0.02 m and 0.03 m grid cases were similar to each other, leading us to choose the 0.02 m grid size for this study to achieve efficient simulation while reducing computational intensity.
Three wave flumes with different beach slopes were simulated to investigate the effect of breaker types on the behaviors of microplastics in nearshore areas.The Iribarren parameter, a dimensionless parameter that categorizes wave breaker types, was used to categorize both beach slope and wave steepness [40].It is calculated as the ratio of the beach slope (β) to the square root of the deep-water wave steepness (H 0 /L 0 ) given by ξ = tan β/ H 0 /L 0 (1) where β, H 0 , and L 0 represent the slope angle, wave height, and wavelength of deep water waves, respectively.This parameter leads to the classification of waves into three distinct categories: spilling, plunging, and surging.Spilling breakers, characterized by a gentle, extensive breaking process, occur when the Iribarren parameter is less than 0.4, typically on flat slopes [41].Plunging breakers, associated with a more abrupt breaking over a shorter distance, are identified when the Iribarren parameter lies between 0.4 and 2.0, often seen on moderately steep slopes [42].Surging breakers, which occur when the Iribarren parameter exceeds 2.0, are typically found on steep slopes, with waves breaking near or at the shoreline, resulting in a swift, rushing bore up the beach [43][44][45].To simulate these various wave breaker types, the numerical experiments were performed using three distinct slope angles (β 1 = 2.5 • , β 2 = 10 • , and β 3 = 35 • ) representing spilling, plunging, and surging breakers, respectively, as classified by the Iribarren parameter.

Simulation Settings
A double-precision, pressure-based solver was utilized to simulate the system's unsteady behavior, employing a bounded second-order implicit scheme for the transient formulation.The coupling of pressure and velocity was achieved through the PISO algorithm in combination with the PRESTO scheme, enabling accurate capture of the dynamics of surface waves, which necessitates efficient handling of rapidly changing pressure and velocity fields.For resolving various quantities, the second-order upwind scheme was applied to momentum, turbulent kinetic energy, specific dissipation rate, energy, and the level-set function.To ensure convergence, integrated quantities, such as bulk flow velocity and turbulence, along with scaled residuals of continuity, x and y velocities, energy, k-ω, and bulk water volume fraction, were monitored.Numerical convergence was determined when these monitored parameters stabilized and the residuals of all quantities became less than 0.001.
The VOF method, a surface-tracking technique that relies on the fluid volume fraction in each computational cell and incorporates open channel boundary conditions, was used to simulate water-air interactions in wave flumes [35] and two-phase flows [46].The densities of water and air were set at 998.2 kg/m 3 and 1.225 kg/m 3 , respectively, in this work.It solves a single set of momentum equations and tracks the volume fraction of each fluid through the computational domain.The fluid volume fraction, represented by α, varies from 0 to 1, where α = 0 denotes the absence of the primary fluid and α = 1 represents the cell completely filled with the primary fluid.The conservation of mass in the VOF framework is governed by the continuity equation as follows: where ρ is the fluid density, u is the velocity vector, and ϕ is the volume fraction of the fluid phase.This equation ensures that the mass is conserved in the fluid flow.The conservation of momentum is described by the Navier-Stokes equations, given by where P, τ, g, F σ , and ⊗ denote the pressure, stress tensor, gravitational force, surface tension, and outer product [30,47].These equations, when solved numerically, enable the detailed modeling of wave dynamics and interactions within a wave flume.The Discrete Phase Model (DPM) was employed in our study to simulate the transport of microplastic particles within a continuous fluid phase.The unsteady Reynolds-averaged Navier-Stokes (uRANS) equations and VOF model are solved based on the finite volume method.
The SST k-ω model [48][49][50] is particularly effective for accurately simulating flow characteristics near walls or boundaries, a critical aspect when studying nearshore environments and marine structures [37,51,52].This turbulence closer model is designed to effectively capture the turbulence behavior in both the near-wall region and the far-field.This hybrid model transitions from a standard k-ω model in the near-wall region to a k-ϵ model in the free stream.The model involves two transport equations.The equation for the turbulent kinetic energy (k) is provided by where P k denotes the production of turbulent kinetic energy, β * is a model constant, µ and µ t are the molecular and turbulent viscosities, respectively, and σ k is the turbulent Prandtl number for k.The specific dissipation rate (ω) is provided by where α, β, σ ω , and σ ω2 can be considered as model constants and F 1 is a blending function.
We simulated the propagation of waves within a transitional gravity wave field, employing Stokes' fourth-order wave theory.The wave height (H) and wavelength (L) were set at 0.25 m and 2.5 m, respectively (Table 1), and the corresponding wave period was 1.45 s.These conditions were imposed at the inlet boundary at x = 0.The resulting Ursell number, HL 2 /h 3 , a dimensionless parameter quantifying the ratio of nonlinearity to shallow water effects [53], was determined to be 3. Hence, the nonlinear features influence the dispersion of buoyant near-surface microplastics.The instantaneous waveforms for different breaker types for each beach slope case are demonstrated in Figure 1.

Microplastics in a Numerical Wave Tank
In the numerical wave tank, microplastic particles were conceptualized as spheres and introduced in a range of sizes (A = 0.2, 1.0, and 5.0 mm) and densities (ρ = 900, 998, and 1100 kg/m 3 ), as specified in Table 1.The particles with these three densities are referred to as buoyant, neutral, and heavy particles, respectively.These particles were deployed near the surface (y = 0.38 m) at four horizontal locations: offshore (red circles), shoaling (blue circles), breaking (green circles), and surfing (black circles), as marked in Figure 1.For each of these positions, a group of 50 particles was released.The 'offshore' location represents the area before the wave starts to interact with the slope.The 'shoaling' position is the point of initial wave-slope interaction.The 'breaking' position represents the zone between the 'shoaling' and 'surfing' areas, while the 'surfing' location is set near the shoreline, where the wave ultimately breaks (Figure 1).The simulations of particle movements began 10 s after the initiation of wave propagation to allow for the establishment of a stable wave field within the computational domain.

Eulerian Mean Characteristics
The Eulerian characteristics were derived from averaged velocity fields throughout the simulations.For the spilling breaker (Figure 2a,d), wave breakage predominantly occurred approximately one wavelength away from the shoreline, specifically near x = 8 m.This location was aligned with both the maximum shoreward horizontal velocity and the maximum wave height (Figure 2a).Across the expansive slope, significant vertical velocity was largely absent (Figure 2d), except for a single column of mean positive vertical velocity situated near x = 2 m outside the sloping region (Figure 2d).For the plunging breaker (Figure 2b,e), wave breaking primarily took place near the shoreline, accompanied by columns of positive vertical velocity outside the sloping region.For the surging breaker (Figure 2c,f), clear wave breaking was absent; instead, the entire wave front exhibited a surging motion along the shore.A pattern of undulations suggested the potential formation of standing waves due to wave reflection from steep beaches.At a distance greater than one wavelength from the start of the slope, a wider column of positive mean vertical velocity was observed (Figure 2f).The mean current revealed that a vertical current connecting the surface layer with the bottom layer was not found over the slope.
J. Mar.Sci.Eng.2024, 12, x FOR PEER REVIEW 6 of 17 breaker (Figure 2b,e), wave breaking primarily took place near the shoreline, accompanied by columns of positive vertical velocity outside the sloping region.For the surging breaker (Figure 2c,f), clear wave breaking was absent; instead, the entire wave front exhibited a surging motion along the shore.A pa ern of undulations suggested the potential formation of standing waves due to wave reflection from steep beaches.At a distance greater than one wavelength from the start of the slope, a wider column of positive mean vertical velocity was observed (Figure 2f).The mean current revealed that a vertical current connecting the surface layer with the bo om layer was not found over the slope.Figure 3 illustrates streamlines calculated from the mean velocity fields.In the spilling breaker, several thin circulations were observed over the sloping region.For the plunging breaker, a single circulation was observed within the slope region, accompanied by additional circulations evenly spread off the slope.In contrast, the surging breaker case lacked circulation over the narrow slope.However, the broadest circulation was detected in the sloping region, with its horizontal reach approximating 1.5 times the wavelength.The vertical columns of positive vertical velocity are aligned with the left side of each circulation zone.The center of circulation was skewed toward the surface, and the circulation zone was elongated in the horizontal direction, resulting in a predominantly negative domain of both horizontal and vertical velocities within the water column (Figure 2).  Figure 3 illustrates streamlines calculated from the mean velocity fields.In the spilling breaker, several thin circulations were observed over the sloping region.For the plunging breaker, a single circulation was observed within the slope region, accompanied by additional circulations evenly spread off the slope.In contrast, the surging breaker case lacked circulation over the narrow slope.However, the broadest circulation was detected in the sloping region, with its horizontal reach approximating 1.5 times the wavelength.The vertical columns of positive vertical velocity are aligned with the left side of each circulation zone.The center of circulation was skewed toward the surface, and the circulation zone was elongated in the horizontal direction, resulting in a predominantly negative domain of both horizontal and vertical velocities within the water column (Figure 2).2b,e), wave breaking primarily took place near the shoreline, accompanied by columns of positive vertical velocity outside the sloping region.For the surging breaker (Figure 2c,f), clear wave breaking was absent; instead, the entire wave front exhibited a surging motion along the shore.A pa ern of undulations suggested the potential formation of standing waves due to wave reflection from steep beaches.At a distance greater than one wavelength from the start of the slope, a wider column of positive mean vertical velocity was observed (Figure 2f).The mean current revealed that a vertical current connecting the surface layer with the bo om layer was not found over the slope.Figure 3 illustrates streamlines calculated from the mean velocity fields.In the spilling breaker, several thin circulations were observed over the sloping region.For the plunging breaker, a single circulation was observed within the slope region, accompanied by additional circulations evenly spread off the slope.In contrast, the surging breaker case lacked circulation over the narrow slope.However, the broadest circulation was detected in the sloping region, with its horizontal reach approximating 1.5 times the wavelength.The vertical columns of positive vertical velocity are aligned with the left side of each circulation zone.The center of circulation was skewed toward the surface, and the circulation zone was elongated in the horizontal direction, resulting in a predominantly negative domain of both horizontal and vertical velocities within the water column (Figure 2).

Lagrangian Mean Characteristics
The average horizontal paths of Lagrangian particles were demonstrated for the spilling breaker (Figure 4), plunging breaker (Figure 5), and surging breaker (Figure 6).Each figure represents the average trajectory of particles originating from four different locations, as shown in Figure 1.For the spilling breaker (Figure 4), buoyant particles were mainly collected at the shoreline, as indicated by the dashed line in Figure 4.Both neutral and heavy particles moved offshore, with heavier particles exhibiting higher speeds, although the differences in retreat speeds among particles of various sizes were not pronounced.Notably, among the buoyant particles, those with a diameter of 1 mm displayed the strongest inclination for accumulation at the shoreline or strong advection toward the shoreline.For instance, the mean forward speeds of buoyant particles released at the 'shoaling' location were 0.12 m/s for those with a diameter of 0.2 mm, 0.16 m/s for those with a diameter of 1 mm, and 0.08 m/s for those with a diameter of 5 mm.

Lagrangian Mean Characteristics
The average horizontal paths of Lagrangian particles were demonstrated for the spilling breaker (Figure 4), plunging breaker (Figure 5), and surging breaker (Figure 6).Each figure represents the average trajectory of particles originating from four different locations, as shown in Figure 1.For the spilling breaker (Figure 4), buoyant particles were mainly collected at the shoreline, as indicated by the dashed line in Figure 4.Both neutral and heavy particles moved offshore, with heavier particles exhibiting higher speeds, although the differences in retreat speeds among particles of various sizes were not pronounced.Notably, among the buoyant particles, those with a diameter of 1 mm displayed the strongest inclination for accumulation at the shoreline or strong advection toward the shoreline.For instance, the mean forward speeds of buoyant particles released at the 'shoaling' location were 0.12 m/s for those with a diameter of 0.2 mm, 0.16 m/s for those with a diameter of 1 mm, and 0.08 m/s for those with a diameter of 5 mm.For the plunging breaker (Figure 5), the buoyant particles also accumulated at the shoreline, but neutral and heavy particles initially moved forward before switching to offshore motion, particularly those particles released from 'shoaling' and 'offshore' locations.The surging breaker (Figure 6) exhibited pa erns similar to those of the other breaker types.However, a distinct oscillatory surging motion was observed, particularly for particles released from the 'surfing' location (Figure 6a).This is likely a ributed to the limited wave breaking, leading to a surging motion with minimal energy dissipation.Similar to the other breakers, neutral particles released from the 'shoaling' and 'offshore' regions take a considerable time to leave the nearshore area.wave breaking, leading to a surging motion with minimal energy dissipation.Similar to the other breakers, neutral particles released from the 'shoaling' and 'offshore' regions take a considerable time to leave the nearshore area.For the plunging breaker (Figure 5), the buoyant particles also accumulated at the shoreline, but neutral and heavy particles initially moved forward before switching to offshore motion, particularly those particles released from 'shoaling' and 'offshore' locations.The surging breaker (Figure 6) exhibited patterns similar to those of the other breaker types.However, a distinct oscillatory surging motion was observed, particularly for particles released from the 'surfing' location (Figure 6a).This is likely attributed to the limited wave breaking, leading to a surging motion with minimal energy dissipation.Similar to the other breakers, neutral particles released from the 'shoaling' and 'offshore' regions take a considerable time to leave the nearshore area.
In general, the mean trajectories indicated that buoyant particles moved forward, while neutral and heavy particles tended to move backward.The average movement of neutral particles was slower compared to heavy particles over a smooth bottom, likely due to the time taken for neutral particles to reach the bottom.Neutral particles released near the surface initially moved forward due to the forward surface current concentrated near the water surface, known as Stokes drift [21,22].Once the neutrally buoyant particles were entrained into deeper waters and distributed throughout the water column via waveinduced mixing, their average horizontal location shifted offshore.This shift occurred because the majority of the particles moved to the deeper layer, where negative horizontal velocity dominates, known as the undertow (Figure 2a-c).
It was also observed that particle advection was selectively influenced by particle density and size.For particles of the same size, mean advection varied dramatically according to particle density; for instance, light particles moved forward while neutral and heavy particles moved backward.For particles of the same density, mean advection was influenced by particle size; small particles (A = 0.2 mm) moved faster, whereas larger particles (A = 5 mm) moved slower.However, the relationship between the strength of advection and particle size, with a fixed particle density, was not linear; medium-sized particles moved the fastest.This suggests that the drag force, influenced by viscosity and the pressure propelling the particle, competes with the inertial force, which can either hinder or sustain the particle's movement.This competition can lead to either stronger or weaker advection at a certain wave period, depending on the particle size.

Reynolds Stresses over the Slopes
We computed the Reynolds stresses (S xx (= u ′ u ′ ), S yy (= v ′ v ′ ), and S xy (= u ′ v ′ )) using the mean current to examine the turbulent flow structure induced by wave motions for each breaker type.These stresses are derived by time-averaging the products of velocity fluctuations (u ′ and v ′ ), which are calculated using Reynolds decomposition (u ′ = u − u), where u represents the raw velocity and the overbar denotes the time-averaged value.For the spilling breaker, strong S xx was observed just below the still water level (SWL, indicated as y ′ = 1) and gradually decreased with depth (Figure 7a).The strong S yy (Figure 7c) signifies robust vertical mixing near the SWL, extending from the wave crest to the trough, except in the bottom layer (Figure 8c).This vertical mixing, however, diminishes near the bottom and past the breaking point.However, beyond the point of maximum wave breaking near the shoreline (x ′ > 0.75), all stresses significantly diminished.Strong Reynolds stresses in the horizontal direction due to vertical velocity fluctuations (S xy ) were evident above the SWL but weakened below it (Figure 7b).In contrast, the surging breaker exhibited pronouncedly high values for all Reynolds stresses near the shoreline (Figure 7g-i).The plunging breaker showed an intermediate distribution between the spilling and surging breakers (Figure 7d-f).The pronounced S xx (Figure 7a) implies that the particles positioned just below SWL, where strong mean shear can be found, undergo significant horizontal dispersion due to wave-induced turbulence.Conversely, particles positioned above the SWL tend to experience strong forward advection, as evidenced by the mean surface current.Smaller particles, which are sufficiently small to be distributed vertically, are more susceptible to layers of intense horizontal dispersion compared to larger, highly buoyant particles, like the 5 mm particles, that are positioned near the surface.Thus, depending on the buoyancy differences due to particle size, the dispersion and advection characteristics of microplastics can vary.

Buoyant Particle Distribution across the Slopes
In order to examine the features of buoyant particle distribution during their phase of accumulation towards the beach, we computed the time series of STD, skewness, and kurtosis of the horizontal positions of particles originating from the 'Shoaling' location (Figure 8).For the spilling breaker (Figure 8a-c), we identified a large deviation in the horizontal particle distribution depending on the particle size while approaching the shoreline.The 1 mm particles with the fastest forward advection (Figure 4) showed the narrowest distribution (small STD) (Figure 8).The 0.2 mm particles showed the broadest distribution (large STD), which aligns with the observations that small particles of order 0.1 mm undergo extensive mixing in all directions, regardless of their density [32].The dispersion was maximized around three-quarters of the distance from the starting point (x ′ = 0.75), analogous to the strong sediment dispersion under conditions of higher breaking wave height [53,54].As particles move toward the coastline, their spatial distribution gradually becomes strongly skewed toward the coastline, which results in high kurtosis and negative skewness (Figure 7b,c), appearing to converge into a narrow region near the coastline.The distribution of buoyant microplastics tended to approximate a Gaussian distribution (kurtosis = 3) in areas with high horizontal stress due to wave motions (Figure 7a-c).However, near the shoreline, where the stress diminished, a noticeable skew towards the coastline was observed.
For the plunging breaker (Figure 8d-f), all particle sizes exhibited similar statistical distributions, though the 5 mm particles showed the slowest advection, and the 0.2 mm particles had the broadest distribution.All buoyant particles were observed moving toward the beach, but their dispersion was weaker than that in the spilling breakers.Interestingly, the particles tended to distribute without a distinct inclination (skewness ~0), with the distribution shape being close to a Gaussian distribution.
For the surging breakers (Figure 8g-i), the weak wave-breaking nature may be responsible for the reduced particle dispersion (lowest STD).The equilibrium center point of particle distribution, in proximity to the coastline, was observed in the order of 0.2 mm, 1 mm, and 5 mm particles.Once near the shoreline, buoyant particles were subjected to strong horizontal stress (Figure 7g-i), leading to their selective positioning based on the relative impact of particle inertia.The 5 mm particles, influenced more by inertia, and the 0.2 mm particles, affected more by viscosity, demonstrated a unique balance between these forces, reflecting distinct movement patterns at equilibrium.

Lagrangian vs. Eulerian Velocities
The inertial effects, depending on the particle size, can be visualized by comparing the Eulerian velocity (u E ) and the Lagrangian velocity (u L ) (Figure 9).The mean squared differences in the horizontal velocity ( (u L − u E ) 2 ) were found to be 0.128, 0.199, and 0.214 m/s for 0.2, 1, and 5 mm particles, respectively.Similarly, the mean squared differences in the vertical velocity ( (v L − v E ) 2 ) were 0.072, 0.186, and 0.311 m/s for 0.2, 1, and 5 mm particles, respectively, as illustrated in Figure 10.The smaller gap between u L and u E for the 0.2 mm particles suggests that small particles of order 0.1 mm tend to follow the background flows, with the observed deviation increasing with particle size (Figure 9). Figure 10 illustrates the differences between the Lagrangian kinetic energy spectrum (LKE) and the Eulerian kinetic energy spectrum (EKE) for both vertical (represented by red lines) and horizontal (black lines) velocities.For small particles (0.2 mm), the spectra show minimal differences, highlighting their tendency to passively follow the background flows.As the particle size increases, a more pronounced deviation between the Lagrangian and Eulerian velocity spectra becomes apparent.In the case of the 1 mm particles, LKE in the horizontal direction (LKE-h) is slightly smaller than EKE in the horizontal direction (EKE-h), while LKE in the vertical direction (LKE-v) is somewhat larger than EKE in the vertical direction (EKE-v).Notably, these differences in the spectrum predominantly occur near the wave period.This suggests that the wave motion tends to suppress horizontal motion while enhancing the vertical motion of these particles (A = 1 mm) moving with the wave period.− indicates the difference between Lagrangian and Eulerian horizontal velocity spectra, and − indicates the difference between Lagrangian and Eulerian vertical velocity spectra.The Lagrangian velocity was measured by the particle released from the 'offshore' location, and the Eulerian velocity was measured at the location of the Lagrangian particle.Figure 10 illustrates the differences between the Lagrangian kinetic energy spectrum (LKE) and the Eulerian kinetic energy spectrum (EKE) for both vertical (represented by red lines) and horizontal (black lines) velocities.For small particles (0.2 mm), the spectra show minimal differences, highlighting their tendency to passively follow the background flows.As the particle size increases, a more pronounced deviation between the Lagrangian and Eulerian velocity spectra becomes apparent.In the case of the 1 mm particles, LKE in the horizontal direction (LKE-h) is slightly smaller than EKE in the horizontal direction (EKE-h), while LKE in the vertical direction (LKE-v) is somewhat larger than EKE in the vertical direction (EKE-v).Notably, these differences in the spectrum predominantly occur near the wave period.This suggests that the wave motion tends to suppress horizontal motion while enhancing the vertical motion of these particles (A = 1 mm) moving with the wave period.Figure 10 illustrates the differences between the Lagrangian kinetic energy spectrum (LKE) and the Eulerian kinetic energy spectrum (EKE) for both vertical (represented by red lines) and horizontal (black lines) velocities.For small particles (0.2 mm), the spectra show minimal differences, highlighting their tendency to passively follow the background flows.As the particle size increases, a more pronounced deviation between the Lagrangian and Eulerian velocity spectra becomes apparent.In the case of the 1 mm particles, LKE in the horizontal direction (LKE-h) is slightly smaller than EKE in the horizontal direction (EKE-h), while LKE in the vertical direction (LKE-v) is somewhat larger than EKE in the vertical direction (EKE-v).Notably, these differences in the spectrum predominantly occur near the wave period.This suggests that the wave motion tends to suppress horizontal motion while enhancing the vertical motion of these particles (A = 1 mm) moving with the wave period.
For the larger 5 mm particles, the difference becomes more distinct but occurs at a different period (Figure 10c).The LKE-h is consistently lower than the EKE-h at the wave period, with a more pronounced difference compared to the smaller particles.Meanwhile, LKE-v significantly exceeds EKE-v, with this marked discrepancy occurring over a longer period of 2.5 s.Consequently, the horizontal motion of these larger particles is noticeably restrained during the wave period, whereas their vertical motion is substantially amplified at periods longer than the wave period, emphasizing the influence of particle size on their dynamic response to wave-induced motions.

Inertial Effects
The discrepancy between Eulerian and Lagrangian velocities can be intricately linked to the Stokes number (St = t p /t f ) [55].This dimensionless parameter is defined as the ratio of the particle relaxation time (t p = ρ p A/18µ, where ρ p , A, and µ are particle density, particle diameter, and fluid dynamic viscosity, respectively) to a characteristic fluid flow timescale (t f = 1/T, where T is a wave period).The particle relaxation time depends on the particle's size and density, while the fluid flow timescale can be considered a wave period.The Stokes number juxtaposes particle inertia against the fluid's viscous forces, helping to understand how particles of varying sizes interact with the fluid flow.The Stokes numbers for buoyant particles of 0.2, 1, and 5 mm were calculated to be 0.001, 0.035, and 0.862 (close to unity), respectively.This suggests that the impact of inertia, and consequently the deviation between Eulerian and Lagrangian velocities, likely increases with particle size.
The inertial effects significantly influence the dispersion and advection of buoyant microplastics.The 5 mm particles, with a higher Stokes number and thus more affected by inertia, exhibited the lowest dispersion and advection, as shown in Figure 11a,b.This is because particles with higher inertia resist changes in their motion due to turbulent fluctuations, making them less responsive to smaller, rapid changes in the fluid's velocity, resulting in less dispersion.Additionally, their higher inertia means these particles maintain their momentum for longer periods, making them less susceptible to being carried by the fluid during advection processes, resulting in lower advection rates as the particles are not easily swept along by wave-induced currents.For the larger 5 mm particles, the difference becomes more distinct but occurs at a different period (Figure 10c).The LKE-h is consistently lower than the EKE-h at the wave period, with a more pronounced difference compared to the smaller particles.Meanwhile, LKE-v significantly exceeds EKE-v, with this marked discrepancy occurring over a longer period of 2.5 s.Consequently, the horizontal motion of these larger particles is noticeably restrained during the wave period, whereas their vertical motion is substantially amplified at periods longer than the wave period, emphasizing the influence of particle size on their dynamic response to wave-induced motions.

Inertial Effects
The discrepancy between Eulerian and Lagrangian velocities can be intricately linked to the Stokes number ( = / ) [55].This dimensionless parameter is defined as the ratio of the particle relaxation time ( = /18 , where , , and are particle density, particle diameter, and fluid dynamic viscosity, respectively) to a characteristic fluid flow timescale ( = 1/ , where is a wave period).The particle relaxation time depends on the particle's size and density, while the fluid flow timescale can be considered a wave period.The Stokes number juxtaposes particle inertia against the fluid's viscous forces, helping to understand how particles of varying sizes interact with the fluid flow.The Stokes numbers for buoyant particles of 0.2, 1, and 5 mm were calculated to be 0.001, 0.035, and 0.862 (close to unity), respectively.This suggests that the impact of inertia, and consequently the deviation between Eulerian and Lagrangian velocities, likely increases with particle size.
The inertial effects significantly influence the dispersion and advection of buoyant microplastics.The 5 mm particles, with a higher Stokes number and thus more affected by inertia, exhibited the lowest dispersion and advection, as shown in Figure 11a,b.This is because particles with higher inertia resist changes in their motion due to turbulent fluctuations, making them less responsive to smaller, rapid changes in the fluid's velocity, resulting in less dispersion.Additionally, their higher inertia means these particles maintain their momentum for longer periods, making them less susceptible to being carried by the fluid during advection processes, resulting in lower advection rates as the particles are not easily swept along by wave-induced currents.Conversely, the 0.2 mm particles, with a lower Stokes number and more influenced by viscosity, exhibited higher dispersion (as seen in Figure 11b).The dominant viscous forces allowed for greater interaction with wave-induced turbulent motions, leading to a Conversely, the 0.2 mm particles, with a lower Stokes number and more influenced by viscosity, exhibited higher dispersion (as seen in Figure 11b).The dominant viscous forces allowed for greater interaction with wave-induced turbulent motions, leading to a broader spread across the water column.The smaller size and lower inertia of these particles make them more responsive to changes in the fluid's velocity, contributing to their extensive dispersal in the nearshore environment.
Interestingly, the 1 mm particles demonstrated the strongest forward advection, as evident in Figures 4g-i and 11a.This pronounced advection is primarily attributed to their optimal balance between inertia and viscosity.The intermediate size of these particles enables them to be significantly influenced by advection forces, such as currents and wave motion, while still being small enough to be impacted by the water's viscous forces.This unique balance facilitates an effective response to the forward momentum of the waves, leading to more pronounced advection.Furthermore, unlike their smaller 0.2 mm counterparts, the 1 mm particles are less susceptible to small-scale turbulent fluctuations, enabling them to maintain a more direct and consistent forward motion, following the prevailing wave and current directions.

Conclusions
This study presents a comprehensive analysis of microplastic behavior in nearshore environments using a 2D CFD framework.We investigated the dynamics of different breaker types-spilling, plunging, and surging-and their influence on microplastic dispersion and accumulation.In spilling breakers, buoyant particles mainly accumulate at the shoreline, while neutral and heavy particles move offshore, influenced by a strong undertow affecting the water column.Plunging breakers show similar shoreline accumulation of buoyant particles, but neutral and heavy particles initially move forward, then offshore, with a more Gaussian distribution and less dispersion compared to spilling breakers.Surging breakers, on the other hand, are marked by a unique oscillatory surging motion, exhibiting the weakest dispersion due to limited wave breaking and minimal energy dissipation.
The Eulerian analysis revealed that the majority of the water column is an undertow zone, except for the near-surface layer, which possesses a strong horizontal forward velocity, impacting the horizontal mean location of particles.The mean velocity field uncovers circulation zones that act as conduits between near-surface and near-bottom transport.From a Lagrangian perspective, particle advection was observed to be selectively influenced by particle density and size.In terms of advection, for a fixed particle size, light particles moved forward while neutral and heavy particles moved backward.Generally, small particles exhibited rapid advection, whereas larger particles (A = 5 mm) showed slower advection.However, the relationship between the strength of advection and particle size, with a fixed particle density, was not linear: the 1 mm buoyant particles demonstrated the strongest forward advection, likely due to an optimal balance between inertia and viscosity.These particles were large enough to be significantly affected by advection forces, such as currents and wave motion, yet small enough to remain influenced by the water's viscous forces.This balance enabled them to effectively respond to the forward momentum of waves, resulting in pronounced advection.Conversely, the 5 mm particles, being more influenced by inertia, tended to settle more rapidly and closer to the point of wave breaking.In contrast, the 0.2 mm particles, more affected by viscosity, displayed higher dispersion, reflecting their greater interaction with wave-induced turbulent motions and responsiveness to fluid velocity changes.
This study sheds light on the complex interplay of physical forces that determine microplastic distribution in nearshore environments.Understanding these dynamics is crucial for developing effective strategies against microplastic pollution in marine ecosystems.Our research highlights the significant impact of inertia on microplastic distribution and accumulation patterns in nearshore areas.This emphasizes the importance of incorporating inertial effects into models to more accurately predict microplastic behavior, especially in zones where varying breaker types influence their movement and final deposition.Future research could extend this study to a 3D wave tank, which would provide a more accurate analysis of particle motion, including the effects of fluctuations in the third dimension.Although our focus was on spherical microplastics, the study could be expanded to include various shapes, such as fibers and plates, to underscore the impact of shape on drag force.This expansion would necessitate an accurate empirical formulation of the drag coefficient corresponding to each shape.

Figure 2 .
Figure 2. Time-average horizontal velocity (a-c) and time-averaged vertical velocity (d-f).Colored contour indicates positive velocity, and the regions of negative velocity were replaced as nulls.Logscale was applied only for v-velocity.

Figure 2 .
Figure 2. Time-average horizontal velocity (a-c) and time-averaged vertical velocity (d-f).Colored contour indicates positive velocity, and the regions of negative velocity were replaced as nulls.Log-scale was applied only for v-velocity.

Figure 2 .
Figure 2. Time-average horizontal velocity (a-c) and time-averaged vertical velocity (d-f).Colored contour indicates positive velocity, and the regions of negative velocity were replaced as nulls.Logscale was applied only for v-velocity.

Figure 4 .
Figure 4. Mean horizontal particle trajectory for a spilling breaker originating from four different release locations: 'surfing' (a-c), 'breaking'(d-f), 'shoaling' (g-i), and 'offshore' (j-l).The x-axis represents the simulation time, and the y-axis represents the horizontal coordinate, as illustrated in Figure1.The four initial release locations are marked with black dots on the y-axis in the left panel.The dashed line indicates the shoreline, and the dashed-dot line indicates the position where the slope starts.So, the region between the dashed and dash-dot lines denotes the sloping region.'A' indicates the particle size, and , , and represent the particle densities of 900, 998.2, and 1100 kg/m 3 , respectively.

Figure 4 .
Figure 4. Mean horizontal particle trajectory for a spilling breaker originating from four different release locations: 'surfing' (a-c), 'breaking' (d-f), 'shoaling' (g-i), and 'offshore' (j-l).The x-axis represents the simulation time, and the y-axis represents the horizontal coordinate, as illustrated in Figure 1.The four initial release locations are marked with black dots on the y-axis in the left panel.The dashed line indicates the shoreline, and the dashed-dot line indicates the position where the slope starts.So, the region between the dashed and dash-dot lines denotes the sloping region.'A' indicates the particle size, and ρ 1 , ρ 2 , and ρ 3 represent the particle densities of 900, 998.2, and 1100 kg/m 3 , respectively.

Figure 5 .
Figure 5. Mean horizontal particle trajectory for a plunging breaker, originating from four different release locations: surfing (a-c), breaking (d-f), shoaling (g-i), and offshore (j-l).The rest of the details are the same as described in Figure 4.

Figure 6 .
Figure 6.Mean horizontal particle trajectory for a surging breaker, originating from four different release locations: surfing (a-c), breaking (d-f), shoaling (g-i), and offshore (j-l).The rest of the details are the same as described in Figure4.

Figure 5 .
Figure 5. Mean horizontal particle trajectory for a plunging breaker, originating from four different release locations: surfing (a-c), breaking (d-f), shoaling (g-i), and offshore (j-l).The rest of the details are the same as described in Figure 4.

Figure 5 .
Figure 5. Mean horizontal particle trajectory for a plunging breaker, originating from four different release locations: surfing (a-c), breaking (d-f), shoaling (g-i), and offshore (j-l).The rest of the details are the same as described in Figure 4.

Figure 6 .
Figure 6.Mean horizontal particle trajectory for a surging breaker, originating from four different release locations: surfing (a-c), breaking (d-f), shoaling (g-i), and offshore (j-l).The rest of the details are the same as described in Figure4.

Figure 6 .
Figure 6.Mean horizontal particle trajectory for a surging breaker, originating from four different release locations: surfing (a-c), breaking (d-f), shoaling (g-i), and offshore (j-l).The rest of the details are the same as described in Figure 4.

Figure 7 .
Figure 7. Momentum fluxes Sxx, Sxy, and Syy for spilling (a-c), plunging (d-f), and surging (g-i) breakers.x' and y' are dimensionless distances.x' = 0 indicates the 'shoaling' location, and x' = 1 indicates the intersection between slope and mean water level.y' = 0 and y' = 1 indicate the bo om and still water level (h = 0.4 m).Negative values were replaced with nulls.All contours were logscaled.

Figure 8 .
Figure 8.Time series of standard deviation (STD), skewness, and kurtosis of buoyant particles ( = 900 kg/m 3 ) distribution for spilling (a-c), plunging (d-f), and surging (g-i) breakers during the simulation time.x-axis indicates the nondimensional mean location of particles released from the 'shoaling' location.x' = 0 indicates the 'shoaling' location, and x' = 1 indicates the intersection between the slope and still water level.The dashed lines for skewness and kurtosis indicate y = 0 and y = 3, respectively.STD and kurtosis are scaled as log scales.

Figure 7 . 17 Figure 7 .
Figure 7. Momentum fluxes S xx , S xy , and S yy for spilling (a-c), plunging (d-f), and surging (g-i) breakers.x ′ and y ′ are dimensionless distances.x ′ = 0 indicates the 'shoaling' location, and x ′ = 1 indicates the intersection between slope and mean water level.y ′ = 0 and y ′ = 1 indicate the bottom and still water level (h = 0.4 m).Negative values were replaced with nulls.All contours were log-scaled.

Figure 8 .
Figure 8.Time series of standard deviation (STD), skewness, and kurtosis of buoyant particles ( = 900 kg/m 3 ) distribution for spilling (a-c), plunging (d-f), and surging (g-i) breakers during the simulation time.x-axis indicates the nondimensional mean location of particles released from the 'shoaling' location.x' = 0 indicates the 'shoaling' location, and x' = 1 indicates the intersection between the slope and still water level.The dashed lines for skewness and kurtosis indicate y = 0 and y = 3, respectively.STD and kurtosis are scaled as log scales.

Figure 8 .
Figure 8.Time series of standard deviation (STD), skewness, and kurtosis of buoyant particles (ρ = 900 kg/m 3 ) distribution for spilling (a-c), plunging (d-f), and surging (g-i) breakers during the simulation time.x-axis indicates the nondimensional mean location of particles released from the 'shoaling' location.x ′ = 0 indicates the 'shoaling' location, and x ′ = 1 indicates the intersection between the slope and still water level.The dashed lines for skewness and kurtosis indicate y = 0 and y = 3, respectively.STD and kurtosis are scaled as log scales.

Figure 9 .
Figure 9. Lagrangian velocity of a buoyant particle (represented by red lines) and Eulerian velocity at the particle's location (represented by black lines) for particles of sizes 0.2 mm (a,d), 1 mm (b,e), and 5 mm (c,f) in the spilling breaker.The upper plots show the horizontal velocity, while the lower plots depict the vertical velocity component.The particle was released from the 'offshore' location.

Figure 10 .
Figure 10.Differences in kinetic energy spectrum between Lagrangian velocity and Eulerian velocity for different particle sizes (A) of 0.2 mm (a), 1 mm (b), and 5 mm (c).−indicates the difference between Lagrangian and Eulerian horizontal velocity spectra, and − indicates the difference between Lagrangian and Eulerian vertical velocity spectra.The Lagrangian velocity was measured by the particle released from the 'offshore' location, and the Eulerian velocity was measured at the location of the Lagrangian particle.

Figure 9 .
Figure 9. Lagrangian velocity of a buoyant particle (represented by red lines) and Eulerian velocity at the particle's location (represented by black lines) for particles of sizes 0.2 mm (a,d), 1 mm (b,e), and 5 mm (c,f) in the spilling breaker.The upper plots show the horizontal velocity, while the lower plots depict the vertical velocity component.The particle was released from the 'offshore' location.

Figure 9 .
Figure 9. Lagrangian velocity of a buoyant particle (represented by red lines) and Eulerian velocity at the particle's location (represented by black lines) for particles of sizes 0.2 mm (a,d), 1 mm (b,e), and 5 mm (c,f) in the spilling breaker.The upper plots show the horizontal velocity, while the lower plots depict the vertical velocity component.The particle was released from the 'offshore' location.

Figure 10 .
Figure 10.Differences in kinetic energy spectrum between Lagrangian velocity and Eulerian velocity for different particle sizes (A) of 0.2 mm (a), 1 mm (b), and 5 mm (c).−indicates the difference between Lagrangian and Eulerian horizontal velocity spectra, and − indicates the difference between Lagrangian and Eulerian vertical velocity spectra.The Lagrangian velocity was measured by the particle released from the 'offshore' location, and the Eulerian velocity was measured at the location of the Lagrangian particle.

Figure 10 .
Figure 10.Differences in kinetic energy spectrum between Lagrangian velocity and Eulerian velocity for different particle sizes (A) of 0.2 mm (a), 1 mm (b), and 5 mm (c).u L − u E indicates the difference between Lagrangian and Eulerian horizontal velocity spectra, and v L − v E indicates the difference between Lagrangian and Eulerian vertical velocity spectra.The Lagrangian velocity was measured by the particle released from the 'offshore' location, and the Eulerian velocity was measured at the location of the Lagrangian particle.

Figure 11 .
Figure 11.Advection and dispersion of buoyant particles in the spilling breaker: (a) mean location of particles ( ), (b) STD of particle locations (STD( )), where is particle locations in the x-direction.

Figure 11 .
Figure 11.Advection and dispersion of buoyant particles in the spilling breaker: (a) mean location of particles (x p ), (b) STD of particle locations (STD(x p )), where x p is particle locations in the x-direction.
administration, Y.-G.P.; Funding acquisition, Y.-G.P. and I.-c.L.All authors have read and agreed to the published version of the manuscript.

Table 1 .
Wave conditions and microplastic particle properties.