Clogging and Unclogging of Fine Particles in Porous Media: Micromechanical Insights From an Analog Pore System

Pore clogging and unclogging in porous media are ubiquitous in subsurface hydrologic processes, which have been studied extensively at various scales ranging from a single pore to porous‐medium samples. However, it remains unclear how fluid flow, particle rearrangement, and the arching effect typical of cone‐shaped pore geometry interact and how they are captured by a pressure drop model at the macroscopic scale. Here, we investigate the pore‐scale feedback mechanisms between fluid flow and pore clogging and unclogging using a fully resolved fluid‐particle coupling approach (lattice Boltzmann method‐discrete element method). We first propose to use a truncated‐cone pore to represent realistic pore geometries revealed by X‐ray images of prepared sand packing. Then, our simulations indicate that the pore cone angle significantly influences the pressure drop associated with the clogging process by enhancing particle contacts due to arching. A modified Ergun equation is developed to consider this geometric effect. At the microscale, clogging can be explained by the interparticle force statistics; a few particles in an arch (or a dome) take the majority of hydrodynamic pressure. The maximum interparticle force is positively proportional to the particle Reynolds number and negatively associated with the tangent of the pore cone angle. Finally, a formula is established utilizing fluid characteristics and pore cone angle to compute the maximal interparticle force. Our findings, especially a modified pressure drop model that accounts for pore geometry resistance, provide guidance for applying pore‐scale models of clogging and unclogging to large‐scale subsurface fines transportation issues, including seepage‐induced landslides, stream bank failure, and groundwater recharge.

Pore clogging has been studied in gravity-driven and fluid-driven systems, respectively.The former system usually consists of a hopper with an opening (constriction) size of D c and particles of diameter D p that pass through the opening under gravity.Pore clogging in this situation is affected by numerous variables, such as the constriction-to-particle size ratio D c /D p (Gerber et al., 2018;Wyss et al., 2006), hopper angle δ (the acute angle between the axial direction and the inclined wall of the hopper) (Lai et al., 2001;López-Rodríguez et al., 2019), particle velocity v p (Gella et al., 2018), molecular force (Jäger et al., 2018;Liu et al., 2019), and particle concentration (Agbangla et al., 2014;Vani et al., 2022).The values of D c /D p and δ characterize the container shape, respectively, and determine the number of particles passing through the constriction and the wall resistance to particle migration.A larger δ enhances pore clogging (Lai et al., 2001;López-Rodríguez et al., 2019).When D c /D p is less than one, size exclusion almost certainly occurs for spherical particles and a circular constriction (Gerber et al., 2018), while when D c /D p is more than one, pore-clogging becomes a probabilistic occurrence with the probability P c decreasing exponentially as D c /D p increases, as stated in Equation 1 (Gella et al., 2018;Valdes & Santamarina, 2008): where v a is the absolute value of the arithmetic mean value of the particles' velocities, a represents the inertial effects, and b describes the impact of particle velocity on clogging.
Equation 1 also demonstrates that P c is inversely proportional to the particle velocity v p , which is strongly affected by D c /D p in a gravity-driven system (Beverloo et al., 1961) and by the fluid velocity in a fluid-driven system (Yin et al., 2021).Pore clogging occurs when particles interacting across the constriction form a two-dimensional arch or a three-dimensional dome.
Unless there is an external disturbance, such as vibration, it seems that pore-clogging is permanent in gravity-driven systems (Caitano et al., 2021;Valdes & Santamarina, 2008).In a fluid-driven system, in contrast, the fluid continues to flow through the interparticle pores and imposes a stronger drag on clogging particles, while the clogging restricts fluid flow, resulting in a pressure drop that provides pressure gradient force for pore unclogging (Bianchi et al., 2018;Jäger et al., 2017;Seybold et al., 2021).According to the simulation investigations conducted by Jäger et al. (2018), pores get unclogged when the rise in pressure drop exceeds the adhesion force between the clogging particles and the pore wall.An equation relating erosion pressure drop ∆p to the YIN ET AL.However, the model above only applies to cylindrical pore geometries but not more realistic pore structures typical of soil and rock mass (Guo et al., 2019;B. Nie et al., 2015;Winterberg & Tsotsas, 2000).Do complex-shaped pores and clogging structures together influence fluid movement and determine the ∆p required for unclogging?Moreover, in a general situation where the pore structure allows particles to form arches and promote jamming, other force terms (in addition to f b ) should also play a role.In particular, the contact force between particles is the primary force responsible for arching clogging.However, it is unclear how to consider this mechanical aspect to establish a more accurate relationship between ∆p and intergranular forces to extend the range of applicability of Equation 2.
To shed light on this issue, we focus on particle-particle interactions, particle-pore interactions, and pore-fluid interactions at the pore scale.Experimentally, it is still challenging to assess the interparticle stress in the three-dimensional structure of particle packings (Shen & Ni, 2017), even though photoelastic methods are often utilized to probe the interaction force between particles in granular mechanics studies (Behringer, 2015).Alternatively, particle-resolved computational methods for fluid-structure coupling enable us to visualize the three-dimensional interactions between particles and flow fields (Jing et al., 2016;Yang et al., 2021).In this regard, numerical methods coupling computational fluid dynamics, lattice Boltzmann method (LBM), or smoothed particle hydrodynamics with the discrete element method (DEM) have been increasingly used to provide fundamental insights into pore clogging (Sun et al., 2019;K. Zhou et al., 2018), landslides (Jing et al., 2018;Yang et al., 2019), internal soil erosion (Xiong et al., 2020), and debris flows (K.F. E. Cui et al., 2021;Ng et al., 2023).
To accurately capture the influence of fluid pressure on each particle during pore clogging and unclogging, it is necessary to fully resolve the flow field surrounding a particle using a finer grid for fluid computation at the pore scale.Therefore, we use the LBM for simulating fluid movement in microscopic pore channels in this paper.Besides, LBM is highly efficient in parallel computing and the calculation of complicated fluid-solid interaction boundaries, which can also guarantee hydrodynamic force accuracy (Neuville et al., 2013;L. Zhang & Wang, 2015).DEM is used to determine the migration of suspended spherical particles and their interactions.
The immersed movement boundary (IMB) is used to calculate the interactions between fluids and particles YIN ET AL.

10.1029/2023WR034628
4 of 29 because this approach assures perfect mass conservation and has a higher sub-grid scale precision (Peskin, 1977).
The LBM-DEM-IMB numerical experiments are conducted using the Palabos and LIGGGHTS open-source programs (Kloss et al., 2012;Latt et al., 2021;Seil & Pirker, 2017;Yang et al., 2019).In the following sections, the methodology, variables, particle and fluid parameters, and model validity verification are first introduced.We then present results, including (a) the movement characteristics of particles and the fluid during pore clogging and unclogging, (b) a modified Ergun model of fluid pressure drop for cone-shaped single-pore clogging, and (c) the characteristics and an empirical relationship for intergranular forces that are responsible for arching effects.Finally, we discuss the significance of the results in the context of engineering practice and geohazard mitigation before conclusions are drawn.

