**3. Non-planar 3D model of hydraulic fractures**

The numerical simulator essentially follows a well accepted numerical methodology [24, 25]. In this numerical framework, the rock mechanics of fracture growth takes place in a linear elastic medium and therefore can be modelled by boundary integral equations based on fundamental solutions of dislocation segments. The fracture growth is controlled by both fracture mechanics criteria and fluid volume balance at the fracture tip. The formulation adequately portrays both viscosity and toughness controlled fracture growth as well as the 'fluid lag' behaviour [25]. The proppant laden frac fluid is idealized as an incompressible power-law fluid commonly characterized by parameters *K* and *n*:

$$
\pi = K \dot{\gamma}^n \tag{7}
$$

Where *τ* is the fluid shear stress and *γ*˙ the fluid shearing rate. *K* and *n* are termed the consistency index and power law index, respectively. They are naturally dependent on the proppant concentration, and this effect of proppant concentration on fluid rheology is accounted for in the simulator using Shah's model [27]. The proppant transport and settlement are computed via the transport equation and empirical equations of particles settling between parallel plates [28]. The fluid loss to the formation is described by a Carter type leak-off model [25].

The ensuing flow of proppant laden frac fluid within a fracture is captured according to the Reynolds lubrication theory, and the derived non-linear fracture growth and fluid flow equations are now solved in a coupled manner in which mass conservation (i.e., frac fluid and proppant) are strictly enforced. The equations are discretized on the same structured grid, and they form a system of moving boundary, transient coupled equations of fracture width and fluid pressure. A new adaptive time integration algorithm has been developed to provide numerically stable solutions for a wide range of power-law fluid properties (especially when the viscosity is low).

Although the proppant transport is calculated at each time step based on a volume conserva‐ tion equation, this computation is decoupled from the process of solving the fluid flow equation for the sake of computational efficiency. Specifically, the fluid pressure and fracture width are first obtained by 'freezing' the proppant concentration associated with each element. The contribution of proppant transport to the volume concentration is then updated based on the calculated fluid velocity field. This proves to be effective as well as efficient since the employed time step is generally small.

#### **3.1. Non-planar fracture propagation**

**Figure 2.** Interaction between two fractures

666 Effective and Sustainable Hydraulic Fracturing

**Figure 3.** Fracture spacing versus stress contrast

In dealing with non-planar fracture growth, the heterogeneous formation is modelled as multiple isotropic parallel layers with heterogeneities characterized by variations in in-situ stress, elastic stiffness modulus, Poisson's ratio, fracture toughness, pore pressure, leak-off coefficients, and spurt-loss coefficients. All fractures are assumed to be vertical but they can turn in any horizontal direction. For each fracture, the growth direction is determined by the stress state of the leading element at the fracture tip. These assumptions exclude the application for certain geological formations. In particular, geological stress regimes as a result of reverse faulting may present difficulties in modelling as the hydraulic fracture would most likely grow in a horizontal plane. However, the assumptions made are sufficient and adequate for typical deep unconventional resource formations where the vertical stress is generally the maximum principal stress, or at least greater than the minimum horizontal principal stress – resulting in vertical fracture growth.

**4. Comparing simulation results with analytical solution and experiment**

The Geomechanical Interaction of Multiple Hydraulic Fractures in Horizontal Wells

http://dx.doi.org/10.5772/56385

669

The complete validation of a hydraulic fracturing simulator is enormously challenging, because of the difficulty in obtaining the general analytical solutions or constructing realistic laboratory experiments. A carefully conducted laboratory hydraulic fracturing experiment has been published [32]. Figure 4 shows our simulation of the experiment. The comparison is

excellent not only in the created fracture geometry, but also the pumping parameters.

**Figure 4.** Comparing with experimental result of [32]

of 0.2 m3

set to 1 with no spurt loss.

**5. Frac fluid viscosity and fractures height growth**

A single stage of four hydraulic fractures is considered. Parameters chosen broadly resemble a typical shale gas development field case. The horizontal wellbore is placed at the relative depth of 0m as shown in Figure 5. Fracture injection points are spaced at 25m (82ft) apart and each fracture has 12 perforations. Fluid is injected into the wellbore at constant injection rate

In this instance, the minimum horizontal stress gradient is set to 14kPa/m (0.62psi/ft) with no upper stress barrier to restrict height growth. The difference between minimum and maximum horizontal stress is chosen to be more than 20%, so that we do not cause the fracture direction to re-orient drastically. The key parameters are shown in Figure 5. The fluid leak-off factor is

/s (~75 bbl/min). Formation heterogeneities are limited to only in-situ stress variation.

Fractures advance when the maximum tensile stress ahead of the crack tip exceeds the intrinsic tensile strength of the rock. In linear elastic fracture mechanics, the stress at the crack tip is singular, therefore the stress intensity factor, which correlates to the strength of such a singular stress field, is generally used to determine fracture tip propagation [29]. In practice, hydraulic fracturing processes take place under relatively large compressive in-situ earth stress that is normally several orders of magnitude higher than the rock strength. Therefore, the process is dominated by the balance of the compressive stress and fluid pressure at the fracture tip region. A group of 'virtual' elements is placed along the crack front and at each time step a check is made on the stress status of these elements. The virtual elements are allowed to be active as part of the new fracture surface if the potential fluid pressure can overcome the compressive stress plus the strength of the rock.

Importantly, the interaction of multiple fractures may result in significant shearing displace‐ ments along the fracture surface, causing the fracture growth to follow a curved path, and the computational technique leads these virtual elements to be oriented according to the stress field so that the local minimum principal stress is perpendicular to the new fracture surface as the fracture front advances.

An outline of the mathematical description of the fracture stiffness matrix, the displacement discontinuity method is given in [19], and will not be repeated here. The numerical approach is robust and efficient but less accurate than the variational boundary integral method which employs second order elements [30, 31]. For the current application, computational efficiency is deemed to outweigh the importance of a moderate improvement in accuracy. [19] also describes the coupling of derived non-linear fracture growth to flow equation and mass conservation.

#### **3.2. Wellbore hydraulic model**

The fractures are assumed to be initiated at the specified locations of a cluster of perforations. Near wellbore formation stress concentration is not considered. The downhole pressure is defined as the fluid pressure just upstream of all injection points. It is obvious that the distribution of the corresponding injection at each fracture follows wellbore fluid mechanics, depending on the injection area connected with the fracture and the fluid pressure within the fracture. An approximation is made by using Bernoulli equation to capture the fluid distribu‐ tion. This wellbore hydraulic model is able to account directly for the empirically derived or calibrated frictional pressure drop at the perforation and along the wellbore. Therefore, it is possible to account for limited entry perforation designs commonly employed in multiple hydraulic fracture design.
