*2.1. Governing Equations*

Based on kinetic theory, continuity and momentum equations were established for fluid and solid phase, respectively [21]. Continuity equation of the solid phase (m = s) and the fluid phase (m = l) is given by:

$$\frac{\partial}{\partial t}(\alpha\_{\\$}\rho\_{\\$}) + \nabla \cdot (\alpha\_{\\$}\rho\_{\\$} \, v \to\_{\\$}) = 0 \tag{1}$$

$$\frac{\partial}{\partial t}(\alpha\_l \rho\_l) + \nabla \cdot (\alpha\_l \rho\_l \,\, \boldsymbol{v} \to\_l) = 0 \tag{2}$$

where αs and αl is the solid and liquid phase volume concentration, dimensionless; ρs and ρl is the density of solid and liquid phase, kg/m3; → *v s* an → *v l* is the velocity vector of solid and liquid phase, m/s. Momentum equation for the fluid phase is given by:

$$\frac{\partial}{\partial t} \left( a\_l \rho\_l \stackrel{\rightarrow}{\boldsymbol{\upsilon}}\_l \right) + \nabla \cdot \left( a\_l \rho\_l \stackrel{\rightarrow}{\boldsymbol{\upsilon}}\_l \stackrel{\rightarrow}{\boldsymbol{\upsilon}}\_l \right) = -a\_l \nabla p\_l + \nabla \cdot \overline{\overline{\boldsymbol{\tau}}}\_l + a\_l \rho\_l \stackrel{\rightarrow}{\boldsymbol{g}}\_l + \beta \left( \stackrel{\rightarrow}{\boldsymbol{\upsilon}}\_s - \stackrel{\rightarrow}{\boldsymbol{\upsilon}}\_l \right) \tag{3}$$

Momentum equation for the solid phase is given by:

$$\frac{\partial}{\partial t}(\boldsymbol{a\_{\delta}}\boldsymbol{\rho\_{\delta}}\stackrel{\scriptstyle \boldsymbol{v}}{\boldsymbol{v}\_{\delta}}) + \nabla \cdot \left(\boldsymbol{a\_{\delta}}\boldsymbol{\rho\_{\delta}}\stackrel{\scriptstyle \boldsymbol{v}}{\boldsymbol{v}\_{\delta}}\stackrel{\scriptstyle \boldsymbol{v}\_{\delta}}{\boldsymbol{v}\_{\delta}}\right) = -\boldsymbol{a\_{\delta}}\boldsymbol{\nabla}\boldsymbol{p\_{l}} + \nabla \cdot \stackrel{\scriptstyle \boldsymbol{\overrightarrow{\nabla}}}{\boldsymbol{\tau}}\_{\delta} + \boldsymbol{a\_{\delta}}\boldsymbol{\rho\_{\delta}}\stackrel{\scriptstyle \boldsymbol{v}}{\boldsymbol{g}} + \boldsymbol{\beta}\Big{(}\stackrel{\scriptstyle \boldsymbol{\overrightarrow{\nabla}}}{\boldsymbol{v}\_{l}} - \stackrel{\scriptstyle \boldsymbol{\overrightarrow{\nabla}}}{\boldsymbol{v}\_{\delta}}\Big{)}\tag{4}$$

where τl and τs is the shear stress tensor of the fluid phase or solid phase, Pa; β is the momentum exchange coefficient between two phases, dimensionless; and g is the acceleration of gravity, m/s2.

The granular temperature Θs for the solid phase represents the kinetic energy of the random motion of the solid particles [26]. The transport equation derived from the kinetic theory is taken as the following form:

$$\frac{3}{2}\rho\_s\Big|\frac{\partial}{\partial t}(a\_s\Theta\_s) + \nabla \cdot \left(a\_s\overrightarrow{\boldsymbol{v}}\_s\Theta\_s\right)\Big| = \nabla \cdot \left(\kappa\_s\nabla\Theta\_s\right) + \tau\_s : \nabla\overrightarrow{\boldsymbol{v}}\_s - \boldsymbol{f}\_s + \Pi\_\Theta \tag{5}$$

where Θ*s* is the granular temperature, m<sup>2</sup>/s2; κ*s* is the diffusion coefficient, kg/(m·s); *Js* is the rate of the dissipation granular energy by inelastic collisions per unit volume, kg/(m·<sup>s</sup>3); ΠΘ is the rate of the dissipation granular energy by inelastic collisions per unit volume, kg/(m·<sup>s</sup>3).

In the artificial fracture, due to the low viscosity and high pumping rate, the flow state of the liquid phase may become turbulent flow. According to the research of Pfleger et al. [29], the predictions of the liquid phase turbulence kinetic energy *k* and its rate of dissipation ε are obtained from the following modified liquid phase *k* − ε turbulence model:

$$\frac{\partial}{\partial t}(\alpha\_l \rho\_l k) + \nabla \cdot \left(\alpha\_l \rho\_l \overrightarrow{\nabla}\_l k\right) = \nabla \cdot \left(\alpha\_l \left(\mu\_l + \frac{\mu\_l^t}{\sigma\_k}\right) \nabla k\right) + G\_{k,l} + \Pi\_k - \alpha\_l \rho\_l \varepsilon \tag{6}$$

