Chapter 28
Dynamics of permeability evolution in stimulated geothermal reservoirs

Joshua Taron, Steven E. Ingebritsen, Stephen Hickman and Colin F. Williams

U.S. Geological Survey, Menlo Park, CA, USA

Abstract

Spatially and temporally evolving permeability fields are fundamentally associated with the operation of enhanced geothermal systems. Indeed, permeability will evolve dynamically in any rock mass subjected to strong alteration of stress and temperature. During nonisothermal fluid injection, thermoelastic stress and fluid pressure changes act upon partially open or hydrothermally altered fracture sets to enhance or degrade formation permeability. The physical processes that drive this behavior are highly nonlinear and interdependent. In this work, we explore the resulting magnitude and patterns of permeability alteration. We utilize a thermal–hydrological–mechanical simulator capable of coupling the dominant physics of shear stimulation. Permeability is allowed to evolve under constitutive models tailored to fractured rock, considering the influence of thermal–hydrological–mechanical stress and elastoplastic shear and dilation in a ubiquitously fractured medium. On this basis, we explore the coupled physical processes that control the evolution of permeability during shear stimulation and long-term evolution of a geothermal reservoir.

Key words: permeability, geothermal, numerical simulation, thermohydromechanical coupling, reservoir stimulations

Introduction

Exchanges between mass, momentum, and energy in hydrothermal systems generate complex physical couplings (Fig. 28.1). In many cases, permeability is the communicative mechanism that determines how these systems evolve and to what extent they are coupled. Thus, while there are a multitude of considerations in the design of enhanced geothermal systems, most relate to the behavior of critically stressed fractures and how their characteristics control the evolution of formation permeability. Minimizing short-circuiting potential (thermal breakthrough), increasing heat-transfer area and fluid residence time, controlling working fluid losses, and mitigating the harmful effects of induced seismicity are primary operational considerations in enhanced geothermal systems, while the mechanisms underlying permeability change are of general scientific interest. Laboratory studies indicate that chemical–mechanical compaction mechanisms (such as pressure solution) may be important at representative enhanced geothermal systems conditions on the order of months (e.g., Polak et al. 2003), while bulk chemical reactions can inhibit or enhance permeability over several years (e.g., Taron & Elsworth 2009). Clearly, shear stimulation and elastic compaction driven by pressure and thermal perturbations are critical concerns at any timescale (e.g., Dempsey et al. 2013).

Fig. 28.1 Physical exchanges between mass, momentum, and energy, and the central communicative importance of permeability in these exchanges.

To address these issues, we perform simulations using a modified version of OpenGeoSys (opengeosys.org). The new simulator is well suited for enhanced geothermal system scenarios. Realistic physics is incorporated with flexible degrees of coupling between the dominant energy conservation, momentum, and mass balance relationships. New fluid equations of state bring the simulation range well above the typical pressure and temperature conditions for enhanced geothermal systems and allow us to consider the influence of the larger scale hydrothermal system, with smooth consideration of liquid/vapor equilibrium conditions at temperatures up to 700 °C and pressures to 100 MPa. Energy, fluid, and solid momentum balance are flexibly combined in any number of “continua” or between materials in local thermal nonequilibrium (i.e., fracture and matrix and/or fluid and solid).

We explore mechanisms involved in the enhancement of fracture permeability, targeting injection pressures that encourage self-propping shear stimulation alongside elastic pressure and thermal dilatation of fracture sets. The elastoplastic response is governed by a Mohr–Coulomb failure criterion and flow rule utilizing a full Newton stress integration scheme. Fracture shear and dilation are linked directly to the failure model with a nonlinear model for dilation angle. Mechanical stiffness of the fractures is modified in response to dilation and relative to a nonlinear fracture compaction model, subsequently impacting bulk stiffness of the system.

Because fractures have a dominant influence on bulk rock properties, several decades of research have focused on their hydrological–mechanical properties. Topics include compressive strength of fractured rock (Goodman 1976; Brown & Scholz 1986; Barton et al. 1985), the effect of normal displacements and fracture roughness on permeability (Goodman 1976; Elsworth & Goodman 1986; Gangi 1978), the effect of shear displacement on permeability (Olsson & Barton 2001; Yeo et al. 1998), and relationships between mechanical and hydraulic aperture (Zimmerman & Bodvarsson 1996; Renshaw 1995; Piggott & Elsworth 1993). Mineral precipitation and dissolution have been studied as well. Two notable works relate to localized dissolution (Detwiler 2008) and heterogeneity of mineral assemblages (Peters 2008a).

