Modeling Fluid Reinjection Into an Enhanced Geothermal System

The manuscript analyzes the stimulation for an Enhanced Geothermal System development in Acoculco, Mexico. Using an analytical penny‐shaped hydraulic fracture model covering different propagation regimes, we computed the final fracture length and width by varying fluid properties with temperature. Our analysis indicates that for the given scenario, the fluid viscosity plays a minor role and instead flow rate and time of the stimulation are the controlling variables. We computed the fracture reopening as a consequence of water reinjection in the second stage of the stimulation through numerical computations based on the enriched discontinuity method. The computation shows that a single isolated fracture will not provide sufficient permeability, as the continuous injection will quickly fill and pressurize the crack. We demonstrate that the fracture needs to be connected to a permeable network to avoid excessive pressurization and achieve a commercially exploitable reservoir for Enhanced Geothermal System.


Introduction
The Acoculco geothermal site, located in the Trans-Mexican Volcanic Belt, is a potential candidate for the development of an Enhanced Geothermal System (EGS). While the rock permeability is very low (∼10 −18 m 2 ), the temperature is sufficient for electricity generation (∼300°C) at a relatively shallow depth (∼2 km) (Lepillier et al., 2019;Peiffer et al., 2015). The heat transfer within the reservoir is considered purely conductive because the two exploration wells measured a linear geothermal gradient of 0.13 K m −1 and retrieved no fluids (Calcagno et al., 2018;Lorenzo-Pulido et al., 2010). An EGS requires advective heat exchange to harvest geothermal energy and water circulation can be achieved by permeability enhancement through propagation of new fractures and reopening of existing ones. Hydraulic stimulation in geothermal reservoirs has certain peculiarities: (i) proppants are almost never employed; (ii) zonal isolation is rarely applied because of either packer's limitations at higher temperature and/or concerns about sealing off existing fractures (feed zones) by cementing; and (iii) stimulation is often carried out within the basement rocks that pose a greater risk of inducing earthquakes.
Hydraulic stimulation is a key component in EGS projects. However, it still remains a rare practice in conventional geothermal reservoirs and sometimes it is serendipitous achievement as a by-product of long-term cold water reinjection (Grant et al., 2013;Yoshioka, Pasikki, & Stimac, 2019). In EGS stimulation, pure shearing (Mode-II) of preexisting fractures is unlikely. The hypothesis of mixed-mechanism stimulation is the propagation of new tensile fractures (Mode-I) interacting with preexisting ones (McClure & Horne, 2014). At Acoculco, we expect Mode-I fracture propagation as the most likely mechanism of stimulation because of a scarcely fractured formation indicated by the low permeability (Lepillier et al., 2020). EGS stimulation is normally performed without pumping proppants or any chemical additives and after the injection, fractures will close up to the minimum aperture that rock asperities can hold. The permeability of fractures will be further decreased by the depleted formation because of fluid extraction during EGS operation, which ultimately hinders fluid flow. Overall, successful heat extraction depends mostly upon the residual permeability that is achieved during stimulation.
In multistage stimulation in oil and gas fields, fracture initiation is enforced at the perforations made through cement and casing which forms complex fracture patterns (Warpinski & Wolhart, 2016) along the horizontal wells. Most geothermal wells, on the contrary, are drilled with some deviation and completed with preperforated liners, aiming to intersect as many preexisting fractures as possible. Zonal isolation is thus rarely employed, with the exception of a few EGS projects (Bendall et al., 2014;Petty et al., 2013;Zimmermann, 2010). Without isolation and/or selective stage control, the state of the stress or weak spots (preexisting fractures, faults, or bedding planes) determine the initiation and propagation of hydraulic fractures. In poorly fractured formations, such as Acoculco, a planar fracture whose normal direction coincides with the least stress is expected to grow.
To maximize the energy output, EGS are usually stimulated in the hotter and deeper part of the reservoir, at depth ranging from 3 to 5 km in the stiff and low permeable crystalline basement. Such depth corresponds to an increase of seismogenic potential and to minimize the risks of inducing large magnitude earthquakes (Ellsworth et al., 2019;Terakawa et al., 2012), relatively low flow rate and limited injection pressure must be applied (Brown, 2009;Chabora et al., 2012;Garcia et al., 2016;Matsunaga et al., 2005;Schill et al., 2017).
The goal of this contribution is to study the stimulation effects at Acoculco. In particular, we aim at characterizing the hydromechanical response of water reinjection in a previously stimulated fracture. We use a combination of analytical and numerical methods: While the former are robust tools to study the involved mechanisms and give indications about preliminary design of stimulation methods (Detournay, 2016;Dontsov, 2016a;Garagash & Detournay, 2000), the latter allow us to explore complex scenarios and reduce the number of hypotheses and simplifications (Hofmann et al., 2016;Lecampion et al., 2018;McClure & Horne, 2014;Yoshioka & Bourdin, 2016). First, we study the propagation of a pure Mode-I fracture in the intact rock and analyze the influence of temperature and flow rate in hydraulic fracture growth. Second, we study the effects of reinjection in the newly created fracture in the reservoir with a finite element based approach with enriched discontinuities (Watanabe et al., 2012). We show that hydraulic stimulation in the low-permeable crystalline basement is likely to produce large fractures with radius ∼200 m, but if such fractures do not intersect preexisting permeable structures, the overall injectivity remains insufficient for successful geothermal exploitation. To overcome the problem, the most promising techniques would be zonal isolation combined with multistage stimulation.

