Variational Phase‐Field Modeling of Hydraulic Fracture Interaction With Natural Fractures and Application to Enhanced Geothermal Systems

In every tight formation reservoir, natural fractures play an important role for mass and energy transport and stress distribution. Enhanced Geothermal Systems (EGS) make no exception, and stimulation aims at increasing the reservoir permeability to enhance fluid circulation and heat transport. EGS development relies upon the complex task of predicting accurate hydraulic fracture propagation pathway by taking into account reservoir heterogeneities and natural or preexisting fractures. In this contribution, we employ the variational phase‐field method, which handles hydraulic fracture initiation, propagation, and interaction with natural fractures and is tested under varying conditions of rock mechanical properties and natural fractures distributions. We run bidimensional finite element simulations employing the open‐source software OpenGeoSys and apply the model to simulate realistic stimulation scenarios, each one built from field data and considering complex natural fracture geometries in the order of a thousand of fractures. Key mechanical properties are derived from laboratory measurements on samples obtained in the field. Simulations results confirm the fundamental role played by natural fractures in stimulation's predictions, which is essential for developing successful EGS projects.


Introduction
Interest in predicting hydraulic fracture propagation is picking up since the Enhanced Geothermal System (EGS) concept could become a competitive solution as a sustainable and essentially carbon-free energy resource. In EGS, the reservoir is stimulated by injecting pressurized fluids in reservoir rock formations with the aim of enhancing permeability. Early application of permeability enhancement were performed in oil and gas reservoirs (Economides & Nolte, 1989), while nowadays the principles of hydraulic fracture mechanics are applied to a broad range of problems, such as nuclear waste disposal (Zoback et al., 2003), carbon-capture storage (Fu et al., 2017), glacier dynamics (Tsai & Rice, 2010), earthquake nucleation (Garagash & Germanovich, 2012), and geothermal systems (Fox et al., 2013;Legarth & Saadat, 2005;McClure & Horne, 2014). Hydraulic fracture propagation is intrinsically a multiscale problem (Garagash et al., 2011), with a wide range of scales of time and length controlling the fluid-driven crack propagation (Detournay, 2016). Under simplified assumptions of problem geometry and physical behavior, analytical solutions (Adachi & Detournay, 2002;Detournay, 2016;Garagash & Detournay, 2000;Savitski & Detournay, 2002) give good predictions of multiscale asymptotic behavior, which has been confirmed by laboratory experiments in highly controlled environments . Numerical methods are more computationally costly but can also overcome the simplifications typical of analytical solutions (Lecampion et al., 2018), such as planar cracks and homogeneous material properties (Bunger et al., 2013). Most numerical methods are based on Linear Elastic Fracture Mechanics (LEFMs) (Adachi et al., 2007), and the problem of hydraulic fracture propagation has been addressed either by (i) assuming planar and single mode crack propagation or (ii) accounting for nonplanar propagating cracks that interact with preexisting natural fractures (Weng, 2014).
The first approach assumes the crack as a planar object splitting the material in two parts with a displacement discontinuity that evolves over time: The dimensions of the hydraulic fracture (such as its length, height, and aperture) change as the fluid is injected. Models of three-dimensional (3-D) bi-wing planar fractures rely upon the known fracture models from Perkins, Kern, and Nordgren (PKN) (Perkins & Kern, 1961), Geertsma and de Klerk (GDK) (Geertsma & De Klerk, 1969), and more generalized three-dimensional models (Nordgren, 1972). The crack propagation criterion is based on the energy release rate (Griffith, 1921) and propagation occurs if the stress intensity factor reaches the critical value (Adachi et al., 2007). Viscous dissipation of fluid is an additional source of energy dissipation in hydraulic fracturing where the LEFM solution is coupled with Poiseuille's flow in the fracture and Carter's equation for leak-off from the fracture to the formation (Detournay & Cheng, 1993). The crack propagates along trajectories in a parametric space whose asymptotic regimes are characterized by a prevailing mechanism among leak-off, toughness, storage, and viscosity (Detournay, 2016). Rock's stiffness, strength, and permeability, fluid's viscosity, and injection rate control the trajectories of the parametric space. Although single Mode-I planar crack models give good estimates of the fracture dimensions whenever the basic assumptions hold valid, they fall short whenever heterogeneities cannot be neglected-a typical occurrence in geothermal reservoirs.
Models of fracture interaction (Jeffrey et al., 1994(Jeffrey et al., , 2009McClure et al., 2015;Renshaw & Pollard, 1995;Warpinski & Teufel, 1987;Weng, 2014) have to account for hydraulic fracture arrest, cross, or branch at the intersection with a natural fracture (Figure 1). Although Yew and Weng (2015) report the Unconventional Fracture Model (UFM) by Weng et al. (2011) as one of the first models of hydraulic fracture propagation that accounts for fluid flow and complex network of natural fractures, several problems regarding the computational mechanics of hydraulic fracture remain unsolved. Two main approaches have so far emerged: (i) the first one employs the Displacement Discontinuity Method (DDM), such as the Unconventional Fracture Model (UFM) or Crack Tip Open Displacement (CTOD), and (ii) the second one uses Finite Elements or Finite Volumes Methods (FEM or FVM), where natural fractures are either smeared using an implicit approach (nonlocal damage or phase field) or embedded into Cohesive Zone Models. The DDM is computationally inexpensive, as it requires discretization of the boundaries only, but cannot handle reservoir heterogeneities. The FEM with explicit embedded discontinuities faces two main drawbacks: (i) it requires fine crack-tip discretization to preserve accuracy, hampering its applicability to real case scenarios where the hydraulic fracture is expected to propagate for several hundreds of meters, and (ii) it suffers from element-distortion issues that generate inaccuracies in crack opening calculations and induce numerical instabilities. The eXtended Finite Element Method (Belytschko & Black, 1999;Belytschko et al., 2001;Gupta & Duarte, 2016;Moës et al., 1999;Wang, 2019;Yazid et al., 2009) overcomes the classical finite elements limitations of resolving field discontinuities by use of enriched shape function, although it is computationally expensive, can hardly handle hydraulic fracture-natural fractures interaction and can on occasions be numerically unstable. The phase-field method of fracture is a valid and promising alternative. Given its success in modeling propagation of brittle fracture, its development has been extended to ductile (Alessi et al., 2017;Ambati et al., 2015;Kuhn et al., 2016;Miehe, Mauthe, & Teichtmeister, 2015), fatigue (Alessi et al., 2018;Carrara et al., 2019;Seiler et al., 2018), and dynamic fractures (Borden et al., 2012;Bourdin et al., 2011;Fischer & Marigo, 2019;Hofacker & Miehe, 2012;Li et al., 2016;Schlüter et al., 2014). The variational phase field (V-pf) is a generalized Griffith criterion (Francfort & Marigo, 1998) numerically implemented using a phase-field variable, which smears the sharp interface fracture with a smooth transition function (Bourdin et al., 2000). The phase-field variable describes the transition from intact to fully damaged state of the material over a specific length scale. Seminal works of the application of the V-pf approach to hydraulic fracture include Bourdin et al. (2012) and Chukwudozie et al. (2013) while following studies addressed problems related to poroelasticity (Mikelić et al., 2015;Miehe, Hofacker, et al., 2015;Santillán et al., 2017;Wilson & Landis, 2016;Wheeler et al., 2014), fracture width computation (Lee et al., 2017;Xia et al., 2017), coupling with the theory or porous media (Ehlers & Luo, 2017;Heider & Markert, 2017), pressure-dependent failure mechanisms (Choo & Sun, 2018), mass conservation (Chukwudozie et al., 2019), and in situ stresses (Shiozawa et al., 2019). The smeared representation can handle complex fracture topology where natural fractures can be represented within nonconforming discretizations, without a priori assumptions on their geometry or restriction on hydraulic fracture growth trajectories (Yoshioka & Bourdin, 2016).
In this study, we solve the toughness-dominated hydraulic fracturing problem where the pressure drop within the fracture is negligible with a V-pf approach formulated with the constitutive model known as no tension or masonry model (Del Piero, 1989;Freddi & Royer-Carfagni, 2011). The main goal is to study the hydraulic fracture interaction with reservoir heterogeneities in the form of preexisting natural fractures with efficient computational V-pf models. We apply the V-pf method to a real case study of a potential EGS system, that is, the Acoculco geothermal field located in Puebla, Mexico. Two exploration wells were drilled within the geothermal field and, through log data analysis, a high-temperature (T ∼ 300°C) and low permeable (k = 1 × 10 −18 m 2 ) reservoir was identified at ∼2 km depth. Natural fractures are modeled as complex Discrete Fracture Networks (DFNs) calculated from outcrop field measurements and mechanical rock properties are derived from laboratory testing on samples collected in the field.
This article is structured as follows: In section 2, we introduce the governing equations of the V-pf model, their implementation in the open-source software OpenGeoSys (Kolditz et al., 2012), the experimental program, and the stochastic method to build DFNs. In section 3 we introduce applications of simple computational scenarios and geothermal reservoir stimulation. Section 4 presents the results of the simulations and contains a wider discussion of our results as well as broader implications of our main findings. Finally, we draw general conclusions of the study in section 6.

Variational Approach to Fracture
LEFM is based on Griffith's observation that the energy dissipation by a propagating crack equals the mechanical energy decay (Griffith, 1921). Thus, the criterion for fracture propagation is given as where G c is the critical surface energy release rate and G is the mechanical energy restitution rate. The energy restitution rate is defined as the derivative of the elastic energy P with respect to a crack increment length a, which is often derived using the concept of stress intensity factor (Irwin, 1957). Griffith criterion reads then as which was revisited by Francfort and Marigo (1998) noticing that it can be described in terms of critical values of the following total energy along a prescribed path as Griffith's criterion was generalized in the variational framework by considering a total energy with any crack set Γ as opposed to a prescribed path a as F |{z} Total energy such that the fracture propagation criterion is obtained by solving for the set of mechanical energy P and the crack geometry Γ that minimizes Equation 4. For a prescribed crack path (e.g., a), the approach converges to Griffith's criterion which can be viewed as a local energy minimum problem.

Governing Equations
The linear elastic constitutive equation of a brittle-elastic porous medium occupying a domain Ω can be expressed as (Biot, 1941) where C is the fourth-order linear elastic tangent operator, α is the Biot's coefficient, p p is the pore pressure, I is the identity matrix, and ε is the linearized strain tensor defined as the symmetric part of the dis- Also, consider crack set Γ filled with a fluid at pressure p f and let ∂Ω N be a portion of its boundary and ∂Ω D : = ∂Ω∖∂Ω N the remaining part, static equilibrium, and continuity of stress at the interfaces mandates that where f ! denotes an external body force and τ is a traction force. Multiplying (7) by a test function δu∈H 1 (Ω∖Γ) that vanishes on ∂Ω D and using Green's theorem, we obtain where N = 2 and N = 3 for 2-D and 3-D, respectively, and κ denotes the material's bulk modulus. We recall that given p p , p f , and Γ, Equation 9 is the unique solution of the minimization among all kinematically admissible displacement of is the poroelastic strain energy density (Yoshioka & Bourdin, 2016).

Phase-Field Approximation
The numerical implementation of the minimization of Equation 4 involves discontinuous deformation across unknown discontinuity surfaces (the cracks), Γ, which pose significant challenges in terms of numerical implementation. Instead, Equation 4 is regularized using the phase-field approach (Bourdin et al., 2000(Bourdin et al., , 2008. Introducing a scalar phase-field variable, v:Ω↦[0,1], which represents a state of the material from intact material (v = 1.0) to fully broken (v = 0.0) with a continuous function and a regularization parameter with the dimension of a length, ℓ s >0, which controls the transition length of the phase-field variable, Equation 4 can be approximated as (Bourdin et al., 2012) where c n is a normalization parameter defined as c n : ¼ ∫ 1 0 ð1−sÞ n=2 ds. Case n = 1 is often referred as AT 1 c n ¼ 3=2 ð Þ and Case n = 2 as AT 2 c n ¼ 1=2 ð Þ (Tanné et al., 2018). It can then be shown (Ambrosio & Tortorelli, 1990, 1992Braides, 1998) that as ℓ s approaches 0, the minimizers of Equation 11 converge to that of Equation 4 in the sense that the phase-field function v takes value 1 far from the crack Γ and transitions to 0 in a region of thickness of order ℓ s along each crack faces of Γ.
We can observe that in Equation 11, the evolution of the phase-field (v) is driven by the strain energy (W) regardless of the deformation direction, which leads to equal strength in tension and compression-a wrong approximation for granular material such as rock. To overcome the limitation, the strain energy can be decomposed into the positive (extension) and negative (shortening) parts Among the several approaches for the strain decomposition (Amor et al., 2009;Freddi & Royer-Carfagni, 2011;Miehe et al., 2010), we employ the so-called masonry model (Freddi & Royer-Carfagni, 2011), in which the material will not withstand tensile stresses.
Open natural fractures can be represented by assigning either the phase-field variable v = 0.0 (Ni et al., 2020;Yoshioka & Bourdin, 2016) or the fracture toughness G int c ¼ 0:0 (v = 0.0 is achieved, where G int c is assigned with 0.0 at the first iteration). In this study we represent discontinuous interfaces by a diffused variable of the phase-field type ( e G int c ) whose fracture toughness or cohesive strength (G c ) is different from the surrounding (Figure 2). Altering

Journal of Geophysical Research: Solid Earth
and ξ is the distance from the crack (v = 0). We built the FEM model containing natural fractures by assigning the equivalent fracture toughness computed in Equation 13 to the region within distance b from the fractures.

Numerical Implementation
We neglect leak-off to the rock mass because the permeability of the rock mass is sufficiently low. The pore pressure p p can be considered invariant and set as p p = 0, and p ′ f ¼ p f − p p in the governing equations. We adopt the notation p ′ f ¼ p . Considering hydraulic fracturing in the toughness dominated regime (Detournay, 2016), the pressure loss within the crack is negligible and p is spatially constant. Equation 11 is solved by the alternate minimization with respect to the displacement u and the phase-field v with a constraint of prescribed time-evolving fluid volume which must be equal to the crack volume, that is, . The minimization problem can be stated as with the constrain The first variation of the energy functional with respect to u is where C ± is the tangent stiffness tensor The first variation of the energy functional with respect to v for AT 1 is and for AT 2 is Equation 7 is linear to p and if we let the displacement solution with p = 1.0 be equal to u 1 , the displacement solution p = p is obtained as u = pu 1 and the crack volume is At a given time step, a volume V inj is injected and Equation 16 yields the mass balance in the porous medium such that the corresponding p is where and the whole solution procedure is described in Algorithm 1.
Parallel simulations were run on the high-performance computing system JUWELS, maintained at the Jülich Supercomputing Centre. The total number of degree of freedom for the Acoculco case scenarios is 513,108 with 170,996 linear quadrilateral elements with a few triangular elements in the mesh resolution transition zone. Domain decomposition was done using METIS (Karypis & Kumar, 1998), and both linear and nonlinear solvers from PETSc Balay et al. (2019) were used. More specifically, the Newton-Raphson solver for the deformation problem and a Newton based variational inequality solver for the phase field, since the phase-field solution is bounded in [0,1] domain and constrained by the irreversibility. The simulations were distributed into 384 cores over eight nodes (Dual Intel Xeon Platinum 8168) with 2×24 cores. While the computation time differs depending on the nonlinearity of each problem, all the simulations shown in the subsequent section were completed within ∼20 hr.

Sensitivity Analyses
We perform a sensitivity analysis to highlight the complex interactions between fluid-driven propagating fractures and existing ones. We analyze simplified models assuming a perfectly homogeneous brittle linear elastic material containing one or two preexisting natural fractures. We evaluate the impact on fracture initiation and propagation given by stiffness (elastic parameters) and strength (fracture toughness) of the bulk rock, existing fractures, state of stress, and orientation angle of the natural fractures. The base case parameters of the sensitivity analysis are in the range of the studied rocks of the Acoculco reservoir, that is, G c = 100 Pa m, E = 25 GPa, and ν = 0.2.
As implemented in the V-pf, the interaction with the preexisting natural fractures is partly controlled by the critical energy release rate of the natural fracture˜G c int . We compare results for G c = 1 with different values of˜G c int spanning 4 orders of magnitude, that is,˜G c int ¼ 0:01; 0:05; 0:1; 10. We analyze the impact of the far-field deviatoric stress by increasing the S Hmax from 21 to 60 MPa with a fixed S hmin = 20 MPa. All case scenarios are presented in Table 1. We finally analyze the influence of a natural fracture oriented 45°from the direction of hydraulic fracture propagation to study the effect of the incidence angle between the direction of propagation and the orientation of a natural fracture.

The Acoculco Geothermal Reservoir
In this manuscript, we analyze the potential permeability enhancement in a fractured reservoir by modeling the fracture growth from a well-bore injection. We apply the developed numerical methodology described in

10.1029/2020JB019856
Journal of Geophysical Research: Solid Earth section 2 on the Acoculco case scenario, considering the local geological features and the rock properties gathered from field campaigns and laboratory measurements. Here we report a brief synthesis of the experimental program, field campaign, DFN construction algorithm, and results, while further details can be found in the original works (Lepillier et al., 2019(Lepillier et al., , 2020. The Acoculco geothermal field, located in Mexico, hosts two vertical exploration wells (named EAC1 and EAC2) drilled at ∼500 m apart horizontally, both reaching a total depth of ∼2 km (Canet et al., 2015;López-Hernández et al., 2009;Weydt et al., 2018). On the one hand, Acoculco is considered a tight reservoir because the rock formations are little permeable (10 −18 m 2 ) and the fractures are scarcely connected (Lepillier et al., 2020); on the other hand, it is a suitable candidate for EGS development because of its high geothermal energy potential given that the geothermal gradient is above average (∼150°C km −1 ). The stratigraphy encountered during drilling is simplified into three lithological units: limestones, marbles, and skarns. Stiffness (E and ν), strength, and critical energy release rate (G c ) of the three lithologies were measured in the rock physics laboratory. Young's modulus and Poisson's ratio were determined by Unconfined Compression Strength (UCS) tests (UCS-20 experiments). Brazilian disc (BD-80 experiments) and Chevron Bend tests (CB-12 experiments) were employed to determine the fracture toughness K Ic of the rock formations, which was later employed to derive the critical energy release rate (G c ). Fracture toughness was determined from the two sets of experiments: (i) from BD tests, it was done following the method proposed by Guo et al. (1993), and (ii) for CB tests, following the method suggested by Franklin et al. (1988). All material parameters employed in the model are summarized in Table 2.
The general DFN is derived from scanline measurements from multiple outcrops analogues of the Acoculco geothermal system (Lepillier et al., 2020) that are later extrapolated using the multiple point statistic method . The method yields three separate DFNs, that is, one per lithology. Each one of the three DFNs is a bidimensional georeferenced section of 600×600 m 2 (Figure 3). Some further processing is necessary to build the FEM models. In the first step, we extracted from each DFN a smaller subdomain of 100×100 m 2 (Figure 3). Each extraction has a specific fracture distribution: to analyze the impact of stimulating one or another specific section of the domain. In the second step we extracted an additional three subdomain from each of the DFNs. The three sub-DFNs, one for each DFN, are then rotated in the third step to align the maximum horizontal stress S Hmax with the x axis and further downscaled in the fourth step to fit the a-dimensional V-pf formulation.
The in situ stress state is believed to be of the normal faulting type and the orientation of the stress tensor is taken from the World Stress Map (Heidbach et al., 2016;Lepillier et al., 2019). Based on this assumption, having S v >S Hmax >S hmin , we defined certain values for S Hmax and S hmin . In normal faulting regime, the hydraulic fracture propagates along the vertical plane oriented perpendicular to S hmin . Because of this, we assume two-dimensional plane-strain conditions were we assign only S Hmax and S hmin . Figure 4 shows the influence of strength (G c ) and stiffness (E and ν) on the internal fracture overpressure and length evolution during hydraulic fracture propagation at constant fluid injection rate. The critical energy release rate G c is the dominant parameter controlling the hydraulic fracture response (Figure 4a): G c represents the resistance to fracture propagation and hence is proportional to the maximum overpressure reached and inversely proportional to the rate of crack length growth during injection. The stiffness parameters play a smaller role on the problem evolution, and, while the influence of Poisson's ratio seems to be negligible over the selected range (Figure 4c), an increase in Young's modulus entails an increase in fracture propagation resistance (Figure 4b). Maximum overpressure is proportional to Young's modulus and inversely proportional to the injected volume at propagation onset.

Sensitivity Analysis
The delay in crack propagation onset is a consequence of lower stiffness: the more the rock is compliant, the larger the volume of fluid needs to be injected before the crack internal pressure reaches the propagation condition and the energy release rate equals its critical value. Globally, it can be interpreted as a higher system compressibility, where more compliant systems require higher volume of injected fluids. Figure 5 shows results of the sensitivity analysis of the interaction between a fluid-driven fracture (the phase field) and two natural fractures at equivalent time steps.
When the natural fracture strength is considerably weaker, the fracture turns at and propagates along the vertical natural fractures (˜G c int ¼ 0:01 and˜G c int ¼ 0:05). These cases exhibit asymmetrical growth when the tips hit the natural fractures by extending only one of the tips rather than the both. Energetically, propagation of one of the tips is as easy as that of the other while propagating the both simultaneously is more expensive. In other words, the energy required to propagate the right tip equals the energy required for the left tip propagation. As a consequence of the energy minimization in such a situation, propagation of one of the tips (not the both) will be eventually chosen as "minimum" energy state since numerically one of the energies is always smaller than the other by some truncation. However, once either one is chosen, this  side will be always the easier direction (will require less energy) in the subsequent propagation. For˜G c int ¼ 0:1, the natural fractures do not have low enough strength and are crossed by the hydraulic fracture without branching ( Figure 5). For˜G c int > G c , the natural fracture act as a barrier to the hydraulic fracture. After the crack hits the natural fracture, it propagates in a path avoiding the natural fracture. In this case, the natural fracture acts as a barrier, shielding the hydraulic fracture propagation. Note that branching in general is energetically more expensive (less favored) as it is avoided in the˜G c int ¼ 10 case but does happen when the surface energy of the natural fractures are so low that crack propagation along them becomes more attractive.
At increasing values of differential stress (Figure 6), and for fixed˜G c int ¼ 0:01, the branching observed at low deviatoric stress disappears for S Hmax ≥ 30 MPa. The critical stress intensity factor at the tip of the natural fracture is proportional to the horizontal stress and propagating a fracture parallel to S Hmax through the bulk rock requires less energy than propagating it through the vertical natural fracture. Therefore, with higher deviatoric stress, considering a natural fracture oriented 90°will not change the propagation direction as the stress dictates the propagation direction.
A 45°oriented fracture has an orientation, which is closer to the critical one for the given state of stress; hence, it influences the propagation and interaction regime differently than vertical natural fracture ( Figure 7). With only one natural fracture present, the problem is intrinsically asymmetric. At˜G c int ¼ 0: 01, the hydraulic fracture first interacts with the natural fracture and later propagates in the direction of S Hmax (Figure 7a) and at˜G c int ¼ 0:1, the hydraulic fracture propagation is still attracted by the inclined natural fracture but not as much as the case with˜G c int ¼ 0:1. For high values of the natural fractures' critical energy release rate, that is, for˜G c int ¼ 10, even though the natural fracture is more favorably oriented, it becomes once again a barrier to fracture propagation (Figure 7a). For˜G c int ¼ 0:01 with varying horizontal stresses S Hmax , the hydraulic fracture propagation along the natural fracture is progressively hindered with increasing S Hmax (Figure 7b). At S Hmax = 40 MPa, the hydraulic fracture shows a small offset at the natural fracture's crossing point while the hydraulic fracture becomes agnostic to the natural fracture with S Hmax = 60 MPa.

Stimulation of the Acoculco Geothermal Reservoir
The natural fractures of the Acoculco reservoir are assumed to be cemented and hydraulically closed prior to stimulation. Following our sensitivity studies, we assume the fracture toughness of the natural fracture is 10% of the surrounding (G int c ¼ 0:1G bulk c ). The simulated domain is discretized using a fixed mesh where the element size is 25 cm. Figures 8-10 show the results of the stimulation scenarios in the Acoculco geothermal reservoir for the different lithologies and for different DFNs. On the left of all figures is plotted the fracture pressure and length with injected volume, while on the right is shown the contour map of the phase-field along with the distribution of natural fractures. For all the cases, the propagation pressure decreases with injected volume as the crack length increases. The pressures started declining rapidly from the onset of the injection/stimulation. This is because the simulations were initiated with a borehole without setting a priori (initial) fracture The stress field is oriented such that S Hmax is aligned along the horizontal direction. The red dots represents the well-bore and initial fracture position (and initial phase-field implementation).

10.1029/2020JB019856
Journal of Geophysical Research: Solid Earth lengths, as often done in practice, which led to the high breakdown pressures. Such high pressure responses may not be observed in fields because (1) the borehole intersects with preexisting fractures or defects or (2) the borehole is completed with perforations or well production packers. However, if fracture is initiated in a intact rock, this level of high pressure should be expected. The fracture length increment with time shows a burst-like behavior: Whenever the hydraulic fracture interacts with a natural fracture, the pressure drops as a consequence of the increase in available fluid storage capacity given by the crack sudden propagation over a finite distance within the natural fracture.
Considering all lithologies, the final fracture length ranges between ∼75 and ∼95 m and a the lowest propagation pressure is observed for the marble stimulation cases (Figure 9), while the highest propagation pressure is observed for the stimulation into the limestone formation ( Figure 8)-a result in agreement with the sensitivity analysis. Figure 11 shows a polar representation of the hydraulic fracture deviation from the x direction during propagation. The limestone simulations show the larger range of fracture lengths spanning from ∼75 to ∼95 m while the marble's one have the smallest range, spanning from ∼78 to ∼79 m. The angular deviation ranges  such that S Hmax is aligned along the horizontal direction. The red dots represents the well-bore and initial fracture position (and initial phase-field implementation).
in an interval of 20°above and below the reference axis given by S Hmax direction. Maximum deviations are observed in marble and skarn simulations, reaching 30°in both simulations, while the deviation angle for the Limestone simulations is contained in a 20°interval.
The asymmetrical propagation of hydraulic fracture from the borehole is a consequence of the intersection angle between natural fracture and the approaching hydraulic fracture. Assuming θ as the angle between a natural fracture and the S Hmax axis, we observed that (i) low-θ natural fracture act as pathways for the hydraulic fracture, which propagates faster along natural fractures; (ii) high-θ natural fracture (∼90°) are bypassed by the hydraulic fracture and no interaction takes place. Intermediate values of θ offer a pathway for hydraulic fracture to propagate along a certain distance, until the pressure buildup is sufficiently high to allow further propagation within in the matrix.

Discussion
The V-pf method presented here is an implicit smeared approach, which represents the fracture with a smoothly transitioning function that spans from intact to fully damaged state of the material. Natural fractures are represented in a non-conforming mesh with the reduced fracture toughness by enforcing energetic Figure 8. Hydraulic fracture models using V-pf with the sub-DFN of the limestone reservoir. The matrix material domain Ω is represented in gray, the natural fracture Γ are discretized in black. Lm01 is composed with 1,483 natural fractures, Lm02 is composed with 709 natural fractures, and Lm03 is composed with 327 natural fractures. The stress field is oriented such that S Hmax is aligned along the horizontal direction. The red dots represent the well-bore and initial fracture position (and initial phase-field implementation).

10.1029/2020JB019856
Journal of Geophysical Research: Solid Earth equivalence, which is one of the advantages of the method since it allows exploring multiple DFNs scenarios with a single discretization. As presented in this study, the ability of the V-pf is to handle complex fracture topologies with unified criteria-energy minimization-that seeks for an admissible displacement and a set of fracture geometry that minimizes the total energy without a need for ad hoc criteria for branching or turning. The model exhibits asymmetric crack growth under some circumstances: the phenomenon is intrinsic to the variational approach, where the energy minimization leads to the occurrence of asymmetric solutions whenever the total energy of the system is smaller than its symmetric counterpart (Tanne, 2017). This instability has been theoretically pointed out by Gao and Rice (1987), and it is observed in experiments for a penny-shaped crack propagation in toughness dominated regime (Bunger et al., , 2013. The interaction behavior between hydraulic fracture and natural fractures depends (i) on the combination of the critically energy release rate ratio between natural fractures and bulk rock (˜G c int =G c ), (ii) on the natural fractures orientation relative to the stress field, and (iii) on the magnitude of principal stress components. Natural fractures can either favor or hamper the propagation of a hydrofracture according to specific combinations of the input parameters. Natural fractures attract hydraulic fractures for relatively low values of Figure 9. Hydraulic fracture models using V-pf with the sub-DFN of the marble reservoir. The matrix material domain Ω is represented in gray; the natural fractures Γ are discretized in black. Ma01 is composed with 295 natural fractures, Ma02 is composed with 215 natural fractures, and Ma03 is composed with 198 natural fractures. The stress field is oriented such that S Hmax is aligned along the horizontal direction. The red dots represent the well-bore and initial fracture position (and initial phase-field implementation).

10.1029/2020JB019856
Journal of Geophysical Research: Solid Earth critical energy release rate ratio, when they have orientations close to the critical ones and for relatively isotropic stress states. Natural fractures can be an obstacle to hydraulic fracture growth whenever the fracture resistance becomes higher than the one of the intact rock. Although counterintuitive, the presence of higher strength discontinuities is a relatively frequent occurrence in deep geothermal systems: The environmental conditions could enhance geochemical reactions of dissolution and precipitation (Singurindy & Berkowitz, 2005;Watanabe et al., 2020), such as silica precipitation (Lu et al., 2018;Scott & Driesner, 2018), and the existence of active volcanism could favor the presence of magmatic intrusions even at shallow depth (Elders et al., 2014), which, if old and cold enough, could represent higher strength and stiffness bodies.
In our analyses, we have assumed a low permeability that is typical of poorly fractured crystalline rock, a hypothesis that entails no leak-off between the fracture and the porous rock. Such an assumption is equivalent to an undrained response where the change in effective stress within the porous rock is null during injection. Although the fracture toughness (critical energy release rate) is more predominant in controlling propagation conditions when compared to stiffness, Young's modulus of the rock also plays a role. In particular, a more compliant rock requires higher injected volumes but overall generates lower overpressure. On Figure 10. Hydraulic fracture models using V-pf with the sub-DFN of the skarn reservoir. The matrix material domain Ω is represented in gray; the natural fractures Γ are discretized in black. Sk01 is composed with 706 natural fractures, Sk02 is composed with 495 natural fractures, and Sk03 is composed with 375 natural fractures. The stress field is oriented such that S Hmax is aligned along the horizontal direction. The red dots represent the well-bore and initial fracture position (and initial phase-field implementation).

10.1029/2020JB019856
Journal of Geophysical Research: Solid Earth the contrary, stiff rocks generate higher overpressure for a lower injected volume. Because of the high fracture strength, high stiffness, and low permeability of basement crystalline rocks, during stimulation of a deep geothermal reservoir high overpressure can be achieved with relatively low values of injected volume (Ellsworth et al., 2019).
The V-pf simulations of the Acoculco reservoir highlighted a fluctuation in the pressure and crack length response in time, with intermittent crack advancement and burst-like behavior-a phenomenon observed during several hydraulic stimulations (Milanese et al., 2016). The V-pf implementation adopted is numerically stable and provides continuous pressure-volume response in absence of viscous flow dissipation. The intermittent advancements are a direct consequence of the interaction between existing fractures with lower crack resistance and the fluid-driven crack: Whenever the hydraulic fracture encounters a natural fracture, if the latter is favorably oriented, the hydraulic fracture encounters almost no resistance and propagates rapidly over a finite length. The pressure drop is associated with a stress release in the rock, which in combination with the crack length increment, can be associated with micro-earthquakes. Microseismicity has been widely observed during hydraulic fracturing (Davies et al., 2013;Lopez-Comino et al., 2017;Schultz et al., 2015) and our results suggest that, in crystalline reservoirs, the phenomenon is associated with hydraulic fractures propagating along preexisting natural fractures.
Results show that the marble formation in the Acoculco reservoir is the optimal one for a potential stimulation because the lowest values of propagation overpressure. The orientation of the natural fractures controls the propagation extent and direction independently of the lithology and the fracture topology dominates the hydraulic fracture response in all cases analyzed. In the present case study we have analyzed homogeneous Figure 11. A comparison of hydraulic fracture simulation V-pf models for (a) the limestones (green color); (b) the marbles (blue color); (c) the skarns (red color); hydraulic fracture lengths are given by the concentric dividers, and hydraulic fracture angles compared to S Hmax 's orientation is given by the radial dividers. rock matrix properties, although a more realistic approach should be based on representing fluctuation of the material properties within the rock matrix. Three-dimensional analyses would be an additional improvement of the current scenarios. Nonetheless, the additional complexity of a three-dimensional reservoir model should be justified by a sufficient knowledge of the reservoir's structure and its property-a current shortcoming for the Acoculco reservoir. Although a normal fault regime is likely at Acoculco reservoir and hydraulic fractures are expected to propagate mainly vertically, there are indications that a strike-slip faulting system could also be active (Liotta et al., 2020), making the full three-dimensional propagation topology rather complex and difficult to estimate a priori.
There is current uncertainty about the in situ state of stress at the Acoculco geothermal reservoir and different values of the stress components could yield a different output in terms of hydraulic fracture propagation.
Although the DFNs that are built from outcrop extrapolations are also a source of additional uncertainty, the small prominence of fractures in the DFNs seems to be in good agreement with the very low permeability that was observed during well logging: small and poorly connected fractures hamper fluid flow in the tight reservoir.
Stimulating a highly fractured zone of the Acoculco geothermal reservoir requires a lower stimulation pressure, therefore reducing the drilling costs. Additionally, according to the well temperature measurements, the marble and skarn formations are more likely to be targeted for stimulation because they are present at a higher depth, and therefore, they are at a higher temperature. The formation breakdown pressure is lower for the marble, which also has a lower density of natural fractures. Nonetheless, the natural fractures in the marble are longer and better connected when compared to the ones in the skarn, which are shorter but more frequent. A trade-off arises between the objective of stimulating the hotter formations of the reservoir on the one hand, and stimulating the formations that would yield a longer propagation of the hydraulic fracture on the other hand. The optimal solution would depend on the ultimate goals of the EGS development and a detailed cost-balance analysis is necessary to optimize the stimulation depth.

Conclusions
We have presented a method for modeling hydraulic fracture propagation and interaction with a network of natural fractures in a geothermal reservoir. The fracture simulations are based on a variational phase-field approach that proved high numerical stability. We have highlighted the main factors controlling the hydraulic fracture propagation and its interaction with natural fractures through sensitivity analyses on simplified models. We have applied the method to model a realistic EGS stimulation scenario of the geothermal reservoir of Acoculco, Mexico. The numerical model is built from field data and model parameters are derived from laboratory experiments.
Building a realistic DFN is an essential piece of the puzzle for numerical analyses of stimulation of complex reservoirs, which can lead to counterintuitive findings of the propagation mechanisms as opposed to simplified models of single-oriented crack families. Pressure fluctuations and burst-like crack propagation are intrinsically connected to the presence of the complex network of natural fractures.