Modeling of complex hydraulic fractures in naturally fractured formation

Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Journal of Unconventional Oil and Gas Resources journal homepage: www.elsevier.com/locate/juogr

Review

Modeling of complex hydraulic fractures in naturally fractured formation Xiaowei Weng Schlumberger, 110 Schlumberger Drive, MD-2, Sugar Land, TX 77478, USA

a r t i c l e

i n f o

Article history: Received 15 February 2014 Revised 14 July 2014 Accepted 21 July 2014 Available online xxxx Keywords: Hydraulic fracture model Complex fractures Fracture network Hydraulic–natural fracture interaction Stress shadow

a b s t r a c t This paper presents a general overview of hydraulic fracturing models developed and applied to simulation of complex fractures in naturally fractured shale reservoirs. It discusses the technical challenges involved in modeling complex hydraulic fracture networks, the interaction between a hydraulic fracture and a natural fracture, and various models and modeling approaches developed to simulate hydraulic fracture–natural fracture interaction, as well as the induced large scale complex fractures during fracturing treatments. Ó 2014 Elsevier Ltd. All rights reserved.

Contents Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Challenges in modeling of complex fractures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Planar hydraulic fracture models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fracture opening displacement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fluid flow in the fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fracture propagation criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . HF–NF interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Analytical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Modeling of large-scale complex fractures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2D Boundary element based models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3D boundary element based models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Coupled geomechanics and reservoir models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Discrete fracture network and distinct element method based models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Lumped elliptic P3D network models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Cell-based P3D complex fracture network model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00

E-mail address: [email protected] http://dx.doi.org/10.1016/j.juogr.2014.07.001 2213-3976/Ó 2014 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

2

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

Introduction It has long been observed in mine-back experiments and corethrough of fractured formations (Warpinski and Teufel, 1987; Jeffrey et al., 1994, 2009; Warpinski et al., 1993) that hydraulic fracture interaction with natural fractures can result in branching and offset at the natural fractures and consequently lead to complex fractures. Fig. 1 shows an example of complex parallel fractures and offsets created as a hydraulic fracture propagates through natural fractures and zone boundaries, observed by Warpinski and Teufel (1987) in a mine-back experiment. However, it was not always clear if these complexities are only small-scale features relative to an otherwise planar fracture at large scale and whether they also occur in formations at greater depth. Limited direct observations available at depth by coring through the hydraulically fractured interval at the GRI/DOEsponsored Multiwell Experiment Site (Warpinski et al., 1993) also revealed multiple closely spaced hydraulically induced fractures, filled with the residue of the fracturing fluid. In spite of the general awareness of potential complexity in the hydraulically induced fractures, based on the observations like these, over the years hydraulic fracturing treatment design continued to be simulated based on the models that assume a planar fracture. Some of the models introduced correction factors to account for increased resistance to fracture propagation (apparent fracture toughness), to fluid flow in the fracture (wall roughness effect), or the effect of multiple parallel fractures that increases apparent fracture stiffness and reduces width. These factors tend to be used as pressure-matching parameters and have limited predictive capability. In the last decade, following the success of horizontal drilling and multistage fracturing in the Barnett shale, exploration and drilling activities in shale gas and shale oil reservoirs have skyrocketed in the US and abroad. Economic production from these reservoirs depends greatly on the effectiveness of hydraulic fracturing stimulation treatment. Microseismic measurements and other evidence suggest that creation of complex fracture networks during fracturing treatments may be a common occurrence in many unconventional reservoirs (Maxwell et al., 2002; Fisher et al., 2002; Warpinski et al., 2005). The created complexity is strongly influenced by the pre-existing natural fractures and in-situ stresses in the formation. However, due to the lack of industry’s modeling

Fig. 1. Complex fractures observed in mine-back experiments, by Warpinski and Teufel (1987).

capability in simulating complex fractures and lack of proper characterization of key reservoir properties, drilling, completion, and stimulation designs for the unconventional reservoirs relied heavily on imprecise estimates of the stimulated reservoir volume (SRV) from microseismic observations and through a highly inefficient trial-and-error approach. Not being able to accurately simulate complex fractures generated during the fracture treatment presented a major limitation that promoted a cookie-cutter completion and fracture design approach rather than one that is based on engineering approach to optimize the fracture parameters and production according to formation properties. Fracture simulation can provide information such as induced overall fracture length and height, propped versus unpropped fracture surface areas, proppant distribution and its conductivity, etc., all of which influence the short- and long-term production from the unconventional reservoir, and cannot be obtained from microseismic measurement alone (Cipolla et al., 2011b). However, modeling the process of hydraulic fracture network creation and interaction between hydraulically induced fractures and natural fractures presents many technical challenges. Significant progress has been made in recent years in the development of complex fracture models to address the needs for more suitable design tools for the unconventional reservoirs than the conventional planar fracture models. However, some aspects of this complex fracturing process are still not fully understood in terms of their impact or importance to the overall fracture geometry creation, or the complexity of simulating them is still beyond the current modeling capabilities or requires computation power or time that is not practical for engineering use. Therefore, these models will continue to evolve in the coming years. One of the difficulties in fracturing design is the lack of clear understanding of the nature of fracture complexity created during fracture treatment. Microseismic monitoring does not provide sufficient resolution to delineate the exact hydraulic fracture planes. Microseismic events are mostly attributed to shear failures along natural fractures or faults surrounding a hydraulic fracture (Rutledge et al., 2004; Williams-Stroud et al., 2012). The events cloud forms a ‘‘halo’’ surrounding the hydraulic fracture. In conventional sandstone formations, the observed events cloud has a relatively narrow width, whereas in unconventional reservoirs, a much wider events cloud is often observed (Fisher et al., 2002). A wide microseismic cloud may possibly be explained by either deep fluid penetration into natural fractures in the shale while the induced hydraulic fracture remains planar or simple (Savitski et al., 2013), or by the creation of complex hydraulic (tensile) fracture network. Although deep fluid penetration into a highly permeable and initially well-connected natural fractures network is certainly possible (Zhang et al., 2013), many unconventional plays have very low effective permeability, and observation of cores shows that most natural fractures in these shales are mineralized (Gale et al., 2007; Gale and Holder, 2008; Han, 2011; WilliamsStroud et al., 2012). Therefore, in very low permeability shale, fluid penetration in the natural fractures network is limited. Fluid penetration into natural fractures can also occur due to dilation of natural fractures as a result of shear, but this typically occurs under the condition of large stress anisotropy and for natural fractures oriented 30 to 60° from the principal stress directions (Murphy and Fehler, 1986). For many shale reservoirs where the tectonic environment is relaxed and the difference between the horizontal stresses is low, a wide microseismic cloud is a strong indication that complex, tensile open hydraulic fracture networks are created, although the hydraulic fractures may follow the paths of the natural fractures. A field case in the Barnett shale presented by Fisher et al. (2002), in which fracturing fluid unexpectedly connected to and brought down the production of several adjacent wells not

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

on the expected fracture plane, provided the supporting evidence of complex hydraulic fracture networks. In the analysis of another Barnett case, Cipolla et al. (2010) showed that the predicted fracture length from a planar fracture model far exceeded the fracture length indicated by the microseismic data, unless a very low fluid efficiency (less than 10%) is assumed in the simulation in order for a planar fracture to accommodate the large volume of fluid injected. Such low efficiency is not consistent with the very slow pressure decline typically observed during the shut-in following the pumping period in most shale formations. In contrast, complex hydraulic fracture networks can explain the much larger fluid volume stored in fracture networks for the same overall fracture network length and yet still high fluid efficiency due to low leakoff. However, in a geological setting with high tectonic stress, shear fracturing may be an important mechanism accounting for the observed wide microseismic cloud and potential permeability enhancement induced by shear. Properly constructed complex fracture models and/or geomechanics models can help answer these questions and provide the tools for optimizing the fracture design and completion strategy. This paper attempts to give an overview of the technical challenges for modeling the complex fractures and various numerical modeling approaches. Some of the advantages and disadvantages of various approaches discussed in the paper are purely based on the author’s personal opinions. Challenges in modeling of complex fractures For a hydraulic fracture propagating in a formation that contains pre-existing natural fractures, or mechanically weak planes relative to the rock matrix, the interaction between the hydraulic fracture (HF) and natural fractures (NF) could cause fluid loss into the NF, dilation of the NF either due to shear or in tension, or even branching or alteration of the HF path, leading to complex fractures. Fig. 2 depicts a hydraulic fracture network created when pumping into a perforation cluster in a horizontal well.

3

Fig. 2 shows some of the scenarios of HF interaction with NF that can lead to fracture branching and complexity. They are discussed below: 1. Direct crossing. When a NF has strong mechanical bonding and/or is subjected to high normal stress, the tensile stress concentration at the tip of approaching HF is readily transmitted across the NF interface to the rock on the opposite side of the hydraulic fracture, causing the rock to fail in tension and allowing the hydraulic fracture to directly propagate through the NF without change of direction. Consequently the HF propagates through the formation as a planar fracture. However, if the fluid pressure can exceed the closure stress acting on the NF, they will open in tension and become a part of a now nonplanar fracture network. 2. HF arrested by NF. This scenario occurs when the NF interface is weaker than the rock matrix and the stress condition is such that the interface fails in shear and slips. Consequently, the tensile stress at the tip of the approaching HF is not sufficiently transmitted to the opposite side of the NF interface to cause the rock to fail in tension. And the HF growth is hence arrested by the NF. If the fluid pressure in the HF continues to increase, it can exceed the closure stress acting on the NF and cause the NF to be opened in tension and become a part of the hydraulic fracture network. 3. Crossing with an offset. It is often observed in laboratory and mineback experiments that when a HF crosses a NF, it can do so with a small offset at the interface as shown in Fig. 2. The offset is typically in the order of one to a few inches. The offset is created due to localized interface separation and shear slip at the point where the HF intersects the NF, which shifts the stress concentration away from the intersection point to the tip of opening/shear slip region (Thiercelin and Makkhyu, 2007). 4. Intersecting natural fractures. Once fluid pressure exceeds the closure stress on the NF, the NF opens up in tension and becomes a part of the HF network. If the NF intersects

Fig. 2. Map view of induced hydraulic fracture network depicting various scenarios of interaction with natural fractures that can lead to fracture branching and complexity.

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

4

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

another NF, when the fluid front reaches the intersection, the HF may branch again at the intersection as long as fluid pressure exceeds the closure stress on NF. 5. Branching or turning of fracture at end of NF. For a HF following the path of a NF to its end, there is no longer a weak plane for fluid to preferably open. Consequently, the fracture either turns itself to align with the preferred fracture direction or creates a T-shaped branch. 6. Shear slip along NF. If the fluid pressure in the NF stays below its closure stress, the fracture interface will not separate in tension. However, it can fail in shear. The shearinduced interfacial slip causes dilation and enhances the permeability of the NF, which can potentially enhance production. The occurrence of shear failure depends on the normal and shear stresses applied on the NF, which depends on the in-situ principal stresses, angle of the NF relative to the in-situ stresses, the fluid pressure (which depends on pressure diffusion in the NF), and interfacial frictional properties. A major challenge in numerical modeling of complex hydraulic fractures is due to the different scales the problem presents. The mechanics that controls the various behaviors at the fracture intersection point is a much localized phenomenon. To solve the problem numerically will require a very fine mesh to compute the local stresses and displacements accurately to capture the correct crossing behaviors. On the other hand, the overall hydraulic fracture network is of much greater scale, in the order of hundreds to a few thousands of feet. Many of the modeling efforts had focused on only one aspect, either looking at the fracture intersection problem to better understand the crossing phenomenon or looking at the fluid flow in the large-scale fracture network with a predefined fracture network pattern. For a comprehensive complex fracture simulator, both the large-scale network and localized crossing processes need to be properly modeled in the simulator. In the large-scale modeling of a fracture network, the mechanical interaction among propagating hydraulic fractures also needs to be taken into account. In the most common horizontal well completion in unconventional reservoirs, a completion stage consists of multiple perforation clusters from which multiple fractures are initiated. Even in the case of noncomplex, planar fractures, the interaction among the fractures can lead to uneven fracture growth. Most conventional planar fracture models do not account for the fracture interaction, or the so-called stress shadow effect. In the case of complex fracture network, the fracture interaction can be even stronger due to typically higher fracture density. Therefore, stress shadowing must be considered in the fracture simulator. In the following, we will start with an overview of various modeling approaches used for simulating the propagation of a single planar hydraulic fracture, followed by separate discussions on the modeling of crossing problem and of large-scale complex fracture networks. Planar hydraulic fracture models To simulate the propagation of a planar hydraulic fracture, whose length, width and height changes with time as fluid and proppant are injected into the fracture, the following modeling elements must be considered. Fracture opening displacement Most models make an assumption that the rock is elastic at large scale. For an elastic medium, the fracture opening width created due to deformation of the rock subjected to the fluid pressure

inside the fracture can be expressed in a general form as (Adachi et al., 2007)

pðx; yÞ  rmin ðx; yÞ ¼

Z

0

0

Cðx; y; x0 ; y0 Þ wðx0 ; y0 Þdx dy

ð1Þ

X

where  p(x, y) = fluid pressure  rmin (x, y) = minimum in-situ stress  w(x0 , y0 ) = fracture width  C(x, y, x0 , y0 ) = a complex stiffness function that contain the information of all the layers of the formation  X = the surface of the fracture. For a homogeneous rock, the elasticity equation is much simplified, with

Cðx; y; x0 ; y0 Þ ¼

G 1 4pð1  mÞ r3

