Numerical simulation of particulate suspension transport and permeability impairment in an actual rough fracture under normal stresses

The transportation and capture of particulate suspension in a rough fracture are widely encountered in hydrocarbon exploitation. According to a full review of the literature, traditional studies mainly focus on the placement and settlement of coarse particles such as proppants and lost circulation materials. Furthermore, a fracture is often considered to be equivalent to smooth parallel plates, whereas actual fractures are roughly walled with microconvex and are affected by normal stresses. The particulate suspension transport mechanism in a fracture under different normal stresses is still unclear. This paper proposes an approach to evaluate the permeability impairment caused by size‐exclusion of solid particles in a fracture under different normal stresses. The morphology of fracture flow space under different normal stresses is obtained by using optical scanning and FEM. Furthermore, the unresolved CFD‐DEM method is used to study the mechanism of the capture and migration of the particulate suspension in a fracture. The simulation results are in agreement with the experimental results. Furthermore, the effects of the particle size, particle concentration, injected particle volume, fracture occurrence, and applied pressure on particle invasion are discussed.


| INTRODUCTION
The transportation and capture of particulate suspension in rough fracture are widely encountered in hydrocarbon exploitation 1 (Figure 1A), such as solid invasion in water injection/drilling process, and fines migration in production. These particulate suspensions remain in the porous media and cause severe permeability impairment 2 ( Figure 1B,C). The fracture aperture distribution under different normal stresses varies, which indirectly leads to a different distribution of the particulate suspension in the fracture. Therefore, a comprehensive understanding of the mechanism of the particulate suspension transport in a rough fracture under different normal stresses is the basis for designing reasonable reservoir protection measures.
Until now, several theoretical models have been devoted to the description of the fluid-particle flow in porous media. Among theoretical studies, Civan 3 provided a detailed summary. At present, the theoretical models are still based on the deep filtration model (DBF) and the subsequent improved model. 4 In addition, numerous experimental methods and conclusions have been proposed, and the core flooding experiments remain the most advantageous. The effects of flow rate (applied pressure), particle concentration, and particle size on the particulate suspension transport and subsequent permeability impairment were studied using core flooding experiments. 5,6 However, it remains challenging to understand the particle distribution in porous media by using the core flooding experiments, as these experiments are still concerned with permeability loss. With the adoption of new technologies, the visualization experiments, computed tomography (CT), nuclear magnetic resonance (NMR), and fluorescence tracking have been introduced in the core flooding experiments. [7][8][9][10][11] Although coarse images can be obtained using these improved experiments, it is still difficult to obtain the quantized relationship between the particles and porous medium.
Considerable research based on numerical simulation methods has been conducted. The numerical simulations can be divided into several types, such as the particle quasi-fluid method, CFD-DEM method (including LBM-DEM, LES-DEM, and DNS-IB), particle tracking method, and Brownian dynamics simulation. [12][13][14][15][16][17][18] The existing studies demonstrated that the numerical simulation is an effective means to study the particulate suspensions transport and permeability damage to the porous media. The particle size-exclusion mechanism in a smooth fracture with a constrict has been extensively studied, and the general conclusion of these studies is that mechanical bridging is more prone to occur when the particle size is larger than one-third of the fracture width. [14][15][16]19 The phenomena of particle penetration and mud cake formation in packed beds have been studied widely. [20][21][22] However, the packed bed cannot mimic the real porous media. The coupled LBM-DEM-IMB method was used to study the particle migration and plugging in real pores. 23 However, the number of simulated particles was small due to the amount of calculation required. In recent years, the phenomenon of particle migration in a real fracture has been studied by combining the optical scanning technology, CT scanning technology, and numerical simulation. The LCM plugging process in the actual fracture obtained by optical scanning was analyzed by CFD-DEM. 24 Koyama et al 12 used optical scanning technology to obtain the fracture surface, and they used a particle tracking method to analyze the particle penetration in a sheared fracture. However, the particle-particle and particle-wall interaction was not extensively considered.
In general, the visualization and quantification of particulate suspension transport and permeability damage in rough fracture is still lacking, and the existing studies do not take into account the effects of normal stresses on the fracture. Generally, the fracture apertures are variable due to the change in the effective stress during drilling, fracturing, and production. 25,26 Moreover, the distribution of the fracture apertures changes significantly when the core is extracted from the reservoir (overpressure state) to the ground (zero pressure state). This phenomenon further results in the change in the solid particle migration and particle capture in a rough fracture.
The purpose of this manuscript is to provide a numerical simulation approach for modeling the particulate suspension flow in an actual fracture under normal stress and the corresponding changes in the permeability. In this paper, the unresolved CFD-DEM approach was used to understand the mechanism of the solid particle capture in an actual rough fracture. To this end, the initial fracture flow F I G U R E 1 (A) Developed fractures in tight sand cores and thin sections; (B) schematic diagram of particle invasion in a fracture; (C) formation damage caused by the particle capture in a rough fracture | 1167 ZHU et al. space distribution was obtained by using the optical scanning technique. In addition, the fracture deformation under different normal stresses was simulated by using the finite element method (FEM). The effect of particle size, particle concentration, injected particle volume, fracture occurrence, particle shape and applied pressure on the particulate suspension transport and subsequent permeability impairment was discussed. The numerical workflow used for this paper was as follows (Figure 2): 1. The profilometry of the actual fracture was carried out on a tight sand fracture sample by using a 3D structured light measurement system. Next, the initial actual fracture flow space was reconstructed using two obtained profilometries by employing the characteristic point cloud searching algorithm. 2. The rock blocks were imported into the FEM software after materialization and meshing. The bottom block was fixed, and normal stresses were applied to the top surface of the upper block. After each step of closure, the distribution of the fracture apertures was output. 3. For each exported fracture flow space, the unresolved CFD-DEM simulation was performed to understand the fracture permeability impairment and particle distribution after the particle invasion.