$$\frac{\partial}{\partial t}(a\mathbf{q}\rho\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}) + \nabla \cdot \left(a\mathbf{q}\rho\boldsymbol{\varepsilon}\overrightarrow{\boldsymbol{\sigma}}\_{1}\boldsymbol{\varepsilon}\right) = \nabla \cdot \left(a\mathbf{q}\left(\mu\_{\mathrm{l}} + \frac{\mu\_{\mathrm{l}}^{\mathrm{f}}}{\sigma\_{\mathrm{c}}}\right)\nabla\boldsymbol{\varepsilon}\right) + a\mathbf{q}\frac{\boldsymbol{\varepsilon}}{\boldsymbol{k}}(\mathbf{C}\_{1\boldsymbol{\varepsilon}}\mathbf{G}\_{\mathrm{k},\mathrm{l}} - \mathbf{C}\_{2\boldsymbol{\varepsilon}}\rho\_{\mathrm{l}}\boldsymbol{\varepsilon}) + \Pi\_{\mathrm{d}}\tag{7}$$

where *k* and ε is turbulence energy and its rate of dissipation, m<sup>2</sup>/s<sup>2</sup> and m<sup>2</sup>/s3; μ*tl* is the turbulence viscosity coefficient and μ*l* is the liquid molecular viscosity, Pa·s, *Gk*,*<sup>l</sup>* is the production term of the turbulence kinetic energy, kg/(m·<sup>s</sup>3). According to Pfleger et al. [29], the source terms Πk and Πε are not considered yet. The constant parameters used in different equations are as follows: σ*k* = 1.0, σε = 1.3, *C*μ = 0.09, *C*1ε = 1.44, *C*2ε = 1.92.

## *2.2. Constitutive Equations*

In the simulation presented here, the accumulation of high-concentration proppant in the proppant dune and the low-concentration on the dune appear at the same time. So, the drag correlation of Gidaspow et al. [30] was used due to the flexibility under different solid-phase concentrations:

$$\beta\_{ls} = \begin{cases} \frac{3}{4} \text{C}\_{D} \frac{\rho\_{l} a\_{l} a\_{s} \left| \stackrel{\rightarrow}{v}\_{l} \right| \stackrel{\rightarrow}{v}\_{l} \left| \stackrel{\rightarrow}{v}\_{s} \right|}{d\_{s}} \alpha\_{l}^{-2.65} & \alpha\_{l} \ge 0.8\\ \frac{150 a\_{s} (1 - a\_{l}) \mu\_{l}}{\alpha\_{l} d\_{s}^{2}} + \frac{1.75 \rho\_{l} a\_{s} \left| \stackrel{\rightarrow}{v}\_{l} \right| \stackrel{\rightarrow}{v}\_{l} - \stackrel{\rightarrow}{v}\_{s}}{d\_{s}} & \alpha\_{l} < 0.8 \end{cases} \tag{8}$$

where αl is the liquid volume fraction and the drag coefficient between phases *CD* is given by:

$$\mathbf{C}\_{D} = \begin{cases} \begin{array}{c} \frac{24}{Re\_{s}} \Big[ 1 + 0.15 \left( Re\_{s} \right)^{0.687} \right] & Re\_{s} < 1000\\ 0.44 & Re\_{s} \ge 1000 \end{array} \tag{9}$$

where *ds* is the particle diameter, m. *Res* is Reynolds number defined by the relative velocity between two phases:

$$Re\_s = \frac{\rho\_l d\_s \alpha\_l \left| \overrightarrow{\upsilon\_s} - \overrightarrow{\upsilon\_l} \right|}{\mu\_l}. \tag{10}$$

The Reynolds stress tensor for the liquid phase is given by:

$$\pi\_1 = \alpha\_l \mu\_l \left(\overrightarrow{\nabla}\overrightarrow{\boldsymbol{\upsilon}}\_1 + \left(\overrightarrow{\nabla}\overrightarrow{\boldsymbol{\upsilon}}\_1\right)^T\right) \tag{11}$$

where μ*t*= <sup>μ</sup>*l*+μ*tl*is liquid phase effective viscosity, Pa·s; μ*tl*is the turbulent viscosity, Pa·s.

 The stress tensor for the solid phase is given by:

$$\tau\_s = \left(-(p\_s + p\_\mathbf{f}) + \eta \mu\_\mathbf{b} \nabla \cdot \overrightarrow{\boldsymbol{\upsilon}}\_s\right) \boldsymbol{I} + \mathcal{2}(\mu\_s + \mu\_\mathbf{f}) \mathcal{S}\_\mathbf{s} \tag{12}$$

where *I* is the second-order identity tensor, dimensionless; *Ss* is the stress strain of solid phase, s<sup>−</sup>1; *p*f is the frictional pressure of the solid phase, Pa.

The strain stress for the solid phase is given by:

$$S\_{\sf s} = \frac{1}{2} \left[ \overrightarrow{\nabla \boldsymbol{v}\_{\sf s}} + \left( \overrightarrow{\nabla \boldsymbol{v}\_{\sf s}} \right)^{T} \right] - \frac{1}{3} \left( \overrightarrow{\nabla \cdot \overrightarrow{\boldsymbol{v}}\_{\sf s}} \right) \boldsymbol{I}. \tag{13}$$

