Sample Preparation

The Yellow River Stone Forest is located upstream of the Heishanxia section along the Yellow River mainstem. Upon completion of the downstream Daliushu High Dam, areas below the 1380 m elevation within the Stone Forest section will be inundated, forming part of the reservoir (Fig. 1b)). The slopes are extensively covered with Loose Gravel Soil (Fig. 1c)), primarily consisting of subangular gravels, sand particles, and clayey silts with well-developed pore structures. Post-impoundment seepage flow pressures can trigger engineering hazards such as bank instability and slope failures. This study targets this region, collecting Loose Gravel Soil samples from five locations within the inundation zone shown in Fig. 1b).

The gravelly sandstone soil displays a loosely compacted granular structure. Prior to sampling, surface aeolian deposits were carefully removed. Undisturbed soil specimens were extracted at a depth of 10 cm using a ring cutter. The specimens were used for porosity and natural density measurements. Then, bulk samples were collected from the corresponding strata for experimental characterization and numerical modeling. The constant-head permeability test was systematically conducted on sandy gravel soils, and the measured hydraulic conductivity values were compared to numerical simulation results to validate the model. The fundamental geotechnical properties of the specimens, determined in accordance with Chinese National Standards GB/T 14684-2022 (Test methods for building sand) and GB/T 34533-2023 (Technical specification for soil tests), are summarized in Table 1.

Particle size distribution analysis was conducted using the sieving method. Previous studies [37, 38] indicate that boundary effects can be avoided when the specimen diameter exceeds five times the maximum particle size. Due to constraints of the experimental apparatus (5-10 cm diameter), scaled specimens were prepared using the equivalent substitution method [39, 40], whereby particles coarser than 10 mm were removed. The resulting gradation curve (Fig. 2) yields a uniformity coefficient C_u = 15.3 and a curvature coefficient C_c = 0.73. This indicates that the overall gradation of the gravelly sandstone crushed stone soil sample is good, with a wide distribution of various grain groups. However, the soil contains a higher proportion of fine particles, and the size variation among these fine particles is relatively small.

Geometric Modeling

The sandy conglomerate gravel soil exhibits a loose structure. For simulations, the discrete element method (DEM) more accurately represents the structural and positional relationships among actual particles, enabling detailed investigation of the pore structure. This study utilized PFC 5.0 (Itasca Consulting Group; http:// www.itascacg.com) to model the particles of sandy conglomerate gravel soil. Following the methodology outlined in the literature [28], the red segment of the continuous gradation range in Fig. 2, after excluding the particle sizes at both ends of the distribution (DEM-CFD), is selected as the modeling basis. The particle diameters smaller than d₁₀ (0.14 mm), d₂₀ (0.20 mm), and d₃₀ (0.28 mm) within this range are generalized as spherical particles, while the larger particles are modeled using aggregate photography techniques. The models containing different proportions of spherical particles are labeled as B-10, B-20, and B-30, respectively.

Fig. 1.

The study area and sampling points. a) location; b) landform of the study area; c) sampling area, where the red dotted line is the Loose Gravel Soil.

Fig. 3. The aggregate imaging protocol involved suspending Loose Gravel Soil samples from the Yellow River Stone Forest area using transparent filaments connected to a low-power adjustable motor. Controlled rotational motion facilitated systematic photographic acquisition across multiple elevation planes for analytical processing [28].