Fluid-Driven Fracture Propagation
The problem of fluid-driven fracture propagation considers a penny-shaped fracture that propagates quasi-statically in a permeable, linear elastic medium and is driven by a pressurized incompressible Newtonian fluid (Detournay, 2016;Garagash & Detournay, 2000;Savitski & Detournay, 2002). The far field stress is σ 0 in Mode-I (Savitski & Detournay, 2002), the well radius is negligible compared to the radius of the fracture R(t) and the porous medium behaves accordingly to linear elastic fracture mechanics (LEFM). The fluid is injected at the center of the fracture at a constant volumetric flow rate Q 0 . The global solution is controlled by the near-tip asymptotic behavior, which involves wide temporal and spatial scale transitions and corresponds to a trajectory in the phase space formed by four vertices (Detournay, 2016): (i) the viscous-storage regime (no leak-off in the matrix and the energy dissipation in the total system is dominated by the fluid's viscous forces); (ii) the toughness-storage regime (no leak-off in the matrix and the dissipation is dominated by the creation of new fracture surface); (iii) the viscous-leak-off regime (presence of leak-off in the matrix and the dissipation is dominated by the fluid's viscous forces); and (iv) the toughness-leak-off regime (presence of leak-off in the matrix and the dissipation is dominated by the creation of the new fracture surface). Following Dontsov (2016a), the lubrication equation describing fluid flow in the fracture writes where w(r) is the fracture opening, t is time, μ is the fluid's dynamic viscosity, μ ′ ¼ 12μ, p w is the fluid pressure within the fracture, C l is the leak-off coefficient, and t 0 (r) is the time at which the fracture is located at point r. The elasticity equation of the porous solid controls the crack opening based on the fluid pressure is where E is the Young's modulus, ν is the Poisson's ratio, where the functions Υ(a) and X(a) are the complete elliptic integrals of the first and second kind, respectively. The fracture propagates if the Mode-I stress intensity factor equals the critical value K ¼ K Ic and the critical stress intensity factor determines the asymptotic limit of the fracture opening at the tip region as . Dontsov (2016a) provided a closed-form approximate solution that covers the whole phase space of propagating regimes. We employ a Matlab © script provided in Dontsov (2016b) to solve the fluid-driven propagation of an initial fracture in the reservoir.

