Modeling Water-Induced Base Particle Migration in Loaded Granular Filters Using Discrete Element Method

: Results are reported from a series of ﬁltration tests simulated using coupled computational ﬂuid dynamics and the discrete element method (CCFD-DEM) to investigate the factors controlling the mechanism of base particle erosion and their subsequent capture in loaded granular ﬁlters. Apart from geometrical factors such as particle and void sizes, the ﬁlter effectiveness was found to be controlled by the magnitudes of the hydraulic gradients and the effective stresses. The results of numerical simulations revealed that the base soils exhibit signiﬁcant stress reduction that reduces further due to seepage, and the base particles migrate into the ﬁlter, bearing very low effective stresses (i.e., localized piping in base soil). Based on the limit equilibrium of hydraulic and mechanical constraints, a linear hydromechanical model has been proposed that governs the migration and capture of base particles in the ﬁlter (i.e., ﬁlter effectiveness avoiding piping) for cases simulated in this study. Nevertheless, the proposed model agrees closely with the simulation results of this study and those adopted from other published works, thereby showing a reasonable possibility of using the proposed model as a measure of retention capacity of loaded protective ﬁlters.


Introduction
In practice, granular filters are installed to protect base soils such as dam cores and subgrades from erosion. Recent advances in geotechnical practices have generated an interest in applying granular filters in geo-environmental and transportation engineering [1][2][3], where filtration occurs in the presence of static and cyclic mechanical loading, e.g., railways and highway filters. Unlike large upstream heads in hydropower dams, the hydraulic excitation in railway sub-structures mainly stems from the development of pore pressure at the interface between natural and engineered fills. The filters within these structures function under complex stress states, which significantly influence the erosion-capture mechanism of fine particles such as base soil migration, suffusion, and internal erosion [4,5]. Thus far, various researchers have studied these processes and proposed numerous concepts and definitions of governing threshold hydraulic gradients, e.g., the critical hydraulic gradient for piping [6,7], segregation piping [8], complete erosion of base soil through filters [9,10], and suffusion or washout failure in granular soils [11][12][13]. The erosion of base soils may occasionally favor the development of self-filtering layers that helps prevent any further erosion and increases the longevity of filters. However, the strong upward seepage could incur significant structural instabilities in the forms of piping, suffusion, and heave development. Alternatively, the reduction in permeability due to base soil eroding into the overlying filter layers would result in progressive clogging, especially under cyclic loading conditions [14].
In this study, the results of filtration tests simulated using a coupled CFD-DEM model are reported, based on which we propose a semi-empirical model governing the migration of base particles into loaded filters to quantify filter effectiveness in retaining base particles while avoiding piping. Following its development, the proposed model is compared with the existing theories, and the results are presented. The results of numerical simulations from the current study plus those adopted from published DEM studies on the topic are used to validate the proposed model. Notably, this study focuses on the analysis of loaded granular filters subject to one-dimensional upward seepage and the erosion of base particles from the horizontal interface between base and filter layers, e.g., drainage layers in transportation embankments (sub-ballast and subbase) and inverted filters at the downstream of embankment dams. The scales of the laboratory and numerical simulations reported here are not comparable with the practical dimensions of the field problems, and this is always a limitation in the modeling of most geotechnical processes.

Numerical Model for Base Soil Erosion
In this study, the base and filter particles are modeled as spheres using the discrete element method (DEM) wherein their motions are governed by the force-displacement law and Newton's second law of motion in tandem [15]. The Hertz-Mindlin non-linear contact model is adopted to simulate particle-particle and wall-particle interactions, while the fluid flow is governed by Navier-Stokes' (N-S) equations incorporating the effect of particles [16]. For brevity, the basic coupled CFD-DEM modeling principals are outlined in Appendix A.

