Numerical study on particle transport and deposition in rough fractures

. The transport and deposition of particulate materials through fractures is widely involved in environmental engineering and resource development engineering. A 3D Computational Fluid Dynamics-Discrete Element Method (CFD-DEM) coupling method was used to investigate the particle and ﬂ uid ﬂ ow. The Gauss Model was applied to construct the rough surfaces. First, the numerical results were compared with the previous results and reasonable agreements were obtained. Second, the results indicated a novel ﬂ ow pattern of particles in rough fractures. Then, a comprehensive particle sedimentary analysis indicated that the deposition distance of particles was inversely proportional to the particle size and density ratio. In addition, the particle deposition rates were increased by the mean roughness and there was an uneven sediment distribution impacted by roughness. Reasons for this uneven sediment distribution were analyzed in detail. More-over, the bridge plugs of particles considering the closure of fractures were simulated as well. A part of particulate materials would be ﬁ ltered at the inlet due to size effect and the transport distance of entered particles decreased signi ﬁ cantly when the particle was large. A critical particle radius ( R < 0.27 mm) that can ﬂ ow through closure fracture in this work was found. This work can provide a clear insight into the migration and deposition characteristics of particles in the rough fractures underground.

Fluid density of phase f, kg/m 3 d n Overlap distance between particles i and j, m xp;i Angular velocity of the particle i, rad/s

