1 Introduction

The propagation direction of hydraulically induced fracture was contested at some point (Harrison et al., 1954; Reynolds et al., 1961) whether it is horizontal along a bedding plane lifting the overburden (Howard & Fast, 1950) or orthogonal to the minimum principal stress (Hubbert & Willis, 1957). However, a number of theoretical analyses and experiments demonstrated that the plane of hydraulic fracture propagation is orthogonal to the minimum principal stress (Harrison et al., 1954; Hubbert & Willis, 1957). Based on this notion, hydraulic fracturing has become an established means to estimate the minimum principal stress as the fracture planes open against the minimum principal stress (Abou-Sayed et al., 1978; Baumgärtner & Zoback, 1989; Haimson & Fairhurst, 1969) or even the intermediate principal stress (Schmitt & Zoback, 1989), though debates on pressure interpretation are yet to be settled (Bröker & Ma, 2022; F. Guo et al., 1993; Jung et al., 2016; Kamali & Ghassemi, 2019). A thorough review on this subject is provided by Schmitt and Haimson (2017).

However, when hydraulic fractures interact with preexisting discontinuities (e.g., natural fractures or defects in the rock mass), the propagation of fractures may deviate from the direction orthogonal to the minimum principal stress. Microseismic records from hydraulic fracturing of tight hydrocarbon formation show that clouds of seismic events extended beyond the plane orthogonal to the principal minimum stress direction (Fischer et al., 2008; Warpinski & Teufel, 1987). Similar observations were also made in geothermal well stimulation (Cladouhos et al., 2016; Cuenot et al., 2008; Kraft & Deichmann, 2014). Furthermore, mine–back experiments (Jeffrey & Weber, 1994; Jeffrey et al., 2009) and retrieved cores from the in situ hydraulic fracturing test (Ciezobka et al., 2018) confirmed such fracture deviation and interaction with discontinuities. Understanding of such complex fracture behaviors is crucial for well stimulation whose main goal is to increase formation permeability in the near wellbore (Riffault et al., 2018). And to achieve that, induced hydraulic fractures are to interact with preexisting discontinuities to increase their morphological complexities (Cipolla et al., 2008; S. Li et al., 2019).

To study the influence of the preexisting fracture geometry and the applied confining pressure on the direction of hydraulic fracture propagation, Lamont and Jessen (1963) conducted hydraulic fracture experiments with a preexisting fracture. They concluded that the fracture propagation direction is mostly orthogonal to the minimum principal stress and deviates only slightly impacted by the preexisting fracture. Their findings are not aligned with many of the studies reported later. One possible explanation could be that the injection rate, which was not reported in Lamont and Jessen (1963), may have been too high for stable fracture propagation (Q. Zhang et al., 2021). Their study was followed by many experimental works that focused on impacts of the approaching angle, differential stress and strength of the preexisting fracture on the hydraulic fracture propagation (Additionally, injection rate and fluid viscosity impact hydraulic fracture interaction with a preexisting fracture for viscous dominated fracture propagation (Chuprakov et al., 2014; Llanos et al., 2017). In this study, however, our focus is on toughness dominated fracture propagation. For fluid effects, see Section 12.) (Blanton, 1982; Dehghan et al., 2015; T. Guo et al., 2014; Lee et al., 2015; Meng et al., 2010; Warpinski & Teufel, 1987; Yushi et al., 2016; Zhao et al., 2022; Zhou et al., 2008). Q. Zhang et al. (2021) provides a comprehensive review on experimental studies. Furthermore, hydraulic fracture interaction has been studied by numerically (Dahi-Taleghani & Olson, 2011; Fatahi et al., 2017; Lepillier et al., 2020; McClure & Horne, 2014; K. Wu & Olson, 2016). We refer to Dahi-Taleghani et al. (2016) and Lecampion et al. (2018) for reviews on numerical aspects on this subject.

Most existing studies have reported that interactions of hydraulic fractures with discontinuities are largely controlled by the in situ stress, and the geometry and strengths of the discontinuity. Using these parameters, studies have formulated criteria for fracture penetration/deflection in terms of the stress (Blanton, 1982; Chuprakov et al., 2014; Fu et al., 2018; Gu et al., 2012; S. Li et al., 2017; Renshaw & Pollard, 1995; Warpinski & Teufel, 1987) or the energy (Dahi-Taleghani & Olson, 2014; X. Li et al., 2020; Zeng & Wei, 2017). These studies have established that a fracture is more likely to deflect at a discontinuity if either of these decreases: (a) the approaching angle to the interface, (b) the interface strength, or (c) the differential stress.