Model and Simulation Scheme
The simulated systems comprised of a 0.03 m thick layer of gap-graded base soil NB (Figure 1), which is a widely used core material in embankment dams that exhibits greater potential for seepage-induced piping and two uniform filters NF1 and NF2 of thickness of 0.06 m each. Two separate base-filter systems consisting of spherical particles are modeled in a cuboid (0.12 m × 0.12 m × 0.09 m) within five fixed walls ( Figure 2). These thicknesses were selected based on a twofold rationale, namely, (1) keeping the size ratio between the smallest cell dimension and the largest particle size well above 5 to avoid potential boundary effects, and (2) in line with the existing studies to accommodate enough particles. For instance, the base layer thickness is 30 mm, while the largest erodible base particle size is 2.8 mm, which yields a size ratio of 10.9 for this study. Note that the conventional definition of relative density (R d ) does not apply to the DEM assemblies [17]; therefore, the base-filter systems were modeled at a target porosity of 50% to indicate a medium dense and stable assembly for which R d ≈ 50% could be assumed. Prior geometrical assessments using retention criteria of Terzaghi [18] ((D 15 /d 85 ) ≤ 4), NRCS [19] (D 15 ≤ 4 × d 85R ), and Indraratna et al. [20] ((D c35 /d 85SA ) < 1) have been made. Here, D 15 and D c35 represent the particle size at 15% finer by mass and controlling the constriction size at 35% finer by surface area for the filter soil, respectively. d 85 , d 85R , and d 85SA represent base soil particle sizes at 85% finer by mass, regraded base particle size at 85% finer by mass, and base particle size at 85% finer by surface area, respectively. It revealed that both the filters NF1 and NF2 were ineffective in retaining the base soil-NB (potential for ineffectiveness; NF2 > NF1), but they were individually internally stable [21,22]. The base particles were assigned different colors according to their group sizes to visualize their positions during the simulations. Similarly, the z-positions of the eroding particles were monitored to determine the particle infiltration depths and hence the retention capacity of the filters.   Figure 3 shows the initial and final effective stress magnitudes in the base and filter at simulation times t = 0, 0.2, and 0.3 seconds. The base soil exhibited significant stress reduction due to seepage that could eventually neutralize to approximately zero (i.e., piping) for the cases with ≤ 30 kPa. This observation endorsed the previous experimental  Table 1 presents the properties and parameters of all model components including balls, walls, and fluid. The simulation domain consisted of 8 × 8 × 8 fluids grid along the x, y, and z-directions, respectively. Keeping the impermeable sidewalls as the slip boundaries, the hydraulic pressure is applied at the bottom (upstream/inflow boundary), while the top of the model (permeable downstream/outflow boundary) is set at zero pressure. A total of 24 filtration cases were simulated for the base-filter systems NF1-NB and NF2-NB under σ vt = 0, 10, 20, 30, 40 and 50 kPa up to a calculation time of 0.3 s. The porosity variations during simulations have been monitored via five measurement spheres: one installed at 0.015 m, three installed at 0.045 m, and only one at 0.06 m from the bottom (Figure 2b). Nevertheless, the calculation time influences the process of erosion; the current time was deemed sufficient when compared with 0.04 s [23] and 0.3 s [24] from published studies. The target σ vt is obtained by a servo mechanism that controls the boundary wall velocities ( . u w ) using the following algorithm [21]:

Results and Discussion
. .
where G is the gain parameter that evolves the following stress increment each cycle: where N c , k w n , ∆t, and A define the number of wall-particle contacts, their average stiffness, time step duration, and the wall area, respectively. A relaxation factor α sets the stability requirements given by Equations (4) and (5), which calculates G and updates Equation (2) during each cycle:  Figure 3 shows the initial and final effective stress magnitudes in the base and filter at simulation times t = 0, 0.2, and 0.3 s. The base soil exhibited significant stress reduction due to seepage that could eventually neutralize to approximately zero (i.e., piping) for the cases with σ vt ≤ 30 kPa. This observation endorsed the previous experimental observations that a relatively higher magnitude of seepage stress would be required to induce any disturbance to the base soil (e.g., piping) when a greater magnitude of σ vt is applied on the upper filter layer. Notably, the stress distributions in base and filter soils for both samples NF1-NB and NF2-NB almost remained unchanged when the simulation time increased from 0.2 to 0.3 s. This clearly depicts that a simulation time of 0.2 s would be sufficient, beyond which the stress conditions become steady.   Figure 4 shows the contact force fc distributions for specimen NF1-NB when subjected to = 10-50 kPa. The initial contact force, fc, and hence the contact stress distribution within the filter layer is almost uniform for the given , but it is reduced significantly inside the base soil layer, although there too it remained uniform. This substantial reduction during transfer from the filter to the base soil clearly indicated that the magnitudes of stress reduction factor ( = ⁄ ) in Appendix B for the current basefilter systems were less than unity. Similarly, the magnitudes of inside both the base and the filter layers markedly increased with the corresponding increase in . The final fc distributions (at t = 0.2 s) in the filter remained almost the same compared to the base, wherein the seepage stresses markedly reduced the base particle contact stresses.