Pore Geometry Simplification
In previous porous media studies, three-dimensional pore channels are usually simplified as cylindrical pores, conical or wedge-shaped pores, parallel plate pores, or narrow bottleneck pores (B.Nie et al., 2015;Y. Zhang et al., 2016).Here, to construct a representative pore structure model that is relevant to typical natural sediments and geomaterials, the X-ray imaging results are used to examine the shape of pore channels in two samples: quartz sand and sandstone, based on which we propose to use a truncated cone-shaped pore to simplify and represent pore structures in these more realistic porous samples.Using X-ray imaging, the first sample we examine is a prepared quartz sand cylindrical sample, which consists of two sieving particle sizes: coarse particles 0.6-0.8mm with a content of 80% (standard quartz sand of China) and fine particles 0.063-0.075mm with a content of 20%.The sample is prepared based on observations that, for loose particle packing in geomaterials, the size ratio of coarse and fine particles is generally less than 10, and the fines content is usually less than 35% in experiments for fines migration in porous media (L.Zhang, Peng, Chang, & Xu, 2016).
As illustrated in Figure 2a, we focus on pores generated by coarse particles without fine particles, which form the channels for fine particle migration.As indicated in the cross-section of the computed tomography (CT) result, the majority of pores are gradually reduced in size from the opening to the pore constriction.The second sample we examine is the CT image of a Bentheim sandstone sample by micro-CT with a voxel resolution of 4.9 μm, which is published in previous open access work (Wetzel et al., 2020), as shown in Figure 2b.Although these researchers acquired CT images to study the effect of mineral growth on pore space and fluid transportation in porous media, their rock CT images indicate that the size of many pores also decreases from openings to constrictions, similar to our observations in Figure 2a.
Using the image processing open-source software Fiji (Schindelin et al., 2012), we conducted a statistical analysis of the angular probability distribution (PD) of 1,154 pore cone angles α in three non-adjacent images of quartz sand samples obtained from CT scans, as illustrated in Figure 2c.The probability of pore cone angle occurrence exhibits a sharp initial increase followed by a gradual decrease as the α increases.The pore cone angles are predominantly concentrated in the range of 10°-50°, with a cumulative probability of 0.9.The highest chance is observed in the range of 15°-20°, reaching 0.14 at this specific angular interval.It is evident that the pore geometry in realistic geomaterial samples is not simply parallel channels but can be better represented by a truncated cone shape.Therefore, in this work, we simplify the three-dimensional pores to a truncated cone shape to carry out numerical simulations for studying pore clogging and unclogging in a more realistic pore geometry, as illustrated in Figure 2d.

Model Configuration and Simulation Conditions
To examine the influence of pore shape parameters on pore clogging and unclogging, the diameter of the pore constriction D c and the length of the pore channel L c can be altered for the truncated cone pore; note that the variation of cone angle α depends on the variations of D c and L c .A decrease in L c increases the cone angle α (Figure 3a), which increases the normal contacts provided by the pore wall in the streamwise direction and hence hinders particle migration more significantly (Figure 3b).On the other hand, an increase in D c leads to a decrease in α and an increase in the number of required particles for pore clogging, jointly reducing the probability of clogging (Figure 3a) (Lai et al., 2001;Valdes & Santamarina, 2008;Zuriguel & Garcimartín, 2020).In this work, YIN ET AL. 10.1029/2023WR034628 5 of 29 to explore the impact of pore shape on pore clogging and unclogging, the pore cone angle α and pore constriction diameter D c are selected as the major geometric parameters.
As shown in Figure 3c, the pore structure is modeled as a cone-shaped wall and interacts with both DEM particles and the flow field solved by LBM.At a distance of 2D p from the pore entrance, 125 particles with a diameter D p of 1.0 × 10 −3 m are produced.The fluid flows along the positive direction of the x-axis.The position of the pore  YIN ET AL.

10.1029/2023WR034628
6 of 29 inlet is kept unchanged, while the position of the pore constriction is changed as we vary the length of the pore channel L c .
Previous research shows that the probability of pore-clogging is proportional to the size-ratio R o = D c /D p of pore constriction to the particle rather than their absolute sizes, and it is defined as a probability event with R o varying from 3.0 to 5.0 (Lai et al., 2001;Valdes & Santamarina, 2008).Moreover, our simulations indicate that the percentage of cases with pore-clogging is low, about 4%, when R o = 3.8, so we do not carry out simulations with R o > 3.8 in this research to focus on clogging behaviors.The values of R o that we use are 3.0, 3.2, 3.5, and 3.8, as listed in Table 1.To systematically vary α, the size-ratio R c = L c /D p of pore channel and particle is set as 2.0, 2.5, 2.9, 3.0, 3.5, 4.0, and 4.5, leading to α values varying from 13.7° to 36.9°.
Equation 1 illustrates that the chance of pore-clogging increases with decreasing particle velocity, and for suspended particles in the fluid, the particle velocity is controlled mainly by the fluid velocity.Coarse particles with diameters ranging from 1 to 5 mm, classified as sand, are commonly used in seepage erosion experiments, and the permeability coefficient for these particle samples typically falls from 1.0 × 10 −4 to 0.1 m/s (Chen et al., 2016;Ouyang & Takahashi, 2015;Won et al., 2019).Therefore, Poiseuille flow with v in as the center velocity is generated and the value of v in is set as 4D p /s, 6D p /s, 8D p /s, and 10D p /s in this study, which is located within the interval of seepage velocity under a conventional hydraulic gradient.As our preliminary simulations indicate that clogging occurs more frequently with R o = 3.0, we include more simulations with v in = 5D p /s, 7D p /s, and 9D p /s to further isolate the influence of fluid velocity.As listed in Table 1, the initial particle Reynolds number Re p = v ave D p /υ (where υ is the fluid kinematic viscosity) is calculated using the average boundary velocity v ave as the velocity difference between a particle and a fluid.Besides, R o , R c , and Re p are used to refer to different cases.