| Fracture reconstruction
The mechanical and hydraulic properties of a fracture under normal compression exhibit a clear dependence on the topography of the contacting surfaces comprising the fracture. In this paper, an advanced optical scanning system, namely, the 3D Eascan-D system ( Figure 3C) was used to scan the fracture surface topography. This system offers a high precision, and the finest sampling interval through each axis is 0.1 mm. The characteristic identification technique was employed in this study to avoid the subjective error in splicing the corresponding fracture surface. A prestep of scanning the outer surface of the fractured rock block is performed in a local coordinate system, and the fracture inner surface was not visible this step ( Figure 3A). In the following scan of the fracture inner surface, the new scanning data at different perspectives of the specimen help search the point cloud with the same distribution characteristics as in the prestep and align the new-found point cloud to the local coordinate system( Figure  3B). Although the two sides of the fracture specimen are scanned separately, the characteristic identification technique ensures that the point cloud of the two corresponding fracture surfaces is in the same coordinate system as the specimen in the prestep. Therefore, the fracture flow space reconstructed using 3D scanning can more accurately represent the conditions of an actual fracture ( Figure 3D). A cylinder tight sand rock core with a natural fracture was collected from the Jurassic formation (7725 m) of the Tarim basin. The fracture was controlled by the tectonic process and reformed by the diagenesis. The scanned fractured rock had a diameter and length of approximately 65 and 70 mm, respectively. The scanned point cloud data required filtering and denoising to improve accuracy. After the abovementioned procedure, the solid model was constructed in FEM software, and the subsequent step was initialized.