Fractures under elevated stress

In enhanced geothermal system and other environments, permeability change is dominated by thermoelastic compression, creep (chemical and/or mechanical), shearing, and mineral infilling or dissolution along fractures. Fracture elastic compaction under stress is generally addressed with (i) a laboratory-derived stress-closure model, or (ii) an elastic contact model based on estimates of real contact area between surfaces of prescribed roughness. Both methods are empirical; the former relies on laboratory information regarding stress and closure and the latter on laboratory information regarding contact area and surface roughness (fracture profiles) and elastic rock properties.

As fractures dilate or compact (by any mechanism), the contact area between opposing fracture surfaces is altered, leading to changes in fracture apparent stiffness. This is accommodated directly in contact models, and indirectly in stress-closure models. Here, a fracture set is considered to respond nonlinearly to an evolving stress field through the often-used relationship for mechanical fracture aperture,

28.1

dating at least to Goodman (1976), in a slightly different form, where is residual aperture, is maximum aperture, is the effective normal stress, and α is the stiffness coefficient. Similar forms have been introduced to include thermal–chemical behaviors (Min et al., 2009). Moreover, Liu et al. (2009, 2013) have derived such an exponential relationship theoretically and shown strong agreement across a wide range of laboratory data.

In addition to Eq. 28.1, another important product of Goodman's (1976) experiments is the observation that “mated” (zero shear offset) and “unmated” (finite shear offset) fractures exhibit the same closure behavior (i.e., can be represented by Eq. 28.1), only differing proportionally with the constants, and (see also Brown & Scholz 1985). This suggests that the contact area relationship for any fracture (plot of contact area versus mechanical closure) after some amount of shear will be an extension of the mated solution. Taron & Elsworth (2010a) performed numerical experiments (utilizing profiles of several laboratory-scale fractures, roughly 20 cm × 15 cm) to quantify this behavior. Fracture surface profiles were separated to zero overlap, numerically sheared at small increments, and closed again to calculate a new relationship for closure at each degree of offset (Fig. 28.2). The closure curve at each offset (as a function of shear distance, x) can be fit with the exponential function, which is identical in form to the relationship presented by Rutqvist et al. (2002) and shown in Eq. 28.1, with effective stress substituted for R_c (ratio of contact area to total square area of the fracture),

28.2

Eq. 28.2 is akin to introducing an angle of dilation for fracture shear, but the correction is nonlinear and bimodal. Laboratory experiments show similar behavior; intuitively, fractures dilate (elastically) under a nonlinear and reproducible rule (e.g., Lee & Cho 2002), while another work (Rutqvist 2015) has successfully matched a similar relationship to in situ data. In this work, we use the direct aperture form of this relationship, which is quite similar to the relationship of Rutqvist (2015),

28.3

Several relationships are available to produce such a response. In this chapter, we utilize a hyperbolic tangent Heaviside type function, with results illustrated in Figure 28.3 for one potential parameter set,

28.4

where x_min is the shear distance to begin dilation (break in), x_max is the shear distance after which dilation no longer occurs, del is a factor controlling the nonlinear angle of dilation during this period, b₀ is the aperture with zero shear, and db_max is the amount of dilation at x_max. An alternative, but similarly behaving, relationship was expressed in Lee & Cho (2002), and, as a simplification, the same dilation is applied to and . Shear slip is computed from inelastic shear strain in the mechanical model. For a discussion of the scale dependence of such behaviors, see Taron & Elsworth (2010a), who utilized fractal upscaling, and Rutqvist (2015) for a more detailed discussion.

Fig. 28.2 Fracture surface profiles separated to zero overlap, numerically sheared, and closed again under stress. The fit lines are produced by Eq. 28.2 (modified from Taron & Elsworth 2010a).

Fig. 28.3 General relationship between fracture dilation and shear displacement used in our model. The dilatation parameters illustrated here are consistent with (but slightly lower than) those obtained in the laboratory experiments of Lee and Cho (2002), and were allowed to vary in our model.