Solid kinetic-collisional pressure is given by:

$$p\_s = \alpha\_s \rho\_s \Theta\_s \left( 1 + 4\eta \alpha\_s g\_0 \right) \tag{14}$$

where η = (1 + *e*)/2 and *e* is the particle-wall restitution coefficient, dimensionless.

The radial distribution function at contact is given by:

$$\text{g}\_0 = \frac{1 - 0.5\alpha\_\text{s}}{\left(1 - \alpha\_\text{s}\right)^3}.\tag{15}$$

Solid-phase shear viscosity model is given by:

$$\mu\_{\rm s} = \left(\frac{2+c}{3}\right) \begin{bmatrix} \frac{\mu^\*}{\mathfrak{z}\eta\eta(2-\eta)} \left(1 + \frac{8}{5}\eta\alpha\_{\rm s}\mathfrak{z}\_0\right) \\ \left(1 + \frac{8}{5}\eta(3\eta - 2)\alpha\_{\rm s}\mathfrak{z}\_0\right) + \frac{2}{5}\eta\mu\_{\rm b} \end{bmatrix} \tag{16}$$

$$\mu^\* = \mu \left[ 1 + \frac{2\beta\mu}{\left(\alpha\_s \rho\_s\right)^2 \mathcal{g}\_0 \Theta\_s} \right]^{-1} \tag{17}$$

$$
\mu = \frac{5\rho\_{\text{s}}d\_{\text{s}}\sqrt{\Theta\_{\text{s}}\pi}}{96} \tag{18}
$$

$$
\mu\_{\rm b} = \frac{256}{5\pi} \mu a\_{\rm s}^2 g\_0 \tag{19}
$$

where μ∗ is the granular viscosity with the effect of the interstitial fluid, Pa·s; μ is the solids phase dilute granular viscosity, Pa·s; and μ*b* is the bulk viscosity of the solid phase, Pa·<sup>s</sup> [22].

Granular energy conductivity coefficient is given by:

$$k\_{\mathfrak{s}} = \left(\frac{\kappa^\*}{\mathcal{S}^0}\right) \begin{bmatrix} \left(1 + \frac{12}{5}\eta \alpha\_{\mathfrak{s}} \mathcal{g}\_0\right) \left(1 + \frac{12}{5}\eta^2 (4\eta - 3)\alpha\_{\mathfrak{s}} \mathcal{g}\_0\right) \\ + \frac{64}{25\pi} (41 - 33\eta)\eta^2 \left(\alpha\_{\mathfrak{s}} \mathcal{g}\_0\right)^2 \end{bmatrix} \tag{20}$$

$$\kappa^\* = \kappa \left[ 1 + \frac{6\beta\kappa}{5(\alpha\_s \rho\_s)^2 g\_0 \Theta\_s} \right]^{-1} \tag{21}$$

$$\kappa = \frac{75\rho\_s d\_s \sqrt{\pi \Theta\_s}}{48\eta (41 - 33\eta)}\tag{22}$$

where κ<sup>∗</sup> is the granular conductivity with the effect of the interstitial fluid, and κ is the solid phase dilute granular conductivity, W/(m·K).

Considering the complexity between the particles in the proppant dune with high volume fraction, solid-phase frictional pressure and viscosity model [22] are given by:

$$\frac{p\_{\rm f}}{p\_{\rm c}} = \left(1 - \frac{\nabla \cdot \vec{\boldsymbol{\upsilon}}\_{\rm s}}{n\sqrt{2}\sin(\rho)\sqrt{\mathbf{S}\_{\rm s} : \mathbf{S}\_{\rm s} + \Theta\_{\rm s}/d\_{\rm P}^2}}\right)^{n-1} \tag{23}$$

$$\mu\_{\rm f} = \frac{p\_{\rm f} \sin(\varphi)}{\sqrt{2} \sqrt{S\_{\rm s} : S\_{\rm s} + \Theta\_{\rm s} / a\_{\rm P}^2}} \left\{ n - (n - 1) \left( \frac{p\_{\rm f}}{p\_{\rm c}} \right)^{\frac{1}{n - 1}} \right\} \tag{24}$$

where

$$p\_c = \begin{cases} \begin{array}{l} 10^{25} \left(\alpha\_s - \alpha\_s^{\text{max}}\right)^{10} & \alpha\_s > \alpha\_s^{\text{max}}\\ \text{Fr} \frac{\left(\alpha\_l - \alpha\_s^{\text{min}}\right)^r}{\left(\alpha\_l - \alpha^\*\right)^s} & \alpha\_s^{\text{max}} \ge \alpha\_s > \alpha\_s^{\text{min}}\\ 0 & \alpha\_s \le \alpha\_s^{\text{min}} \end{array} \tag{25}$$

$$m = \begin{cases} \frac{\sqrt{3}}{2} \sin(\delta) & \nabla \cdot \overrightarrow{\boldsymbol{v}}\_s \ge 0 \\\ 1.03 & \nabla \cdot \overrightarrow{\boldsymbol{v}}\_s < 0 \end{cases} \tag{26}$$

where *Pc* is the solid critical pressure, Pa; *Pf* is the solid friction pressure, Pa. δ is the internal frictional angle of the solid phase; αmax *s* is the maximum packing concentration of the solid phase, and αmin *s* is the minimum frictional solid concentration, dimensionless.