where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  r ¼ ðx  x0 Þ2 þ ðy  y0 Þ2  G, m = shear modulus and respectively.

Poisson’s

ratio

of

the

rock,

An equation equivalent to Eq. (1) in boundary element formulation is given by Bui (1977):

pðx; yÞ þ rmin ðx; yÞ ¼

G 4pð1  mÞ    Z    @ 1 @w @ 1 @w 0 0 dx dy þ  0 0 @x @y r @y X @x r

ð2Þ

Therefore, the homogeneity approximation is made in many models to simplify the problem, and average moduli of the layered medium are used. The boundary integral formulation given above reduces the problem domain from 3D to 2D so the problem can be numerically solved much more efficiently. Fluid flow in the fracture Fluid flow in the fracture must satisfy the equation for conservation of mass, which reduces to Eq. (3) for an incompressible fracturing fluid:

@w @qx @qy þ þ þ ql ¼ 0 @t @x @y

ð3Þ

where qx and qy are components of flow rate per unit length in the fracture, ql is the fluid leakoff rate defined as

2ct ql ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t  sðx; yÞ

ð4Þ

ct is fluid leakoff coefficient, and s(x, y) is the time at which the fluid leakoff begins at location (x, y) on the fracture surface. Fluid flow must also satisfy the equation for conservation of momentum, which, when neglecting the inertia effect, reduces to the following relationship between the flow rates and pressure gradients (Ouyang, 1994; Yew, 1997):

"   2 #n1 2n 2 @p @p @p þ þ qg @x @y @x n1 " #  2  2 2n   2nþ1  n @p @p @p 1 w n k n nþ1 þ þ qg þ qg qy ¼  2n þ 1 @x @y @y 2n 2nþ1

n 1 w n qx ¼  k n nþ1 2n þ 1 2n

ð5Þ

where n and k are power-law and consistency indices of the fracturing fluid that in general has a power-law rheological behavior, q is fluid density, and g is gravity.

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

5

Fracture propagation criterion Along the fracture front, the stress intensity factor is equal to the rock critical stress intensity factor KIc,

K I ¼ K Ic

ð6Þ

While Eq. (6) is theoretically the valid condition at the fracture tip, the classic square root singularity at the fracture tip based on linear elastic fracture mechanics (LEFM) may only exist in a very small region near the tip, due to the large fluid pressure gradient in the fracture near the tip region and existence of a fluid lag between the fracture tip and fluid front. At a small distance into the fracture from the tip, the width profile no longer follows the square root distribution, since it is often dominated by coupled fluid flow effect rather than just elastic deformation. The near-tip pressure and width distribution can be described by a 2D KGDtype model (Geertsma and de Klerk, 1969). General 2D solutions have been obtained for various parametric limits, including stress intensity factor, viscosity, and leakoff-dominated regimes (Detournay et al., 2002; Adachi and Detournay, 2002; Garagash and Detournay, 2002; Detournay and Garagash, 2003; Detournay, 2004). Due to large grid size that is often used in the numerical models, significant error can be introduced if care is not taken to accurately model the tip region. Instead of using a very fine grid near the tip region, the tip asymptotic solutions based on 2D models are often used for the tip element. This allows an accurate solution to be obtained without expensive fine gridding. From the above equations, the fluid pressure p and width w within the planar fracture, and the fracture front evolution with time, can be numerically solved, in conjunction with appropriate boundary conditions along the fracture front and at the wellbore. Additionally, the proppant transport in the fracture can be solved from the transport equation (Yew, 1997):

@ðcwÞ þ r  ðcwv p Þ ¼ 0 @t

(a)

Planar 3D model using moving grid.

(b)

Planar 3D model using fixed grid.

ð7Þ

where c is the volume concentration of the proppant and vp is the velocity of proppant. The majority of commercial planar hydraulic fracture models can be categorized either as a planar 3D (PL3D) model, or a pseudo-3D (P3D) model. Fig. 3, taken from Adachi et al. (2007), illustrates various numerical schemes used for discretizing the fracture plane. Fig. 3a shows a moving boundary element mesh used by Clifton and Abou-Sayed (1981) and Clifton and Wang (1991) in their planar 3D model for one-half of the fracture plane (the fracture is assumed symmetrical with respect to the wellbore). The fracture plane is discretized into triangular elements using an automated grid generation scheme. As the fracture shape evolves at each time step, the mesh is regenerated. A similar moving mesh scheme was used by Gu (1987) and Gu and Yew (1988). The mesh evolution is based on a node convection scheme and adjusted as triangular elements become distorted or unevenly distributed. Fig. 3b shows a fixed grid used in the planar 3D model by Peirce and Siebrits (2001) and Siebrits and Peirce (2002). The fracture front is interpolated from the width solution obtained at the interior nodes of the grid at each time step based on the use of the asymptotic solution at the fracture front and mass balance. A fixed grid of equal grid size in x and y directions is used by Barree (1983), in which the fracture front is not specifically computed. The fracture element opens when fluid pressure at tip overcomes a ‘‘process zone stress’’ (Green et al., 2007). Besides the differences in the meshing scheme, the solution approaches are also quite different among the different PL3D models, primarily in the solution of the elasticity equation. Clifton and Abou-Sayed (1981) solves the integral equation Eq. (2) using

(c)

(d)

Cell-based pseudo-3D model

Lumped elliptical pseudo-3D model.

Fig. 3. Various numerical schemes for planar 3D and pseudo-3D models.

expansion in Chebyshev polynomials. Gu (1987) solves the variational form of Eq. (2) using the finite element approach. Barree (1983) computes fracture opening using an approximate solution

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

6

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

based on superposition of the deformation of elastic half-space subjected to the pressure applied at each square element of the fracture plane. The above models solve the elasticity equation assuming a homogeneous medium or using averaged elastic moduli. However, when there is very large variation in the elastic moduli from layer to layer, error in the computed fracture width can be introduced (Smith et al., 2001). Siebrits and Peirce (2002) adopted an approach that solves the elasticity problem for a layered medium rigorously. The kernel function C in Eq. (1), which corresponds to the normal stress r(x, y) induced by a unit displacement at point (x0 , y0 ), is obtained by applying Fourier transform to the 3D elasticity equations in x and z (normal to fracture plane) directions, solving the 1D ordinary differential equation of the transformed quantities, and then applying the inverse transform to obtain the kernel function. The boundary integral Eq. (1) is discretized and solved numerically. Although planar 3D models provide the most accurate solutions, they were computationally expensive and simulations typically took a long time to run on a personal computer. Only in recent years with the significant improvement in computer hardware, have planar 3D models become more widely used. For this reason, simpler P3D models have historically been more widely used and adopted in most commercial fracture design software. P3D models are computationally much faster than planar 3D models, well suited for day-to-day fracture design, simulation, post-treatment analysis, and design optimization. Fig. 3c and d illustrate two types of P3D models, cell-based and lumped models. In a cell-based P3D model (Settari and Cleary, 1982; Mack et al., 1992), the fracture is divided into a number of elements, or cells, along the length of the fracture. Each cell can have different height, as shown in Fig. 3c, and may cover multiple formation layers. 2D plane strain assumption is made for each vertical cross section. With this assumption, the fracture opening width is directly related to the local fluid pressure and decoupled from the pressure in the neighboring cells, greatly simplifying the elasticity equation. Neglecting vertical flow, the stress intensity factors at fracture top and bottom tips KIu and KIl and fracture width w can be directly computed analytically for a layered medium with piecewise constant stress (but not layered moduli), as given below (Fung et al., 1987; Mack and Warpinski, 2000):

rffiffiffiffiffiffi   ph 3 pcp  rn þ qf g hcp  h 2 4 rffiffiffiffiffiffi n1    qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2X h h  2hi þ ðriþ1  ri Þ arc cos  hi ðh  hi Þ 2 ph i¼1 h rffiffiffiffiffiffi   ph h K Il ¼ pcp  rn þ qf g hcp  2 4 rffiffiffiffiffiffi n1    qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 2 h h  2hi þ hi ðh  hi Þ ð8Þ ðriþ1  ri Þ arc cos þ 2 ph i¼1 h

K Iu ¼

  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 h z wðzÞ ¼ 0 pcp  rn þ qf g hcp   zðh  zÞ 4 2 E h2h  2 3 i þh 1 z i h n1  zÞ cos h ðh i 4 X 6 7 jzhi j ðr  ri Þ4 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ 0

5 pE i¼1 iþ1 i þ zðh  zÞarc cos h2h h

ð9Þ

where hi is the distance from top of the ith layer to fracture bottom tip, pcp is the fluid pressure at a reference (entry) depth hcp measured from the bottom tip, and qf is fluid density. The fracture height can be determined at each position of the fracture by matching KIu and KIl, given by Eq. (8), to the fracture toughness KIc of the corresponding layer containing the tips. The fracture height computed directly using Eqs. (8) and (9) is referred to as equilibrium

height. When vertical fluid flow effect is considered, the fracture height would be less than the equilibrium height at a given reference pressure. The fluid flow effect can be computed numerically (Settari and Cleary, 1982; Weng, 1991) or by using the analytical tip asymptotic solution. In a lumped model (Cleary, 1980; Cleary et al., 1983), the governing equations are nondimensionalized, with spatial variation in the fracture represented by integral coefficients. These coefficients can be either computed with simplifications with self-similar assumptions or specified, e.g., based on the results from numerical simulation such as the cell-based model discussed above. Meyer (1989) presented a lumped model in which the fracture plane is approximated by an upper and a lower half ellipse as shown in Fig. 3d. The governing equations reduce to a set of ordinary differential equations in time in terms of variables including wellbore pressure, width, half-length, upper and lower heights, etc., by assuming certain analytical functional forms for the variables in the fracture and carrying out the integrals along the fracture length. The lumped models reduce the order of governing equations and can be solved very efficiently. In Eqs. (8) and (9), a homogeneous medium is assumed which, as mentioned above, could introduce errors for a layered formation with large modulus contrast between layers. Siebrits et al. (2001) developed an enhanced P3D model using a rigorous 2D solution for a layered-modulus medium by adopting the Fourier transform technique as in the PL3D model discussed earlier. The effect of modulus contrast on fracture height containment is investigated by Gu and Siebrits (2008). Enhanced P3D model using the finite element approach for layered modulus is also available in a commercial design software. The above presents a brief overview of various conventional planar 3D and P3D fracture models. For more comprehensive reviews of these models, refer to Mack and Warpinski (2000) and Adachi et al. (2007). As may be evident from the discussions above, almost all planar hydraulic fracture models are based on the boundary element approach by reducing the elasticity equation to an integral equation and discretizing only the fracture surface, hence achieving computational efficiency. The conventional finite element based approach has seen only limited use in simulation of hydraulic fracturing problems, mainly because of the challenge of requiring a very fine mesh near the fracture front in order to achieve sufficient accuracy and having to constantly remesh as the fracture front evolves with time; both operations are computationally expensive and very time consuming. However, the boundary element approach also becomes much more complex when the fracture becomes nonplanar and/or when many fracture planes coexist and intersect as is the case for complex fractures.

HF–NF interaction As discussed earlier, the interaction between HF and NF plays a critical role in creating fracture complexity during hydraulic fracturing treatments in formations with preexisting NF. Understanding and proper modeling of the mechanisms governing HF–NF interaction are keys to explain fracture complexity and microseismic events observed during hydraulic fracturing treatments, and, ultimately, to be able to properly predict fracture geometry and reservoir production. When a HF intersects a NF it can cross or be arrested by the NF, and can subsequently open (dilate) the NF, as discussed earlier. If the HF crosses the NF, it remains planar, with the possibility to open the intersected NF if the fluid pressure at the intersection exceeds the effective stress acting on the NF. If the HF does not cross the NF, it can dilate and eventually propagate into the NF, which leads to a more complex fracture network. So, the crossing

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

7

criterion, in general, strongly influences the complexity of the resulting fracture network. The interaction between the HF and the NF depends on the insitu stresses, mechanical properties of the rock, properties of natural fractures, and the hydraulic fracture treatment parameters, including fracturing fluid properties and injection rate. During the last decades, extensive theoretical and experimental work has been done to investigate, explain, and develop the rules controlling HF–NF interaction (Warpinski and Teufel, 1987; Blanton, 1982, 1986; Renshaw and Pollard, 1995; Hanson et al., 1982; Beugelsdijk et al., 2000; Potluri et al., 2005; Zhao et al., 2008; Gu and Weng, 2010; Gu et al., 2011) as well as numerical simulations to model fracture behaviors at the intersection (Heuze et al., 1990; Zhang and Jeffrey, 2006, 2008; Thiercelin and Makkhyu, 2007; Zhang et al., 2007a, 2007b, 2009; Zhao and Young, 2009; Chuprakov et al., 2010, 2013a; Meng and de Pater, 2010; Dahi-Taleghani and Olson, 2011; Sesetty and Ghassemi, 2012). Analytical models Based on relatively simple analytical considerations of the stress field at the propagating fracture tip, analytical crossing criteria have been developed (Warpinski and Teufel, 1987; Blanton, 1986; Renshaw and Pollard, 1995; Gu and Weng, 2010). Renshaw and Pollard (1995) developed a simple criterion for predicting if a fracture will propagate across a frictional interface orthogonal to the fracture. The criterion is based on the linear elasticity fracture mechanics solution for the stresses near the fracture tip. It determines the stresses required to prevent slip along the interface at the moment when the stress on the opposite side of the interface is sufficient to reinitiate a fracture. The Renshaw and Pollard crossing criterion is given as

rH T 0 þ rh

