Rock properties in simulation

## 1. Introduction

Hydraulic fracturing has been recognized as the most effective technique for economic recovery in tight oil and gas formations in North America [30], [36]. Hydraulically induced fractures increase well-reservoir contact area enormously; hence well productivity improves greatly after stimulation. During a typical hydraulic fracturing treatment, a mixture of proppant and viscous fluids is injected into the formation to create a fracture. The main mechanism responsible for fracturing the rock is the generation of tensile stresses ahead of pressurized fracture. The direction of the fracture will be perpendicular to the direction of minimum principal in-situ stress [23]. Well-testing analysis done at the early production life of these wells provide estimations for hydraulically induced fracture surface areas which are much larger than the fracture dimensions estimated in fracturing design or predicted areas constrained by the scattering domain of microseismic events. Presence of microcracks might be indicated by increased pore volume and compressibility, as well. This finding is speculated to be related to microcracking [42] i.e. a large population of microcracks could essentially explain this result. It is notable that microcracks are not necessarily micron size. We call them microcrack because they are much smaller than the major hydraulic fractures (millimetres in size). This hypothesis becomes more plausible by considering the fact that a large number of these tight sand and shale gas reservoirs [17] are naturally fractured. Presence of natural fractures and their fractal distribution is a widely observed fact in various tight sand and shale formations. The significance and role of these pre-existing natural fractures on the performance of fracturing treatments and post-frac production are not well-understood; consequently, most analysis is mainly descriptive rather than quantitative. In summary, there is no model to predict the likelihood of opening these fractures in different scales. For instance, few models have been introduced to predict interaction of hydraulic fractures with large size natural fractures [18], [31]. Here, large size natural fractures are fractures with the lengths and heights comparable with the size of hydraulic fractures. Laboratory experiments [10] have confirmed the influence of these large fractures in changing the direction of fracture propagation, and earlier shallow depth mineback experiments have shown similar outcomes [45]. However, there is no similar study about the role of microfractures. Almost all published models in the literature are limited to the cases in which natural fractures have the same height as that of hydraulic fractures. Considering the fact that power-law distribution of natural fractures implies population of small size fractures to be orders of magnitudes more than that of the large size fractures, it is not surprising that induced large fractures are intersecting thousands of these fractures. Due to their small sizes, small natural fractures cannot be propped by proppants; their aperture and therefore their hydraulic conductivity is a function of the inner-fracture fluid pressure. Due to their presence in large numbers, only partial reactivation of these natural microfractures may affect fluid flow pattern near the major fracture. These effects could be in the form of increasing the total effective wellbore-formation contact area and consequently improving hydrocarbon production, or oppositely, these microfractures could act as capillary traps for the fracturing fluid. The entrapped water, which is essentially part of the leakoff volume that will never produce, could hinder hydrocarbons flow from the formation into the major hydraulic fracture.

Low required energy for the re-opening of natural fractures makes them susceptible to re-opening if large enough tensile or shear stresses are somehow generated on the surface of major fractures. Then, depending on the distribution of natural fractures and the strength of their digenetic cements, their possible reactivation may influence hydrocarbon flow considerably. Despite the predominantly compressive stress regime around pressurized fractures under certain circumstances, it is possible to have tensile stresses. Two main mechanisms responsible for inducing tensile and/or shear forces on the surface of major fractures are thermal stresses and residual stresses due to the plastic deformation of the rock during hydraulic fracturing.

Figure 1 shows a typical response of the bottomhole pressure and temperature measurement during a fracturing treatment. Fluid and proppants have been pumped for a period of time, and the termination is marked by a red line and followed by an extended period of shut-in that lasts much longer than the pumping time [23]. Of particular interest here is that minimum temperature, minimum fracturing fluid pressure and maximum leakoff fluid pressure occurs almost simultaneously within a short period of time after shut-in. Minimum downhole temperature and maximum pore pressure due to leakoff could be essential factors in reducing rock effective stress. The red mark also indicates the onset of depressurization which also locally develops tensile stresses.