The generation term of the turbulent kinetic energy is given by

$$G\_{\mathbf{k},\mathbf{l}} = \varepsilon\_{\mathbf{l}\dagger} \mu\_{\mathbf{l}} \left( \overrightarrow{\nabla \upsilon\_{\mathbf{l}}} + \overrightarrow{\nabla \upsilon\_{\mathbf{l}}}^{\mathbf{t}} \right) \colon \overrightarrow{\nabla \upsilon\_{\mathbf{l}}}.\tag{27}$$

Collisional dissipation term of the granular energy is given by

$$J\_s = \frac{48}{\sqrt{\pi}} \eta (1 - \eta) \frac{\rho\_s \alpha\_s^2 g\_0}{d\_s} \Theta\_s^{3/2}. \tag{28}$$

The interaction model of turbulence [31] is given by:

$$
\Pi\_{\Theta} = -3\beta\Theta\_{\sf s} + 81 \frac{\alpha\_{\sf s} \mu\_1^2 \left| \overrightarrow{\boldsymbol{\upsilon}}\_1 - \overrightarrow{\boldsymbol{\upsilon}}\_s \right|^2}{g\_0 d\_s^3 \rho\_s \sqrt{\pi \Theta\_s}}.\tag{29}
$$

In this study, the parameters used in upper equations are set as: *e* = 0.9, *c* = 1.6, *Fr* = 0.05, *r* = 2 *s* = 5 αmax *s* = 0.63, αmin *s* = 0.5, δ = π/6.

## *2.3. Boundary Conditions*

The values of the liquid phase velocity *vl*,*in* the solid phase velocity *vs*,*in*, and the volume concentration were given as the inlet conditions. The pressure value of 0 Pa was set as the outlet conditions.

The nonslip condition was used for the liquid phase, and Johnson-Jackson model [20] was used for the solid phase.

$$\rho \left(\mu\_s + \mu\_f\right) \frac{\partial v\_{sl}}{\partial \mathbf{x}} = -\frac{q \pi \rho\_s \alpha\_s g\_0}{2\sqrt{3} \alpha\_s^{\text{max}} s!} \tag{30}$$

$$\kappa\_{\rm s} \frac{\partial \Theta\_{\rm s}}{\partial \mathbf{x}} = \frac{\phi \pi \mathbf{v}\_{\rm s1} \mathbf{p}\_{\rm s} \alpha\_{\rm s} \mathbf{g}\_{0} \sqrt{\Theta\_{\rm s}}}{2 \sqrt{3} \alpha\_{\rm s}^{\rm max}} - \frac{\sqrt{3} \pi \mathbf{p}\_{\rm s} \alpha\_{\rm s} \mathbf{g}\_{0} \left(1 - \mathbf{e}\_{\rm w}^{2}\right) \sqrt{\Theta\_{\rm s}}}{4 \alpha\_{\rm s}^{\rm max}} \Theta\_{\rm s} \tag{31}$$

where *vs*,*l* is the slip velocity between the solid particles and the wall, m/s; φ is the specularity coefficient at the wall, dimensionless; *ew* is the particle–wall restitution coefficient, dimensionless. The parameters used in the upper equations are as follows: φ = 0.001, *ew* = 0.9.

The use of wall function allows the computation of turbulent flows without the need to resolve the turbulent boundary layer. The boundary layer for liquid phase [31] is defined as

$$\frac{\partial \upsilon\_l}{\partial \mathbf{x}} = -\frac{\rho\_l \kappa\_v \mathbf{C}\_{\mu}^{1/4} \sqrt{\mathbf{k}}}{\left(\mu\_l + \mu\_l^t\right) \ln(E \mathbf{x}^\*)} \upsilon\_l \tag{32}$$

$$\mathbf{x}^\* = \frac{\rho\_l \mathbf{C}\_{\mu}^{1/4}}{\mu\_l} \frac{\sqrt{k} \Delta \mathbf{x} / 2}{\partial \mathbf{x}}.\frac{\partial \mathbf{k}}{\partial \mathbf{x}} = 0 \,\frac{\partial \varepsilon}{\partial \mathbf{x}} = 0\tag{33}$$

The production *Pk* and dissipation *Dk* of *k* as well as ε at the fluid cells adjacent to the wall, is given by:

$$P\_k = \alpha\_l \rho\_l \sqrt{\mathbf{C}\_\mu k} \frac{v\_l}{\Delta x / 2 \ln(E x^\*)} + \beta k\_{12} \tag{34}$$

$$D\_k = \alpha\_l \rho\_l \varepsilon \tag{35}$$

where *E* is a Constant in wall function formulation, 9.81; κ*v* is Von Karmen constant of value, 0.42; Δ*x* is width of computational cell next to the wall, m.

#### *2.4. Geometric Model and Solution Strategy*

As shown in Figure 2, the length, width, and height of the slot were 4 m, 0.01 m, and 0.4 m, respectively. Three types of inlet positions of the top inlet, the middle inlet, and the bottom inlet were chosen to discuss the influence of the perforation positions in vertical hydraulic fractures. Eight outlets with a height of 0.02 m and a width of 0.01 m are uniformly distributed at the right wall.

**Figure 2.** Schematic diagram of slot geometry.

