Quantitative Study on Bed Load Proppant Transport during Slickwater Hydraulic Fracturing


 Bed load proppant transport is a significant phenomenon during slickwater hydraulic fracturing. However, the mechanism of bed load proppant transport is still unclear. In the present study, the proppant transport process during slickwater hydraulic fracturing was simulated with a coupled computational fluid dynamics- (CFD-) discrete element method (DEM) model, and the mechanism of the bed load proppant transport was analyzed. A model for calculating the mass flux of the bed load layer was proposed and verified with experimental results from the literature. The results show that bed load migration is an essential mechanism of proppant transport. When the shear force of fluid acting on the surface of the sand bank reaches the critical Shields number, the proppant in the upper layer of the sand bank begins to migrate in the form of bed load. The movement of the bed load layer increases the time for the sand bank to reach the equilibrium height. In addition, the mass flux of the bed load layer significantly affects the equilibrium height of the sand bank. The mass flux of the bed load layer decreases, and the equilibrium height increases as the proppant density, proppant diameter, or rolling friction coefficient and static friction coefficient of the proppant increase, but the mass flux of the bed load layer increases, and the equilibrium height decreases as the fluid viscosity increases.


Introduction
Slickwater hydraulic fracturing is a crucial technique for unconventional reservoir development. Due to the low viscosity of slickwater, the proppant settles out of suspension quickly to form a sand bank, and then bed load proppant transport becomes a significant mechanism of proppant transport [1]. Therefore, the study of bed load proppant transport is of great significance to analyze the proppant distribution in the fractures.
Experimental and numerical simulations are the key ways to study proppant distribution in fractures. Kern et al. [2] studied the proppant transport in fractures through experiments. They proposed that only when the sand dune reaches the equilibrium height can the proppant be further transported to the fracture in the form of fluidization. Woodworth and Miskimins [3] carried out the erosion slot flow experiments and performed a sensitivity analysis to show the sensitivities of various parameters, such as prop-pant density and fluid viscosity. Mack et al. [1] indicated that the static friction coefficient and coefficient of restitution of proppant influence the bed load proppant transport and further influence the proppant distribution in fractures. The advanced ceramic proppant can be transported more indepth into the fracture due to its higher coefficient of restitution and lower friction coefficient compared to sand and ceramic proppant. Tsai et al. [4] analyzed the proppant transport in fractures with the Eulerian-Lagrangian model. They found that the proppant settling behavior is influenced by proppant density, diameter, and injection rate. Blyton et al. [5] simulated the sand-carrying fluid flow in fractures with the CFD-DEM model. The simulation results indicated that the relative velocity of proppant and fluid is affected by sand ratio and the ratio of proppant diameter to fracture width. Zhang et al. [6] analyzed the impacts of various factors such as injection velocity and proppant size and density on proppant transport by simulation with the CFD-DEM model. Baldini et al. [7] presented similar to those of Zhang et al. [6]. Li et al. [8] studied the influence of interaction parameters between proppants and walls on proppant transport with a CFD-DEM model. They found that the friction coefficient of the proppant affects the bed load proppant transport. A lower friction coefficient is beneficial for bed load proppant transport. In addition, many scholars have carried out experiments or numerical simulations of proppant transport in complex fracture networks. Sahai et al. [9] experimentally studied the proppant distribution in fracture networks. Li et al. [10] researched the laws of proppant distribution in complex fracture networks through many experiments and the Eulerian-Eulerian simulation method. They proposed that using a large injection rate, low particle diameter, and sand ratio can improve the effective proppant placement in fracture networks. Kong and McAndrew [11] studied the effects of fracture width, fluid viscosity, sand load, particle size, and particle density on proppant distribution in fracture networks with the Eulerian-Eulerian model. These studies mainly focus on the controversial issue of whether the proppant can enter secondary fractures.
From the literature review, one can observe that the proppant transport in the fractures has two mechanisms: suspension and bed load. Most research is about the influence of various factors on proppant distribution in fractures, while research on bed load proppant transport is limited. In this paper, we intend to conduct quantitative research on the transport of bed load proppants in fractures. This study is crucial for clarifying the proppant distribution in fractures.