Fracturing fluid is frequently pumped with the temperature very close to the surface temperature; hence its temperature at the bottomhole usually differs from the initial temperature of the reservoir, especially in the case of deep and hot formations. The temperature gradient between the fracturing fluid and formation is a function of formation temperature, injection rate, casing/tubing diameter, the distance from perforations to the surface, heat capacity of fluids, fracture width and treatment pressure [8]. For most cases, fracturing fluid does not have enough time to reach the formation temperature due to its high velocity in the tubing. Because of the fluid migration and heat transfer in the reservoir, such differential temperature induces thermal stresses. The tensile and shear stresses induced by this temperature difference could be large enough to initiate small cracks on the fracture surface or in the case where pre-existing natural fracture are present, these stresses may open them. Thermal cracking happens when induced stresses inside the rock due to cooling exceed the in-situ stress of the formation, this phenomenon is well-documented in waterfloodings of brittle hot rocks [39] and geothermal systems with cold water circulation [46]. Thermal cracking may lead to the formation of clusters of small cracks, or so-called secondary fractures, which are very similar to pavement cracks but on the surface of the main hydraulic fracture.

As mentioned earlier, when the induced stresses inside the rock overcome formation in-situ stresses, re-opening of natural fractures may also occur. However, the spacing and geometry of opened cracks, in addition to previously mentioned parameters, are also functions of natural fractures distribution. Although these thermal induced cracks and re-opened parts of the pre-existing natural fractures have small size in comparison to the main hydraulic fracture, they can tremendously increase the well-formation contact area. For the case of no capillary trapping, the fluid flux inside these secondary fractures is roughly proportional to the cube of the fracture width and to the inverse of spacing length. Based on thermoelasticity analysis for closely spaced fractures, the fracture width is proportional to fracture spacing. Therefore, the fluid flux inside the thermal induced fracture is a quadratic function of spacing length [5]. Moreover, Bazant et al. [5] showed that the ratio of crack depth-to-spacing in pavements (elastic half-space) is a sensitive function of temperature profile inside the crack. Heat transfer for hydraulic fracturing has been studied in the last couple of decades. For instance, Biot et al. [8] proposed a one-dimensional analytical solution for heat transfer in the plane strain geometry. The fundamental solution for a centre of dilation and a point source fluid injection was provided earlier by Cleary [13]. Clifton and Wang [14] utilized this fundamental solution for a pseudo-three dimensional hydraulic fracturing simulations. However, these models are mainly investigating local changes of in-situ stresses rather than the likelihood of initiating secondary fractures. Study on the effect of stress redistribution on fractures due to thermal gradient of rock mass and fracturing fluid received more attention for geothermal reservoirs due to the presence of large temperature differences [3], [5], [46]. Zhou et al. [46] adapted this problem in the context of initiation of secondary fractures from a hydraulic fracture in hot dry geothermal systems with brittle rocks. Dahi-Taleghani et al. [19] considered the effect of induced thermal stress during hydraulic fracturing on opening of cemented natural fractures. They used the concept of cohesive interfaces in the framework of three dimensional finite element methods to show how thermal conductivity of the rock mass could make the population of opened natural fractures clustered rather than uniform. Additionally, their model considered interaction between propagating fractures.

Thermal stresses are not necessarily the only driving force behind formation of microfractures or opening of pre-exiting natural fractures. Plastic deformation induced during fracture pressurization results in tensile residual stress upon reduction of fracturing fluid pressure. Therefore, microcrack initiations can be enhanced upon unloading, as long as the pressurization at the pumping stage induces plastic deformation in rock. Cracking due to stress release resulted from unloading is a well-established mechanism in indentation experiments [16]. Choi et al. (2012) explains that plasticity is playing the main role in nucleation of microfractures during unloading. They showed the nucleation of microfractures from microscopic voids during unloading of hydraulic fracture. Plastic deformation induced during pressurization of main hydraulic fracture creates a tensile residual stresses during depressurization of hydraulic fracture. Therefore, these tensile residual stresses initiate the nucleation of microfractures; but compared to microcracks induced by thermal gradient, the effect of tensile residual stresses due to plasticity has not been studied so much. In this paper, the effect of plastic deformation on opening the natural fractures has been studied. In terms of methodology for modelling natural fracture reactivation, this paper is an extension of the work done by Dahi-Taleghani et al. [19] regarding the effect of plastic residual stresses.

