Interactions of Hydraulic Fractures With Grain Boundary Discontinuities in the Near Wellbore Region

Hydraulic fractures often turn or branch, interacting with preexisting discontinuities in the rock mass (e.g., natural fractures or defects). The criteria for fracture penetration or deflection are typically based on the in situ stress, and the angle and strength of discontinuities. However, in hydraulic fracture experiments on carbonate rocks (Naoi et al., 2020, https://doi.org/10.1093/gji/ggaa183), small scale analyses show that the fractures deflected more frequently at discontinuities (grain boundaries) as they propagated farther from the wellbore, a finding not explained by the conventional criteria. Here, we demonstrate that the energy dissipation of a deflecting fracture increases with the distance from the wellbore, such that a propagating hydraulic fracture more easily deflects at a discontinuity from an energetic standpoint. This tendency was confirmed by hydraulic fracture simulations based on a successive energy minimization approach. Our findings, which show that wellbores appreciably affect the behavior of hydraulic fractures, highlight the importance of energetic stability analysis for determining fracture paths.

. 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 5.) (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.

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.

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  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.

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 Sample #1, #2, and #3 were called EFS1701, EFS1704, and EFS1706 respectively by Naoi et al. (2020). "Upper" represent the upper half of a sample and "lower" the lower half of a sample.
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).

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, 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: 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 ̇≥ 0 , we can rewrite the Griffith criterion as, 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 where the total energy  is the sum of the potential energy and the fracture surface energy: 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: Then we can recover the Griffith criterion from the first-order minimization condition,

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, Γ = [− ] × {0} , 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, where E ′ = E/(1 − ν 2 ) 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), Accordingly, Griffith's fracture propagation criterion in Equation 1 can be written using the stress intensity factor as 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 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, 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.

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 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 where int c and bulk c 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 (0) of a straight fracture propagation. The normalized energy release rate ( )∕(0) depends only on the deflection angle α: 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 ((0) ≥ ( )) . 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: This result leads us to two possibilities, α = β or 0. And whether the energy dissipation along the interface, ( )∕ int c , is greater than that of straight propagation, (0)∕ bulk c , depends on int c and bulk c . 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, Using Figure 6b, we can determine this graphically. Given an interface angle β and fracture toughness int c , if the point corresponding to

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 Figure 6. (a) A fracture arriving at an interface at angle β with a possible propagation angle α. (b) Energy release rates of a deflecting fracture normalized by the energy release rate of a penetrating fracture, plotted against the propagation angle α. The shaded region corresponds to fracture deflection and the unshaded region corresponds to fracture penetration; the line between them corresponds to the He and Hutchinson criterion (see the text for details). The point marked by "x" represents an interface angle of 20° and an interface fracture toughness of 30% of the bulk fracture toughness, which lies in the deflection region.
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.
In the absence of an analytical expression for this particular setting, we numerically approximated energy release rates Equation 2 as, 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 (d) against the energy release rate of a penetrating fracture ( p ) . 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).

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 −1 . (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, Γ = {0} × [ + 0.25 ] , was prescribed. We represented the grain boundary by assigning a weaker fracture toughness to the interface than that of the bulk rock (i.e., int c < bulk c ) while keeping their elastic properties identical.
We considered three different fracture toughnesses at the grain boundary: strong ( . 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.

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, 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 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 3.3.
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 4.1 are based on linear elastic fracture mechanics, which is applied to describe fracture propagation in a wide range of scale sizes from 10 −6 m (Shimada et al., 2015) to 10 2 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  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.

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.

Appendix A: Fracture Path Quantification
From the thin section images (Figure 3), we determined the fracture paths and lengths as follows ( Figure 4): 1. Register the main and branch fractures in sequences of pixels. 2. Set the beginning of the main fracture as the origin and assign the pixels to Cartesian coordinates (x j , y j ) by multiplying them by the image resolution (2.2 μm). 3. Compute the overall propagation direction vector for the main fracture, w = (cos θ, sin θ), by the least squares method, that is, where (x 1 , y 1 ) = 0 and N corresponds to the point when √ 2 + 2 equals the end of the image (18 mm).
4 In each 10 μm interval along the propagation direction w, compute the discrete fracture length l j assuming that the fractures are linear within each interval. 5 Integrate the discrete fracture lengths as = ∑ .
Steps 2-5 were repeated for L m and L t except that the propagation directions (θ in step 3) were only calculated for the main fractures.

Appendix B: Numerical Energy Release Rate Analysis
For computation of energy release rates, several integral-based approaches are available such as J-integral (Rice, 1968) or the G θ method (Dubois et al., 1998). However, for non-smooth deflection (kinking), these methods do not appear to compute the energy release rate accurately, even with reasonable element refinement. Therefore, we approximated the energy release rates using Equation 18 by evaluating potential energies from a linear elastic finite element analysis with perturbations in the fracture length. Theoretically, we can estimate  by running two linear elastic analyses with, say, Δl 1 and Δl 2 (Δl 2 > Δl 1 ) and computing the corresponding potential energies (Π 1 and Π 2 ) as: In this study, we computed the potential energy with five different fracture tip perturbation lengths for each geometrical configuration to ensure that (1) the perturbation of l is small enough and (2) the element discretization is adequate by checking whether the Π versus Δl curve falls on a straight line ( Figure B1).
The potential energy is written according to the Clapeyron theorem (Fosdick & Truskinovsky, 2004) with stress σ and strain ɛ as where ɛ is the linearized strain ( ∶= 1∕2 ( ∇ + ∇ T )) and u is the displacement. Figure B1 shows the computation of  for 80° fracture deflection in the He and Hutchinson problem (Figure 6a). To verify our computational procedure for  , we compared the numerically obtained values of ( )∕(0) with the analytical expression of He and Hutchinson (1989) in Figure B2a. The energy release rates computed with the G θ method are also plotted. Both approaches follow a similar decreasing trend, but as the deflection angle grows, the G θ method deviates from the analytical expression.
Furthermore, we verified computations of the absolute  value for a line fracture ([−a, +a] × {0}) pressurized by a constant pressure, p, in an infinite domain. Sneddon and Lowengrub (1969) derived an analytical solution to this problem: The computed energy release rates are in close agreement with the analytical solution with various fracture lengths ( Figure B2b).

C1. A Variational Phase-Field Model for Fracture
Although the Griffith criterion is presented in a generalized form in Equation 4, its numerical implementation is not straightforward because of the discontinuous and evolving crack set Γ. Among available methods (Schmidt et al., 2009;Wang & Sun, 2017), a regularization method via Gamma-convergence (Bourdin et al., 2000; Figure B2. Comparisons of the computation of an energy release rates with results from the solutions of (a) He and Hutchinson (1989) and (b) Sneddon and Lowengrub (1969). Chambolle, 2004) has become the standard implementation method for Equation 4 in the last decade (e.g., Ambati et al., 2014;Borden et al., 2012;Freddi & Royer-Carfagni, 2010) and is now typically referred to as the phase-field model.

C2. Extension to Hydraulic Fracturing
To The stress-strain relation is where ℂ is the tangential stiffness tensor, p p is the pore-pressure and α B is Biot's coefficient. In this study, because of the tight permeability of the sample and the short duration of injection, we assumed no leak-off. Thus p p is a temporary constant (p p (x, t) → p p (x)). Additionally, because of the small sample size, the pressure loss within the fracture is negligible. Thus we assume that p f is spatially uniform (p f (x, t) → p f (t)) and the pore-pressure is ambient (p p (x) → 0). Denoting p f (t) with p, the regularized total energy now reads: This model does not differentiate fracture evolution with respect to the sign of strain; that is, the fracture evolves identically under tension and compression (Amor et al., 2009). To avoid such unphysical behaviors, a number of different models that split the strain energy, the first term in Equation C4, have been proposed (Freddi & Royer-Carfagni, 2010;Miehe et al., 2010;Steinke & Kaliske, 2019;J.-Y. Wu, 2017). In this study, we adapted an approach proposed by Freddi and Royer-Carfagni (2010) in which the strain is split into positive and negative parts as where ɛ + is a positive definite tensor that meets the following orthogonality condition (T. Li et al., 2016;Ziaei-Rad et al., 2023): With this split, we rewrite the first integrand so that only the positive part of the strain energy contributes to fracturing: Equivalently, we can write where ℂ± = 2 2 ( 1 2 ℂ ∶ ± ∶ ± ) . (C9)

C3. Numerical Implementation
We control the injected volume in a quasi-static setting (t = t 1 , t 2 , …, t n ) as q = q 1 , q 2 , …, q n . In this simplified setting or the so-called toughness dominated regime (Detournay, 2016), we can simply enforce mass conservation at time step i as where V i is the fracture volume at time step i and is given as (Bourdin et al., 2012) With this mass conservation constraint, our minimization problem Equation C4 now reads where  is the kinematically admissible displacement set: and H 1 (Ω) is a Sobolev space of order 1 on Ω. The kinematically admissible set of v requires the irreversible condition (Bourdin, 2007;Burke et al., 2013) and v ir is the irreversible threshold ∈ [0, 1] (e.g., 0.05).
The current implementation guarantees to find local minima by considering the first-order optimality condition where  ′ ( , )(̃,̃) is the Gâteaux derivative of  at (u, v) in the direction of (̃,̃) and is given by: where ℂ( ) = 2 ℂ+ + ℂ− .
The minimization is performed via an alternate minimization scheme (Bourdin et al., 2000), taking advantage of the convexity of  with respect to u and v.  is first minimized with respect to u while fixing v and then minimized with respect to v while fixing u with the irreversible condition. As the solution for v requires an inequality constraint, 0 ≤ v ≤ η, we employed a variational inequality solver from the PETSc library (Balay, Abhyankar, Adams, Benson, et al., 2020; Balay, Abhyankar, Adams, Brown, 2020).

C4. Interface Model
Interfaces such as a grain boundary occupy negligible physical space, and thus it is challenging to account for their presence in numerical simulations. We follow an approach proposed by Hansen-Dörr et al. (2019 and Yoshioka et al. (2021) in which an interface is numerically diffused over a certain length b ( Figure C1). Within the diffused zone, the material posses an effective interface fracture toughness and outside of the diffused zone it posses the bulk fracture toughness. Toughness profiles in the sharp and diffused interface representations are presented in Figure C2.
We can estimate the effective fracture toughness ̃i nt c by equating the energies from the sharp and diffused representations. In the sharp representation, the fracture surface energy is int c ∫ Γ dΓ . This is diffused using the effective fracture toughness ̃i nt c :  where ξ is the shortest distance from the interface and ( ∇ ) ∶= 3 8 ( 1 − + |∇ | 2 ) . (C18) In earlier work (Yoshioka et al., 2021), Equation C17 was solved analytically in 1D and verified numerically in 2D and 3D. In this study, we computed G c profiles from the interface toughness int c and the toughness of bulk rock, bulk c .