Introduction
With the development of unconventional reservoirs, the flow mechanism in pores and fractures have attached the attention of scholars (Li and Voskov, 2018;Sun et al., 2017;Yang et al., 2019c;Zhu et al., 2016).Numerous publications have been reported focusing on pore-scale numerical investigations of the gas flow in porous media and many of the results are helpful to study reservoir permeability (Song et al., 2016;Wang et al., 2019a;Zhang et al., 2019;Zhao et al., 2020).Such as Tianyi Zhao et al. indicated gas flow characteristics in the shale matrix by Lattice Boltzmann Method (LBM) and the effects of slippage, adsorption layer, residual water saturation and surface area all were considered (Zhao et al., 2018).
Especially the dynamic sorption was included in LBM schemes which are robust and efficient by Zhang and Sun (2019).However, except for the fluid flow mechanism the permeability damage results from the plugging of particulate materials in pore or fractures are studied as well because the loss of the particulate material also is a thorny issue for underground engineering and resource development engineering (Lijun et al., 2019;Oliveira Jr. et al., 2012;Wang et al., 2019c).The improper operation or engineering accidents will cause particles to leak into fractures underground and then result in significant economic losses and damage to the reservoirs (Sun et al., 2017;Yang et al., 2019b).As the basis of this particulate pollutant control, it is necessary to investigate particle flow and deposition characteristics in fractures.In addition, particle migration and deposition in fractures are common challenges in many other environmental and industrial projects such as the transport of proppant particles in fracturing engineering, treatment and plugging of slurry losses and structure optimization of granular flowing experiments (Haddadi and Di Carlo, 2017;Patankar et al., 2000;Yan et al., 2019;Zeng et al., 2019).
There are many efforts have been made to prevent this hazard and many research groups investigated the particle flow underground by experimental and numerical methods (Haddadi and Di Carlo, 2017;Hemmati and Rafee, 2018;Li et al., 2020;Yang et al., 2019a).It was indicated that the fractures were the important transport channels of particulate materials (Wang et al., 2019b;Yan et al., 2016).The previous results also showed that the characteristics of particle flow were complicated and involved the inherent structure complexities of fracture surfaces, multiple properties of particles, and the difficulties of coupling calculation (Han and Cundall, 2017;Wang et al., 2017;Zhang et al., 2017).However, not many field measurements and visual analysis for particle transport and deposition through fractures underground have been reported, since it is infeasible to detect directly in real-time.Fortunately, the numerical method has been considered as an effective method to further study these complex processes that cannot be realized by existing experimental techniques.Generally, the simulation strategies for particle flow are divided into the Eulerian-Eulerian coupling method and the Eulerian-Lagrangian coupling method (Basu et al., 2015;Han et al., 2007b).A Eulerian-Eulerian coupling method is a macroscopic approach treating fluid and particle as continuum media which is not suitable for tracking movement characteristics.In contrast, the particle is regarded as a discrete phase in the Eulerian-Lagrangian coupling strategy.In this strategy, particle flow is computed by the Discrete Element Method (DEM), fluid flow is computed by the Lattice Boltzmann Method (LBM) or Finite Volume Method (FVM).In addition, the Immersed Boundary or Moving Immersed Boundary Method is introduced in the resolved coupling method and drag force formula is used in unresolved coupling conditions (Han et al., 2007a;Kuruneru et al., 2017;Zeng et al., 2016;Zhang et al., 2017).Some semi-resolved Computational Fluid Dynamics-Discrete Element Method (CFD-DEM) approaches with advanced particle volume fraction calculation methods are introduced recently as well (Cheng et al., 2018;Wang et al., 2019d;Zhang et al., 2018).The CFD-DEM coupling method is adopted in this manuscript because it has both computational efficiency and particle computing accuracy.
In the field condition, there are many factors that should be taken into consideration for particle flow simulation.The collisions between particles, the non-uniform sizes and densities of mixed particles, and the different amounts of every kind of particle all would affect the distribution of deposited particles in fractures (Huang et al., 2018;Zou et al., 2015).In relevant work, DEM was proposed as a suitable method to fully consider these factors (Wu et al., 2017;Zeng et al., 2016Zeng et al., , 2018;;Zhou et al., 2018).Compared with other particle tracking methods, such as DPM (Discrete Phase Model), the DEM has more advantages for its richness of contact force model and its accuracy of computing motion of particles with multi sizes, shapes, and densities.Especially, the DEM can compute particle deposition reasonably from the perspective of energy losses rather than using rigid boundary judgment conditions.
Random and complex morphology is one inherent feature of underground fractures, especially the surface roughness is increased due to acid corrosion or treatment when the calcite-enrichment area is dissolved in priority (Li et al., 2016;Ogilvie et al., 2006;Yang et al., 2018).Because of the fluctuation characteristics of fracture surface accord with normal distribution, the Gauss Model has been widely accepted for constructing surface profiles that more tally with the actual situation (Guo et al., 2015).The roughness has an important influence on the flow and deposition pattern of dispersed particles in fracture and extensive relative study has been done in experiments and numerical simulation aspects.On the one hand, the visual experiments were conducted.It was found that particle size, mass loading, and the roughness degree play roles on the particle velocity profiles (Sommerfeld and Kussin, 2004).The uneven distribution of particles in the rough fracture was studied in experiments as well (Huang et al., 2018;Tan et al., 2017;Zou et al., 2015).On the other hand, various numerical simulation works obtained more detailed information on particle deposition characteristics, such as the effects of the micro shape of the rough channel (Dorari et al., 2015;Hemmati and Rafee, 2018;Hong et al., 2018).However, the compressive numerical investigation is still worth implementing.
The aim of this work is to investigate the deposition characteristics of particulate materials in rough fractures underground.There are following three sections in this manuscript.Section 2 presents the details of the CFD-DEM coupling model for the particle-water flow and Gauss model for surface roughness.Section 3 validates the model compared with previous work and discusses the results of this manuscript.The influence of multiple factors on particle transport and deposition is analyzed comprehensively such as the surface roughness, fracture closure, particle injection height, particle density, particle size, and particle number difference of mixed particles.The conclusions are summarized in Section 4 eventually.

Methodologies 2.1 CFD-DEM coupling model
The particle and water flow in the rough fracture are investigated by the CFD-DEM coupling method.The CFD equation set is computed in a commercial software ANYSY Fluent to obtain the information of fluid field containing particles.The movement of particles is computed by the Discrete Element Method in the EDEM software.The coupling between these two phases is implemented by a User Defined Functions (UDFs) code in Fluent.Through this UDF, the numerical simulation results between two software can be shared and invoked in real-time.The whole flowchart about the coupling method is presented in Figure 1.
The movement of particles is computed by the Discrete Element Method and the interaction between contacted particles is characterized by the overlap as presented in Figure 2. The particle translation and rotation are computed by equations ( 1)-(3): where m p;i g represents the gravity; ṽp;i represents the translational velocity of the particle i; xp;i represents the angular velocity of the particle i; P FC P;i represents the total contact force on particle i; ñij represents the normal direction between the center of particle i and j; I p;i represents the particle rotational inertia; P T i represents the vector sum of the torques on the particle i.
The governing equations of fluid flow based on Navier-Stokes description are: where a f represents the fluid volume fraction; q f represents the density of fluid; ũ represents the velocity of fluid; sf represents the viscous stress tensor; F * pf f represents the source term of the external forces; FrP represents the pressure gradient force; FD represents the drag force and FL represents the lift force.
The fluid forces on the particles are calculated in the coupled CFD-DEM model and the detailed fluid force calculation formulas are presented in the Appendix.

