## Abstract

We implement proppant transport in a three-dimensional hydraulic fracturing simulator, including proppant settlement due to gravity, tip screen-out, and fracture closure. Constitutive equations are used that account for processes that can cause the ﬂowing fraction of proppant to be different from the volumetric fraction of proppant. The constitutive equations capture the transition from Poiseuille ﬂow to Darcy ﬂow as the slurry transitions from dilute mixture to packed bed. We introduce new constitutive equations that allow the simulator to seamlessly describe the process of fracture closure, including a nonlinear joint closure law expressing fracture compliance and roughness and accounting for the effect of proppant accumulation into a packed layer between the fracture walls.

#### Sogo Shiozawa^{a,n}, Mark McClure^{a,b}

^{a}Department of Petroleum and Geosystems Engineering, The University of Texas at Austin, Austin, TX, USA. ^{b}McClure Geomechanics LLC, Palo Alto, CA, USA

*Received 6 October 2015, Received in revised form 16 December 2015, Accepted 5 January 2016.*

© 2016 The Authors.

We perform sensitivity analysis simulations to investigate the effect of ﬂuid viscosity, proppant density, proppant size, and formation permeability. The simulations conﬁrm that tip screen-out can limit fracture length, cause proppant banking, and increase injection pressure. Sensitivity analysis indicates that reasonably accurate results can be achieved without excessive mesh reﬁnement. We also perform a simulation of hydraulic fracture propagation through a complex natural fracture network. In this simulation, proppant tends to accumulate at the intersections between natural and hydraulic fractures. Overall, the results suggest that in very low permeability formations, proppant settling is a major problem for proppant placement because proppant tends to gravitationally settle before fracture closure can occur. Because leakoff is so slow, proppant immobilization through bridging is critical for vertical proppant placement. Bridging can occur at aperture approximately three times greater than particle diameter, which will occur much sooner after shut-in than full mechanical closure. Even though larger diameter proppant settles more rapidly, it may lead to better proppant placement because it will bridge sooner, at a larger fracture aperture. These results also suggest that it is critical to optimize injection schedule in order to avoid tip screen-out, which leads to a shorter, wider fracture in which bridging is less likely to occur. Our modeling approach can be used practically for optimization of proppant placement through selection of ﬂuid properties, proppant properties, and injection schedule.

## 1. Introduction

Hydraulic fracturing is performed by injecting ﬂuid into the subsurface at high rate and pressure, opening and propagating fractures through the formation. In the majority of fracturing treatments, particulate matter called proppant is pumped in a slurry with the injection ﬂuid. After injection is stopped, ﬂuid pressure decreases, and the fractures close. The proppant holds the fractures open and increases their ability to conduct ﬂuid after closure.

Several approaches have been used for numerical simulation of ﬂuid–solid two-phase systems, such as proppant slurry. The two most common frameworks are Eulerian–Eulerian and Eulerian–Lagrangian (Hu et al., 2001; Zhang and Chen, 2007; Tsai et al.,2012). In the Eulerian–Eulerian technique, the particles and ﬂuid are both treated with an Eulerian framework.

Each component is governed by conservation equations in stationary control volumes (Clifton and Wang, 1988; Ouyang et al., 1997; Mobbs and Hammond, 2001; Adachi et al., 2007; Weng et al., 2011; Dontsov and Peirce, 2015). In the Eulerian–Lagrangian technique, proppant transport is described with a Lagrangian framework, which tracks the locations of individual particles or groups of particles, and ﬂuid ﬂow is described with an Eulerian framework (Tsai et al., 2012; Tomac and Gutierrez, 2015). For describing slurry ﬂow, it is necessary to calculate an effective ﬂuid viscosity (Adachi et al., 2007).

The earliest major contribution on this topic was the theory of dilute suspensions of particles (Einstein, 1905). For concentrated suspensions of particles, one of the simplest expressions was introduced by Mooney (1951). For the modeling of proppant transport, an expression similar to the Krieger–Dougherty equation (Krieger and Dougherty, 1959) is usually used (Adachi et al., 2007). In this study, we follow the method of Dontsov and Peirce (2014), who used the constitutive model introduced by Boyer et al. (2011), which is described below.

The slip velocity vector expresses the difference in average velocity between the particles and ﬂuid. There is a tendency for transverse particle migration away from the fracture walls, where shear stress is maximum, to the center of the ﬂow channel, where shear stress is lowest. This phenomenon causes higher proppant concentration at the center of channel, where ﬂuid velocity is highest (Constien et al., 2000). Some models assume that proppant distribution is uniform across the aperture, and so the velocity difference between ﬂuid and proppant is caused only by gravity (Adachi et al., 2007). Other models account for proppant migration away from the fracture walls to the center of the ﬂow channel. Mobbs and Hammond (2001) performed simulations of proppant transport taking into account the migration effect with an assumed proppant distribution across the aperture. Boronin and Osiptsov (2014) performed a similar analysis with a different assumed particle distribution and achieved good agreement with experiment results. Eskin and Miller (2008) developed a model accounting for micro-level particle dynamics from kinetic theory. Dontsov and Peirce (2014) derived an expression for the dis- tribution of proppant velocity across the fracture aperture as a function of average proppant concentration using the empirical constitutive model introduced by Boyer et al. (2011). Taking into account the slip velocity and performing a boundary layer calculation, Dontsov and Peirce (2014) described the transition from Poiseulle to Darcy ﬂow that occurs as proppant concentration increases. Using their model, Dontsov and Peirce (2015) performed simulations of proppant transport with Khristianovich–Geertsma–de Klerk (KGD) and pseudo-three dimensional (P3D) hydraulic fracture models using the Carter leakoff model (Howard and Fast,1957).

In this study, using the approach introduced by Dontsov and Peirce (2014), we perform simulations of proppant transport in a fully three-dimensional hydraulic simulator, CFRAC (McClure and Horne, 2013; McClure et al., 2016). We perform simulations of the entire injection and post-injection period, simulating fracture propagation and closure. We extend the framework of Dontsov and Peirce (2014) to describe the process of fracture closure after the end of injection. At zero or low proppant concentration, the algorithm allows ﬂuid storage and conductivity of closed fractures to be described by a nonlinear joint closure law related to the stiffness of the asperities in the fracture walls. At high proppant concentration, ﬂuid storage and conductivity of closed fractures are primarily controlled by the properties of the packed proppant bed. The algorithm seamlessly handles the transition between these two end members during fracture closure, for all possible values of proppant concentration.

This paper provides the details of our model, including governing equations, the method for handling fracture closure with proppant, and the method of solution. Simulations of injection into a single planar fracture are provided, with key parameters varied for sensitivity analysis. The results demonstrate that the model is capable of describing the tip screen-out (TSO) process, in which proppant accumulates at the tip and slows or prevents further fracture propagation. Finally, we describe a proppant transport simulation performed in a large, complex discrete fracture network model.

## 2. Methodology

### 2.1. Model setup and assumptions

In this study, a fully three-dimensional hydraulic fracturing simulator, CFRAC, is extended to perform proppant ﬂow calculations, including: pressure-driven convection, gravitational settling, and fracture closure. The details of CFRAC have been described in previous publications (McClure and Horne, 2013; McClure et al., 2016). The governing equations, constitutive equations, and methods of solution are summarized in the following sections, along with the modiﬁcations to the code that are performed for this study. The simulator couples unsteady state mass balance equations for ﬂuid and proppant with mechanical calculations for the stresses induced by fracture opening and sliding. Fully implicit timestepping is used, guaranteeing numerical stability.

The simulations are fully three-dimensional, and so the fractures are meshed in both the vertical and horizontal directions. CFRAC can simulate ﬂow in individual hydraulic fractures or in large discrete fracture networks involving both hydraulic fractures and natural fractures. In all cases, each individual fracture is individually discretized, and ﬂuid ﬂow and deformation of each fracture is calculated fully numerically, without upscaling to a continuum approximation. For example, Section 3.5 shows a simulation of hydraulic fracture propagation through a naturally fractured reservoir.

