Femi Oyedokun | Freelancer Portfolio Item #252981

Journal of Natural Gas Science and Engineering 44 -e364 Contents lists available at ScienceDirect Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse A quick and energy consistent analytical method for predicting hydraulic fracture propagation through heterogeneous layered media and formations with natural fractures: The use of an effective fracture toughness Oluwafemi Oyedokun*, Jerome Schubert Department of Petroleum Engineering, Texas A&M University, College Station, TX, USA a r t i c l e i n f o a b s t r a c t Article history: Received 30 August 2016 Received in revised form 3 May 2017 Accepted 4 May 2017 Available online 8 May 2017 An effective fracture toughness, based on equivalent energy release-rate hypothesis is presented for homogenizing heterogeneous layered media. Using crack closure method, the energy released when a mode-1 fracture propagates through an equivalent homogenized layer is equated to the sum of the energies released in the heterogeneous layered media. And from extensive numerical experiments, the predictions of fracture tips' positions through this proposed method are of the same range of accuracy as the known linear blend rule; the weakest link arguments technique performed poorly when compared to the other two methods. Therefore, homogenizing heterogeneous layered media with energy consistent approach will reduce the complexities associated with modeling fracture propagation in multi-layer without losing accuracy. Furthermore, an effective fracture toughness, based on equivalent energy release rate, equivalent strain energy, and modiﬁed Kachanov's damage theory is presented in this study. The proposed approach will reduce the computation time required in predicting fracture containment potential in formations with opened or sealed natural fractures. And comparing the proposed phenomenological model with the rigorous solution provided by Mori-Tanaka, it was observed that the margin of error was negligible. The beneﬁt of the proposed phenomenological model over the MoriTanaka's effective shear modulus model is the ease of estimating the effective fracture toughness. Thus, these models can be applied to quickly estimate the potential for hydraulic fracture containment or broaching possibilities during oil and gas blowouts and fracturing operations. Published by Elsevier B.V. Keywords: Effective fracture toughness Heterogeneous layered media Hydraulic fracture propagation Crack closure method Continuum damage mechanics Natural fractures 1. Introduction “Pseudo-3D” (P-3D) models are often used for quick hydraulic fracturing design or prediction of hydraulic fracture containment in the petroleum industry; these models are less computationally expensive compared to the fully 3D and planar 3D models (Morales et al., 1993; Clifton and Abou-Sayed, 1979, 1981; Hirth and Lothe, 1968; Bui, 1977). In the P-3D models, the ﬂuid ﬂow is assumed to be one-dimensional along the length of the fracture, while the pressure proﬁle in each elliptic cross sections is constrained to be linear. And as assumed in the PKN model (Perkins and Kern, 1961; Nordgren, 1972) that the cross sections are independent of each * Corresponding author. E-mail addresses:--(J. Schubert). http://dx.doi.org/10.1016/j.jngse-/Published by Elsevier B.V. (O. Oyedokun), other, so that the plane strain assumption decouples the ﬂuid ﬂow and solid mechanics; similarly, this constraint is applied in the P3D models. The linear pressure proﬁle in the cross sections is based on the assumption that the vertical fracture extension is slow such that the dynamic pressure gradient in the vertical sections is negligible. Many authors have contributed to the development of the equilibrium height problem (P-3D), starting from the work of Simonson et al. (1978), who developed the model for symmetric 3layer formation, with constant internal applied pressure. Ahmed (1984), and Newberry et al. (1985) solved an asymmetric 3-layer equilibrium height problem; while Fung et al. (1987) and Economides (2000) extended the asymmetric equilibrium height problem to multi-layer formations with constant and linearly varying pressure proﬁles in the vertical sections respectively. But as pointed out by Liu and Valko and Liu (2015), the model developed by Economides and Nolte has some intrinsic errors; and sometimes 352 O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 yield unacceptable results. Other notable contributions to the development of the P-3D models include Palmer and Carroll 1983a,b, Palmer and Craig (1984), Settari and Cleary (1986), Meyer (1986), Advani et al. (1990), and Adachi et al. (2010). Analytical P-3D model for asymmetric multi-layer formations with linearly varying pressure proﬁle in each cross section is very complicated, and difﬁcult to tract. To reduce the complexity, an effective fracture toughness, based on equivalent energy release rate, developed for the upper and lower layers will condense the multi-layer problem to the classical three-layer problem, which algebra is tractable. In the same vein, we also proposed in this paper the use of an effective fracture toughness, based on the same hypothesis mentioned earlier, which can approximately describe hydraulic fracture propagation in formations with natural fractures. The interactions between hydraulic and natural fractures have been studied by many authors in the past (Olson and Taleghani, 2009, Zhou et al., 2015; Renshaw and Pollard, 1995; Dahi-Taleghani and Olson, 2011, Gu et al., 2012; Gu and Weng, 2010; Chuprakov et al., 2014, Zhou et al., 2008; Warpinski and Teufel, 1987; Gao and Rice, 1989; Wu, 2014). From the various studies, there are three possibilities that may occur when a hydraulic fracture intersect a natural fracture: (1) the hydraulic fracture may cross the natural fracture without change in propagation direction (2) the hydraulic fracture is arrested by the natural fracture; consequently, the natural fracture, becoming part of the hydraulic fracture network and propagates from its end(s) or (3) from a weak point along its length. In this study, the complex hydraulic fracture network created by the interactions of the parent hydraulic fracture with ordered and disordered natural fractures will be simpliﬁed such that an effective fracture toughness, dependent on area of the discontinuities of a representative volume element and the mechanical properties of the formation matrix can suggest the containment potential of the hydraulic fracture in the host formation. In past studies, the use of linear blend rule (Atkins, 1979, Atkins and Mai, 1985; Eriksson and Atkins, 1995; Eriksson, 1998) and weakest link arguments (Landes and Shaffer, 1980; Wallin et al., 1984; Slatcher, 1986; Iwadate et al., 1983; Beremin et al., 1983) have been used to homogenize heterogeneous layered media. However, Heerens et al. (1994) showed that the predictions through the weakest link method are inaccurate. Fig. 1. Asymmetric multilayer hydraulic fracture propagation. having mode-1 fracture toughnessKc;1 ðzÞ. Therefore, the multilayer equilibrium-height problem reduces to the classical threelayer problem (Fig. 2). Assuming the formations have a linearly elastic response, the energy released per unit length of the fracture is P¼ ZHN ZH1 G1þ ðdzÞ þ 0 G1 dz (2) 0 Using superposition theorem, the normal stresses at the tips of the vertical fracture plane (in the y z plane) at arbitrary positions z0 and z1 are K I syy ðz ¼ z1 Þ ¼ syy ¼ Shmin ðz ¼ z1 Þ þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2pðzÞ 2. Mathematical formulations 2.1. Effective fracture toughness for layered media Fractures tend to grow in the direction that maximizes the potential energy released. And as a result, most hydraulic fractures propagate in mode-1; although diversion may occur at the boundary of two formations or when intercepted by natural fractures. In this study, the hydraulic fracture is assumed to propagate through boundaries in mode-1 only. As the hydraulic fracture propagates through each formation, there is an increase in the total energy released; and considering Fig. 1, the total energy released is PT ¼ ZHN ZL G1 dxdz H1 (1) 0 where G1 is the energy release rate, as deﬁned by Grifﬁth (1921). Layers n þ 1 to N can be lumped together as a heterogeneous formation having a varying fracture toughnessKc;3 ðzÞ, while layers 1 to n 1 are also lumped together as a heterogeneous formation 1, Fig. 2. Asymmetric three-layer equilibrium height problem. (3) O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 353 propagates into the abutting layer (lumped), is K Iþ ﬃ syy ðz ¼ z0 Þ ¼ syy þ ¼ Shmin ðz ¼ z0 Þ þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2pðzÞ (4) þ with the extension of the upper tip byz0, new fracture surfaces are created in 0 z z0 and the displacement of the upper-right face is uyy þ k þ 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2pðz0 zÞ ¼ K 4mp Iþ (5) Using crack closure method (Irwin, 1958; Broek, 1991), the work done by the normal stress in the fracture extension is equal to the energy required to close the fracture after opening Pþ ¼ 2 Zz0 1 þ syy uyy þ dz 2 (6) 0 P ¼ hN1 Z þh 0 syy ðzÞ k þ 1 þ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ K ðzÞ 2pðh zÞdz 4mp Ic (8) Hence, the effective fracture toughness for the lumped upper layer, þ þ KIC can be estimated from the hypothesis that P ¼ Pþ . Supposing the in-situ stresses are uniformly distributed throughout each layer, and the upper layers are lumped as a homogeneous þ formation with fracture toughness, KIC , therefore, the effective fracture toughness of the upper layer can be determined from the algebraic expression below. It should be noted that a Heaviside function was used to represent the uniformly distributed in-situ stresses in each layer in Eq. (8) above. And k and m are averaged properties (Poisson ratio function and shear modulus respectively) of the formations; for plane strain approximation, k ¼ 3 4n and plane stress k ¼ ð3 nÞ=ð1 þ nÞ. #) " N 3 3 2 X 2 2 ðh S þ h hn1 Þ Shmin;n ðhN1 þ h hn Þ 3 n¼1 hmin;n N1 " # pﬃﬃﬃﬃﬃﬃ N X 3 k þ 1 þ 2 2 2p kn þ 1 kn þ 1 2 2 Shmin;n KIC ðhN1 þ hÞ ¼ K ðhn hn1 Þ þ K ðhn hn1 Þ þ 8m 4mn p IC;n 8mn IC;n 3 n¼1 pﬃﬃﬃﬃﬃﬃ k þ 1 þ K 2p 4mp IC ( The factor 2 in Eq.6 accounts for the right and left surfaces of the fracture. And2uyy þ ðzÞ ¼ wþ ð0; zÞ. Therefore, the total energy released when the upper fracture tip propagates from the bottom of formation nþ1 to point h from the bottom of layer N is Alternatively, the shear modulus and Poisson ratio can also be represented with Heaviside functions, which will further make Eq. (9) more complicated. Similarly, the effective fracture toughness for the formations below the host formation can be estimated from Eq. 10 " #) M 3 3 2 X Shmin;n ðhM1 þ h hn1 Þ2 Shmin;n ðhM1 þ h hn Þ2 3 n¼1 " # pﬃﬃﬃﬃﬃﬃ M X 3 k þ 1 2 2 2p kn þ 1 k þ 1 n 2 Shmin;n KIC ðhM1 þ hÞ ¼ K ðhn hn1 Þ2 þ K ðhn hn1 Þ þ 8m 4mn p IC;n 8mn IC;n 3 n¼1 pﬃﬃﬃﬃﬃﬃ k þ 1 K 2p 4mp IC Pþ ¼ hnþ1 Z 0 ( hnþ2 Z syy;nþ2 hnþ1 hN1 Z þh …þ hN1 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ knþ2 þ 1 KIc;nþ2 2pðhnþ2 hnþ1 zÞ dz þ … 4mnþ2 p syy;N (10) 2.2. Effective fracture toughness derivation without tensile stresses qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ k þ1 syy;nþ1 nþ1 KIc;nþ1 2pðhnþ1 zÞ dz 4mnþ1 p þ (9) qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ kN þ 1 KIc;N 2pðh zÞdz 4mN p As fracture propagates through each medium, the tensile forces acting on the fracture edges, away from the tips reduces with crack length. Thus, neglecting the tensile force contributions and following the same procedure in Section 2.2 above, the effective fracture toughness for the homogenized upper and lower barriers, based on equivalent energy-release rate hypothesis are derived as þ (7) Where, the Irwin (1957) criterion, KI ¼ KIC was implied in Eq. (7). From Fig. 2, the equivalent energy released, as the upper tip K IC ¼ PN 3 kn þ1 2 n¼1 mn Shmin;n KIC;n ðhn hn1 Þ n 3 3 kþ1 PN S 2 2 n¼1 hmin;n ðhN1 þ h hn1 Þ ðhN1 þ h hn Þ m o (11) 354 K IC O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 3 PM kn þ1 2 n¼1 mn Shmin;n KIC;n ðhn hn1 Þ n o ¼ 3 3 P M kþ1 2 2 n¼1 Shmin;n ðhM1 þ h hn1 Þ ðhM1 þ h hn Þ m (12) The fracture extension pressure may not only be affected by gravity, but the relative positions of the fracture tips can affect the value. Considering the three different cases in Fig. 3. The critical extension pressures for the lower and upper tips in case 1 respectively, as derived in the Appendix are 0 1 1 BK C ﬃﬃﬃﬃﬃﬃﬃﬃ þ S2;hmin ghx rA PFD ¼ @qIC 2 phx 2.3. Minimum fracture extension pressures When the pressure driving the hydraulic fracture is assumed to be constant, the fracture extension pressure was derived by Simonson et al. (1978) as KIC PF ¼ S2;hmin ðzÞ þ qﬃﬃﬃﬃﬃﬃﬃﬃ phx (13) where, hx is the half thickness of the host formation. But when the displacements of the upper and lower tips of the fracture are signiﬁcant, the hydrostatic pressure can signiﬁcantly affect the extension potential. 0 PFU þ (14) 1 1 BK C ﬃﬃﬃﬃﬃﬃﬃﬃ þ S2;hmin þ ghx rA ¼ @qIC 2 phx (15) Noting that S2;hmin is the minimum principal in-situ stress at the host formation. Cases 2 and 3 are not very common; but could be relevant to the extension of existing fractures in formations. In case 2, the critical extension pressures for the upper and lower tips are respectively derived in the Appendix as sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! #) 1 þ k3 1 þ k3 2KIC þ k3 þ p S2;hmin 1 k3 1 k3 sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! qﬃﬃﬃﬃﬃ 1 þ k3 1 þ k3 þ k3 2 hx þ 1 k3 2 þ p 1 k3 1 k3 "rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ # qﬃﬃﬃﬃﬃ 1 k3 4 hx S3;hmin S2;hmin arcsin 2 sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ þ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! qﬃﬃﬃﬃﬃ 1 þ k3 1 þ k3 2 þ k3 þ 1 k3 þ p 2 hx 1 k3 1 k3 " ( qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃ p þ hx ghx pr þ 2 S3;hmin 1 k3 2 þ þ pﬃﬃﬃ PFU Fig. 3. Possible fracture tips positions in three-layer media. (16) O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 PFD sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ( " ! #) pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 k3 1 k3 1 k3 k3 þ 1 ghx pr þ 2 S3;hmin ð1 þ k3 Þ k þ þ p S2;hmin 1 þ k3 1 þ k3 3 1 þ k3 sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ ! ) ( pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 k3 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 k3 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 k3 1 þ k3 þ 1 k3 1 þ k3 p 1 þ k3 2 k3 1 þ k3 1 þ k3 sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ "rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ # pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 k3 pﬃﬃﬃ 1 þ k3 4 k3 þ 1 S3;hmin S2;hmin arcsin 2KIC p hx 2 sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ( ! ) pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 k3 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 k3 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 k3 1 þ k3 þ 1 k3 1 þ k3 p 1 þ k3 2 k3 1 þ k3 1 þ k3 h x Noting that, k3 ¼ h þh , where hu is the displacement of the u x upper tip. For case 3, the critical extension pressures for the upper and lower tips are 355 (17) Sometimes the natural fractures may be ﬁlled with materials that have higher modulus than the rock matrix; on the other hand the inﬁll materials may have lower Young's modulus than the sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! #) 1 k1 1 k1 2KIC þ k1 p S1;hmin 1 þ k1 1 þ k1 sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! qﬃﬃﬃﬃﬃ 1 k1 1 k1 2 k1 þ 1 k1 þ p 2 hx 1 þ k1 1 þ k1 "rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ # qﬃﬃﬃﬃﬃ 1 þ k1 4 hx S2;hmin S1;hmin arcsin 2 sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ þ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! qﬃﬃﬃﬃﬃ 1 k1 1 k1 2 k1 þ 1 k1 þ p 2 hx 1 þ k1 1 þ k1 (18) sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ rﬃﬃﬃﬃﬃ p 1 þ k1 1 þ k1 KIC ghx pr S2;hmin þ k1 S2;hmin þ S1;hmin 1 k1 2 þ pS1;hmin hx 1 k1 1 k1 sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! 1 þ k1 1 þ k1 2 þ k1 þ 1 k1 þ p 1 k1 1 k1 "rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ # "rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ # 1 k1 1 þ k1 ghx r arcsin þ ghx r þ 2S2;hmin 2S1;hmin arcsin 2 2 sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ þ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! qﬃﬃﬃﬃﬃ 1 þ k1 1 þ k1 þ k1 þ 1 k1 2 þ p 2 hx 1 k1 1 k1 (19) " ( qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃ p þ hx ghx pr þ 2 S2;hmin 1 k1 2 þ pﬃﬃﬃ PFU PFD And k1 ¼ h hx . d þhx 2.4. Effective fracture toughness for formations with disordered natural fractures The growth of hydraulic fractures in formations with natural fractures depends on many factors, including the orientation of the micro cracks, the inﬁll materials in the cracks, the tensile strength of the rock, treatment pressure, ﬂuid viscosity, and many other factors as studied by Zhou et al., 2015. To describe the random path of the hydraulic fracture in disordered natural fracture network will be very complicated. But the use of damage theory with the crack closure method (described above) can approximately suggest the containment potential of the hydraulic fracture. matrix. The surface density of the cracks and its mechanical strength affect the bearing capacity of the crack-matrix system. If the cracks are ﬁlled with higher modulus materials, they help to reinforce the strength of the crack-matrix system; but when the cracks are empty or ﬁlled with lower modulus materials, there is a reduction in the bearing capacity of the crack-matrix system. This is the crux of the effective stress introduced by Kachanov (2013); in the classical continuum damage theory, the inclusions are empty; thus, this study will cover both opened and ﬁlled discontinuities. For disordered natural fracture network, the impact of the inclinations of the natural fractures to the propagating hydraulic fracture cannot be included in the model development in a deterministic way, except by the use of random theory, which we have decided not to pursue in this paper. As mentioned in the previous section, the effective fracture 356 O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 Fig. 4. Region map of the tips positions when the pressure inside the fracture is (a) 53.09 MPa (7700 psi), (b) 53.78 MPa (7800 psi), and (c) 54.47MPa (7900 psi) for 11-EHP. toughness that describes the containment potential of the hydraulic fracture, propagating in the rock with natural fractures, is based on the use of equivalent energy released rate and equivalent strain energy (by compressional loading) hypotheses. In the equivalent system, the hydraulic fracture network is represented by a mode-1 fracture propagating in a homogeneous medium; and having the same energy released as the fracture propagates to the same height or length. In the equivalent system, the energy released due to mode-2 propagation of the hydraulic fracture is disregarded: when hydraulic fracture is arrested by a natural fracture and consequently the natural fracture propagates from its ends, if the natural fracture is inclined to the direction of maximum in-situ stress in the plane, the propagation will be mixed mode. The model development are based on the following assumptions: Assumption 1: The surface discontinuities do not evolve with the applied loadings. This suggest that there is no need to develop an evolution equation for the damage variable, which is the ratio of the effective surface area of the discontinuities in a representative volume element (RVE) and surface area of the RVE. Assumption 2: The surface discontinuities are isotopically distributed. Assumption 3: The rock behaves as a linearly elastic material. Assumption 4: The natural fractures are ﬁlled with a material that also has a linearly elastic response. Table 1 Description of the formation properties and in-situ stress proﬁles for the 11- eleven problem. Layer - Top Thickness Stress KIC E m. MPa pﬃﬃﬃﬃﬃ MPa m n m. MPa e e - - - - 29,510 29,510 29,510 29,510 29,510 29,510 29,510 29,510 29,510 29,510 29,510 - Shale Shale Sand Shale Sand Shale Sand Shale Sand Shale Shale Lithology O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 357 Fig. 5. Region map of the tips positions when the pressure inside the fracture is (a) 53.09 MPa (7700 psi), (b) 53.78 MPa (7800 psi), and (c) 54.47MPa (7900 psi) for 3ER-EHP. Table 2 Description of the formation properties and in-situ stress proﬁles for the equivalent three-layer problem. 358 O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 Fig. 6. Region map of the tips positions when the pressure inside the fracture is 54.47MPa (7900 psi) for 3BR-EHP. The effective stress tensor deﬁned by Kachanov for a material body having empty inclusions is given as s0 ¼ s (20) 1f But assuming the ﬁlling of the discontinuities reduces its effective area relative to the matrix, then, the effective area of the discontinuities is proposed as Sf;e ¼ Sf ð1 lÞ b Fig. 8. Region map of the tips positions when the pressure inside the fracture is 54.47MPa (7900 psi) for 11-EHP, Case 2. material has higher elastic modulus than the matrix, then the overall elastic modulus of the inclusion-matrix system will be higher; a reinforcement of the matrix. But the works by Eshelby (1957), Hashin and Shtrikman (1963), Hill (1965), Mori and Tanaka (1973), Willis (1981), Walpole (1981), Weng (1984), and Chen et al. (1992) provide more rigorous studies on the impact of inclusions on the elastic modulus of the host matrix. However, the computation of the effective fracture toughness is demanding (Li and Zhou, 2013). As stated in Assumption 1, the damage variable is deﬁned as (21) where l is the ratio of elastic moduli of the inﬁll material and the matrix, and b is a constant parameter. In a case where the inﬁll Fig. 7. Region map of the tips positions when the pressure inside the fracture is 54.47MPa (7900 psi) for 3WL-EHP. D¼ Sf;e ¼ fð1 lÞb S (22) Fig. 9. Region map of the tips positions when the pressure inside the fracture is 54.47MPa (7900 psi) for 3ER-EHP, Case 2. O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 359 Fig. 10. Region map of the tips positions when the pressure inside the fracture is 54.47MPa (7900 psi) for 3BR-EHP, Case 2. Fig. 12. Region map of the tips positions when the pressure inside the fracture is 54.47MPa (7900 psi) for 3EER-EHP, Case 2. and using the crack closure method, the energy released in the damaged rock material is equivalent to the energy released in a virgin rock with the normal stress acting on the face of the fracture being the effective stress, i.e. u0yy ¼ k0 þ 1 0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2pðz0 zÞ K 4m0 p IC (24) and, K0 IC s0yy ¼ S0hmin þ pﬃﬃﬃﬃﬃﬃﬃﬃ 2pz P¼2 Zz0 0 1 syy uyy dz ¼ 2 2 Zz0 1 0 0 s u dz 2 yy yy (23) uyy ¼ 0 Where, (25) k þ 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ K 2pðz0 zÞ 4mp IC (26) K IC syy ¼ Shmin þ pﬃﬃﬃﬃﬃﬃﬃﬃ 2pz (27) by substituting equations 22e25 into (21), the effective fracture toughness can be determined from the algebraic expression K 0IC p 3 k0 þ 1 0 2 z0 K IC S0hmin z0 2 þ pﬃﬃﬃﬃﬃﬃ 0 3 m 2p 2 3 kþ1 2 KIC p z0 ¼ KIC Shmin z0 2 þ pﬃﬃﬃﬃﬃﬃ 3 m 2p 2 = = (28) Table 3 Description of the formation properties and in-situ stress proﬁles for the 11-EHP, Case 2. Layer Fig. 11. Region map of the tips positions when the pressure inside the fracture is 54.47MPa (7900 psi) for 3WL-EHP, Case 2. - Top Thickness Stress KIC E m. MPa pﬃﬃﬃﬃﬃ MPa m n m. MPa e e - - - - 29,510 29,510 29,510 29,510 29,510 29,510 29,510 29,510 29,510 29,510 29,510 - Shale Sand Shale Shale Shale Shale Sand Shale Sand Shale Shale Lithology 360 O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 2.5. Effective fracture toughness for formations with ordered natural fractures When the natural fractures are oriented in the same direction, the path of the main hydraulic fracture stem can be easily traced compared to the disordered pattern. In this case the energy released as the fracture propagates in both modes 1 and 2 is P¼2 Z c 1 syy uyy ds þ 2 Z c2 1 2N sxy uxx ds2 ¼ 2 2 Zz0 1 0 0 s u dz 2 yy yy 0 (32) Fig. 13. Variations of the normalized shear modulusmD , and normalized fracture toughness with damage in any material displaying a linear elastic response. l ¼ 0 in this case. Since damage prior to loading is under consideration and not damage induced by the loading; thus, considering a compression loading of the RVE. The elastic strain energy (compressional loading) and the effective elastic energy are equivalent, i.e. 1 2 Z sεdv ¼ B 1 2 Z s0 ε0 dv0 (29) B0 Hence, the relationship between the damaged and virgin shear moduli is h i2 m ¼ m0 1 fð1 lÞb (30) The constant parameterb, can be determined by curve-ﬁtting the effective shear modulus in Eq. (30) with Mori-Tanaka's model or from experiment. When the contributions of the tensile stresses along the fracture edges are disregarded, the damaged fracture toughness of the formation becomes KIC ¼ K 0IC qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ b 1 fð1 lÞ (31) Where, c is the path of the main hydraulic fracture, and c2 is the distance the fracture propagates in mixed mode. N is the number of repeated pattern of the natural fractures; in other words, the number of rows. Eq. (32) is derived based on the assumption that the natural fracture will propagate from one of its ends, when it arrests the hydraulic fracture. uxx ¼ k þ 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2pðz0 zÞ K 4mp IIC K IIC sxy ¼ pﬃﬃﬃﬃﬃﬃﬃﬃ 2pz (33) (34) Depending on the projected path, the steps in Section 2.3 can be followed in estimating the effective fracture toughness from equation (32). 3. Numerical examples A hydraulic fracture stands in layer 6, with its tips at the top and bottom of the formation, which is 9 m (30 ft) thick. Using the full multi-layer equilibrium height model, the tips displacements at different internal pressures are shown in Fig. 4a b, and c. Table 1 shows the description of the layers. While the use of the energyreleased-based effective fracture toughness for the same range of pressures are shown in Fig. 5a, b, and c. The equivalent three-layer system, based on equivalent energy release rate hypothesis, (3ER-EHP) to the eleven-layer problem (11EHP) is shown in Table 2 below. Comparing the region maps in Fig. 4a,b and c and Fig. 5a and b, and 5c, it is evident that the two solutions are equivalent. Points A, B, and C in Fig. 4a, b, 5a, and 5b Fig. 14. (a) Variation of the normalized shear modulusmD , considering different inﬁll materials (b) Variation of normalized fracture toughness, KIC;D , with damage considering different inﬁll materials. O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 Table 4 Elastic properties of the matrix and inclusions. Case 1 2 3 Matrix Poisson Ratio Inﬁll Poisson Ratio Matrix Young's Mod. Inﬁll Young's Mod. GPa. GPa e e - 0 5.0 25.0 - - are not the practical solutions to the nonlinear equation; at these pressure values, the initial upper and lower tips are stable. In Figs. 4c and 5c Point B is the pair solution to the practical solution A. From the map, the upper tip is located at depth 2790 m, while the lower tip is located at depth 3005 m. Using the linear blend rule and weakest link arguments, the effective fracture toughness for the upper and lower barriers are 1.11 MPa and 1.1 MPa respectively. The locations of the tips in these equivalent systems are similar to the 3ER-EHP (Fig. 6 and 7); for shorthand writing, the equivalent three layer problem based on blend rule is 3BR-EHP, and that based on weakest link arguments is 3WL-EHP. In this case, both 3ER-EHP and 3WL-EHP yielded similar results, while 3BR-EHP and 3ERR-EHP gave more accurate positions of the tips relative to the other homogenization methods (Figs. 8e12). From these examples, it is evident that the use of energy-release 361 rate hypothesis for homogenizing multilayer problems may not be accurate (to a small degree), except when the contributions of the tensile stresses are neglected. Table 3 shows the description of eleven layers with signiﬁcant modulus contrasts. And the effective fracture toughness for the upper and lower barriers using the equivalent energy release rate hypothesis are 1.44MPa and 1.59 MPa respectively. While the effective fracture toughness for the upper and lower homogenized barriers when using the linear blend rule are 4.17 MPa and 5.46 MPa respectively; and using the weakest point arguments, the effective fracture toughness are 1.6MPa and 3.28MPa respectively. In the same vein, the effective fracture toughness using the reduced equivalent energy release rate hypothesis (3ERR-EHP) are 3.38MPa and 4.15 MPa respectively. Alternatively, the use of equivalent strain energy hypothesis, which is mathematically tedious, may be used in lieu of the energy released rate; this comparison will not be pursued in this paper. The tensile stress is greatest at the tips of the fracture, and reduces along the length of the fracture. Along the edges of the fracture, away from the tips, the Irwin criterion, KI ¼ KIC , does not apply; but KI < KIC . Thus, neglecting the tensile stress contribution in the crack closure energy function will yield more accurate results; and this is evident in the two example cases. Fig. 13 shows that the impact of damage on the elastic properties of the formation matrix do not follow the same path. The Young's and shear moduli are more severely damaged than the fracture toughness. And based on isotropic damage assumption, the Poisson Fig. 15. Comparing the performance of the proposed model for effective shear modulus of a material body having (a) empty or open natural fractures, Case 1, (b) micro-fractures ﬁlled with lower modulus material, Case 2, and (c) micro-fractures ﬁlled with higher modulus material, Case 3. 362 O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 Appendix. Derivation of minimum fracture propagation for three-layer problem ratio will not be affected by the presence of the discontinuities (damage), based on the continuum damage theory; the impact of damage on other elastic properties are beyond the scope of this article. And Fig. 14a and b, show the variations of the effective shear modulus and fracture toughness when the discontinuities (inclusions, or natural fractures) are ﬁlled with material having lower or higher Young's modulus. To compare the predictions of the proposed effective shear modulus for the damaged/reinforced material with Mori-Tanaka's model, Table 4 shows the three cases: As observed in Fig. 15a, the prediction from the proposed phenomenological model mimics the rigorous Mori-Tanaka model. In the second and third cases, the values of b are 2.5 and 1.2 respectively. Sneddon and Elliott (1946) derived the stress intensity factor at the tip of a fracture subjected to a constant internal pressure. Modifying the equation by including the effect of gravity, the stress intensity factor at the bottom tip, when the fracture is contained in formation n only: KI;D 1 ¼ pﬃﬃﬃﬃﬃ pl Let z* ¼ 4. Summary Using equivalent energy release-rate hypothesis, an effective fracture toughness was derived to homogenize heterogeneous layered media. Homogenizing heterogeneous layered media into a single layer can reduce multi-layer equilibrium-height problem to the classical three-layer equilibrium-height problem. And the reduction of the multi-layer problem to the classical three-layer problem reduces the model complexity. The predictions from the proposed model have the same range of accuracy as the wellknown linear blend rule; while the predictions from the weakest link arguments method were also observed to be inaccurate. Furthermore, an effective fracture toughness formulation for predicting the growth of hydraulic fracture in a formation with micro-cracks or any other form of ordered or disordered inclusions was developed. The formulation is based on the use of damage theory, equivalent energy-release rate, and equivalent strainenergy hypotheses. The presence of inﬁll materials in the microcracks (closed natural fractures) will affect the overall elastic properties of the damaged/reinforced rock-fracture system. And using a phenomenological approach, which is based on the use of an effective area, the proposed model performs very well as the rigorous Mori-Tanaka model, when the fractures are open. Then using the Mori-Tanaka model for the effective shear modulus, as reference, the constant parameter in the proposed model can be determined; subsequently, estimating the effective fracture toughness for the rock-closed fracture system. The beneﬁt of the proposed model over the Mori-Tanaka is its simplicity in computing the effective fracture toughness. And at Zz l rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ lþz pðz; tÞ sn dz lz z l (A.1) (A.2) z ¼ z; z* ¼ z* z ¼ l; z* ¼ 1 dz* 1 ¼ dz l (A.3) (A.4) Substituting (A.3) and (A.4) into (A.1), KI;D ¼ rﬃﬃﬃ Zz* rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 þ z* l dz* pQ ðtÞ þ rglz* sn 1 z* p (A.5) 1 Solving (A.5) and substituting (A.2) back yields n pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 KI;D ¼ pﬃﬃﬃpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ðl zÞ l þ zðgð2l þ zÞr 2sn Þ 2 p lðl zÞ 0qﬃﬃﬃﬃﬃﬃ1 lþz pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ l 1 @ pﬃﬃﬃ A þ 2l l z ðglr 2sn Þsin 2 0 0qﬃﬃﬃﬃﬃﬃ119 lþz = pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ l 1 þ pQ @2ðl þ zÞ l þ z þ 4l l z sin @ pﬃﬃﬃ AA ; 2 (A.6) Acknowledgement lim KI;D ¼ z/l This research is primarily supported by the Blowout Risk Assessment Joint Industry Projects. pﬃﬃﬃﬃﬃ lp glr þ 2pQ 2sn 2 (A.7) Therefore, the critical pressure at the mid position Q to cause the lower tip to move is Nomenclature D GIþ GI hu hd KIþ KI K IC K IC PFU þ damage variable energy release rate in upper barriers energy release rate in lower barriers displacement of the upper tip displacement of the lower tip stress intensity factor at upper tip stress intensity factor at lower tip effective fracture toughness for upper homogenized upper barriers mn mu md l effective fracture toughness for homogenized lower barriers f minimum fracture extension pressure for upper tip PFD Sh;min uyy þ uyy minimum fracture extension pressure for lower tip minimum horizontal in-situ stress fracture half width in the upper barriers fracture half width in the lower barriers formation, n, shear modulus Averaged shear modulus for upper barriers Averaged shear modulus for lower barriers ratio of inclusion and undamaged matrix elastic moduli Ratio of the surface area of discontinuity to surface area of representative volume element O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 pQ ¼ PFD ¼ KIC;D 1 pﬃﬃﬃﬃﬃ þ sn glr 2 pl (A.8) k3 approaches unity, the fracture extension pressure for the upper tip corresponds to case 3 in Fig. 3. Similarly, the stress intensity factor at the lower tip is At the upper tip, the stress intensity factor is KI;U 1 ¼ pﬃﬃﬃﬃﬃ pl Zz l sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ lz dz pðz; tÞ sn lþz KI;D (A.9) And at (A.10) (A.11) (A.12) Substituting (A.11) and (A.12) into (A.9), rﬃﬃﬃ Z1 rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 þ z* l dz* pQ ðtÞ rglz* sn 1 z* p (A.13) z* Solving (A.13) and substituting (A.2) back yields 1 KI;U ¼ pﬃﬃﬃpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 p lðl zÞ ðl zÞ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ðl þ zÞðgð2l þ zÞr þ 2sn Þ 0qﬃﬃﬃﬃﬃﬃ1 lþz qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ l 1 @ pﬃﬃﬃ A 2l ðl zÞ ðglr þ 2sn Þcos 2 0 0qﬃﬃﬃﬃﬃﬃ119 lþz qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ = l þ 2pQ @ðl zÞ ðl þ zÞ þ 2l ðl zÞ cos1 @ pﬃﬃﬃ AA ; 2 (A.14) lim KI;U ¼ z/l pﬃﬃﬃﬃﬃ lp glr þ 2pQ 2sn 2 (A.15) Therefore, the critical pressure at midpoint Q to cause the upper tip to move is KIC;U 1 pﬃﬃﬃﬃﬃ þ sn þ glr 2 pl (A.16) Noting that l ¼ hx . Hence, extending the derivation to three layer: the stress intensity factor at the upper tip is KI;U 8 rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ > Z1 1 þ z* 1 < dz* ¼ pﬃﬃﬃﬃﬃ pQ ðtÞ rglz* s3 1 z* pl > : k3 Zk3 þ k1 Zk1 þ z* k3 ZZ z ¼ z; z* ¼ z* z ¼ l; z* ¼ 1 pQ ¼ PFU ¼ 1 þ dz* 1 ¼ dz l KI;U ¼ 8 rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ > Zk3 1 þ z* 1 < dz* ¼ pﬃﬃﬃﬃﬃ pQ ðtÞ þ rglz* s3 1 z* pl > : Zk1 z Let z* ¼ l rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 þ z* dz* pQ ðtÞ rglz* s2 1 z* (A.17) 9 rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ > = 1 þ z* dz* pQ ðtÞ rglz* s1 > 1 z* ; z/l Therefore, the minimum fracture extension pressure for the upper tip at any location z is pQ in (A.17); since the algebraic expression is complicated, it will not be displayed in this paper. In the same vein, when k1 approaches unity, the fracture extension pressure for the upper tip corresponds to case 2 in Fig. 3. And when 363 þ k1 pQ ðtÞ þ rglz* s2 rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 þ z* dz* 1 z* (A.18) 9 rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ > = 1 þ z* dz* pQ ðtÞ þ rglz* s1 > 1 z* ; z/l and when k1 approaches unity, the fracture extension pressure for the lower tip corresponds to case 3 in Fig. 3. While when k3 approaches unity, the fracture extension pressure for the lower tip corresponds to case 2 in Fig. 3. References Adachi, J.I., Detournay, E., Peirce, A.P., 2010. Analysis of the classical pseudo-3D model for hydraulic fracture with equilibrium height growth across stress barriers. Int. J. Rock Mech. Min. Sci. 47 (4), 625e639. Advani, S.H., Lee, T.S., Lee, J.K., 1990. Three-dimensional modeling of hydraulic fractures in layered media: part Idﬁnite element formulations. J. Energy Resour. Technol. 112 (1), 1e9. Ahmed, U., 1984, January. A practical hydraulic fracturing model simulating necessary fracture geometry, ﬂuid ﬂow and leakoff, and proppant transport. In: SPE Unconventional Gas Recovery Symposium. Society of Petroleum Engineers. Atkins, A.G., 1979. Cracking of layered structures: the suppression of yielding by cast hardening. Mech. Behav. Mater. 3, 341e349. Atkins, A.G., Mai, Y.W., 1985. Elastic and Plastic Fracture: Metals, Polymers, Ceramics, Composites, Biological Materials. Halsted Press, Ellis Horwood. Beremin, F.M., Pineau, A., Mudry, F., Devaux, J.C., D'Escatha, Y., Ledermann, P., 1983. A local criterion for cleavage fracture of a nuclear pressure vessel steel. Metall. Trans. A 14 (11), 2277e2287. Broek, D., 1991. Elementary Engineering Fracture Mechanics, 4th Revised Edition. Kluwer Academic Publishers, Dordrecht, Netherlands. Bui, H.D., 1977. An integral equations method for solving the problem of a plane crack of arbitrary shape. J. Mech. Phys. Solids 25 (1), 29e39. Chen, T., Dvorak, G.J., Benveniste, Y., 1992. Mori-Tanaka estimates of the overall elastic moduli of certain composite materials. J. Appl. Mech. 59 (3), 539e546. Chuprakov, D., Melchaeva, O., Prioul, R., 2014. Injection-sensitive mechanics of hydraulic fracture interaction with discontinuities. Rock Mech. Rock Eng. 47 (5), 1625e1640. Clifton, R.J., Abou-Sayed, A.S., 1979, January. On the computation of the threedimensional geometry of hydraulic fractures. In: Symposium on Low Permeability Gas Reservoirs. Society of Petroleum Engineers. Clifton, R.J., Abou-Sayed, A.S., 1981, January. A variational approach to the prediction of the three-dimensional geometry of hydraulic fractures. In: SPE/DOE Low Permeability Gas Reservoirs Symposium. Society of Petroleum Engineers. Dahi-Taleghani, A., Olson, J.E., 2011. Numerical modeling of multistrandedhydraulic-fracture propagation: accounting for the interaction between induced and natural fractures. SPE J. 16 (03), 575e581. Economides, M.J. (Ed.), 2000. Reservoir Stimulation, vol. 18. Wiley, Chichester. Eriksson, K., 1998. The effective fracture toughness of structural components obtained with the blend rule. Nucl. Eng. Des. 182 (2), 123e129. Eriksson, K., Atkins, A.G., 1995. The effective through crack fracture toughness of structural components of older inhomogeneous steels. Mater. Ageing Compon. Life Ext. 1, 147e154. Eshelby, J.D., 1957, August. The determination of the elastic ﬁeld of an ellipsoidal inclusion, and related problems. In: Proceedings of the Royal Society of London a: Mathematical, Physical and Engineering Sciences, vol. 241. The Royal Society, pp. 376e396. No. 1226. Fung, R.L., Vilayakumar, S., Cormack, D.E., 1987. Calculation of vertical fracture containment in layered formations. SPE Form. Eval. 2 (04), 518e522. Gao, H., Rice, J.R., 1989. A ﬁrst-order perturbation analysis of crack trapping by arrays of obstacles. J. Appl. Mech. 56 (4), 828e836. Grifﬁth, A.A., 1921. The phenomena of rupture and ﬂow in solids. Philosophical Trans. R. Soc. Lond. Ser. A, Contain. Pap. a Math. or Phys. character 221, 163e198. Gu, H., Weng, X., 2010, January. Criterion for fractures crossing frictional interfaces at non-orthogonal angles. In: 44th US Rock Mechanics Symposium and 5th USCanada Rock Mechanics Symposium. American Rock Mechanics Association. Gu, H., Weng, X., Lund, J.B., Mack, M.G., Ganguly, U., Suarez-Rivera, R., 2012. Hydraulic fracture crossing natural fracture at nonorthogonal angles: a criterion 364 O. Oyedokun, J. Schubert / Journal of Natural Gas Science and Engineering 44 -e364 and its validation. SPE Prod. Operations 27 (01), 20e26. Hashin, Z., Shtrikman, S., 1963. A variational approach to the theory of the elastic behaviour of multiphase materials. J. Mech. Phys. Solids 11 (2), 127e140. Heerens, J., Zerbst, U., Bauschke, H.M., Schwalbe, K.H., 1994. On the Application of the Weakest Link Model in the Lower Shelf and Transition, 2013, February. In ECF10, Berlin. Hill, R., 1965. A self-consistent mechanics of composite materials. J. Mech. Phys. Solids 13 (4), 213e222. Hirth, J.P., Lothe, J., 1968. Theory of Dislocations, 780 pp. Irwin, G.R., 1957. Analysis of stresses and strains near the end of a crack traversing a plate. Spie Milest. Ser. MS 137 (167e170), 16. Irwin, G.R., 1958. Fracture in “Handbuch der Physik,”, vol. V. Iwadate, T., Tanaka, Y., Ono, S., Watanabe, J., 1983, January. An analysis of elasticplastic fracture toughness behavior for J IC measurement in the transition region. In: Elastic-plastic Fracture: Second Symposium, Volume II Fracture Resistance Curves and Engineering Applications. ASTM International. Kachanov, L., 2013. Introduction to Continuum Damage Mechanics, vol. 10. Springer Science & Business Media. Landes, J.D., Shaffer, D.H., 1980. Statistical characterization of fracture in the transition region. In: Fracture Mechanics. ASTM International. Li, Y., Zhou, M., 2013. Prediction of fracture toughness of ceramic composites as function of microstructure: I. Numerical simulations. J. Mech. Phys. Solids 61 (2), 472e488. Meyer, B.R., 1986, January. Design formulae for 2-D and 3-D vertical hydraulic fractures: model comparison and parametric studies. In: SPE Unconventional Gas Technology Symposium. Society of Petroleum Engineers. Morales, R.H., Brady, B.H., Ingraffea, A.R., 1993, January. Three-dimensional analysis and visualization of the wellbore and the fracturing process in inclined wells. In: Low Permeability Reservoirs Symposium. Society of Petroleum Engineers. Mori, T., Tanaka, K., 1973. Average stress in matrix and average elastic energy of materials with misﬁtting inclusions. Acta metall. 21 (5), 571e574. Newberry, B.M., Nelson, R.F., Ahmed, U., 1985, January. Prediction of vertical hydraulic fracture migration using compressional and shear wave slowness. In: SPE/DOE Low Permeability Gas Reservoirs Symposium. Society of Petroleum Engineers. Nordgren, R.P., 1972. Propagation of a vertical hydraulic fracture. Soc. Petroleum Eng. J. 12 (04), 306e314. Olson, J.E., Taleghani, A.D., 2009, January. Modeling simultaneous growth of multiple hydraulic fractures and their interaction with natural fractures. In: SPE Hydraulic Fracturing Technology Conference. Society of Petroleum Engineers. Palmer, I.D., Carroll Jr., H.B., 1983a. Three-dimensional hydraulic fracture propagation in the presence of stress variations. Soc. Petroleum Eng. J. 23 (06), 870e878. Palmer, I.D., Carroll Jr., H.B., 1983b, January. Numerical solution for height and elongated hydraulic fractures. In: SPE/DOE Low Permeability Gas Reservoirs Symposium. Society of Petroleum Engineers. Palmer, I.D., Craig, H.R., 1984, January. Modeling of asymmetric vertical growth in elongated hydraulic fractures and application to ﬁrst MWX stimulation. In: SPE Unconventional Gas Recovery Symposium. Society of Petroleum Engineers. Perkins, T.K., Kern, L.R., 1961. Widths of hydraulic fractures. J. Petroleum Technol. 13 (09), 937e949. Renshaw, C.E., Pollard, D.D., 1995, April. An experimentally veriﬁed criterion for propagation across unbounded frictional interfaces in brittle, linear elastic materials. Int J. Rock Mechanics. Min Sci. Geomechanics. Abstr 32 (No. 3), 237e249. Pergamon. Settari, A., Cleary, M.P., 1986. Development and testing of a pseudo-threedimensional model of hydraulic fracture geometry. SPE Prod. Eng. 1 (06), 449e466. Simonson, E.R., Abou-Sayed, A.S., Clifton, R.J., 1978. Containment of massive hydraulic fractures. Soc. Petroleum Eng. J. 18 (01), 27e32. Slatcher, S., 1986. A probabilistic model for lower-shelf fracture toughnessdtheory and application. Fatigue & Fract. Eng. Mater. Struct. 9 (4), 275e289. Sneddon, I.N., Elliott, H.A., 1946. The opening of a Grifﬁth crack under internal pressure. Q. Appl. Math. 4 (3), 262e267. Valko, P.P., Liu, S., 2015. February. An improved equilibrium-height model for predicting hydraulic fracture height migration in multi-layered formations. In: SPE Hydraulic Fracturing Technology Conference. Society of Petroleum Engineers. €rro €nen, K., 1984. Statistical model for carbide induced brittle Wallin, K., Saario, T., To fracture in steel. Metal Sci. 18 (1), 13e16. Walpole, L.J., 1981. Elastic behavior of composite materials: theoretical foundations. Adv. Appl. Mech. 21, 169e242. Warpinski, N.R., Teufel, L.W., 1987. Inﬂuence of geologic discontinuities on hydraulic fracture propagation (includes associated papers 17011 and 17074). J. Petroleum Technol. 39 (02), 209e220. Weng, G.J., 1984. Some elastic properties of reinforced solids, with special reference to isotropic ones containing spherical inclusions. Int. J. Eng. Sci. 22 (7), 845e856. Willis, J.R., 1981. Variational and related methods for the overall properties of composites. Adv. Appl. Mech. 21, 1e78. Wu, K., 2014. Mechanics Analysis of Interaction between Hydraulic and Natural Fractures in Shale Reservoirs. Unconventional Resources Technology Conference (URTEC). Zhou, J., Chen, M., Jin, Y., Zhang, G.Q., 2008. Analysis of fracture propagation behavior and fracture geometry using a tri-axial fracturing system in naturally fractured reservoirs. Int. J. Rock Mech. Min. Sci. 45 (7), 1143e1152 (SS). Zhou, J., Huang, H., Deo, M., 2015, November. Modeling the interaction between hydraulic and natural fractures using dual-lattice discrete element method. In: 49th US Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association.