Consider one further numerical experiment. As a fracture is offset, its contact area is decreased, leading to a decrease in apparent stiffness and, in the presence of chemical and/or mechanical creep (i.e., Polak et al. 2003; Tada & Siever 1989; Beeler & Hickman 2004; Taron & Elsworth 2010a,b), an increased tendency for irreversible closure due to increased contact stress. In Figure 28.4, a single representative fracture is offset discretely and allowed to creep under its new conditions of contact area (using the model from Taron & Elsworth (2010b). Because the contact area decreases with each offset, it becomes more susceptible to compaction and ultimately achieves a new equilibrium closure. While creep is not considered here, the results are directly analogous to elastic closure, where each exponential closure segment represents a range of stress, rather than equilibration time. In other words, as the reservoir undergoes shear stimulation, compliance and general characteristics of the responsible fracture sets are both reversibly and irreversibly altered.

Fig. 28.4 Representative fracture profile offset in increments and allowed to creep under new conditions of contact area. Results are directly analogous to elastic closure for each exponential closure segment (modified from Taron & Elsworth 2010a).

Under an anisotropic stress state, the resulting aperture may have directional components of permeability. Under such anisotropy, the cubic law may be expressed as (Bear 1972),

28.5

for fracture spacing s and the diagonal components of the permeability tensor k. The derivative of Eq. 28.1 with respect to stress gives fracture stiffness. Here, b should be the hydraulic, rather than mechanical, aperture (e.g., Zimmerman & Bodvarsson 1996). The subscript of b and s refers to the direction orthogonal to the respective quantity, and the subscript on k refers to the direction of flux. Thus, permeability in x depends on the apertures on the xy () and xz () planes. Note that for the isotropic case this reduces to the cubic law, . Later, we assume isotropic spacing in an orthogonal network, with apertures derived from the stress acting orthogonal to their respective plane.

Methodology

Conservation of mass

Here, we present, for brevity, only the basic forms for the case of single porosity. The general mass balance for a fluid of multiple phases is, for phase ,

28.6

for bulk fluid density , and density, , saturation, S, porosity, , and fluid velocity , and where superscripts l and v refer to the liquid and vapor phases, respectively, and q_p is a pressure source term. Considering Darcy's law for flow and substituting the solid mass balance into the aforementioned relationship, we can write

28.7

28.8

where and , for the Biot–Willis coefficient , grain bulk modulus , bulk linear coefficient of thermal expansion , and granular thermal expansion , where is relative fluid velocity with respect to the solid for solid velocity . The form of the density relationship is important, especially when working with enthalpy, h, whose derivative is not merely temperature dependent,

28.9

Note that in multiphase conditions (i.e., if the current state of the system is that of liquid–vapor equilibrium) the density derivative must include the derivative of vapor saturation with respect to pressure and enthalpy. All fluid properties and their derivatives in water systems are calculated directly from the IAPWS industrial formulation (IAPWS 1997).

Conservation of energy

Both enthalpy (h)- and temperature (T)-based derivations begin with the general balance of specific internal energy (e), where superscript s refers to the solid, considering both viscous dissipation and pressure–volume work terms (e.g., Garg & Pritchett 1977; Faust & Mercer 1979),

28.10

Now taking and canceling terms,

28.11

Lastly, by writing , and substituting the fluid mass balance Eq. 28.8 into Eq. 28.8 to remove many expansions,

28.12

which maintains thermal equilibrium between fluid and solid. Alternatively, the solid temperature balance may be extracted and treated separately. Here, is solid specific heat capacity at constant pressure, and similarly for the fluid.

The temperature formulation is obtained from Eq. 28.11 by recognizing . This leads also to the fact that, in enthalpy form, must be numerically discretized with respect to both enthalpy and pressure. Note the left-hand side term that arises in Eq. 28.12 accounting for any fluid mass sources, . This term is highly nonlinear and must be linearized when pressure and energy are solved in the same global equation. To consider multiple continua or local nonequilibrium, transfer functions must be included in the fluid and energy balance relationships. For energy, the exchange is governed by the interphase or intercontinuum heat-transfer coefficient, , where in Eq. 28.12 is the result of an additional energy balance partial differential equation for that continuum or phase. With enthalpy as a primary variable, since the interphase transfer is temperature based, a linearization is required for this term, unless it is to be placed on the right-hand side. Both of these terms are strongly coupled between respective balance equations, and so left-hand side coupling (monolithic) is important in most hydrothermal cases. See, for example, Galet (2011), for a thorough discussion of interphase transfer.

Conservation of momentum

In the usual manner (i.e., Zienkiewicz & Taylor 2005; Jaeger et al. 2007), we can write solid momentum equilibrium as,

28.13

where and determines the effective stress relationship, where for isotropic elasticity using the Biot–Willis coefficient along with moduli for the bulk solid, , and grains, . Also introduced is the stress–strain constitutive relationship and the strain–displacement matrix B, for the total, , and effective, , stress, and the total, , thermal, , and plastic, , strain. Note that we solve this system as incremental in strain and displacement (in keeping with small strain approximations), but that thermal strains are relative to a fixed thermodynamic initial state, plastic strain to an undamaged state, and elastic strain relative to a state of zero stress. In other words, the stress–strain relationship is not incremental, which is important when moduli evolve with strain and which we find introduces less numerical error. For a Mohr–Coulomb plasticity model, the basic relationships are nonassociated in shear and given by

28.14

Tensile flow follows the usual associated form (i.e., Jaeger et al., 2007), while for the case when two flow functions are simultaneously active (i.e., two shear yield criteria are met, which is common in two dimensions), we adopt the Koiter (Koiter 1953; De Borst 1987) composite flow rule:

28.15

In these relationships, is the failure criterion (for a friction angle ) and the flow rule (nonassociated reproduction of , with N_d in place of N_f, and with dilation angle ), while is the “plastic” stress change required for return to the yield surface, is plastic strain, are cross-coupling coefficients for return to the α and β surfaces, is the matrix of elastic coefficients, and is the flow multiplier resulting from the plasticity integration. The Drucker–Prager form exhibits only a single smooth yield surface, amenable to an elegant solution in deviatoric stress invariant space (Westergaard (1952) coordinates), while we thus far find the multiple surface Mohr–Coulomb (MC) model to be preferable to a smooth surface or corner smoothing MC model (see Zienkiewicz & Taylor 2005; Sloan & Booker 1986) for smoothing model suggestions). All results here use the MC form shown earlier.