**Fig. 1. **Fracture coordinate system used in this study.

The fractures are not meshed across their aperture, though cross-aperture variation in ﬂow velocity is implicitly considered by the constitutive equations. All fractures in the model are assumed to be vertical. The coordinate system is shown in Fig. 1.

Both ﬂuid and proppant particles are described with an Eulerian framework (Eulerian–Eulerian approach), with ﬂuid and proppant mass conservation equations solved in stationary control volumes with the ﬁnite volume method (Karimi-Fard et al., 2004). The mass balance equations are solved simultaneously with mechanical equilibrium equations using iterative coupling. The primary variables at each element are (1) ﬂuid pressure in fracture, (2) fracture aperture, (3) fracture sliding displacement, and (4) volume fraction of proppant. The simulations are isothermal.

It is assumed that the proppant consists of non-colloidal spherical particles and all of the particles have the same size. Proppant particles are assumed to be incompressible in the ﬂow calculation, but when the proppant packs into a porous bed, the compressibility of the porosity is taken into account (Section 2.3). The carrier ﬂuid is slightly compressible and Newtonian. We use constitutive relations for proppant and ﬂuid ﬂow in the fracture that were developed by Dontsov and Peirce (2014). Their relations assume laminar ﬂuid ﬂow and negligible Brownian motion of the proppant, implying:

where P_{e} and R_{e} are Péclet and Reynolds numbers, respectively, μ is ﬂuid viscosity, γ is the shear rate of a simple-shear ﬂow, d is proppant diameter, k is the Boltzman constant, T is temperature, and ρ_{f} is ﬂuid density (Morris and Boulay, 1999).

### 2.2. Flow equations

The unsteady state ﬂuid mass balance equation and ﬂuid ﬂux are given as

where E is the aperture, ρ_{f} is the ﬂuid density, φ is the volume fraction of proppant, ∇ is the gradient operator, qf,flux is ﬂuid mass ﬂux (mass ﬂow rate per cross-sectional area), μ is ﬂuid viscosity, q_{leak} is ﬂuid mass leakoff rate per fracture surface area, k is fracture permeability, P is ﬂuid pressure, s_{f} is a ﬂuid source term, Q_{s} and Q_{p} are dimensionless functions of proppant concentration and aperture, and χ is a blocking function (described later in this section). In some versions of CFRAC, a distinction is made between the void aperture E (ﬂuid volume per surface area) and the hydraulic aperture (used for calculating the fracture transmissivity), but in the present work, these apertures are assumed the same.

**Fig. 2. **The dimensionless functions Q_{s}, Q_{p} and G_{p} introduced by Dontsov and Peirce (2014) as functions of normalized proppant concentration φ for different values of the parameter E/d. Normalized proppant concentration is equal to the volumetric fraction of proppant divided by φm. The black, blue, and red lines represent cases with E/d = {50, 5, 3}, respectively. (For interpretation of the references to color in this ﬁgure caption, the reader is referred to the web version of this paper.)

The model calculates leakoff of ﬂuid from the fractures into the surrounding rock using a one-dimensional leakoff model (Vinsome and Westerveld, 1980). The model solves the one-dimensional diffusivity equation, which implies single phase, single component Darcy ﬂow in a porous media with constant ﬂuid viscosity and constant total compressibility. This is a simpliﬁcation because leakoff is actually a multiphase, multicomponent process. Despite (or perhaps because of) the complexity of the leakoff process, in hydraulic fracturing simulators, leakoff is nearly always calculated with a highly simpliﬁed model. The most common leakoff model is Carter leakoff (Howard and Fast, 1957), which is derived based on the same simplifying assumptions used by the Vinsome and Westerveld (1980) model, as described above.

However, Carter leakoff assumes constant ﬂuid pressure in the fracture, which can lead to unacceptable inaccuracy in problems involving fracture closure after shut-in, when ﬂuid pressure decreases signiﬁcantly. The Vinsome and Westerveld (1980) model solves the one-dimensional diffusivity equation using a semi-analytical technique that accounts for changing ﬂuid pressure in the fracture over time. The poroelastic stresses induced by the ﬂow of ﬂuid into the matrix and the poroelastic stress changes induced in the matrix due to fracture deformation are neglected, as is conventional in hydraulic fracturing simulators. Wallace et al. (2014) found that these stresses are typically negligible.

Mack and Warpinski (2000) describe why the use of Carter leak off is considered acceptable. Carter leakoff uses a leakoff coefﬁcient, an effective parameter that lumps together the effects of different complex processes, and which can be estimated from fracture calibration tests (Nolte, 1979; Mayerhofer and Economides, 1993). In CFRAC, to account for the complexity of the leakoff process, the leakoff rate can be calculated using a special leakoff viscosity, approximately equal to the viscosity of the clean ﬂuid, and not equal to the viscosity of the ﬂuid in the fracture, which is elevated due to the inclusion of gelling and cross linking agents (Mack and Warpinski, 2000). These large molecular weight viscosifying agents are too large to enter into the tiny pores of shale formations. The use of a lower leakoff viscosity in the calculations increases the overall leakoff rate.

The normalized proppant concentration is expressed as

where φ_{m} is the maximum possible volume fraction of proppant, which is set to be 0.585 in this study. When the volume fraction of proppant reaches its maximum, proppant forms immobile cluster with porosity of 1−φ_{m}. The normalized proppant concentration ranges from 0 to 1. The unsteady state proppant mass balance equation and proppant ﬂux are given as

respectively, where ρ_{p} is proppant density, q_{p},flux is proppant mass ﬂux, s_{p} is a proppant source term, d is a proppant diameter, g is the gravitational acceleration, and G_{p} is a dimensionless function of proppant concentration and aperture.

Fracture permeability is deﬁned as

Transmissivity is deﬁned as a product of permeability and aperture, which yields the cubic law for fracture transmissivity (Witherspoon et al., 1980):

Arithmetic averaging is used for calculating the transmissivity for ﬂow between elements. Harmonic averaging is not used because it can lead to strong mesh dependency in cases with strong transmissivity contrast and coupling between ﬂuid pressure and transmissivity. Q_{s}, Q_{p}, and G_{p} are the dimensionless functions introduced by Dontsov and Peirce (2014) and expressed as

where Q_{s}, Q_{p}, G_{s}, and G_{p} are dimensionless functions of normalized proppant concentration and D_{p} is a constant related to packed particles (=8(1-φ_{m})^{β}/(3φ_{m})). The expression (1-φ_{m})^{β} comes from the normalized settling rate, which can be hindered due to interaction between particles (Morris and Boulay, 1999). In this study, φm and β are chosen to be 0.585 and 4.1, respectively, following Dontsov and Peirce (2014). These functions can be calculated numerically. For this study, tables providing the function values were provided by Egor Dontsov (personal communication).

Qs captures the transition from Poiseuille to Darcy ﬂuid ﬂow as φ increases from 0 to 1.0. The ﬁrst term of Eq. (9) represents the inverse of the effective viscosity of slurry. The viscosity of slurry is higher than the viscosity of pure ﬂuid because of interactions between particles and interactions between particle and the ﬂuid. The second term describes Darcy ﬂow in porous media, and its effect is only signiﬁcant when the normalized proppant concentration is close to 1. Qp and Gp control the ﬂowing volume of proppant due to convection and gravitational settling, respectively.