Results and Discussion
In this study, the hydraulic gradients = 20 and 35 were applied on both specimens NF1-NB and NF2-NB. Notably, these high hydraulic gradients were chosen to ensure that the base particles are fully mobilized. For no external loading, the magnitude of minimum hydraulic gradient for the complete erosion of base soil through an ineffective granular filter could be computed using model of Indraratna and Radampola [10] and the current magnitudes have been kept slightly larger. The erosion ratio was computed to quantify the extent of base soil erosion and hence the filter effectiveness [23]: where , , and , define the number of eroded and original base particles and their respective radii. The filter was divided into six equal layers of 0.01 m each, and the quantities of eroded base fines and their respective values were calculated. Based on the published DEM studies, a filter exhibiting < 2.25% could be considered effective [24]. Figure 5a presents the correlation between and , NF1-NB exhibited relatively smaller compared to NF2-NB at = 0 and their magnitudes decreased significantly with the increase in the magnitudes of . For instance, at = 35, the specimen   Figure 4 shows the contact force f c distributions for specimen NF1-NB when subjected to σ vt = 10-50 kPa. The initial contact force, f c , and hence the contact stress σ c distribution within the filter layer is almost uniform for the given σ vt , but it is reduced significantly inside the base soil layer, although there too it remained uniform. This substantial σ c reduction during transfer from the filter to the base soil clearly indicated that the magnitudes of stress reduction factor (β = σ base /σ f ilter ) in Appendix B for the current base-filter systems were less than unity. Similarly, the magnitudes of σ c inside both the base and the filter layers markedly increased with the corresponding increase in σ vt . The final f c distributions (at t = 0.2 s) in the filter remained almost the same compared to the base, wherein the seepage stresses markedly reduced the base particle contact stresses.
In this study, the hydraulic gradients i a = 20 and 35 were applied on both specimens NF1-NB and NF2-NB. Notably, these high hydraulic gradients were chosen to ensure that the base particles are fully mobilized. For no external loading, the magnitude of minimum hydraulic gradient for the complete erosion of base soil through an ineffective granular filter could be computed using model of Indraratna and Radampola [10] and the current magnitudes have been kept slightly larger. The erosion ratio R e was computed to quantify the extent of base soil erosion and hence the filter effectiveness [23]: where n j , n k , and r j , r k define the number of eroded and original base particles and their respective radii. The filter was divided into six equal layers of 0.01 m each, and the quantities of eroded base fines and their respective R e values were calculated. Based on the published DEM studies, a filter exhibiting R e < 2.25% could be considered effective [24]. Figure 5a presents the correlation between R e and σ vt , NF1-NB exhibited relatively smaller R e compared to NF2-NB at σ vt = 0 and their R e magnitudes decreased significantly with the increase in the magnitudes of σ vt . For instance, at i a = 35, the specimen NF1-NB exhibited reductions in R e from 6.3% (ineffective) at σ vt = 0 to 2.0% (effective) at σ vt = 40 kPa, which reduced further to almost 0.95% (effective) at σ vt = 50 kPa. Similarly, at i a = 20, NF2-NB showed reductions in R e from 7.5% (ineffective) at σ vt = 0 to 1.4% (effective) at σ vt = 50 kPa.
Notably, when sample NF2-NB was subjected to i a = 35, a subtle increase in the magnitude of R e was observed. However, subjecting NF1-NB to a relatively smaller i a = 20 at the same stress levels could substantially reduce the R e values to almost half of that at i a = 35. It clearly depicted that the erosion and hence the effectiveness of a base-filter system would be significantly influenced by the magnitude of applied hydraulic gradient.
er 2021, 13, x FOR PEER REVIEW 6 of 18 NF1-NB exhibited reductions in from 6.3% (ineffective) at = 0 to 2.0% (effective) at = 40 kPa, which reduced further to almost 0.95% (effective) at = 50 kPa. Similarly, at = 20, NF2-NB showed reductions in from 7.5% (ineffective) at = 0 to 1.4% (effective) at = 50 kPa. Notably, when sample NF2-NB was subjected to = 35, a subtle increase in the magnitude of was observed. However, subjecting NF1-NB to a relatively smaller = 20 at the same stress levels could substantially reduce the values to almost half of that at = 35. It clearly depicted that the erosion and hence the effectiveness of a base-filter system would be significantly influenced by the magnitude of applied hydraulic gradient.    Figure 5b,c show the spatial distributions of eroded base fines captured at various depths inside the filters at t = 0.2 s and 0.3 s, respectively. In both cases, the computed R e values were higher near the interface, because the hydraulic gradients were higher at the interface than elsewhere, and R e gradually decreased with the filter depth. Notably for the given i a values, a significant amount of base fines could reach the downstream filter boundary at relatively smaller magnitudes of applied effective stresses, i.e., σ vt ≤ 20 kPa. These fines could have eroded further had either the magnitude of i a increased or that of σ vt reduced. Interestingly, the magnitudes of eroded fines (R e values) for effective basefilter systems were insensitive to the simulation time because the steady-state conditions could be reached until t = 0.2 s. In contrast, the continuous erosion from ineffective basefilter systems at t > 0.2 s resulted in markedly increased R e values, as shown by red lines in Figure 5b,c.
The expressions of hydraulic gradients for a unit displacement of a base particle inside the filter layer and for complete erosion through the entire filter layer read [10]: where γ and γ w represent unit weights of soil and water, respectively, ∅ is the soil's angle of internal friction in degrees, k f c is the particle contact factor, d b is the base particle size, δ z is the unit displacement of base particle, D cm is the mean constriction size, and h f is the filter depth. Figure 6a,b plot the critical hydraulic gradient values for unit displacement and complete erosion of base particles versus those predicted by the above Equation (7), respectively. Not surprisingly, the magnitudes of i c0 computed through Equation (7) agree closely with those recorded during the numerical simulations of this study, and it may be because both models consider the particles to be spheres. Nevertheless, the values of the critical hydraulic gradients would be expectedly higher for actual granular soil particles (i.e., non-spheres). It clearly establishes that the values of i c0 obtained from Equation (7) would be conservative when applied to the actual field conditions. Therefore, this study adopts both Equation (7) to capture the unit displacement and the complete erosion of base particles through the filter layer (i.e., unit displacement approaching h f ), respectively. Furthermore, the additional hydraulic gradient (i cp ) arising due to the normal loading needs to be overcome for once by hydrodynamic forces to dislodge an erodible fine and induce unit displacement (i.e., initiation of erosion) and can be given by (Appendix B): Thus, the expression for the complete erosion of base particles (i.e., both initiation and development) from the filter layer can be given by: Figure 7 presents the base-filter particle distributions for select specimens at the end of simulations i.e., at t = 0.2 s for filter NF2 subjected to 0, 10 kPa, 30 kPa, and 50 kPa. As shown, the magnitude of erosion has been significantly reduced with the increase in the magnitude of normal stress. Similarly, the filter's ability of capturing the eroded fines has been markedly enhanced with the increase in stress magnitude. This clearly depicts that the filter effectiveness could be increased under higher stress magnitudes. Nevertheless, it is fully consistent with our understanding of increasing the stability and effectiveness of inverted filters by constructing berms and providing additional surcharge at downstream of hydraulic structures.
Water 2021, 13, x FOR PEER REVIEW 9 of 18 Figure 6. (a) Plot of hydraulic gradient for unit displacement of base particle observed [10] and that estimated by Equation (7) Figure 6. (a) Plot of hydraulic gradient for unit displacement of base particle i c0 observed [10] and that estimated by Equation (7) versus mean constriction size D cm values, and (b) comparison between experimentally observed and theoretically estimated i cr data for the complete erosion of base soil through filter.
Water 2021, 13, x FOR PEER REVIEW 9 of 18 Figure 6. (a) Plot of hydraulic gradient for unit displacement of base particle observed [10] and that estimated by Equation (7) Table 2 presents the results of numerical data of this study plus those adopted from the published DEM studies for validation. In Figure 8, the magnitudes of maximum applied hydraulic gradients during DEM simulations, i a, max , and those predicted by the proposed model, i ct (Equation (8)), were plotted. The line of equality demarcates a visible boundary between effective and ineffective base-filter systems, whereby the ones plotting on its right were deemed ineffective and vice versa. The assessment results of filter effectiveness shown in Figure 8 have been quantified in Table 2 in the form of effectiveness of a filter in fully protecting the base soil. For brevity, the current model assesses a filter as ineffective should the applied hydraulic gradient be greater than its anticipated i ct value (Equation (8)). Then, this assessment is compared with the actual DEM response observed in this study or that reported in the adopted works. As Table 2 shows, out of 40 simulation results from three independent DEM studies including the current work, only five incorrect assessments (i.e., two conservative assessments including N3 and N10 and three unsafe including N11, N25, and N39) are obtained by the currently proposed theoretical model. This acceptably smaller discrepancy in numerical results and theoretical assessments of this study could be because the proposed model is an explicit solution using a continuum approach, while CCFD-DEM simulations of this study capture micro-mechanical responses at the fluid-particle level. Note: Here n f , i a,max , and i ct define filter porosity, maximum applied hydraulic gradient, and hydraulic gradient governing particle migration estimated from Equation (8), respectively. observed in this study or that reported in the adopted works. As Table 2 shows, out of 40 simulation results from three independent DEM studies including the current work, only five incorrect assessments (i.e., two conservative assessments including N3 and N10 and three unsafe including N11, N25, and N39) are obtained by the currently proposed theoretical model. This acceptably smaller discrepancy in numerical results and theoretical assessments of this study could be because the proposed model is an explicit solution using a continuum approach, while CCFD-DEM simulations of this study capture micromechanical responses at the fluid-particle level.  Table 2.