>

0:35 þ 0:35 k 1:06

ð10Þ

where rH and rh are maximum and minimum horizontal stresses, respectively, T0 is the rock tensile strength, and k is the coefficient of friction of the NF interface. The criterion was validated by laboratory experiments for dry cracks. Natural fractures are often not aligned with the contemporary principal in-situ stress directions in the rock formation, in which case the intersection angle of a HF approaching a NF is between 0° and 90°. Because the intersection angle has a significant effect on crossing, the Renshaw and Pollard criterion has been extended to fracture intersection at nonorthogonal angles by Gu and Weng (2010). It shows that it becomes more difficult to cross an interface when the intersection angle decreases from 90°. This crossing criterion was compared to the laboratory results from block experiments and good agreement was obtained (Gu et al., 2011; also see Fig. 4). With their relative simplicity, these analytical criteria do not take into account the detailed evolution of stress field and rock deformation after the initial contact by the tip. As a result, they are insensitive to parameters of fluid injection affecting the hydraulic fracture geometry and the fluid infiltration into the natural fracture after the contact. These criteria capture the first-order effects of fracture intersection angle, NF friction coefficient, and in-situ stresses. However, field and laboratory observations showed that fluid properties are important and should also be accounted for. The experimental study by Beugelsdijk et al. (2000) showed that flow rate and fracturing fluid viscosity strongly influence hydraulic fracture complexity in a prefractured block. With low value of the product of the injection rate Q and fracturing fluid viscosity l (Ql), fluid tends to leak into the preexisting

Fig. 4. Comparison of laboratory ‘‘crossing-arresting’’ data with previous analytical models (Blanton, eRP), new analytical model OpenT, and MineHF2D simulations (Chuprakov et al., 2013a).

discontinuities in the rock and creates tortuous fracture paths following the discontinuities. With a large Ql value the hydraulic fracture tends to cross most discontinuities and the overall fracture path is nearly straight. A similar observation is made based on the microseismic monitoring of a treatment using gel fracturing fluid and a treatment using slick water in the same well in Barnett shale (Warpinski et al., 2005). To improve the description of HF-NF interaction, a much more sophisticated analytical model that takes into account the mechanical influence of the HF opening and the hydraulic permeability of the NF, was recently developed by Chuprakov et al. (2013a). The model (referred to as the OpenT model) solves the problem of elastic perturbation of the NF at the contact with a blunted HF tip, which is represented by a uniformly open slot (thereby, giving it the name OpenT). The opening of the HF at the junction point wT (blunted tip) develops soon after contact, and approaches the value  defined by the of the average opening of the hydraulic fracture w; injection rate Q and the fluid viscosity l. The model computes the elastic stress field along the interface and in the vicinity of the activated NF, and determines the size of open and shear sliding zones. The analysis of the generated stress field gives the sites of sufficient tensile stress concentration where the new fractures can be nucleated. A combined stress and energy criterion determines if the fracture can be initiated. If the fracture initiation criterion is satisfied, a secondary fracture is initiated on the opposite side of the interface, leading to either direct crossing (if the initiation site is right at the intersection point) or crossing with an offset (if initiation site is located away from the intersection point). The solution shows that the spatial extent of the open and sliding zones strongly depends on the fluid pressure inside the activated part of the NF. The larger the inner fluid pressure, the larger the open and sliding zones at the NF. Consequently, it is expected that after the HF contacts the NF, the injected fracturing fluid will gradually penetrate the NF with finite hydraulic permeability j and thus enhance the fluid pressure within the NF. The average pressure inside the NF is computed approximately based on the fluid diffusion into the fracture. As a result, the fluid pressure builds up much more easily in a highly permeable NF or with a low-viscosity fracturing fluid, leading to greater likelihood of HF being arrested. Conversely, it is more difficult to penetrate the NF with a high-viscosity fluid, leading to greater likelihood of HF crossing the NF.

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

8

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

The predictions of the analytical models for crossing behavior generally agree with the experimental results for various fracture intersection angles, stress contrasts. Fig. 4 shows the comparison between different analytical models by Blanton (1986), Gu and Weng (2010, also referred to as extended Renshaw and Pollard, or eRP model in Fig. 4), OpenT model (Chuprakov et al., 2013a), and the experimental results from Gu et al. (2011), as well as numerical simulation using the MineHF2D code (Zhang and Jeffrey, 2006, 2008; Zhang et al., 2007b). It shows that both the analytical and the numerical models agree with the experiments, suggesting they capture the first-order crossing-arresting behavior. Note that the discrimination between the different models would require additional data points in the transitional region near the crossing-noncrossing boundary. Additionally, the injection rate and viscosity were not changed in this series of experiments, so it was not possible to assess their effect on the fracture interaction outcome. Future experimental verification aimed to further narrow the probed parameter range in the transitional region and also incorporating variation of fluid properties and flow rate is highly desirable, though these experimental programs tend to be very expensive. To compensate for this lack of laboratory data for full validation, numerical experiments were conducted using MineHF2D code to assess the sensitivity of the injection rate on fracture crossing. The results for one set of parametric runs are shown in Fig. 5. The OpenT model agrees well with the numerical results in the sense that it captures the crossing-arresting transition as flow rate changes. The model has been implemented in a P3D-based complex fracture model called the UFM model which is to be discussed later in this paper. Numerical models For more accurate modeling of HF–NF interaction, numerical models are often employed. Zhang et al. (2009) give a generalized 2D elasticity equation for arbitrary fracture surfaces X in an infinite elastic medium:

rn ðxÞ  r1 n ðxÞ ¼ sðxÞ  s1 ¼

Z

Z

½G11 ðx; sÞ wðsÞ þ G12 ðx; sÞ v ðsÞds

X

ð11Þ

½G21 ðx; sÞ wðsÞ þ G22 ðx; sÞ v ðsÞds

X

where x = (x, y) is a point on the fracture surface X, ds is an infinitesimal length increment along the fracture, rn(x) and s(x) are the 1 normal and shear stresses acting on the fracture face, r1 are n and s

the normal and shear stress induced by the remote in-situ stresses on the fracture at point x, w and v are opening and shear displacements along the fracture, and Gij are hyper-singular Green’s functions given by Zhang et al. (2005). Note that Eq. (11) not only contains the normal components of the stress and displacement as in the planar equation Eq. (1), it also contains the shear components. It applies not only to the open hydraulic fracture, where rn(x) is equal to p(x) and s(x) = 0, but also to the closed fracture where shear stress s(x) may be nonzero. The mass conservation equation (for incompressible fluid) can be written as

@q @ðw þ -Þ þ þ qL ¼ 0 @s @t

ð12Þ

where - is aperture of closed fracture due to surface roughness, and the equation for fluid flow in the fracture becomes

 n @p 4n þ 2 qn ¼ 2k @s n ðw þ -Þ2nþ1

ð13Þ

For the part of the fracture system that is closed, the Coulomb frictional law applies, which gives

jsj 6 kðrn  pÞ

ð14Þ

The aperture of a closed fracture, x, is a function of both effective normal stress and shear displacement. When the shear stress acting on the fracture interface reaches the maximum shear dictated by Coulomb frictional law (equality in Eq. (14)), the interface slips. This shear slip causes dilation of the interface, i.e., an increase in x and consequently its conductivity to fluid flow. In the Zhang et al. formulation, the mechanical aperture that gives rise to an increase in porosity as implied in Eq. (12) and the hydraulic aperture that dictates fracture conductivity are considered the same for simplicity. These two quantities can be significantly different, especially when mechanical aperture is small. The shear slip induced dilation is also not considered. The description of natural fracture dilation will not be given here; refer to Zhang et al. (2009) and Yew and Weng (2014) for further detailed discussions on this subject. Dividing the fractures into small equal-size elements, and discretizing the governing equations by using the displacement discontinuity method (DDM) (Crouch and Starfield, 1983), the resulting coupled nonlinear system of equations can be solved for pressure, opening width, and stresses along the fractures. Using this model, Zhang and Jeffrey (2006, 2008), Zhang et al. (2007a, 2007b, 2009), and Jeffrey et al. (2009) investigated various

Fig. 5. Comparison of HF–NF crossing-arresting behavior between analytical models and numerical model using MineHF2D code. The red crosses and squares respectively indicate crossing and arresting behavior from MineHF2D code, solid green curves correspond to analytical predictions using the OpenT criterion, dashed yellow curve corresponds to the Blanton criterion, and the eRP criterion is given by dashed blue lines. The interaction is studied for various injection rate and relative stress difference for two different HF-NF contact angles, b = 90° (left) and b = 60° (right) (Chuprakov et al., 2013a). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

problems related to a hydraulic fracture crossing a natural fracture or a bedding plane, effect of fracture offset, and interaction of multiple fractures reinitiated from the interface. In Zhang et al. (2007a), the effect of fracture arrest at a formation interface on the pressure response was studied, for both situations of stiff-tosoft and soft-to-stiff interfaces. The effect of fracture offset on the fluid pressure and evolution of fracture width were studied in Zhang et al. (2007a, 2009) and Jeffrey et al. (2009). The offsets cause increase in the fluid pressure, due to higher frictional pressure drop through the offsets that are subjected to higher normal stress and have narrow width. Cooke and Underwood (2001) numerically investigated three types of fracture intersection with bedding contacts: fracture propagation through bed contacts, termination (abutment) at contacts, and step-over of fractures at bedding contact. Their numerical code is also based on 2D DDM, and results suggest that local interface opening near the tip of approaching fractures, rather than sliding, is responsible for fracture termination and step-over at bedding contacts. Similarly, Chuprakov et al. (2010) used the 2D DDM method (without full fluid coupling) to study the detailed stress distribution in the rock on the opposite side of the intersecting hydraulic fracture tip. The results showed high stress concentration at the intersection point as well as at the tip of the opening zone along the natural fracture near the intersection point. Sesetty and Ghassemi (2012) also used the 2D DDM method coupled with fluid flow to simulate the interaction of a hydraulic fracture with a natural fracture and the pressure response, as well as fracture growth after reaching the ends of natural fracture. A different numerical approach was adopted by Dahi-Taleghani and Olson (2011). They used the 2D extended finite element method (XFEM) to simulate the HF growth when intersecting a NF. They investigated the effect of fracture toughness of the cementing material in the NF relative to the rock matrix on the crossing and arrest behaviors. They also showed asymmetry growth in the NF when the HF intersects the NF at an angle. Whereas the numerical models have been successfully used to simulate the detailed stress field and obtaining the insights to the complex process of hydraulic fracture–natural fracture interaction, these simulations are quite computationally intensive, even though the problem is only 2D, because of the fine numerical elements needed to produce accurate results and the complexity of the coupled solid–fluid problem. It is difficult to apply these models to simulate the reservoir-scale problems that have much larger length scale, with 3D geometry (having finite fracture height), and possibly involving a large number of fractures. In contrast, the analytical crossing models discussed earlier can be computed very fast, which make them much more suited for integration in large-scale complex fracture models.

Modeling of large-scale complex fractures In recent years, new hydraulic fracture models have been developed, or existing geomechanics models adapted, for simulating induced complex fracture networks in a hydraulic fracture treatment of an unconventional reservoir (Weng et al., 2011; Xu et al., 2009; Meyer and Bazan, 2011; Nagel et al., 2011; Fu et al., 2011; Dershowitz et al., 2010; Savitski et al., 2013; Wu and Olson, 2013; McClure, 2012). As discussed in the previous sections, the process of hydraulic fracture propagation in a formation with preexisting natural fractures is very complex. Simulation of this process requires proper consideration of the key physical elements governing the process that are important for ultimate reservoir production, including rock deformation, fracture propagation, fluid flow in the complex fracture networks, interaction between the hydraulic fractures and natural fractures, interaction among

9

different hydraulic fractures, fracture height growth, and proppant transport in the fracture networks. Although complete and rigorous simulation of all the aspects is technically very challenging and most models have simplifying assumptions, it is important that the models capture the most essential elements so that simulation reasonably represents the real process. Today, many complex fracture models are still evolving and are applicable only under the limited conditions where the model assumptions are valid or are missing important elements to fully simulate the complete fracturing process. A critical consideration of a complex fracture model is the interaction between the hydraulic fracture and the natural fracture. For a formation that initially contains a large number of well-connected, highly permeable natural fractures, and is subjected to low in-situ stress anisotropy, fracturing fluid has a tendency to follow the existing natural fracture network. In that case, the induced hydraulic fractures follow the natural fracture network, and relatively few new fracture paths are created. Under this scenario, modeling of fracture propagation and intersection with natural fractures are not as critical. A static numerical grid properly modeling the natural fracture system can be used to simulate the problem. A coupled geomechanics-reservoir model that is capable of solving the coupled fracture (or joint) deformation and fluid flow in the fracture networks can be well suited for this kind of simulation. On the other hand, in formations that contain a large number of isolated or poorly connected natural fractures, the natural fractures are only hydraulically activated when intersected by hydraulic fractures, the consequence of which is to alter the path or cause branching of the hydraulic fractures. To simulate the fracture treatment in this type of reservoir, the model needs to properly simulate hydraulic fracture propagation and interaction with the natural fractures. In the following, various modeling approaches to simulate large-scale complex fracture networks are discussed. 2D Boundary element based models The 2D models already discussed in the previous section can theoretically be used to simulate large-scale fracture networks. However, it becomes challenging when it comes to selecting the proper grid size, one that allows the simulator to model the reservoir-scale problem and the potentially large number of fractures, while still not compromising the key physical process. As discussed before, if the fracturing process to be simulated involves interaction between hydraulic and natural fractures, the grid size required is very small in order to accurately simulate the crossing process, which would result in very large number of elements needed to simulate the propagation of a full-scale complex fracture network. Consequently, it limits the applicability of the models for such simulations. Olson (2008) and Olson and Dahi-Taleghani (2009) used an enhanced 2D DDM method that is based on the classic 2D DDM approach (Crouch and Starfield, 1983) but with a correction factor to account for the effect of finite fracture height to simulate complex fracture pattern development. The model was based on a model developed earlier by Olson (1993, 2004) for predicting natural fracture formation over geological time as a result of subcritical crack growth under in-situ stresses. The model does not specifically couple the fluid flow in the fracture network. Instead, a static pressure is assumed, and hence the model is very efficient and can simulate the growth of very large number of fractures and complex fracture network patterns. The model can be useful as a research tool for obtaining insights into the fracture interactions and pattern development, but cannot be used as a fracture design simulator since fluid flow is not considered.

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