Q_{s} simpliﬁes to the cubic law when proppant concentration is low and simpliﬁes to Darcy's law when proppant concentration is high. Q_{p} and G_{p} describe pressure-driven proppant convection and gravity settling, respectively. At low concentration, the value of Q_{p} is such that the ﬂowing fraction of proppant is about two times larger than the volume fraction of proppant. In the nomenclature of multiphase ﬂow theory, this is equivalent to saying that the fractional ﬂow is twice as large as the saturation. Both Q_{p} and G_{p} become zero when proppant volume fraction reaches the maximum allowed value, indicating that proppant cannot ﬂow because it has formed a packed, immobile bed. At dilute proppant concentration, the settling velocity is approximately given by Eq. (31) in Section 3.2 (Dontsov and Peirce, 2015).

The blocking function is expressed as

The blocking function accounts for proppant bridging. The function allows proppant to ﬂow into or out of a fracture element only when a fracture is open and its aperture is greater than Nmin times the proppant diameter. To avoid a discontinuity in the ﬂuid ﬂow derivatives, the function linearly increases from 0 to 1 as aperture changes from N_{min} d to N_{max} d. Nmin and Nmax are chosen to be 3 and 4, respectively, in this study. The value of χ is determined based on Eq. (12) at the beginning of a time step and kept constant during a single timestep, as shown in Fig. 3.

**Fig. 3.** Summary of the iterative coupling approach during a timestep.

In order to obtain a numerically stable scheme for the proppant transport, it is necessary to choose either upwinding (calculating ﬂuxes using values from the element ﬂuid is ﬂowing out of) or downwinding (using values from the element ﬂuid is ﬂowing into). The winding for each connection between adjacent fracture elements is determined by the sign of the derivative of proppant ﬂux with respective to the product of aperture and normalized proppant concentration (Dontsov and Peirce, 2015). If the signs of both derivatives are positive, the ﬂux is upwinded. If the sign of both derivatives is negative, the ﬂux is downwinded. If the signs are different, a shock velocity is calculated and used to determine the wind direction. If the sign of shock velocity is positive, up-winding is used. If the sign of shock velocity is negative, down-winding is used. In this case, the shock velocity is deﬁned as

where i and j are two adjacent elements and located upstream and downstream, respectively. For simplicity, the wind direction is determined for each connection at the beginning of each timestep and ﬁxed during the timestep. The treatment of winding is critical to implement properly. If upwinding/downwinding is not determined appropriately, the model can yield physically unrealistic results, such as dimensionless proppant concentration greater than 1.0.

### 2.3. Fracture aperture

We deﬁne an “open” fracture as a fracture where the ﬂuid pressure has reached the normal stress, the fracture walls have come out of contact, and compressive stress is not transmitted across the fracture by the proppant. The aperture (volume of ﬂuid and proppant stored in the fracture per surface area) of an open element is decomposed into three components: E_{0}, representing the contribution of the roughness of the fracture walls, E_{p}, the hypothetical ﬂuid and proppant volume per surface area if ﬂuid was drained from the fracture until φ equaled 1.0, and Eopen, the additional separation of the fracture walls. The aperture of an open element is

With these deﬁnitions in place, it is possible to consistently deﬁne aperture through fracture opening and closure for any proppant concentration, either very high or very low.

A fracture is deﬁned as open if E_{open} is greater than 0. Some amount of proppant could become lodged within the “roughness” dominated portion of the aperture, represented by E_{0}. If the fracture contains less than that maximum capacity, then E_{p} is set to zero. This is calculated by evaluating the following:

In this case, the aperture at the transition from open to closed is E_{0}, and the volumetric fraction of proppant at closure is equal to the maximum volume fraction of proppant φ_{m} or less. Eq. (15) relies the simplifying assumption that the proppant would be able to ﬁll the roughness-generated aperture of the fracture up to a volumetric fraction of exactly φ_{m}.

If a fracture element contains a sufﬁciently large amount of proppant, there will be a layer of proppant separating the fracture walls when the fracture closes, and the aperture at closure will be greater than E_{0}. In this case, the aperture at the transition between open and closed is equal to E_{0}+Ep. Along with Eq. (15), Ep is updated according to

Eqs. (15) and (16) are applied explicitly by updating Ep only at the beginning of every timestep. When Ep is updated, the value of Eopen is changed by the opposite amount, ensuring that E remains constant. Because E_{p} is updated only at the beginning of the timesteps, the conditions in Eqs. (15) and (16) can sometimes be slightly violated at the end of a timestep, but this has a very minor effect.

The aperture of a closed element is deﬁned as

where σ_{n}′ is the effective normal stress, σ_{n,ref} is a constant deﬁned as the effective normal stresses required to cause a 90% reduction in aperture, and c_{p} is the compressibility of the part of the aperture ﬁlled with proppant. The effective normal stress σ_{n}′ is deﬁned as Jaeger et al. (2007):

where σ_{n} is the normal stress. The E_{0} term in Eq. (17) represents the natural compressibility of the fracture due to deformation of the asperities in the fracture wall. The E_{p} term in Eq. (17) represents the compressibility of the porosity of the proppant bed. Eq. (17) is somewhat ad hoc, but it is physically plausible. Because the proppant is assumed incompressible, but the aperture of a closed fracture element is assumed compressible, the normalized volume fraction of proppant in a closed element can exceed 1.0. This approximates the effect of proppant embedment in the fracture walls.

The E_{0} term in Eq. (17) was ﬁrst used by Willis-Richards et al. (1996) to describe joint closure, based on the work of Barton et al. (1985). In their work, the value of σn,ref is a constant, considered a property of the fracture. However, if proppant occupies the fracture, we would expect that the fracture would become stiffer, due to the presence of a stiff, bridged, immobile solid phase wedged between the fracture walls. Therefore, we chose to make σ_{n,ref} a function of normalized proppant concentration at closure. We use a simple expression relating these values to proppant concentration:

where σ_{n,ref}, max and σ_{n,ref},min are the maximum and minimum values of the effective normal stresses required to cause a 90% reduction in aperture. If Eq. (19) is not used to increase σ_{n,ref} as a function of φ, physically unrealistic behaviors can result, such as aperture lower than the aperture value required to contain the proppant grains in the element.

Because of the bridging of the proppant between the fracture walls, proppant is not permitted to ﬂow into or out of a closed fracture element. For preexisting fractures, E_{0} is treated as a constant. However, a special algorithm must be used to deﬁne E_{0} for hydraulic fracture elements (McClure, 2014). When a hydraulic fracture element is initiated, it is given the very small initial value of E_{0} of 0.1 microns.

The newly initiated element is ﬁlled with ﬂuid, and water is not taken from an adjacent element to compensate. This does not strictly conserve mass, but the mass balance error is slight because the initial aperture is so small. The element cannot be initialized with a larger E_{0} because it would result in a more signiﬁcant mass balance error. But the residual aperture of a typical joint, E_{0}, is on the order of 100s of microns. Therefore, an algorithm is needed to allow E_{0} to increase as the element progressively opens. The algorithm sets E_{0} to be equal to 90% of E up until a maximum value, E_{hf max,resid}. E_{0} can only increase when the element is open. When E_{0} is updated, Eopen is decreased by an equal amount in order to maintain constant E. E_{0} is not permitted to decrease. This algorithm mimics the natural process of fracture roughness development as a fracture forms and opens. The algorithm in Eqs. (15) and (16) for updating E_{p} is applied after the updating of E_{0}.

### 2.4. Wellbore calculations

Injection is performed with a speciﬁed maximum rate and pressure. The simulator enforces both conditions simultaneously, maintaining either constant rate or constant pressure, depending on which boundary condition will permit both the maximum pressure and rate conditions to be enforced. The ﬂuid ﬂow calculation does not include frictional pressure drop in the wellbore.

When proppant is injected at the surface, it does not enter the formation until all of the ﬂuid in the well has ﬁrst ﬂowed into the formation. This process is included in the simulator. The wellbore is meshed into a series of elements of constant volume. The ﬂuid velocity is assumed constant along the entire well and the ﬂowing volume fraction of proppant is assumed equal to the actual volume fraction of proppant. The mass balance equations for proppant are solved in each wellbore element simultaneously with the mass balance equations for proppant in the fracture elements (Eq. (5)).