In order to verify the independence of the mesh, Figure 3 shows the proppant placement profile with different kinds of mesh conditions at T = 15 s. Obviously, mesh sizes will affect the final calculation results. But when the grid size is changed from 200 × 40 × 4 to 200 × 40 × 6 in the fracture length, height, and width directions, the proppant placement profile will no longer change. In order to obtain more detailed flow field information, we chose 200 × 40 × 6 as our computing grid in this study.

**Figure 3.** Numerical simulation result of proppant placement profile under different mesh conditions at T = 15 s.

Multiphase solver MFIX was used for the simulation. MFIX (multiphase flow with interphase exchanges) is an open-source code used for multiphase computational fluid dynamics simulations that were developed at National Energy Technology Laboratory (NETL, Pittsburgh, United States) [14,21,24,26]. The governing equations were discretized with the finite volume method and solved on a staggered grid. A second-order super bee discretization scheme was used for all variables. The transient numerical solution was obtained within a residual tolerance of less than 10−4. It takes about 74 h to complete the calculation of a single working condition in 60 s, with i7-7700 cpu, calling 4 threads. The plane at the center along the slot width was used to display the simulation results.

#### **3. Results and Discussions**

#### *3.1. Verification of Simulation on Proppant Transport*

To verify the Eulerian two-phase simulation on proppant transport, the proppant transport experiments on laboratory scale were conducted. The experiment was conducted in a fracture with a length, width, and height of 4 × 0.006 × 0.6 m, which comprised two pieces of parallel transparent plexiglass (shown in Figure 4). The proppant-water mixture was injected into a slot from perforation to mimic an artificial fracture. The parameters for experiment and simulation were set as the same, and the average Reynolds number is 3000 for both cases. The specific details were shown in Table 1.

**Figure 4.** Proppant transport experimental system.


**Table 1.** Parameters set for experiment and simulation.

Figure 5 shows the comparison of the proppant placement structure between experiment and simulation in the whole transportation process. In the early stage, the experiment result showed the proppant settled near the entrance and formed two dunes (Figure 5a), which was coinciding with the simulation case (Figure 5e). As proppant continued to be injected, the proppant dune first grew vertically and then grew horizontally (Figure 5b,c), and these characteristics were also can be demonstrated through numerical simulation (Figure 5f,g). Finally, the whole process had reached the equilibrium state for experiment and simulation cases. The final placement of proppant obtained by numerical simulation was basically consistent with the experimental results (Figure 5d,h).

**Figure 5.** Experimental and numerical simulation results of proppant transport and settlement process. (**<sup>a</sup>**–**d**): Proppant placement in the experiment process; (**<sup>e</sup>**–**h**): Corresponding proppant placement in the simulation process.

Specifically, Figure 6 shows the comparison of the proppant placement pattern obtained from the experiment and the simulation at the final equilibrium state. Due to the fracture height were different, the relative height was used to compare the equilibrium height in the experiment and simulation, which was defined as the height of the proppant bank to the fracture height. The background picture

was obtained from the experiment result, while the white dot line was extracted from the proppant placement profile obtained from the simulation result. The result shows that the profile of the proppant bank for the simulation case is generally consistent with the experiment result, even though some small differences exist in the inlet and outlet place. The reason for this difference may be: (1). The proppant particle size used in the experiment is nonuniform (0.4–0.8 mm), while the proppant size for simulation is set as constant at 0.8 mm; (2). the fracture outlet in the experiment is a single outlet at the top of the wellbore, while six outlets were uniformly distributed along the wellbore in numerical simulations. Hence, the placement near the outlet of the two cases was slightly different from each other. Despite these differences, the Eulerian two-phase model can effectively simulate proppant dynamic behavior and static placement pattern under HRNCs, which proves that this model can be used to simulate and predict the proppant transport and settlement in the fracture.

**Figure 6.** Comparison of proppant placement pattern obtained from experiment and simulation at the final equilibrium state.

#### *3.2. Proppant Transport Process*

In this section, by choosing the same simulation parameters in Table 1, proppant transport characteristics and mechanisms under HRNCs were comprehensively studied by analyzing the solid phase velocity field. As we discussed above, the typical proppant transportation process was divided into four typical main stages; those are: Initial settlement stage, vertically grow stage, horizontally grow stage, and equilibrium state stage.

#### 3.2.1. Stage1: Initial Settlement Stage

In the initial settlement stage, as the water–proppant mixture was injected form the middle inlet, two opposite vortexes were generated at the upper and lower ends of the slurry jet, and the lower one was smaller because of the inhibition generated by the slurry jet, as shown in Figure 7a. The injected proppant tended to quickly settle to the bottom and formed a proppant dune with a distance from the left boundary due to the gravity effect and the energy dissipation generated by the vortex. We defined this dune as the primary dune. Meanwhile, the small part of the proppant was dragged by the lower vortex move to the opposite direction and settle near the left boundary. We defined the small dune generated here as the secondary dune. Due to the slurry flush, a none proppant placement area was created between the two dunes. In this stage, this process was governed by the slurry vortex and the settlement.

**Figure 7.** Four stages for proppant transport under high Reynolds number condition. (**a**): initial settlement Stage; (**b**) vertically grow Stage; (**c**) horizontally grow Stage; (**d**) equilibrium State Stage.