10

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

simulate fractures with finite height. The model is used to study the effect of stress shadowing on the fracture growth pattern among fractures propagating from multiple perforation clusters in a horizontal well. One general limitation of 2D models is not being able to simulate fracture height. This is not only an important consideration in fracture designs, but can also drastically affect fracture length growth. Therefore, a general-purpose fracture design model needs to be able to model fracture height growth.

3D boundary element based models

Fig. 6. An Example of a fracture network with potential hydraulic fractures by McClure (2012). The black line is the wellbore; the blue lines are preexisting fractures, and the red lines are the potentially forming hydraulic fractures. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

McClure (2012) developed a comprehensive 2D DDM based model fully coupled with fluid flow, stresses induced by fracture opening and sliding, friction evolution, and fracture propagation in a pre-existing discrete fracture network (DFN). An efficient numerical scheme is developed to allow simulation of very complex fracture networks. Four different stimulation mechanisms, including pure opening mode, pure shear stimulation, mixedmechanism stimulation, and primary fracture with shear stimulation, are simulated and investigated. Fig. 6 shows a simulation of an induced hydraulic fracture network in a formation with preexisting natural fractures by McClure. One drawback in McClure’s approach is that the path of any new tensile fracture, e.g., generated from interaction with natural fractures, must be predefined. In other words, certain crossing behaviors at NFs are predefined, rather than predicted. The locations of potential new fractures are predefined and discretized before the simulation, which may be activated during the simulation. Wu and Olson (2013) developed a model based on the extended 2D DDM approach coupled with fluid flow to more realistically

When there are multiple hydraulic fractures propagating simultaneously, or when fracture surface deviates from a plane, the problem becomes fully 3D, and the corresponding equations are much more complicated to solve. The 3D boundary element approach extends the elasticity equation as given earlier in Eq. (1) for planar 3D to nonplanar fractures or multiple isolated or interconnected fracture surfaces. The flow equation for the planar 3D still applies but is based on the local coordinates along the fracture surface. Fracture front growth becomes much more complicated since the fracture tip may be subjected to not only mode I, but possibly also mode II and III stress singularities. Lam et al. (1986) presented one of the earliest fully 3D hydraulic fracture models and used it to simulate a curved fracture development and compared it to the laboratory experiment. A similar model was developed by Vandamme et al. (1988). Sousa et al. (1993) and Carter et al. (2000) developed a 3D fracture model that had been used for studying fracture initiation from a deviated wellbore and its interaction with casing and borehole. Yamamoto et al. (2004) developed a full 3D model to simulate interaction of multiple nonplanar fractures (see Fig. 7a). All the models mentioned are based on the boundary element approach to achieve computational efficiency by solving the integral equation over the fracture surfaces. The fracture surface is typically discretized into a triangular or quadrilateral mesh that is adjusted as fracture front evolves. However, even with the boundary element scheme, the numerical solution is extremely time consuming and is often plagued by numerical instability. Therefore, these models are only used as research tools rather than for engineering design purposes.

Fig. 7. Multiple fractures simulated using 3D boundary element based models. (a) Yamamoto et al. (2004); (b) Xu and Wong (2013).

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

11

Fig. 8. Multiple initially parallel fractures from a horizontal well simulated by Castonguay et al. (2013): (a) three fractures, isotropic stress, equal wellbore pressure, 0.25 Pa-s viscosity; (b) same condition but equal injection rate; (c) four fractures, higher vertical stress.

In recent years, with the multistage fracturing in horizontal wells becoming the most preferred completion method in unconventional reservoirs, interest in simulating multiple, simultaneously growing fractures and stress shadowing revived new development in more robust 3D models. Xu and Wong (2013) and Wong et al. (2013) presented a 3D model and its use for simulating multiple fractures propagating from perforation clusters and effect of stress shadowing on sequentially placed fractures. Fig. 7b shows four initially parallel fractures propagating from a horizontal well, with fractures contained in height, by Xu and Wong (2013). Rungamornrat et al. (2005) developed a 3D model for nonplanar fractures, which has been extended by Castonguay et al. (2013) to simulate multiple initially parallel fractures from a horizontal well. Fig. 8a and b show an example by Castonguay et al. (2013) of three fractures in an isotropic stress field under equal wellbore pressure and equal flow rate conditions, respectively. Fig. 8c shows four fractures subjected to isotropic horizontal stresses but a higher vertical stress. These examples demonstrated several important effects of stress shadowing that influence the fracture geometry. In Fig. 8a, under equal wellbore pressure condition, the width of middle fracture is suppressed by the outer fractures, leading to the outer fractures having greater length than the middle fracture. Under the same formation condition, with equal flow rate injected into each fracture, all fractures have approximately equal dimension (Fig. 8b). This is because they must contain equal volumes of fluid volume injected. This case represents the case of limited-entry completion in which a limited number of perforations are shot in each perforation cluster, resulting in very high perforation friction that forces approximately equal flow rate distribution among the perforation clusters. Fig. 8c shows that when closely spaced fractures are forced to propagate parallel to each other, the system is inherently unstable, resulting in ‘‘petals.’’ Fractures grow in an asymmetric way to avoid each other, attaining a lower energy state with wider fracture width, as opposed to parallel fractures that have narrower width and require higher energy to propagate. While fully 3D boundary element based models are capable of simulating complex nonplanar fractures, they are still limited in

their ability to efficiently simulate interaction with natural fractures and complex fracture networks. Future improvements in these models are expected to extend their use in the field applications. Coupled geomechanics and reservoir models Geomechanics models are widely used in the oil industry in a variety of problems involving drilling, well construction, fracture stimulation, reservoir compaction and subsidence, sanding, fault activation, etc. (see, e.g., Koutsabeloulis and Hope, 1998; Koutsabeloulis and Zhang, 2009). Related to hydraulic fracturing design, geomechanics models can be used to compute the initial stress state in the formation of interest, based on the seismic and log-derived geological surfaces, faults, and formation properties, as well as knowledge of regional tectonic strains. They can also be used to compute the evolution of stress state resulting from hydrocarbon production and reservoir pressure depletion, which can change both magnitude and orientation of the stresses and is an important consideration especially for wells drilled in a mature field. Geomechanics models can be coupled with reservoir simulators to take into account the effect of the stress changes on the rock properties and reservoir transmissibility which impacts production or water injection. A geomechanics model is typically based on finite element or finite difference methods. It solves the equations governing the deformation and stresses of a poroelastic or poroplastic medium. By coupling with reservoir simulation, the pore pressure change in the reservoir changes the effective stress acting on the rock matrix, which in turn affects the permeability and porosity of the rock matrix and natural fractures, and subsequently reservoir fluid flow. The effect of pore pressure change on permeability is relatively weak in reservoir engineering applications where pore pressure is below the minimum in-situ stress. This weak dependency allows the coupling be accomplished in the numerical algorithms more easily. However, in the case of hydraulic fracturing, where the fluid pressure exceeds the minimum insitu stress and creates tensile fractures, the fluid transmissibility

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

12

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

changes drastically for a small change in pressure, since the fracture aperture changes linearly with net pressure (fluid pressure minus minimum in-situ stress), while fracture conductivity increases rapidly with fracture aperture according to the cubic law. This presents a major challenge when using conventional geomechanics models for simulating hydraulic fractures. Special algorithmic treatments are often necessary to achieve convergence and stability. Ji et al. (2009) presented a model that couples a 3D finite element geomechanics model with a conventional finite difference reservoir flow simulator. The geomechanics module implicitly models fracture propagation via displacements on the fracture face. The flow and geomechanics models are coupled in an iterative manner. The 3D planar fracture geometry and the pressures in it are the common dynamic boundary conditions for the flow and stress modules. The iterative process yielded smooth fracture propagation which was validated against the analytical 2D models. The coupled geomechanics and reservoir model allows simulation of fracture propagation, post-fracturing multiphase cleanup, and well performance in a changing stress and pressure environment, all within the same system. The model is applied in a field example to simulate complex injection history and calibrate the stress dependency of formation permeability. Dean and Schmidt (2008) presented a fully coupled geomechanical reservoir simulator, in which the poroelastic/poroplastic equations are discretized over the 3D reservoir grid. The geomechanical equations are solved together with the reservoir flow equations using a Newton–Raphson technique in which a Jacobian is generated for the entire system of equations to increase robustness and convergence of the nonlinear iterations. The model explicitly simulates fracture propagation, and fracture growth computation can be based on either critical stress intensity factors or use of cohesive elements that exhibit strain-softening behavior. Chen (2012) used a standard finite element modeling (FEM) code, ABAQUS, to simulate the rock deformation and incorporate separate fluid flow calculations for the fluid pressure inside the fracture as a traction applied on the fracture surfaces in the FEM grid. A cohesive-zone fracture tip model is used instead of a stress intensity based tip model. A special meshing scheme with fine grids in the cohesive zone and adjacent elements is used to obtain sufficient accuracy. The FEM simulation was compared to 2D KGD and PKN solutions and showed good agreement. Similarly, Shin and Sharma (2014) also used ABAQUS to model the reservoir as a poroelastic medium, coupled with fracture planes represented by pore pressure cohesive elements, and they simulated the initiation and propagation of multiple parallel fractures. The fractures are subjected to equal pressure at the wellbore; no perforation friction is considered. Similar to the equal pressure case by Castonguay et al. (2013), the results show the fracture interference caused the middle fracture to be suppressed by the outer fractures due to stress shadow effect. Singh et al. (2014) presented a coupled reservoir-fracture flow model that accounts for varying reservoir geometries and complexities including nonplanar fractures, faults, and barriers. It considers multiphase flows in a fractured poroelastic reservoir where fractures are modeled as surfaces. A multipoint flux mixed finite element (MFMFE) method is used for spatial discretization of the reservoir and fracture domains. An iterative coupled implicit pressure explicit saturation (IMPES) scheme is used for solving the system. The demarcation of reservoir and fracture as separate domains allows for special treatment of computationally challenging regions. Examples of water injection and production simulation in fractured reservoirs were presented. Fracture simulation based on a geomechanics model coupled with reservoir simulation is most versatile in terms of simulating complex reservoir structures and flow behaviors, as well as

heterogeneous and nonlinear mechanical behaviors of the rock such as plasticity. However, due to highly localized stress concentration, rock deformation, and fluid pressure gradient near fracture tips, a fine numerical mesh is generally required near the tip of the fracture (while coarser away from the tip to reduce the number of grids) to obtain a reasonably accurate solution. As fractures grow, remeshing is needed to move the fine mesh with the fracture tips. Use of special tip element will help reduce the grid size requirement but will require tracking the fracture tip movement as well as remeshing. An alternative is to put a fine mesh along the entire expected path of fracture growth, which significantly increases the number of grids. In general, the finite element or finite difference based approaches are much less efficient than the boundary element based approaches and require much longer computation time. Because of this, the extended finite element method (XFEM) has been explored as an alternative method that can potentially be much more efficient than standard FEM. In the XFEM approach, a fixed mesh is used. Discontinuities such as fractures are allowed to propagate independently of the mesh configuration. The elements occupied by the fractures are enriched by introducing additional shape functions that have the form of analytical or asymptotic fracture tip solutions. Doing so introduces additional degrees of freedom and complexity in the algorithm for the enrichment elements that represent the fractures, but eliminates the need for remeshing. Dahi-Taleghani and Olson (2011, 2013) used a 2D XFEM model to simulate a HF intersecting a NF, and fracture propagation in a medium containing many pre-existing natural fractures. Weber et al. (2013) presented a general formulation of XFEM and simulation of a 2D fracture using a linear elastic fracture mechanics tip solution as an enrichment function. Chen (2013) achieved an XFEM implementation using special user elements based on XFEM incorporated in ABAQUS and was able to obtain good agreement with the 2D KGD model. Gordeliy and Peirce (2013) presented an in-depth description of XFEM formulation using a hydraulic fracture asymptotic function as the enrichment function for the fracture tip elements and implicit level set schemes for determining the fracture front. Abbas et al. (2014) used the XFEM model to study the effect of bedding plane and fracture offset at the bedding plane on the vertical fracture height containment and fracture width. Application of XFEM in modeling hydraulic fracturing problems is still in its infancy, but it has shown great promise. Discrete fracture network and distinct element method based models For simulation of fluid flow in a reservoir containing a large number of natural fractures, two general approaches had been taken in reservoir engineering. One is based on the continuum approach and the other based on modeling fluid flow in discrete fracture network (DFN). Continuum reservoir simulation is based on the idea that a heterogeneous reservoir containing fractures can be represented by an equivalent homogeneous medium. Warren and Root (1963) proposed to describe the fractured reservoir with two types of porosity: primary or matrix porosity, and secondary or fracture porosity (dual porosity). The two systems are linked through a transfer function (depending on the shape and dimension of matrix blocks, transmissibility of the block, relative permeability, capillary pressure, etc.), which represents the exchange between them. The discrete fracture network approach simulates fluid flow in a reservoir containing discrete fractures, taking into account their size, shape and distribution, as well as their own flow properties. It starts with spatial statistics associated with a fracture network (fracture orientation, size, density, transmissivity, etc.) measured from cores, well tests, surface outcrop, etc. These statistics are used to generate realizations of fracture network with the same spatial properties. Reservoir flow through these fracture network

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