### 2.5. Mechanical calculations

The displacement discontinuity method, a boundary element method, is used for the calculation of stresses induced by fracture deformation, assuming an elastically homogeneous and isotropic formation, linear elastic deformation, and small strain (Okada, 1992). The Okada (1992) method simultaneously satisﬁes the equations of quasistatic stress equilibrium, the compatibility equations, and the constitutive equations of linear elasticity, providing solutions that are convergent with mesh reﬁnement to exact analytical solutions from classical continuum mechanics.

The Okada (1992) method assumes that the fractures are embedded in a semi-inﬁnite half-space, enforcing a traction free boundary condition at the edge of the half-space (representing the Earths surface) and enforcing the boundary condition that induced stresses and displacement go to zero as distance goes to inﬁnity within the semi-inﬁnite domain.

To satisfy force equilibrium, the effective normal stress of open elements is enforced to be equal to zero (Crouch et al., 1983):

To satisfy force equilibrium, the effective normal stress of open elements is enforced to be equal to zero (Crouch et al., 1983):

The shear stress on open elements is also enforced to be equal to zero:

where τ is shear stress, v_{s} is the sliding velocity, η is the radiation damping coefﬁcient (Rice, 1993). For closed elements, the Coulomb failure criterion with a radiation damping term is enforced (Jaeger et al., 2007):

where μ_{f} is the coefﬁcient of friction and S_{0} is fracture cohesion. If the left-hand side of Eq. (22) is less than the right-hand side, the fracture sliding velocity is assumed to be zero because the shear stress is less than the frictional resistance to slip. If the fracture is sliding, then equality is enforced in either Eq. (21) or (22), depending on whether the element is open or closed.

Fracture propagation is predicted using linear elastic fracture mechanics. The stress intensity factor is calculated at elements along fracture tips according to the method described by Olson (2007). If the stress intensity factor exceeds the fracture toughness, K_{Ic}, the fracture is extended by creating a new adjacent element. A limitation of CFRAC is that the location and propagation path of potentially forming fractures must be speciﬁed in advance. For simulations with a single hydraulic fracture, this is not a major limitation because the fracture can reasonably assumed to be planar in most cases. But for simulations in a network, such as shown in Section 3.5, this assumption may result in some loss of realism, because in reality, hydraulic fractures may curve in response to stress heterogeneity, and new fractures may form off natural fractures.

### 2.6. Solving the coupled equations

The system of equations is coupled either with sequential coupling or explicit coupling with adaptive timestepping used to control error (Kim et al., 2011). In sequential coupling, the code sequentially solves (1) the shear traction equations, (2) the ﬂuid ﬂow and normal traction equations, and (3) the proppant transports equations (Fig. 3). In the scheme, each system of equations is solved while holding the primary variables from the other systems of equations constant.

The process is repeated until all equations are simultaneously satisﬁed within a certain tolerance (Settari and Mourits, 1998; Dean and Schmidt, 2009; Kim et al., 2011; Mikelić and Wheeler, 2013). Each of the systems of nonlinear equations is solved with a modiﬁed version of Newton–Raphson iteration. This approach is an extension of the method described by McClure and Horne (2013). After each system of equations is solved, the residuals of each element are calculated again and made dimensionless.

where R_{A}, R_{B1}, R_{B2}, and R_{C }are the (dimensional) residuals of the discretized forms of Eqs. (21) (or 22), (20), (2) and (5), respectively, in the ﬁnite volume scheme, and so their units are MPa, MPa, kg/m^{3}, and kg/m^{3}, respectively. R_{d,A}, R_{d,B1}, R_{d,B2}, and R_{d,C} are the dimensionless forms of those residuals. ς_{A} and ς_{B} are user-deﬁned reference stresses, and A_{s} is the area of fracture surface of an element. The value of E used for this nondimensionalization is the value from the previous timestep. In calculations A–D in Fig. 3, the convergence criteria are: A:max (R_{d,A,i)}<ε_{A}, B:max (R_{d,B1,i})<ε_{B1} and max (R_{d,B2,i}) < ε_{B2}, C: max (R_{d,C ,i})< ε_{C}, and D: max (R_{d,A,i})<ε_{D1}, max (R_{d,B1,i}) < ε_{D2}, and max (R_{d,B2,i}) < ε_{D3} where i is element number and ε_{A}, ε_{B1}, ε_{B2}, ε_{C,} ε_{D1}, ε_{D2}, and ε_{D3} are tolerances of each calculation. The tolerances used in the ﬁnal residual check for iterative coupling, ε_{D1}, ε_{D2}, and ε_{D3} are set to be looser than those used in the original system of equations.

When proppant concentration becomes very high in an element, the coupling between proppant concentration and the ﬂuid ﬂow equations becomes very strong, and convergence of the iterative coupling scheme becomes poor. To avoid this problem, when dimensionless proppant concentration goes above 0.8 anywhere in the simulation, the code switches to “explicit coupling”, in which the cycle shown in Fig. 3 is terminated after a single iteration (Kim et al., 2011). An alternative would be to solve all equations in a monolithic scheme with one large system of equations. However, this would be complex to implement and has not been attempted.

The danger to explicit coupling is that it could lead to signiﬁcant coupling error. For example, this could occur if solving the proppant system of equations introduces large error into the residuals of Eqs. (2) and (20)–(22). To minimize this problem, adaptive timestepping is used in which the residuals to Eqs. (2) and (20)–(22) are enforced to be below a certain tolerance. If they are too high, the timestep is aborted and repeated with a smaller dt. If the timestep is accepted, adaptive timestepping is performed to keep the coupling error near a certain target value. The new timestep size is selected using the method suggested by Grabowski et al. (1979):

where δni is the normalized residual of Eq. (2) or (20) (R_{d,B1} or R_{d,B2}) at element i, η_{targ} is a user speciﬁed target for the maximum residual, which is one fourth of a user speciﬁed tolerance, and ω is a factor that can ranges from zero to one (ω is set to one in this study). If δ_{i} > 4η_{targ} for any element, the timestep is discarded and repeated with a smaller value of dt. When the iterative coupling scheme is converging efﬁciently, coupling error can be driven down to very low levels within a few iterations. But with the explicit scheme, coupling error is difﬁcult to drive to zero without taking extremely small timesteps. Therefore, the convergence tolerances ε_{D2} and ε_{D3} are loosened to 10 times larger than their original values in Table 1.

**Table 1** Simulation settings for the base case.

We test the effect of this loosening by performing global mass balance calculations on the ﬂuid in the problem domain. This is performed by evaluating:

## 3. Results

First, we describe a base case simulation using the settings shown in Table 1 and the injection schedule shown in Table 2. Then, sensitivity analysis is performed by changing variables such as proppant density, proppant diameter, and viscosity of ﬂuid injected. One simulation is described with severe tip screen-out (TSO). Finally, a simulation is demonstrated with hydraulic stimulation of a complex fracture network. To validate the accuracy of the code, Shiozawa (2015) performed a simulation of hydraulic fracture propagation and proppant transport designed to closely imitate a published result from Dontsov and Peirce (2015). The simulation results showed good agreement.

**Table 2** Injection schedule.

### 3.1. Base case

The base case simulation is performed with the settings and pumping schedule shown in Tables 1 and 2, respectively. pad of clean ﬂuid is injected for 1,000s. After the pad injection,slurry (ﬂuid mixed with proppant) is injected for 2,000 s. Then clean ﬂuid is injected again for another 1,000 s to sweep the proppant out of the wellbore. Finally, the well is shut-in. For the entire injection period, the volumetric injection rate is constant, 0.04 m^{3}/s. The amount of clean ﬂuid injected in Stages 1 and 3 is same as the volume of the wellbore (40 m^{3}). The normalized proppant concentration, φ, is 0.2 during the slurry injection period.