Lower Interface Elements for Reinjection Modeling
The hydraulic fracture theory provides an estimate of the fracture length and width at the end of the stimulation treatment. Once the injection stops, the newly created fracture starts to close unless proppants are used during the stimulation. We analyze fracture reopening by continuous water reinjection and assess the newly created EGS. We perform hydromechanical finite element analyses using the open source platform OpenGeoSys (www.opengeosys.org). The fracture created in the first stage (stimulation) is modeled as a local enrichment at the element boundaries applying a Lower-dimensional Interface Element (LIE) (Watanabe et al., 2012;Yoshioka, Parisio, Naumov, et al., 2019).
The governing equations of the problem are derived from the linear momentum and mass balance equations. The linear momentum balance writes where σ ′ is the effective stress tensor in the porous matrix, α is the Biot's coefficient, p is the pore pressure, I is an identity tensor, ρ m ¼ ϕρ w þ 1 − ϕ ð Þρ s is the density of the porous medium as a function of porosity ϕ, water density ρ w and solid density ρ s , and g is the gravity acceleration vector. The elastic constitutive equation is σ ′ ¼ C: ϵ, where C is the fourth-order elastic stiffness tensor and the linearized elastic strain tensor is defined as with u the displacement vector. The fluid mass balance in the porous medium writes 10.1029/2020GL089886

Geophysical Research Letters
where β w is the fluid's compressibility, K s is the bulk modulus of the solid phase, k is the intrinsic permeability tensor, and Q m is a source/sink term. The relative displacement between the two faces of a fracture has a normal w n and a tangential component w t as w ¼ fw n ; w t g T which, in a local coordinate system, are where the superscripts + and − indicate the values of the displacement field components at each side of the fracture and j refers to either normal j = n or tangential j = t components. The relative displacement is a jump in the global displacement field u. The total stress acting on the fracture plane σ f is a function of the effective stress σ ′ f and the pore pressure in the fracture p w as where m ¼ f0; 1g T and σ′ f ¼ fσ′ f; t ; σ′ f ;n g T splits into normal and tangential component. The constitutive relation in the fracture writes with the fracture stiffness matrix with k t and k n the tangential and normal fracture stiffness, respectively. The mass balance equation in the fracture writes where the fracture aperture b ¼ m T w has been introduced along with the fluid fluxes from each side of the fracture into the matrix, that is, q + and q − . The fracture permeability is based on the lubrication flow in fracture and is a function of the aperture k f = b 2 /12.
We assumed two stimulation scenarios for the fluid temperature at downhole: (i) the fluid reaches a thermal equilibrium with the rock (hot fluid) and (ii) the fluid remains at wellhead temperature (cold fluid). The fluid properties for the hot injection are taken at temperature and pressure at 1,900 m depth, that is, p ¼ 16:7 MPa and T ¼ 298°C so that the fluid viscosity is μ ¼ 9:1 × 10 −5 Pa s, the specific density ρ w ¼ 793 kg m −3 , enthalpy h ¼ 1397:3 kJ kg −1 and the compressibility β w ¼ 2:46 −9 Pa −1 . The properties of the fluid for the cold injection scenario are computed for T ¼ 20°C and p ¼ 16:7 MPa. The two cases aim to investigate the effects of the fluid rheology on hydraulic fracture propagation: we are not accounting for the effects induced by temperature variations such as thermo-elastic strains or thermal expansion of water. The injection flow rate for the base case is Q 0 ¼ 0:1 m 3 s −1 and a stimulation time of T ¼ 20 min is considered for a total injection volume of 120 m 3 .

Fluid Reinjection After the Initial Stimulation
Given a fixed size fracture, we investigated the consequences of fluid reinjection in terms of fracture reopening and pressurization with the finite element method with lower interface element. We analyzed three 10.1029/2020GL089886