realizations is computed by discretizing the fracture network and solving a very large system of equations governing the flow in the discretized system. A DFN model generally needs extensive computational resources. Discrete fracture models, in which the fractures are represented individually, can be used both as standalone tools as well as for the evaluation of transfer functions for dual-porosity models. Such models can also be used in combination with the dual-porosity approach. There is a large body of literature related to simulation of fluid flow in naturally fractured formation, for both the continuum and discrete fractures approaches. However, many of these models either do not consider fracture permeability change with pore pressure or not to the extent of tensile opening of the fractures, which is an essential feature of hydraulic fracturing. Therefore, an extensive review of these models will not be given here. To model hydraulic fracturing at a pressure above the in-situ stress, the model must consider the geomechanical deformation of the fractures and its effect on fluid flow in response to the pressure changes. Some of the models that couple fracture deformation and elasticity and fluid flow can be found in Bagheri and Settari (2008), Bai et al. (2002), Sanyal et al. (2000), Hossain and Rahman (2008), Karimi-Fard et al. (2004), Kazemi et al. (2005), Koudina et al. (1998), Nakashima et al. (2000, 2001), Rijken (2005), Sarda et al. (2001), Willis-Richards et al. (1996) and Zhang et al. (2009). Much of the earlier interest in modeling of injection into fractured reservoirs stemmed from the study of enhanced geothermal systems (EGS), under the premise that the behavior of EGS will be dominated by fracture flow. The paper of Sanyal et al. (2000) reviewed the special features that would be required of any practical numerical simulator for EGS. These features, required in addition to the basic features of conventional geothermal simulators (ability to handle two-phase fluid flow, heat transfer, and tracer transport in porous and fractured media) are: explicit representation of the fractures, change in fracture aperture due to effective stress and shear, thermo-elastic effects, relation between fracture aperture and conductivity, and channeling of fluid flow within fractures. In EGS systems, the main challenges are improving permeability through enhancement of natural fractures or creation of induced fractures, and optimizing heat recovery through injection. For simulating fluid flow and rock deformation in a fractured system, the discrete element method (DEM) is widely used. The discrete element method is a family of numerical methods for computing the motion of a large number of particles like molecules or grains of sand. One of the branches of DEM is termed as distinct element method by Cundall (1971, 1988) and Cundall and Strack (1979). In this method, the rock mass is represented as an assembly of discrete blocks connected by joints or interfaces. It assumes

13

that the positions of major faults and fractures are predefined and then solves the resulting system as a set of interacting cracks and blocks that is elastic or plastically deformable and moves according to Newtonian mechanics. This relatively robust explicit method has proved to be able to solve a variety of geotechnical problems. A similar method called discontinuous deformation analysis (DDA) has also been used for modeling the fractured rock medium (Shi, 1988; Jing et al., 2001). The DEM model has been used to analyze the injection treatment and interpret microseismic observations in hot dry rock (HDR) as an EGS by Murphy and Fehler (1986). Similar to the observations in injections into shale formations, a wide cloud of microseismic events was detected. Murphy and Fehler attributed the microseismic response to shear slip along the joints in the HDR. In recent years, a number of studies were carried out using DEM to simulate hydraulic fracturing injection into a naturally fractured medium. Nagel et al. (2011, 2012) and Savitski et al. (2013) carried out simulations of a planar hydraulic fracture cutting through a pre-existing DFN, and investigated the extent of shear failure in the natural fractures system as a function of various fracturing and rock parameters. As expected, the induced shear failure in the rock joints is highly dependent on the friction angle. The extent of shear failure can be very large for low friction angles (5° and 15°), but becomes more localized near the propagating fracture tip for friction angle of 30° (Nagel et al., 2012). Savitski et al. (2013) studied the effect of DFN on fluid loss from the main hydraulic fracture, as well as fracture apertures of the hydraulic and natural fractures and treating pressure. Fig. 9 shows the simulated pressure distribution in the DFN for a well-connected and a sparsely connected DFN, respectively. In the parametric study shown, the natural fractures were initially conductive, which allows fluid pressure to diffuse into the DFN and causes dilation. However, observation of shale cores suggests that most natural fractures in many shales are mineralized and few are open (Gale et al., 2007). This raises a question of whether an initially mostly nonconductive fracture systems can develop into a more connected and conductive DFN during stimulation. Savitski et al. further evaluated the potential for shearing of nonconductive fractures and concluded that high pressure fluid percolating through the DFN is the major cause of shear failure, and the stress perturbation alone accounts for only a small fraction of the failure. Riahi and Damjanac (2013) used a 2D DEM model to study the DFN connectivity on the induced hydraulic fracture network. Fig. 10 shows the predicted hydraulic fracture network for various DFN patterns. The results seem to suggest that well-connected DFN result in a longer hydraulic fracture system. However, this may be biased by a

Fig. 9. Fluid pressure distribution on a horizontal cross-section cutting through the injection point after 50 min of injection for two DFN realizations, by Savitski et al. (2013). A hydraulic fracture is located at the center of the network in the east–west direction.

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

14

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

Fig. 10. Effect of DFN connectivity, by Riahi and Damjanac (2013). The maximum horizontal stress is in the east–west direction.

major limitation of the DEM model, that is, the fluid can only flows along the predefined connected DFN paths. The model does not simulate creation of new tensile fractures through the rock blocks if the fracture paths are not predefined as block boundaries, artificially restricting the growth of hydraulic fractures. An enhanced DEM model that contains an FEM grid embedded within the DEM blocks was presented by Fu et al. (2011). The embedded FEM analysis allows initiation of new fractures and propagation of fractures through the DEM blocks. Fig. 11 shows a simulation using the DEM-FEM model.

While the DEM-FEM method overcomes the limitation of DEM of not being able to create new fracture paths, it is still quite challenging for it to be able to predict fracture initiation accurately. As mentioned before, the interaction between a hydraulic fracture tip and a natural fracture interface is a much localized process. To model it accurately requires a fine numerical grid, which is very difficult to incorporate in a large-scale problem, without requiring an enormous amount of computational resources. Recently, a numerical approach using the synthetic rock mass (SRM) concept has been developed (Pierce et al., 2007) for

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

15

Fig. 11. The original fracture network and the stimulated fracture network at three levels of stress anisotropy, by Fu et al. (2011).

simulation of fractured rock mass. The SRM consists of two components: (1) the bonded particle model (BPM) to model deformation and fracturing of intact rock and (2) the smooth joint model (SJM) to model mechanical behavior of discontinuities. Damjanac et al. (2013) presented a hydraulic fracture simulator based on implementation of the SRM in the lattice, a quasi-random array of nodes (with given masses) in 3D connected by springs. The lattice nodes are connected by two springs, one representing the normal and the other shear contact stiffness. The springs represent elasticity of the rock mass. The calibration factors for spring stiffness are built in the model based on the user-provided macroscopic elastic properties. The tensile and shear strengths of the springs control the macroscopic strength of the lattice. The fluid flow occurs through the network of pipes that connect fluid elements, located at the centers of either broken springs or springs that represent pre-existing joints (i.e., springs intersected by the surfaces of pre-existing joints). The mechanical and flow models are fully coupled through the following conditions: 1. Fracture permeability depends on aperture or on the deformation of the solid model. 2. Fluid pressure affects both deformation and the strength of the solid model. 3. The deformation of the solid model affects the fluid pressures.

Fig. 12 shows a simulation of hydraulic fractures propagating from two perforation clusters in a horizontal well.

Lumped elliptic P3D network models Even though 3D models presented above are capable of simulating more complex fracture geometry or fluid interaction with natural fractures and shear failure, they are missing some important features required for simulating a hydraulic fracturing treatment, such as proppant transport, non-Newtonian fluid behavior, temperature modeling, interaction with perforations, fracture conductivity, etc. The computation time is typically far too long for them to be suitable for fracture design and analysis purposes. Therefore, their use is limited to that of a research tool for obtaining insight to the fracturing process in a complex reservoir and as benchmark for simpler models. For engineering design application, simpler P3D based complex fracture models have been developed. Xu et al. (2009) and Meyer and Bazan (2011) presented network models in which the induced hydraulic fracture network consists of two orthogonal sets of parallel uniformly spaced fractures, in the directions of maximum and minimum horizontal stresses, respectively. The assumed fracture network geometry can be a close representation of the actual fracture networks, if the natural fractures system consists of two

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

16

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

Fig. 12. Hydraulic fractures generated in a medium with three pre-existing joints (blue disks are microcracks) by Damjanac et al. (2013). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

orthogonal sets of fractures that crisscross each other as observed in some shale outcrops (Meyer and Bazan), or can be thought as an idealized representation of the induced hydraulic fracture network which may have different fracture orientations and uneven fracture distribution, much like the dual porosity concept of Warren and Root (1963). With the assumed fracture network pattern, the coupled elasticity and flow equations are solved, and overall mass balance in the fracture system is observed to yield the network dimensions as a function of injection time or fluid volume. The effect of interference among parallel fractures is also taken into account in the fracture stiffness. Proppant transport in the fracture network is also modeled. Fig. 13 shows a simulated fracture network by Meyer and Bazan (2011). Cell-based P3D complex fracture network model In the lumped elliptic P3D network model presented above, the induced hydraulic fracture network pattern is predefined. With exception of the case in which the pre-existing DFN consists of two orthogonal sets of fractures and the induced fractures simply open along the existing DFN, the assumed fracture geometry may

not resemble the actual fracture network created. The model does not consider interaction between the hydraulic fracture and natural fractures, and is hence unable to predict the influence of initial DFN on the final induced hydraulic fracture network. Weng et al. (2011) presented a P3D based complex fracture network model, referred to as unconventional fracture model (UFM). It simulates the propagation, deformation, and fluid flow in a complex network of fractures. The model solves the fully coupled problem of fluid flow in the fracture network and the elastic deformation of the fractures, and has similar assumptions and governing equations as found in conventional P3D fracture models. But instead of solving the problem for a single planar fracture, the UFM model solves the equations for the complex fracture network. Fracture height growth is modeled in the same manner as in a conventional pseudo-3D model. A three-layer proppant transport model, consisting of a proppant bank at the bottom, a slurry layer in the middle, and clean fluid at the top, is adopted for simulating proppant transport in the fracture network. Transport equations are solved for each component of the fluids and proppant pumped. A key difference between the UFM model and the conventional planar fracture model is being able to simulate the interaction of

Fig. 13. DFN 2D and 3D aerial views (x–y plane) by Meyer and Bazan (2011).

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

hydraulic fractures with preexisting natural fractures, using the analytical OpenT crossing model discussed earlier in this paper. Additionally, UFM model also considers the interaction among hydraulic fracture branches by computing the stress shadow effect on each fracture exerted by the adjacent fractures. The basic governing equations to be solved include the equations governing fluid flow in the fracture network, mass conservation, fracture deformation, and the fracture propagation/ interaction criterion. The mass conservation equation, the equivalent of Eq. (3) in 1D along any branch of the fracture network, is given as

 @q @ðHfl wÞ þ þ qL ¼ 0; qL ¼ 2hL uL @s @t

ð15Þ

 is an average width or opening of where q is the local flow rate, w the fracture at position s = s(x, y), Hfl(s, t) is the local height of the fracture occupied by fluid, and qL is the leakoff volume rate through the wall of the hydraulic fracture into the rock matrix per unit length (leakoff height hL times leakoff velocity uL) which is expressed through Carter’s leakoff model, Eq. (4). The fracture tips propagate as a sharp front, and the total length of the entire hydraulic fracture networks at any given time t is defined as L(t). The fluid flow in the fracture could be laminar, turbulent, or Darcy flow through the proppant pack and is described correspondingly by different laws. For the case of laminar flow along any given fracture branch, Eq. (5) becomes

 n1 @p 1 q  q  ¼ a0 2nþ1  @s Hfl Hfl  w

ð16Þ

with

a0 ¼

2K 0 /ðn0 Þn

0



 0 n0 4n þ 2 ; n0

1 Hfl

/ðn0 Þ ¼

 2n0 þ1 wðzÞ n0 dz  w Hfl

Z

ð17Þ

Here w(z) represents fracture width as a function of depth z at the current position s(x, y). In addition to the equations described above, the global volume balance condition must be satisfied:

Z 0

t

QðtÞdt ¼

Z 0

LðtÞ

 tÞds þ hðs; tÞwðs;

Z HL

Z 0

t

Z

LðtÞ

2uL dsdtdhL

ð18Þ

0