| Fracture measurement under normal stresses
In this work, ABAQUS (version 6.14) was used to analyze the static stress problems during fracture closure. An unstructured mesh of tetrahedral elements was generated in both the blocks. A mesh refinement was performed on the fracture surface in FEM to obtain a more accurate fracture aperture after the deformation. The following boundary conditions were applied ( Figure 3E): (a) a normal constraint on the bottom surface of the bottom block; (b) a normal displacement with an interval of 0.02 mm on the top surface of the upper block; (c) a hard contact between the corresponding fracture surfaces. After the calculation was completed, the average stress on the top surface of the upper block was extracted as the normal stress in each step, and the fracture aperture distribution was output.
It should be noted that the rock was set as a hard rock, which is linearly elastic and does not exceed the yield point; consequently, the deformation was reversible. In addition, the shedding of debris on the fracture surface was not considered. The Young's modulus and Poisson's ratio of the rock matrix were set as 40 GPa and 0.25, respectively.
The analysis results for the change in the averaged fracture aperture under different normal stresses are shown in Figure 4A and Table 1. The average fracture aperture decreased with an increase in the normal stress. Furthermore, the average fracture aperture decreased rapidly at the beginning, and the drop rate reduced after the normal stress reached a certain value. Before the normal stress reached 10 MPa, the average fracture aperture decreased from 0.218 to 0.1298 mm, and the fracture aperture drop rate was close to 40%. When the normal stress increased from 10 to 30 MPa, the fracture aperture reduced from 0.1298 to 0.0817 mm, and the fracture aperture drop rate was only 20%. Similarly, the contact area of corresponding fracture surfaces increased with the increase in the normal stress. Under a normal stress of 10 MPa, the contact rate of the corresponding fracture surfaces was 12.55%. When the normal stress increased to 30 MPa, the contact rate of the corresponding fracture surfaces increased to 22%. This Because the presence of a highly channelized fluid flow affects the particle transport, it is especially important to study the fracture aperture distribution instead of the average fracture aperture. To this end, the variation in the fracture aperture distribution contour map, distribution histogram, fractal dimension, and deviation with the normal stresses is also discussed in this paper (Figures 4B and 5; Table 1). The triangular prism method was used to calculate the fractal dimension. It can be seen from the contour map that as the normal stress increases, the contacted position continues to close, while the new contact occurs at some narrow positions. At the same time, the fracture aperture at the wide position gradually decreases.  It can be seen from the histogram that as the normal stress increases, the frequency of the large apertures decreases, while the frequency of the small apertures increases. As a consequence, the frequency histogram of the aperture becomes increasingly sharper and narrower, and the deviation and fractal dimension of the fracture aperture are also gradually reduced. The permeability change in a single-phase flow under different normal stress conditions assuming no particle injection occurs was calculated by using the calculation method defined in Section 3.1. It can be seen from the Figure 4A that when the normal stress increases to 10 MPa, the permeability is reduced by nearly 80%. When the normal stress increases to 30 MPa, the permeability is reduced by 95%. The downtrend of the average fracture aperture and permeability was successfully reproduced, and it demonstrated agreement with the experimental results obtained using the core from the same block in the Tarim basin. 28,29 3 | NUMERICAL METHOD FOR

PARTICULATE SUSPENSION TRANSPORT
Traditional methods such as the pseudofluid approach cannot analyze the particle movement accurately and intuitively as the particle-particle and particle-fluid collisions are not considered. Compared with the Eulerian-Eulerian approach, the individual motion of a solid particle is analyzed using Newton's motion theorem in CFD-DEM simulation. The particle capture, accumulation, and plugging can be observed. Therefore, the CFD-DEM coupling was used in this paper to simulate the movement and distribution of the particulate suspension in an actual rough fracture. Combined with the flow characteristics, the particle invasion performance was analyzed to provide ideas for the principle research and optimization of formation protection. The CFDEM coupling library (https ://www.cfdem. com), which is a bridge between the fluid and particle calculations, was used. This library calculates and transfers the exchange momentum terms between the fluid motion calculated using the OpenFOAM with the particle motion calculated using LIGGGHTS. 30

| Governing equations
The equilibrium equation of the solid particles can be expressed as: In Equation (1), F p,n and F p,t denote the normal and tangential contact forces, respectively; F p,g is the gravity, F p,e denotes the fluid-exerted force, such as drag force; r p,c and T p,r denote the particle radius and particle torque, respectively; and u p is the particle velocity.
According to the relationship between the particle size D p and CFD mesh size L, the resolved (D p ≫ L, such as immersion boundary method) and unresolved (L ≈ 3-10 D p , such as particle centroid method and divided particle volume method) CFD-DEM coupling methods are commonly used. The flow details around the particles can be resolved in the former method. However, this method is only used for the calculation of a few particles due to the expensive calculation cost. Only the force exerted on the particle mass center is considered in the unresolved CFD-DEM method, and thus, it is widely used in calculations involving several particles. 30 The continuity equation and momentum equation for particle-fluid flow in the unresolved CFD-DEM method can be written as: where l is the fluid volume fraction calculated using the divided particle volume method in this paper, l is the drill-in fluid density, u l is the fluid velocity, p is the pressure, l represents the liquid-stress-tensor, and R l,s is the fluid-particle momentum exchange. The PISO coupling is used to solve the abovementioned equations.
In this study, only the drag force, pressure gradient force, and buoyancy force were considered. Considering the large range of Reynolds numbers considered in this paper, the Di Felice 31 model was used to calculate the drag force: where C d is the solid-fluid drag coefficient, which can be calculated using Reynolds number Re s . According to previous studies, a better performance can be obtained when the CFD-DEM coupling interval is 50-100 DEM time steps.