Construction of fracture with rough surfaces
Roughness often was ignored by the traditional fracture model using two parallel plates due to the complexity of geometry but it influenced the micro flow characteristics in fracture a lot.Most of the morphology information of rough surfaces such as peak size and distributions of roughness peaks are consistent with normal distribution, thus rough surface can be computed and reconstructed by Gauss model statistically.In other words, these peaks map information is concluded in one mathematical model equivalently for fracture roughness description.
In the Gauss (GS) Model a Fourier analysis is used for generating the surface roughness z (x, y) in equations ( 7) and ( 8): where x and y represent the number of nodes in the direction X and Y of fracture geometry region respectively; l x is the surface length in the direction X and l y is the surface length in the direction Y; cl x is the correlation length in the direction X and cl x is the correlation length in the Y direction; z represents the height of every rough peak;  the ift is the function that returns the two-dimensional inverse Discrete Fourier Transform (DFT) computed by the Fast Fourier Transform (FFT) algorithm and F g is the Gaussian filter: The geometry information of rough fracture (Fig. 3), the properties of the fluid and particle used in this work are summarized in Table 1.

Validation
To verify the reliability of this model to simulate the dense particle and fluid two-phase flow in rough fracture the profile of packed particles in fracture simulated by the model in this work was compared with the reference result (Zeng et al., 2016).The simulation setting and parameters of cases for comparison are presented in Figure 4 and Table 2.The deposited particles are in red color and the fluid is in blue color in Figure 4.The particle dune in the reference is outlined by the dotted line in black.The height and shape of the particle dune in both results are constant with a slight difference.
In addition, to validate the accuracy of particle tracking of the CFD-DEM coupling model the particle flow through the channel was simulated and compared with previous work (Haddadi and Di Carlo, 2017).A group of particles is injected from the inlet and flow out from the outlet of the micro channel as shown in Figure 5.The trajectory of a representative particle is tracked and compared with that in reference.Thus, the particle tracking results computed by the CFD-DEM method have a reasonable match with the reference results.

Effects of the rough bottom
The influence of fracture surface features on particle flow patterns and deposition characteristics has attracted the attention of many scholars recently.In this paper, some novel migration characteristics of suspending particles (density ratio was 1) in rough fracture were investigated.
The particle trajectories in the rough fracture (R-P1 and R-P2) are different from those in the smooth parallel plates (S) as shown in Figure 6.The vibration amplitude on x, y and z directions (AX, AY and AZ) of the particle trajectory was computed to describe the particle flowing fluctuation.It was indicated that a certain particle moves forward in the three-dimensional channel and experiences up and down motion during swaying around.The overall movement is a spiral trajectory that cannot be visualized  by both two-dimensional simulations and three-dimensional simulations without considering rough surfaces before.The velocity vector of particles shown in Figure 7 indicates that the particles would wind around the peak of the rough bottom, while the other particles with less kinetic energy deposit in the valley of the rough bottom after a period of spiraling motion.In addition to the effects on the particle trajectories, the uneven distribution of deposited particles results from the surface roughness was found as well.The roughness increased the probability of particle-particle and particle-wall collisions.After a series of collisions, the kinetic energy of some particles was reduced to a level which was not enough to escape vortexes.Eventually, particles deposit (Fig. 7) and the uneven distribution of particles is   formed that more particles tend to deposit on the valleys of the rough bottom.
Furthermore, the effect of the mean surface roughness on the particle deposition ratio was computed.The mean roughness of fracture surfaces in this work is 0.264 mm, 0.539 mm, 0.787 mm and 1.080 mm, respectively.The length (l) of the investigated area is 80 mm and the particle deposition rate is defined in equation ( 11).It was found that fluid maximum velocity and deposition rate increased as mean roughness grown as shown in Figure 8.
Particle deposition rate ¼ Deposited particle number Total particle number :

Transport and deposition of mixed particles
Considering the actual mixture condition of particulate materials, the effects of density ratio and size on particle transport and deposition were investigated.

