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.
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 . 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.  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  carried out the erosion slot flow experiments and performed a sensitivity analysis to show the sensitivities of various parameters, such as proppant density and fluid viscosity. Mack et al.  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 in-depth into the fracture due to its higher coefficient of restitution and lower friction coefficient compared to sand and ceramic proppant. Tsai et al.  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.  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.  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.  presented similar to those of Zhang et al. . Li et al.  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.  experimentally studied the proppant distribution in fracture networks. Li et al.  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  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.
2. Numerical Simulation of Bed Load Proppant Transport in Fractures
The Eulerian-Eulerian [11–14], Eulerian-Lagrangian [4, 15, 16], and coupled CFD-DEM [5–7, 17–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.
2.1.1. Particle Control Equations
2.1.2. Fluid Control Equations
2.2. Initial and Boundary Conditions
Figure 1 illustrates the computational domain of the fracture, where , , and . 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 (), 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 , 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 .
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 (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.
3. Results and Discussion
Figure 4 illustrates the result of proppant transport in the fracture simulated with the parameters in 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.
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 , , and denote the heights of the fluid flow area, the bed load layer, and the immobile sand bank, respectively.
4. Mass Flux of the Bed Load Layer
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.
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.
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.
4.2. Model Validation
The experimental data of Wang et al.  were used to verify the built model. The size of the fracture used in their experiments was , and the other parameters in their experiments were , , and . The static friction coefficient and rolling friction coefficient between the proppants (sand) in their experiments are not given. We took larger values than ceramic proppant , where and . The results with various injection rates of fluid and proppant by the built model are shown in Table 4. Figures 6 and 7 show the comparison of and 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%.
In this study, bed load proppant transport was studied. From the research, the following conclusions can be obtained:
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
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
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%
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 .
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
We acknowledge the funding by the National Natural Science Foundation of China (Grant No. 51804175).