| Simulation conditions
The CFD mesh quality is key to the computation accuracy and convergence in the CFD-DEM simulation. However, the geometry of the fracture flow space under different normal stresses is extremely complicated, which makes it difficult to satisfy the mesh quality requirement. Consequently, a reasonable simplification of the directly exported fracture flow space is required. Two simplified methods were adopted in this paper. The first stage of the operation was to reduce the simulation area (reduced to 25 mm in the length direction). The next step was to deduct the contact area (d = 0 μm) and an extremely narrow region (d < 0.01 mm) near the contact area. It should be noted that this deducted area has a minimal contributes to the fluid-particle flow in the complete fracture. Performing the abovementioned operation helps avoid the occurrence of certain abnormal meshes. Continuous particles at the desired concentration were created in the simulated borehole to mimic a steady particle inflow. A constant pressure boundary condition was applied on the inlet (top of the borehole) and outlet (end of fracture). The no-slip wall condition was applied on the fracture surface. The fracture model was placed horizontally or vertically to simulate the different fracture occurrences ( Figure 6). In this paper, a Newtonian flow was assumed to be laminar. The fluid parameters were derived from the drill-in fluid parameters of the tight sandstone reservoirs in the Tarim basin. The particles were assumed to be calcium carbonate, which is the most commonly used solid particle in the drilling operations in this area ( Table 2).
The determination of the CFD mesh size depends on the balance among the fracture surface detail, fluid detail, and computational stability of the fluid volume fraction. To check the mesh independence, fluid cells with different sizes were generated. The maximum particle diameter used in the simulations of this study was 0.25 mm, and thus, five different fluid cell sizes (0.25, 0.5, 0.75, 1, 1.25, and 1.5 mm) were generated. The applied pressure was 0.1 MPa, and the particle concentration was 1%. The change in the mass of the retained particles (particle velocity < 0.1 m/s) with the machine time was monitored, as shown in Figure 7. The mass of the retained particles was noted to be stable when the fluid cell size was >1 mm. Therefore, a mesh size of 1 mm was chosen for the simulation.  Figure 8 shows retained particle distribution and the fluid streamlines in the fracture under the different particle sizes for a constant total injected particle volume (0.05 PV). Figure 9 compares the dimensionless remaining permeability under the different particle sizes. The term '0.05 PV' denotes that the volume of the cumulative injected particles is 5% of the fracture flow space volume. The dimensionless remaining permeability is used to evaluate the severity of the permeability damage after particle suspension invasion, and it can be defined as where K r , K, and K 0 denote the dimensionless remaining permeability, permeability before the particle invasion, and permeability after the particle invasion, respectively. As the normal loading increases, the fracture aperture gradually decreases, and the contact area increases, resulting in a more complicated void space of the fracture. The particles must change their paths to bypass the contact area and the positions in which the aperture is smaller than the particle size. Therefore, the particle path becomes more tortuous and variable. 12 Not all the particles can migrate smoothly in the fracture. Some particles will be retained around the contact area or in certain low permeability areas due to the settlement, hydrodynamic effects, and mechanical bridging. Due to the fracture closure, the particles gather closer to the inlet as the normal stress increases at the same particle size ( Figure  8). Once the particles cannot effectively pass through the fracture, the dimensionless permeability retention considerably changes.

| Effect of particle size
In the case of a fine particle size, only a few particles are retained in the nonflowing area, and the permeability is not considerably affected (Normal stress = 0, 2.8, and 6.9 MPa in Figure 1A). As the particle size increases, numerous particles remain in the shallow part of the fracture, and the dimensionless permeability remaining rate is significantly reduced (Normal stress = 0, 2.8, and 6.9 MPa in Figure 8B). As the particle size continues to increase, most particles will collect near the inlet (mouth plugging) and block all the flow paths, preventing other particles from entering the fracture (Normal stress = 14.1 MPa in Figure  8B). However, the permeability of the fracture is still high because of the large particle size, according to the Kozeny-Carman equation: where a is the unity factor, and the is the void fraction.