As it mentioned before, thermal stresses and plasticity induced residual stresses may generate some microfractures or reactivate pre-existing natural fractures, but activation of these fractures does not necessarily lead to production enhancement due to the increase in contact area. If microfractures act as capillary traps, contact area and productivity index can be considerably influenced. Capillary trapping occurs when hydraulic pressure cannot overcome the capillary entrance pressure of microfracture to open it, and it’s a function of pore geometry, rock-fluid interaction and fluid flow inside the pores; therefore, considering capillary pressure effect and trapping mechanism is quite important to achieve a realistic prediction of fractured well productivity and the amount of producible leakoff fluid volume. Pore geometry and rock-fluid interaction control capillary trapping. Capillary trapping effect can become a quite interesting topic in hydraulic fractured reservoirs and naturally fractured reservoirs. To activate natural fractures, fracturing fluid pressure should go beyond the in-situ rock stresses; however, due to small aperture size of these fractures, if the hydraulic pressure cannot overcome the capillary entry threshold pressure of microfracture, formation fluid may not flow via the microfracture to reach the main fracture.

Due to the limited knowledge about the presence of natural fractures and their potential distribution in different formations, their contribution has been ignored or at least has not received enough attention. Only recent advances in characterization of natural fractures and verification of power-law distribution of fractures in different length scales [35], as well as the development of more sophisticated mechanistic models for fracture initiation and propagation such as cohesive crack models, made the investigation about the role of these natural fractures possible.

The remainder of the paper is organized as follows. In the next section, we talk about distribution of natural fractures, which is followed by sections about rock plasticity and a section regarding cohesive interface constitutive equations to model mechanical behaviour of pre-existing cemented natural fractures. At the end, numerical results of implementing this model for several examples will be presented to examine the significance of induced thermal stress in different situations.

## 2. Natural facture distribution

Fracture is a mechanical discontinuity in the rock mass formed due to the presence of stress fields in earth’s crust (Figure 2). There are wide scale ranges for fractures from micrometre (microfractures) to kilometres (lineaments). Presence of fractures in earth’s crust can influence underground fluid flow and physical properties of rock like rock strength. Fractures can influence the velocity of elastic waves and rock elastic moduli [41]. Natural fractures are categorized into four groups [40] based on their genesis : (1) tensile fractures due to compressive stresses, (2) shear fractures due to compressive stresses, (3) tensile fractures due to unloading of compressive stresses, (3) natural hydraulic fractures. Despite indeterministic nature of the aforementioned mechanisms, a large number of outcrop studies have revealed pattern and identifiable organization in fracture orientation and spacing. Due to the limited access to the subsurface to map fractures and limited precision of seismic techniques, outcrops are the main source to speculate fracture’s geometry in the subsurface. There are different distribution models used to describe fracture size like fracture length, aperture and tangential or perpendicular displacement due to fracture. Scale-limited laws (lognormal, exponential, gamma and power law) are methods in literature to characterize fracture systems [9], but it should be mentioned that scaling exponents alone cannot act as good criterion to define the whole pattern of fracture networks. Moreover, Bonnet et al. [9] showed that there is a linear relationship between rupture area and frequency scale of tensile fractures in seismometers acting. Field studies have confirmed the existence of a critical threshold that cracks with aperture less than this threshold are fully filled with digenetic materials [34]. Although microfractures are filled with calcite or quartz cements, laboratory measurements have proved that these filled natural fractures may still act as weak surfaces, or in other words, potential paths for rock failure. For instance, lab measurements for Barnett shale samples have shown tensile strength of cemented cracks to be about 10 times less than the tensile strength of intact rocks [27]. There exist some integrated models in the literature that can be utilized for this purpose [33]. By combining the knowledge of natural fracture patterns, cement properties and current in-situ stresses, it is possible to build a model to make a realistic prediction about the distribution of natural fractures in the case of limited core and outcrop data.

Proppants cannot move into microfractures opened during hydraulic fracturing due to their small aperture, which is less than a couple of microns. However, hydraulic pressure can open the microfractures if it goes beyond the local closure stress; therefore, activation of microfractures is function of confining pressure and pore pressure. As it mentioned earlier, contact area between rock-fluid can be considerably affected by the presence of microfractures in large quantities despite their small aperture and depth.

## 3. Elastoplastic effect in fracturing

The mechanical behaviour of quartz or calcite is essentially identified as elastic and brittle, however, clay/organic dominated regions can undergo significant plastic strains. Hence, it is not surprising that excessive fluid pressure present during hydraulic fracturing treatments may induce plastic deformations. This issue has been the subject of several studies in the literature [37], [38], [44]. For instance, it has been shown that plasticity causes shorter and wider fractures. However, most of these plastic deformations are due to high stress near the tip of the hydraulic fracture. The excess pressure in the main fracture may be only 1 or 2 MPa higher than the minimum in-situ stress, and this amount of additional stress may not cause a considerable plasticity unless in very weak formations