After the well is shut-in, the ﬂuid pressure in the fracture decreases due to ﬂuid leakoff, and eventually the fracture closes. The fracture is meshed into elements of size one square meter.

#### 3.1.1. Simulation result

Figs. 4 and 5 show the results from the base case simulation at different points in time (1,000, 2,000 3,000, 4,000, 6,000, 10,000,15,000, and 1,000,000s), showing proppant distribution (with two different color scales), aperture, and ﬂuid pressure. The fracture is vertical and planar, and its height is limited to 100 m. The wellbore connection to the fracture is located at (x, z)^{¼ }(0, 0). A bi-wing fracture geometry forms in the fracture, but only one wing is shown in the ﬁgures.

**Fig. 4.** Base case simulation during pumping with the element size 1m 0.99m.

Fig. 4 a shows the fracture after 1,000 s of injection of clean ﬂuid. At this point, proppant slurry begins to be injected at the surface. Proppant does not enter the fracture immediately, because ﬁrst the clean ﬂuid in the wellbore must be displaced into the formation. Ideally, proppant should begin to enter the formation at 2,000 s, after 40 m^{3} of slurry have been injected at the surface.

However, as seen in Fig. 4b, a small amount of proppant has already entered the formation at 2,000s because of numerical dispersion in the calculation of proppant advection through the wellbore (Section 2.4). This could be reduced further by reﬁning the discretization of the wellbore. As can be seen in Fig. 4d, proppant accumulates near the fracture tip. The proppant cannot ﬂow all the way to the fracture tip because of proppant bridging. The aperture must be at least three times greater than the proppant diameter for proppant to ﬂow (Eq. (12)). Because the proppant tends to ﬂow in the center of the aperture, where the ﬂowing velocity is greatest, the ﬂowing fraction of proppant is typically about double the volumetric concentration of proppant unless the concentration is very high (Fig. 2). Comparing Fig. 4a–d, we can see proppant ﬂowing rapidly through the fracture, faster than the rate of fracture propagation, and accumulating near the tip. Fluid leakoff into the formation is another process that causes proppant concentration to increase.

At shut-in, the effect of gravity settling is subtle due to the high viscosity (100 cp) of the ﬂuid (Fig. 4d). The concentration of the accumulated proppant in Fig. 4d reaches its maximum (dimensionless proppant concentration equal to 1.0) in some elements near the fracture tip, an example of tip screen-out (TSO). In these elements, proppant is packed and immobile. The obstruction to ﬂow created by the TSO causes a discontinuity to develop in the distribution of ﬂuid pressure (Fig. 4d). The fracture continues propagating for a signiﬁcant time after shut-in. The screen-out bank is immobile, and so even though it is located near the tip at shut-in, it is located a signiﬁcant distance from the ﬁnal location of the crack tip. The obstruction to ﬂow created by the TSO bank causes the signiﬁcant time delay in crack tip extension, because ﬂow through the bank is relatively slow. This is apparent in the discontinuity in pressure distribution that occurs at the screen-out bank in Figs. 4d and 5a–c.

**Fig. 5.** Base case simulation after the well is shut-in with the element size 1m 0.99m.

Because the matrix permeability is so low (1000 nd), full mechanical closure does not occur for a substantial period of time after shut-in. The delay in mechanical closure gives the proppant time to settle to the bottom of the fracture. Settling is slowed by the high ﬂuid viscosity (100 cp), but is not slowed enough to prevent eventual settling. Some fracturing ﬂuids exhibit a yield strength, in which proppant ﬂow is zero unless shear stress exceeds a certain value. This could theoretically prevent proppant settling in the period before mechanical closure, but ﬂuid yield stress is not included in the model implemented for this study.

The ﬁnal distribution of proppant is seen in Fig. 5d. Proppant has only been emplaced into a fraction of the total fracture surface area. The total mass of ﬂuid and proppant injected is 150,640 and 23,400 kg, respectively. The total mass of ﬂuid leakoff is just 3955 kg at the end of injection (Fig. 4d) and 100,848 kg at the end of the simulation (Fig. 5d).

**Fig. 6.** Base case simulation with the element size 5m 4.76m.

#### 3.1.2. Mesh dependency

The base case simulation is repeated with a coarser simulation mesh in order to investigate mesh dependency and convergence.

Fig. 6a and b shows the results of the base case with the element size 5m 4.76m at the end of pumping and after closure, respectively. It is seen that proppant concentration at the proppant bank in Fig. 6a is lower than that in Fig. 4d. This may be due to greater numerical dispersion in the simulation with a coarser grid. This is also due to the difﬁculty of capturing highly localized accumulations of proppant with the coarse mesh. As described in Section 3.1.1, once proppant concentration reaches its maximum, it creates a ﬂow obstruction.

Thus, with a coarser mesh, the TSO bed takes longer to form, which has some effect on the results. In the comparison of Figs. 6b and 5d, even though they have the same trend in proppant distribution, there are subtle differences: (1) the former has one or two proppant banks while the latter one has three, (2) the fracture is shorter in the ﬁne-mesh simulation probably due to the earlier formation of the proppant bank, and (3) the propped aperture is smaller in coarse mesh simulation.

**Fig. 7. **The transition of proppant concentration at the speciﬁc points. Δx represents element length in the x-direction.

Fig. 7 shows the normalized proppant concentration versus time elapsed for simulations with different mesh sizes at four speciﬁc points: (x, z)=(50, 0), (100, 0), (150, 0), and (200, 0). The chosen mesh sizes are 10 m 9.09 m, 5 m 4.76 m, 2.5 m 2.44 m, and 1m 0.99 m. Generally, the simulation result should converge to the true solution as the mesh size decreases (as the number of elements increases). We can see convergence in Fig. 7 at x=50, 100, and 150 m. However, at x=200 m, the area where proppant accumulates and the proppant bank forms, the transition of the concentration is highly mesh dependent. These results suggest that coarser scale models are adequate for modeling the overall behavior of the system, but will vary signiﬁcantly in predicting the exact geometry of a tip screen-out proppant bank.

**Fig. 8.** Simulation results when ﬂuid viscosity is 10 cp.

#### 3.2. Sensitivity analysis

Four sensitivity analysis simulations are performed by changing ﬂuid viscosity, proppant density, proppant size, and formation permeability. Other settings are the same as in Table 1. Fig. 8a and b shows the results at the end of pumping (4,000 s) and after closure (1,000,000 s) for the simulation with ﬂuid viscosity of 10 cp.

Note that in all simulations, the leakoff viscosity is 1 cp (Section 2.1), and so leakoff occurs at roughly the same rate. The lower ﬂuid viscosity affects only ﬂuid and proppant ﬂow in the fracture itself. Because viscosity is lower than in the base case, the effect of gravity settling is more significant. Proppant does not accumulate at the tip because the proppant settles before it can reach the tip. The fracture is longer, and aperture is smaller than in the base case. The total mass of ﬂuid and proppant injected is 150,640 and 23,400kg, respectively. The total mass of ﬂuid leakoff is 5433kg at the end of injection (Fig. 8a) and 108,275kg at the end of the simulation (Fig. 8b).

*Fig. 9. Simulation results when proppant density is 1054 kg/m ^{3}.*

Fig. 9a and b shows simulation results with a proppant density of 1054 kg/m^{3}. Proppant with such low density is not common but is sometimes used (The Hole Solution Company, 2015). Because the proppant is only slightly denser than the ﬂuid, gravitational settling is negligible at the end of injection (Fig. 9a). Proppant settling is apparent after closure (Fig. 9b) but is still subtle. The total mass of ﬂuid and proppant injected is 150,640 and 9865 kg, respectively. The total mass of ﬂuid leakoff is 3956 kg at the end of injection (Fig. 9a) and 102,201 kg after closure (Fig. 9b).

**Fig.10.** Simulation results when proppant diameter is 200 microns.