Fluid and Particle Properties
The LBM-IMB (D3Q19 model) with 19 components f i (x, t) and 19 lattice vectors e i is used in this study to simulate three-dimensional (3D) fluid-solid interactions, as described in greater detail in Appendix A. The fluid properties (density and viscosity) are set according to those of water at room temperature (Table 2).In DEM simulations, the DEM particles are assumed to be suspended (neutrally buoyant) during migration in the soil matrix (L.Zhang, Peng, Chang, & Xu, 2016), as we assume that the fluid injection induced particle migration is much faster than the particle settling due to gravity according to the Stokes formula (Stokes, 1850).Moreover, the particle radius we use is 5.0 × 10 −4 m, greater than typical clay particles (diameter usually is less than 7.5 × 10 −5 m) migrating inside soils (P.Cui et al., 2014), because using larger particles has the advantage of higher computational efficiency.Nevertheless, since the effect of the pore-to-particle size ratio on pore clogging and unclogging is prominent, the absolute value of the particle size does not affect our results.
Overall, the results of 110 simulations (Table 1) are reported, which were run on three high-performance computers with CPUs with 16 cores and 64 GB of RAM, with two controlling variables including the pore cone angle (adjusted by constriction size and pore length) and the velocity at the fluid inlet.

Result Validation
To examine the feedback process between the fluid pressure drop and the particle behaviors, including pore clogging and unclogging, Bianchi et al. (2018) experimented with the circulating injection of particle suspension with a constant concentration into a granular filter.As seen in Figure 4, the pressure difference between the inlet and outlet of the granular filter, ∆p, exhibits a rapid decline in the early stage and then transitions to a moderate decline and finally to a nearly constant value throughout pore unclogging.The experimental data are widely utilized to validate the accuracy of numerical simulations, including the LBM-DEM (Figure 4) and the coupling of Navier-Stokes equations and convection-diffusion equations for exploring the process of pore-clogging at the pore scale (Jäger et al., 2017(Jäger et al., , 2018)).Because the change in macroscopic ∆p can be considered as a collection of multiple pores clogging or unclogging at the pore scale, and its magnitude is proportional to the number and size of pores clogging or unclogging, macroscopic monitoring data can be used to compare with the microscopic results of the changes in ∆p.Based on this understanding, we use the experimental data set of Bianchi et al. (2018) to validate our LBM-DEM-IMB simulations for the single-pore configuration.
The pressure difference ∆p between two pressure monitoring cross sections S1 and S2 (Figure 3) is measured for R o = 3.0, R c = 3.5, and Re p = 1.9, where ∆p is determined by averaging the fluid pressure over all lattices on the two cross sections.Moreover, the pressure difference ∆p and the time t are normalized as ∆p* and t* using Equation 3: YIN ET AL.

10.1029/2023WR034628
8 of 29 where the subscripts max and min represent the maximum and minimum values of a variable ψ, respectively, and ψ* is the normalized version of ψ.
As shown in Figure 4, following unclogging, ∆p* initially decreases and then steadily stabilizes, which follows a power function with an index of −1.96.Our simulation results of ∆p* exhibit fluctuations (which are obtained from single-pore simulations), whereas the curves of the experiment and microscale simulation results are smoother.The fluctuations are likely associated with the random action of particle collisions and particle-wall interactions.Nevertheless, regarding the average behavior, Figure 4 shows that the fitting curve of ∆p* agrees well with the results of other researchers, which validates the numerical method used in this work in simulating the interactions between particles, fluids, and pores.

Particle Migration Phenomenology and Clogging Regimes
Under continuous fluid injection, particles may continuously flow through the pore, or the pore may become clogged by particles.However, pore-clogging does not necessarily result in the shutdown of the fluid-particle system; unclogging may occur as a result of an increase in the fluid pressure drop.Furthermore, clogging may reoccur after unclogging, leading to consecutive clogging and unclogging sequences.In the following, two typical particle migration sequences are presented before we identify several clogging regimes for different simulation conditions.The time t is nondimensionalized by the ratio of the x-length of the fluid domain to the average velocity of the inlet boundary v ave as  t , which represents the proportion of the fluid volume in the fluid domain that is replaced.The particle velocity v p is nondimensionalized by v ave as    , which represents the hydrodynamic force exerted on the particle as a result of the velocity difference between the particle and the driving fluid.In addition to the clogging and unclogging behaviors shown below, particle motion characteristics during continuous washout (no clogging) are depicted in Figure B1.

A Single Clogging-Unclogging Event
As illustrated in Figure 5, the scenario where R o = 3.0, R c = 3.5, and Re p = 1.9 is chosen to examine the characteristics of particle motion during particle migration, pore-clogging, and unclogging.The particle flow is inter- YIN ET AL. 10.1029/2023WR034628 9 of 29 rupted due to weak inter-particle support at  t = 0.27 (Figures 5a1 and 5a2).In contrast to the occurrence where particle velocity increases as the particles are closer to the outlet at  t = 0.27 in the case of continuous particle loss with α = 18.4° (Figure B1), the particle velocity in the pore is often lower than the v ave (    < 1 ) in this case with α = 23.2°.This suggests that the pore-clogging effect is enhanced at higher α.Furthermore, the percentage P e of particles washed out versus  t plot in Figure 5f shows that particle flow interruption occurs at least three times during  t from 0.2 to 0.3, and this increased frequency after  t = 0.1 seems to imply that pore-clogging is going to occur.
However, the particle structure causing the interruption is not stable because it cannot withstand the hydrodynamic action of the fluid.As the particles are washed out, the remaining particles become spatially sparse, and the particle velocity increases and eventually decreases with increasing    (Figures 5b1 and 5b2).The creation of a particle structure strong enough to defy hydrodynamics leads to pore clogging, as depicted in Figures 5c1  and 5c2, where the particle velocity approaches zero near the pore constriction.However, particles close to the pore inlet still travel at velocities much slower than the mean fluid velocity.Figure 5h shows that the    , the ratio of the average particle displacement in the pore to the particle diameter, increases slowly during pore clogging, indicating a compressed clogging structure.As can be seen in Figures 5d1 and 5d2, the majority of the particles come to a stop, but the pore-clogging is about to burst as the particles at the front of the clogging structure are washed out with    ≥ 6.0 .
As seen in Figure 5e1, the pore eventually becomes unclogged.Figure 5e2 also demonstrates that, in contrast to the gradual increase in    when    is less than 0.6 at  t = 0.33 , the particle velocity increases linearly along the direction of fluid motion at  t = 3.95 .The ratio    , the average velocity of particles in the pore to the v ave , is higher after unclogging than before clogging (Figure 5g).These findings suggest that the particles experience a greater hydrodynamic impact following pore unclogging, resulting in a higher average particle velocity.Particles are eventually washed away quickly.

