Abstract
Massive quantities of fluid are injected into the subsurface during the creation of an engineered geothermal system (EGS) to induce shear fracture for enhanced reservoir permeability. In this numerical thermoelasticity study, we analyze the effect of cold fluid injection on the reservoir and the resulting thermal stress change on potential shear failure in the reservoir. We developed an efficient methodology for the coupled simulation of fluid flow, heat transport, and thermoelastic stress changes in a fractured reservoir. We performed a series of numerical experiments to investigate the effects of fracture and matrix permeability and fracture orientation on thermal stress changes and failure potential. Finally, we analyzed thermal stress propagation in a hypothetical reservoir for the spatial and temporal evolution of possible thermohydraulic induced shear failure. We observe a strong influence of the hydraulic reservoir properties on thermal stress propagation. Further, we find that thermal stress change can lead to induced shear failure on nonoptimally oriented fractures. Our results suggest that thermal stress changes should be taken into account in all models for long-term fluid injections in fractured reservoirs.
1. Introduction
One of the primary driving mechanisms for permeability creation in engineered geothermal systems (also known as enhanced geothermal systems), or EGS, involves shear failure induced by fluid injection at high pressures. In environments with low differential stress, tensile fractures may develop if the injection pressure exceeds the minimal principal stress (e.g., fracking). The injection of cold water into a reservoir at substantially higher temperature also induces thermal stress changes that contribute to the overall evolution of the local stress and failure potential. This rapid cooling of the reservoir can lead to thermal cracking and thus further enhance permeability [1–4], but thermal fracturing, fracture propagation, or fracture reactivation may also contribute to premature cold water breakthrough into producing wells.
The basis for EGS is usually geothermal plays of the “hot dry rock” type where the available water in the porous medium is considered negligible. These conditions are found primarily in metamorphic or igneous terrains with low permeability and porosity containing fractures and faults that provide the major pathways for fluid flow. In geothermal energy systems, the fracture’s surfaces serve as the main heat exchanger. Clearly, preexisting, critically stressed, and optimally oriented fractures provide the most favorable conditions for enhancing permeability of EGS [5, 6].
In this paper, we focus on the role of thermal stresses during cold fluid injection and stimulation of an EGS site. Of special interest, here is the interplay of hydraulic and thermally induced stresses. The processes involved in permeability creation during hydraulic stimulation act on different timescales. While the poromechanical coupling is active throughout the injection, its dominance over thermomechanical effects depends on the state of the injection. Thermomechanical coupling plays a particularly important role during prolonged periods of injection (weeks to years) because the variation of injectivity with injection water temperatures can be attributed to thermal stress [7].
A major concern in EGS is induced seismicity at levels above that tolerated by the local population, in either frequency or magnitude. Usually induced seismic events are attributed to the change in effective stress due to the change in fluid pressure [8–10]. However, thermal stress also significantly contributes to induced seismicity in petroleum and geothermal fields [7, 11, 12]. Stark [12] found that, in the Geysers geothermal field in northern California, USA, half of the measured earthquakes appear associated with cold water injection.
In this paper, we present the theoretical basis for thermal stresses, evaluate the temperature distribution during injection in a borehole, and determine a modeling framework for evaluating the influences of thermal stress generation and propagation in a hypothetical reservoir. We describe the numerical method and present results of numerical experiments focusing on the influence of thermal stress on permeability, fracture orientation, and failure potential. We discuss the results in terms of thermal influence on induced seismicity and reservoir characteristics.
2. Theory
Here we present the mathematical basis for the formulation of thermal stresses. In addition, the temperature profile in an injection well is considered because this has a significant impact on the initial conditions of the numerical simulations.
2.1. Mathematical Description of Thermal Stress
A body will change its shape and/or volume when exposed to a temperature change . This change is called thermal strain and can be expressed aswhere is the coefficient of linear thermal expansion in 1/K. In most materials, is positive and on the order of . In isotropic materials, the thermal strain acts only on the normal strains with the same magnitude. If the body’s deformation is restricted, as it would be the case for a small volume inside a rock mass, the strain results in thermal stress.where is the elastic stiffness tensor of the material. If we restrict ourselves to isotropic conditions, the thermal strain has only normal components with equal magnitude in which case we can simplify the previous expression towhere is Young’s modulus (Pa) and Poisson’s ratio (—). It is important to note that (3) is only nonzero for the three normal stresses . It is immediately obvious that larger temperature differences will result in higher thermal stress changes. Additionally, the thermal stress is positive (relative compression) if the temperature difference is positive (), and if the temperature difference is negative, the thermal stress is negative (relative tension). The magnitude of thermal stress can change widely depending on the material. Assuming a constant thermal expansion coefficient and fixed temperature difference of , the resulting thermal stress for an elastic sandstone ( = ~20 GPa) is smaller than the resulting thermal stress in a typical granite ( = ~50 GPa). The granite would undergo a stress change of 1 MPa in this case compared to a 0.4 MPa stress change in the sandstone.
Thermal expansion coefficients are well-constrained by experiments and show only minor influences of temperature and pressure on the thermal expansion coefficients [13, 14]. Cooper and Simmons [15] attributed some of the change in the thermal expansion coefficient to the formation of microcracks by differential expansion of mineral grains. Considering the small magnitude in the change of the thermal expansion coefficient compared with the order of magnitude expected in temperature and pressure change, it is a valid assumption that the thermal expansion coefficient is constant.
In the following, we assume that the thermal stress is independent of the fluid pressure and the in situ stress state of the rock. Thus, the resulting stress can be obtained by superposition of the effective stress () and the thermal stress. Changes in the in situ stress of the rock are negligible on the timescales of interest for hydraulic stimulation. Considering only stress changes resulting from pore pressure and thermal expansion, we can formulate the total stress change asClearly, other stress contributions as slip-induced stresses and stresses induced by chemical reaction have to be considered in a general case. However, for reasons of simplicity, we restrict ourselves to only pore pressure and thermally induced stress changes.
2.2. Induced Shear Failure Potential
Induced shear failure potential is estimated by a Mohr-Coulomb failure condition. We restrict ourselves to a cohesion-less material with a friction coefficient of . A fracture segment is able to slip and is thus categorized as “potential slip” if the following condition is met:The effective normal stress is defined aswhere is the normal stress, is the fluid pressure, and the thermal stress as introduced in the previous section. The effective normal stress and shear stress acting on a fracture segment arewhere and are the maximum and minimum principal stresses acting in the far field and is the angle between and the fracture segment measured from the normal to the plane.
The poroelastic deformation of the fractured reservoir during the injection as well as the deformation due to fracture slip is not included in the present model. Consequently, this simple model does not predict magnitudes or estimation of the amount of slip but rather identifies when a frictional failure condition is met, similar to other approaches in modeling EGS [9, 10].
2.3. Heat Distribution in a Geothermal Well
The temperature inside an injection well is not constant with depth. The injected water is heated by the rock mass surrounding the borehole while moving downwards through the borehole. Although the heat distribution in geothermal boreholes can be measured, in many cases it is still useful to describe it mathematically. Such calculations can be used in numerical simulations and aid the drilling crews during their operation. In most cases, description of well bore heat transmission is based on a well bore heat balance equation. Most of the literature is based on the initial work of Ramey [16] in which they derived the temperature distribution in a well used for hot fluid injection. The work was later enhanced by the rate of heat loss from the well to the formation by Ramey [17]. Recent work from Hagoort [18] reevaluated Ramey’s classical work and found that it is an excellent approximation.
Following Satman and Tureyen [19], the temperature in an injection well into which a single-phase fluid is injected is given byHere, is the distance downwards from surface in meters, is the geothermal gradient in , is the surface temperature in , and is the temperature in of the injected fluid. is a variable defined aswhere is the mass flow rate in ks/s, the heat capacity of the fluid (assumed constant), and the thermal conductivity of the formation (also assumed constant with depth). The dimensionless function describes the transient heat transfer to the formation. There are a number of different formulations for function available. Kutun et al. [20] provide a simple formulation:which is based on a best curve fit of the data provided by Ramey [16] and Ramey [17, 21]. According to Satman and Tureyen [19], it is accurate within 1% for the relevant timescales. In (10), is a dimensionless time defined bywhere is the mean thermal diffusivity of the formation in and the well radius in meters.
A detailed review on the methods to describe well bore heat distributions, prevalent assumptions, and different formulations for the transient heat transfer function can be found in Satman and Tureyen [19].
2.4. Heat Distribution in a Geothermal Well
Based on (8) to (11), the general characteristics of fluid injection in the well are presented and examined in detail. The model parameters are given in Table 1. Figure 1 shows the temperature distribution with depth of well bore temperature during prolonged fluid injection. The profiles are plotted for different injection times: 1 d, 10 d, 30 d, and 365 d. The geothermal gradient shown as a dashed line represents also the static (long term) temperature distribution in the well bore. It is clearly visible that it is not a static profile but contains a dynamic evolution with time. During the first days of injection, the temperature profile changes most rapidly. After about 10 days of continuous injection, the rates of the temperature profile change decrease. After 365 days of continuous injection, the bottom hole temperature has decreased from after 1 day of injection to in the case of an injection rate of . When the injection rate is lower at , the temperature decreases from initially after 1 day of injection to after one year. This demonstrates the highly dynamic temperature distribution during injection with respect to time and injection rates.
(a)
(b)
The dynamic behavior is driven by advective heat transfer during injection of the fluid at the well head and conductive heat exchange with the formation. Initially heat transfer from the formation to the fluid is high, leading to a rapid increase in the fluid temperature. During the course of injection, the formation is cooled by the injected fluid leading to reduced conductive heat transfer from the formation to the fluid. This subsequently lowers the overall temperature distribution in the well bore.
Figure 1 shows that the assumption of a constant temperature boundary in numerical modeling of fluid injection is not valid either in space or in time. Especially for very low injection rates (not shown here), the thermal gradient in the borehole is comparable to the static geothermal gradient. In long duration high injection rate scenarios, a constant temperature at the borehole becomes more acceptable as the formation is cooled quickly and the well bore temperature approaches the injection temperature. From the findings of this section we can conclude that the temperature change in the well bore is significant and must be taken into account when thermal stresses are to be evaluated accurately.
3. Methods
We developed an embedded discrete fracture model (EDFM) using the phenomenology described in other studies [22–24]. The conceptual idea of the EDFM is the distinct separation of a fractured reservoir into a fracture and a damaged matrix domain. We introduce a transfer function to account for coupling effects between the two domains (Figure 2), so the fracture and matrix domains are computationally independent except for the transfer function. As the fractures are generally very thin and highly permeable compared to the surrounding matrix rock, the gradient of fracture pressure normal to the fracture is negligible. This allows for a lower dimensional representation of fractures (i.e., 1D objects within a 2D reservoir).
3.1. Governing Equations
The flow in naturally fractured reservoirs is often described by the equations for nearly incompressible single-phase flow. We include gravity effects in the formulation because gravity can play an important role in the flow field evolution of a reservoir fluid with variable densities. We note that the following equations are assumed valid in both the damaged matrix and the fractures. From mass balance and single-phase fluid flow, the pressure equation iswhere [—] is the porosity, is the fluid density, and [Pa] is the fluid pressure. Compressibilities [] are denoted with subscripts for fluid and for rock, respectively. Moreover is the permeability and the fluid viscosity. In general, the permeability is an anisotropic 2nd-order tensor. However, only an isotropic permeability is used hereafter and in the implementation.
From the fluid pressure , the fluid velocity is calculated using Darcy’s law; that is,
The total mass balance equation derived above is separated into the matrix and the fracture domains, that is,where and are the flux transfer functions between the damaged matrix and the fractures. Superscripts and denote matrix and fracture quantities, respectively.
The velocities can subsequently be used in the heat transport equation. The heat transport equation is derived similarly to the continuity equation by the balancing the heat transport mechanisms and described asif local thermal equilibrium is assumed. Overlined properties denote volume averaged mean values for the porous medium (i.e., ). The heat transport equation is separated into matrix and fracture parts according to the same procedure as for the fluid pressure.
3.2. Fracture-Matrix Coupling
We apply a transfer function governing the mass and heat exchange between the two domains since flow in the damaged matrix is treated separately from flow in the fractures. The transfer function is treated as a source/sink term in the pressure and heat transport equations for damaged matrix and fracture, respectively, similar to classical well models (Peaceman [25]), that is,with being the mean total mobility of the fluid, defined as the fraction of permeability and viscosity. CI is the connectivity index between matrix and fracture that is grid dependent and defined based on the linear pressure distribution assumed within a grid cell intersected by a fracture [22]. The connectivity index is defined as the length fraction of fracture segment inside matrix cell divided by the average distance between matrix cell and fracture segment :The average distance can be calculated aswhere is the distance from the fracture within the matrix cell and the volume of the matrix cell. This allows properly accounting for the reduced influence of a fracture segment on a matrix cell if the fracture segment does not cross the matrix cell through its center. In many cases, (18) has to be evaluated by numerical integration. For rectangular grids, however, there exists an analytical solution [22] for fracture intersections horizontally, vertically, or on the diagonal to a grid cell. For enhanced efficiency, the analytical expressions are used in the present implementation. From the separated mass balance equations, it becomes immediately clear that the total flux between matrix and fracture has to be conserved: Obviously, is only nonzero in matrix cells that are actually intersected by at least one fracture segment.
3.3. Fracture Intersections
Fractures often intersect other fractures in naturally fractured reservoirs, which potentially significantly impact flow dynamics in the reservoir. Therefore we must also consider fracture-fracture coupling in the model. The additional transmissivity at a fracture intersection can be obtained similarly to the approach used in electrical engineering known as the star-delta transformation in circuits (Karimi-Fard et al. [26]). The additional fracture-fracture transmissivity can be calculated aswhere denotes the fracture aperture, the total mobility, and the numerical discretization spacing in the fracture. This approach can be generalized to more than two fractures intersecting in a single point on a fracture segment but is omitted here due to its rare occurrence.
3.4. Discretized EDFM Equations
The EDFM equations in two dimensions are discretized by the Finite Volume Method (FVM). Both the pressure equation and the heat transport equation are discretized by a two-point flux approximation scheme. The time derivatives in (12) and (15) are treated by an implicit time discretization. The advection term in the heat transport equation is treated by the QUICK scheme [27]. Although an explicit coupling also exists between the pressure and transport equations themselves, the current implementation uses a serial scheme to solve the coupled problem. Instead of assembling and solving one very large system for pressure and transport, the problem is divided into two parts. In a first step, the pressure system is assembled and solved. Using Darcy’s law, the fluid velocities can be calculated in an intermediate step. Once the fluid velocities are found, the transport system can be assembled and subsequently solved. In strongly coupled flow and transport problems, an iterative scheme must be used to capture any arising nonlinearities. In most cases, the flow and transport exhibit rather loose coupling in which only few iterations are needed to converge to the solution. It is restated that both equations are discretized by an unconditionally stable implicit time-discretization scheme. In case the problem shows nonlinear behavior, issues with nonconvergence might appear and place an indirect restriction for the time step. Nonetheless, much larger time steps are allowed in the implemented approach when compared to explicit schemes.
4. Results and Discussion
We present the results of three numerical experiments that provide insight into the effects of thermal stress on the stimulation of a geothermal reservoir. First we evaluate the influence of contrasts in the fracture permeability and the fracture-matrix permeability. We then investigate the influence of the fracture orientation on potential slip in combination with thermal stress. In a final numerical experiment, we simulate fluid injection into a complex fracture network over a prolonged period of time to evaluate the temporal evolution of thermal stress and failure potential. In all experiments, the pressure and/or temperature dependence of the fluid density and viscosity is taken into account. The underlying equation of state is given by Sun et al. [28] for density and Al-Shemmeri [29] for the viscosity of water. In this section, we combine results and their discussion to emphasize each experiment’s outcome. A more general discussion of the results as well as a conclusion is provided in Section 5.
4.1. Influence of Fracture and Matrix Permeability on Thermal Stress
Fracture permeability is one of the most important variables to be determined in EGS reservoirs to accurately predict flow in the reservoir. Here we evaluate to what degree the fracture permeability has an impact on thermal stress. To this end, we model the fluid injection in a well that is intersected by a single fracture as depicted in Figure 3.
Table 2 lists the physical parameters for the matrix and fracture used in this study. In this experiment, we investigate the results after a continuous injection of 30 days. In contrast to the previously presented temperature distribution in the well, we assume a homogeneous temperature distribution in the open hole section throughout the injection as the open hole section is short (100 m) compared to the length of the borehole. The injection temperature is variable with time and obtained by a least-squares fit of the data presented in Figure 1 for a depth of 5 km. The initial temperature is set at . The injection pressure is constant at 25 MPa. Gravity is neglected in this experiment. We vary the fracture permeability in five steps over a range within . Further we evaluate two matrix permeabilities and .
Figure 4 shows the temperature distribution and resulting thermal stress after 30 days of fluid injection for a matrix permeability of . The temperature profiles along the fracture show significant differences based on fracture permeability. For the two lowest fracture permeabilities (), the profiles are very similar and show temperature variations only within the first 10 meters of the fracture. At the injection point, the temperature dropped down to approximately . The resulting thermal stress is approximately −70 MPa. An intermediate fracture permeability of shows a further propagation of the thermal front within the 30 days. We observe temperature change on almost 60 m of the fracture resulting thermal stress perturbations along the same length. At an approximate distance of 30 m from the injection point, the thermal stress reaches −9 MPa. The two highest tested permeabilities () show distinctively different profiles compared to the lower permeability values. Here the temperature propagates the full length of the fracture with an almost linear profile. The temperature at the injection point however is greater than for the other permeabilities with approximately . Due to the temperature profiles, the thermal stress is significantly larger along the whole length of the fracture for these very permeable fractures.
We further evaluated the mean velocity in the fracture throughout the 30 days of injection. The result is shown in Figure 5. Overall the mean velocity varies over four orders of magnitude due to the wide range of fracture permeabilities. For all fracture permeabilities, an initial decline in the fracture velocity is visible. This is more pronounced for high fracture permeabilities () with a “stabilization time” of only a couple of days. Lower fracture permeabilities show less variability in the mean fracture velocity.
Figure 6 shows the temperature distribution and resulting thermal stress after 30 days of fluid injection for a matrix permeability of . The temperature profiles along the fracture show significantly less differences compared to higher matrix permeability. The three highest permeabilities show nearly identical behavior. The temperature front advanced as far as 25 m which is the maximum for a matrix permeability of . The two lower fracture permeabilities () show less temperature change after 30 days compared to the other permeabilities as well as the higher matrix permeability. All profiles exhibit a temperature at the injection point of .
A comparison of thermal stress distribution for the fracture permeability with both matrix permeabilities is shown in Figure 7. As already discussed earlier, the different propagation depths of the temperature and thus thermal stress are visible. Additionally the matrix effect is also shown. For the higher matrix permeability, the thermal stress front propagates as far as 6 m into the reservoir, and we observe a wedge-shaped thermal stress disturbance close to the fracture. In the case of the lower matrix permeability, the thermal stress front barely penetrates the matrix. The 1 MPa isoline is already at a distance of approximately 2 meters from the injection well. Close to the fracture a more pronounced thermal stress, alteration zone is visible. In contrast to higher matrix permeability this zone is relatively small.
(a)
(b)
Figure 8 shows that the pressure diffusion for the higher matrix permeability is much more homogeneous than for a lower matrix permeability, leading to a significant pressure gradient in the fracture. If the matrix permeability is very low, a high pressure zone develops around the fracture causing a smaller pressure gradient in the fracture. This is also observable in the mean fracture velocities. For high fracture permeabilities, the difference in fracture velocity increases one order magnitude with increasing matrix permeability. For low fracture permeabilities, there is no significant change in observed mean fracture velocity with respect to matrix permeability. The observed “stabilization time” in the fracture velocities is caused by the initial pore pressure diffusion through the matrix and fracture.
(a)
(b)
The differences in the pressure field and the fracture velocities can, however, not fully explain the observed thermal stress distributions. In order to explain these differences, we examine the heat transport in more detail. It is clear that the matrix-fracture interaction plays an important role in the thermal stress propagation. With increasing matrix and fracture permeabilities, also the interface permeability between matrix and fracture is increased allowing heat transfer by advection between matrix and fracture. Only by advection has the fluid in the fracture a chance to cool down the surrounding matrix. Once the surrounding matrix is cooled to a certain extent, the temperature front can also advance within the fracture. In case of a low matrix permeability, the main heat transfer mechanism between matrix and fracture is heat diffusion. Here heat diffusion from the matrix into the fracture dominates, leading to a heating of the fracture and thus significantly slowing down the thermal stress propagation. The second scenario is favorable in terms of sustainability of an geothermal reservoir. Here the heat is efficiently extracted from the matrix rather than prematurely cooling down the matrix surrounding the fractures.
4.2. Influence of Thermal Stress on Shear Failure Potential
In this section, we investigate the thermal stress influence on shear failure potential. As mentioned in the introduction and methods, shear failure in the reservoir occurs mainly on optimally oriented fractures. Here we investigate whether thermal stress leads to earlier onset of slip and how this is influenced by fracture orientation and model the fluid injection into a single fracture (Figure 9).
All parameters except for the matrix and fracture permeabilities are identical to the previous experiment and shown in Table 2. We apply a fracture permeability of with a matrix permeability of . The injection pressure is constant at 9 MPa. The minimum and maximum principal stresses are 45 MPa and 20 MPa, respectively. Stresses are oriented in alignment to the coordinate axis as shown in Figure 9. We vary the fracture orientation (measured in degree from the normal to the maximum principal stress) from to in steps of . The injection pressure was chosen so that only very close to optimally oriented fractures at are eligible for slip. In this experiment, we investigate the results after a continuous injection of 10 days. The injection temperature is variable with time and obtained as presented in the previous experiment. The initial temperature is set again at and we neglect gravity.
Figure 10 shows the result for a rotation of after 10 days of injection. It shows a pressure ellipse with its principal axis aligned with the fracture. The initial point source is not recognizable in the matrix pressure after 10 days of continuous injection. This can be explained by the fast pressure diffusion in the fracture and consequential pressure diffusion from the fracture to the matrix. The thermal stress is distributed symmetrically around the injection point in the fracture. The temperature front propagates roughly 10 m in each direction. The maximum cooldown and consequently maximum thermal stress are found at the injection point. Thermal stress in the matrix is propagated again from the fracture. The propagation direction is normal to the fracture leading to an elliptic thermal stress alteration zone surrounding the injection point. The rotation of the fracture barely influences the pressure and thermal stress distributions. Minor differences were observed but can be attributed to numerical effects.
(a)
(b)
We use the Mohr-Coulomb diagram to evaluate the failure potential due to thermal stress with respect to fracture orientation. Figure 11 shows the failure potential for each fracture orientation, where each point in the diagram is the mean normal and shear stress calculated over the whole fracture length. We show the effective stress modified solely by the fluid pressure as well as modified by thermal stress and fluid pressure combined. Figure 11 shows that, accounting only for fluid pressure in the effective normal stress, only fractures oriented at are able to slip. All other orientations do not meet the failure condition. If we further consider the effects of thermal stress, we observe a wide range of fracture orientations that fulfill the failure condition and are able to slip. The results after 10 days of injection show that fractures oriented up to from the optimum () are critically stressed. As Figure 11 only considers the average stress state in the fracture, we evaluate the failure potential in greater detail with a histogram. Figure 12 shows the percentage of fracture segments for each orientation that could demonstrate slip. Considering only fluid pressure effects, we find that the injection pressure of 9 MPa propagated as far as necessary to induce slip over the whole length of the optimally oriented fracture. Fractures oriented at from the optimum already show significantly less segments with failure potential (). Adding the thermal stress effect, we observe a broad spectrum of potential slip in up to of the fracture segments (at ). However, also fractures with an orientation further from the optimum show a significant ability to slip with ~20% at and ~5–10% for fracture oriented at from the optimum.
The observed asymmetry in the failure potential is of numerical origin because as the fracture is rotated, the number of matrix grid cells intersected by the fracture is changing. Consequently, this leads to small differences in the solution that propagate to the failure potential.
This experiment shows that thermal stress can greatly enhance the range of fracture orientations eligible for slip. This has a potentially large impact for the stimulation of geothermal reservoirs with a complex tectonic history and multiple fracture sets. If one or more fracture sets are nonoptimally oriented and unsuitable for hydraulic stimulation, they might still be suitable for thermal stimulation. However, due to different propagation speeds of the fluid pressure and thermal front, thermal stimulation might be challenging for short stimulation scenarios. The significantly slower propagation of the thermal front suggests that enhanced slip on close to optimally oriented fractures is unlikely. Nevertheless, thermal stresses can be expected to add to slip in the reservoir on nonoptimally oriented fractures and on fractures with highly heterogeneous frictional properties in the long run. The role of the different timescales of fluid pressure and thermal stress will be examined in more detail, the third numerical experiment.
4.3. Cold Water Fluid Injection into a Complex Fracture Network
In the third numerical experiment, we investigate the spatial and temporal evolution of potential failure due to fluid pressure and thermal stress. We model the fluid injection into a complex fracture network with a range of fracture orientations. The initial setup is shown in Figure 13. The domain represents a fractured reservoir at a depth of 5 km. We used the fracture network generator FracSim3D [30] to create the fracture network used in this study. It consists of a total of 310 large-scale fractures within a damaged rock matrix. The borehole is located in the middle of the domain with an open hole length of 50 m (cf. Figure 13).
All parameters except for the matrix and fracture permeabilities are identical to the previous experiments (cf. Table 2). We apply a fracture permeability of with a matrix permeability of . The injection pressure is constant at 35.1 MPa. The minimum and maximum principal stresses are 113 MPa and 175 MPa, respectively. The in situ pore pressure is assumed to be hydrostatic at 50 MPa. We neglect the hydrostatic gradient within the reservoir due to the relatively small vertical extent of the model. Stresses are oriented with rotation to the coordinate axis as shown in Figure 13. In this experiment, we investigate the results during a continuous injection of 30 days. The injection temperature is variable with time and obtained as discussed previously. The initial temperature is set to , and we include gravity effects.
An important addition to this experiment is a stepwise change in fracture permeability when a fracture segment reaches the failure condition [31]. That is, when the failure condition is reached, then the permeability adopts , where is a multiplication factor. For the simulations presented here, the enhancement factor is set to . The geological basis for stepwise change in permeability [31] rests with the strong aperture dependence of permeability (e.g., Nemčok et al. [32]), where small changes in aperture result in very large changes in permeability. The proposed model has been used to successfully describe the distribution of the induced seismicity in the Basel EGS site and for modeling fluid-driven aftershock sequences [10, 33].
Figure 14 shows the pressure distribution after 30 days of injection. Due to the heterogeneity of the fracture network, the pressure distribution is complex. As is also visible in Figure 13, the fracture density on the left side of the injection well is much higher than on the right side. This leads to more flow towards the left side of the reservoir as clearly indicated by the pressure field. The pressure distribution is dominated by the fracture network. After 30 days of injection, a pressure change of at least 1 MPa is measured in the whole domain. The maximum extent of the area with at least 30 MPa is approximately 30 m on the right side of the injection and up to 75 m on the left side of the well. A distinct anisotropy is visible between the horizontal and vertical extents of the high pressure zone. Vertically the pressure propagates slower compared to the horizontal direction and is caused by the main direction of the fracture set.
Figure 15 shows the thermal stress caused by the fluid injection into the reservoir. As predicted by the results of the previous experiments, the thermal stress alteration is concentrated to a region close to the injection well. Thermal stress evolves around the fractures intersecting the well as most of the fluid enters the reservoir here. The color scale in the figure starts at 0.65 MPa with the darkest blue. Everything below is neglected in the graphical representation and shown in the background color. This low cutoff shows the extent of the thermal stress alterations in the domain. The propagation of thermal stress in the matrix is rather homogeneous and mainly driven by heat diffusion. Additional thermal stress of 0.65 MPa occurs in a distance of 6–10 m from the borehole. Close to the borehole, thermal stress in the matrix reaches up to 20 MPa. The thermal stress in the fractures differs from the thermal stress in the matrix. Maximal thermal stress perturbation is measured at the fracture intersections at the borehole with up to 40 MPa additional stress. Due to the complex flow pattern, the thermal stress propagation in the fractures is equally strong. Some fractures show stress perturbations up to 10 m from the borehole, while, for others alterations, they do not reach 5 m into the reservoir. This shows that heat in the fractures is transported also by advection. A distinctive separation of the heat transfer mechanisms as previously performed for a single fracture is not possible here due to the interplay of the various effects.
We further evaluate failure potential in the reservoir due to the fluid injection. Figure 16 shows the distribution of potentially slipping failure segments in the reservoir at the end of the injection period. Each fracture segment is tested for potential failure based on (5). If the failure criterion is met, a star is shown in the figure to indicate potential slip on this segment. The color indicates the failure mechanism. Symbols in red indicate potential shear failure solely due to the fluid pressure. Blue symbols, on the other hand, signal failure due to the combined effects of increased fluid pressure and thermal stress. Pressure induced shear failure (or hydraulic failure) occurs mainly on the left side of the injection well and up to 65 m away from the well. Note that not all fractures are hydraulically stimulated as their orientation is not suitable for slip, although the fluid pressure inside the fractures is comparable (cf. Figure 14). Close to the well thermohydraulic failure is visible. The extent corresponds to the size of the thermal stress alteration zone. These fractures are not aligned optimally for slip and thus do not show potential failure with hydraulic stimulation alone. Once thermohydraulic stimulation is taken into account, the reduction in effective normal stress on these fractures is high enough to make them eligible for slip as indicated by the blue symbols in Figure 16.
The temporal evolution of the failure potential has been neglected so far. In Figure 17, we evaluate the simulated seismic events to gain more insight into the mechanisms at play and their evolution. An event is counted if the corresponding fracture segment is eligible for slip. So-called H-events indicate failure only due to fluid pressure. The TH-events are all events including both hydraulic failure and thermohydraulic failure. Initially, during the first days of injection, both event curves are close to each other. Later in the injection cycle, the difference between the curves grows larger, indicating the increased significance of the thermohydraulic failure events. Throughout the 30 days of fluid injection, the difference between - and TH-events is about 15%. With continued injection, this difference is expected to increase as the fluid pressure propagation slows down the higher the radial distance from the well and thermal stress propagation continues. Note that both curves share many of the sharp gradients and subsequent plateaus. This indicates that hydraulic failure dominates the shape of the measured event curves. Steep gradients indicate abrupt failure in many connected fracture segments possibly leading to microseismicity. Overall, the microseismicity is likely to increase due to combined thermohydraulic failure. However, large felt microseismic events are not a necessary consequence as there are no sharp gradients in the difference between the two curves (black dashed line in Figure 17) during the 30 days of injection into the reservoir.
This third numerical experiment shows that thermal stress can play and important role in shear stimulation of a fractured reservoir. Figure 16 showed that two distinct stress domination regimes can be defined. In one region, only hydraulic failure is observed, while, in the other region (more or less confined close to the well), thermohydraulic failure occurs. We showed that these regions act on different timescales. Fractures that are not optimally oriented for failure and are not prone to slip early in the injection cycle might be triggered due to the combined thermohydraulic stress changes later during the injection. This could help to explain the sometimes nonintuitive temporal evolution of microseismicity observed during fluid injection.
We can compare our results to typical changes in the static Coulomb stress transfer model that is widely used. Typical changes in Coulomb stress during fluid injection are observed in the range of −2 MPa to 1 MPa [34, 35]. In our experiments, we showed that thermally induced stress changes can easily exceed these changes in close proximity to the fractures and the injection point. This further emphasizes the importance of thermal stress during reservoir stimulation.
5. Conclusion
We developed and implemented a fast and efficient methodology for investigating the role of thermal stress in a geothermal reservoir. The numerical experiments presented in this study give insight into the role of thermal stress during hydraulic stimulation over short and long periods. We showed that thermal stresses can facilitate slip on nonoptimally oriented fractures and that thermal stress propagation is largely influenced by the hydraulic properties of both fracture and matrix. We found that thermal stress change is spatially more concentrated and generally propagates much slower than fluid pressure. Nevertheless thermal stress changes can exceed fluid pressure especially close to the well.
Some of our simulations showed that the combination of pore pressure and thermal stress can exceed the minimum principal stress, which would thermally induce tensile cracks. However, these tensile cracks are not self-propping, so hydraulic pressurization would be needed to facilitate flow. This is not considered in our current model but could be introduced as shown, for example, in Ghassemi [36]. Figure 6 uses a damage methodology for the matrix blocks surrounding the fracture. Including these effects would lead to an enhanced heat exchange between matrix and fracture in addition to that discussed previously. We also identified two different failure regimes that act on different timescales. The hydraulic failure regime acts on the timescale of days to weeks, whereas the thermohydraulic regimes acts on the scale of weeks to years. Future studies should aim to determine whether this can be observed in microseismic data from real-world fluid injection experiments.
Finally, we showed that thermal stress changes can be significantly larger than the injected fluid pressure as well as the slip-induced Coulomb stress change in some situations. This is especially important in long-term injection scenarios where the thermal stress changes become more significant with time and during the production phase of high flow rate geothermal systems.
Although we assumed water as the working fluid, it would be interesting to observe how the model behaves by implementing an equation of state for as the working fluid. Since is typically at supercritical conditions in a similar reservoir, we would expect similar hydraulic outcomes. However, specific heat capacity, density, and viscosity of supercritical differ significantly from those of water, which could affect the spatial evolution and heat flow, consequently influencing thermal stress propagation within the reservoir.
Our model using simplified mechanics provided important insight into the dominant thermoelastic effects and helped to identify some challenges and opportunities for future studies. More advanced models currently under development will consider both preexisting fractures as in the present work but also the generation of new fractures in response to the evolving stress state from both thermo- and hydraulic perturbations. Future models might also include fracture roughness and morphology and solving the full equilibrium equations to estimate aperture changes that influence permeability. Finally, including a method to estimate magnitude from computed slip [37], not currently modeled, is essential for a mechanistic assessment of seismic hazard associated with injection. Induced seismicity cannot be quantified in terms of magnitude as fracture slip in the fractures is not computed. These are all areas that we are currently pursuing.
We conclude that thermal stress changes in the reservoir should be incorporated in models that seek to fully understand the processes at play during fluid injection (covering initial stimulation as well as the production phase) in fractured reservoirs over long periods.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors appreciate the helpful discussions with B. Galvan and R. Sohrabi. This work was supported by the Swiss National Fond (SNF), for the financial support through Grant “NFP70: Energy Turnaround” under Project no. 153971.