These papers were mainly focused on plastic deformations induced at the tip of fractures due to stress concentration at the tip of fractures, while plastic deformation of the surrounding rocks and its possible effects was out of the scope of these papers. Irreversible strain characterizes the plasticity when stress reaches a certain point. After this yielding point, the material shows elastoplasticity, which means its behaviour is somewhat plastic and also elastic. Equations (1) to (3) show general elastoplastic behaviour in three-dimensional problems for a strain increment

For flow rule and yield criterion, we used Drucker-Prager criterion, which is a well-known model to describe plastic behaviour of rocks. The Drucker-Prager criterion is an adjusted version of the Von Mises criterion for granular materials like soils and soft rocks. The yield surface for the Drucker-Prager criterion is a circular cone with the form shown in equation (4) where α and k are constants related to internal friction and cohesion of material. The stress at any point can be represented by the vector (σ_{1}, σ_{2}, σ_{3}). This vector can be shown by a corresponding stress point on the π-plane which is constituted of vector *s* (σ_{1}- σ_{m}, σ_{2}- σ_{m}, σ_{3}- σ_{m}) and ρ (σ_{m}, σ_{m}, σ_{m}), where σ_{m} equals to (σ_{1}+σ_{2}+σ_{3})/3 (Figure 3a). The constants can be obtained from the plot of failure in

## 4. Initiation and propagation of cracks

There is a considerable number of publications for modelling hydraulic fracturing treatments published since 1955, these solutions are varying from analytical and asymptotic solutions [21], [22], [32] to finite element or boundary element numerical schemes [11], [18], [31]. A comprehensive review of these models is given by Adachi et al. [1]. Griffith’s criterion is the common method to model fracture propagation in all of these techniques. Fracture propagation in Griffith’s criterion is a function of stress intensity factor and rock toughness. Griffith’s criterion presumes the presence of an initial fracture and predicts its propagation, hence it is an appropriate method to model major hydraulic fracture propagation, but it cannot predict fracture nucleation. Since, we are interested in predicting fracture initiation on the surface of intact rock or along the cemented natural fractures; we cannot limit our analysis to Griffith’s criterion. We use the cohesive interface technique to model reopening of cemented fractures. The cohesive interface model is a constitutive equation to model deformation of discontinuities, which can be easily applied to multiple cracks or incorporates their coalescence. Cohesive interfacial model could also be used to simulate fracture propagation with the advantage of removing stress singularity at the fracture tips [2]. Later laboratory experiments showed that nonlinear region added to cohesive crack models provides better prediction for fracture growth in granular cementious materials like rock and concrete [6]. The cohesive interfacial model considers a cohesive crack of zero width with traction transferring capacity, thus eliminates the stress singularity problem at the crack-tip. Additionally by nature, cohesive interface concept is the best fit for the problems with predefined fracture propagation paths like this problem; however, some sophisticated algorithms has been invented to adaptively add or remove cohesive elements in the computational model upon necessity [47]. In addition, cohesive interfacial models, despite their nonlinear nature, are easy to implement and we will see in the next section how we used this capability to model the initiation of microcracks during hydraulic fracturing.

Cohesive interface model is mainly a nonlinear constitutive equation between the traction and jump in displacement between two bodies. Cohesive interface starts to fail when the applied traction exceeds a critical value and followed by softening, and both are described by traction separation law [43]. Traction-separation law has the flexibility to tune the parameters to find potential function solutions for fracture propagation in different materials [24], [29]. Bilinear law is the simplest form of traction separation law composed of two piecewise linear sections for pre- and post-failure situations (Figure 4). Despite its simplicity as constitutive equation, it has proven capabilities to model fracture behaviour in cementious materials [6]. Quadratic stress law is good candidate for mixed mode condition. In this model, damage initiates the quadratic interaction function involving nominal stress ratios reaching the value of one (Equation 7), while t_{n} and t_{s} represents the real values of normal and tangential tractions along the interface, respectively.

Fracture energy release, cohesive strength, initial cohesive stiffness, critical separation gap at complete failure and critical separation at damage initiation are key parameters to specify irreversible fracturing based on bilinear cohesive law. To model the mixed mode fracture propagation, Benzeggagh-Kenane (BK) fracture criterion is used [7].