Fig. 10a and b shows simulation results with a proppant diameter of 200 microns. This simulation is especially CPU intensive because the small proppant diameter forces the simulator to use very small timesteps when it switches to explicit coupling. To perform the simulation in a reasonable period of time, slightly larger elements are used than in the other simulations, 2.5m 2.44 m.

In this case, proppant bridging is less severe, and the proppant reaches the fracture tip (Fig. 10a). The fracture does not propagate much after the well is shut-in (Fig.10b) because the effective permeability of the proppant bank is lower with a lower proppant diameter (Equation (6)). The small size of the proppant allows it to settle all the way to the bottom of the fracture. The total mass of ﬂuid and proppant injected is 150,640 and 23,400 kg, respectively. The total mass of ﬂuid leakoff is 3,946 kg at the end of pumping (Fig. 10a) and 83,736 kg after closure (Fig. 10b).

**Fig. 11.** Simulation results when formation permeability is 1mD.

Fig. 11a and b shows simulation results with formation permeability of 1 mD. The fracture is much shorter because the ﬂuid leakoff rate is larger due to high permeability. The total mass of ﬂuid and proppant injected is 150,640 and 23,400 kg, respectively. The total mass of ﬂuid leakoff is 74,275 kg at the end of pumping (Fig. 11a) and 133,903 kg after closure (Fig. 11b), respectively. The high permeability simulation is the only simulation where nearly the entire hydraulic fracture remains propped after mechanical closure. Because closure occurs much earlier, the proppant has much less time to settle to the bottom of the fracture before it becomes immobile. These results underscore the importance and difﬁculty of proppant placement in very low permeability formations, where very low permeability can delay fracture closure and encourage settling.

#### 3.3. Factors affecting the timing of fracture closure and proppant settling

The simulation results demonstrate that in very low permeability formations, proppant settling before closure may be a major issue. Proppant is only placed across most of the created fracture surface area in the simulation with 1mD matrix permeability. In the simulation with ultralight proppant, proppant is successfully placed across most of the fracture height, but only to a point where screenout occurred about half-way along the eventual fracture length where a screenout occurs. In the other simulations, proppant settles to the bottom of the fracture before closure.

Proppant bridges and becomes immobile when the fracture aperture isless than approximately three times greater than the proppant diameter. Whether proppant will settle to the bottom of the fracture before closure depends on whether the time to immobilization, t_{immob}, deﬁned as the time when nearly all the average fracture aperture has become less than 3d, is less than the timescale of settling t_{settling}, the time when nearly all the injected proppant will have settled into a bank at the bottom of the fracture. The time to immobilization depends on two processes: ﬂuid leakoff and post-injection fracture propagation, both described below.

The 1D leakoff model (Vinsome and Westerveld, 1980) takes into account the changing ﬂuid pressure in the fracture over time. However, prior to closure, ﬂuid pressure in the fracture is fairly constant, and so it is roughly acceptable to use the Carter leakoff model, which assumes constant ﬂuid pressure in the fracture (Howard and Fast, 1957). The leakoff coefﬁcient is

where C is the Carter leakoff coefﬁcient, ΔP is differential pressure (the difference between the ﬂuid pressure in the fracture and the initial formation ﬂuid pressure), ϕ is the porosity of the formation around the fracture, c_{t} is the total compressibility calculated as c_{t}=c_{ϕ}+c_{f}, and μleak is the ﬂuid viscosity used for leakoff. ϕ and c_{ϕ} are assumed to be 0.03 and 0.00145 MPa^{-1}, respectively, and c_{t} is equal to 0.00185 MPa. In the simulations, ΔP is approximately equal to 4 MPa. Using these values, we can calculate C=1.68x10^{-5}m/s^{1/2} for the simulation with higher permeability and 5.32x10^{-7}m/s^{1/2} for the other simulations. The total volume of ﬂuid and proppant injected is Q_{tot}=0.04x4000=160 m^{3}.

The Carter leakoff model states that cumulative leakoff volume equals 4A_{s}C(t)^{1/2}, Therefore, if the fracture surface area is constant, the time to proppant immobilization due to leakoff can be calculated as

where N equals 3 and A_{s} is the fracture area. The numerator represents the volume of ﬂuid that must leakoff for closure to occur. It is equal to the total volume injected minus the ﬂuid remaining in the fracture when proppant can no longer ﬂow. The equation uses the simplifying assumptions that proppant cannot ﬂow once the average aperture is less than 3d and that the entire fracture surface area forms instantaneously at the beginning of injection and is then held constant.

A problem with applying Eq. (30) is that fracture surface area may continue to grow after closure. For example, in the base case simulation, the fracture half-length at shut-in is about 200 m, As around 40,000 m^{2}, and the ﬁnal fracture half-length is about 325m, A_{s} around 65,000 m^{2} (Figs. 4d and 5d). The changing fracture surface area hugely impacts the timing of closure, not only because this affects leakoff rate (the denominator in Eq. (30)), but also because it changes the volume of ﬂuid that needs to leak off for closure to occur (the numerator in Eq. (30)). For a sufﬁciently large fracture surface area, proppant could become immobile even without any leakoff. This would occur when As is large enough that the numerator of Eq. (30) is zero. For Qtot equal to 160 m^{3}, the average aperture is less than 3d for A_{s}=Q_{tot}/(3d), or 66,667 m^{2} for the values of Q_{tot} and d used in the base case simulation, which is close to the actual ﬁnal fracture surface area in the base case simulation. This suggests that in the base case simulation, the timing of proppant immobilization is controlled more by the rate of post-injection propagation than by ﬂuid leakoff. Eq. (30) yields values that differ by orders of magnitude, depending on whether As is assumed to be 40,000 m^{2} (the surface area at shut-in) or 65,000 m^{2} (the ﬁnal surface area). The base simulation shows a strong degree of post-injection propagation because the proppant bank that forms from TSO causes a screen that elevates net pressure during pumping and delays fracture propagation. Fluid can ﬂow quickly through an open fracture due to the cubic law, but ﬂow ﬂuids relatively much more slowly through a packed bed of proppant that behaves like a layer of porous media. This discussion suggests that effective vertical proppant placement in very low permeability formations may depend on forming a high surface area, low aperture fracture where proppant bridges soon after closure. This could be accomplished by avoiding TSO by using a sufﬁciently large pad or gradually increasing proppant concentration during a stage. TSO increases the net pressure during pumping and delays fracture growth, which delays the immobilization of proppant in the fracture and increases the probability of excessive settling. Larger diameter proppant will settle more rapidly but counterintuitively may be better for vertical proppant placement because it will bridge more readily. Higher viscosity will lead to better proppant carrying capacity but may reduce leakoff rate. For settling, using the asymptotic behavior of the function G_{p}, Dontsov and Peirce (2014) show that the settling velocity can be estimated as

Note that the unit of μ is MPa s. The time for settling can be calculated by dividing the fracture height h by the velocity plus the time when proppant starts to enter the fracture after the beginning of injection t_{p}, and so:

For all the simulations, t_{p} is approximately equal to 2,000. The ratio between the time for closure and settling time is

S is a dimensionless number that describes whether proppant will settle to the bottom of the fracture before mechanical closure. If it is more than approximately 1.0, the proppant will settle to the bottom of the fracture prior to closure. If less than approximately 1.0, proppant will be placed across the fracture height because closure will occur before the proppant has time to settle to the bottom.

Eq. (33) cannot be easily used if the fracture surface area is changing strongly after shut-in, as in most of the simulations in this paper, because timmob becomes too difﬁcult to estimate. The value of timmob depends on the rate of post-injection fracture propagation. However, it is a useful heuristic for understanding vertical process placement. In the higher permeability simulation, limited fracture propagation occurred after shut-in. In this case, the value of S can be estimated as 39, 281 s/97, 663 s, around 0.40. As suggested by the value of S less than one, there is good vertical proppant placement in this simulation (Fig. 11b).