However, these criteria may not apply straightforwardly to the behaviors of hydraulic fractures in the near wellbore region. Grain scale images from hydraulic fracturing experiments on carbonates (Naoi et al., 2020) have shown that hydraulic fractures deflect more frequently at discontinuities (grain boundaries) as they propagate farther from the wellbore (Figure 1). We can postulate that weakly bonded grain boundaries act as defects and cause the fractures to deflect, as reported in other experimental studies on granites (Chen et al., 2015), but we cannot deduce that grain boundaries are coincidentally weaker farther from the wellbore.

Though Naoi et al. (2020) remarked the increasing fracture complexity away from the wellbore, their study focused on the failure mechanisms from moment tensor analyses and no quantitative analyses of fracture morphology or computation on fracture propagation was performed. In this study, our objective is to investigate their experimental observations from an energetic perspective. Our contributions are three-hold. First, we analyzed the fracture paths quantitatively using digitized images from the experiments in Naoi et al. (2020) including unpublished results. Second, we computed energy release rates of deflecting and penetrating fractures with varying distances from the wellbore using the potential energy differences induced by fracture tip perturbations. Third, we numerically simulated hydraulic fracture interactions with a discontinuity (grain boundary) using an approach based on successive energy minimization (Bourdin et al., 2012). The results of the energy release rate analysis and the numerical simulations hydraulic fracture propagation agree well with the experimental observations, indicating that the fracture propagation path is affected not only by the in situ stress, and the geometry and strength of discontinuities (grain boundaries), but also by the distance from the wellbore.

2 Fracture Path Analysis

We analyzed fracture paths observed in three samples from the hydraulic experiments conducted in Naoi et al. (2020). For experimental and property measurement details, we refer to Supporting Information S1 and only brief experimental procedures are provided here.

2.1 Experimental Procedures

Hydraulic fracture experiments were conducted on outcrop samples of Eagle Ford formation, a carbonate rock from Langtry, Texas, USA. The samples have a representative Young's modulus (E) of 50 GPa and a Poisson's ratio (ν) of 0.33 measured perpendicular to the direction of sedimentary plane. The porosity is estimated to be lower than 1%. X-ray diffraction analysis indicates that they consist of 96.3% calcite and 3.7% quartz. Although their mineral composition is nearly homogeneous, the samples include textural heterogeneity in the form of calcite microfossils, as seen in thin sections (Figure 2a). These inclusions are expected to have similar mechanical properties to the bulk rock, except that the grain boundaries may not be bonded perfectly.

The total of three samples were prepared with dimensions of 65 mm × 65 mm × 130 mm ([−32.5 mm, +32.5 mm] × [−32.5 mm, +32.5 mm] × [−65 mm, +65 mm]) with a wellbore 3 mm in radius drilled in the center (Figure 2b). The top and bottom surfaces were subjected to an axial loading of 5 MPa, and the other surfaces were not loaded. The wellbore was sealed with a packer which has a 30 mm pressurizing section (x = −15 to +15 mm) sealed with two O-rings. Through this open section, a thermosetting acrylic (methyl methacrylate) resin with a viscosity (It is comparable to the water viscosity at 20°C (1.0 mPa⋅s).) of 0.8 mPa⋅s was injected to induce fractures. The resin contained a fluorescent compound so that hydraulically induced fractures can be distinguished from preexisting fractures by their brightness in a thin section under ultraviolet light (Chen et al., 2018; Ishida et al., 2016) (Figure 1). And the limit of detection is reported to be less than a few micrometers (Nishiyama & Kusuda, 1994).