Geophysical Research Letters
PARISIO AND YOSHIOKA different scenarios with a bidimensional model containing a 180 m long radial crack: (i) a short-term (∼10 3 s) reinjection in an isolated crack; (ii) a medium-term (∼2 × 10 4 s) reinjection in an isolated crack; (iii) a medium-term reinjection (∼2 × 10 4 s) in a crack connected at its tip with a highly permeable fracture network (simulated as a constant pore pressure boundary equal to the reservoir initial pore pressure, i.e., free-flow conditions). To study the impact of the crack size, we simulated two additional scenarios of long-term (∼2 × 10 6 s) reinjection in an isolated crack of 180 and 930 m length. Due to the size effects, the fracture stiffness decreases with size and for fractures ∼100 m of length, the shear stiffness k s should be in the order of 10 −2 to 10 −1 GPa m −1 (Bandis et al., 1983). At σ f; n ′ > 1 MPa, k n =k s ¼ 2 ÷ 20 so that k n ¼ 0:02 ÷ 2 GPa m −1 . In this case, we assume k s ¼ 0:01 GPa m −1 and k n ¼ 0:1 GPa m −1 for the 180 m long fracture k s ¼ 0:06 GPa m −1 and k n ¼ 0:6 GPa m −1 for the 930 m long fracture. Figure 1 shows evolution of width (a) and permeability (b) for the two scenarios (hot and cold fluids). Injecting hot fluid generated a smaller fracture width at the injection point compared to the cold fluid case. The fracture radius is slightly larger in early times and became almost identical after 20 min of injection (Figure 1a). The difference between the two cases reduces with time, since the hydraulic fracturing conditions transition from a viscous to a toughness dominated regime, in which the viscous forces are a secondary control with respect to the fracture resistance. The reduced opening in the lower viscosity (hot fluid) case is related to the higher matrix leak-off (more fluid is transferred from the fracture to the porous matrix), since Carter's leak-off coefficient is inversely proportional to the viscosity. The fracture permeability increases with time ( Figure 1b) and the difference between the two scenarios is marginal. The increase of permeability at the injection point is about 5 times between 1 and 20 min of injection for both cases. As permeability is a quadratic function of the width, the value drops by 3 to 4 orders of magnitude at the fracture tip where the width tends to 0 following a power law asymptote of exponent 1/2 for the toughness regime and 2/3 for the viscous regime. At the end of the stimulation, the maximum fracture radius achieved is in the order of ∼180 m.

Initial Stimulation
Following the base case scenario, we investigated the evolution of the fracture radius and the maximum fracture width at the injection point (Figures 2a and 2c) in time and the effects of different flow rate (Figures 2b  and 2d) at a fixed time (20 min). The fracture radius is insensitive to the temperature of the fluid. It increases rapidly in the early time and tends to stabilize with lower growth rate toward the end of the injection: Propagating a radial fracture becomes more energetically expensive at longer time. Similarly, the fracture width at the injection point increased quickly and stabilized in the later time. On the contrary, the width is sensitive to the fluid's temperature. The final fracture radius is proportional to the flow rate with a higher proportionality up to roughly 100 L s −1 and a lower thereafter. The fracture width at the injection point is sensitive to the combination of flow rate and fluid viscosity: The difference between the two curves in Figure 2d increases at higher injection rate.