Two Sequential Clogging-Unclogging Events
Clogging may reoccur after unclogging.Generally, after the first clogging event, particles continue to migrate slowly to the outlet under the influence of hydrodynamics.During this interval, individual particles may be washed out if they are incapable of resisting hydrodynamic forces due to the structural adjustment during clogging (Figures 6a1 and 6a2), which might be accompanied by a slight increase in P e as shown in Figure 6f.However, this rearrangement of particles does not necessarily cause unclogging.As the particles in front of the clogging become less tightly packed, localized unclogging then occurs, as depicted in Figure 6b1.At the same time, the particle velocities increase, but only a handful of particles have    reaching a higher value between 2.0 and 4.0, while the remaining particles all have    lower than 2.0 (Figure 6b2).Then, as these front particles are washed out, the remaining particles are not accelerated due to a stable structure gradually forming again in the pore (Figures 6c1 and 6c2), leading to a second clogging.
Figure 6h demonstrates that the average particle displacement    does not change much during the second clogging and is significantly lower than the first clogging.This finding suggests that the second clogging structure is conducive to the fluid flowing through and experiences less hydrodynamic pressure.Before the second unclogging, the particle velocity    increases along the    , indicating that the clogging structure is becoming unstable (Figures 6d1 and 6d2).Subsequently, the unclogging occurred again.Although the average particle velocity    is close to the first unclogging (Figure 6g),    of most particles at this time is almost larger than zero (Figure 6e1), and the closer to the outlet, the higher the particle velocity (Figure 6e2), which indicates that the possibility of reoccurrence of clogging is very low.
These results indicate that the first clogging occurs until the particle flow is interrupted and a stable intergranular structure is formed to resist the hydrodynamic force.The lower the particle velocity after the first unclogging, the more favorable the formation of the second clogging.Moreover, during the clogging, the particle structure of the second clogging is more difficult to compress than that of the first clogging.

Occurrence of Clogging and Unclogging in Parameter Spaces
In Figure 7, we present the occurrence of clogging and unclogging in a parameter space defined by the pore geometry (R o and tanα) and the hydrodynamic condition (Re p ).Among the 110 reported numerical simulations, clogging occurs in 37 cases, and pore clogging-unclogging occurs twice in five cases.In four cases, however, no unclogging is detected; yet, unclogging is possible if sufficient time is provided for these simulations.The 38 incidents, including one whole process of pore clogging and unclogging, have provided sufficient samples for our subsequent analysis.
As shown in Figure 7a, the probability of clogging P c decreases as the size ratio R o increases.At R o = 3.8, the P c reaches its lowest value of 0.04.P c is plotted in Figure 7b as a function of tanα, which is calculated by counting clogging occurrences for tanα ranging between 0.25 and 0.75 with an interval of 0.1.Similar to the relationship between P c and tanα in dry hopper flow (López-Rodríguez et al., 2019), P c increases with increasing tanα.Alternatively, the cases of pore-clogging with Re p equal to 1.9, 2.8, 3.8, and 4.7 are counted to calculate the clogging probability P c (Figure 7c).Cases with Re p = 2.3, 3.3, and 4.2 are omitted from the calculation due to insufficient data.The value of P c fluctuates obviously with increasing Re p because the velocity of particles, which directly determines the clogging probability according to Equation 1, is influenced not only by the fluid inlet velocity but also by the concentration of particles before clogging, tanα, and the strength of particle collisions (Jäger et al., 2018;Zuriguel & Garcimartín, 2020).1, considering the effect of particle velocity and R o .For suspended particles driven by low-velocity fluids, the inertia effect is small, so parameter a in Equation 1 can be set to zero.When a = 0, the constant b does not affect the ratio or relative magnitude of clogging probability for cases with different v m and R o , so b is set to 1.In our simulation, the particle velocity is determined by the average boundary velocity v ave , and the Re p is proportional to v ave .Therefore, we use  − 2  ln(Re) as a variable to show the likelihood of pore-clogging under fluid driving.Subsequently, we present the distribution of cases with and without clogging in a parameter space defined by tanα and  − 2  ln(Re) as shown in Figure 7e.The results indicate the existence of a boundary: when  −5 < − 2  ln(Re) ≤ −14.5 , the proportion of clogging cases reaches 42.1%, while when  −23 < − 2  ln(Re) < −14.5 , the proportion of clogging cases is 14.7%.This interval with  −23 < − 2  ln(Re) < −14.5 seems to be a transition region from a high probability of clogging to no clogging.Nevertheless, a study with a wider range of parameters of the pore and fluid is in the future needed to determine the position of these boundary lines accurately.

Fluid Pressure Drop During Particle Migration
The spatial distribution of particles changes the flow field which in turn impacts particle movement.In this section, we investigate the interaction between the fluid pressure drop ∆p and the particle migration dynamics.In the flow through porous media, a widely used model for the pressure drop is the Ergun equation, as shown in Equation 4, which indicates that ∆p is primarily dependent on the porosity ε of the particle packing, the superficial fluid velocity U, and the length L s of the seepage path (Ergun, 1952).
where μ is the fluid dynamic viscosity.
We utilize Equation 4 to determine the initial fluid pressure drop ∆p i when the particles are packed on a cubic lattice with ε = 0.476, L s = 5D p , and U = v ave in the initial state, and ∆p i is then used to normalize ∆p in subsequent time instances.As illustrated in Figure 8  The results of scenario (1) are presented in Figure C1, where ∆p/∆p i fluctuates with time and remains below 0.87 throughout the process of particle loss due to loose particle packing and the combined action of pore obstruction and fluid driving force on particle movement.
In scenario (2), the change in ∆p/∆p i is similar to that in scenario (1) before particle entry into the pore (Figure 8a).When pore clogging occurs, ∆p/∆p i increases sharply to 1.42 and pressure p increases at the pore inlet but declines dramatically at the outlet, as shown in the diagram of fluid pressure (  t = 0.36 ) in Figure 8a.The findings demonstrate that a clogging structure effectively blocked the free flow of fluids.The ∆p/∆p i continues to increase by 26% to 1.79 at  t = 3.87 as a result of particle migration, which further reduces the space available for fluid movement.After pore unclogging, the ∆p/∆p i declines significantly and the fluid pressure is rapidly discharged to the outlet, leading to a relatively uniform distribution of p at  t = 3.94 as illustrated in Figure 8a.Eventually, ∆p/∆p i approaches zero.
For scenario (3), the change of ∆p/∆p i with time is similar to the case with one pore clogging before the first unclogging, as shown in Figure 8b.However, the clogging period is shorter in this case, and p does not change considerably according to the pressure field in Figure 8b.After the first unclogging, which releases fluid pres-  The findings show that the fluid pressure drop is both a sensitive indicator of pore-clogging and one of the requirements for pore unclogging.To evaluate fluid conditions for pore unclogging, it is crucial to determine the extent to which pore clogging raises the fluid pressure drop.The Ergun equation (Equation 4) is widely applicable for determining the pressure drop of a fluid passing through a particle packing.In this study, the particle packing contains a tiny number of particles, and their spatial distribution cannot be considered to be uniform.There are frequently only a few particles in a loose state at the front and back of the clogging, resulting in a spatially substantial difference in porosity.Therefore, the system mean packing porosity is not well defined.To solve this issue, we calculate the porosity ε and the average surface velocity U for each section by dividing the flow field and particles into segments of 1.0% of the particle diameter perpendicular to the direction of fluid flow.Then, we calculate the fluid pressure drop ∆p s of each segment using Equation 5 and sum up the ∆p s from the position where the particle is present to the absence of the particle using Equation 6.Finally, the total drop in fluid pressure ∆p e of the particle packing is obtained.
where x s is the x-position of the cross-section.
where L p is the length of particle packing along the direction of fluid flow.
The fluid pressure drop ∆p si is calculated by monitoring the results of cross-sections S1 and S2 just before the imminent unclogging to verify the accuracy of the Ergun equation.In Figure 9, ∆p si is set as the x-coordinate, and YIN ET AL. 10.1029/2023WR034628 14 of 29 the fluid pressure drop estimated from the models above is plotted as the y-coordinate.The line y = x indicates perfect agreement between the model predictions and the simulation results.It can be found that Ergun's pressure drop ∆p e differs significantly from the simulation data.This is attributed to the fact that the Ergun equation, which ignores the wall boundary effect, is only suitable for the fluid pressure drop of the particle packing in the cylindrical container with a large R o (Winterberg & Tsotsas, 2000).
By contrast, the pore diameter gradually decreases along the fluid flow direction in this research, and a larger pore cone angle α means a decrease in the length and volume of the pore as well as in the outlet diameter, leading to enhanced resistance to the fluid flow because the wall provides streamwise resistance to particle motion and promotes particle jamming (Figure 3).Therefore, the pore cone angle α should be incorporated into the Ergun equation to improve how the pressure drop is calculated for this specific pore geometry.Based on the above analysis, α should generally have a positive impact on the fluid pressure drop, and when α = 0, the modified model should reduce to the classic Ergun equation.
To explore an empirical relationship between the enhanced pressure drop and the pore cone angle, we examine raw data from 38 clogged incidents as shown in Figure 10 and Figure C2 of Appendix C, and find the following relation to account for the influence of tanα on ∆p e and A.
where A is a fitting parameter that defines the sensitivity of the fluid pressure drop to α.As shown in Figure 10, the fluid pressure drop ratio ln(∆p e /∆p si ) increases with 1 − tanα, and the majority of the results lie within the region bounded by two linear functions y = 0.78x − 0.2 and y = 0.78x − 0.27.The result indicates that A is close to 0.78 in most simulations, as detailed in Figure C2 of Appendix C.However, when tanα = 0.75, A significantly deviates from the bounded region, as indicated in the blue box in Figure 10, which indicates that a larger pore cone angle α further enhances the resistance of clogging to the fluid.