#### 3.2.2. Stage2: Vertically Grow Stage

Figure 7b shows the vertically grow stage. The flow region in the fracture was divided into three main regions; those are free clean fluid flow area in the upper side of the primary dune, transition flow area on the surface of the dune, and the immobile area beneath the surface. In this stage, the liquid phase velocity near the primary dune was too small to carry proppant to a distant place. Therefore, most of the injected particles settled quickly on the proppant dune, and only a small part of them was carried to the deep part of the fracture. With a continuous injection of slurry, more proppant settled on the primary and secondary dunes, and the flow path between the dune and fracture upper surface became narrower. Therefore, the fluid velocity became larger, which increased the fluidization of the settled particles and reduce the proppant settlement. Besides, the primary dune reached a specific height at 0.28 m and stopped growing vertically, and we defined the height as the primary equilibrium height (PEH). Additionally, the secondary dune grew gradually near the entrance. The lower vortex was getting weaker in this stage due to the collective effects caused by the boundary, secondary dune, and the inject slurry flow. Consequently, the lower vortex disappeared while the secondary dune stopped growing. The valley between the two dunes was still existed due to the flushing. In this stage, the proppant transport process was governed by the fluidization, settlement, and suspension.

#### 3.2.3. Stage3: Horizontally Grow Stage

Figure 7c shows the proppant bank development process in a horizontal direction. In this process, because of the presence of the narrow gap between the proppant bank and the upper boundary, water had enough speed to fluidize and suspend the proppant. Therefore, a grea<sup>t</sup> amount of proppant transported through the currently generated proppant bank rather than retarded on the surface. However, once the proppant passes through the channel between the primary dune and wall, the velocity decreased dramatically because of the increased flow cross-section. Besides, another vortex was generated just behind the dune. Due to these two reasons, water cannot carry proppant into the deep of the fracture. Proppant gradually accumulated in the front of the formed proppant dune until reaching the equilibrium height, in which the water can suspend the proppant again. The settlement-reach equilibrium height process will repeat again and again, and from a macro perspective, the proppant dune will gradually grow horizontally. In this stage, the proppant distribution is determined by the settlement, fluidization and saltation.

#### 3.2.4. Stage4: Equilibrium State Stage

As the proppant dune grew horizontally and reached the outlet boundary, the proppant dune remained unchanged and reached a comprehensive equilibrium height (CEH), as shown in Figure 7d. The CEH was defined as the average equilibrium height of the primary dune when the proppant bank was stable during the injection process. The suspension and saltation were the main mechanisms to control the transport.

Figure 8 shows the proppant dune growth trajectory from 5 s to 65 s, in which the characteristic of the four main stages during the injection can be obviously obtained. Once the primary dune was formed, the dune will gradually grow vertically, and the front settlement angel almost remains constant at 20◦ until reaching the primary equilibrium height (PEH) by the end of the second stage (at around T = 20 s). The third stage began at T = 25 s, in which the proppant mainly grows horizontally, and the settlement angle of the upper side becomes more steep form 40◦ to nearly 90◦. At T = 60 s, the proppant dune reaches the CEH.

**Figure 8.** Proppant dune growth trajectory from for the simulation of the control sample.

## *3.3. Parametric Study*

In this section, the parametric study for di fferent injection velocity, inject position, particle diameter, density, and concentration were conducted for a better understanding of the e ffects caused by these parameters on the proppant transport process under HRNCs. The parameter sets are listed in Table 2. The bold items are set as the control sample in the parametric study. In each section, besides the investigated parameter, other parameters were set as same as the control sample.


**Table 2.** Parameter sets for parametric study and primary data for simulation.

The bold items are the setting value for the control sample.

#### 3.3.1. Inlet Velocity Effect

Figure 9 shows the distribution of the liquid volume concentration during the transport process for different injection velocity. In this study, the injection velocity of fluid and proppant was set to be the same, and the velocity uniformly expressed as v. The parameters except velocity were set as the same as the control sample. The average Reynolds number are 2000–5000 for v = 2–5 m/s, respectively. Compared with the control sample, the lower velocity case for *v* = 2 m/s (Figure 9A) showed a closer distance from the primary dune to the entrance and more particle placement between the dunes, which was mainly because the vortex is weaker than that of the control sample. Meanwhile, the proppant dune grew rate for v = 2 m/s is a little bit slower than that of the control sample. On the contrary, as the injection velocity became faster, the vortex generated near the entrance became larger and stronger, resulting in a larger distance for primary dunes and disappeared secondary dune. As the velocity of the slurry increases from 2 m/s to 5 m/s, the non-propped area near the inlet increase by 5.3 times.

**Figure 9.** Volume concentration of liquid phase during the transport process for different injection velocity cases: (**A**) v = 2 m/s, (**B**) v = 3 m/s, (**C**) v = 4 m/s, and (**D**) v = 5 m/s.

The e fficiency of proppant placement is also an essential factor in evaluating the e ffect of proppant placement. Therefore, we defined the proppant occupation as a stationary proppant bed area over the total fracture area. The proppant occupation in equilibrium stages was defined as equilibrium proppant occupation (EPO), which indicated the maximum proppant accumulation in the fracture under the specific condition. Proppant occupation, together with the EPO and CEH, were used to characterize the proppant transport and settlement for each case.