Coupling strategy

A fully coupled thermal–hydrological–mechanical–chemical simulator is often discussed as a critical step toward reliable prediction both of the results of enhanced geothermal system stimulation and of long-term enhanced geothermal system reservoir performance. However, “fully coupled” should refer to the physics of the problem (i.e., strong adherence to coupled terms between partial differential equations, as presented earlier) rather than to the numerical scheme used to combine necessary equations or the number of simulators that must be combined to achieve this. Although there are certain advantages to a single, self-contained simulator, in most scenarios we have explored, a full monolithic solution (all equations solved by a single call to a linear solver) is required only between those equations that are most strongly coupled in a given scenario (typically fluid mass and energy balance). Often a sequential approach is valid or even preferable for some of the equations (i.e., mixed sequential approach). To allow for these various approaches, we allow flexibility in how and where coupling occurs. In a general residual form of finite element method temporal discretization,

28.16

or

28.17

where M is the mass accumulation matrix, J the Jacobian (not present if using a Picard iterative solution), K the stiffness matrix, and R any right-hand-side terms (not residual) of all active equations in a local element. Note that M and K must be evaluated at time . In this formulation, a particular partial differential equation coefficient can be applied either on the left- or right-hand side of the balance equations presented earlier. For instance, if fluid and heat flow are to be solved in the same global equation, the term from Eq. 28.8 can either be placed in M (left-hand side) or in R (right-hand side). There is no additional computational overhead in choosing one or the other, other than the potential need for a greater number of iterations due to the nonlinearity that is now weakly coupled if placed in R. Furthermore, any number of equations can be solved in an iterative loop. Thus, any coupling strength can be dealt with in the appropriate manner. A full Newton–Raphson finite element method solution is rarely necessary in our experience, provided the left-hand side expansion includes at least all first derivatives. Instead, a Picard iteration works quite well provided that the time step control considers the error in an appropriate manner and the full forms of the equations in “Methodology” Section are utilized. Equation 28.16 can obviously be found in many sources, but Bergamaschi and Putti (1999) provide a good discussion.

Results