For an incompressible and Newtonian fluid, equation (9) represents the continuity equation where the first term shows fracture capacity based on its width change and the second term represents the cross-sectional flow rate of the fracture (no leak-off from fracture into the formation). The tangential flow along the gap between two cohesive walls base on momentum equation for Poiseuille’s flow pattern is represented in equation 10 where q, w, p and μ are local flow rate, local crack width, fluid pressure inside the fracture and fracturing fluid viscosity, respectively [4].

Analytical solutions for coupled thermo- poroelasticity and only restricted to simple geometries. Therefore, they are not pertinent for a system with numerous factures. Therefore, we used a commercial finite element package, ABAQUS, for modelling purposes.

## 5. Results and discussion

Commercial finite element software, ABAQUS (Dassault Systèmes Simulia Corp.) [20], is chosen for implementing cohesive crack methods to model secondary fractures initiation and propagation. We begin with the simplest possible geometry for the hydraulic fracture, i.e. a planar fracture. Extension of utilized techniques to non-planar hydraulic fractures does not require introducing any new concept and should be straightforward. Due to the symmetry of the problem with respect to the fracture plane, we only need to simulate half of the geometry. Figure 5 shows the numerical grid with the blue zone showing fracture surface, and the surrounding red zone showing the intact rock. To model pre-existing cemented natural fractures, cohesive elements have been embedded as parallel planes perpendicular to the fracture surface (Z-direction) with 5 cm spacing for this example. Hence, fluid pressure during pumping stage will be introduced only to the fracture surface (blue zone), and the rest of the model will be under the effect of in-situ stress only. In case that any part of the natural fractures (cohesive elements) reaches failure threshold, following the opening of the crack, fracturing fluid is supposed to reach the opened part of the fracture and pressurized it, which is considered by removing failed cohesive elements and adding gap flow to our model to include fracturing fluid pressure and their cooling effect in natural fractures. Our preliminary models showed the significance of this effect on clustering of secondary fractures and their depth of penetration.

The model assumed homogeneous and isotropic properties for mechanical and hydraulic properties of the rock. Additionally, fracturing of the rock is assumed as an irreversible process. Furthermore, we considered development of major fracture parallel to maximum principal stresses, and natural fractures are assumed to be fully cemented with digenetic cements like quartz or calcite. We used Drucker-Prager model to describe plastic behaviour of rock. A summary of mechanical properties considered for the rock is provided in Table1.

To estimate the parameters of the bilinear cohesive law, no new experimental step is needed. The classic lab tests to measure tensile strength and critical energy release rate should be enough to derive bilinear cohesive law parameters. Among various fracture testing techniques available for homogeneous materials, those which facilitate the stable advance of a fracture are more preferred for interface toughness (G_{c}) measurements. Examples include the double cantilever crack specimen [26] and Brazilian disk. The tensile strength (σ_{max}) may also be measured using common Brazilian beam tests. Based on the definition of crack energy release rate and bilinear softening, the maximum separation at failure can be determined:

Some parametric studies have been done on both cohesive parameters to observe their effects on the pattern of reactivated fractures; however, both parameters are basically a function of the composition of digenetic materials and environmental conditions at the time of their precipitation. Obviously, the values of tensile strength and energy release rate for fracture cements are much lower than the values for rock matrix. Strength of a cemented natural fracture is a function of cement type (composition) and its continuity [27]. Gale et al. [28] evaluated the rupture strength in Barnett shale at different depths for different lithofacies. Based on their published laboratory tests, we chose cohesive properties and other mechanical properties of rock. The values of 12 MPa and 3.2 Pa.m are considered for rupture strength and fracture toughness of cemented fractures, respectively. We assumed that net fracture fluid pressure, applied on blue zone in figure 5, is gradually increased to reach 2 MPa and then slowly bleed off to become equal to reservoir pressure. The induced tensile stress due to loading and unloading process during hydraulic fracturing can reactivate pre-existing fractures. Figures 6 and 7 show that the induced tensile stresses due to fluid pressure decline inside fracture can be as important as induced stresses in the loading process. Figure 6 shows reactivated fractures at different depth from the fracture surface at the peak of fracturing fluid pressure and before pressure decline due to leakoff; main mechanism for failure is shear associated with compressive stresses. However, the main mechanism for failure in Figure 7 is associated with residual tensile stresses induced in the unloading due to plasticity. Fractures in Figure 8 are reactivated due to not only plastic effect but also considering a temperature difference between the rock matrix and the hydraulic fracturing fluid.