Figure 10 shows the proppant occupation of di fferent time steps and comparison of equilibrium height for di fferent injection velocity cases. The case for v = 5 m/s reached the equilibrium stage at T = 22 s, which was the fastest among all instances. As velocity decreased, it takes longer for the proppant dune to reach equilibrium in the fracture. However, the v = 2 m/s cases yielded the highest CEH at 35.9 cm, compared with other cases, which were at 35.2 cm, 34.6 cm, and 34.1 cm, respectively. For the low-velocity instance, the gap between the primary dune needed to be narrower to obtain enough speed to balance the amount of settled and rolled-up proppant, resulting in higher CEH. Besides, due to the smaller (none) secondary dunes and the relatively shorter primary dune, it was observed that the EPO for large velocity is lower than those of slower cases during the transportation process, which are 80.5%, 74.6%, 66.4%, and 56.5% for the velocity from 2 to 5 m/s. The small EPO, especially caused by the no proppant placement neat the entrance, will cause worse conductivity in artificial fracture, which should be avoided in the field fracturing design.

**Figure 10.** Proppant occupation percentage in fracture change with time a comprehensive equilibrium height (CEH) for di fferent velocity cases.

#### 3.3.2. E ffects of Inlet Position

Figure 11 shows the distribution of the proppant placement in the transport process under di fferent inlet positions. According to Figure 11A, when the slurry was injected from the bottom, the distance between the initial settlement location and left wall was 1.1 m, which indicated that the flow velocity at this position was not able to overcome the frictional force between the proppant particles and the bottom wall. Since no flow vortex was formed upon the bottom wall, no proppant particles were situated in the region near the left wall. As a result, only the primary dune was generated in the fracture. Proppant particles injected later would ge<sup>t</sup> over the dune and settle at the rear of the dune. Finally, the proppant dune gradually grew until it reached the CEH.

**Figure 11.** Proppant placement during the transport process for different inlet position cases: (**A**) Inject from the bottom inlet; (**B**) inject from the middle inlet; and (**C**) inject from the top inlet.

Figure 11C shows the proppant transport process of the top inlet condition at different time points. When the mixture was injected into the slot, the proppant quickly settled to the bottom at the position about 0.65 m from the left wall. A big clockwise vortex was generated near the entrance and brought proppant back to the inlet side, forming a sizeable secondary dune. As the proppant was injected continuously, the heights of both dunes increase, and the trough between two dunes gradually filled with the proppant. Consequently, the primary and secondary dunes gradually became an integrated dune, and the CEH was obtained almost alone the whole length of the fracture.

According to Figure 12, it took almost the same time (55 s) to reach CEH for the three cases with different inlet positions. Injecting proppant from the top could produce the highest equilibrium height at 36.1 cm, compared with that of the middle case at 35.2 cm and the bottom situation at 34.3 cm. Besides, the proppant occupation curve indicated that the proppant injected through the top yielded the highest proppant placement percentage, and the EPO reaches 82.3% compared with the middle case at 74.6% and the bottom case at 67.8%. The result identified that the vortex at the front end of the fracture is the dominant factor for proppant placement near the entrance. The large vortex could carry the proppant back to the inlet place. The large non-propped area in the fracture may cause fracture closure after pumping, resulting in the permeability significantly decreased. Therefore, choosing the relatively high inlet positions will result in better patterns of the proppant placement and improve fracture conductivity.

**Figure 12.** The curve for proppant occupation percentage in fracture change with time and CEH for different inlet position cases.

#### 3.3.3. Effects of Particle Diameter

Figure 13 shows the effects of particle diameter on transport processes. The accumulation curve demonstrated that the larger particles settle more in the fracture at any time during the transport process because smaller particles were easier carried by the fluid. The smaller particle was carried into the far end of the fracture rather than accumulated during the transportation. For the same reason, it took a long time for smaller particles to reach CEH. However, the CEH difference is tiny among different cases, which indicated the proppant diameters had little effect on the primary dune for a relatively long time in the transport process. At last, EPO for 1 mm to 0.6 mm are 67.6%, 68.8%, 71.0%, 73.1%, and 74.6% respectively.

**Figure 13.** The curve for proppant occupation percentage in fracture change with time and comprehensive equilibrium H = height (CEH) for different particle diameter cases.

#### 3.3.4. Effects of Particle Density

Figure 14 shows the detailed data extracted from the simulation results for different particle density. The results indicated proppant density has a significant influence on the proppant settlement in the fracture. The curve shows that the lighter particle has the least amount of settlement in the fracture during the injection process. The reason for that phenomenon was that the lighter particles were more likely to be fluidized and dragged by the fluid to the further side of the fracture. All the cases almost reached the equilibrium height at the same time, around 45 s. The EPO was 60.39% when the density was 1800 kg/m<sup>3</sup> is, which was 21.3% smaller than that of 3200 kg/m3. For the CEH study, it was interesting to obtain a proppant density threshold for the particle dune growth. When the density is small, the equilibrium height increases rapidly as the proppant density increases, which changed from 32.8 cm in the 1800 kg/m<sup>3</sup> case to 35.8 cm in the 2500 kg/m<sup>3</sup> case. However, when the density was larger than 2500 kg/m3, the CEH almost unchanged.