Utilizing this system of finite element method discretized equations and the permeability constitutive models discussed previously, we now explore the permeability response of a critically stressed, fractured geothermal reservoir. We take a single geometric conceptualization for all scenarios (Fig. 28.5). A radial symmetry model is adopted, with a 300-m-tall open-hole injection interval represented as a line source. Owing to vertical variations in porosity, with constant solid density and state-dependent fluid density, the initial stress magnitudes and fluid pressures are nonlinear functions of depth and must be equilibrated numerically.

Fig. 28.5 Simulation domain: Reservoir is more highly fractured than caprock or basement. Parameters are shown for the previously stimulated reservoir (“Long-term observations” Section).

A lateral constraint boundary is used to first determine a vertical stress profile and, simultaneously, a pore-pressure distribution with depth assuming rock matrix density of 2600 kg m⁻³, depth-varying porosity, an initial temperature gradient of 100 °C km⁻¹, and properties of pure water calculated at Gaussian quadrature points from the IAPWS 1997 Industrial Formulation equation of state (IAPWS 1997). Then, a friction coefficient of 0.6 is used to apply horizontal stress to the outer boundary consistent with zero cohesion and stress limited by frictional failure on optimally oriented normal faults (e.g., Hickman & Davatzes 2010), wherein both horizontal principle stresses are assumed equal in magnitude. The lateral constraint is then removed, while the equilibrated fluid pressure is also applied to the outer boundary and the surface pressure removed to allow open flow at the surface. Because we are using an enthalpy formulation for the energy balance, after setting up an initial temperature profile the geothermal gradient is recalculated at each equilibration step to ensure an accurate inversion from temperature to enthalpy. This geothermal gradient is maintained at the outer boundary as well. The rather high initial geotherm of 100 °C km⁻¹ is used to generate temperatures similar to enhanced geothermal systems in the reservoir zone (250 °C at the center of injection), resulting in temperatures and pressures in the lower basement that are supercritical (Fig. 28.5).

Fracture parameters are highly variable and should be chosen with reference to field data if a particular site is to be examined. Our baseline fracture parameters were selected to be representative of very weathered surfaces – not highly compliant initially, but susceptible to large (relative) aperture increases when subjected to shear dilation. Unstimulated reservoir characteristics include low compliance fractures (∼15–20 µm max aperture and potential closure of ∼10 µm) on a 1 m orthogonal spacing to produce a desired in situ permeability of ∼ at the center of the injection interval. Laboratory evidence suggests that natural features of this scale remain permeable and compliant to stresses well above those examined here (Pyrak-Nolte et al., 1987). Alternative combinations are possible, such as larger fractures occluded by mineral deposition or simply with larger spacing, but these alternate parameterizations are not explored here.

Short-term stimulation behavior

Here, we explore the near-well region during a high-injection-rate shear stimulation exercise. The fluid injection temperature for these tests is 70 °C and fractures are cohesionless with a friction angle of 32° (coefficient of friction of 0.6). A constant pressure/temperature boundary is used at injection. This allows tight control on the fluid pressure, while staying below the least principal stress, and mimics what is sometimes attempted during a controlled shear stimulation enhanced geothermal system exercise (Chabora et al. 2012). However, this approach introduces a discontinuity in the pressure field, which can cause stability problems with plasticity. The injection pressure is applied on a mild Gaussian taper over the injection zone to reduce the discontinuity at each end of the open-hole interval. Injection pressure is set at 6 MPa above the in situ pore-fluid pressure at the center of the injection zone, or roughly 2 MPa below the least principal stress. The initial, unstimulated target formation permeability is roughly , about 5 times larger than the prestimulation values at Soultz and Basel (Evans et al. 2005; Haring et al. 2008).

Dual temperature control is used, consisting of a fracture and intact rock domain. Because the timescale is short, heat transfer between fracture and rock becomes important. For pressure, we assume zero storage in the matrix: either because the pressure diffusion into the matrix is rapid relative to other processes, and thus equilibrates quickly with fracture pressure, or because the matrix is of such low porosity and/or permeability that it does not contribute strongly as a storage mechanism on short timescales. To do this, we introduce an additional term in the pressure equation that considers the different temperatures, and thus different fluid properties, in the fractures relative to the matrix; this is a mass accumulation term. This is a hybrid approach falling between a local thermal nonequilibrium simulation and a full dual-porosity simulation. The thermal–mechanical response to cooling in such a simulation is controlled by the temperature of the matrix, not the temperature of fluid in the fractures.