Fluid Reinjection Into the Stimulated Fracture
Following the stimulation, the newly created fracture will reduce its opening as the fracture pressure reaches equilibrium with the reservoir pressure. In oil and gas industry, this is prevented by enriching the injection fluid with solid proppants that serve the purpose to keep the fracture open and permeable after the stimulation job. In geothermal operations, because of the very high temperature and the necessity to keep the fracture open for longer periods of time (years to decades), proppants are hardly employed and the fracture will close after stimulation. We have investigated the reopening of a fracture created by a stimulation job during reinjection. The fracture is assumed to be 180 m of radius with an initial opening of 100 μm. Figure 3a shows the case of rapid reinjection with stepwise increase of the injected flow rate in an isolated crack: (i) the pressure rapidly increases during reinjection and as a consequence the crack reopened; (ii) the crack reopening increases the permeability and reduces the overpressure by increasing storage. During continuous and constant flow rate reinjection (dotted line) the pressure monotonically decreases and reaches a plateau relatively quickly, which is temporary. If the injection continues for a prolonged time (Figure 3b), the opening reaches an equilibrium with the deformation of the surrounding rock contrasting the pressure inside the crack. The pressure response shows a minimum at which the crack is completely filled with fluid.
Once the crack is filled, continuous injection pressurizes the crack further, with a linear slope that is a function of the whole system compressibility and a stepwise increase in the flow rate from 25 to 40 L s −1 has a minimum effect in terms of pressure response. If the crack is connected with a permeable fracture network that offers no resistance to flow, no pressurization occurs as the injected fluid can escape with no resistance other than frictional viscous forces (dotted line in Figure 3b). Figure 3c (light blue line) shows that after ∼4 day, the pressure within the fracture reached the total far field stress (σ 0 ¼ 50 MPa), a tensile state developed at the fracture tip, which would eventually overcome the tensile strength, at which point the fracture would repropagate. Comparing the result against the case of a 930 m long crack (obtained after approximately 10 days of initial stimulation) shows a short delay of the pressurization response for the longer initial crack (dotted black line in Figure 3c). Pressurization increases linearly in the log-log plot with a slope almost independent on the fracture length.
The crack width profile during continuous reinjection shows no difference between the nonconnected case ( Figure 3d, solid line) and the connected case (Figure 3d, dots) in the early times, that is, when the crack was not filled and pressurized yet. Once the crack approaches the beginning of repressurization, the responses from the two cases start to diverge, with the width of the connected crack being smaller than the nonconnected one. The connected crack is more difficult to pressurize once the fluid reopens the crack entirely, as evidenced by the width (Figure 3d) and the pressure (Figure 3e) profiles along the crack.

Geophysical Research Letters
During short-time reinjection, stepwise increase in the flow rate generates additional overpressure and crack reopening, which then reduces accordingly the pressure with time. Figures 4a-4h show on the left the evolution of the deformed FEM model (10 4 times amplified) with crack width and pore pressure at different instants of time; the different time instants of the contour map are shown on the right as a point in the time versus pressure curve from Figure 3a (stepwise injection increase). The crack initially opens as a consequence of the overpressure, which is also transmitted to the rock matrix. The pressure decreases to values lower than the initial pressure at the crack tip because of equilibrium conditions stemming from LEFM. The crack tip advances and the width increases all along its profile during the injection (Figures 4a-4d); after shut-in at 60 L s −1 (Figures 4e-4h), the crack tip continues to advance as the pressure quickly decays with time and the fluid is redistributed toward the crack tip. As a consequence of the fluid redistribution and pressure signal traveling along the crack, the width at the injection point is lower than the tip (crack reopening).

