Introduction

Tailings impoundments are designed to accumulate solid waste generated during ore refining processes. These impoundments are characterized by high potential energy, significant pollution risks, and elevated hazards. If a dam failure or leachate leakage occurs, it can result in substantial property damage and environmental contamination downstream. The filter, composed of coarse-grained soil materials such as sand, gravel, or pebbles arranged according to a specific grading, ensures efficient drainage, effectively reduces the phreatic line within the dam body, and enhances the permeability stability of the dam structure. However, over time, fine particles progressively enter the filter, leading to siltation both on its surface and within its structure. This process significantly reduces the filter’s permeability, causing a rise in the water level within the dam, potential leachate leakage, and a corresponding risk to both the dam’s structural integrity and the surrounding environment. According to statistics1, 42.86% of dam failures in tailings impoundments are attributed to the loss of permeability stability, which is one of the most critical factors affecting dam safety. Therefore, investigating the siltation characteristics of the filter and understanding the evolution of its internal pore structure and permeability is of significant engineering importance.

The clogging behavior of filters is fundamentally a process of fine particle migration and deposition within the filter under seepage forces, representing a typical unstable seepage problem in variable porous media2,3,4,5. During this process, physical parameters such as soil mass, pore structure, and hydraulic gradient undergo significant temporal variations, leading to a significant reduction in soil porosity and hydraulic conductivity6,7. Based on the theory of muddy water seepage under the context of high-sediment-content Yellow River water infiltration, this problem is characterized as one involving variable seepage paths and variable hydraulic conductivity2. The primary factors influencing filter clogging8,9 can be categorized into three major types: physical clogging10,11,12, chemical clogging12,13, and biological clogging13,14,15,16. Among these, physical clogging is the most common form and can be further subdivided into three mechanisms based on the geometric relationship between fine particles and filter pore sizes: surface interception, internal blockage, and particle attachment17,18. The physical factors influencing filter clogging are complex and can be broadly classified into internal and external factors.

The internal factors influencing filter clogging primarily include particle gradation19, particle morphology20,21,22, and pore geometric characteristics23. Among these, particle gradation directly affects the distribution and deposition depth of clogging materials within the filter24. When the coefficient of uniformity of the soil is large, the pore size distribution is wide, and fine particles in muddy water are more likely to fill the voids between larger particles, leading to localized internal clogging. Conversely, when the soil particles are uniform, the regular pore structure is more prone to forming surface bridging clogging25. Particle morphology significantly influences the clogging process26: rounded particles tend to pass through the pore structure more easily, while angular particles, due to their tendency to interlock, have a significantly higher clogging probability than rounded particles. Studies show that the clogging probability of non-spherical particles (e.g., disk-shaped, rod-shaped) can be more than three times that of spherical particles25,27,28. The size of the pore diameter directly determines the extent and spatial distribution of fine particle clogging within the pores29,30,31. Conical pore structures in sandy soils can promote close contact between migrating particles, forming stable arch-like or dome-shaped clogging structures23. Additionally, pore tortuosity is another critical factor. Pores with low tortuosity and high connectivity are less prone to clogging, whereas complex pore structures with high tortuosity and numerous dead-end pores accelerate the clogging process32.

The external factors influencing filter clogging primarily include hydraulic boundary conditions as well as the concentration and size distribution of fine particles in muddy water. The magnitude of the hydraulic gradient directly affects particle migration behavior. Higher hydraulic heads accelerate the transport velocity of particles within the soil column while promoting particle deposition at the top, forming a filter cake and significantly increasing the reduction in filter permeability33,34. Studies have shown that under high-pressure and high-flow-rate conditions, fine particles are more likely to penetrate the filter, resulting in a slower decline in filter permeability34. In contrast, under low-flow-rate conditions, fine particles tend to deposit within the filter35,36. For a specific filter structure, the particle size distribution in muddy water influences the clogging pattern through mechanisms such as bridging, migration, and filling. Coarse particles primarily cause rapid surface clogging, while fine particles lead to deep clogging37. Small-sized particles form mixed clogging due to Brownian motion and aggregation effects, whereas large-sized particles cause surface clogging due to gravity and interception38. Additionally, higher concentrations of fine sediment particles in water facilitate the formation of clustered structures, which block pore channels and significantly reduce permeability29.

The migration of particles within soil mass under hydraulic action represents a classic fluid–solid coupling problem. This research focus has led to the development of numerous numerical computation methods. Among these, the most prevalent is the CFD-DEM fluid–solid coupling method21,29. This approach simulates fluid flow by solving the macroscopic Navier–Stokes equations and achieves bidirectional coupling through momentum exchange between fluid and particles. The method enables high-precision capture of particle–fluid interaction processes, extensively supports complex particle shapes and non-Newtonian fluids, and has been widely applied in soil internal erosion and slurry clogging30,39,40. Secondly, the LBM-DEM fluid–solid coupling method41,42, based on mesoscopic particle collision and migration mechanisms, can capture microscopic flow details43. It demonstrates inherent advantages in handling complex boundaries (e.g., porous media) and non-Newtonian fluids. Additionally, the DEM-DFM method7 utilizing the dynamic fluid mesh (DFM) approach dynamically updates tetrahedral fluid meshes to capture synergistic effects between pore deformation and seepage. It calculates porosity and permeability in real-time based on particle positions (via the Kozeny-Carman equation) and solves the Darcy flow field (pressure–velocity field).

This study employs the CFD-DEM fluid–solid coupling approach to simulate the clogging behavior of filters. By controlling the size of fine particles entering the filter, the pressure difference between the fluid inlet and outlet, and the concentration of fine particles in the water, the spatiotemporal evolution of void ratio, hydraulic conductivity, and dry density during the clogging process is investigated. The study aims to reveal the evolution patterns of macroscopic physical and mechanical parameters from the perspective of microscopic particle behavior, providing a theoretical foundation for the rational design of filters.

Discrete element theory of fluid–solid coupling and method of programme implementation

Control equations for fluids