**Figure 14.** The curve for proppant occupation percentage in fracture change with time and CEH for different inlet particle density.

#### 3.3.5. Effects of Solid-Phase Volume Concentration

Figure 15 shows the particle concentration effects on proppant transport. From the curves, higher proppant concentration resulted in a greater amount of settlement in the fracture at the same time. Besides, the lowest concentration (10%) case was the last one to reach the CEH at T = 76 s, which almost double that of 25% concentration This phenomenon was caused by two reasons; those are fewer injection particle amounts and easy fluidization for proppant when the concentration was low. Moreover, the CEH for lower concentration cases was also shorter, which may result from more liquid fluid flowing past over the top of the proppant bed or the higher ability of lower solid phase volume concentration to fluidize the proppant particles. In conclude, the lower concentration in the transport process lead to a more proppant place in the further part of the fracture, which was good for longer effective propped fractures. However, less proppant in the fracture may result in poor conductivity. The proppant particle concentration selection in the field should be considered comprehensively for both sides.

**Figure 15.** The curve for proppant occupation percentage in fracture change with time and CEH for different inlet position cases.

#### *3.4. Equilibrium Height Prediction Model*

The height of the equilibrium gap H, which was defined as the distance between the top of the proppant bed and the upper wall boundary, was used to represent the equilibrium height. J. Wang et al. [13] established a Bi-power law correlation for the equilibrium gap of the proppant bed, and it was given by:

$$\frac{\mathcal{H}\_1}{\mathcal{W}} = \left[ -2.3 \times 10^{-4} \ln(R\_G) + 2.92 \times 10^{-3} \right] \mathcal{R}\_f^{1.2-1.26 \times 10^{-3} \lambda^{-0.428} \times [15.2 - \ln(R\_G)]} R\_p^{[-0.0172 \ln(R\_G) - 0.12]}.\tag{36}$$

The gravity Reynolds number *RG* is defined as:

$$R\_{\mathbb{G}} = \rho\_l (\rho\_s - \rho\_l) \lg d\_s^2 / \mu\_l^2. \tag{37}$$

The gravity Reynolds number for the fluid λ is defined as:

$$
\lambda = (\mu\_l / \rho\_l) / \left(\mathcal{W}^{3/2} \sqrt{\mathfrak{g}}\right). \tag{38}
$$

The fluid Reynolds number *Rf* is shown as:

$$R\_f = \left(\rho\_l \overline{v} \mathcal{W}\right) / \mu\_l. \tag{39}$$

The proppant Reynolds number *Rp* is defined as:

$$\mathcal{R}\_p = \left(\rho\_p \overline{v} \mathcal{W}\right) / \mu\_l \tag{40}$$

where H1 is the height of equilibrium gap, ml; W is the fracture width, m; *H*/*W* is gap height over fracture width, dimensionless;

Along the altitude direction of the fracture, the flow field can be divided into three layers. The bottom of the proppant bank is an immobile layer in which the velocity of the proppant particles is approaching zero. The middle layer is a mobile zone that is above the stationary bed, in which the proppant particles move by sliding and rolling or a combination of both. The top layer is a clean liquid zone. As the proppant concentration increases gradually from the fresh liquid layer to the mobile layer and then to the immobile layer, the point with the proppant concentration equal to 0.5 was chosen to be the top of the stationary bed.

Figure 16 shows the comparison of dune height between the simulation result and Wang's model under different conditions. As the parameters change, they had the same trend. However, compared with Wang's model, the equilibrium height obtained by the simulation was always higher. The mean deviation for the case of density, diameter, velocity, and solid concentration was 45.8%, 36.5%, 57.7%, and 59.4%, respectively, which indicated it is difficult to accurately predict the equilibrium height by using this model.

**Figure 16.** Comparison of the equilibrium dune height between simulation result and Wang's model. (**a**) Comparison of particle density; (**b**) comparison of particle diameters; (**c**) comparison of injection velocity; and (**d**) comparison of particle concentration.

In this paper, some modification about Wang's model was conducted based on the simulation result, to accurately predict the equilibrium height of the proppant dune in the fracture under HRNCs. The new prediction model is shown as:

$$\frac{\mathcal{H}\_1}{\mathcal{W}} = \left[ -5.76 \times 10^{-3} \ln(R\_{\mathcal{G}}) + 1.87 \times 10^{-3} \right] \mathcal{R}\_f^{1.2-1.008 \times 10^{-3} \lambda^{-0.428} \times [15.2 - \ln(R\_{\mathcal{G}})]} R\_p^{[-0.04 \ln(R\_{\mathcal{G}}) - 0.12]}.\tag{41}$$

Figure 17 shows the relationship between the new prediction model and simulation results. It can be seen from the comparison that this prediction model can fit the simulation result very well. Through calculation, the mean deviation is only 3.8%. The applicable range for this prediction model is the case where the Reynolds number is 2000–5000. This correlation can be used to quickly predict proppant placement in a fracturing simulator, which can greatly improve the simulation accuracy.

**Figure 17.** Comparison of the equilibrium dune height between simulation result and new prediction model. (**a**) Comparison of particle density; (**b**) comparison of particle diameters; (**c**) comparison of injection velocity; and (**d**) comparison of particle concentration.