Discussion
Fluid rheology plays a minor role in fracture propagation by decreasing viscosity and increasing leak-off with the temperature increase. Flow rate and total time of the stimulation procedure (i.e., injected volume) instead play a dominant role in determining the final fracture length and its opening profile.
After a stimulation performed at high injection rate, the fracture closes and the permeability decreases by several orders of magnitude. Finite element analyses have shown that fluid reinjection reopens the Figure 3. Fracture reopening during water injection in the stimulated area. A rapid injection case (a) shows the initial pressure decay as a consequence of reopening the existing crack, while continuous long-term injection will later pressurize the crack up to failure (c), unless the crack is connected to a permeable fracture network (b). The width profile (d) along the crack at different time shows the different pressurization efficiency for the nonconnected (solid lines) and connected (dots) case. Similar argument can be drawn for the pore pressure profile (e).
existing fracture and pressure decays with time as the fracture reopens. When the fracture is filled with fluid, continuous injection further pressurizes the fracture at a rate that is a function of the system compressibility (fracture + water + rock matrix). If the fracture is nonconnected, the pressurization continues and the injectivity is a decreasing function of time. When the pressure reaches the fracture propagation pressure again, the fracture will further grow. Nonetheless, the longer fracture (∼5 times longer) with a higher storage capacity reaches the propagation pressure delayed by a few days only when compared to a shorter one. If the fracture is connected to a highly permeable network, the repressurization does not occur and the injectivity reaches a constant and time-independent value.
Our analyses have shown that isolated fractures are unlikely to create the suitable conditions for geothermal development: during injection, because of limited storage, the fracture quickly fills in the order of hours and is continuously pressurized. For a successful exploitation of the geothermal reservoir, it is imperative that the newly created fracture intersects a permeable network such that fluid can circulate between wells. One way to achieve a maximum connectivity is zonal isolation (Yoshioka et al., 2015). Injection into a long section of open hole or slotted liner will likely initiate a fracture at the least resistance point (either the least stress or preexisting weakness) and every subsequent fluid reinjection will redirect the fluid to the already permeable fracture. The fracture is then quickly filled and the injectivity will decrease with time.
The Acoculco rock mass is poorly fractured and heat transfer is mainly conductive with an estimated small penetration depth (<1 m) of the thermal front in the rock mass within the considered time frame (∼10 days). Additionally, since the fluid does not escape the fracture because of the impermeable rock mass, the fluid velocity in the fracture is almost null and thermally induced opening would translate into an increase in permeability that will have no effect on the rate of pressurization. For these reasons, we have neglected the effects of thermal strains. Nonetheless, thermally induced fractures can provide additional surface area and connectivity (Joulin et al., 2020), which may impact the economical viability of EGS. Future studies should investigate the impact of thermally induced deformation in low-permeable formations.
Prolonged stimulation is not an effective strategy unless the fracture finds a permeable fracture network, because the gain in additional storage space (longer fracture and wider opening) marginally shifts the pressurization time and tensile fracturing is delayed most likely by days, according to our analyses. Furthermore, prolonged stimulation in a closed fracture drastically increases the overpressure, posing a great risk of triggering unabated dynamic ruptures in critically stressed faults (Ellsworth et al., 2019;Garagash & Germanovich, 2012). Zonal isolation will instead impose to create multiple fractures within a given section of the borehole, increasing the probability of intersecting and connecting to an existing network of fractures.

Conclusions
We have studied the propagation of a fluid-driven fracture in an ultratight and high-temperature potential geothermal reservoir in Acoculco, Mexico, and the hydromechanical response of the newly created fracture during water reinjection. In this scenario, a large fracture can be achieved as long as injection is provided and the temperature of the injected fluid plays a minor role in the viscous forces that develop during propagation. Fracture connectivity is crucial to avoid excessive pressurization and is necessary to create an adequate flow path in the reservoir by connecting the stimulated fractures with the existing ones. Stimulation protocols of EGS should employ zonal isolation to improve fluid injectivity and obtain economically viable geothermal exploitation. All data are generated through publicly available software, which are referenced accordingly in the text, along with proper links to download the full source codes. The calculations of radial fracture propagation are performed through the Matlab © script provided in Dontsov (2016b). The finite element analyses of water reinjection into a fracture are performed with OpenGeoSys-6, an open source computational platform for the simulation of thermo-hydromechanical-chemical (THMC) processes in porous and fractured media (project's website: www. opengeosys.org; source code: https:// github.com/ufz/ogs). The contribution of F. P. is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) project number PA 3451/1-1. K. Y. acknowledges the GeomInt project funded by the Federal Ministry for Education and Research (BMBF) under grant 03G0866A.