Numerical Simulation of Bed Load Proppant Transport in Fractures
The Eulerian-Eulerian [11][12][13][14], Eulerian-Lagrangian [4,15,16], and coupled CFD-DEM [5][6][7][17][18][19][20] models are commonly used for simulating proppant distribution in fractures. The Eulerian-Eulerian model treats particles as another fluid phase, and the interaction between the fluid and particles is calculated using the momentum exchange technique. Compared to the Eulerian-Eulerian model, the Eulerian-Lagrangian and CFD-DEM models can calculate the trajectory of each particle. However, the Eulerian-Lagrangian model is not suitable for simulating particle accumulation. Therefore, the coupled CFD-DEM model was chosen to simulate the bed load proppant transport during slickwater hydraulic fracturing, which is helpful for us to observe and analyze the bed load proppant transport. In addition, during slickwater hydraulic fracturing, the fracturing fluid leak off will lead to the reduction of the fluid velocity. However, as the fracture model used in the simulation is small, the influence of the fracture fluid leak off on fluid velocity is little. Thus, the fracturing fluid leakage was ignored in the simulations.

Models
2.1.1. Particle Control Equations. By Newton's second law, the particle movement is predominantly translational and rotational, and the expressions are as follows [8,21]: where v i refers to the translational velocity of particle i (m/s); w i refers to the angular velocity of particle i (rad/s); m i refers to the mass of particle i (kg); I i refers to the moment of inertia of particle i (kg m 2 ); F ci refers to the contact force acting on particle i (N); T ij refers to the torque acting on particle i by particle j or wall (N m); n refers to the number of total contacts for particle i (dless); F di refers to the drag force of the fluid acting on the particle (N); and F b refers to the buoyancy of the particle (N).
In the simulations, the particle-moving pattern was calculated with the Hertz-Mindlin (no slip) model [22,23].
where F cni and F cti are the contact forces in the normal and tangential directions (N), respectively, and F dni and F dti are the damping forces in the normal and tangential directions, respectively (N).
where E ′ refers to the effective Young's modulus (N/m 2 ); R ′ refers to the effective particle radius (m); α refers to the displacement of the amount of normal overlap (m); E refers to Young's modulus (N/m 2 ); v refers to the Poisson ratio (dless); and R refers to the particle radius (m).
where m ′ refers to the equivalent mass (kg); S n refers to the normal stiffness (N/m); v rel n refers to the normal relative velocity component (m/s); m refers to the particle mass (kg); e refers to the coefficient of restitution (dless); v ! refers to the particle velocity before collision (m/s); n ! refers to the 2 Lithosphere normal unit vector at the time of collision (dless); and r ! denotes the position vectors (m/s).
where S t is the tangential stiffness (N/m); ζ is the overlap along the tangent direction (m); G ′ is the effective shear modulus (Pa); and G refers to the shear modulus of the particle (Pa).
where v rel t denotes the relative velocity along the tangent direction (m/s).
where μ r is the rolling friction coefficient (dless); R i is the length from the centroid to contact point (m); and w i ! refers to the unit angular velocity vector at the contact point. The CFD and the DEM modules are coupled by the forces of fluid applied on the proppant.
where μ f is the fluid viscosity (Pa s); ε s denotes the volume fraction of the solid phase (dless); d s is the proppant density (kg/m 3 ); and Re s is the particle Reynolds number, which can be calculated with Equation (11) (dless).
2.1.2. Fluid Control Equations. In the CFD part, the fluid phase is solved directly by the Eulerian-Eulerian model [8,11,14].
where ε f refers to the volume fraction of the fluid (dless); ρ f refers to the density of the fluid (kg/m 3 ); v f refers to the average velocity of the fluid in a cell (m/s); t refers to the time (s); g refers to the gravitational acceleration (m/s 2 ); τ f refers to the stress-strain tensor (Pa); m refers to the number of particles in one fluid cell (dless); and F di refers to the momentum exchange between the fluid and solid phases (N). The turbulence fluctuation was accounted for a standard k-ε turbulence model, as follows [14]: where k refers to the turbulent kinetic energy (m 2 /s 2 ); ε refers to the turbulent dissipation rate (m 2 /s 3 ); G k represents the generation of turbulent kinetic energy (kg/(m 3 s)); χ k and χ ε is the turbulent exchange terms between the fluid phase and particles (kg/(m 3 s)); μ t refers to the turbulent eddy viscosity (Pa s); σ k = 1 and σ ε = 1:3 denote the Prandtl numbers of turbulent kinetic energy and turbulent dissipation rate (dless), respectively; and C 1ε = 1:44 and C 2ε = 1:92 are empirical constants (dless).

