• A hydro-mechanical model of effective stress changes from injection and production is calibrated using surface deformation from InSAR • Seismicity, located predominantly in the basement, is primarily induced by poroelastic stress changes from cold water re-injection • Similar stress changes trigger four times more earthquakes per unit volume in the basement than in the sedimentary cover
Supporting Information may be found in the online version of this article.
or due to thermo-or poroelastic stress changes (e.g., Goebel & Brodsky, 2018;Segall & Fitzgerald, 1998;Segall & Lu, 2015;Segall et al., 1994). The effect of pressure diffusion might, in principle, be identified from the migration pattern of seismicity (e.g., Shapiro et al., 1997). By contrast, poroelastic stresses due to a well operation are imparted almost instantaneously and should result in a spatially stationary pattern, which can reach large distances shortly after injection (Goebel et al., 2017;Zhai et al., 2019). In reality, disentangling the relative role of these mechanisms is often challenging. One reason is that the temporal evolution of seismicity also depends on the earthquake nucleation process, which may result in a time lag between fluid injection and the seismicity response (e.g., Alghannam & Juanes, 2020;Dieterich, 1994;Zhai et al., 2019).
Another reason is that field operations often involve a complex set of injecting and extracting wells (e.g., Hornbach et al., 2015) in addition to the regional tectonic stresses (Zhai & Shirzaei, 2018); assigning seismicity to particular sources of stresses is therefore often ambiguous in such a context.
Here, we study the case-example of seismicity induced at the Raft River geothermal field, which is particularly appropriate to gain insight into induced earthquake triggering mechanisms (Figure 1). The number of wells is relatively small, and the geological setting is simple and well documented (e.g., Bradford et al., 2013;Nash & Moore, 2012). In addition, a prominent geodetic signal has been previously observed using Synthetic Aperture Radar Interferometry, InSAR (Ali et al., 2018;Liu et al., 2018). Such measurements can indeed provide important constraints on the spatio-temporal evolution of the pressure field and fluid flow in the sub-surface reservoirs (Ali et al., 2018;Chaussard et al., 2014;Hoffmann et al., 2001). Finally, and most importantly for this study, local seismic monitoring has revealed seismicity that clearly correlated with the geothermal field (Figure 1a). A specific feature of the Raft River geothermal field is that hot water is produced from a ∼100-m thick reservoir at a depth of ∼1,500 m below surface, and the cold water is reinjected at a shallower depth. These injecting wells extend to the basement (∼1,550 m depth) but are open at depths greater than 500 m. The majority of the reinjected cold water flows into a "thief zone" consisting of ∼300 m thick permeable layer at 500 m depth (Ahmed et al., 1979;Bradford et al., 2013;Liu et al., 2018;Spencer, 1979).
Given that the production reservoir directly overlies the seismogenic basement, it is unlikely that pore pressure diffusion into the basement is responsible for the observed seismicity. If the zone of seismicity was hydraulically connected to the reservoir, fluid extraction should have decreased the pore pressure, leading to fault strengthening. As a result, the thermo-and poroelastic contributions from production and injection are perhaps the primary mechanisms behind the induced seismicity. Hereafter, we first give an overview of the setting of the Raft River geothermal field. We next present the methods used to model surface deformation and effective stress changes in the sub-surface. In the results section, we investigate the relative importance of injection and production on the spatial and temporal distribution of induced seismicity.
The Raft River geothermal site was identified in the 1950s, first developed as a United States Department of Energy (DOE) demonstration site in the 1970s, and significant production started by the end of 2007. The geology and geomechanical properties of the subsurface is relatively well documented in publicly available documents thanks to investigations conducted by the USGS and recent studies, and in relation to the recent EGS demonstration (Bradford et al., 2014;Nash & Moore, 2012, 2015;Yuan et al., 2020). The Miocene to Pliocene sedimentary cover is about 1.5-2 km thick and consists of sub-horizontal interbedded volcaniclastics, volcanics, sandstones, and siltstones. It overlies a proterozoic metamorphic basement consisting of quartzite, schist, and quartz monzonite. The area is part of the Basin-and-range zone of ∼East-West active extension. Historical seismicity rates were however very low prior to production, with normal-faulting focal mechanisms broadly consistent with E-W extension (Zandt et al., 1982). Well logs show steeply dipping fractures trending approximately N-S (Bradford et al., 2013;Nash & Moore, 2012).
Hot water is extracted at a depth of ∼1,500 m below surface from a ∼100 m thick quartzite layer through four wells (Figure 1b). One well, RRG-9 (Figure 1) was found to be poorly connected to the geothermal reservoir and selected for an Enhanced Geothermal System (EGS) stimulation. To avoid thermal drawdown, the cold water is injected through three wells and, as noted above, flows mostly into a shallower reservoir at a depth of ∼500 m. This reservoir consists of Late Miocene tuff, about 300 m thick.
In this study, we used the baseline mechanical, hydrological, and poroelastic parameters of the reservoir rocks determined in previous hydro-thermal-mechanical studies (Bradford et al., 2014;Liu et al., 2018;Yuan et al., 2020) (Table 1). The reservoir permeability and compressibility were then calibrated to the geodetic data (Text S1 in Supporting Information S1) to better constrain the pore pressure diffusion and poroelastic stress change models, respectively.
Significant production started only by the end of 2007 and ramped up gradually until 2010. The site was selected for an Enhanced Geothermal System (EGS) demonstration by DOE. The EGS demonstration was conducted at well RRG-9 over three phases between February 2012 and April 2014 with the goal of using hydraulic stimulation to improve the injectivity of the surrounding reservoir (flow rates are shown in Figure 3b). Each phase consisted of cold fluid injection followed by a shut-in period to assess the pressure falloff characteristics. The well has since been used for continuous commercial use.
The Lawrence Berkeley National Laboratory deployed a local seismic network consisting of 10 stations, which started operating in early 2010 (Figure S2). Their seismicity catalog, which is not de-clustered, shows seismicity clearly concentrated around the geothermal field and correlated with the surface deformation pattern measured from InSAR (Figure 1a). Seismicity is located primarily in the metamorphic basement immediately underlying the production reservoir. It started in late 2010, lagging the onset of large-scale production by about 2 years and pre-dating the EGS stimulation. We note that the magnitude of completeness appears to decrease from M C ∼ 0.5 prior to the EGS stimulation to M C ∼ 0 during and after the EGS stimulation. The seismicity cluster along a steeply dipping zone trending ∼N60°E coincides with the Narrows Fault zone, a basement structure which had been inferred from geophysical and hydrological studies (Bradford et al., 2013;Dolenc et al., 1981).
In our formulation, we account for the effects of pore pressure diffusion and its resultant poroelastic stress changes. We do not consider thermal effects given that the thermal contrasts are generally most important for injection of cold water, and the injection reservoir is very shallow (∼500 m depth) so that the poroelastic effects would dominate over the thermal effects as observed at the Salton Sea Geothermal Field (Barbour et al., 2016). In the production reservoir, fluid extraction should only cause minor thermal effects as heat is advected horizontally, and so, the horizontal temperature gradients should be small. This is corroborated by Liu et al.'s study (2018), where their coupled hydro-thermal-mechanical model of the Raft River geothermal field predicts that the poroelastic effects of fluid extraction and injection would dominate for the first 10 years, and the long-term effects of cooling in the injection reservoir would only emerge after 20 years of production.
In our workflow, we calculate reservoir pressure change resulting from fluid injection and extraction, surface displacements for comparison with the InSAR measurement, and stresses in the basement and overlying sedimentary cover for comparison with the observed seismicity. We employ analytical solutions as they
Injection reservoir Production reservoir Narrows fault Fluid density ρ, (kg/m 3 ) 998 998 Fluid viscosity, (Pa s) 8.9 × 10 -4 8.9 × 10 -4 Fluid compressibility, C f (Pa -1 ) 4.6 × 10 -10 4.6 × 10 -10 Skempton's coefficient, B 0.75 0.75 Permeability, k (m 2 ) 5.9 ± 1.5 × 10 -13 (4.7 × 10 -13 ) 5 ± 1.3 × 10 -13 (4 × 10 -13 ) Porosity 0.15 0.15 Compressibility, C m (Pa -1 ) 6.1 ± 0.3 × 10 -9 (5.1 × 10 -9 ) 6.3 ± 1 × 10
-10 (2.04 × 10 -9 ) Shear modulus, μ (GPa) 0.833 2.08 Thickness, (m) 300 100 Biot coefficient, α 0.31 0.31 Poisson's ratio, ν 0.2 0.2 Flow anisotropy factor, ξ 1.4 1.4 Flow anisotropy angle, 20° (slow diffusion N70°E) 20° (slow diffusion N70°E) Friction angle, (°) 31 Strike, (°) N60°E Dip, (°) 76 Sense Sinistral Note. The parameters which are calibrated using the surface deformation measurements are shown in bold. The values in parentheses indicate a priori values, and uncertainties at the 95% confidence level are estimated from bootstrap sampling (Figure S1).
minimize the number of model parameters while providing an excellent fit to the geodetically measured surface deformation data. In addition, they are computationally inexpensive, which can enable, in principle, scaling to larger systems involving thousands of wells.
Given that the local geology is to first order horizontally layered (Liu et al., 2018), we employ a 2D axisymmetric model for pore pressure diffusion. The production and injection reservoirs are assumed to be unconnected as suggested by geochemical studies (Ayling & Moore, 2013), and the fact that the permeability of the Salt Lake formation underlying the injection reservoir is significantly smaller than the reservoir permeability (see Figure S5 for an estimate of the pore pressure diffusion from the injection reservoir to the basement). In this formulation, considering a single well, the pore pressures at a given time, t, and a given location, x = (x 1 , x 2 ), from the well is,
where q(t') is the mass flow rate, c is the diffusivity, ρ 0 is the reference density, H is the thickness of the layer, κ = k/η, k is the permeability, and η is the fluid viscosity. The diffusivity c is k c
(2
where C m is the solid compressibility, n is the porosity, and C f = 4.6 × 10 -10 Pa -1 is the fluid compressibility. This solution is taken from Rudnicki (1986), and modified for anisotropy by transforming the radius r into an elliptical anisotropic radius r' defined as
where ξ = c fast /c slow is the anisotropy factor, is the anisotropy angle of the direction of low permeability, and
. This formulation follows the functional form analytically derived by Carcione and Gei (2009), and is included since the measured surface displacement (Figure 1a) suggests higher permeability in the ∼N-S direction. The pressure changes are calculated independently for each injection and production well and superposed for all wells within the same reservoir.
Given the pore pressure distribution in the production and injection reservoirs, the induced poroelastic stress and displacement fields can be approximated by gridding the pore pressure solution following Kuvshinov (2008). The solutions for the entire stress tensor and displacement vector can be found in the supplementary materials.
We evaluate the displacement at the surface (Text S1 in Supporting Information S1) and calibrate the model parameters using the geodetic measurements. We also evaluate the stress changes in the basement and the sedimentary cover to compare with seismicity. The stress changes are rotated to obtain normal and shear components on a strike-slip fault oriented at N60°E with 76° dip corresponding to the Narrows fault zone (Nash & Moore, 2012) and consistent with the orientation of the seismicity cluster. We calculate the Coulomb failure stress (CFS) change
where Δ is the change in shear stress on the fault, Δσ n is the change in normal stress, ΔP is the change in pore pressure, and μ = 0.6 is the friction coefficient.
We use a volume rate of 0.1 m 3 /s for injection wells, and 0.075 m 3 /s for production wells following Liu et al. (2018), with a ramp-up time of 2 years. We also include the effect of the EGS injection at well RRG-9 ST1 (Figure 1) using the data reported in the thesis by Bradford (2016), discretized into monthly rates. The best-fit parameters for permeability and compressibility are given in Table 1, please see Text S3 in Supporting Information S1 for additional details on the calibration procedure.
Figure 1a shows the spatial distribution of vertical displacement rate averaged from 2007 to 2011, which shows a clear signal of uplift around the injection wells. These wells are injecting at 500 m below the surface and thus have a stronger signature across a smaller footprint compared to the production which occurs at a deeper depth of ∼1.5 km below the surface. The elliptical bulb of uplift is well captured by our model (Figure 1b) estimated over the same time period of 2007-2011, which similarly shows highly localized uplift that is elongated in the NW-SE direction and generally captures the transition distance from uplift to subsidence thanks to the inclusion of permeability anisotropy. The isotropic model (Figure S3) predicts an E-W elongated uplift pattern due to the distribution of wells that fails to match this feature. The best-fitting permeability for the injection reservoir is 5.9 ± 1.5 × 10 -13 m 2 in the direction of slow diffusion, which is trending N70°E, and 8.3 ± 2 × 10 -13 m 2 in the direction of fast diffusion, which is trending N20°W. The N70°E permeability is within error from the value of 4.7 × 10 -13 m 2 used by Liu et al. (2018), while the N20°W permeability indicates faster diffusion than predicted in their study. While we are unable to directly determine the specific physical mechanism behind our observed anisotropy, we hypothesize that it arises from pre-existing fractures, which are oriented ∼N-S (Liu et al., 2018;Nash & Moore, 2012). This interpretation would however not be consistent with the northeast-oriented fracture zones and northeast-southwest direction of maximum permeability reported by (Yuan et al., 2020). The N20°W direction of fast diffusion and the anisotropy factor of 1.4 are also consistent with the analysis of Yuan et al., (2020), who found the permeability to be larger by a factor 1.2 to 2 along the NNE-SSW direction compared to the orthogonal direction. Note that some previous works on the Raft River site have modeled the injection wells at the same depth as the production reservoirs (e.g., Bradford et al., 2017). We have modeled this scenario in Figure S6 and have found the shallow injection scenario to represent a better fit (40% lower Wasserstein distance) to the observed surface displacements.
Our model also reproduces the temporal behavior shown in Figure 1d, with initial rapid uplift at location "A" before 2010, which then transitions to steady-state with no further uplift. The seismicity is located within the general vicinity of injection and production wells as shown in Figure 1 and begin in the basement in mid-2010 then continue in swarms of activity.
Figures 2a and 2b show the increase and decrease in reservoir pressures resulting from injection and production respectively. Note that the effect of the stimulation of well RRG-9 is not visible due to the very small injection volume (Figure 3b). The resultant poroelastic stress changes in the basement, shown for April 2014 in Figures 2c and 2d for injection and production respectively, indicate that the poroelastic contributions to Coulomb failure stress (CFS) changes affect fault stability in the same manner as the pore pressures, i.e., injection causes increases in CFS due to both pore pressure and poroelastic effects, and vice versa for the production. The poroelastic stress changes from production are smaller than those from injection, owing to the differences in compressibility and thickness. Specifically, the compressibility in the injection reservoir is 10× higher than in the production reservoir, and the injection reservoir is 3× as thick as the production reservoir. Given that these factors all affect the final poroelastic solutions linearly, a unit volume of injected fluid has 3.33× the effect of an equivalent extracted unit volume of fluid. As a result, the total ΔCFS in the basement (Figure 3a) suggests that injection is primarily responsible for the observed induced seismicity, which occurs almost exclusively in the regions of positive ΔCFS. To explore the sensitivity to the assumed fault plane orientation, we test different receiver fault orientations and find that the zone of increased ΔCFS around the injection and production wells is largely unchanged (Figure S4).
The time evolution of CFS at the location of maximum uplift (point A) and the edge of the seismicity cloud (point B) suggests that the eastern edge of the seismicity is dominated entirely by stresses from the injection, while stresses on the western edge of the seismicity are inhibited significantly by the production (Figure 3c). Note that the pressure evolution at point B, located near well RRG-9, shows a very short-lived pressure increase. This figure also shows that although the simulation assumes a simple ramp to constant flow rate from 2008 to 2010, the zone of Coulomb stress increases due to poroelastic loading of the basement keeps increasing until the end of the simulation as the zone of high pore pressure in the shallow reservoir expands away from the injection wells. This can explain the sustained seismic activity. Additionally, the ΔCFS at "A" as shown in Figure 3c indicates that there may be a critical Coulomb stress increase of approximately 50 kPa required to trigger the seismicity. That would represent the initial strength excess, that is, the initial distance from the failure criterion, within the Narrows fault zone. To further investigate this, we calculate the ΔCFS at the specific origin time and hypocentral location of the observed earthquakes (Figure 4a). In this calculation, we assume the same hypocentral depth of 500 m for events in the basement and -500 m for events in the overlying sedimentary cover. We can see that nucleation stresses range from 0 to 70 kPa for events in the basement and the overlying sedimentary cover, and the histograms shown in Figures 4b and 4c show a possibly bimodal distribution of nucleation stress for the basement with a main mode at LI ET AL. 10.1029/2021GL095108 8 of 12
19 kPa. It is interesting that the primary modes of the basement and sedimentary cover earthquake ΔCFS histograms are actually quite similar despite the fact that very few earthquakes occurred in the sedimentary cover. According to our analysis, a given volume of basement rocks submitted to a given ΔCFS produces four times more earthquakes than the same volume of sedimentary rocks submitted to the same ΔCFS. This is because the earthquake productivity, representing the density of possible nucleation points of earthquake with magnitude larger than the detection threshold, must be four times larger in the basement than in the sedimentary cover (Figure 4d).
Given the magnitude of the pore pressure changes in the production reservoir, which are about 100 times larger than the magnitude of the poroelastic stress changes, we believe it is justified to consider the production reservoir to be uncoupled from the basement, otherwise we would not expect any seismicity whatsoever since the pore pressure reduction from extraction would completely dominate the state of stress in the basement. We account for the hydraulic stimulation at RRG-9-ST1 as an additional injection well located in the production reservoir, although we note that its corresponding pore pressure contributions are small relative to the standard production wells given the low injection volumes (Figure 3b). This can be seen in the stress changes at point "B" (Figure 3c), which is located close to well RRG-9 but only experiences small stress changes on the order of ∼1 kPa as a result of the stimulation. Considering the spatial distribution of seismicity prior to the first stimulation phase, we see that the majority of the central seismicity cloud is already activated before March 2012 even with the lower magnitude of completeness (Figure S2), and so, the hydraulic stimulations likely did not significantly affect the hydraulic connection between the basement and the production reservoir. This is supported by the observation that we do not generally see spatial diffusion of earthquakes away from well RRG-9 after the three stimulation events (Figure S2c).
Together with previous studies which have adopted a similar strategy (Juanes et al., 2016;Shirzaei et al., 2016), our study demonstrates the of combining observation of surface deformation and simple fluid flow and geomechanical modeling to analyze induced seismicity due to injection and extraction of fluids from the subsurface. We acknowledge that our evaluation of uncertainties incorporates only the sensitivity of our analytical solution and does not account for potential sources of correlated error in the InSAR measurements (e.g., tropospheric corrections, ionospheric effects) or epistemic uncertainties relating to our assumption of homogeneous reservoir properties in a homogeneous elastic half-space, the depth of injection relative to production, and no vertical diffusion of pore-pressure. In particular, the shear modulus at the shallow depth of the injection layer is lower than that of the basement. Making it the same as half-space may bias our estimate of the effect of production. Nevertheless, we believe that our main findings, at least qualitatively, are not dependent on these uncertainties, given the similarity of our modeled and measured surface displacements with previous works (Ali et al., 2018;Liu et al., 2018), as well as our ability to resolve the permeability anisotropy suggested by Yuan et al. (2020).
The spatial distribution of ΔCFS contributions (Figures 2c and 2d) in the basement strongly suggests that poroelastic stress changes from shallow injection are responsible for the timing and location of observed induced seismicity. The pore pressures from injection are not connected to the basement, and the poroelastic stress changes resulting from fluid extraction are smaller in magnitude compared to the injection. Figure 3a, which shows the ΔCFS in the basement, indicates that all earthquake events occur in the region of positive slip potential, suggesting that while the magnitude of the contribution from production is small, it nevertheless has an inhibiting effect on fault slip potential, for example to the north of well RRG-9 ST1. In terms of the distribution of ΔCFS at the time and location of earthquake nucleations (Figure 4), we also find that the basement has an initial strength excess of approximately 15-20 kPa and an excess of 20-25 kPa in the overlying sedimentary cover, which is required to activate the fault, similar to the observations at Groningen by Smith (2019). This strength excess may be considered analogous to a distance to failure as defined in Rate-and-State formulations (Dieterich, 1994;Heimisson, 2019). Although minimal further uplift of the ground surface over the geothermal field occurred past 2010, Coulomb stress continued to increase slightly, and seismicity continued in the basement. This is consistent with our modeling results because the competing effects of injection and production contribute differently to surface deformation and to basement deformation owing to the different functional forms of the displacements and stress components.
In terms of the timing of seismicity in relation to the DOE hydraulic fracturing project at well RRG-9 ST1, our results show that the seismicity begins in late 2010 whereas the first phase at well RRG-9 ST1 begins in early 2012. This, in conjunction with the minimal ΔCFS at "B" (located close to well RRG-9 ST1) in Figure 3c, show that while the project did indeed trigger some induced events (Bradford, 2016), the timing and location nevertheless correspond to the larger-scale spatial distribution of ΔCFS shown in Figure 3a.
Overall, our results show that observed spatial and temporal distribution of surface deformation and induced seismicity at the Raft River geothermal site may be reasonably explained by simple analytical solutions for axisymmetric pore pressure diffusion and its resultant poroelastic stress changes. Our model uses only a simplified 1D geological model and injection and production rates as inputs, where geomechanical parameters are constrained by fitting the modeled surface displacement to the true surface displacement as observed by InSAR. This simple framework is scalable and easily calibrated given the small number of parameters and can be readily applied to investigate large multi-well systems involving combined injection and production such as Oklahoma and Texas (Langenbruch & Zoback, 2016;Walter et al., 2018;Zhai & Shirzaei, 2018), where existing studies are primarily focused on the effects of fluid injection. Our framework may be used to inform geothermal and carbon storage strategies, where the seismicity in the basement and elsewhere could in principle be controlled by regulating the injection rates (Birkholzer et al., 2012;Cihan et al., 2015;Kwiatek et al., 2019). We conjecture that, using the modeling workflow presented in this study, the location and time-evolution of seismicity induced by fluid extraction could in principle be controlled by adjusting well flow rates.