Effects of density ratios
The particles with various density ratios were injected into the fracture in the same way and each kind of particle has the same size and amount.The density ratio is 1.5, 2.5, 3.5 and 4.5 respectively.The particle radius is 0.25 mm and the number of every kind of particle is 100.It can be found that particles experienced separation due to the particles with a larger density ratio had smaller horizontal deposition distance.Furthermore, the relationship between the maximum deposition distance and density ratio is consistent with the parabola law as shown in Figure 9.Additional simulations with the density ratio of 2, 3 and 4 are implemented for validating this interpolation pattern and these data (in red) is fitting well with the dotted line in Figure 9.
According to the law of deposition distance for particles with different density in Figure 9, the amount adjustment or difference of each kind of particle would change the eventual distribution of deposited particles in fracture.The effects of the so-called amount difference were investigated in this work.The detailed particle number of every kind of density is listed in Table 3 and Case 1 is regarded as the reference case that all kinds of particles have the same   As shown in Figure 10c, there are much fewer lowdensity particles deposited at area C in Case 2 compared with reference Case 1. Due to the blue particles is the only kind in area C, the amount difference leads to a reduction in the total number in this area as well.In Figure 10b, Case 3 has fewer particles in green and Case 4 has fewer particles in red which both lead to sedimentation reduce in middle section B. Case 5 has fewer particles in the near entry area A than other four cases as shown as Figure 10d due to the dominance of the dense density particle in black.

Effects of particle sizes
Because it is hard to keep the radius of particle the same and many kinds of sizes were always mixed during the transport of particulate materials in fractures.The uneven size becomes one significant factor that will affect the transport process and distribution results in a fracture.The effects of particle sizes and their injection heights were investigated in this paper.The particles were injected randomly from the inlet area, whose size had a normal distribution (Fig. 11).The horizontal migration distances of every kind of particle are recorded in Figure 12 and a previous separation phenomenon is found.It was indicated that the horizontal transport distance is inversely proportional to the particle diameter but it is proportional to the injection heights.
The row of particles with various radius injected from high position experienced sorting.Smaller particles broke away from the cluster of particles gradually during the sedimentation and had farther deposition distance.Particles from medium and low injection height had a similar law of sorting during the settlement process.However, because of the shorter settling time from lower injection height, the separation was not obvious.Besides, big particles injected from a higher position and small particles from the lower position would land on the bottom at the same time and had close horizontal migration distance.Consequently, the particles from different initial layers were mixed together after settlement.

Effect of rough sides
Particle migration and deposition characteristics in fracture with rough top and bottom walls were investigated and potential influence from uneven density, amount, size and injection height was analyzed in detail.However, tortuosity of fracture, in other words, the rough sides should be considered as well.It is a significant feature of fractures and has effects on the flowing and deposition patterns of particulate materials.The particle migration and deposition process through fracture with rough surface and sides were simulated.It was predicted that a part of particles would deposit in the corner of rough sides in priority.
On the one hand, the zigzag corner of rough sides affected the fluid flow such as inducing some local vortices.The motion of inertia particles was closely related to the fluid flow characteristics and a part of particles would go through these corners.Under the effects of fluid forces and collisions, the kinetic energy of particles experiences a fluctuation and decreases to zero eventually.The particles deposited when kinetic energy reduced to a level that was not enough to escape vortexes.To analyze the particle movement in this fracture a single particle trajectory from the top view was computed.As presented in Figure 13, in the region inside the blue dashed box corresponding to region 1, the particle dives along the channel under the drag of fluid; in the green region corresponding to region 2, particle still has enough kinetic energy to escape the sharp corner of the rough side in a parabolic path; in the region 3 corresponding to the region with orange color, particle swirls in the vortices of the fluid then it still moves forward; in region 4 that inside the red dashed box, the kinetic energy of the particle has been consumed a lot after collisions and the particle has a tendency to stop at the corner.Particles finally reach the region 5 indicated by the purple box where kinetic energy consumption further and the  particle velocity gradually approached 0 after swirling and collisions.The particle deposit in region 5 eventually in Figure 13.
Moreover, the impacts of rough bottom and sides on the particles were investigated compressively.The roughness induced local fluid vortices, increased particle collisions, and kinetic energy dissipation.Driven by the fluid, most of the particles flowed into the valleys of the bottom and zigzag corners of sides and lost kinetic energy after collisions.Consequently, a number of particles tended to deposit in the valley and zigzag corners in priority.This uneven sedimentary distribution can be observed clearly from simulation results in Figure 14.