Initial and Boundary
Conditions. Figure 1 illustrates the computational domain of the fracture, where L = 0:24 m, W = 0:006 m, and H = 0:032 m. The sand-carrying fracturing fluid was injected into the fracture at the perforation and flowed out from the right outlet. In the simulations, the perforation was simplified to a small rectangle (0:004 × 0:006 m), which was located at the center of the left boundary of the fracture.
In the CFD part, we used a structured mesh that was generated by hexahedral elements, as illustrated in Figure 2. The cell size was 0:002 × 0:002 × 0:002 m, which ensured that the mesh size was larger than the particle diameter. The simulations applied a velocity boundary condition and a constant pressure boundary condition to the inlet and outlet, respectively. The timestep of the CFD solver was 8 × 10 −4 s.
In the DEM part, a particle factory was built at the position of the perforation. Randomly located particles were generated from the factory at a fixed rate. The generated rate of the particles was calculated based on the sand ratio (particle concentration) in the fracturing fluid. In addition, the velocity of the generated particles was the same as that of the injected fluid. The timestep of the DEM solver was 8 × 10 −6 s (29% of the Raleigh times), which is one hundredth of that in the CFD solver. The mesh size used in the DEM part is 2.5 times the particle radius.
The CFD and DEM modules are coupled through the drag force of fluid acting on particles, and the schematic of the CFD-DEM coupling calculation process is shown in Figure 3.  Table 1. The proppant settles out of suspension rapidly due to the effect of gravity when it is injected into the fracture, and a dune forms near the wellbore. The flow channel of the fluid above the sand bank decreases gradually as the height of the sand bank increases, which increases the shear force acting on the dune generated by the fluid. A bed load layer is generated when the shear force reaches the critical Shields number. However, the height of the dune increases gradually due to the continuous injection of proppant. When the shear force is enough to transport the injected proppants to move forward, the height of the dune stops increasing and reaches an equilibrium height. Due to the shear effect of fluid acting on the sand bank, there is an unsupported area at the front of the sand bank. The size of the unsupported area is related to the injection rate, the settlement velocity of the proppant, and the position of the perforation.