Theoretical Envelope
A graphical illustration of the proposed theoretical hydro-mechanical envelopes (Equation (8)) governing the migration of base particles in loaded filters is presented in Figure 9. It followed a straight-line relationship with the dimensionless mechanical constraint Ɲ (= ′ ∆ ⁄ ), having a slope equal to the unit (Equation (7)). This definition of slope is more realistic because the proposed envelope for a given base-filter system tends to be unique by the virtue of its physical properties. For an internally stable filter under ′ = 0 (i.e., with no additional surcharge on the filter layer except the self-weight), the can be estimated by the ordinate intercept (OA) of the theoretical envelope in Figure 9. Almost all the inter-particle contacts are mechanically active in an internally stable soil [25]; thus, any value of ′ > 0 can help a filter resist the erosion of particles by virtue of the increased mechanical and frictional constraints. In turn, this requires larger values to trigger the migration of base particles (path BCD), thus imitating Terzaghi's [6] Estimated (Eq. 8) Applied hydraulic gradient, Qing-fu et al. [24] Current study Figure 8. The comparison between observed and estimated i cr values by the proposed model for numerical modeling data in Table 2.

Theoretical Envelope
A graphical illustration of the proposed theoretical hydro-mechanical envelopes (Equation (8)) governing the migration of base particles in loaded filters is presented in Figure 9. It followed a straight-line relationship with the dimensionless mechanical constraint three unsafe including N11, N25, and N39) are obtained by the currently proposed theoretical model. This acceptably smaller discrepancy in numerical results and theoretical assessments of this study could be because the proposed model is an explicit solution using a continuum approach, while CCFD-DEM simulations of this study capture micromechanical responses at the fluid-particle level.  Table 2.