In the CFD-DEM fluid–solid coupling method used in this study, the fluid is primarily analyzed using the Dense Discrete Phase Model (DDPM), which has been extensively applied in the study of dense particle motion in fluidized beds. The advantage of the DDPM lies in its ability to account for the effects of particle volume on the fluid flow, as well as its capability to simulate particle collisions during motion. The governing equations for the fluid are the continuity equation and the momentum equation, as shown in Eqs. (1) and (2).

$$\frac{{\partial \left( {\alpha_{f} \rho_{f} } \right)}}{\partial t} + \nabla \left( {\alpha_{f} \rho_{f} {\varvec{u}}_{f} } \right) = 0$$
(1)
$$\frac{\partial }{\partial t}(\alpha_{f} \rho_{f} {\varvec{u}}_{f} ) + \nabla \cdot (\alpha_{f} \rho_{f} {\varvec{u}}_{f} {\varvec{u}}_{f} ) = - \alpha_{f} \nabla {\varvec{p}} + \alpha_{f} \nabla \cdot \left( {\mu_{f} \nabla {\varvec{u}}_{f} } \right) + \alpha_{f} \rho_{f} {\varvec{g}} + f_{f,s}$$
(2)

where \(\alpha_{f}\) is the fluid volume fraction; \(\rho_{f}\) is the fluid density; t is the time; \({\varvec{u}}_{f}\) is the fluid velocity; \(\mu_{f}\) is the continuous phase fluid viscosity; \({\varvec{p}}\) is the dynamic water pressure; \(f_{f,s}\) is the drag force of the fluid on the particles; \({\varvec{g}}\) is the gravitational acceleration.

Equation of motion for solid particles

The motion of particles is modeled using the Discrete Element Method (DEM), which is governed by equations of motion and rotational equations. The equation of motion accounts for various forces acting on the particles, including physical forces (such as gravity), the interaction forces between particles, and the combined drag force exerted by the fluid on the particles’ velocity and acceleration. The equation of motion for any given particle can be derived from Newton’s second law, as shown in Eq. (3).

$$m_{i} \frac{{{\text{d}}{\varvec{u}}_{i}^{{\text{p}}} }}{{{\text{d}}t}} = {\varvec{F}}_{i}^{{\text{g}}} + {\varvec{F}}_{i}^{d} + \mathop \sum \limits_{{\text{C}}} {\varvec{F}}_{ij}^{{\text{C}}}$$
(3)

where \(m_{i}\) is the mass of particle i; \({\varvec{u}}_{i}^{{\text{p}}}\) is the velocity of particle i; \({\varvec{F}}_{i}^{{\text{g}}}\) is the gravitational force acting on the particle; \({\varvec{F}}_{ij}^{{\text{C}}}\) is the interaction force between particle i and particle j; \({\varvec{F}}_{i}^{d}\) is the drag force of the fluid on the particle.

The equation of rotation of the particle when subjected to an external moment is presented in Eq. (4).

$$I_{i} \frac{{{\text{d}}{\varvec{\omega}}_{i} }}{{{\text{d}}t}} = \mathop \sum \limits_{j = 1}^{{n_{i}^{c} }} \left( {{\varvec{M}}_{{{\text{t}}ij}}^{c} + {\varvec{M}}_{{{\text{r}}ij}}^{c} } \right)$$
(4)

where, \(I_{i}\) is the moment of inertia of the ith particle; \({\varvec{\omega}}_{i}\) is the angular velocity of the particle; \(n_{i}^{c}\) is the number of particles in contact with particle i; \({\varvec{M}}_{{{\text{t}}ij}}^{c}\) is the moment generated by the tangential force (including inter-particle and fluid forces); \({\varvec{M}}_{{{\text{r}}ij}}^{c}\) is the rolling moment generated by rolling friction.

Modelling of particle–fluid interactions

The interaction force \({\varvec{F}}_{i}^{f,p}\) between the fluid and the particles is comprised of several components, including drag \({\varvec{F}}_{i}^{d}\), buoyancy force \({\varvec{F}}_{i}^{b}\), and pressure gradient force \({\varvec{F}}_{i}^{g}\), as illustrated in Eq. (5):

$${\varvec{F}}_{i}^{f,p} = {\varvec{F}}_{i}^{d} + {\varvec{F}}_{i}^{b} + {\varvec{F}}_{i}^{g}$$
(5)

The drag force \({\varvec{F}}_{i}^{d}\) experienced by a single particle in the flow field is a complex fluid–solid coupling problem, and currently, only empirical models are available for its calculation. In this paper, we use the free-sphere drag computational model to compute the force acting on the particles, as shown in Eq. (6)