∆𝑝𝑝𝑚𝑚 = ∆𝑝𝑝𝑒𝑒
(1 − tan )  = ∫ Figure 9 shows that ∆p m predicted using the modified Ergun equation agrees excellently with the measured data (red symbols fall closely on the y = x line).Note that ∆p m is computed by averaging all the A values for tanα = 0.75 and with A = 0.78 for other tanα.These findings show that when calculating the fluid pressure drop in porous media with many particles deposited and adsorbed on the pore wall or pore clogging, it is quite important to account for the pore geometry, particularly the pore cone angle.

Arch Formation and Micromechanical Analysis
The unclogging of pores requires not only the hydrodynamic force provided by the fluid but also the resistance of the clogged structure provided by inter-particle forces, f c .As shown in Figure D1 in Appendix D, the ratio f c /w p (where w p denotes particle weight) during clogging increases along the direction of fluid flow and can be separated into two clusters.The particles in the cluster with large values create an arch at the front of the clogging and carry the majority of the fluid load.This condition occurs in all 38 instances of clogging, as shown in Figure 11a.
Along the direction of the fluid motion, the normalized interparticle contact force   *  (using Equation 3) initially increases slowly and then abruptly.K-means cluster analysis is used to partition with a small value of less than 0.27.These data imply that a small number of particles bear the majority of the fluid load during the clogging event, which determines the unclogging.By isolating particles with a high   *  from the clogging structure, it is found that these particles exhibit four types of spatial distribution: a single arch (Figure 11c), a branch arch (Figure 11d), two separated arches (Figure 11e), and a dome structure (Figure 11f).And these particles with a high   *  in 38 instances of clogging are depicted in Figure D2.In 3D clogging, a dome is produced at the front of the clog, but it may be the 3D arch that is essential for clogging stability, with the surrounding particles bearing weaker forces but allowing the arch to limit deformation and avoid unclogging.When an arch or a dome is damaged by hydrodynamic forces, the clogging is freed.Thus,   *  in domes and arches should likewise have a significant impact on the unclogging.
In the current fluid-driven system, it is the hydrodynamic force that gives rise to the interparticle force, so the interparticle force should scale somewhat with the fluid pressure gradient.We calculate the fluid pressure ∆p si s p on the cross-section with the area s p through the center of a particle, followed by the ratio of ∆p si s p to w p to demonstrate that the ∆p si s p is equivalent to the dry particle weight.Figure 12 displays the results for the maximum f m and average f a of f c when unclogging impends.Both f m and f a increase linearly with ∆p si s p /w p , but the data for f m fluctuates more significantly with a smaller deterministic coefficient R 2 around the fitted linear curve.According to these findings, it is evident that the inter-particle force is highly dependent on the macroscopic fluid pressure, although f m maybe affected by other factors, such as fluid boundary conditions and pore geometry.Note We further consider the connection between f m and the two major simulation conditions, Re p and tanα, in Figure 13.The value of f m /w p increases with increasing Re p , which can be fitted by y = 1.42x − 0.96 with R 2 = 0.5.The relationship between its average f ma /w p and Re p can be described by y = 1.46x − 1.2 with R 2 = 0.9 (Figure 13a).Clogging requires a particle structure that can bear the drag force of the fluid, and this force increases with increasing Re p .At low f m , if an arch is produced, it will only create a temporary particle flow interruption.The relationship between f m /w p and tanα follows a power function y = 1.32x −1.2 with R 2 = 0.27, and a similar relationship can be observed between f ma /w p and tanα with higher R 2 (Figure 13b).Recall that Figure 7b demonstrates that larger tanα is more prone to pore clogging.These results indicate that the more favorable the conditions are for pore clogging, the less stable the particle structure is and the lower the inter-particle force before unclogging.The impact of the fluid pressure drop on f m /w p is also considered, as shown in Figure 14b.Significant correlation (R 2 = 0.85) is found between f m /w p and (tanα/Re p )/(∆p si s p /w p ), and their relationship can be described as a power law y = 1.84x −0.57 .Therefore, an empirical function of f m is obtained as follows: where a m and b m are constant (a m = 0.015, b m = 0.57).