Results and Discussion
As seen from the process of proppant transport, the movement of the bed load layer only increases the time for the sand bank to reach the equilibrium height. The dune still reaches the equilibrium height near the wellbore rapidly, and the newly injected proppant moves forward with the bed load mechanism. The final morphology of the sand bank is less affected by the migration of the bed load layer. Figure 5 illustrates the mechanism of proppant migration in the fracture after the sand bank reaches the equilibrium height, in which H 1 , H 2 , and H 3 denote the heights of the fluid flow area, the bed load layer, and the immobile sand bank, respectively. (14) [1, [25][26][27].

Model. The Shields number is a common parameter that governs bed load proppant transport, as shown in Equation
where τ is the shear stress of fluid acting on the surface of the dune (N); ρ s refers to the proppant density (kg/m 3 ); and S is the Shields number (m 2 ). When the Shields number is low, the sand bank is immobile, and when the Shields number reaches the critical level, the proppant up the sand bank begins to move in the form of a bed load under the action of the fluid drag force. As the shear stress can be calculated with Equation (15) in a fracture [1], the critical average velocity can be expressed as Equation (16), which denotes the average velocity of the fluid up the sand bank at which bed load begins.
where V f refers to the average velocity of fluid up the sand bank (m/s); W refers to the fracture width (m); η refers to a coefficient, which is a constant (dless); V f c is the critical average velocity (m/s); and S c is the critical Shields number, which varies between 0.03 and 0.06 [27]. Wang and Zheng [28] established the equation of mass flux of bed load layer as Equation (17), and as seen from it, the mass flux of the bed load is affected by the density of the fluid and proppant, the proppant diameter, the shear   Lithosphere velocity of fluid acting on the bed, the static friction coefficient, and the rolling friction coefficient.
where c ≈ 1:5, μ m refers to the static friction coefficient between proppants (dless); μ d refers to the rolling friction coefficient of the proppant (dless); U * refers to the shear velocity of fluid acting on the bed (m/s); and ε s denotes the average volume fraction of proppant in the bed load layer (dless). The shear velocity can be calculated with Equation (18) [28,29]. Then, the mass flux of the bed load layer can be calculated with Equation (19).
As the coefficient η is unknown, the simulation results of the equilibrium height with various parameters were used to regress the coefficient. Table 2 shows the simulation results of the equilibrium height with various parameters.
In the simulations, assume that all the injected proppants migrate in the form of bed load when the sand bank reaches the equilibrium height. The mass flux of the bed load layer can be calculated with where N s is the number of injected proppants per unit time (1/s) and m s is the mass of proppant (kg). The thickness of the bed load layer can be calculated as where u creep is the average velocity of the bed load layer (m/s) and H 2 is the thickness of the bed load layer (m).

Lithosphere
The velocity distribution of the creep layer is complex [30,31]. The Saint-Venant model was used to calculate the velocity of the bed load layer along the Z-direction [28] to simplify the calculation.
where u creep is the velocity of the bed load layer (m/s) and y is the coordinate along the Y-direction (m). And the coordinate system is shown in Figure 5. Then, the average velocity of the bed load layer can be expressed as Þy + H 2 ð Þdy Assume that the velocity of the fluid can be approximated as the proppant velocity in the bed load layer. Then, where Q fin is the volume of injected fluid per unit time (m 3 /s). According to the simulation results, the coefficient η can be calculated by combining Equations (19)- (24). Table 3 shows the results of the coefficient η based on the simulation results in Table 2 when the average volume fraction of proppant in the bed load layer is 0.5. As shown, the average result of η is 268.
During actual slickwater fracturing, the mass flux of the bed load layer can be calculated with Equation (25). The shear velocity, the thickness of the bed load layer, and the height of the fluid flow area can be calculated with Equations (26)- (28). Then, the equilibrium height of the sand bank can be calculated with Equation (29).
where Q pin is the packing volume of injected proppant per unit time (m 3 /s); ε s ′ is the volume fraction of proppant dune, which is generally taken as 0.6 (dless); H EH is the equilibrium height of the sand bank (m); and H is the height of the fracture (m). According to the simulation results shown in Table 2 and the equation of mass flux of the bed load layer, the mass flux of the bed load layer decreases, and the equilibrium height increases as the proppant density, proppant diameter, or the rolling friction coefficient and static friction coefficient of proppant increase, but the mass flux of the bed load layer increases, and the equilibrium height decreases as the fluid viscosity increases.  Table 4. Figures 6 and 7 show the comparison of H 1 and H 2 between the experimental data and the calculated results by the built model, respectively. As shown, the calculated results by the built model fit well with the experimental results when the sand ratio is larger than 0.5%.

Conclusions
In this study, bed load proppant transport was studied. From the research, the following conclusions can be obtained: (1) Bed load proppant transport is an essential mechanism during slickwater hydraulic fracturing. The mass flux of the bed load layer mainly affects the equilibrium height of the sand bank (2) The mass flux of the bed load layer increases with an increase in the fluid viscosity and with a decrease in the proppant density, proppant diameter, and friction coefficients of the proppant (3) A model was built to calculate the mass flux of the bed load layer and the equilibrium height of the sand bank. The calculated results by the built model fit well with the experimental results when the sand ratio is larger than 0.5%

Data Availability
The results of the simulations or calculated by the built model (sand bank) used to support the findings of this study are included within the article. Previously reported (experimental results) data in Figures 5 and 6 were used to support this study and are available at doi:10.1016/ S0301-9322(02)00152-0. This study is cited at relevant places within the text as reference [23].

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.   Lithosphere