That is, the total volume of fluid pumped is equal to volume of fluid in fracture network and volume leaked from the fracture up to time t. The boundary conditions require the flow rate, net pressure, and fracture width to be zero at all fracture tips. The total fracture network system consists of not only fractures but also the perforations and wellbore. The fracture networks communicate through injection elements to account for perforation friction, and perforation clusters are connected through wellbore elements to account for the friction in the casing. The system of equations, Eqs. (15)–(18), together with the elasticity equations Eqs. (8) and (9), and the initial and boundary conditions, plus the equations governing fluid flow in the wellbore and through the perforations, represent a complete set of governing equations. Combining these equations and discretizing the fracture network into small elements leads to a nonlinear system of equations in terms of fluid pressure p in each element, simplified as f(p) = 0. At each time step, each propagating fracture tip extends an incremental distance, according to the propagation criterion based on the local fluid velocity and fracture tip stress intensity factor, in the direction of local maximum horizontal principal stress (accounting for the stress shadow effect). The intersection of fracture tip with any natural fracture is checked, and if intersection with natural fracture occurs, the crossing model is applied to determine whether the tip crosses the natural fracture

17

or is arrested by it, and adjustment in fracture grid is applied accordingly. The system equation f(p) = 0 is then solved by using the damped Newton–Raphson method to obtain the new pressure and flow distribution in the fracture networks. Proppant transport equations are solved to update the proppant movement and settling in the fractures, and fracture height and stress shadow are also updated. A more detailed description of the model can be found in Weng et al. (2011). Fracture network growth pattern is affected by the mechanical interaction among the adjacent fractures. Generally known as the ‘‘stress shadow’’ effect, the stress field of one fracture is perturbed by the opening and shearing displacements of other nearby fractures. In a 2D plane-strain displacement discontinuity solution, Crouch and Starfield (1983) described the normal and shear stresses (rn and rs) acting on one fracture element induced by the opening and shearing displacement discontinuities (Dn and Ds) from all fracture elements as following:

rin ¼

N X j¼1

r

i s

Aij C ijns Djs þ

N X Aij C ijnn Djn j¼1

N N X X ¼ Aij C ijss Djs þ Aij C ijsn Djn j¼1

ð19Þ

j¼1

where Cij are the 2D, plane-strain elastic influence coefficients. This method, referred to as the 2D displacement discontinuity method (2D DDM) is used in the model to compute the additional stresses induced on each fracture element from the displacements of adjacent elements. In addition, a 3D correction factor Aij suggested by Olson (2008) (referred to as enhanced 2D DDM) was further introduced to modify the influence of coefficients Cij in the above equation to account for the 3D effect due to finite fracture height that leads to diminishing interaction between any two fracture elements when the distance increases (see Warpinski and Branagan, 1989). This additional normal stress due to the stress shadow is computed at each time step in the UFM model and is then added to the initial in-situ stress field on each fracture element during pressure and width iteration. The effect of stress shadow, including the shear stress, on the directional change of propagating fracture tips is essential for proper modeling of fracture network propagation pattern. The UFM model that incorporates the enhanced 2D DDM approach is validated against full 2D DDM simulator incorporating a full solution of coupled elasticity and fluid flow equations by CSIRO (Zhang et al., 2007) in the limiting case of very large fracture height (because the 2D DDM approach does not consider the 3D effect of finite fracture height). The comparison of influence of two closely propagating fractures on each other’s propagation paths has been provided in Kresse et al. (2012) and is shown in Fig. 14. The effect of stress shadow on fracture geometry is highly influenced by many parameters. Fig. 15, for example, shows the fracture geometry predicted by the UFM model for five fractures propagating from a horizontal well with perforation cluster spacing of 10, 20, and 40 m, respectively. When fracture spacing is large (larger than fracture height), the effect of stress shadow dissipates, and fractures have approximately the same dimensions; as the distance between fractures is reduced, the effect of stress shadow becomes greater and is noticeably observed in the reduced width and length for the inner fractures. These observations are consistent with those from 3D simulations such as those by Castonguay et al. (2013). The UFM model also considers the stress shadow contribution from the still open or propped fractures generated in the previous pumping stages on the current treatment stage. In the case of complex fractures, the existing fractures not only affect the current treatment stage through stress shadow effect, they also present as open fracture planes that can terminate and communicate with

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

18

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

Fig. 14. Comparison of propagation paths for two initially parallel fractures, with and without offset, in isotropic and anisotropic stress fields, by Kresse et al. (2012).

Fig. 15. Fracture geometry and fluid pressure for five fractures propagating from a horizontal well when distance between perforation clusters is equal to 10, 20, and 40 m; by Kresse et al. (2012).

the newly created fractures in the current stage. Examples of stress shadow effect on complex fractures can be found in Kresse et al. (2012). The UFM model has been used to simulate complex fractures in unconventional reservoirs. Some field examples can be found in Cipolla et al. (2011b), Kennaganti et al. (2013), and Liu et al. (2013). Kresse et al. (2013) presented UFM simulation of a field case in Barnett shale by Warpinski et al. (2005), as shown in Fig. 16. The well was first treated with a crosslinked gel, and then refractured several months later with slickwater. The wellbore is approximately aligned with the direction of maximum horizontal stress (which favors longitudinal fractures, untypical of most current transverse stimulation orientations). The microseismic events of the two treatments, also shown in Fig. 16, indicated a much narrower induced network for the crosslinked gel treatment than for the slickwater treatment. Using a DFN model based on the

dominant set of fractures in the north 70° west direction by Gale et al. (2007), the simulated hydraulic fracture network in the two treatments closely agree with the microseismic footprint, as shown in Fig. 16. The example presented in Fig. 16 showing the difference in induced fracture networks from two treatments with different types of fluid matches the microseismic cloud trend observed in Warpinski et al. (2005) and shows that the UFM model with the OpenT crossing model is able to correctly replicate the observed hydraulic fracture patterns. The UFM model has been integrated with a reservoir simulator using an automatically generated unstructured grid that has very fine grid next to the fractures and a logarithmically increasing grid size away from the fracture (Cipolla et al., 2011a) to properly simulate reservoir flow in an ultralow-permeability formation. The integrated fracture-reservoir model can be used to evaluate

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

19

Fig. 16. Hydraulic fracture network simulated by the UFM model for a Barnett case by Kresse et al. (2013); the thin blue lines indicate the traces of the natural fractures on a horizontal plane. The insets of the microseismic map are from Warpinski et al. (2005). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

production performance from hydraulic fracturing stimulation of naturally fractured formation and optimize completion parameters as presented in Cohen et al. (2014). Summary This paper presents a general overview of the state-of-art of hydraulic fracturing models applied to simulate complex fractures created when stimulating naturally fractured shale reservoirs. Even though many modeling approaches have been explored and

models are developed to simulate complex fractures in the naturally fractured reservoir, most have some limitations, have limited focus, or lack full functionalities to simulate the entire fracturing process. Therefore, it is important for the engineers to understand the limitations and applicability of the models and apply them judicially. The advancement in complex fracture modeling progresses rapidly, driven by the need for predictive fracturing design tools for the rapidly growing development in unconventional reservoirs. Vastly improved models are expected to come in the nottoo-distant future.

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

20

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

Even though there has been significant progress made in understanding the mechanisms influencing the interaction between hydraulic fracture and natural fracture, precise prediction of detailed fracture geometry is still very challenging. This is especially true for small-scale natural fractures and weak planes that are difficult to characterize to be included in the fracture network modeling, and yet they could cause offsets/kinks in fracture path that affect pressure loss in the fracture and proppant transport. The impact of natural fractures on the stimulation outcome also presents new challenges to the operators and engineers who design completion and evaluate fracturing treatments. Due to the high uncertainty and stochastic nature of the natural fractures as input for fracture simulation, there is inherent large uncertainty in the model predictions. Obtaining good understanding of the reservoir and properly characterizing the natural fracture system provide a good baseline for model simulation and an important first step towards reducing uncertainty. Additional measurements such as microseismic monitoring and treating pressure data will greatly help calibrate the model and uncertain parameters to further reduce the uncertainty. Even so, the uncertainties in the predicted results need to be considered in well performance forecasts and engineering decision making. Acknowledgements The author would like to thank Schlumberger for permission to publish this paper. He would like also to thank his current and former colleagues in Schlumberger—Olga Kresse, Charles Cohen, Hongren Gu, Ruiting Wu, Wenyue Xu, Dimitry Chuprakov, Ed Siebrits, Mark Mack, Craig Cipolla, and Jose Adachi—whose work had directly or indirectly contributed to this paper. Appreciation also goes to many colleagues in the industry who have done ingenious and highly inspiring work in this fascinating field of hydraulic fracture modeling, of which only limited number are referenced here. The author owes a sincere apology for any important work overlooked. References Abbas, S., Gordeliy, E., Peirce, A., Lecampion, B., Chuprakov, D.A., Prioul, R., 2014. Limited height growth and reduced opening of hydraulic fractures due to fracture offsets: an XFEM application. In: Paper SPE 168622, SPE Hydraulic Fracturing Technology Conference, The Woodlands, Texas, USA, 4–6 February. Adachi, J.I., Detournay, E., 2002. Self-similar solution of a plane-strain fracture driven by a power-law fluid. Int. J. Numer. Anal. Meth. Geomech. 26, 579–604. Adachi, J.I., Siebrits, E., Peirce, A., Desroches, J., 2007. Computer simulation of hydraulic fractures. Int. J. Rock Mech. Min. Sci. 44, 739–754. Bagheri, M., Settari, A., 2008. Modeling of geomechanics in naturally fractured reservoirs. SPE Res. Eval. Eng. 11 (1), 108–118, Paper SPE 93083-PA. Bai, T., Maerten, L., Gross, M.R., Aydin, A., 2002. Orthogonal cross joints: do they imply a regional stress rotation? J. Struct. Geol. 24, 77–88. Barree, R.D., 1983. A practical numerical simulator for three-dimensional fracture propagation in heterogeneous media. In: SPE 12273, SPE Symposium Reservoir Simulation, San Francisco, California, USA, 15–18 November. Beugelsdijk, L.J.L., de Pater, C.J., Sato, K., 2000. Experimental hydraulic fracture propagation in a multi-fractured medium. In: SPE 59419, SPE Asia Pacific Conference in Integrated Modeling for Asset Management, Yokohama, Japan, 25–26 April. Blanton, T.L., 1982. An experimental study of interaction between hydraulically induced and pre-existing fractures. In: SPE 10847, SPE/DOE Unconventional Gas Recovery Symposium, Pittsburgh, Pennsylvania, USA, 16–18 May. Blanton, T.L., 1986. Propagation of hydraulically and dynamically induced fractures in naturally fractured reservoirs. In: Paper SPE 15261, SPE Unconventional Gas Technology Symposium, Louisville, Kentucky, USA, 18–21 May. Bui, H.D., 1977. An integral equation method for solving the problem of a plane crack of arbitrary shape. J. Mech. Phys. Solids 25, 29–39. Carter, B.J., Desroches, J., Ingraffea, A.R., Wawrzynek, P.A., 2000. Simulating fully 3D hydraulic fracturing. In: Zaman, M., Gioda, G., Booker, J. (Eds.), Modeling in Geomechanics. Wiley, pp. 525–557. Castonguay, S.T., Mear, M.E., Dean, R.H., Schmidt, J.H., 2013. Prediction of the growth of multiple interacting hydraulic fractures in three dimensions. In: Paper SPE 166259, SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, USA, 30 September – 2 October.