Table 3 presents the effect of fracture cement strength, fracture toughness and cement resilience on the pattern of opened fractures i.e. failed cohesive elements. By increasing cement strength, the number of initiated cracks on the cohesive layer decreases (Figure 9 (a) and (b)); in addition, the decrease in cement toughness by nearly one-fifth makes considerable increase in the initiated cracks (Figure 9 (a) and (h)). The values for the failed cohesive elements in Table 3 show that flaws initiation has noticeable dependency on cohesive stiffness. Small decrease in the cohesive stiffness leads to considerable increase in the cohesive damage, indicating increase of contact area between rock-fluid. As it showed in Figures 6 and 7, loading-unloading residual stress may be considerably effective in reactivating pre-existing natural fractures; moreover, it may have a significant influence on fracture reactivation when assisted by induced thermal stresses due to temperature difference between hydraulic fracturing fluid and formation rock (Figure 9).

In Figures 8 and 9, initial temperature of 200°C for the rock is selected, and the expansion coefficient of rock matrix is assumed to be 1.5^{x}10^{-5}. The surface of hydraulic fracture is exposed to hydraulic fracturing fluid with temperature of 150°C. The temperature difference between hydraulic fracturing fluid and rock induced tensile stresses on the rock matrix. This tensile stress can be intensified by the induced tensile stresses due to elastic unloading of hydraulic fracture plastic deformation as shown in Figure 8. Figure 9 presents the pattern of reactivated natural fractures of Table 3 at the depth on-third of cohesive layer depth far away from the surface of hydraulic fracture.

Young’s modulus | 26 GPa |

Rock density | 2100 Kg/m^{3} |

Rock friction angle | 30° |

Rock dilation angle | 20° |

Poisson’s ratio | 0.27 |

Rock yield stress | 30MPa |

Case | Tensile strength (MPa) | Cement toughness (Pa.m) | Hydraulic pressure (MPa) | Cohesive stiffness (GPa) |

1 | 12 | 3.2 | 2 | 6.4 |

2 | 1.2 | 3.2 | 2 | 6.4 |

3 | 12 | 3.2 | 2 | 0.64 |

4 | 12 | 32 | 2 | 6.4 |

5 | 20 | 3.2 | 2 | 6.4 |

6 | 12 | 3.2 | 2 | 1 |

7 | 12 | 3.2 | 2 | 3.4 |

8 | 12 | 15 | 6.4 | 2 |

Case | Total cohesive elements | Failed cohesive elements after loading | Failed cohesive elements after unloading | Failed cohesive elements by Thermoplasticity | Studied parameter |

1 | 14700 | 3406 | 4692 | 5511 | Reference case |

2 | 14700 | 0 | 0 | 7923 | Tensile strength |

3 | 14700 | 0 | 0 | 0 | Cohesive Stiffness |

4 | 14700 | 0 | 0 | 1 | Cement toughness |

5 | 14700 | 415 | 1247 | 1683 | Tensile strength |

6 | 14700 | 0 | 0 | 0 | Cohesive Stiffness |

7 | 14700 | 18 | 74 | 154 | Cohesive Stiffness |

8 | 14700 | 298 | 484 | 2440 | Cement toughness |

## 6. Conclusion

Induced tensile stresses are expected to occur in rocks with significant plastic behaviour during loading and unloading of hydraulic fractures. These induced stresses can open pre-existing natural fractures in the formation and even open the cemented natural fractures. These activated fractures can provide more rock-fluid contact area. The size of these natural fractures is much less than the main hydraulic fracture but presence of these fractures in considerable numbers can significantly increase the contact area between the wellbore and formation. The path of initiation and propagation of these fractures can be induced by natural fractures. To study their effect, a three-dimensional finite element model with cohesive interfaces embedded in the rock. The effect of energy release rate and cohesive tensile strength investigated and it was shown that the decrease of these parameters can activate more natural fractures but not in a uniform pattern. Moreover, the simulation showed that the effect of plasticity can be more considerable when it would be helped by thermal stress induced by temperature difference between rock matrix and fracturing fluid.