**Fig. 12. **Comparison of simulation results. TSO signiﬁcantly limit fracture propagation. (a) With proppant. (b) Without proppant.

#### 3.4. Tip screen-out (TSO)

In order to create a scenario in which severe TSO occurs, a schedule is designed with injection of slurry for 4,000 s without a pad of clean ﬂuid injection. The normalized proppant concentration in the ﬂuid injected is 0.4, which is higher than the base case simulation, and the volumetric injection rate of slurry is constant, 0.04 m^{3}/s. Proppant density and diameter are 1054 kg/m^{3} and 200 microns, respectively. Other settings are the same as in Table 1. Fig. 12a shows the simulation result at the end of injection. Fig. 12b shows the result of a simulation with identical settings, except that clean water is injected with no proppant. The simulation shows that TSO signiﬁcantly affects fracture propagation, aperture, and ﬂuid pressure. The proppant accumulates at the fracture tip and forms a low permeability bank that prevents further fracture propagation and causes a sharp buildup of ﬂuid pressure (Nolte and Smith, 1981).

**Fig. 13. **Time elapsed and pressure above closure in log–log scale.

Fig. 13 shows time elapsed and bottomhole pressure minus minimum horizontal stress ( BHP − σ_{yy}) in log–log scale. The TSO causes increasing net pressure during pumping. The total mass of ﬂuid and proppant injected is 122,560 and 39,462 kg, respectively, for the TSO simulation (Fig. 12a), and 160,000 kg of ﬂuid for the simulation without proppant injection (Fig. 12b). The total amounts of ﬂuid leakoff for the simulations shown in Fig. 12a and b are 3235 and 4045 kg respectively.

**Fig. 14.** Geometry of natural fractures (gray) and path taken by the hydraulic fracture (red). For visibility, fractures are partially transparent. (For interpretation of the references to color in this ﬁgure caption, the reader is referred to the web version of this paper.)

#### 3.5. Complex fracture network

A simulation is performed in the complex natural fracture geometry shown in Figs. 14 and 15. It is assumed that only a single hydraulic fracture will form, and it will be planar and height conﬁned. The natural fractures are assumed to be vertical.The natural fractures are meshed into rectangular elements and the governing equations (described in Section 2) are solved in the same way as in the hydraulic fracture elements. Fluid and proppant ﬂow between hydraulic fracture elements and natural fracture elements is calculated in the same way as any other pair of elements. Thus, the model does not account for the possibility that proppant has difﬁculty turning the corner into natural fractures due to inertial effects. For natural fracture elements, E_{0} is constant and set equal to E_{hf max,resid }.

**Fig. 15. **Natural fracture distribution (viewed from the top). The black line at the center represents the wellbore.

The simulation settings and pumping schedule are otherwise identical to the base case simulations, as shown in Tables 1 and 2. Figs. 16 and 17 show the proppant distribution at the end of pumping (4,000s) and after closure (1,000,000s) viewed from the side of the hydraulic fracture. In Fig. 16, it can be seen that proppant accumulates in the hydraulic fracture near intersections with natural fractures. The natural fractures do not open sufﬁciently to allow proppant to ﬂow into them. However, ﬂuid ﬂow slows down at these natural fracture intersections because the ﬂow is diverted into the natural fractures. This reduces the proppant-carrying ability of the ﬂuid and removes ﬂuid, causing an accumulation of proppant. The asymmetric distribution of proppant and fracture propagation after closure in Fig. 17 is caused by the interactions with the natural fractures.

**Fig. 16. **Proppant distribution at the end of pumping (4,000 s), viewing the hydraulic fracture from the side. Proppant accumulates at intersections with natural fracture elements and near the fracture tip.

**Fig. 17. **Proppant distribution after fracture closure (1,000,000 s), viewing the hydraulic fracture from the side. The natural fractures causes asymmetric fracture propagation and proppant distribution.

## 4. Conclusions

Simulations of proppant transport are performed in an Eulerian–Eulerian framework in a three-dimensional hydraulic fracturing simulator, CFRAC. The ﬂow equations (ﬂuid and proppant) and the mechanical equations are coupled with either iterative coupling or explicit coupling with adaptive timestepping. The approach introduced by Dontsov and Peirce (2014) is used to capture the transition from Poiseuille to Darcy ﬂow, proppant settling due to gravity, and the processes that cause the ﬂowing proppant fraction to be different from the volumetric fraction of proppant. Except at high proppant concentration, the ﬂowing proppant fraction is greater than the volumetric fraction of proppant. This can cause tip screen-out (TSO), even if the leakoff rate is very low. We introduce a framework that allows the simulator to seamlessly describe the process of fracture closure, including a nonlinear joint closure law for describing fracture compliance and roughness and including the effect of proppant accumulation into a packed layer between the fracture walls.

Iterative coupling convergence is poor when dimensionless proppant concentration is greater than approximately 0.8. In these cases, the simulator uses explicit coupling, in which iteration is terminated after a single iteration, and timestep duration is adaptively varied to minimize coupling error. When using explicit coupling, it is necessary to use a relatively loose error tolerance to achieve reasonable performance. In our simulations, this leads to a cumulative global mass balance error of around 1–3%, which is adequate for most practical applications.

Our sensitivity analysis simulations show the effects of ﬂuid viscosity, proppant density, proppant size and formation permeability. The simulations show that proppant settling and tip screen-out can cause serious problems for proppant placement. In low permeability formations, mechanical closure may take hours or days, giving proppant time to settle to the bottom of the fracture even if a high viscosity fracturing ﬂuid is used. Proppant immobilization can occur due to bridging long before the fracture fully mechanically closes. Bridging is especially critical for effective vertical proppant placement in very low permeability formations because the aperture change due to leakoff is very slow. This suggests, counterintuitively, that smaller diameter proppant may be worse for vertical proppant placement. Smaller proppant will settle more slowly but will require a much smaller aperture to become immobilized through bridging. This also suggests that it is critical to optimize injection schedule to avoid tip screen-out. TSO increases net pressure and delays fracture propagation, resulting in a shorter, higher aperture fracture that is not conducive to proppant bridging.

We perform one simulation with leakoff into discrete natural fractures. In this simulation, proppant does not ﬂow into the natural fractures because their aperture is too low. However, due to ﬂuid leakoff into the natural fractures, proppant accumulates at the intersections.

##### Acknowledgment

The authors would like to thank Jian Huang and Uno Mutlu for helpful advice, Weatherford International Inc. and Cockrell School of Engineering at The University of Texas at Austin for ﬁnancially supporting this research, and Egor Dontsov at The University of Houston for providing tables with the dimensionless functions controlling slurry ﬂow.

##### Corresponding author.

E-mail addresses: sogo@utexas.edu (S. Shiozawa), mark@mccluregeomechanics.com (M. McClure).

http://dx.doi.org/10.1016/j.petrol.2016.01.002

0920-4105/& 2016 The Authors. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

## References