Chen, Z., 2012. Finite element modeling of viscosity-dominated hydraulic fractures. J. Pet. Sci. Eng. 88–89, 136–144. Chen, Z., 2013. An ABAQUS implementation of the XFEM for hydraulic fracture problems. In: International Conference for Effective and Sustainable Hydraulic Fracturing, Brisbane, Australia, 20–22 May. Chuprakov, D.S., Akulich, A.V., Siebrits, E., Thiercelin, M., 2010. Hydraulic fracture propagation in a naturally fractured reservoir. In: SPE 128715, SPE Oil and Gas India Conference and Exhibition, Mumbai, India, 20–22 January. Chuprakov, D., Melchaeva, O., Prioul, R., 2013a. Injection-sensitive mechanics of hydraulic fracture interaction with discontinuities. In: ARMA 47th US Rocks Mechanics/Geomechanics Symposium, San Francisco, California, USA, 23–26 June. Cipolla, C.L., Williams, M.J., Weng, X., Mack, M., Maxwell, S., 2010, Hydraulic fracture monitoring to reservoir simulation: maximizing value. In: Paper SPE 133877, SPE Annual Technical Conference and Exhibition, Florence, Italy, 19–22 September. Cipolla, C.L., Fitzpatrick, T., Williams, M.J., Ganguly, U.K., 2011a. Seismic-tosimulation for unconventional reservoir development. In: Paper SPE 146876, SPE Reservoir Characterization and Simulation Conference and Exhibition, Abu Dhabi, UAE, 9–11 October. Cipolla, C.L., Weng, X., Mack, M., Ganguly, U., Gu, H., Kresse, O., Cohen, C., 2011b. Integrating microseismic mapping and complex fracture modeling to characterize fracture complexity. In: Paper SPE 140185, SPE Hydraulic Fracturing Technology Conference and Exhibition, The Woodlands, Texas, USA, 24–26 January. Cleary, M.P., 1980. Comprehensive design formulae for hydraulic fracturing. In: Paper SPE 9259, SPE Annual Technical Conference and Exhibition, Dallas, Texas, USA, 21–24 September. Cleary, M.P., Keck, R.G., Mear, M.E., 1983. Microcomputer models for the design of hydraulic fractures. In: SPE 11628, SPE/DOE Symposium on Low Permeability Gas Reservoirs in Denver, Colorado, USA, 14–16 March. Clifton, R.J., Abou-Sayed, A.S., 1981. A variational approach to the prediction of the three-dimensional geometry of hydraulic fractures. In: Paper SPE 9879, SPE/DOE Low Permeability Gas Reservoirs Symposium, Denver, Colorado, 27–29 May. Clifton, R.J., Wang, J.J., 1991. Adaptive Optimal Mesh Generator for Hydraulic Fracturing. In: J.C. Roegiers (Ed.), Rock Mechanics as a Multidisciplinary Science, Rotterdam: Balkema, pp. 607–616. Cohen, C.E., Kamat, S., Itibrout, T., Onda, H., Weng, X., Kresse, O., 2014. Parametric study on completion design in shale reservoirs based on fracturing-toproduction simulations. In: Paper IPTC 17462, International Petroleum Technology Conference, Doha, Qatar, 20–22 January. Cooke, M.L., Underwood, C.A., 2001. Fracture termination and step-over at bedding interfaces due to frictional slip and interface opening. J. Struct. Geol. 23, 223– 238. Crouch, S.L., Starfield, A.M., 1983. Boundary Element Methods in Solid Mechanics, first ed. George Allen & Unwin Ltd, London. Cundall, P.A., 1971. A computer model for simulating progressive large scale movements in blocky rock systems. In: Proceedings of the Symposium of the International Society for Rock Mechanics, 1, paper II-8. Cundall, P.A., 1988. Formulation of a three-dimensional distinct element model— Part I: a scheme to detect and represent contacts in a system composed of many polyhedral blocks. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 25, 107–116. Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Geotechnique 29, 47–65. Dahi-Taleghani, A., Olson, J.E., 2011. Numerical modeling of multistrand hydraulicfracture propagation: accounting for the interaction between induced and natural fractures. SPE J. 16 (3), 575–581, SPE 124884-PA. Dahi-Taleghani, A., Olson, J.E., 2013. How natural fractures could affect hydraulic fracture geometry. In: SPE J (preprint). SPE 167608-PA. Damjanac, B., Detournay, C., Cundall, P.A., Varun, 2013. Three-dimensional numerical model of hydraulic fracturing in fractured rock masses. In: International Conference for Effective and Sustainable Hydraulic Fracturing, Brisbane, Australia, 20–22 May. Dean, R.H., Schmidt, J.H., 2008. Hydraulic fracture predictions with a fully coupled geomechanical reservoir simulator. In: Paper SPE 116470, SPE Annual Technical Conference and Exhibition, Denver, Colorado, USA, 21–24 September. Dershowitz, W.S., Cottrell, M.G., Lim, D.H., Doe, T.W., 2010. A discrete fracture network approach for evaluation of hydraulic fracture stimulation of naturally fractured reservoirs. In: Paper ARMA10-475, 44th US Rock Mechanics Symposium, Salt Lake City, Utah, USA, 27–30 June. Detournay, E., 2004. Propagation regimes of fluid-driven fractures in impermeable rocks. Int. J. Geomech. 4, 1–11. Detournay, E., Garagash, D., 2003. The tip region of a fluid-driven fracture in a permeable elastic solid. J. Fluid Mech. 494, 1–32. Detournay, E., Adachi, J.I., Garagash, D., 2002. Asymptotic and intermediate asymptotic behavior near the tip of a fluid-driven fracture propagating in a permeable elastic medium. In: Dyskin, A., Hu, X., Sahouryeh, E. (Eds.), Structural Integrity and Fracture. Swets & Zeitlinger, Lisse, pp. 9–18. Fisher, M.K., Davidson, B.M., Goodwin, A.K., Fielder, E.O., Buckler, W.S., Steinberger, N.P., 2002. Integrating Fracture Mapping Technologies to Optimize Stimulations in the Barnett shale. In: Paper SPE 77411, SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, 29 September – 2 October. Fu, P., Johnson, S.M., Carrigan, C.R., 2011. Simulating complex fracture systems in geothermal reservoirs using an explicitly coupled hydro-mechanical model. In: Paper ARMA 11–244, 45th US Rock Mechanics Symposium, San Francisco, California, 26–29 June.

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx Fung, R.L., Viayakumar, S., Cormack, D.E., 1987. Calculation of vertical fracture containment in layered formations. SPE Form. Eval. 2 (4), 518–522. Gale, J.F.W., Holder, J., 2008. Natural fractures in the barnett shale: constraints on spatial organization and tensile strength with implications for hydraulic fracture treatment in shale-gas reservoirs. In: Paper ARMA 08–096, 42nd US Rock Mechanics Symposium and 2nd Canada Rock Mechanics Symposium, San Francisco, California, USA, 29 June and 2 July. Gale, J.F.W., Reed, R.M., Holder, J., 2007. Natural fractures in the Barnett shale and their importance for hydraulic fracture treatment. AAPG Bull. 91 (4), 603–622. Garagash, D.I., Detournay, E., 2002. Near tip processes of a fluid-driven fracture. ASME J. Appl. Mech. 67, 183–192. Geertsma, J., de Klerk, F., 1969. A rapid method of predicting width and extent of hydraulically induced fractures. J. Pet. Tech. 21 (12), 1571–1581, SPE 2458. Gordeliy, E., Peirce, A., 2013. Implicit level set schemes for modeling hydraulic fractures using the XFEM. J. Comput. Meth. Appl. Mech. Eng. 266, 125–143. Green, C.A., Barree, R.D., Miskimins, L., 2007. Development of a methodology for hydraulic fracturing models in tight, massively stacked, lenticular reservoirs. In: Paper SPE 106269, Hydraulic Fracturng Technology Conference, College Station, Texas, USA, 29–31 January. Gu, H., 1987. A Study of Propagation of Hydraulically Induced Fractures (Ph.D. dissertation), Department of Engineering Mechanics, The University of Texas at Austin, Austin, TX. Gu, H., Siebrits, E., 2008. Effect of formation modulus contrast on hydraulic fracture height containment. SPE Prod. Oper. (May), 173–176. Gu, H., Weng, X. 2010. Criterion for fractures crossing frictional interfaces at nonorthogonal angles. In: 44th US Rock Mechanics Symposium and 5th US–Canada Rock Mechanics Symposium, Salt Lake City, Utah, USA, 27–30 June. Gu, H., Yew, C.H., 1988. Finite element solution of a boundary integral equation for mode I embedded three-dimensional fractures. Int. J. Num. Meth. Eng. 26, 1525–1540. Gu, H., Weng, X., Lund, J.B., Mack, M., Ganguly, U., Suarez-Rivera R., 2011. Hydraulic fracture crossing natural fracture at non-orthogonal angles, a criterion, its validation and applications. In: Paper SPE 139984, Hydraulic Fracturing Conference and Exhibition, Woodlands, Texas, 24–26 January. Han, G., 2011. Natural fractures in unconventional reservoir rocks: identification, characterization, and its impact to engineering design. In: Paper ARMA 11–509, 45th US Rock Mechanics/Geomechanics Symposium, San Francisco, California, 26–29 June. Hanson, M.E., Anderson, G.D., Shaffer, R., Thorson, L.D., 1982. Some effects of stress, friction, and fluid flow on hydraulic fracturing. SPE J. (June), 321–332. Heuze, F.E., Shaffer, R.J., Ingraffea, A.R., Nilson, R.H., 1990. Propagation of fluid-driven fractures in jointed rock. part 1 – development and validation of methods of analysis. Int. Rock Mech. Min. Sci. Geomech. Abstr. 27 (4), 243–254. Hossain, M.M., Rahman, M.K., 2008. Numerical simulation of complex fracture growth during tight reservoir stimulation by hydraulic fracturing. J. Pet. Sci. Eng. 60, 86–104. Jeffrey, R.G., Weber, C.R., Vlahovic, W., Enever, J.R., 1994. Hydraulic fracturing experiments in the great northern coal seam. In: Paper SPE 28779, SPE Asia Pacific Oil & Gas Conference, Melbourne, Australia, 7–10 November. Jeffrey, R. G., Bunger, A., Lecampion, B., Zhang, X., Chen, Z., As, A., Allison, D.P., de Beer, W., Dudley, J.W., Siebrits, E., Thiercelin, M., Mainguy, M., (2009). Measuring hydraulic fracture growth in naturally fractured rock. In Paper SPE 124919, SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 4–7 October. Jeffrey, R.G., Zhang, X., Thiercelin, M., 2009. Hydraulic fracture offsetting in naturally fractured reservoirs: quantifying a long-recognized process. In: Paper SPE 119351, SPE Hydraulic Fracturing Technology Conference, Woodlands, Texas, 19–21 January. Ji, L., Settari, A., Sullivan, R.B., 2009. A novel hydraulic fracturing model fully coupled with geomechanics and reservoir simulation. SPE J. (September), 423–430. Jing, L., Ma, Y., Fang, Z., 2001. Modeling of fluid flow and solid deformation for fractured rocks with discontinuous deformation analysis (DDA) method. Int. J. Rock. Mech. Min. Sci. 38, 343–355. Karimi-Fard, M., Duflovsky, L.J., Aziz, K., 2004. An efficient discrete-fracture model applicable for general-purpose reservoir simulators. SPE J. 9 (2), 227–236. Kazemi, H., Atan, S., Al-Matrook, M., Dreier, J., Ozkan, E., 2005. Multilevel fracture network modeling of naturally fractured reservoirs. In: Paper SPE 93053, SPE Reservoir Simulation Symposium, The Woodlands, Texas, 31 January – 2 February. Kennaganti, K.T., Grant, D., Oussoltsev, D., Ball, N., Offenberger, R.M., 2013. Application of stress shadow effect in completion optimization using a reservoir-centric stimulation design tool. In: Paper SPE 164526, SPE Unconventional Resources Conference USA, The Woodlands, Texas, April 10–12. Koudina, N., Gonzalez Garcia, R., Thovert, J.-F., 1998. Permeability of threedimensional fracture networks. Phys. Rev. E 57, 4466–4479. Koutsabeloulis, N.C. Hope, S.A., 1998. Coupled stress/fluid/thermal multi-phase reservoir simulation studies incorporating rock mechanics. In: Proceedings of SPE/ISRM EUROCK-98 Symposium, Norway. Koutsabeloulis, N., Zhang, X., 2009. 3D reservoir geomechanical modeling in oil/gas field production. In: Paper SPE 126095, SPE Saudi Arabia Section Technical Symposium and Exhibition, Al Khobar, Saudi Arabia, 9–11 May. Kresse, O., Weng, X., Wu, R., Gu, H., 2012. Numerical modeling of Hydraulic fractures interaction in complex naturally fractured formations. In: ARMA-292, 46th US Rock Mechanics/Geomechanics Symposium, Chicago, Illinois, USA, 24– 27 June.

21