Theoretical Envelope
A graphical illustration of the proposed theoretical hydro-mechanical envelopes (Equation (8)) governing the migration of base particles in loaded filters is presented in Figure 9. It followed a straight-line relationship with the dimensionless mechanical constraint Ɲ (= ′ ∆ ⁄ ), having a slope equal to the unit (Equation (7)). This definition of slope is more realistic because the proposed envelope for a given base-filter system tends to be unique by the virtue of its physical properties. For an internally stable filter under ′ = 0 (i.e., with no additional surcharge on the filter layer except the self-weight), the can be estimated by the ordinate intercept (OA) of the theoretical envelope in Figure 9. Almost all the inter-particle contacts are mechanically active in an internally stable soil [25]; thus, any value of ′ > 0 can help a filter resist the erosion of particles by virtue of the increased mechanical and frictional constraints. In turn, this requires larger values to trigger the migration of base particles (path BCD), thus imitating Terzaghi's [6] Estimated (Eq. 8) Applied hydraulic gradient, Qing-fu et al. [24] Current study (= σ vt /γ w ∆y i ), having a slope equal to the unit i c0 (Equation (7)). This definition of slope is more realistic because the proposed envelope for a given base-filter system tends to be unique by the virtue of its physical properties. For an internally stable filter under σ vt = 0 (i.e., with no additional surcharge on the filter layer except the selfweight), the i c0 can be estimated by the ordinate intercept (OA) of the theoretical envelope in Figure 9. Almost all the inter-particle contacts are mechanically active in an internally stable soil [25]; thus, any value of σ vt > 0 can help a filter resist the erosion of particles by virtue of the increased mechanical and frictional constraints. In turn, this requires larger i cr values to trigger the migration of base particles (path BCD), thus imitating Terzaghi's [6] recommendation of using inverted loaded filters to avoid quicksand conditions at the toe of embankment dams.  Nevertheless, a different response would be expected with internally unstable soils whereby the inter-particle contacts are not mechanically significant, and hence, the load transfer is not efficient. This would not increase the inter-particle friction sufficiently, and Nevertheless, a different response would be expected with internally unstable soils whereby the inter-particle contacts are not mechanically significant, and hence, the load transfer is not efficient. This would not increase the inter-particle friction sufficiently, and hence the theoretical envelope tends not to obey the proposed law. The seepage forces accompany erodible finer fractions from the filter itself, which progressively increases its porosity (thus, the constriction sizes), and then, the base particles require less i cr for erosion through the filter. Hence, a prior assessment of internal stability through an acceptable geometrical criterion is imperative to select internally stable and effective filters.
In Figure 10, the applied i a and the Figure 8. The comparison between observed and estimated values numerical modeling data in Table 2.