A first analysis shows the relative changes of parameters at two locations on a centerline outward from the injection well (Fig. 28.6). As expected, the fluid pressure increase and resulting dilation is rapid and significant while the thermal response is delayed and clearly distinguishable following passage of the pressure front. Importantly, a second, smaller pressure increase slightly precedes the second dilation phase (at roughly 1.2 days 2 m from injection and 1.5–2 days at 4 m). This indicates that permeability is increasing nearer the well, allowing the advance of a higher-pressure front, and corresponds closely with the arrival of colder injection fluids. Permeability increase is allowing greater rates of heat and pressure diffusion, which feeds back to higher rates of permeability increase due to both thermal and hydrological–mechanical effects.

Fig. 28.6 Parameter evolution over time for the shear stimulation (short-term) case. At 2 m: temperature (250–70 °C), pressure (21.4–27.3 MPa), slip (0–3.0 mm), and permeability (). At 4 m: temperature (250–98 °C), pressure (21.4–27.2 MPa), slip (0–2.6 mm), and permeability ().

Next, in Figures 28.7 and 28.8, we compare two cases with differing thermal transfer efficiencies between fracture and rock. Both have identical bulk permeability, but in one case heat transfer from fracture to matrix is more efficient and in the other the heat-transfer coefficient is reduced by a factor of 10. This latter case is similar to what would occur with larger fracture spacing, creating more inefficient or slower heat transfer.

Fig. 28.7 Efficiency of cooling, expressed in terms of fluid and rock temperature very near the injection point.

Fig. 28.8 Shear stimulation at two heat-transfer efficiencies, 2 m from injection. Rapid transfer (frequent spacing or higher efficiency) is capable of creating secondary slip through rapid cooling. Slow transfer (wide spacing or lower efficiency) cannot create the necessary thermal gradient for secondary slip.

Figure 28.7 shows the difference in temperature between fracture and matrix at the injection point. Note that the fluid temperature in the fractures is immediately constant at the injection value, while slow diffusion into the matrix introduces a delayed thermal response. Predictably, the more efficient simulation allows for much more rapid cooling of the matrix than the inefficient simulation.

Figure 28.8 explores the inelastic slip and permeability response in these two cases. Both cases display immediate shear failure and resulting permeability gain from the hydraulic reduction in effective normal stress, amounting to about 60% of the total permeability change realized during the simulation with better heat-transfer efficiency. The thermally inefficient case, however, experiences no additional dilation after this initial shear failure. The rate of heat transfer to the matrix is not sufficiently rapid to induce an inelastic response. The efficient heat-transfer case (corresponding to tighter fracture spacing) experiences an additional 40% gain in permeability due to the combined effects of shear failure and thermal contraction. It has a secondary slip mechanism as a result of the sharp thermal gradient in the rock, despite the identical initial permeability field.

Long-term observations

Now, the baseline values are made more compliant by introducing a small amount of initial shear, effectively assuming that the entire reservoir has been previously stimulated. Initial bulk permeability at the center of the injection interval is roughly . This value is near the lower end of commercial geothermal reservoirs, which range from roughly to , with permeability generally decreasing with temperature over an observed temperature range of 200–350 °C (Bjornsson & Bodvarsson 1990). This value is also similar to the poststimulation permeability at 2.85–3.4 km depth at Soultz (: Evans et al. 2005) and at 4.6–5.0 km depth at Basel (: Haring et al. 2008). Fractures are considered cohesionless with a friction angle of 34°, corresponding to a coefficient of friction of 0.7.

Fluid pressures are kept below the minimum principle stress and not greatly in excess of values required for Coulomb (frictional) failure, at an injection rate of 10 kg s⁻¹. This injection rate is maintained at a temperature of 80 °C (not as a thermal Dirichlet boundary, but as a corresponding energy flux). Local thermal equilibrium is assumed for the simulations in this section, in that the same temperature is applied to fracture and adjacent matrix and fluid and rock. The injection is a line source that is slightly smoothed to zero flux at the ends (on a mild Gaussian taper), which removes some of the tendency for shear slip to concentrate at the lower and upper ends of the well; some failure does initiate here nonetheless. The overall behavior is visible in Figure 28.9, where arrows show the behavior of solid displacement throughout the domain as the injection region cools and contracts.