Kresse, O., Weng, X., Chuprakov, D., Prioul R., Cohen, C., 2013. Effect of flow rate and viscosity on complex fracture development in UFM model. In: International Conference for Effective and Sustainable Hydraulic Fracturing, Brisbane, Australia, 20–22 May. Lam, K.Y., Cleary, M.P., Barr, D.T., 1986. A complete three dimensional simulator for analysis and design of hydraulic fracturing. In: Paper SPE 15266, Unconventional Gas Technology Symposium, Louisville, Kentucky, 18–21 May. Liu, H., Luo, Y., Zhang, N., Yang, D., Dong, W., Qi, D., Gao, Y., 2013. Unlock shale oil reserves using advanced fracturing techniques: a case study in China. In: IPTC 16522, International Petroleum Technology Conference, Beijing, China, 26–28 March. Mack, M.G., Warpinski, N.R., 2000. Mechanics of hydraulic fracturing. In: S. Economides, K. Nolte (Eds.), Reservoir Stimulation, third ed., John Wiley & Sons (Chapter 6). Mack, M.G., Elbel, J.L., Piggott, A.R., 1992. Numerical representation of multilayer hydraulic fracturing. In: ARMA Paper 92–0335, 33rd U.S. Symposium on Rock Mechanics, Santa Fe, New Mexico, 3–5 June. Maxwell, S.C., Urbancic, T.I., Steinsberger, N.P., Zinno, R., 2002. microseismic imaging of hydraulic fracture complexity in the Barnett shale. In: Paper SPE 77440, SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 29 September – 2 October. McClure, M.W., 2012. Modeling and Characterization of Hydraulic Stimulation and Induced Seismicity in Geothermal and Shale Gas Reservoirs (Ph.D. dissertation), Stanford University. Meng, C., de Pater, C.J., 2010. Hydraulic fracture propagation in pre-fractured natural rocks. In: ARMA 10–318, 44th US Rock Mechanics Symposium and 5th US–Canada Rock Mechanics Symposium, Salt Lake City, Utah, 27–30 June. Meyer, B. R., 1989. Three-dimensional hydraulic fracturing simulation on personal computers: theory and comparison studies. In: SPE 19329, SPE Eastern Regional Meeting in Morgantown, West Virginia, October. Meyer, B.R., Bazan, L.W. 2011. A discrete fracture network model for hydraulically induced fractures: theory, parametric and case studies. In: Paper SPE 140514, SPE Hydraulic Fracturing Conference and Exhibition, Woodlands, Texas, USA, 24–25 January. Murphy, H.D., Fehler, M.C., 1986. Hydraulic fracturing of jointed formations. In: Paper SPE 14088, SPE 1986 International Meeting on Petroleum Engineering, Beijing, China, 17–20 March. Nagel, N.B., Gil, I., Sanchez-Nagel, M., Damjanac, B., 2011. Simulating hydraulic fracturing in real fractured rock – overcoming the limits of pseudo 3D models. In: Paper SPE 140480, SPE Hydraulic Fracturing Conference and Exhibition, Woodlands, Texas, 24–26 January. Nagel, N.B., Sanchez-Nagel, M., Lee, B., 2012. Gas shale hydraulic fracturing: a numerical evaluation of the effect of geomechanical parameters. In: Paper SPE 152192, SPE Hydraulic Fracturing Technology Conference and Exhibition, Woodlands, Texas, 6–8 February. Nakashima, T., Sato, K., Arihara, N., Yazawa, N., 2000. Effective permeability estimation for simulation of naturally fractured reservoirs. In: Paper SPE 64286, SPE Asia Pacific Oil and Gas Conference and Exhibition, Brisbane, Australia, 16– 18 October. Olson, J.E., 1993. Fracture pattern development: the effects of subcritical crack growth and mechanical interaction. J. Geophys. Res. 93, 12,251–12,265. Olson, J.E., 2004. Predicting fracture swarms – the influence of subcritical crack growth and crack-tip process zone on joint spacing in rock. In: Cosgrove, J.W., Engelder, T. (Eds.), The Initiation, Propagation, and Arrest of Joints and Other Fractures, vol. 231. Geological Society of London Special Publication, pp. 73– 87. Olson, J.E., 2008. Multi-fracture propagation modeling: applications to hydraulic fracturing in shales and tight sands. In: 42nd US Rock Mechanics Symposium and 2nd US–Canada Rock Mechanics Symposium, San Francisco, California, 29 June – 2 July. Olson, J.E., Dahi-Taleghani, 2009. Modeling simultaneous growth of multiple hydraulic fractures and their interaction with natural fractures. In: Paper SPE 119739, Hydraulic Fracturing Technology Conference, Woodlands, Texas, 19–21 January. Ouyang, S., 1994. Propagation of Hydraulically Induced Fractures with Proppant Transport (Ph.D. dissertation), The University of Texas at Austin, Austin, Texas. Peirce, A.P., Siebrits, E., 2001. The scaled flexibility matrix method for the efficient solution of boundary value problems in 2D and 3D layered elastic media. Comput. Methods Appl. Mech. Eng. 190 (45), 5935–5956. Pierce, M., MasIvars, D., Cundall, P.A., Potyondy, D.O., 2007. A synthetic rock mass model for jointed rock. In: Eberhardt, E. et al. (Eds.), Rock Mechanics: Meeting Society’s Challenges and Demands, 1st Canada–U.S. Rock Mechanics Symposium, Vancouver, May 2007, Fundamentals, New Technologies & New Ideas, vol. 1. Taylor & Francis Group, London, pp. 341–349. Potluri, N., Zhu, D., Hill, A.D., 2005. Effect of natural fractures on hydraulic fracture propagation. In: Paper SPE 94568, SPE European Formation Damage Conference, Scheveningen, Netherlands, 25–27 May. Renshaw, C.E., Pollard, D.D., 1995. An experimentally verified criterion for propagation across unbounded frictional interfaces in brittle, linear elastic-materials. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 32 (3), 237– 249. Riahi, A., Damjanac, B., 2013. Numerical study of interaction between hydraulic fracture and discrete fracture network. In: International Conference for Effective and Sustainable Hydraulic Fracturing, Brisbane, Australia, 20–22 May. Rijken, P., 2005. Modeling Naturally Fracturing Reservoirs: From Experimental Rock Mechanics to Flow Simulation (Ph.D. theses), University of Texas at Austin.

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001

22

X. Weng / Journal of Unconventional Oil and Gas Resources xxx (2014) xxx–xxx

Rungamornrat, J., Wheeler, M.F. Mear, M.E. 2005. A numerical technique for simulating nonplanar evolution of hydraulic fractures. In: Paper SPE 96968, SPE Annual Technical Conference and Exhibition, Dallas, Texas, 9–12 October. Rutledge, J.T., Phillips, W.S., Meyerhofer, M.J., 2004. Faulting induced by forced fluid injection and fluid flow forced by faulting: an interpretation of hydraulicfracture microseismicity, Carthage cotton valley gas field, Texas. Bull. Seismol. Soc. Am. 94, 1817–1830. Sanyal, S.K., Butler, S.J., Swenson, D., Hardeman, B., 2000. Review of the state-of-art of numerical simulation of enhanced geothermal systems. In: Proc. World Geothermal Congress, Japan, pp. 3853–3858. Sarda, S., Jeannin, L., Bourbiaux, B., 2001. Hydraulic characterization of fractured reservoirs: simulation of discrete fracture models. In: Paper SPE 66398, SPE Reservoir Simulation Symposium, Houston, Texas, 11–14 February. Savitski, A.A., Lin, M., Riahi, A., Damjanac, B., Nagel, N.B., 2013. Explicit modeling of hydraulic fracture propagation in fractured shales. In: IPTC 17073, International Petroleum Technology Conference, Beijing, China, 26–28 March. Sesetty, V. Ghassemi, A., 2012. Simulation of hydraulic fractures and their interactions with natural fractures. In: ARMA 12–331, 46th US Rock Mechanics/Geomechanics Symposium, Chicago, Illinois, 24–27 June. Settari, A. Cleary, M.P., 1982. Development and testing of a pseudo-threedimensional model of hydraulic fracture geometry (P3DH). In: Paper SPE 10505, Sixth SPE Symposium on Reservoir Simulation, New Orleans, Louisiana, 31 January – 3 February. Shi, G., 1988. Discontinuous Deformation Analysis: A New Numerical Model for the Static and Dynamics of Block System (Ph.D. theses), University of California, Berkeley. Shin, D.H., Sharma, M.M., 2014. Factors controlling the simultaneous propagation of multiple competing fractures in a horizontal well. In: Paper SPE 168599, SPE Hydraulic Fracturing Technology Conference, The Woodlands, Texas, 4–6 February. Siebrits, E., Peirce, A.P., 2002. An efficient multi-layer planar 3D fracture growth algorithm using a fixed mesh approach. Int. J. Numer. Meth. Eng. 53, 691–717. Siebrits, E., Gu, H., Desroches, J., 2001. An improved pseudo-3D hydraulic fracturing simulator for multiple layered materials. In: Proceedings of 10th International Conference on Computer Methods and Advances in Geomechanics, Tucson, Arizona, 7–12 January. Singh, G., Pencheva, G., Kumar, T., Wick, T., Ganis, B., Wheeler, M.F., 2014. Impact of accurate fractured reservoir flow modeling on recovery predictions. In: Paper SPE 168630, SPE Hydraulic Fracturing Technology Conference, The Woodlands, Texas, 4–6 February. Smith, M.B., Bale, A.B., Britt, L.K., Klein, H.H., Siebrits, E. Dang, X., 2001. Layered modulus effect on fracture propagation, proppant placement and fracture modeling. In: Paper SPE 71654, SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September – 3 October. Sousa, J.L.S., Carter, B.J., Ingraffea, A.R., 1993. Numerical simulation of 3D hydraulic fracture using Newtonian and power-law fluids. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 30, 1265–1271. Thiercelin, M., Makkhyu, E., 2007. Stress field in the vicinity of a natural fault activated by the propagation of an induced hydraulic fracture. In: Proceedings of the 1st Canada–US Rock Mechanics Symposium, vol. 2, pp. 1617–1624. Vandamme, L., Jeffrey, R.G., Curran, J.H., 1988. Pressure distribution in threedimensional hydraulic fractures. SPE Prod. Eng. 3, 181–186. Warpinski, N.R. Branagan, P.T., 1989. Altered-stress fracturing. In: SPE 17533, JPT, September, pp. 990–997. Warpinski, N.R., Teufel, L.W., 1987. Influence of geologic discontinuities on hydraulic fracture propagation (includes associated papers 17011 and 17074). SPE J. Pet. Technol. 39 (2), 209–220. Warpinski, N.R., Lorenz, J.C., Branagan, P.T., Myal, F.R., Gall, B.L., 1993. Examination of a cored hydraulic fracture in a deep gas well. J. SPE Prod. Facil. (August), 150– 164. Warpinski, N.R., Kramm, R.C., Heinze, J.R., Waltman, C.K., 2005. Comparison of single- and dual-array microseismic mapping techniques in the Barnett shale. In: Paper SPE 95568, SPE Annual Technical Conference and Exhibition, Dallas, Texas, 9–12 October. Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. SPE J. 3, 245–255. Weber, N., Siebert, P., Willbrand, K., Feinendegen, M., Clauser, C. Fries, T.P., 2013. The XFEM with an explicit–implicit crack description for hydraulic fracture

problems. In: International Conference for Effective and Sustainable Hydraulic Fracturing, Brisbane, Australia, 20–22 May. Weng, X., 1991. Incorporation of 2D fluid flow into a pseudo 3D hydraulic fracturing simulator. In: Paper SPE 21849, Rocky Mountain Regional Meeting and Low Permeability Gas Reservoirs Symposium, Denver, Colorado, April. Weng, X., Kresse, O., Cohen, C., Wu, R., Gu, H., 2011. Modeling of hydraulic fracture network propagation in a naturally fractured formation. In: Paper SPE 140253, SPE Hydraulic Fracturing Conference and Exhibition, Woodlands, Texas, USA, 24–26 January. Williams-Stroud, S.C., Barker, W.B., and Smith, K.L. 2012. Induced hydraulic fractures or reactivated natural fractures? modeling the response of natural fracture networks to stimulation treatments. In: ARMA 12–667, 46th US Rock Mechanics/Geomechanics Symp., Chicago, Illinois, 24–27 June. Willis-Richards, J., Watanabe, K., Takahashi, H., 1996. Progress toward a stochastic rock mechanics model of engineered geothermal systems. J. Geophys. Res. 1010 (B8), 17,481–17,496. Aug. 10. Wong, S.W., Geilikman, M., Xu, G., 2013. Interaction of multiple hydraulic fractures in horizontal wells. In: Paper SPE 163982, SPE Middle East Unconventional Gas Conference and Exhibition, Muscat, Oman, 28–30 January. Wu, K., Olson, and J.E., 2013. Investigation of Critical In situ and injection factors in multi-frac treatments: guidelines for controlling fracture complexity. In: Paper SPE 163821, SPE Hydraulic Fracturing Conference, Woodlands, Texas, 4–6 February. Xu, G., Wong, S.W. 2013. Interaction of multiple non-planar hydraulic fractures in horizontal wells. In: Paper SPE 17043, International Petroleum Technology Conference, Beijing, China, 26–28 March. Xu, W., Calvez, J.L., Thiercelin, M. 2009. Characterization of hydraulically-induced fracture network using treatment and microseismic data in a tight-gas formation: A geomechanical approach. In: Paper SPE 125237, SPE Tight Gas Completions Conference, San Antonio, Texas, USA, 15–17 June. Yamamoto, K., Shimamoto, T., Sukemura, S., 2004. Multiple fracture propagation model for a three-dimensional hydraulic fracturing simulator. Int. J. Geomechanics, 46–57. Yew, C.H., 1997. Mechanics of Hydraulic Fracturing. Gulf Publishing Company. Yew, C.H., Weng, X. 2014. Mechanics of Hydraulic Fracturing, second ed., Butterworth-Heinemann and Gulf Professional Publishing (to be published). Zhang, X., Jeffrey, R.G., 2006. The role of friction and secondary flaws on deflection and re-initiation of hydraulic fractures at orthogonal pre-existing fractures. Geophys. J. Int. 166 (3), 1454–1465. Zhang, X., Jeffrey, R.G., 2008. Reinitiation or termination of fluid-driven fractures at frictional bedding interfaces. J. Geophys. Res.-Sol. Ea. 113 (B8), B08416. Zhang, X., Jeffrey, R.G., Detournay, E., 2005. Propagation of a fluid driven fracture parallel to the free surface of an elastic half plane. Int. J. Numer. Anal. Meth. Geomech. 29, 1317–1340. Zhang, X., Jeffrey, R.G., Thiercelin, M., 2007. Effects of frictional geological discontinuities on hydraulic fracture propagation. In: SPE 106111, SPE Hydraulic Fracturing Technology Conference, College Station, Texas, 29–31 January. Zhang, X., Jeffrey, R.G., Thiercelin, M., 2007b. Deflection and propagation of fluiddriven fractures as frictional bedding interfaces: a numerical investigation. J. Struct. Geol. 29, 390–410. Zhang, X., Jeffrey, R.G., Thiercelin, M., 2009a. Mechanics of fluid-driven fracture growth in naturally fractured reservoirs with simple network geometries. J. Geophys. Res. 114, B12406. Zhang, Y., Sayers, C.M., Adachi, J., 2009b. The use of effective medium theories for seismic wave propagation and fluid flow in fractured reservoirs under applied stress. Geophys. J. Int. 177, 205–221. Zhang, F., Nagel, N., Lee, B., Sanchez-Nagel, M., 2013. The influence of fracture network connectivity on hydraulic fracture effectiveness and microseismicity generation. In: Paper ARMA 13–199, 47th American Rock Mechanics Symposium, San Francisco, California, USA, 23–26 June. Zhao, X.P. Young, R.P., 2009. Numerical simulation of seismicity induced by hydraulic fracturing in naturally fractured reservoirs. In: Paper SPE 124690, SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 4–7 October. Zhao, J., Chen, M., Jin, Y., Zhang, G., 2008. Analysis of fracture propagation behavior and fracture geometry using tri-axial fracturing system in naturally fractured reservoirs. Int. J. Rock Mech. Min. Sci 45, 1143–1152.

Please cite this article in press as: Weng, X. Modeling of complex hydraulic fractures in naturally fractured formation. J. Unconventional Oil Gas Resourc. (2014), http://dx.doi.org/10.1016/j.juogr.2014.07.001