- Adachi, J., Siebrits, E., Peirce, A., Desroches, J., 2007. Computer simulation of hydraulic fractures. Int. J. Rock Mech. Min. Sci. 44 (5), 739–757.
- Barton, N., Bandis, S., Bakhtar, K., 1985. Strength, deformation and conductivity coupling of rock joints. Int. J. Rock Mech. Min. Sci. Geomech. Abst. 22 (3),121–140.
- Boronin, S., Osiptsov, A., 2014. Effects of particle migration on suspension ﬂow in a hydraulic fracture. Fluid Dyn. 49 (2), 208–221.
- Boyer, F., Guazzelli, É., Pouliquen, O., 2011. Unifying suspension and granular rheology. Phys. Rev. Lett. 107 (18), 188301.
- Clifton, R., Wang, J., 1988. Multiple ﬂuids, proppant transport, and thermal effects in three-dimensional simulation of hydraulic fracturing. In: SPE Annual Technical Conference and Exhibition, Houston, Texas, 2-5 October.
- Constien, V.G., Hawkins, G.W., Prud'homme, R.K., Navarrete, R., 2000. Performance of Fracturing Materials. In: Economides, M.J., Nolte, K.G. (Eds.), Reservoir Stimulation, third ed. Wiley & Sons, Chichester UK, Chapter 8.
- Crouch, S.L., Starﬁeld, A.M., Rizzo, F., 1983. Boundary element methods in solid mechanics. J. Appl. Mech. 50, 704.
- Dean, R.H., Schmidt, J.H., 2009. Hydraulic-fracture predictions with a fully coupled geomechanical reservoir simulator. SPE J. 14 (4), 707–714.
- Dontsov, E., Peirce, A., 2014. Slurry ﬂow, gravitational settling and a proppant transport model for hydraulic fractures. J. Fluid Mech. 760, 567–590.
- Dontsov, E., Peirce, A., 2015. Proppant transport in hydraulic fracturing: crack tip screen-out in KGD and P3D models. Int. J. Solids Struct. 63, 206–218.
- Einstein, A., 1905. Eine neue bestimmung der moleküldimensionen (Ph.D. thesis). Buchdruckerei KJ Wyss.
- Eskin, D., Miller, M., 2008. A model of slurry ﬂow in a fracture. J. Can. Pet. Technol. 47 (3).
- Grabowski, J.W., Vinsome, P.K., Lin, R.C., Behie, G., Rubin, B., 1979. A fully implicit general purpose ﬁnite-difference thermal model for in situ combustion and steam. In: SPE Annual Technical Conference and Exhibition, Las Vegas, Nevada, 23-26 September.
- Howard, G.C., Fast, C., 1957. Optimum ﬂuid characteristics for fracture extension. Drill. Prod. Prac. 24, 261–270.
- Hu, H.H., Patankar, N.A., Zhu, M., 2001. Direct numerical simulations of ﬂuid–solid systems using the arbitrary Lagrangian–Eulerian technique. J. Comput. Phys. 169 (2), 427–462.
- Jaeger, J.C., Cook, N.G., Zimmerman, R., 2007. Fundamentals of Rock Mechanics, fourth ed.. Blackwell Publishing, Malden, MA.
- Karimi-Fard, M., Durlofsky, L., Aziz, K., 2004. An efﬁcient discrete fracture model applicable for general purpose reservoir simulators. SPE J. 9 (2), 236–277.
- Kim, J., Tchelepi, H., Juanes, R., 2011. Stability and convergence of sequential methods for coupled ﬂow and geomechanics: drained and undrained splits. Comput. Methods Appl. Mech. Eng. 200 (23), 2094–2116.
- Krieger, I.M., Dougherty, T.J., 1959. A mechanism for non-newtonian ﬂow in suspensions of rigid spheres. Trans. Soc. Rheol. 3 (1), 137–152.
- Mack, M.G., Warpinski, N.R., 2000. Mechanics of hydraulic fracturing. In: Economides, M.J., Nolte, K.G. (Eds.), Reservoir Stimulation, third ed. Wiley & Sons, Chichester UK, Chapter 6.
- Mayerhofer, M., Economides, M., 1993. Permeability estimation from fracture calibration treatments. In: SPE Western Regional Meeting, Anchorage, Alaska, 26-28 May.
- McClure, M., 2014. The potential effect of network complexity on recovery of injected ﬂuid following hydraulic fracturing. In: SPE Unconventional Resources Conference, The Woodlands, Texas, 1-3 April.
- McClure, M., Horne, R.N., 2013. Discrete Fracture Network Modeling of Hydraulic Stimulation: Coupling Flow and Geomechanics. Springer International Pub- lishing, Cham, Switzerland.
- McClure, W.M., Babazadeh, M., Shiozawa, S., Huang, J., 2016. Fully coupled hydro- mechanical simulation of hydraulic fracturing in 3D discrete fracture networks. SPE J. http://dx.doi.org/10.2118/173354-PA.
- Mikelić, A., Wheeler, M.F., 2013. Convergence of iterative coupling for coupled ﬂow and geomechanics. Comput. Geosci. 17 (3), 455–461.
- Mobbs, A., Hammond, P., 2001. Computer simulations of proppant transport in a hydraulic fracture. SPE Prod. Facil. 16 (2), 112–121.
- Mooney, M., 1951. The viscosity of a concentrated suspension of spherical particles. J. Colloid Sci. 6 (2), 162–170.
- Morris, J.F., Boulay, F., 1999. Curvilinear ﬂows of noncolloidal suspensions: the role of normal stresses. J. Rheol. 43 (5), 1213–1237.
- Nolte, K.G., 1979. Determination of fracture parameters from fracturing pressure decline. In: SPE Annual Technical Conference and Exhibition, Las Vegas, Nevada,23-26 September.
- Nolte, K.G., Smith, M.B., 1981. Interpretation of fracturing pressures. J. Pet. Technol. 33 (9), 1767–1775.
- Okada, Y., 1992. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 82 (2), 1018–1040.
- Olson, J., 2007. Fracture aperture, length and pattern geometry development under biaxial loading: a numerical study with applications to natural, cross-jointed systems. Geol. Soc. Lond. Spec. Publ. 289 (1), 123–142.
- Ouyang, S., Carey, G., Yew, C., 1997. An adaptive ﬁnite element scheme for hydraulic fracturing with proppant transport. Int. J. Numer. Methods Fluids 24 (7), 645–670.
- Rice, J.R., 1993. Spatio-temporal complexity of slip on a fault. J. Geophys. Res. 98 (B6), 9885–9907.
- Settari, A., Mourits, F., 1998. A coupled reservoir and geomechanical simulation system. Spe J. 3 (03), 219–226.
- Shiozawa, S., 2015. Designing enhanced geothermal and hydraulic fracturing systems based on multiple stages and proppant (Master's thesis). The University of Texas at Austin.
- The Hole Solution Company, 2015. FracBlack HT, URL 〈http://www.sundrilling.com/documents/FracBlackHTProductBulletin-March2015PDF.pdf〉.
- Tomac, I., Gutierrez, M., 2015. Micromechanics of proppant agglomeration during settling in hydraulic fractures. J. Pet. Explor. Prod. Technol. 5 (4), 417–434.
- Tsai, K., Fonseca, E., Lake, E., Degaleesan, S., 2012. Advanced computational modeling of proppant settling in water fractures for shale gas production. Spe J. 18 (1), 50–56.
- Vinsome, P., Westerveld, J., 1980. A simple method for predicting cap and base rock heat losses in' thermal reservoir simulators. J. Can. Pet. Technol. 19 (3), 87–90.
- Wallace, J., Kabir, C., Cipolla, C., 2014. Multiphysics investigation of diagnostic fracture injection tests in unconventional reservoirs. In: SPE Hydraulic Fracturing Technology Conference, The Woodlands, Texas, 4-6 February.
- Weng, X., Kresse, O., Cohen, C.-E., Wu, R., Gu, H., 2011. Modeling of hydraulic- fracture-network propagation in a naturally fractured formation. SPE Prod. Oper. 26 (4), 368–380.
- Willis-Richards, J., Watanabe, K., Takahashi, H., 1996. Progress toward a stochastic rock mechanics model of engineered geothermal systems. J. Geophys. Res. 101 (B8), 17481–17496.
- Witherspoon, P., Wang, J., Iwai, K., Gale, J., 1980. Validity of cubic law for potential size effects in experimental determination of the hydraulic properties of a fracture. Water Resour. Res. 16 (6), 1–16.
- Zhang, Z., Chen, Q., 2007. Comparison of the Eulerian and Lagrangian methods for predicting particle transport in enclosed spaces. Atmos. Environ. 41 (25), 5236–5248.