Models of Fluid Pressure Drop and Critical Interparticle Force for Pore Unclogging
Consistent with previous studies (Bianchi et al., 2018;P. Cui et al., 2014;Jäger et al., 2017Jäger et al., , 2018;;Kermani et al., 2020Kermani et al., , 2021;;Yin et al., 2021), we find that unclogging is associated with an increase in fluid pressure drop following clogging.More importantly, the influence of pore geometry on fluid pressure drop ∆p, which cannot be described by the Ergun equation, has not been well studied; hence, we make further modifications to the Ergun equation to more precisely compute the ∆p taking into account the effect of pore shape, as shown in Figure 9.
When the pore size is modest in comparison to the fine particle size, as it is in rock and soil, this adjusted equation works better.When this ratio is sufficient to disregard the non-cylindrical pore geometry, the Ergun equation should be more applicable (Ergun, 1952;Winterberg & Tsotsas, 2000).When tanα is larger, though, the channel length may decrease or the diameter of the outlet may become smaller, and the regulation of the constriction size on the pressure drop becomes more significant.Modifying the parameter A to adjust the influence of tanα can make the modified Ergun equation applicable.In this study, the pore cone angles of all numerical simulations are restricted to a range of 13.7°-36.8°due to the extensive computation required.However, the value of parameter A is obscured by the infinity-trending expression of 1/(1 − tanα), which occurs when the pore cone angle is 45°.Therefore, more research is needed to determine how to calculate fluid pressure drop with a larger pore cone angle during pore clogging.
From the viewpoint of the resistance of the particles to the fluid flow, the spatial distribution of inter-particle force reveals that only a small number of particles bear the majority of the hydrodynamic force.Therefore, the critical condition that best characterizes pore unclogging is the maximum interparticle force f m , which could be a key point for us to study pore unclogging under a driving fluid.Furthermore, we find that f m is governed by the fluid flow velocity and the pore shape and decreases with the increase of tanα/Re p , indicating that when conditions are conducive to clog formation, they can also be easily unclogged with low hydrodynamic force.
By factoring in the fluid pressure drop, we can obtain the equation , which we may use to determine f m .The link between inter-particle force and fluid pressure is described by both this equation and YIN ET AL. 10.1029/2023WR034628 18 of 29 Equation 2, but the underlying physical mechanisms are very different.Equation 2 highlights that the pressure condition for unclogging is defined by the adhesion force (Jäger et al., 2017(Jäger et al., , 2018)).The clogging corresponding to Equation 2 is forcibly generated by adhesion and cohesion forces, and there is no requirement for the pore shape.In our equation, clogging is a probability associated with the establishment of a stable particle packing arrangement, which is difficult to form in cylindrical pores.In addition, our equation highlights the fact that the unclogging process is governed by several particles in front of the clogged structure.Nevertheless, the two equations represent two primary reasons for clogging: molecular forces and arch bridge structure, and they do not contradict one another.