Theoretical Envelope
A graphical illustration of the proposed theoretical hyd (Equation (8)) governing the migration of base particles in load Figure 9. It followed a straight-line relationship with the dimen straint Ɲ (= ′ ∆ ⁄ ), having a slope equal to the unit (E tion of slope is more realistic because the proposed envelope for tends to be unique by the virtue of its physical properties. For under ′ = 0 (i.e., with no additional surcharge on the filter laye the can be estimated by the ordinate intercept (OA) of the th ure 9. Almost all the inter-particle contacts are mechanically acti soil [25]; thus, any value of ′ > 0 can help a filter resist the ero of the increased mechanical and frictional constraints. In turn, th ues to trigger the migration of base particles (path BCD), thus values for the numerical simulation results of this study are plotted with the corresponding theoretical envelopes for the base-filter systems NF1-NB and NF2-NB. Apart from a single incorrect assessment NF1-NB-30 (i.e., test N20 in Table 2), the rest of the base-filter systems are correctly plotted in their respective effective and ineffective zones by the proposed model. Nonetheless, based on the geometrical and mechanical factors, the proposed theoretical envelopes governing base particle erosion for a given base-filter system should be unique.  Table  2 for data).

Conclusions
The following conclusions could be drawn from this study. The results of CCFD-DEM numerical modeling of this study plus those adopted from the published literature showed good agreements with both the experimental observations and the model predictions. Unlike the available particle size distribution (PSD)based geometrical criteria for assessing the effectiveness of filters, a performance-based hydro-geo-mechanical approach that considers the level of hydro-mechanical excitation (ia and ′ ) and the geometrical and physical characteristics of filter media (e.g., PSD, Rd, and ∅' etc.) was presented.
No significant erosion occurs from an effective base-filter system at extended simulation times (t = 0.2 s to t = 0.3 s), whereby steady-state conditions prevail due to the formation of self-filtering layers inside the filter. In contrast, the erosion of base particles continues through an ineffective filter should the simulation time increase. Furthermore, the changes in the magnitudes of applied hydraulic gradients markedly affect the magnitude of erosion.
An enhanced alternative method for quantifying filter effectiveness through the factor of safety against base particle erosion instead of traditional "effective or ineffective" bifurcation may be more efficient and cost effective in practice; for instance, designing filters for properly monitored low-risk structures such as drainage layers subject to low   Figure 10. Validation through coupled CFD-DEM simulation results (test series numbers in Table 2 for data).

Conclusions
The following conclusions could be drawn from this study. The results of CCFD-DEM numerical modeling of this study plus those adopted from the published literature showed good agreements with both the experimental observations and the model predictions. Unlike the available particle size distribution (PSD)-based geometrical criteria for assessing the effectiveness of filters, a performance-based hydrogeo-mechanical approach that considers the level of hydro-mechanical excitation (i a and σ vt ) and the geometrical and physical characteristics of filter media (e.g., PSD, R d , and ∅' etc.) was presented.
No significant erosion occurs from an effective base-filter system at extended simulation times (t = 0.2 s to t = 0.3 s), whereby steady-state conditions prevail due to the formation of self-filtering layers inside the filter. In contrast, the erosion of base particles continues through an ineffective filter should the simulation time increase. Furthermore, the changes in the magnitudes of applied hydraulic gradients markedly affect the magnitude of erosion.
An enhanced alternative method for quantifying filter effectiveness through the factor of safety against base particle erosion instead of traditional "effective or ineffective" bifurcation may be more efficient and cost effective in practice; for instance, designing filters for properly monitored low-risk structures such as drainage layers subject to low magnitude hydraulic excitation in slopes, embankments, levees, railway, and highway sub-structures. A novel theoretical model has been proposed to quantify filter effectiveness under static loading condition that could also be presented graphically. Nevertheless, the current simulation results and those adopted from the published literature could sufficiently verify the proposed model, thus enhancing the user confidence in adopting it for various preliminary analyses and assessment purposes.
Needless to mention, the findings from this study only govern internally stable and non-cohesive base-filter systems, but the authors do envisage extending the scope to internally unstable granular filters and cohesive base soils.
where , = volume fractions and , = particle radii of filter and base, respectively, and f , f = normal and tangential vector components of contact forces. Figure A1. Vertical effective stresses on soil volume subject to the upward seepage flow.

Appendix B. Equivalent Hydraulic Gradient for the Mechanical Constraint
For a granular medium of arbitrary depth ∆ (from to ℎ ) under an effective stress ′ and subject to upward flow due to ( Figure A1), the bottom effective stress ′ reads: where; (A14 For a granular medium of arbitrary depth ∆y i (from D f to h f ) under an effective stress σ t and subject to upward flow due to i a (Figure A1), the bottom effective stress σ b reads: where; At the onset of erosion into a filter, the stress σ b on the eroding base particles should be zero: Equation (A15) governs the uniform soils where all particles participate equally in transferring σ vt . However, with a non-uniform soil, the share of coarse particles due to their larger surface areas is more than with fine particles. The contribution of fines declines further as the seepage forces increase, and it becomes zero when they begin to erode (quick condition). Thus, a stress reduction factor β (=i c0 /i cr,0 )) is introduced with i c0 denoting the actual hydraulic gradient governing the particle erosion at σ b = 0, and i cr,0 (= γ /γ w ) is the critical hydraulic gradient for quick condition in base soil [6]. The factor β has physical meanings, whereby the particles will not erode individually until β→1 (i c0 ≈ i cr,0 ). However, for the ineffective base-filter systems, the base particles may start to erode even when β < 1 and the filter particles are still intact, i.e., at σ base ≈ 0, it is likely that σ f ilter > 0.
i cp = β i cr,0 + σ vt γ w .∆y i = i c0 i cr,0 i cr,0 + σ vt γ w .∆y i (A16) At the inception of erosion in internally stable base soils, the term i cr,0 ≈ 1 (fines in their loosest state). Similarly, the term i cr,0 in Equation (A16) represents the contribution of buoyant weight that could be relevant for full-scale problems with non-negligible filter thickness and is negligible in laboratory tests. Therefore, ignoring the first term in Equation (A17) when compared with the magnitude of the second term for a 150 mm thick sub-ballast layer or a 500 mm thick filter layers under more than 30 kPa [28] and 100 kPa surcharges, respectively: i cp = i c0 σ vt /γ w .∆y i (A17)