Fig. 28.9 Simulation domain: Injection process showing fluid density within the domain of Figure 28.5 after 20 years of injection; note the transition to lower supercritical densities beneath the reservoir. Arrows represent displacement of the solid (elastic and inelastic); rock is contracting toward the cool injection area creating areas of shear.

The first areas to fail in shear are the ends of the injection interval, albeit over only a small area, followed by failure along the full length of the injection interval, and then mild shear advancing radially away from the well. Over time, the cold zone becomes large enough to disturb the stress field in a more significant way and plasticity begins to concentrate at the outer edge of the thermally perturbed zone (Fig. 28.10B). There is no gap in plasticity from the well to the outer edge of the advance (Fig. 28.10), only a larger accumulation at greater distance. Slip occurs in areas exhibiting the highest thermal gradient: the advancing cold boundary. The shape of the thermal plume influences the spatial concentration of frictional failure; a point injection with radial advance will look different from the line injection with linear advance shown here.

Fig. 28.10 Total inelastic slip distance (m) on the fractures within the domain of Figure 28.5 after 6 months (A) and 20 years (B). Area shown is roughly 500 m in height and 650 m in width.

Figure 28.11 shows reservoir response in the center of the reservoir at two different distances from the injection well. The pressure pulse has a minimal impact on total permeability change at any considerable distance from injection. The thermal/plastic response is delayed, but leads to increases in permeability starting at about 1–2 years. In the most active areas, permeability increases by a factor of ∼3. Adjustment of fracture stiffness and the fracture dilation curve would, of course, modify this value. In this stimulation, the fractures are initially sheared to represent an already stimulated system; the plastic dilation response would be different if this were not the case.

Fig. 28.11 Evolution of parameters with time at two radial distances from the center of the injection interval. At 25 m: temperature (250–85 °C), pressure (21.4–21.7 MPa), slip (0–0.6 mm), and permeability (). At 60 m: temperature (250–136 °C), pressure (21.4–21.6 MPa), slip (0–0.8 mm), and permeability ().

We highlight a few primary observations regarding the permeability response. As seen in Figures 28.10 and 28.11, shear failure concentrates at the boundary of the advancing cool region. Thus, the onset of shear failure and dilation occurs before a particular point cools significantly. The ensuing permeability change is associated with the magnitude of temperature change, not with the onset of fracture dilation. While dilation opens the fracture, it also increases its compliance, resulting in minimal net observable increase in permeability. However, once the cold front arrives fully, the fracture is more complaint due to its previous slip, and thus dilates more easily as the thermal stress declines. It is safe to assume that larger scale fractures will behave differently due to their tendency to dilate at a greater rate.

Conclusion

Simulation of a short-term, high-injection-rate enhanced geothermal system stimulation explored the importance of hydraulic versus thermal effects and the impact of a temperature differential between fluid in the fractures and the adjoining rock matrix. Under less efficient thermal transfer (corresponding to larger fracture spacing), thermal effects are relatively unimportant during short stimulation operations, inducing no significant additional shear-enhanced permeability gain. Alternatively, if the matrix cools more rapidly (due to tighter fracture spacing), additional shear stimulation and stress relaxation from thermal effects add another order of magnitude to the total permeability change. This is a two-stage slip mechanism. Importantly, these two very different results occurred within a reservoir represented by an identical initial permeability field. Thus, the effects of variations in fracture spacing on the evolution of permeability is at least as important as the magnitude of initial permeability, and quantifying the driving mechanisms is critical to understanding reservoir response. Thermal effects are probably always important. Even without significant mechanical influence, they drive large changes in fluid properties; their impact on fluid density can introduce buoyancy effects, promoting downward migration of colder injection fluids, while viscosity differences influence the inelastic response via pressure diffusion.

Simulation of long-term enhanced geothermal system operation demonstrates the ability of sustained pressure/thermal perturbation to drive fracture slip and corresponding permeability gain. The most active inelastic areas exist on the boundary of the advancing cold injection fluid, at the location of the highest thermal gradient. In the case of small-scale fractures, inelastic slip does not significantly impact permeability, at least directly, at any significant distance from the injection point. However, such slip may weaken fractures in shear by altering their contact area, allowing them to respond significantly to sustained cooling by relaxing the thermal stress field.