| Effect of particle concentration
In this study, four particle volume concentrations of 1%, 3%, 5%, and 10% were simulated. The total injected particle volume remained unchanged (0.05 PV) during the simulation. The retained particle distribution and the streamlines under a particle concentration of 5% are shown in Figure 10. Figure 11 compares the dimensionless remaining permeability under different particle concentrations. It can be seen from Figure 11 that the permeability remaining rate decreased significantly with an increase in the concentration. However, the trend of decrease in the permeability remaining rate varied with the changes in the normal stress. When the normal stress was 0 MPa, the permeability remaining rate decreased slowly. This phenomenon occurred because a considerable particle retention did not occur when the particle concentration was <10%. However, some slight plugging was noted at some narrow positions, which led to a slight decrease in the permeability ( Figure 10A).
When the normal stress values were 2.8 and 6.87 MPa, the permeability remaining rate decreased rapidly as the particle concentration increased from 1% to 5%. This phenomenon occurred because the number of particles entering the fracture at the same time increased significantly with the increase in the particle concentration, which resulted in the increase in the mechanical bridging probability. The most significant reflection of the increase in the mechanical bridging probability was that the particle retention became more severe (as noted by comparing Figure 10B/C with Figure 8A). However, the permeability remaining rate decreased slowly when the particle concentration increased from 5% to 10%. This phenomenon occurred because the continued increase in particle concentration did not lead to a continuous increase in the mechanical bridging probability. 15 When the normal stress was 14.1 MPa, the particle retention could occur at the same position regardless of the particle concentration. Therefore, the increase in the particle concentration did not lead to a significant change in the permeability remaining rate.

| Effect of injected particle volume
As seen from Figure 12, the change in the permeability damage with the cumulative injection particle volume can be classified into two types. (a) In the case in which the particle retention is slight, the permeability damage does not change significantly with the accumulated injected particle volume (normal stress = 0, 2.8, and 6.9 MPa in Figure 8A; Figure 10A). (b) In the case in which serious particle retention can occur when the particles are just injected, the permeability impairment becomes more severe with continuation of the particles injection. Therefore, a larger number of retained particles under the same applied pressure correspond to a greater seepage resistance.

| Effect of fracture occurrence
In geology, fractures can be divided into horizontal and vertical fractures according to the inclination angle. Considerable research on the particle migration in vertical fractures has been conducted by many researchers by using experiments and numerical simulations; however, these studies focused mainly on the placement of proppant. The widely accepted conclusions are that narrower fracture, higher particle concentrations (or frequent particle-particle interaction), smaller proppant size, and rough fracture surface can hinder the particle settlement. [32][33][34][35] However, the particulate suspension migration in the vertical fracture is considerably different from the proppant placement due to the difference in the particle size and injection velocity.
The effect of fracture occurrences on the permeability impairment after the particle invasion was investigated in this paper. The characterization of fracture occurrence was achieved by adjusting the direction of gravity. Figure 13 shows the retained particle distribution and streamline in the horizontal fracture, and these aspects can be compared with the corresponding characteristics for the vertical fracture shown in Figure 10. It can be seen from the comparison that the retained particle distribution and the streamline in the horizontal fracture and the vertical fracture are almost identical except in the lower right corner of Figure 13B. The comparison between the dimensionless remaining permeability in the horizontal and vertical fractures is shown in Figure 14. It can be seen that the permeability impairment in the vertical fracture is relatively more critical when the normal stress is 2.8 MPa. However, when the normal stress values are 0, 6.9, and 14.1 MPa, the permeability impairment in the horizontal and vertical fractures is nearly equivalent.
Overall, the particle settlement does not play a dominant role in the particulate suspension transport in a fracture due to the following three main reasons: (a) Because the pressure difference between the wellbore and the formation is high during the drilling process, the drag force applied to the particles is considerably greater than the gravity of the particles (especially that of the fine particles). (b) The fracture aperture is small, which hinders the particle settlement due to the wall-effect. 17 (c) The surface roughness hinders the particle settlement, such as in the form of mechanical bridging.

| Effect of applied pressure
Under the premise of the total amount of the cumulative injected particles, the particle retention and permeability loss at different applied pressures (0.1, 0.5, and 1.0 MPa) were studied. Figure 15 shows the retained particle distribution and streamlines under 0.5 MPa. The dimensionless remaining permeability and particle retention rates at different applied pressures are shown in Figure 16. The particle retention rate is defined as the ratio of the retained particles to the total injected particles. According to Equation (4), the drag force applied to the particles is proportional to the slippage velocity between the fluid and the particles. That is, when the particles are retained, the drag force applied to the particles is proportional to the fluid velocity. As the fluid velocity increases, the drag force is increased, which may result in the destruction of the mechanical bridging, especially for multiparticle bridge. 15 Therefore, the particles can migrate to the fracture depth or flow out from the fracture outlet.
It should be noted that no significant change occurs in the retained particle distribution under a normal stresses of 14.1 MPa because no particles can pass through the fracture smoothly even under a high applied pressure. However, in the fracture in which the particle bridge is unstable (normal stresses values are 0, 2.8, and 6.9 MPa), the mechanical bridge is easily destroyed by the pressure difference (as noted by comparing Figure 15B with Figure 10B). Consequently, the particle retention rate is significantly reduced, and in the final state, the dimensionless remaining permeability is high (Figure 16). It can also be seen from Figure 16 that the influence of the applied pressure on the permeability remaining rate and the particle retention rate is limited. As the applied pressure increases from 0.5 to 1 MPa, the permeability remaining rate and the particle retention rate remain stable.