Implications for the Calculation of Subsurface Fines Transported by Water
Studies in the past typically treated pore clogging and erosion as separate processes, paying less attention to the evolution of clogging, unclogging, and erosion (Bernatek-Jakiel & Poesen, 2018;Beverloo et al., 1961;Gerber et al., 2018;Lai et al., 2001).However, erosion will be finished rapidly if the rock and soil have good erosion conditions like wide pores, low confining pressure, and a high hydraulic gradient (L.Zhang, Peng, Chang, & Xu, 2016), but erosion may not occur at all if pore clogging is the dominant process.However, the permeability coefficient of rock and soil mass is greatly reduced by pore-clogging (Han & Kwon, 2023;Parvan et al., 2020;Wang et al., 2022), which both impedes fluid movement and causes soil subsurface erosion to develop relatively slowly, suggesting that the evolution process may account for the majority of the time during soil erosion.And because pore clogging easily occurs with a size ratio of the pore constriction to the particle of 3-5 (Gella et al., 2018;Valdes & Santamarina, 2008), the geometry of pores has an impact on fluid flow that can no longer be ignored (Dincau et al., 2022(Dincau et al., , 2023;;Jäger et al., 2017Jäger et al., , 2018;;Vani et al., 2022).The modified Ergun equation obtained at a single pore with a cone angle is the foundation for assessing the pressure response for pore clogging in porous media with complex pore structures.The simulation determines the fluid pressure drop just before pore unclogging occurs, which can become indicative of the fluid condition to eliminate clogging.Furthermore, preventing or encouraging clogging requires an appreciation of the underlying causes and structural mechanics of pore clogging.Our equation of critical interparticle force provides the conditional relationship between the clogging formation and the unclogging, and can preliminarily judge the stability of the clogging structure according to the given hydrodynamic conditions and pore characteristics.
More notably, it is essential to figure out how to adapt the modified Ergun equation (Equation 8) derived in a single pore to small-scale porous media elements to use this equation to tackle real-world challenges.The calculation of fine particle migration and clogging in larger-scale soil, such as the representative volumetric element scale, often neglects the influence of pore morphology.The intricate pore structure in porous media results in a complex flow field, causing clogging to be more stochastic (Dincau et al., 2023;Parvan et al., 2021).Therefore, certain assumptions are requisite prior to the extended application of Equation 8 for calculating clogging-erosion in porous media.It is commonly assumed that the migration and clogging of fine particles in the soil are calculated based on fine particle concentration, porosity, and Darcy's velocity.It is also assumed that the pore structure formed by the coarse particles is permanently fixed.We propose a parametric study framework considering pore morphology and using fluid pressure drop as the critical condition to conduct a saturated and larger-scale soil internal fine particle migration and clogging analysis as shown in Figure 15.
We assume the actual large-scale saturated soil is divided into many computational grids with size L g .Given the content and particle size of mobile fine particles in the soil and the spatial distribution of the pore structure of the skeleton particles obtained through CT scanning, we statistically analyze the size distributions of pore throats, pore inlets, the length L t between the pore inlet and pore throat, and pore cone angles.Based on pore morphology parameters, the pore cone angle can be computed.We specifically focus on pores whose extension direction aligns closely with the direction of fluid movement.The clogging probability is determined based on the clogging probability model for hopper flow under dry conditions.
At the initial moment, fluid pressure acts at the boundaries and varies along the flow direction according to Darcy's law.The permeability coefficients of each grid are calculated based on the Kozeny-Carman equation (Carrier, 2003;Parvan et al., 2021).It is assumed that fluid pressure at the inflow boundaries of individual grids decreases along the flow direction, resulting in a pressure drop ∆p g between the inflow and outflow boundaries.
YIN ET AL.The parametric study framework adopts numerous assumed conditions; however, the actual pore distribution within the soil is considerably more complex.There is still a long way to go before effectively applying models for fluid pressure drop in pores to real soil, considering its intricate nature.
The two formulas can be used as a benchmark for fluid pressure drop to unclog pores in the aquifer to achieve high-efficiency artificial groundwater recharge.They can provide theoretical references for the quick removal of pollutants from the soil.Moreover, a pressure drop threshold can be calculated to mitigate internal erosion in hillslopes, stream banks, and earth-rockfill dams, and a more accurate prediction of the stability of these large-scale rock and earth masses can be carried out.

Conclusions
Extensive research has been conducted on the clogging of particles in a flow through a constriction driven by gravity, but the interplay between the fluid flow and the clogging or unclogging of the pore in a fluid-driven system has received less attention.Here, we simulate particle migration through a single pore driven by a continuous fluid flow at the pore scale.The effects of the pore cone angle and the particle Reynolds number are systematically studied.The numerical simulation results are validated by comparing the change in fluid pressure drop over time between the simulation and experimental data following unclogging.In addition, we conducted a comprehensive analysis of the macroscopic fluid pressure drop caused by particle clogging and the microscopic inter-particle force statistics in the clogging structure.In summary, three conclusions are drawn carefully: 1.The pore cone angle α not only influences the probability of pore clogging but also makes the Ergun equation no longer suitable for the calculation of the fluid pressure drop induced by the clogging structure.The measured pressure drop ∆p si in the numerical simulation is always larger than the Ergun equation calculation unless the calculated result is divided by (1 − tanα) A to account for the amplification of the pressure drop by the pore However, many additional factors affecting clogging and unclogging are not included in this study, such as non-spherical particles, irregularly shaped pores, and a wider range of fluid velocity.Furthermore, our equations are all derived based on the single-pore model.How to extend the model on the single-pore scale to porous media is still a challenging question to be addressed in future work.
the particles in the central region is higher than that of the surrounding particles.Therefore, the particle velocity exhibits an oblique distribution along    (    is the ratio of the x-position of a particle to the length of the fluid domain).
As the pore diameter decreases gradually in the direction of fluid motion, the fluid velocity near the outlet increases.Consequently, more than half of the particles have    greater than 2.0 while entering the pore, as shown in Figures B1b1 and B1b2.In addition, the pore hinders the migration of particles on the outside of particle pack ing, resulting in a very low    , even near zero.
Continuous fluid begins to wash away the particles.As illustrated in Figures B1c1 and B1c2, the    along    initially climbs and then falls, with the turning point with the particle velocity being greater than five times the v ave as this particle approaches the outlet location.However, the    of the vast majority of particles is below 2.0 and even near zero.The data reveals that particles tend to develop an arch-bridge structure in the pore; however, this structure is unstable, and particles are eventually washed away.This process diminishes the continuity of the particles being washed away, causing the particles to appear spatially separated along the flow direction, as seen by the blue circle in Figures B1c1 and B1c2  22 of 29 constant, indicating that particle flow interruption is common throughout the particle loss process, analogous to particle flow in a hopper under dry conditions (Wu et al., 1993).
Figures B1d1 and B1d2 illustrate that as more particles are washed out, fluid resistance diminishes, and particles near the outlet obtain higher velocities with    larger than 6.0.Massive particles have been washed out, and the rate at which particles are washed out gradually decreases in the final stage, as shown in Figure B2, but a handful of particles remain at the pore inlet or on the pore wall without attaining a higher velocity due to pore obstruction, as shown in Figures B1e1 and B1e2.In addition, the current value of  t equals 0.54, indicating that half the volume of the fluid domain is sufficient to wash away the majority of particles.

Appendix C: Fluid Pressure Drop for Particles Continually Washed Out
For scenario (Equation 1) with particles continually washed out, ∆p/∆p i fluctuates with decreasing amplitude when particles migrate to the pore inlet (  t ≤ 0.14 ) (Figure C1).Because there is a substantial difference in fluid velocity at the inlet and outlet of the fluid domain, the particles gradually approach the fluid velocity under the action of hydrodynamics and their inertia.And when ∆p/∆p i ≈ 0, the particles follow with the fluid, and their effect on the fluid pressure becomes insignificant.Subsequently, as the particles aggregate in the direc tion of decreasing pore diameter, the space within the pore to transport fluid decreases, and the pores hinder particle movement, resulting in a decrease in particle velocity and, consequently, an increase in ∆p/∆p i during  0.14 < t ≤ 0.5 (Figure C1).But these particles are accelerated again and washed out, resulting in a drop in ∆p/∆p i accordingly.Therefore, the ∆p/∆p i remains below 0.87 throughout the particle loss process, suggesting that the particle packing is even looser than the cubic packing.
Parameter   and its average value as a function of the pore cone angle

Appendix D: Inter-Particle Force for Clogging Structure
Figure D1a shows that when the particles are constantly washed out and pore unclogging occurs with R o = 3.0, R c = 3.5, and Re p = 1.9, the ratio f c /w p (w p represents particle weight) is close to the zero axis, indicating that the particle packing is loose.After pore clogging, however, f c /w p climbs as    increases, and it exhibits two clusters with two distinct regions: f c /w p < 1.25 and f c /w p > 1.25 (marked by the red rectangular box).
Figure D1b shows the PD of the f c /w p by combining the data of f c /w p every 500 ∆t d for the 5,000 ∆t d before the pore unclogging occurs.When f c /w p > 1.25 is considered, the corresponding PD value is only 0.02 (marked by the blue rectangular box).A closer examination reveals, however, that f c /w p > 1.25 forms at the forefront of clogging and particles form an arching structure (Figure D1a particles marked by red) close to the pore constriction.These findings suggest that just a small number of particles with f c /w p > 1.25 carries the bulk of the load and determines whether or not an unclogging happens.
The interparticle force f c is divided into two clusters by K-means, and the particles bearing the main load are obtained by extracting the cluster with the larger f c .The spatial distribution of these particles can be divided into four types: a single arch, a branch arch, two separated arches, and a dome structure for 38 clogging incidents.A single arch and branch arch are more frequent, at 34.2% and 39.5%, respectively, while the type with two separated arches only accounts for 7.9% (Figure D2).

List of Notations
f b , clogging volume V c , and constriction diameter (cylindrical pore) D c was developed as given in Equation2

Figure 1 .
Figure 1.Particle migration and pore clogging are driven by the fluid in many hydrologic processes.(a) Scanning electron microscopy images of pores clogged by silicate fines during water injection in an aquifer's limestone sample (modified from Wang et al. (2022)).(b) Fine sand infiltrates into static gravel beds in flume experiments (modified from Gibson et al. (2009)).(c) Latex particles clog pores in several simplified pores (modified from Liu et al. (2019)).

Figure 2 .
Figure 2. Simplification of the three-dimensional pore of numerical simulations based on micropore CT images.(a) Binary CT image of quartz sand particle packing and extraction of two-dimensional pore structure with a spatial resolution of 9 μm obtained by the Shanghai Synchrotron Radiation Facility's BL13W1 X-ray imaging and biomedical applications beamline.(b) Two-dimensional pore structure in the binary CT image of sandstone (modified from Wetzel et al., 2020).(c) Probability distribution P a of pore cone angle α.(d) Simplified three-dimensional pore (D in is the diameter of the inlet, D c is the diameter of the constriction, L c is the length of the pore channel, and α is the pore cone angle).

Figure 3 .
Figure 3. (a) Effect of the length of the pore channel L c , the diameter of the pore constriction D c , and the pore cone angle α on particle migration.(b) An increase in the pore cone angle α strengthens the resistance of the pore wall to particles.(c) Placement of pore, particles, and fluid pressure monitoring cross-sections S1 (x = 0.1D p ) and S2 (x = 15.9Dp ) in the flow field.

Figure 4 .
Figure 4. Numerical model verification by comparing the change of fluid pressure difference ∆p* of the case R o = 3.0, R c = 3.5, and Re p = 1.9 over time t* after pore unclogging with the results of previous experiments and microscale simulations.

Figure 5 .
Figure 5. (a 1 -e 1 ) Spatial distribution of particles with a different color representing    and (a 2 -e 2 ) distribution of    along    for the simulation case with R o = 3.0, R c = 3.5, and Re p = 1.9.(f) The percentage P e of particles washed out, (g) the mean velocity    , and (h) the average displacement    of particles not washed out versus  t .

Figure 6 .
Figure 6.(a 1 -e 1 ) Spatial distribution of particles with a different color representing    and (a 2 -e 2 ) distribution of    along    for the simulation case with R o = 3.2, R c = 3.5, and Re p = 4.7.(f) The percentage P e of particles washed out, (g) the mean velocity

Figure 7 .
Figure 7. (a) Probability P c of clogging versus R o , (b) P c versus tanα, and (c) P c versus Re p .Distribution of cases with and without clogging in a parameter space defined by (d) Re p and   −1  tan  , and (e) tanα and  − 2  ln ( Re ) .
and FigureC1in Appendix C, we compare the change of ∆p/∆p i with the particle migration state for three distinct scenarios: (1) particles are continually washed out with R o = 3.0, R c = 4.5, and Re p = 1.9, (2) pore clogging and unclogging occur once with R o = 3.0, R c = 3.5, and Re p = 1.9, and (3) pore clogging and unclogging occur twice with R o = 3.2, R c = 3.5, and Re p = 4.7.The negative values in the pressure field in Figure8are derived relative to the initial pressure setting and do not represent absolute

Figure 8 .
Figure 8. Changes of fluid pressure with particle state.(a) Particle migration with one occurrence of pore clogging and unclogging with R o = 3.0, R c = 3.5, and Re p = 1.9.(b) Particle migration with two occurrences of pore clogging and unclogging with R o = 3.2, R c = 3.5, and Re p = 4.7.Note that the particles in our simulations have a uniform particle size, and the varying particle size in the snapshots are due to clipping effects.
/∆p i sharply decreases to 0.42, climbs to 0.82, and then decreases to about 0.71 slowly.The ∆p/∆p i in the first clogging is 1.82 times the ∆p/∆p i in the second clogging, as the clogged structure contains fewer particles, as depicted in the fluid pressure diagram (  t = 0.78 ) (Figure8b).Eventually, a second unclogging occurs, and ∆p/∆p i sharply decreases to zero.

Figure 9 .
Figure 9.Comparison between the theoretical value of fluid pressure drop calculated by the equations of Ergun (∆p e ) and modified Ergun (∆p m ), respectively, versus the simulation results (∆p si ).

Figure 11 .
Figure 11.(a) Normalized interparticle contact force   *  versus normalized x-position   *  of contacts.(b) Probability distribution (PD) of   *  for 38 clogging events.Spatial distribution forms of particles with large   *  : (c) a single arch, (d) a branch arch, (e) two separated arches, and (f) a dome structure.

Figure 13 .
Figure 13.The f m /w p and its average f ma /w p change with (a) particle Reynolds number Re p and (b) pore cone angle tanα.

Figure 12 .
Figure 12.The ratio of maximum f m and average f a of inter-particle force to particle weight w p versus fluid pressure drop on the cross-section through the center of the sphere particle ∆p si s p /w p when unclogging occurs.

Figure 14 .
Figure 14.(a) Relationship between the f m /w p and the ratio of tanα to Re p and (b) relationship between the f m /w p and (tanα/ Re p )/(∆p si s p /w p ) follow the negative power functions with different exponents.
of particles within the pores based on pore length, and calculate the length of fine particles packing obstructing fluid flow in the direction of motion, taking into account their diameter, the pore packing fraction at close particle contact, and pore morphology.Ultimately, ∆p e and ∆p m can be determined according to Equations 6 and 8. Comparing ∆p m with ∆p g L t /L g helps to ascertain variations in the PD of clogging.If ∆p m is smaller than ∆p g L t /L g , the fine particles in the pore are transported into the adjacent grid in the direction of flow.If ∆p m is smaller than ∆p g L t /L g , the eroded mass of the fine particles in the pore is inversely proportional to the clogging probability.

Figure 15 .
Figure 15.Application of the modified Ergun equation for the fluid pressure drop in a single pore considering pore morphology in porous media.
Figure B1.(a 1 -e 1 ) Spatial distribution of particles with a different color representing    and (a 2 -e 2 ) distribution of    along    at five stages for the simulation case with

Figure B2 .
Figure B2.The ratio P e of the number of particles washed out to their initial amount changes with time  t .Particle separation occurs at locations marked by blue circles.

Figure C1 .
Figure C1.Changes in fluid pressure with particles continually washed out with R o = 3.0, R c = 4.5, and Re p = 1.9.

Figure C2 .
Figure C2.Parameter A and the average of A versus the pore cone angle tanα.

Figure D1 .
Figure D1.Inter-particle force characteristics with R o = 3.0, R c = 3.5, and Re p = 1.9.(a) Distribution of the ratio of the inter-particle force f c to particle weight w p along the direction of fluid flow    .(b) Probability distribution (PD) of the f c /w p .

Figure D2 .
Figure D2.Spatial distribution forms of particles with large   *  for 38 clogging incidents (① and ② are used separately to mark the first and second clogging for one case).

Table 1
Parameters of the Pore Geometry and Fluid in 110 Numerical Simulations c is negatively correlated with R o and positively correlated with tanα, indicating a positive correlation between P c and   −1  tan  .Therefore, to clarify the influence of pore geometry and fluid characteristics on clogging, we present the distribution of cases with and without clogging in a parameter space defined by Re p and   −1  tan  as shown in Figure 7d.The results show that there seems to be a boundary separating clogging and no clogging.On this boundary line, when 2.3 ≤ Re p ≤ 4.7,   −1  tan  decreases with increasing Re p , and when Re p ≤ 2.3,   −1  tan  may increase with increasing Re p .The probability of clogging can be calculated by Equation