$${\mathbf{F}}_{i}^{d} = 0.5C_{D} \rho_{f} A_{p} \left( {{\varvec{u}}_{f} - {\varvec{u}}_{p} } \right)\left| {{\varvec{u}}_{f} - {\varvec{u}}_{p} } \right|$$
(6)
$$C_{D} = \left\{ \begin{gathered} \frac{24}{{Re}}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \;Re \le 0.5 \hfill \\ \frac{{24\left( {1.0 + 0.15Re^{0.687} } \right)}}{Re}\quad \quad \;0.5 < Re \le 1000 \hfill \\ 0.44\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad Re > 1000 \hfill \\ \end{gathered} \right.$$
(7)
$$Re = \frac{{\alpha_{f} \rho_{f} L\left| {{\varvec{u}}_{f} } \right|}}{{\mu_{f} }}$$
(8)

where CD is the drag force coefficient; \(\rho_{f}\) is the density of the fluid; Ap is the projected area of the particles, \({\varvec{u}}_{f}\),\({\varvec{u}}_{p}\), is the velocity of the fluid and the particles; \(Re\) is the Reynolds number; L is the diameter of the particles.

An algorithmic implementation of discrete elements for fluid–solid coupling

In this paper, the CFD-DEM fluid–solid coupling analysis is carried out using ANSYS Fluent 2021R1 and EDEM 2022.2 software. The coupling between the two is achieved through the UDF (User-Defined Function) interface in Fluent, which is based on the Eulerian multiphase flow model. This allows for bidirectional transmission of the interaction forces between the fluid and the particles. The CFD-DEM solver flow process is illustrated in Fig. 1. At the beginning of the simulation, Fluent calls the UDF to obtain the particle position data from EDEM. This data is then used to compute the mesh porosity and the interaction forces between the fluid and the particles. The velocities and positions of the particles are updated according to the governing Eqs. (3) and (4). Once the particle data is updated, it is transferred back to Fluent, where the momentum and volume fraction of the fluid phase are solved. The fluid mesh and interaction forces are then automatically updated at the end of each convergence iteration, after which the computation proceeds to the next step.

Fig. 1
figure 1

CFD-DEM solving process.

Fluid–solid coupling calculation model

In this study, the methodology involves first generating a cylindrical filter specimen with a length of 10 cm and a diameter of 5 cm using discrete element modeling, as shown in Fig. 2. The soil particles of the filter are modeled as spherical particles with a uniform diameter of 5 mm. These particles are generated within the calculation area based on the specified mass and are then extruded simultaneously from the top and bottom rigid planes toward the center, forming the 10 cm high filter specimen. The specimen density is set to 1.5 g/cm3, and the initial porosity is 0.4. The filter is divided into ten layers, each 1 cm thick, labeled from the top as Layer 1 (L1) to Layer 10 (L10). This division facilitates the statistical analysis of fine particle deposition within each layer (L1–L10) and allows for tracking changes in the void ratio, hydraulic conductivity, and dry density across the layers.

Fig. 2
figure 2

Model schematic.

Secondly, a muddy water seepage model is applied above the filter specimen to simulate the inflow of fine particles from the dam body into the filter. In the DEM framework, a fine particle generation surface with the same diameter as the filter is positioned 5 mm above the top filter screen. Fine particles enter the filter under the combined action of gravity and fluid drag forces and disappear upon contacting the bottom surface if they penetrate the filter. To enhance computational efficiency, all particles are assumed to be spherical. For the CFD setup, the mesh size is set to 2–3 times the particle diameter (approximately 10–15 mm), resulting in a total of 832 grids. The specimen top is defined as a water pressure inlet, the bottom as a zero-pressure outlet, and the surrounding walls as impermeable. Notably, the term “muddy water” in this study refers to water containing fine particles, while the fluid phase is assumed to be Newtonian clear water (dynamic viscosity = 0.001 Pa·s) to isolate particle–fluid interactions. Based on the Rayleigh criterion, the DEM time step is set to 1 × 10−6 s, and the CFD time step to 1 × 10−4 s.

To investigate the influence of the particle size ratio (Ra) on filter clogging during muddy water seepage, this study focuses on the geometric compatibility between the filter and fine particles in muddy water. The filter particle size was standardized to a uniform gradation of 5 mm, eliminating interference from complex gradations on pore structure. By varying three parameters—fine particles size, hydraulic head at the specimen inlet, and fine particles concentration-we analyze the mechanistic roles of particle size ratio (Ra), seepage head, and particle concentration in governing the spatiotemporal evolution of void ratio, hydraulic conductivity, and dry density within the filter.

The particle diameter and density of the filter were kept constant, while the particle size of the fine particles in the muddy water, the top water pressure, and the concentration of fine particles in the muddy water were varied. This approach was used to study the effects of the particle size ratio between the filter material and the protected soil particles, the infiltration head, and the concentration of fine particles on the void ratio, hydraulic conductivity, and dry density of the filter. Additionally, the temporal and spatial evolution of these properties during the muddy water seepage process were investigated.

The numerical simulation conditions for the coupled CFD-DEM in this study are summarized in Table 1, and the material properties of particles and fluid are detailed in Table 2. All parameters were adopted from previous studies29,39,40,44,45. The concentration of fine particles is denoted by Cc (in g/s). The ratio of the particle diameter in the filter to the average particle size of the fine particles in the muddy water is referred to as the particle size ratio, denoted by Ra. When the term "particle size ratio" is mentioned in the following sections, it refers to this specific definition. Additionally, the water pressure at the top of the filter is denoted by Pr (in kPa). For example, AC-Ra6.25-Pr5 represents a working condition where the particle size ratio is 6.25 and the water pressure at the top of the filter is 5 kPa. C-Ra6.25-Cc25 refers to a condition where the concentration of fine particles in the muddy water is 25 g/s. The "A", "C", and "B" labels correspond to working conditions related to changes in the particle size ratio, concentration, and water pressure at the top of the filter, respectively.

Table 1 CFD-DEM numerical calculation working condition.
Table 2 Model parameter.

The effect of characteristic parameters of muddy water on the clogging characteristics of the filter

Effect of particle size ratio (Ra) on the filter

Analysis of the siltation state of the filter with different particle size ratios (Ra)

The particle size ratio (Ra) directly influences the siltation distribution of fine particles in the filter during the infiltration of muddy water. The siltation distribution of fine particles for different particle size ratios (Ra = 2.5, 3.3, 5.0, 6.25, 8.3, 10.0, 16.7) is shown in Fig. 3. From the early stage of infiltration (t = 0.1 s), as shown in Fig. 3(a), the depth of infiltration of the fine particles in muddy water varies with the particle size ratio. Specifically, the larger the particle size ratio, the deeper the fine particles penetrate into the filter. At an infiltration time of t = 1.0 s, as shown in Fig. 3(b), the fine particles with smaller particle size ratios (Ra = 2.0, 3.3) have essentially stabilized and no longer migrate downstream. In contrast, fine particles with larger particle size ratios (Ra = 10.0, 16.7) continue migrating downwards. After t = 4.0 s, as shown in Fig. 3(c), the fine particles for Ra ≤ 10.0 have stabilized, while for Ra = 16.7, a significant number of fine particles directly penetrate through the filter and continuously and stably move downward, eventually disappearing at the bottom.

Fig. 3
figure 3

Siltation process of fine particles in muddy water with different particle size ratios.

When the infiltration of muddy water stabilizes, the siltation mode is classified into surface siltation, internal siltation, and internal penetration based on the ___location of the fine particles in the filter. When the particle size ratio (Ra) is small (e.g., Ra = 2.5), the fine particles are unable to penetrate the pore channels between the particles in the filter. Most of the particles will accumulate and silt on the surface and top of the filter, which is referred to as surface siltation. As the particle size ratio increases (Ra = 3.3, 5.0, 6.25, 8.3), the fine particles can move freely through the pore spaces of the filter. However, most of the particles are retained within the filter’s pores. As the accumulation of fine particles increases, the seepage channels become increasingly blocked, and the infiltration of fine particles into the lower layers decreases. Eventually, the infiltration process reaches a stable state, which is known as internal siltation. When the particle size ratio (Ra) increases significantly (Ra = 10.0, 16.7), the pore diameter of the filter becomes relatively larger compared to the fine particles in muddy water. Under hydraulic transport conditions, a portion of fine particles can penetrate through the filter to reach the bottom escape layer, while others are intercepted in narrower pore constrictions. Nevertheless, interconnected seepage channels remain preserved, defining this state as internal penetration.

The effective pore diameter formed by the skeleton particles of the filter was calculated as follows46:

$$D_{e} = \frac{2}{{3\alpha_{1} }}\frac{n}{1 - n}d_{e}$$
(9)

where De is the effective pore size corresponding to the equivalent particle size of the skeleton particles; α1 is the shape correction coefficient of soil particles and capillary channels, which is the ratio of the surface area of the natural soil particles to the surface area of the sphere of the same volume, and in this paper, we take α1 = 1; n is the void ratio; de is the equivalent particle size of the skeleton particles.

Theoretically, when the ratio of the pore diameter of the filter skeleton particles to the fine particles in the muddy water is less than or equal to 1, the fine particles cannot pass through the pores of the filter, resulting in surface siltation. Conversely, when this ratio is greater than 1, the fine particles can penetrate through the pores. However, in practice, even when the pore diameter of the spherical particles forming the filter is large, the effective pore size is only an idealized average value. Therefore, even when the ratio is greater than 1, internal blockage can still occur. As a result, internal blockage is observed for ratios greater than 1, while a small number of fine particles may still pass through the filter when the ratio is less than 1. Table 3 outlines the siltation patterns for the filter. For different particle size ratios, the calculation of De/det under various particle size ratios helps to identify the corresponding siltation pattern, as shown in Table 4. Combined with the results from Fig. 3, it is evident that the siltation patterns for different particle size ratios are consistent with the calculation results in Table 3. According to Table 4, the internal siltation pattern is further subdivided into the surface-internal double-layer siltation pattern (I-S) and the internal siltation pattern (I). For the case of Ra = 3.3, it falls into the I-S pattern, which is consistent with the observation that fine particles in muddy water silt up mainly in the shallow layer under these conditions. Through inverse calculation using Eq. (9) and the critical values of the siltation states provided in Table 3, the upper and lower thresholds of Ra for fine particles to either enter the soil body and cause siltation or not (I-S, I) are 8.8 and 2.2, respectively. This means that for a filter with a particle size of 5 mm, fine particles with a diameter smaller than 0.57 mm will penetrate internally, while particles with a diameter greater than 2.2 mm will be retained on the surface, leading to surface siltation. When the particle size ratio (Ra) approaches critical thresholds, the clogging behavior of the filter exhibits a mixed state of two clogging modes near the threshold values. For instance, when Ra approaches 2.2, fine particles may form surface siltation while partially migrating into the filter to create localized bridging (resembling internal siltation). Similarly, when Ra approaches 8.8, internal siltation predominates, yet some particles may penetrate the filter. For example, at Ra = 10.0, particles should theoretically penetrate the filter, but during actual infiltration, most particles clog internally, with only a few penetrating through preferential flow channels (as shown in Fig. 3). This hybrid behavior indicates that the clogging mechanism in transitional zones results from competition between geometric constraints and hydrodynamic forces.

Table 3 Determination of the siltation pattern of the filter47.
Table 4 Calculation of siltation pattern discrimination for different particle size ratios.

The effect of the particle size ratio (Ra) on the porosity of the filter

After computational processing, DEM software quantifies the volume of fine particles in each layer and coarse particles in filter at discrete time intervals. This dataset enables precise determination of void ratio distribution surfaces across temporal and spatial domains, as depicted in Fig. 4. The red color in the figure indicates a larger decrease in the void ratio, while blue indicates a smaller decrease. The void ratio forms a quadrilateral surface within the spatiotemporal ___domain, exhibiting a progressive decrease toward the tail of the time axis and the origin of the spatial axis. It can be observed that the temporal and spatial distribution of the void ratio is similar when Ra ≤ 10.0, with the largest decrease in the void ratio occurring near the top of the sample as the infiltration time increases. Additionally, as the particle size ratio (Ra) decreases, the void ratio exhibits a steeper slope in both the time and depth directions. This suggests that a smaller particle size ratio results in a quicker attainment of siltation stability and a shallower penetration of fine particles into the internal filter. When Ra = 16.7, a significant number of fine particles penetrate the filter, and the changes in its internal pore structure become more complex.

Fig. 4
figure 4

Spatial and temporal distribution of void ratios for different particle size ratios.

By projecting the void ratio surface in Fig. 4 along the depth direction, the variation of void ratio with time at different depths can be obtained. Taking the conditions of Ra = 10.0 and Ra = 16.7 as examples, as shown in Fig. 5(a). When the particle size ratio Ra is small, with the continuous entry of fine particles into the filter, the void ratio in each layer starts to decrease monotonically to varying extents. The decrease becomes progressively smaller from the top layer to the bottom layer, and there is no intersection between the layers. Additionally, the pattern remains consistent when Ra ≤ 10.0. Combining the particle movement state from Fig. 3, it can be observed that the fine particles in the muddy water first start forming clusters at the top of the filter, causing clogging in some of the pore channels of the filter. Subsequent particles, if they fall above the cluster, will cause clogging; if they enter the pore channels, they will continue moving downward until being intercepted within the filter. After several layers of fine particles are captured by the filter, the number of pore channels leading to the lower layers progressively decreases, resulting in a stable clogging state within the filter. As a result, the void ratio inside the filter decreases progressively with each layer, never intersecting. When the particle size ratio Ra is large (Ra = 16.7), as shown in Fig. 5(b), the pore channels in the filter are ineffective at intercepting the fine particles in the muddy water, and the fine particles pass directly through the entire filter under hydraulic transport. The void ratio of the L1–L2 layers begins to decrease and stabilize first, with the smallest decrease in void ratio. The stabilization of the void ratio in layers L1–L2 is dynamic, meaning that the fine particles entering and exiting the layer have reached a balance. Subsequently, the L3–L5 layers begin to activate and gradually stabilize, with the void ratio decreasing further, indicating that they contain not only the fine particles in the muddy water but also some fine particles that have flowed downward from the upper L1–L2 layers. The reason for the further decrease in void ratio for the L6–L8 layers is the same, but their activation occurs later compared to the L1–L5 layers. For the L9–L10 layers, due to their proximity to the bottom outlet, the number of fine particles flowing out is greater than in the L6–L8 layers, so the void ratio slightly increases. Figure 5(c) illustrates the temporal evolution of the void ratio in the first filter layer under different particle size ratios (Ra). During the initial infiltration stage, the void ratio decreases rapidly across all Ra values, characterized by steep curve slopes. Over time, the slope gradually diminishes until stabilizing horizontally, indicating the formation of a stable clogging layer at the top. For Ra ≤ 10.0, larger Ra values correspond to smaller void ratios in the first layer, as most fine particles in muddy water form clogging near the filter surface. For example, at Ra = 10.0, the void ratio in Layer L1 decreases by 27.2%, representing a 14.3% greater reduction compared to Ra = 3.33. However, at a particle size ratio Ra = 16.7, the reduction magnitude of the void ratio in the first layer becomes smaller. This occurs because the fine particles in the muddy water are too small to be effectively intercepted by the filter, resulting in partial particle migration to the middle and lower layers.

Fig. 5
figure 5

Variation of void ratio with time.

Figure 5d shows the variation trend of the total void ratio of the filter with time under different particle size ratios. The results demonstrate that during the initial clogging stage, the total void ratio of the specimen decreases sharply, with the rate of decline gradually slowing over time. A larger particle size ratio (Ra) corresponds to a delayed entry into the stable clogging phase and results in a smaller final total void ratio. For instance, when Ra = 3.33, the void ratio stabilizes at approximately 1 s with a final value of 0.657, whereas for Ra = 16.7, stabilization occurs after 2.5 s with a final void ratio of 0.580. The following generalized equation is proposed to describe the variation relationship of the void ratio over time under different particle size ratios (Ra) :

$$e(t) = e_{\infty } + (e_{0} - e_{\infty } )e^{ - \lambda t}$$
(10)

where e0 is the initial void ratio, which is 0.667 in this study; e is the void ratio at steady state; λ is the clogging rate; t is time.

The fitting parameters are shown in Table 5. As the particle size ratio increases, both the clogging rate (λ) and the steady-state void ratio (e) gradually decrease. This can be attributed to the fact that the pore structure of the filter has a weaker interception capability for smaller particles, which reduces the clogging efficiency within the filter. Consequently, more fine particles penetrate into the filter, leading to a decrease in the void ratio.

Table 5 Fitting parameters for the variation of the void ratio with time.

By projecting the surface of the void ratio in Fig. 4 along the time axis, the trend of the void ratio with respect to depth is obtained. This is illustrated at the time when silt plugging has stabilized (i.e., t = 4s), as shown in Fig. 6. Overall, when the particle size ratio is relatively small (Ra ≤ 10.0), the void ratio of the filter exhibits significant anisotropy along the depth direction. The reduction in void ratio becomes more pronounced near the top of the filter, with most particles accumulating within a 4 cm range from the filter surface. For Ra = 10.0, the pore ratios of the first and tenth layers are 0.485 and 0.664, respectively, indicating a decrease of 27.3% in the first layer and only 0.3% in the tenth layer. The decrease in void ratio at the same position is significantly influenced by the particle size ratio Ra; for example, the void ratio decrease in the first layer when Ra = 10.0 is 18.6% larger than when Ra = 3.3. In the case of internal penetration with Ra = 16.7, the variability of the void ratios in each layer during infiltration stabilization is smaller, and the siltation is more concentrated in the middle and lower layers, with no clear pattern. Additionally, although the siltation pattern for Ra = 10.0 is classified as internal penetration, its De/det value is close to the critical threshold, meaning this condition exhibits characteristics of internal siltation. In this case, most of the particles are trapped, with only a few particles escaping through the soil layer, resulting in an evolution pattern for its void ratio that still follows an exponential function distribution.

Fig. 6
figure 6

Change in void ratio along depth at infiltration stabilization.

When Ra ≤ 10.0, the variation of the void ratio with depth for different particle size ratios is described by the following equation:

$$e(z) = e_{0} + \varepsilon e^{ - \beta z}$$
(11)

where ε is the fitting coefficient; β is the interception efficiency; z is the infiltration depth.

The fitting parameters are presented in Table 6. As the particle size ratio Ra increases, the interception efficiency β progressively decreases. This occurs because the filter’s efficiency in intercepting finer particles diminishes, allowing more fine particles to move deeper into the filter, which leads to a greater decrease in the void ratio of the deeper layers as Ra increases. Furthermore, when Ra ≤ 10.0, finer particles fill the internal pore spaces of the filter more thoroughly, which causes a more pronounced decrease in the void ratio of the first layer as the particle size ratio increases.

Table 6 Fitting coefficient of void ratio along depth direction.

Effect of muddy water particle size ratio on the hydraulic conductivity of the filters

The following analysis is based on the Kozeny-Carman and the Darcy hydraulic conductivity formula:

$$k_{k} = \frac{{c_{2} \rho_{{\text{w}}} g}}{{s^{2} \mu }}\frac{{e^{3} }}{1 + e}$$
(12)
$$k_{d} = \frac{{\beta d^{2} \gamma_{{\text{w}}} }}{\lambda \mu }\frac{{e^{2} }}{1 + e}$$
(13)

where, k is the hydraulic conductivity (m/s); c2 is the coefficient related to the shape of the particles and the actual flow direction of water, and Carman suggests to take 0.2; ρw is the density of water (kg/m3); g is the acceleration of gravity (m/s2); s is the specific surface area of the particles (m-1); μ is the coefficient of the kinetic viscosity of the water (Pa·s); e is the void ratio; β is the volume coefficient of the soil particles (β = π/6 for the sphere); d is the equivalent diameter of the particles (m); λ is the influence coefficient of the neighbouring particles (3π for the sphere).

The numerical results allow for the convenient calculation of particle specific surface area and equivalent diameter, with the hydraulic conductivity computed using Eqs. (12) and (13) under the aforementioned conditions, as shown in Figs. 7, 8, 9 and 10. Figures 7 and 8 show the variations over time of the Kozeny-Carman hydraulic conductivity and the Darcy hydraulic conductivity for the entire filter under different particle size ratios. As the particle size ratio (Ra) increases, the fine particles entering the filter cause a decrease in porosity, while simultaneously increasing the specific surface area and reducing the equivalent diameter of the filter, leading to a gradual reduction in permeability. For the case of Ra = 10.0, the initial Kozeny-Carman hydraulic conductivity is 0.24 m/s, which decreases to 0.15 m/s at stability, a reduction of 38%. The initial Darcy hydraulic conductivity is 3.63 m/s, which decreases to 3.20 m/s at stability, a reduction of 12%.

Fig. 7
figure 7

Kozeny-Carman equation as a function of time.

Fig. 8
figure 8

Darcy’s law as a function of time.

Fig. 9
figure 9

Variation of Kozeny-Carman equation along depth.

Fig. 10
figure 10

Variation of Darcy’s law along depth.

Figures 9 and 10 show the trend of the hydraulic conductivity along the depth direction. The dashed line labeled "E" in the legend represents the change in specific surface area (s) and equivalent particle size (d) during the seepage process, without accounting for changes in these parameters. This serves to highlight the difference in the hydraulic conductivity when considering muddy water seepage versus clear water seepage. From the graphs, it can be observed that when the particle size is relatively small (Ra = 3.3, 5.0), the number of fine particles entering the filter is low, and the changes in specific surface area and equivalent particle size are minimal. As a result, the impact on the hydraulic conductivity is relatively small, with the pore size ratio playing a dominant role. However, when the particle size is larger (Ra = 6.25, 8.3, 10.0), the infiltration of a larger number of fine particles significantly impacts the hydraulic conductivity. For instance, when Ra = 10.0, the Kozeny-Carman hydraulic conductivity in the first layer was 0.026 m/s without considering the changes in specific surface area and equivalent particle size. When these changes were accounted for, the hydraulic conductivity increased to 0.105 m/s, which is about four times higher. Similarly, the Darcy hydraulic conductivity considering the change in equivalent particle size was 2.16 m/s, compared to 1.68 m/s when this factor was ignored, representing a 28.6% increase. Therefore, in seepage calculations, when water contains a significant amount of fine particles much smaller than the filter’s particle diameter, the permeability damage to the filter caused by these fine particles in muddy water must be considered. Additionally, the Kozeny-Carman equation relies on pore geometric parameters (e.g., specific surface area). When fine particles fill pores and increase the specific surface area, this equation can more sensitively capture permeability changes, making it suitable for describing clogging processes with significant pore-scale variations (e.g., internal bridging or particle filling). In contrast, Darcy’s law is better suited for steady-state seepage in structurally stable conditions, where permeability calculations depend on the equivalent particle diameter derived from geometric averaging—an approach less advantageous for clogging processes dominated by coarse particles in filters. Thus, the Kozeny-Carman equation proves more appropriate for dynamic pore-scale analysis during muddy water infiltration.

Effect of muddy water particle size ratio on the dry density of the filter

The dry density of each layer can be calculated by measuring the mass of the sample layers in the DEM software. Figure 11 illustrates the time-dependent change in the dry density of the filter. As the particle size ratio increases, more fine particles enter the filter, resulting in a gradual increase in dry density. As clogging enters a stable phase, the changes in dry density become more stabilized. Figure 12 illustrates the relationship between the dry density of the filter and depth under stable muddy water seepage conditions. As the depth increases towards the bottom, the increment in dry density becomes smaller. Most of the fine particle accumulation occurs in the middle and upper sections of the filter, with the largest increase in dry density occurring in the L1 layer. Compared to the initial dry density of the sample, the dry density increases by 5.44%, 7.38%, 9.26%, 10.66%, and 12.17% for particle size ratios Ra = 3.3–10.0, respectively. The increase in dry density of the filter indicates an increase in its compaction, caused by the clogging of fine particles during the infiltration of muddy water.

Fig. 11
figure 11

Dry density as a function of time.

Fig. 12
figure 12

Dry density variation along depth.

Influence of hydraulic pressure difference on changes in the pore structure of the filter

Figure 13 presents the temporal evolution of the total void ratio in the filter under particle size ratio Ra = 10.0 with varying top water pressures (2, 5, and 10 kPa). During the initial infiltration stage, the total void ratio showed negligible differences across pressure conditions. Over time, the void ratio gradually stabilized, with marginally greater reductions observed at higher pressures. At clogging stability, the final void ratios reached 0.646, 0.645, and 0.644 for 2, 5, and 10 kPa, respectively.

Fig. 13
figure 13

Variation of total void ratio with time at different pressures.

Figure 14 displays the depth-dependent void ratio distribution under identical conditions. Increased top water pressure caused slight void ratio reductions across all layers. For instance, Layer L1 exhibited a 2.49% greater void ratio decrease at 10 kPa compared to 2 kPa. This indicates that with increasing fluid pressure, a small amount of fine particles were pushed deeper into the filter.

Fig. 14
figure 14

Variation of void ratio with depth at different pressures.

However, while elevated fluid pressure moderately influenced particle penetration depth, its impact on the ultimate clogging pattern remained limited34. This behavior stems from the geometric dominance (e.g., De/det) in particle interception, where hydrodynamic factors (e.g., pressure-driven flow velocity) play a secondary role. Fine particles must overcome geometric constraints (e.g., extrusion deformation or rotational resistance) at pore throats to penetrate the filter. However, the existing hydraulic potential energy proves insufficient to provide particles with adequate energy for overcoming these constraints.

Influence of muddy water particle concentration on the silting effect of the filter

Influence of muddy water particle concentration on the void ratio of the filter

To investigate the influence of fine particle concentration in muddy water on filter void ratio during seepage, four mass flow rates of fine particles (Cc = 5, 10, 25, and 50 g/s) were controlled to represent different sediment concentrations. The variation in the void ratio of the filter with time and depth at different concentrations is shown in Fig. 15. It can be observed that under different infiltration concentrations, the filter’s void ratio exhibits a quadrilateral surface distribution in the spatiotemporal ___domain. As the concentration of fine particles in muddy water increases, the void ratio reaches clogging stability at an earlier time, undergoes a more significant reduction in magnitude, and develops steeper gradients along the depth direction.

Fig. 15
figure 15

Spatial and temporal distribution of void ratios at different concentrations.

Figure 16 shows the relationship between the total pore ratio and time, revealing that concentration significantly influences the pore ratio. As the mass of fine particles flowing into the system increases within a unit time, the rate of pore ratio decrease accelerates in the initial stage, and the system reaches the stability stage earlier. This indicates that higher concentrations of infiltration lead to an earlier establishment of the siltation stability stage. At the same time, higher concentrations of fine particles result in a smaller decrease in the pore ratio after siltation stabilizes, suggesting fewer fine particles enter the filter. Eq. (10) was employed to characterize this process, with fitting parameters summarized in Table 7. The results reveal that increasing fine particle concentration in muddy water enhances the clogging rate λ while elevating the stabilized void ratio e. For instance, at Cc = 5 g/s, the clogging rate λ is 9.3 times greater than at Cc = 50 g/s, while the void ratio reduction is 2.4 times more pronounced. This concentration-dependent behavior directly governs particle retention efficiency and consequently determines the final void ratio. Mechanistically, high concentrations (e.g., Cc = 50 g/s) promote the formation of bridging clogging at the top of the filter. This mechanism restricts deeper particle migration (analogous to a "self-healing" effect), thereby preventing further void ratio reduction.

Fig. 16
figure 16

Variation of total void ratio with time at different fine particle concentrations.

Table 7 Fitting parameters for total void ratio over time.

Figure 17 presents a cloud view of particle positions projected along the negative direction of the Z-axis within the range of -0.5 cm < Z < 0.5 cm for the filter. Red particles represent those near Z = -0.5 cm, indicating the formation of a silt plug at the top of the filter, while blue particles represent those near the Z = 0.5 cm cross-section. Green particles indicate those that have formed a shelf structure within the pore structure of the filter. As the concentration of fine particles increases, certain areas transition through three stages: siltation, semi-siltation, and hollowing. The internal siltation of fine particles at the same cross-section gradually decreases, while the hollowing zone begins to expand. The formation of these hollow zones is the primary factor impeding the downward migration of fine particles.

Fig. 17
figure 17

Fine particle clogging at different concentrations.

Figure 18 illustrates the relationship between the total porosity ratio and depth at different fine particle concentrations during siltation stabilization, with the data fitted using Eq. (11). The fitted parameters are shown in Table 8. It can be observed that the fine particle concentration affects the variation of the porosity ratio along depth. At the first layer, the porosity ratios for concentrations of 5 g/s and 50 g/s are 0.518 and 0.552, respectively, with the former showing a decrease of 5.1% compared to the latter. The siltation efficiency (β) shows a gradual increase with the concentration, although its impact is relatively weaker compared to the effect of the particle size ratio (Ra) on β.

Fig. 18
figure 18

Variation of total void ratio along depth direction.

Table 8 Fitting parameters for the variation of total void ratio along the depth direction.

Influence of muddy water particle concentration on the hydraulic conductivity of the filters

The hydraulic conductivity under different concentrations of fine particles were calculated using Eqs. (12) and (13). Figures 19 and 20 present the Kozeny-Carman and Darcy hydraulic conductivity as a function of time for the entire filter, while Figs. 21 and 22 show the corresponding hydraulic conductivity for the first layer. It can be observed that the variation in the hydraulic conductivity over time follows a trend similar to that of the void ratio discussed in Section “The effect of the particle size ratio (Ra) on the porosity of the filter”. For the average hydraulic conductivity of the entire specimen, the results calculated using Eqs. (12) and (13) at different concentrations show only minor differences. For example, at a concentration of 50 g/s, the Kozeny-Carman hydraulic conductivity at final stabilization is 5.1% higher than at 5 g/s, and the Darcy hydraulic conductivity is 2.4% higher. In contrast, the hydraulic conductivity for the first layer exhibit greater variability. At final stabilization, the Kozeny-Carman and Darcy hydraulic conductivity at 50 g/s are 21.6% and 13.2% higher, respectively, compared to those at 5 g/s. Overall, the effect of fine particle concentration on the hydraulic conductivity is smaller than that of the particle size ratio.

Fig. 19
figure 19

Kozeny-Carman hydraulic conductivity as a function of time.

Fig. 20
figure 20

Darcy’s hydraulic conductivity as a function of time.

Fig. 21
figure 21

Kozeny-Carman hydraulic conductivity of layer L1.

Fig. 22
figure 22

Darcy hydraulic conductivity of L1 layer.

Figures 23 and 24 illustrate the relationship between the two hydraulic conductivity, calculated using Eqs. (12) and (13), along the depth direction when the final siltation has stabilized. The dotted line labeled "E" in the legend carries the same meaning as described in Section “Influence of muddy water particle concentration on the void ratio of the filter”. It can be observed that the influence of the fine-scale parameters in Eqs. (12) and (13) on the hydraulic conductivity is less significant than that of the particle size ratio at different concentrations. For instance, the Kozeny-Carman hydraulic conductivity of the L1 layer, when the change in specific surface area is not considered, is 1.98 times larger than when the change is considered at a concentration of 5 g/s. Similarly, the Darcy hydraulic conductivity of the L1 layer, when the change in equivalent diameter is not considered, is 16.8% larger than when it is taken into account.

Fig. 23
figure 23

Variation of Kozeny-Carman hydraulic conductivity along depth.

Fig. 24
figure 24

Variation of Darcy’s hydraulic conductivity along depth.

Effect of muddy water particle concentration on the dry density of the filter

Figure 25 illustrates the relationship between the dry density of the specimen and time at different fine particle concentrations. Similar to the changes in the void ratio, the dry density increases rapidly during the initial phase of silt plugging. At higher concentrations, the growth rate is faster, and the specimen reaches the stabilization stage earlier. Figure 26 presents the distribution of dry density along the depth direction. It can be observed that the dry density of the first layer (L1) undergoes a more significant change compared to the other layers. Specifically, the dry density of the L1 layer increases by 9.76%, 9.26%, 7.82%, and 7.37% for concentrations ranging from 5 g/s to 50 g/s, respectively. When comparing the effect of fine particle concentration on dry density to the effect of particle size ratio, it is evident that the influence of concentration at the same gradient is smaller.

Fig. 25
figure 25

Dry density versus time.

Fig. 26
figure 26

Distribution of dry density along depth.

Analysis of the effect of siltation on the internal pore structure of the filter

Fine particles in muddy water can alter the original pore characteristics when they are silted up within the filter. The three-dimensional pore reconstruction method, utilizing digital core technology, provides an accurate representation of the pore structure inside the soil matrix. In this study, the changes in the pore structure of the filter, both before and after infiltration, are examined by considering different working conditions with varying void ratios as examples.

The DEM-calculated particle model was imported into CT image processing software. Axial slicing along the Z-axis enabled porosity calculation of the filter under different particle size ratios (Ra) during clogging stabilization, as demonstrated in Fig. 27. The horizontal axis slices numbered 0 and 2000 correspond to Z = 10 cm and Z = 0 cm, respectively. It can be observed that the initial porosity of the filter was approximately 0.4, and the internal porosity of the sample was relatively uniform. As fine particles of varying sizes silt up in the filter, the porosity decreases to different extents. Along the Z-axis, as shown in Fig. 27(a), the porosity change is more pronounced near the top of the specimen, indicating that for particle size ratios of Ra = 3.3–10.0, the topmost layers experience the most significant siltation. Furthermore, comparing the changes in porosity across different particle size ratios, it can be seen that a larger particle size ratio leads to a greater reduction in porosity, with the decrease becoming more noticeable as the depth increases, as illustrated in Fig. 27(b) and( c).

Fig. 27
figure 27

Porosity of slices along the depth direction.

The pore network model (PNM) simplifies the pore structure of porous media into a combination of pores and pore throats. In this model, expanded cavities within the pores are defined as pores and are represented by spheres, while narrow channels are considered pore throats and are represented by cylinders48. Figure 28 illustrates the pore network model of the filter in its initial state. The effective pore diameters, calculated using the pore grid model, were then used to generate both the relative frequency distribution and the cumulative frequency distribution graphs, as shown in Fig. 29. It can be observed that as fine particles of varying sizes silt up the filter, the pore structure changes, causing the distribution of the equivalent pore diameter to gradually deviate from the initial curve. Specifically, in Fig. 29, the equivalent pore radius of the filter exhibits two peaks at 1.1 mm and 1.7 mm, which closely match the effective pore diameter De = 2.2 mm calculated using Eq. (9). The results from the PNM calculations provide a comprehensive representation of the distribution of the equivalent pore diameter within the pore space, further validating the accuracy of Eq. (9). When comparing the relative frequency of the equivalent pore radius across different particle size ratios, it is evident that the frequency of larger equivalent pore radii decreases with an increase in particle size ratio, while the frequency of smaller equivalent pore radii increases. This suggests that fine particles in muddy water progressively fill the larger pores within the filter, while also creating many smaller pores. Notably, the effect is more pronounced at higher particle size ratios.

Fig. 28
figure 28

Filter and its pore network model.

Fig. 29
figure 29

Relative Frequency Distribution of PNM Pore Radii.

Summary

This study employs the CFD-DEM fluid–solid coupling method to investigate the siltation process of the filter under different conditions. A quantitative analysis was conducted on the spatiotemporal evolution of physical parameters such as porosity, hydraulic conductivity, and dry density. The following conclusions are drawn:

  1. (1)

    CFD-DEM simulations of muddy water seepage demonstrated strong consistency with filter clogging mode criteria derived from the effective pore diameter formula, while also showing high agreement with pore network model (PNM) computational results. Surface siltation occurs when the particle size ratio (Ra) is less than 2.2, internal penetration occurs when it exceeds 8.8, and internal siltation occurs when the particle size ratio falls between these values.

  2. (2)

    When the particle size ratio (Ra) approaches the internal siltation, both the temporal evolution of void ratio and its spatial variation along the depth under stabilized conditions demonstrate strong conformity with exponential functions.

  3. (3)

    The particle size and concentration of fine particles in the muddy water significantly influence both the clogging process and the final stable state of the filter. A larger particle size ratio leads to a greater influx of muddy water particles into the filter, while a smaller particle size ratio makes it easier for surface siltation to occur. Additionally, a higher concentration of fine particles in the muddy water results in a shorter time for blockage stabilization, making it more prone to quicker siltation. The fluid pressure difference has a minimal effect on the siltation pattern.

  4. (4)

    Changes in the internal structure of the filter during the siltation process significantly impact the hydraulic conductivity. The CFD-DEM numerical method enables convenient calculation of the equivalent average particle size and specific surface area of the specimen, which are difficult to obtain through experimental methods.