Figure 2c shows the response of pressures and acoustic emissions from one of the experiments (sample #1). The hydraulic pressure increased almost linearly and dropped abruptly at the onset of fracturing. This brittle and linear elastic behavior is consistent with the concentration of acoustic emissions around the point of failure. The peak (breakdown) pressures from sample #1, #2, and #3 were 11.79, 10.30, and 10.71 MPa respectively. The fractured samples were heated immediately to harden the resin (Chen et al., 2018). The thin sections were then prepared at the central cross-section (i.e., x = 0.0 mm) and we analyzed fractures propagating roughly along the z-axis.

2.2 Fracture Propagation Paths

We examined the thin sections of the experimental specimens under ultraviolet light and studied the fractures connected to the wellbore, classifying long fractures with bright fluorescence as main fractures and short, deflected fractures as branch fractures (Figure 3). When the whole fracture is connected, we picked the brighter fluorescence (wider fracture aperture) as main. Finally, all the identified fractures were visually confirmed as some calcites also have fluorescent properties.

In these images, the fractures appear to have increased in topological complexity as they propagated farther from the wellbore. To investigate this apparent change, we calculated the main fracture lengths (L_m) and the total fracture lengths (L_t) against the distance from the wellbore (d). L_t is the sum of the branch and main fracture lengths and thus L_t ≥ L_m. As shown in Figure 4, we divided the images laterally into 10 μm intervals in the direction of fracture propagation, measured the fracture lengths within each interval, and then determined the total fracture lengths by summing all of the intervals (see Appendix A for the detailed procedure). For a measure of fracture complexity, we normalized the total fracture length against the corresponding main fracture length as L_t/L_m. This measure has the value 1 (L_t/L_m = 1) until the first deflection, after which the fracture complexity increases if branch fractures take longer paths than the main fracture and decreases otherwise.

When the main fracture lengths L_m are plotted against the normalized distance from the wellbore (d/r_o, where r_o is the wellbore radius of 3 mm), the main fractures all fall on a nearly identical linear trend (Figure 5a). When the normalized total fracture lengths L_t/L_m are plotted against d/r_o (Figure 5b), the fracture complexities initially increase with d/r_o to values as high as about 4.5, but they reach a plateau of complexity at a distance between one-half and twice the wellbore radius (d/r_o = 0.5 − 2 or d = 1.5 − 6 mm).

3 Fracture Propagation Mechanism

3.1 Revisiting the Griffith Criterion

The Griffith criterion (Griffith, 1921) states that a line fracture of length l propagates when the energy release rate reaches a critical value,

(1)

where G_c is the critical energy release rate and is the energy release rate, defined as the released potential energy (Π) with respect to fracture growth:

(2)

While G_c is a measure of material resistance against fracturing (material property), is dependent on the fracture geometry and applied loadings. The Griffith criterion is complete with the irreversibility of fracture, meaning that the fracture length only increases. With more recent mathematical tools, we can view Griffith's arguments as a constraint minimization problem known as the Karush-Kuhn-Tucker conditions (Nocedal & Wright, 2006). By denoting the rate of length growth as , we can rewrite the Griffith criterion as,

(3)

Francfort and Marigo (1998) recast the Griffith criterion Equation 1 as a minimization problem where the displacement u and the crack set Γ are to be determined from

(4)

where the total energy is the sum of the potential energy and the fracture surface energy:

(5)

From Equation 4, we can view the Griffith criterion as a special instance of the case in which the energy is minimized only along the prescribed fracture path l:

(6)

Then we can recover the Griffith criterion from the first-order minimization condition,

(7)

3.2 The Griffith Criterion in Terms of Stress

The Griffith criterion Equation 1 is often applied in terms of stress using the concept of stress intensity factor, K, through Irwin's formula (Irwin, 1957). Here we note a few identities between the energy-based (Griffith's) and stress-based (Irwin's) approaches.

For a line crack, , subjected to a tensile (mode-I) loading in an infinite 2D domain, the energy release rate defined in Equation 2 can be equivalently expressed using the stress intensity factor as,

(8)

where E^′ = E/(1 − ν²) for plane strain and E^′ = E for plane stress, and K_I is the (mode-I) stress intensity factor, defined in polar coordinates (r, θ) as (Anderson, 1995),

(9)

Accordingly, Griffith's fracture propagation criterion in Equation 1 can be written using the stress intensity factor as

(10)

where K_Ic is the critical stress intensity factor for mode-I fracture. As the energy release rate is expressed in terms of the mode-I stress intensity factor (K_I) in Equation 8, the critical energy release rate (G_c) can be expressed in terms of the mode-I stress intensity factor (K_Ic) as

(11)

Both G_c and K_Ic represent the material's resistance to fracturing and are related through Equation 11. And, interchangeably, they are often referred to as “fracture toughness.”

In the following, we analyze fracture behaviors using the energy release rate rather than the stress intensity factor(s), but the energy release rate can be converted from the stress intensity factors. More generally, when all modes of loadings (I, II, and III) are present, the energy release rate is,

(12)

where K_II and K_III are the stress intensity factors for shear (mode-II) and tearing (mode-III) loadings, respectively, and μ is the shear modulus. From Equation 12, we can see that an increase in the energy release rate means an increase in the stress intensity factor and vice versa.

3.3 Fracture Path Determination

In which direction will a fracture propagate under multidirectional loadings? The two main approaches to this question are the principle of local symmetry (Cotterell & Rice, 1980; Gol'dstein & Salganik, 1974) and the maximum energy release rate (Hussain et al., 1974; Nuismer, 1975). Under mode-I (tensile) and mode-II (shear) loadings, for example, the principle of local symmetry mandates model-I propagation, where the mode-II stress intensity factor is zero (K_II = 0), whereas the maximum energy release rate approach seeks to determine the direction that maximizes the energy release rate.

These two approaches predict a similar propagation direction (identical for smoothly growing cracks and almost identical for kinking cracks) in isotropic, homogeneous materials (Amestoy & Leblond, 1992; Cotterell & Rice, 1980; Hutchinson & Suo, 1991; Nuismer, 1975), but the principle of local symmetry does not correctly predict propagation paths in anisotropic or heterogeneous materials (Brach et al., 2019; Marder, 2004; Takei et al., 2013). In contrast, the maximum energy release rate approach can predict the right direction if we consider the energy dissipation associated with a material heterogeneity (Chambolle et al., 2009; Gurtin & Podio-Guidugli, 1998; Hakim & Karma, 2009; He & Hutchinson, 1989; B. Li & Maurini, 2019; Takei et al., 2013). This modified (or generalized) energy release rate approach seeks a direction that maximizes the energy dissipation as

(13)

where α is the propagation angle we wish to find.

Let us consider a fracture arriving at an inclined interface (e.g., a natural fracture) with an angle of β (Figure 6a), where the interface has a different fracture toughness from the bulk material given by

(14)

where and are the fracture toughness of the interface and the bulk material, respectively.

He and Hutchinson (1989) computed the energy release rate of a deflecting fracture under mode-I loading and normalized it against the energy release rate of a straight fracture propagation. The normalized energy release rate depends only on the deflection angle α:

(15)

which is plotted in Figure 6b against α. We can see from Equation 15 that the energy release rate is always the greatest with straight fracture propagation . Thus, the energy dissipation of straight propagation (α = 0) is always greater than those of other propagation angles except for the direction of the interface (α = β), which leads to:

(16)

This result leads us to two possibilities, α = β or 0. And whether the energy dissipation along the interface, , is greater than that of straight propagation, , depends on and . Using this criterion, we can determine whether the fracture deflects or penetrates by simply comparing the respective energy dissipations of deflection and straight propagation (Xu et al., 2003). In other words,

(17)

Using Figure 6b, we can determine this graphically. Given an interface angle β and fracture toughness , if the point corresponding to falls below the He and Hutchinson criterion Equation 15, then the fracture deflects; otherwise, the fracture penetrates.

4 Results

4.1 Energy Release Rates of a Deflecting Fracture

To examine how energy release rates change with the distance from the wellbore, we computed the energy release rates of penetrating and deflecting fractures at a grain boundary in our experimental setting. In a half-domain of the sample in plane strain, we considered a straight hydraulic fracture arriving at the center of a grain located at a distance d away from the wellbore (Figure 7). Because the fluid viscosity was low (0.8 mPa⋅s) and the fracture scale was small (∼6 mm), we neglected the pressure loss within the fracture and applied the same hydraulic pressure (15 MPa) to the wellbore's internal wall and the fracture surfaces. Because all mineral grains were calcite, Young's modulus and Poisson's ratio of the inclusion were considered identical to those of the bulk rock. Table 1 lists the physical and geometrical parameters used in the analyses. All the parameters were taken from Naoi et al. (2020) except for the fracture toughness, which was not measured, and we used a typical fracture toughness value for carbonates (10 Pa⋅m). The inclusion radius (0.1 mm) is a typical size observed in the samples.

Table 1. Material Properties and Geometrical and Loading Parameters

Name	Symbol	Value	Unit
Young's modulus	E	50	GPa
Fracture toughness of bulk		10	Pa⋅m
Poisson's ratio	ν	0.33	−
Axial loading	–	5	MPa
Wellbore pressure	–	15	MPa
Width	W	65	mm
Height	H	65	mm
Wellbore radius	ro	3	mm
Inclusion radius	–	0.1	mm

In the absence of an analytical expression for this particular setting, we numerically approximated energy release rates Equation 2 as,

(18)

We perturbed the crack tip by a small increment (Δl), and then computed the corresponding potential energy (Π) at each increment. To represent a penetrating crack, we extended the tip straight, and for a deflecting crack, we extended the tip by splitting it into two tips deflected by 90° in opposite directions (Figure 7) (see Appendix B for the computational details and verification examples).

In accordance with the deflection/penetration criteria in Figure 6, we normalized the energy release rates of a deflecting fracture against the energy release rate of a penetrating fracture . This normalized energy release rate represents the likelihood of fracture deflection; the higher it is, the more likely a fracture deflects for the same fracture toughness. Figure 8 plots the normalized energy release rates against the normalized distance from the wellbore (d/r_o). The energy release, or the likelihood, increases with the distance from the borehole. Recalling the relationship between energy release rate and stress intensity factor Equation 8, this trend is equivalent to an increase in stress intensity factors with the distance. This increasing trend continues to about twice the wellbore radius (d/r_o ∼ 2).

4.2 Simulation of Hydraulic Fracture Interactions at a Grain Boundary

To examine the effect of grain boundary strength and location on hydraulic fracture behavior, we simulated the propagation of a hydraulic fracture and its interaction with a mineral grain using a variational phase-field model (Bourdin et al., 2012). The model is based on minimization of the Francfort-Marigo energy Equation 4 with phase-field regularization. Details are given in Appendix C.

Again, given the low fluid viscosity and small fracture scale, we neglected pressure loss within fractures. Furthermore, because of the low sample permeability (around 1 nanodarcy, see Supporting Information S1), we neglected fluid leak-off to and poro-elastic deformation in the rock formation. The injection rate was set at 1 cc min⁻¹. (Because pressure loss and leak-off are ignored, the injection rate is simply the fracture volume increase in one time step and has no other physical impact.) A mineral grain was placed at a distance d from the wellbore and, for numerical simplicity, a short fracture, , was prescribed. We represented the grain boundary by assigning a weaker fracture toughness to the interface than that of the bulk rock (i.e., ) while keeping their elastic properties identical.

We considered three different fracture toughnesses at the grain boundary: strong , moderate , and weak . We also considered three different grain positions relative to the wellbore: near (d = 1/3r_o), middle (d = 2/3r_o), and far (d = 3/3r_o). Our simulations show that the fracture penetrates the grain without deflecting at the grain boundary for all toughness cases if the grain is in the near position (Figure 9a). With the grain in the middle position (Figure 9b), the fracture deflects in the weak boundary toughness case, but not in the other cases. With the grain in the far position (Figure 9c), the fracture deflects in the moderate and weak cases, but not in the strong case.

5 Discussion

Our work shows that the energy release rate of a deflecting fracture increases with the distance from the wellbore, which means that the farther a fracture propagates away from the wellbore, the more likely it is to deflect at a defect in the rock mass. According to our analysis, the likelihood of fracture deflection continues to increase up to about twice the wellbore radius (2r_o) if all other things are equal (Figure 8). Because grains and defects in real rocks are not distributed uniformly and their boundary toughnesses vary, our experiments exhibited varying degrees of complexity, but deflections peaked somewhere between d/r_o = 0.5 and d/r_o = 2 (Figure 5b), a result remarkably close to our analysis. Simulations with various grain locations and boundary toughnesses confirmed two things: (a) for grains at a given distance, the fracture deflects more easily at a weaker grain boundary, and (b) for grains with the same boundary toughness, the fracture deflects more easily at a grain more distant from the wellbore.

The first claim seems intuitive, but the second may not. To rationalize the second claim, let us consider the disturbance that the fluid pressure imposes on the stress field. Fluid pressure acts on both the wellbore walls and the fracture surfaces. While the fluid pressure inside the fracture imposes a tensile loading at the crack tip, the fluid pressure on the wellbore induces a loading parallel to the main fracture as depicted in Figure 7. This parallel loading to the propagation direction then suppresses the tendency of the fracture to deviate from the propagation direction. Considering Kirsch's solution in an infinite domain under uniaxial loading of σ_H (Kirsch, 1898) and Lame's solution for an infinitely thick wall cylinder (Timoshenko & Goodier, 1970, p. 69), the radial stress in the direction parallel to σ_H is,

(19)

where p_w is the wellbore pressure. The radial stress σ_r acts against the fracture deflection near the fracture tip. In accordance with our analysis and experiment, let us set p_w = 3σ_H (p_w = 15 MPa, σ_H = 5 MPa). Then we have

(20)

The formulation in Equation 20 shows that the radial stress σ_r decreases rapidly with the distance from the wellbore; from the wellbore r = r_o (d = 0) to r = 3r_o (d = 2r_o), the normalized radial stress (σ_r/p_w) drops by 60%. Thus, a fracture will be less suppressed against deflection as it propagates farther from the wellbore. This is an overly simplified analysis of our experiments ignoring the finite domain and pressurized fractures, but it helps understand the implications of our energy release rate analyses. We also note that a rigorous analysis of stresses around the fracture tip with a discontinuity is a daunting and, unfortunately, futile effort because approaches based on stress intensity factors will fall short of predicting the fracture path in anisotropic or heterogeneous media, as discussed in Section 7.

Most previous studies of fracture propagation criteria have focused on the interplay between the in situ stress and the geometry and strength of natural fractures. However, our study demonstrates that the presence of a pressurized wellbore appreciably affects the fracture interactions with grain boundary discontinuities in the near wellbore region. Although the pressurized wellbore impact has been considered in the past in the context of tortuous fracture propagation when a fracture is initiated in a direction not orthogonal to the minimum in situ stress (Feng & Gray, 2018; X. Zhang et al., 2011), to the best of our knowledge, this study is the first to consider the effects of pressurized wellbore on fracture penetration/deflection at a grain boundary discontinuity in the near wellbore region. Our findings also underscore that analysis should address not only the local stress stability around the fracture tip but also the global energetic field which takes into account the presence of other fractures or a wellbore. Thus, we would caution against uncritically transferring criteria for fracture deflection/penetration from one setting to another, especially when other wells or fractures are present.

Our energy release rate analyses in Section 9 are based on linear elastic fracture mechanics, which is applied to describe fracture propagation in a wide range of scale sizes from 10⁻⁶ m (Shimada et al., 2015) to 10² m (Detournay, 2016). At our experimental scale, a grain boundary acts as a discontinuity. At larger scales, however, discontinuities may be defects in the rock mass or preexisting fractures, and our findings on near wellbore fracture behaviors may apply not just to grain boundaries, but to these discontinuities as well. We also note that our analysis considered a pressurized wellbore, which is applicable to open-hole completions but does not apply to cased and perforated completions, in which the casing will prevent the wellbore pressure force from directly acting on the rock formation. However, if the cementing of the casing is not perfect, then seeping fluid may pressurize the interface between the cement and the rock formation, bringing the conditions closer to those of this study. Conversely, if a rock–cement interface is tightly sealed, induced fractures will likely interact more strongly with grain boundary discontinuities in the near wellbore region shown in Figure 8 (d = 2r_o). Alternatively, if perforations extend past d = 2r_o, then fractures issuing from the perforations may become complex. Both situations should increase the fracture complexity immediately around the wellbore and result in improved well productivity.

Among all the fractures analyzed (Figure 4), the upper fracture in sample # 1 exhibits a much greater normalized length (∼4.5) than the others. In terms of the mineral grain distribution, we did not observe anything particular compared to the other samples. On the other hand, the recorded breakdown pressure of sample #1 was higher (11.09 MPa) than the other 2 samples (10.30 and 10.71 MPa). One possible explanation for the high degree of complexity is that the rock mass had a higher fracture toughness and induced a higher breakdown pressure which led to faster fracture propagation. This then may have caused more unstable and more complex fracture propagation. To investigate this claim, however, dynamic fracture mechanics needs to be considered and may be a topic of future study.

Furthermore, in this study, we focused only on the very first deflection of a fracture, when a straight fracture arrives at a grain boundary discontinuity. Fracture interactions with discontinuities after the first deflection (non-straight fracture) may be important subjects for future studies. Additionally, we did not consider fluid effects, although we can no longer neglect pressure loss within the fracture or fluid leak-off to the rock formation under some circumstances. For example, experimental studies have reported increasing hydraulic fracture complexities with the use of a low–viscosity fluid (Ishida et al., 2016; Watanabe et al., 2017). These effects too should be considered in future studies.

6 Conclusions

We have studied a phenomenon repeatedly observed in our hydraulic fracturing experiments: the increasing tendency of fracture deflection with the distance from the wellbore. We were able to quantitatively explain this increasing complexity through an analysis on the energy release rates. Our numerical simulation results also report an increasing likelihood of fracture deflection at grain boundaries with the distance from the wellbore. Both computations were based on energy minimization and agreed with our experimental observations. From these findings, we conclude that fracture paths are better determined on the basis of global energy dissipation rather than local stresses around the fracture tip, especially in anisotropic or heterogeneous formations.