In 3DF Zephyr (3DFLOW; http://www.3dflow.net), multi-view stereo matching technology is applied to analyze the photographs captured from different angles. Identical feature points in the images are identified and matched to calculate their positions in three-dimensional space. The depth value of each pixel, representing the distance between each pixel and the imaging device, is determined through stereo matching. Combined with the focal length and optical center coordinates, the corresponding three-dimensional coordinates are derived using Equations (1) and (2). This set of coordinates forms a sparse point cloud.

Fig. 2.

Grading distribution of Loose Gravel Soil.

Fig. 3.

Geometric modeling process.

where Z is the depth value; f_x and f_y are the focal lengths of the imaging device; c_x and c_y are the coordinates of the optical center.

Cloth Simulation Filter (CSF) [41] algorithm optimization enhanced point cloud fidelity by outlier removal and noise reduction. Feature points with close depth values can be converted into three-dimensional coordinates to construct a dense point cloud, followed by Poisson surface reconstruction generating watertight mesh models.

The import of 72 clump templates into PFC 5.0 required internal vectorization via the Bubble Pack algorithm [42], with surface topography regulated through sphere centroid spacing adjustments using the “distance” command [43]. Virtual specimen generation replicated actual gradation profiles, with spherical fines occupying interstitial voids between realistic aggregates. Initial particle overlaps were removed through servo-controlled pre-compaction [44], resulting in a target porosity of 30% in the final virtual specimens. Due to the high porosity exhibited by the Loose Gravel Soil sample, variations in aggregate and pore structures are not highly sensitive as long as the total porosity of the model remains comparable [28]. Under these conditions, the virtual specimen is constructed as a cube with a 2 cm edge length.

Since the model relies only on mechanical contact among particles and between particles and boundaries, the mechanical properties of the system have not been evaluated. Mechanical interactions were modeled exclusively using linear contact formulations (Equation (3)) [45, 46].

where F_n is the normal contact force, K_n is the normal contact stiffness coefficient, and R_n is the normal overlap amount. The calculation of shear force and tangential overlap is shown in Equation (4).

where ΔF_s is the shear force; K_s is the tangential contact stiffness coefficient; ΔR_s is the tangential overlap amount.

Pore Structure 3D Reconstruction

Direct export of the three-dimensional model data of virtual specimens from PFC 5.0 can result in excessive memory consumption. A sequential sectioning method [47] was applied to the B-30, B-20, and B-10 virtual specimens to mitigate this. Model slices were generated along the Z-axis at 0.016 mm intervals. After trimming the edges, 1250 high-resolution slices (1250×1250 pixels) were processed to reconstruct the 3D pore architecture using digital image analysis.

The AVIZO platform (FEI Visualization Sciences Group; https://www.fei.com/software/amira-avizo) performed threshold segmentation using Otsu’s method to determine pore-grain interfaces. Interpolated slices were then stacked to reconstruct the 3D geometry. The Separate Objects module was employed to isolate regions with inadequate resolution or ambiguous segmentation. This preprocessing step ensured precision in the PNM construction and in subsequent fluid flow simulations through pores and throats, as illustrated in Fig. 4.

Fig. 4.

Pore 3D reconstruction process.

Fig. 5.

REV extraction schematic.

REV Extraction

The Representative Elementary Volume (REV) is defined as the minimum volumetric unit that statistically captures the microstructural characteristics of porous media when analyzing their macroscopic physical properties [32]. Conventional 2×2×2 mm³ models used in pore-scale seepage flow simulations are inadequate for specimens containing heterogeneously distributed pores >2 mm [48-50]. This necessitates careful REV selection to balance physical representativeness and computational feasibility. The B-10 specimen serves as the primary case for subsequent REV determination.

REV determination for geomaterials requires consideration of both geometric configuration and spatial anisotropy of pore structures. Following established methodologies, a nine-point sampling scheme (Fig. 5) is implemented to characterize virtual pore structures [32, 51, 52]. The protocol sequentially selects eight vertices (A–H) and the centroid (O) of the cubic structure. Starting with reference points (H and O), spatial increments (L, 2L, ...) expand sub-cube dimensions until encompassing the entire structure, generating multiple scaled sub-volumes per reference point. A total of 180 sub-cubes are created for comprehensive parameter correlation analysis.

Following Chu et al.’s [53] approximation criterion, maximum pore diameter emerges as the critical parameter for porosity-based REV determination, governing key transport properties, including permeability and thermal/electrical conductivity. Therefore, the REV dimension analysis focuses on porosity characterization.

Porosity-size relationships derived from nine-point sampled sub-cubes (Fig. 6) indicate optimal REV edge lengths of 400-600 voxels for the 3D pore architecture. Quantitative REV determination is achieved through porosity coefficient of variation (C_V) analysis (Fig. 7), eliminating subjective interpretation. Adopting Min’s [54] 5% C_V threshold, 540 voxels (8.64×8.64×8.64 mm³) is established as the optimal REV dimension for the 3D Loose Gravel Soil pore model.

Seepage Flow Simulation

PNM Single-Phase Flow Calculation

The computation of microscopic pore flow in porous media correlates with intra-pore flow dynamics. This study conducts pore-scale flow simulations using the XLabSuite Extension module integrated in Avizo. The fluid flows from high-pressure regions to low-pressure regions, where the pressure gradient demonstrates a linear correlation with flow quantity, as illustrated in Fig. 8. The flow rate q_ij between pores i and j is determined by Equation (5).

(5)

qij=gij(Pi−Pj)

where g_ij is the hydraulic conductivity of the throat between pore i and pore j, calculated from Equation (6); P_i – P_j is the pressure difference between the two pores.

(6)

gij=π8μrij4lij

(7)

Q=∑(Pi−Pj)gij

where μ is the dynamic viscosity of the fluid; r is the throat radius; l is the throat length. The total flow rate Q can be calculated, as shown in Equation (7), by summing the flow rates of all interconnected pores.

The pore spaces in macropore-containing network models can be resolved by establishing linear equation systems through Equations (5) and (7). Permeability, which quantifies the fluid transmission capacity of soil pores, is typically computed via Darcy’s law.

(8)

K=QμLΔPA

Fig. 6.

Variations in the relationship between porosity and REV scale for Loose Gravel Soil.

Fig. 7.

Porosity coefficient of variation versus subsample size for Loose Gravel Soil.

Fig. 8.

Pore and throat connection unit model.

In accordance with Darcy’s law, Equation (8) enables the determination of the permeability of the virtual specimen REV. The accuracy of the model is verified by comparing it to experimental permeability measurements. In this context, ΔP represents the imposed pressure gradient at the boundaries; A is the cross-sectional area; L is the characteristic flow length.

Boundary condition initialization: Slip phenomena at fluid-solid interfaces are neglected, as they rarely occur in pore structures. The watertight geometric model establishes inlet and outlet pressures along the seepage flow direction, ensuring a quasi-static pore pressure distribution and facilitating fluid migration through pore networks. With fluid properties maintained at 20℃ and a dynamic viscosity of 0.001 Pa·s, the boundary pressures are specified as 0.15 kPa (inlet) and 0 kPa (outlet).

Navier-Stokes Flow Calculation

In pore-scale seepage flow scenarios, the fluid exhibits a low flow velocity and minimal Reynolds numbers, justifying the omission of inertial forces. The Creeping Flow interface in COMSOL (COMSOL Inc.; http://www.comsol.com) applies the reduced Navier-Stokes equations to characterize mass-momentum dynamics in fluid motion, corresponding to Stokes flow [55]. The mathematical formulation of Stokes flow is presented in Equation (9).

(9)

ρ(∂u∂t+u·Vu)⏟1=−∇p⏟2+∇·(μ(∇u+(∇u)T)−23μ(∇u)I)⏟3+F⏟4

where u is fluid velocity; p is fluid pressure; ρ is fluid density; ∇ is the Laplace operator; μ is fluid dynamic viscosity; I is the identity matrix. The terms 1, 2, 3, and 4 correspond to inertial forces, pressure, viscous forces, and external forces acting on the fluid, respectively.

Defined as the ratio of inertial forces 1 to viscous forces 3, the Reynolds number remains exceptionally small in pore-scale flows, indicating laminar flow regimes and quasi-incompressible behavior within microscopic pores. The divergence-free constraint for incompressible fluids is expressed in Equation (10).

(10)

∇·u=0

(11)

0=−∇p+∇·(μ(∇u+(∇u)T)

At ultralow Reynolds numbers, inertial forces become negligible in comparison to viscous forces 3, permitting their exclusion from computational models. In the absence of external forces (external force 4 = 0) and by substituting Equation (10), the reduced form of the Stokes equations is obtained as Equation (11). When the Reynolds number is extremely small, inertial forces are negligible relative to viscous forces and can be omitted from the calculations. In addition, as no external force is applied, the force acting on the fluid is zero. Substituting Equation (10) into the governing equation yields Equation (11), the simplified Stokes equation.

Based on the temporal evolution of the hydraulic gradient during reservoir impoundment, the erosion model validation framework established by Zhang et al. [56, 57] is implemented to examine depth-dependent flow velocity responses to hydraulic gradient variations in Loose Gravel Soil.

The test was conducted using the apparatus illustrated in Fig. 9. The soil column had a height of 20 cm, with piezometers BWK-1, BWK-2, and BWK-3 installed at depths of 15 cm, 10 cm, and 5 cm from the base, respectively. The upper end of the column was submerged to maintain a constant hydraulic head. Pore water pressure measurements were recorded at 30-second intervals, while the hydraulic gradient was incremented by 0.4 every 5 minutes through controlled adjustments of the free surface level in the water tank. Electrical signals acquired by the piezometers were processed via a data logger and transmitted to a computational unit to record temporal pressure variations at each monitoring point. Seepage flow velocities at each cross-section were subsequently calculated using Equation (12).

(12)

q=−kdpdx

where q is the water flow velocity; k is the permeability measured in the experiment; dp is the pore water pressure; dx is the height of the permeable soil column.

In Stokes flow simulations, initial inlet and outlet boundary pressures were defined using pore water pressures measured per minute at each piezometer depth, with the outlet configured as a free surface. Since pore water pressure dynamically varied with hydraulic gradient changes during the test, the inlet pressure in Stokes flow calculations was proportionally adjusted to match hydraulic gradient variations. Quantitative comparison between seepage test results and Stokes flow simulation outcomes further validated the model’s accuracy.

Fig. 9.

Test apparatus for seepage.

Verification of the Numerical Model

Connectivity of Pore Structure

The topology of the pore model describes the connectivity relationships between pores and throats. The spatial topology of pore networks is quantitatively characterized through pore coordination number distributions.

REV exhibits comparable throat parameters across spatial positions, where the primary pores show equivalent radii of 0.50-0.61 mm and surface areas between 9-13 mm². Distinct maximum coordination numbers are observed across REVs with varying spherical particle contents: B-30 specimens demonstrate peak values of 24, 22, 23, 22, 21, 22, 26, 27, and 28, in contrast to B-10’s 30, 27, 29, 30, 28, 29, 27, 30, and 30.

Fig. 11 exhibits that the coordination number distributions of REVs are compared across B-30, B-20, and B-10 configurations at various spatial locations. The results indicate that while spatial heterogeneity exists in coordination number distributions, the cumulative probability curves exhibit similar patterns under constant spherical particle content. For B-30 specimens, 80% of coordination numbers fall within the 3-10 range, while values above 20 constitute less than 2%, indicating attributes of a normal distribution.

A comparison of pore coordination number evolution across B-10, B-20, and B-30 demonstrates that reduced spherical particle content leads to a sequential decrease in REV pore counts (approximately 1900, 1400, and 1000). The modeled coordination number frequency declines progressively, with pores possessing coordination numbers ranging from 1-10 comprising 77%, 82%, and 87% of the total populations. Among these, B-10 exhibits enhanced pore connectivity. The underlying mechanism is attributed to the spatial architecture of spherical particles. Fig. 12 illustrates that higher concentrations of spherical particles increase the likelihood of particle-to-particle interactions. Under identical spatial constraints, the regular geometry of spherical particles results in homogeneous pore structures during packing, characterized by morphological and structural consistency. In contrast, the geometric variability of irregular particles leads to complex and disordered pore architectures during packing, where such structural irregularities favor the formation of increased porosity, particularly within the micropore domain.

Fig. 10.

Comparison of the numerical results and the test results.

Fig. 11.

Coordination number distribution. (a-i) is REV taken at points A, B, C, D, E, F, G, H, and O.

Fig. 12.

Comparison of pore structure between real aggregate model and spherical particle model. (1-7) is the contact between particles; (a-d) is the pores between particles.

Pore Diameter Analysis

Pore diameters are obtained by measuring the circumscribed sphere diameters within the three-dimensional reconstructed pore network model. All pore diameters within each REV are statistically analyzed to assess pore size distribution patterns. Fig. 13 presents the frequency distribution histogram of pore diameters.

The frequency curves reveal that pore diameters in B-10, B-20, and B-30 specimens predominantly range from 0.10-1.0 mm, displaying characteristics of a normal distribution. The most frequent pore diameters fall within the 0.26-0.7 mm interval, with mean values ranging from 0.42-0.55 mm and standard deviations between 0.23-0.38.

Current literature classifies pores into four categories based on size: micropores (<10 nm), small pores (10-100 nm), mesopores (100-1000 nm), and macropores (>1000 nm) [61]. Based on this classification, the reconstructed pores in Loose Gravel Soil from this study primarily fall within the macropore range.

Characteristics of Pore Morphology Distribution

Loose Gravel Soil exhibits considerable pore arrangement diversity and randomness due to irregular aggregate particles and an inherently loose structure, with pores displaying highly irregular morphologies [62]. Although simple geometric shapes cannot fully represent pore complexity, shape factors are utilized to quantify irregularity and better approximate actual pore geometries [63]. The formula for calculating the shape factor is provided in Equation (13).

where V is the pore volume; L is the pore length; A is the pore surface area; C is the pore shape factor.

A larger pore shape factor indicates a more regular pore shape, while a smaller shape factor corresponds to a flatter pore shape. The shape factor of a cube is 0.03, that of a sphere is 0.05, and the shape factor of a tetrahedron, which varies in form, ranges from 0 to 0.07.

Fig. 13.

Distribution curve of pore diameter frequency. (a-i) refers to REV taken at points A, B, C, D, E, F, G, H, O.

Fig. 14.

Distribution curve of pore shape factor frequency. (a-i) is REV taken at points A, B, C, D, E, F, G, H, and O.

Using the Label Analysis module in AVIZO software, pore volume, length (derived from circumscribed sphere diameters), and surface area (calculated through voxel quantification) were obtained for each pore. These parameters were then applied in Equation (13) to determine individual pore shape factors. Statistical analysis of shape factors reveals characteristics of pore regularity, as shown in the frequency distribution curves in Fig. 14.

Analysis of shape factor distributions reveals unimodal, near-symmetric patterns across B-10, B-20, and B-30 REVs. Most values fall within the 0.01-0.05 range (normal distribution), with mean values between 0.025-0.035. This indicates that cubic or rectangular prism geometries dominate, while spherical or tetrahedral pores constitute minor components. The structural instability of cubic or rectangular shapes increases their susceptibility to hydraulic-induced damage during seepage flow processes.

Seepage Characteristics

PNM-based seepage flow simulations were conducted for models B-10, B-20, and B-30 in three orthogonal directions (X, Y, Z). Fig. 15 and Fig. 16 illustrate the pressure field configuration and streamline distribution of the velocity field during fluid migration. As the REV pore architecture and connectivity display limited spatial heterogeneity within individual models, the analysis focused solely on REVs located at position O.

Fig. 15 illustrates progressive pressure attenuation along flow paths. Enlarged pore throat dimensions with increased branching are associated with reduced pressure gradients and smoother transitions, whereas geometrically constrained pores retain elevated pressure levels accompanied by intensified pressure oscillations. The pressure field exhibits consistent evolutionary patterns across flow directions, indicating similar behavioral characteristics between transverse and longitudinal flows. The velocity field depicted in Fig. 16 reveals significant flow heterogeneity within the pore structure. Central regions exhibit more continuous and densely clustered streamlines, indicating enhanced connectivity. Both X- and Y-direction flows exceed Z-direction flow regarding global connectivity metrics. In addition, across all spherical particle content ratios, REV simulations indicate localized velocity amplifications at throat constrictions, along with concentrated streamline distributions in all flow directions. Increased spherical particle fractions facilitate streamlined flow trajectories with greater stability, improving the global permeability characteristics of the system. In contrast, lower spherical particle concentrations (B-10) produce geometrically complex pore networks due to enhanced aggregate representation, resulting in tortuous flow paths with proliferated micro-channels that substantially reduce permeability efficiency.

Fig. 15.

Pressure field distribution of REVs under different flow directions. (a-c) is the pressure field of REV at point O of B-30 in the X, Y, and Z directions; (d-f) is the pressure field of REV at point O of B-20 in the X, Y, and Z directions; (g-i) is the pressure field of REV at point O of B-10 in the X, Y, and Z directions.

Fig. 16.

Velocity field distribution of REV under different flow directions. (a-c) is the flow velocity field of REV at point O of B-30 in the X, Y, and Z directions; (d-f) is the flow velocity field of REV at point O of B-20 in the X, Y, and Z directions; (g-i) is the flow velocity field of REV at point O of B-10 in the X, Y, Z directions.

Fractal Dimension and Heterogeneity of Pore Structure

The fractal dimension serves as a measure of the complexity and dynamic evolution of the internal structure of porous media. The structures of soil, rock, and similar porous materials frequently exhibit distinct fractal characteristics. Among the various methods available for calculating fractal dimensions, the box-counting method is most widely used in geotechnical engineering [64, 65]. This method involves covering the fractal target with small cubes of side length a and counting the number of cubes, N(a), needed to cover the fractal. The fractal dimension D of the porous medium is then calculated using Equation (14). Using AVIZO’s 3D fractal module, REV-specific fractal dimensions were obtained, with the B-10 model serving as an example. In this study, AVIZO’s 3D fractal module was applied to determine the fractal dimension of the pores at each REV location. For the B-10 model, the fractal dimension was computed to assess the complexity of the pore structure.

Table 3 lists the calculated fractal dimensions of the pores. The three-dimensional fractal dimension (D_3D) ranges from 2.5882 to 2.6101, indicating minimal overall variation. In contrast, the two-dimensional fractal dimension (D_2D) ranges from 1.5147 to 1.6006, with a maximum difference of 0.086. Based on Feng et al. [66], the ideal difference between the three-dimensional and two-dimensional fractal dimensions should be 1.0. In this study, the average difference between the three-dimensional and two-dimensional fractal dimensions is approximately 1.05. This indicates that the pore development and structure are relatively consistent across different REV locations. Although the uneven pore distribution at various locations leads to a slightly greater difference than 1.0, it still conforms to the fractal self-similarity characteristics described by Feng et al.

It is necessary to introduce the concept of pore structure heterogeneity to examine the causes of the non-uniform pore distribution. Previous numerical simulation methods often assume that pore distributions within geotechnical materials are homogeneous, which can result in errors when applied to practical engineering problems. It is essential to fully consider the heterogeneity of pore structures in geotechnical materials to reduce these errors [67].

Porosity, pore diameter, and pore shape, as essential parameters characterizing pore structures, can also be utilized to assess pore structure heterogeneity. The heterogeneity coefficient is defined using the coefficient of variation, calculated as the ratio of the standard deviation of the pore structure characterization parameters to their mean values. The magnitude of the heterogeneity coefficient reflects the non-uniformity of these parameters. The primary factors influencing pore structure heterogeneity can be identified by comparing the heterogeneity coefficients.

Fig. 17 shows that the heterogeneity coefficients of the pore structure characterization parameters are compared. Among the three parameters, the heterogeneity coefficient of pore shape is the highest, making it the dominant factor influencing pore structure heterogeneity in this study. In contrast, Wang et al. [61] in their research on coarse-grained sandstone, identified pore diameter as the primary factor affecting pore distribution heterogeneity, followed by porosity, with pore shape exerting a lesser influence. This variation can be related to the porosity of the studied materials. The Loose Gravel Soil structure, compared to sandstone reservoir structures, is more loosely packed and exhibits greater porosity. In addition, the virtual sample generation process leads to more uniform particle deposition than natural conditions, reducing the sensitivity of pore size to spatial structure. As a result, the influence of pore diameter on pore structure heterogeneity is slightly less significant than that of pore shape.

Permeability Anisotropy

The effect of non-uniform pore structure distribution on permeability can be evaluated through permeability anisotropy. The seepage flow behavior of pores at various locations can be analyzed by comparing the permeability ratios in different directions. Typically, the horizontal permeability of a sample exceeds its vertical permeability [68, 69], implying that the anisotropy indices X/Z and Y/Z should be greater than 1. In addition, the greater the deviation of the anisotropy index from 1, the higher the permeability anisotropy.

Table 4 indicates that all anisotropy indices in the REV pore model of B-10 exceed 1. This outcome is attributed to the use of the pre-compaction method, which effectively simulates the deposition of actual samples under natural conditions. During the pre-compaction process, the movement and rotation of aggregate particles result in the displacement of fine particles within the pores, modifying the pore structure. This observation is consistent with the findings of Masad et al. [68] on the anisotropic distribution of permeability in asphalt concrete, indicating that the permeability of Loose Gravel Soil in different directions adheres to objective principles.

In addition, permeability anisotropy does not exhibit a consistent trend across different spatial locations. Several factors, including the non-uniform distribution of particle sizes, the complexity of pore structures, and the irregularity of fractures and bedding, can influence this randomness. The inclusion of additional variables and uncertainties further complicates the prediction and simulation of this randomness.

During actual seepage flow in Loose Gravel Soil bank slopes, the shallow subsurface exhibits permeability anisotropy slightly greater than 1, aligning with the patterns summarized in Table 4 of this study. In deeper slope strata, however, the internal structure of the gravelly sandstone matrix becomes progressively denser due to sedimentation and stacking processes, resulting in lower porosity compared to shallow zones. Concurrently, fine particle migration blocks partial longitudinal seepage paths, weakening vertical seepage behavior (Fig. 18). Driven by hydraulic head differences and porosity gradients, water continuously migrates from high-potential to low-potential zones, redistributing saturation and thereby altering the slope’s permeability coefficient and reducing its shear strength. Concurrently, seepage-induced displacements develop within the slope. Although constrained by horizontal permeability, these displacements remain minor. Nevertheless, progressive water infiltration saturates the lower slope strata, increasing pore water pressure and diminishing slope strength. As the soil’s mechanical properties degrade, the upper slope becomes susceptible to sliding failure under gravitational loading.

Fig. 17.

Comparison of non-uniformity coefficients of pore structure characterization parameters.

The preceding discussion demonstrates the validity of the proposed modeling method from the perspectives of fractal dimension and permeability anisotropy. Compared to traditional DEM spherical modeling and nondestructive electronic detection techniques for extracting pore structures, the approach of scanning and modeling actual aggregate particles employed in this study better simulates the real structural state of geotechnical materials. In addition, extracting pore structures using this method is more convenient than techniques such as CT scanning. The pore structure can be altered, resulting in more accurate outcomes in seepage flow simulations by adjusting the DEM particle gradation and the precision of aggregate modeling. This provides a practical and effective methodology and framework for future investigations into microscopic pore-scale seepage flow problems, such as reservoir slope instability. However, during the accumulation of Loose Gravel Soil, aggregate particles form a specific inclination angle with the ground due to gravitational forces. Variations in this inclination angle inevitably affect the stability of slope seepage flow to varying extents [70], a factor not yet considered in this study. Therefore, it is recommended that future applications of the proposed modeling method incorporate the effects of aggregate particle inclination to more accurately simulate seepage flow conditions in slopes under different accumulation processes. In addition, the primary causes of engineering problems, such as dam breaches and slope instability following reservoir impoundment, are the heterogeneity of internal pore spaces and internal erosion caused by the migration of fine particles [25]. This study currently focuses on single-phase flow in seepage flow paths and does not thoroughly examine internal erosion. Investigating fluid-solid coupling issues from the perspective of microscopic pores using the DEM-CFD method is a critical direction for future engineering research. Therefore, future studies on pore-scale seepage flow should highlight the influence of internal erosion. In addition, to address the computational challenges associated with DEM, machine learning and deep learning techniques can be integrated to predict and monitor internal erosion [24], enabling early identification of potential risks in geotechnical slopes and reservoir dams and supporting timely preventive actions.

Fig. 18.

Influence of bank slope pore structure on seepage flow behavior.