| Damage mechanism and experimental validation
In this section, to clearly explain the mechanism of particle invasion damage, the distribution of the retained particles was projected onto the contour map of the fracture aperture ( Figure 17). It can be found that the particle retention can be divided into two processes: single-double particle bridging and multiparticle accumulation. In particular, as the fine particles move along the preferential tortuous flow pathways existing in the fracture, they are captured and retained by the size-exclusion at the positions in which the fracture aperture is equal to the injected particle diameter (single-double particle bridging). Thereafter, the flow pathway is truncated, and the subsequent particles are either captured or they bypass the zone occupied by the retained particles. After all the pathways in front of the flow direction are jammed, the subsequent particles can only be captured. The particle size and the fracture aperture are no longer required to be equivalent in this stage; this stage is called multiparticle accumulation. The particle size, particle concentration, applied pressure, and retained particle distribution of all the simulated data were placed on a figure to determine the intrinsic relationship among them ( Figure 18). Since the fracture aperture distribution is not uniform, a few large aperture channels extend considerable contributions to the particle-fluid flow. Specifically, when 1 2 W f mean < D p < W f mean , most of the retained particles are distributed near the fracture inlet and are shaped like islands at the leading edge of the fracture. These island-shaped particles tend to remain in location in which the fracture microconvex body suddenly increases. Therefore, the fracture roughness critically affects the distribution of the retained particles. When D p > W f mean , numerous particles accumulate near the fracture inlet, while the fracture leading edge remains clean. When D p < 1 2 W f mean , almost no serious particle retention occurs. Furthermore, the applied pressure influences the particle retention when W f mean ∕D p is close to 2.
The particle-fluid mixed flow in rough fracture depends on the distribution of the fracture aperture. It is difficult to directly verify the numerical simulation results in this paper with the experimental results because no identical fractures are present. However, many experiments pertaining to the particle invasion in the fracture have been performed, including the experiments with tapered slots and a fractured core. The experimental results are compared with the numerical simulation results of this paper to verify the F I G U R E 1 7 Projecting the retained particles (black balls) onto the fracture aperture contour map reliability of the proposed numerical method. Comparing the experimental results reported in the literature with the numerical simulation results (Table 3; Figure 4A), it can be noted that the proposed method provides a reliable tool for studying the particulate suspension transport in a fracture under different normal stresses and is expected to yield highly reliable results.

| CONCLUSION
In this paper, the particulate suspension transport and permeability impairment in a fracture under different normal stresses was simulated by coupling the unresolved CFD-DEM approach and FEM. The fracture mechanical behaviors under normal loading and simulated particle-fluid flow results are noted to be in qualitative agreement with the results of the published experimental studies. The following conclusions can be drawn from this study.
1. When normal loading increased, the frequency histogram of the aperture becomes increasingly sharper and narrower, and the deviation and fractal dimension of the fracture aperture are also gradually reduced. These findings indicate that the asperities are smoothed out during the normal loading. 2. The influences of some key model parameters such as the particle size, particle concentration, injected particle volume, fracture occurrence, and applied pressure were investigated. The increase in the particle size, particle concentration, and total number of injected particles will lead to a significant increase in the particle retention in the fracture, which in turn leads to an increase in the permeability impairment. The change in the fracture occurrence does not significantly affect the permeability damage and particle retention. 3. By projecting the distribution of the retained particles onto the contour map of the fracture aperture, it can be found that the process of particle retention can be divided into two processes: single-double particle bridging and multiparticle accumulation 4. The aperture characterization, including its evolution with the normal loading, is the key parameter for the reliable simulation of the particle-fluid mixed flow in the fractures. When 1 2 W f mean < D p < W f mean , most of the retained particles are distributed near the fracture inlet and are shaped like islands at the leading edge of the fracture. When D p >W f mean , most of the particles accumulate near the fracture entrance while the fracture leading edge remains clean. When D p < 1 2 W f mean , almost no severe particle retention occurs.