Effect of closure of the rough fracture
When considered the closure of rough fracture, the openings of different parts of the inlet are different and there are some contact spots due to the initial apertures between rough surfaces are uneven.A part of particulate materials would be filtered due to size effect.For the particles that entered fracture, some flowed around the rough protruding peaks while others formed the bridge plugging.The fingerlike particle distribution was formed in rough fractures rather than the particle layers.It was found that both the number of entering particles and transport distance decrease as particle size increases in Figure 15.When the particle is rather small (R < 0.27 mm in this work), the closure just has a weak impact on the number of particles entered fractures.The impact becomes significant for cases with bigger particles (R > 0.27 mm in this work) and particle transport distance is short due to the bridge plugs.

Conclusion
As the basis of particulate pollutant control, it is necessary to investigate particle flow and deposition characteristics in fractures.The CFD-DEM coupling numerical method was used to compute the migration and deposition process of particles in rough fractures.In addition, the Gauss Model was applied to construct the rough walls of fracture.These numerical simulation results are summarized as followed: 1.The simulated trajectories showed that the particles within the rough fracture held different flow patterns from those in smooth condition.The particles experienced up and down motion, swaying around and overall trajectories were spiral.2. A comprehensive sedimentary analysis indicated that the deposition distance of particles was inversely proportional to the particle size and density ratio.The separation and deposition of mixed particles were simulated efficiently as well.3. It was found that both rough surfaces and sides of fracture have impacts on the deposition and distribution of particles.Bottom mean roughness increased particle deposition rates.Particles tended to settle in the valleys of the rough bottom and the zigzag corners of rough sides.Possible reasons for this uneven sediment distribution were analyzed: firstly, local fluid vortices were induced by surface roughness; then particles flow into the valleys and corners driven by fluid and lost kinetic energy gradually after collisions; consequently, most of the particles deposited in the valley and corners in priority.4. The bridge plugs in closed rough fractures were simulated efficiently.A part of the particle would be filtered and stopped form entering fracture due to size effect and the transport distance decreased significantly when the particle was large.A critical particle radius (R < 0.27 mm) that can flow through closure fracture in this work was found.
Overall, the trajectory of sparse particles and the distribution of deposited particles were influenced by wall roughness while the law of particle separation was more related to the particle parameters.This work gave a clear insight into the transport and deposition process of particles in fractures which could provide guides for hazardous and pollutant particle treatments underground.To investigate the particle transport and deposition in the environment and engineering further, the improved related numerical model and experiments are still necessary.Such as the model of this work is in the lab-scale and more efforts are needed to achieve the simulation in the field scale of hydraulic fracturing.In addition, the fracture propagation and the leakoff are not considered in this study and more effects should be done in the future.
Drag coefficient of the particle lFluid viscosity, Pa Á s q f

Fig. 3 .
Fig. 3.The rough fracture geometry and the investigated region.

Fig. 4 .
Fig. 4. The particle dune profile in the Eulerian view.The dune outline of the reference result (Zeng et al., 2016) was outlined by a dotted black line.

Fig. 5 .
Fig.5.The particle trajectory comparison between the simulation and reference results(Haddadi and Di Carlo, 2017).

Fig. 6 .
Fig. 6.(a) Comparison of particle movement characteristics in rough and smooth fracture.(b)-(d) Particle trajectory on the X-Z section (view along the direction of movement), Y-Z section (side view), and X-Y section (top view).

Fig. 7 .
Fig. 7.The velocity vector of particles near the rough bottom and the decrease of kinetic energy due to collisions.

Fig. 8 .
Fig. 8.The velocity of fluid and particle deposition rates under different roughness.

Fig. 9 .
Fig. 9. Maximum deposition distance of particles with various density ratios.

Fig. 10 .
Fig. 10.Particle amount distribution in areas A, B, and C.

Fig. 12 .
Fig. 12.The influence of particle diameter and injection position on the deposition distance of particles (the three rows of particles in orange were observed).
Fig. 13.(a) The trajectory of certain selected particles.(b) The decrease process of average kinetic energy.

Table 1 .
Simulation settings and parameters.

Table 2 .
The properties of particle and fluid.