The spatial footprint of injection wells in a global compilation of induced earthquake sequences

See allHide authors and affiliations

Science  31 Aug 2018:
Vol. 361, Issue 6405, pp. 899-904
DOI: 10.1126/science.aat5449

Seismic limits for hard and soft rock

Induced earthquakes from oil, gas, and geothermal energy exploration projects can damage infrastructure and concern the public. However, it remains unclear how far away from an injection site an earthquake can still be triggered. Goebel and Brodsky looked at 18 different earthquake-producing injection sites around the world to address this issue. Injecting fluid into softer layers increased the range for seismic hazard, whereas harder basement rock better confined the fluid. These findings should be considered when regulating and managing projects with the potential to induce seismicity.

Science, this issue p. 899


Fluid injection can cause extensive earthquake activity, sometimes at unexpectedly large distances. Appropriately mitigating associated seismic hazards requires a better understanding of the zone of influence of injection. We analyze spatial seismicity decay in a global dataset of 18 induced cases with clear association between isolated wells and earthquakes. We distinguish two populations. The first is characterized by near-well seismicity density plateaus and abrupt decay, dominated by square-root space-time migration and pressure diffusion. Injection at these sites occurs within the crystalline basement. The second population exhibits larger spatial footprints and magnitudes, as well as a power law–like, steady spatial decay over more than 10 kilometers, potentially caused by poroelastic effects. Far-reaching spatial effects during injection may increase event magnitudes and seismic hazard beyond expectations based on purely pressure-driven seismicity.

Human-induced seismicity close to geothermal, hydraulic fracturing, and wastewater disposal wells is a source of substantial seismic hazard. Such injection activity has led to many moderate-magnitude earthquakes and an exceptional increase in earthquake rates in parts of North America and Central Europe after ~2006 (13). The hazard from injection-induced earthquakes is particularly difficult to manage because the earthquakes frequently occur at large distances (>10 km) from the targeted injection zones (46).

A better understanding of the driving mechanisms of injection-induced seismicity is vital for improving seismic hazard assessment and mitigation. Traditionally, such hazard assessment has concentrated on pore-pressure increase in a volume hydraulically connected to the injection wells (7, 8). Pore-pressure increase, which we here call the direct pressure effect, is thought to reduce the normal load on locked faults, thereby allowing sliding to occur.

As induced seismicity cases and studies have proliferated, the importance of additional mechanisms such as elastic and fully coupled poroelastic stresses have become increasingly clear (911). The pore-pressure increase in the injection zone is expected to load the surrounding rock matrix and result in a fully coupled fluid-solid stress field. Elasticity is an effective means of transmitting forces to great distances, and therefore the fully coupled poroelastic stress field can extend well beyond the fluid pressure increase in the hydraulically connected region. For example, large-scale, field-wide injection can perturb faults and induce earthquakes more than 30 km away (11). In controlled injection experiments that resolve coupled aseismic and seismic processes during fluid injection, induced earthquakes are absent within the pressurized zone but occur as a result of elastic stress changes in the surrounding rock volume (12). These observations require a reexamination of the controls on induced seismicity and suggest that the spatial reach of the earthquakes may be a useful discriminant of triggering mechanisms. The observations also suggest that induced earthquakes may extend farther than previously thought [i.e., at least to 30 km from wells (4, 11)], and empirically measuring the spatial extent is necessary to develop appropriate mitigation and regulatory approaches.

Here we examine the distance of induced earthquakes from injection wells to understand the mechanical controls on the spatial extent of induced seismicity. We start by quantifying the shape of spatial decay and find two groups with distinct decay patterns. We then analyze migratory behavior within each group, as well as the relationship between operational parameters and the resultant spatial pattern of earthquakes. We examine differences in observed maximum magnitudes, which appear to be linked to the type of spatial pattern. Lastly, we discuss possible physical mechanisms that control the distinct behavior within the two groups, with a focus on potential differences in poroelastic properties between sedimentary and basement units.

We compiled measurements of the spatial decay of induced earthquakes from causative wells worldwide (1323). We focused on isolated injection sites with well-recorded induced sequences, which are clearly connected to a single geothermal, scientific, or waste-fluid injection well (fig. S1; section S2 describes two exceptions with more than one well). This strategy eliminates some major regions, such as Oklahoma, where effects of field-wide injection are difficult to unravel. All of the data used are in the public domain (24). The initial 27 selected injection sites are predominantly located in intraplate regions with low seismic activity throughout the United States, central Europe, and Australia (Fig. 1A). We included one well from the northwestern edge of the Geysers geothermal field because of its isolated location and several-kilometer distance to other active injection wells during the analyzed time period from 2007 to 2009 (16).

Fig. 1 The spatial decay of induced sequences can be classified as abrupt or steady.

(A) Map of the injection sites (blue triangles) in the United States, Australia, and Europe (fig. S1). (B) Seismic density of all studied induced sequences, normalized by number of events above completeness. We show cases with abrupt decay in shades of blue and steady decay in shades of red. We shifted abrupt decays vertically (translucent blue dots show one example of the original vertical position; the blue arrow shows the shift) to allow for easier visualization. (C) Theoretical expectation of spatial density fall-off using Eq. 2; the blue curve shows the abrupt decay of pressure-dominated sites, and the red curve shows a power law decay for sites with strong poroelastic coupling. (D) Merged densities above rc for sequences with steady decay (black markers) and power law fit (red dashed line) with ~r−1.8, which is more gradual than spatial aftershock decay in California (gray markers and black dashed line).

Before determining the spatial seismicity decay from wells, we assessed the quality of the seismic record, focusing on spatial and magnitude information. We first selected events above the magnitude of completeness (Mc), determined by minimizing the misfit between the observed distribution and the Gutenberg-Richter relation (fig. S2). The datasets for Soultz-Sous-Foretz, France, and Fenton Hill, New Mexico, do not include magnitude information. For the Soultz site, magnitude estimates are available from a surface array for a subset of events, showing consistent results between borehole and magnitude-corrected surface catalogs. For the Fenton-Hill site, the entire seismic record was used.

In a second quality-control step, we tested whether the observed distributions exhibit significant deviations from random Gaussian location uncertainty and show significant spatial clustering close to wells at rates above the background activity, using a two-sided Kolmogorov-Smirnov test (fig. S4). The quality-control steps eliminated nine induced cases, leaving 18 sequences for further analysis.

For the high-quality sequences, we computed two-dimensional (2D) distances between wells and earthquakes at the average depth of the injection interval, taking into account the well trajectories. Relative horizontal location uncertainties of seismic events are on the order of tens of meters, whereas absolute location uncertainty ranges from 100 to 500 m (section S4), and average vertical location errors are more than 1 km. Because the vertical uncertainties are large, we initially focused on 2D distances and later compared the results with 3D distances. Seismicity distance fall-off from wells was calculated from areal densities determined by nearest-neighbor binning of the k closest events with a moving window of k/2, resulting in density estimates, ρ = kr, where Δr is the distance between the first and kth event.

The seismicity density distributions can generally be described by two types of spatial decay: (i) sequences with an extended plateau close to the well, followed by an abrupt decay within less than 1 km, and (ii) sequences with steady, power law–like decay out to distances of more than 10 km (Fig. 1B). The shapes of the two types of spatial decay remain consistent across a range of spatial scales and durations of injection operations, spanning hours (e.g., Rustrel) to years (e.g., Paradox), both for the complete datasets and for the subsets of events from individual sites (fig. S8). The consistency in the shapes of spatial decay suggests an underlying time-invariant process. We distinguished the two types of distance fall-offs quantitatively by fitting with the functional formEmbedded Image(1)where ρ is seismic density, ρ0 describes the short-distance density plateau, rc is the corner distance, and γ is an exponent describing the abruptness of the decay at larger distances. We chose this functional form to capture both differences in near-well density plateaus and the abruptness of distance fall-offs. We determined the three free parameters by fixing ρ0 using the average density of the first five sample points and then inverted for rc and γ using a maximum likelihood approach that assumes Poissonian uncertainty in each distance bin (24). The distribution of the resulting decay exponents exhibits a natural break between γ = 3.1 and 4.3, separating sequences into steady decay with γ = 1.5 to 3.1 and abrupt decay with γ = 4.3 to 5.9 (fig. S5). The difference between these two populations exceeds the maximum expected uncertainty of ±0.34 based on the 95% confidence intervals of the maximum likelihood fits. In addition to having larger γ values than steady-decay sites, sites with abrupt decay are characterized by corner distances that are substantially closer to the maximum extent of the sequence (fig. S6).

Several observations indicate that the result is not an artifact of specific instrumental, statistical, or measurement practices: (i) The maximum distance of the earthquakes is smaller than array aperture in all cases, so that the spatial decay is not a product of limited array extent, which may truncate the data. We further tested the influence of array geometry by computing Mc as a function of distance from wells and found no systematic bias. (ii) Where a single site allows for comparison between instrumentation (i.e., surface versus borehole arrays), the results are identical. Specifically, the spatial decay at the Soultz site shows abrupt fall-offs in both the borehole and the wide-aperture surface array catalogs when corrected for Mc. (iii) In addition to fitting Eq. 1, the two types of decay can be identified by using pure power laws. A power law is well fit over a large distance range for sites with steady decay, whereas abrupt-decay sites show power law–like behavior only over a limited range of distances that are near the maximum distance of the datasets (fig. S6). (iv) Lastly, stacking seismicity density estimates on the basis of 2D or 3D distance leads to a visible separation of steady and abrupt decay (Fig. 1B and fig. S11).

To connect the spatial decay to perturbing stresses, we considered induced pore-pressure changes and poroelastic coupling. We expressed the observed seismicity fall-off within a Coulomb failure framework as the product of stress perturbation and the number of faults sufficiently close to failure (25)Embedded Image(2)where ρ is seismicity density, Δσ is the induced stress perturbation, and Nfl/(2πrΔr) is the density of faults per area. To build intuition for the influence of changes in Δσ, we first assumed that the fault density is constant and investigated the distance fall-off expected for the direct pressure effect and poroelastic coupling. For the direct pressure case, we can write the distance fall-off of the first term in Eq. 2 in a vertically confined reservoir as Δσ = ΔPp ~ W(u), for which ΔPp is pore pressure, W(u) is the exponential integral, and the argument u = r2/(4Dt); r is distance, t is time from injection, and D is hydraulic diffusivity (26). The corresponding distance fall-off is substantially faster than poroelastic stress decay in the far field of injection, which decays as 1/r2 (26). The expected distance decay based on pure pressure and poroelastic models approximately matches the shapes of observed density fall-offs for abrupt and steady decays (Fig. 1, B and C).

Although the similarity between the simple models and the observed spatial decay is compelling, additional processes are expected to contribute. Such processes include the coupling between pressure and permeability, as well as heterogeneity in the permeability structure in the presence of faults (8, 14, 27). In addition, elastic stress transfer and event-event interaction of both seismic and aseismic ruptures can increase the extent of induced seismicity sequences (12).

We further analyzed sequences with extended, power law–like decay by fitting the joined densities above corner distances (rc) of individual sequences (Fig. 1D). The corresponding linear least-square fit of the log-transformed data has a value of γ = 1.8. This value is substantially smaller (i.e., more gradual) than the spatial decay of aftershocks from mainshock epicenters in California, for which γ = 2.4 (25) (Fig. 1D). The latter value was determined by identifying mainshock-aftershock clusters by using a nearest-neighbor distance in a space-time-magnitude domain (28). This difference in γ exponents indicates that spatial seismicity decay contains information on the distinctive forcing stresses specific to injection-induced sequences. The stresses from a mainshock that produce aftershocks cannot explain the data, so we must invoke additional processes such as fluid migration and poroelasticity. The smaller γ value for induced sequences points to a more gradual spatial relaxation of the underlying stress field.

To test the influence of location uncertainty on the shape and extent of spatial decay, we performed Monte Carlo simulations of random power law–distributed data convolved with random normal distributions with standard deviations that correspond to the observed location uncertainty (section S4). Our simulations reveal that decay exponents can in principle be inflated owing to location uncertainty. However, our location uncertainties are small enough that the effect is negligible for this study (fig. S9). We can robustly identify sites with abrupt decay in the 2D data because both relative and absolute uncertainties are smaller than the determined rc values (fig. S10).

We next investigated potentially systematic migratory behavior of seismicity in each sequence. Migration is particularly interesting because the most commonly invoked mechanism for induced seismicity—direct pore-pressure diffusion from a nearly instantaneous injection—has a well-understood square-root dependence of distance on time (9). Other cases worth considering are no migration, as might be expected for a rapidly applied elastic stress perturbation, and linear migration. These three models were fitted to the seismicity envelope using a maximum-likelihood method, and the preferred model was chosen on the basis of a Bayesian information criterion (Fig. 2 and section S3).

Fig. 2 Sequences with abrupt spatial decay are dominated by square-root migration, an indication of pressure diffusion.

Examples of (A) linear and (B) square-root migration. Gray markers show event time and distance from injection; white markers show the 95th percentile of distances in specific time bins. (C) Histogram of the number of abrupt- and steady- decay sites with square-root (blue), linear (orange), or no (gray) migration.

Square-root migration is most common at abruptly decaying sites, consistent with direct pore-pressure effects (Fig. 2). Cases with steady decay, on the other hand, show linear migration (five sites; Fig. 2C) or a lack of migration (six sites). The latter may be characteristic for poroelastic stress–dominated sequences during which earthquakes are triggered beyond the pressure-dominated region even shortly after injection starts, leading to a breakdown of square-root migration (1). Processes that may contribute to linear migration are related to elastic stress transfer, including event-event triggering, slow slip, and aseismic slip (12, 18, 29). Additional mechanisms that likely affect seismicity migration include fracture creation, changes in permeability structure, and thermally induced stresses. These mechanisms are expected to be most pronounced during geothermal injection and close to injection wells (16).

Distance decay and space-time migration highlight two distinct groups of induced seismicity sequences, consistent with direct pressure effects on the one hand and far-field elastic stresses on the other hand. To better understand the physical controls on this distinct behavior, we investigated operational and geologic parameters at all injection sites. We evaluated the ability of injection volume, average flow rate, peak well head pressure, injectivity (which is approximated by flow rate divided by well head pressure), injection depth, distance to basement (i.e., vertical distance to the crystalline basement rock), reservoir temperature, stress state, and injection duration to discriminate the two populations defined on the basis of spatial decay (fig. S12) (24). The most important governing parameter for the separation of abrupt and steady decay is distance to basement: All sites with abrupt decay were injecting below the upper basement surface, and sites with steady decay (except for Basel) were located in sedimentary rocks above the basement (Fig. 3).

Fig. 3 Above-basement injection commonly results in larger spatial footprints and a higher probability of inducing larger-magnitude earthquakes.

(A) The spatial decay, separated into abrupt (blue) and steady (orange), is controlled by the distance between injection and crystalline basement. (B) Maximum magnitude of each sequence as a function of total injected volume for steady (orange) and abrupt (blue) decays. The dashed line is the theoretically expected maximum moment based on (35), for which G is shear modulus and V is total injected volume. (For Fenton Hill, we report the largest recorded magnitude during all stimulations. No stimulation-specific fluid volume is available for Raft River.)

This difference in behavior based on lithology at injection depth can potentially be explained by a difference in the ability of the rocks to transmit fluid stresses into the solid material. The Biot-Willis α coefficient (26) captures the poroelastic coupling so that larger α values correspond to more effective couplingEmbedded Image(3)where K is the bulk modulus and Ks is the modulus of the solid, so that α is expected to be close to 0 for stiff, low-porosity rocks in the basement and close to 1 for soft sediments (30). Experimental results confirm that bulk strains in low-porosity rocks are less sensitive to fluid pressure changes (31). In other words, pressure perturbations within the larger, interconnected pore space in sedimentary rocks are much more efficient in changing bulk stresses than pressure changes in smaller, isolated pores in crystalline rock. Sedimentary units may also have higher permeability, which influences the spatial extent of injection but is less consequential for the observed separation into steady and abrupt decay.

A notable exception to the lithology-controlled seismic response is the Basel injection site. At this site, steady spatial decay of earthquakes up to local magnitude (ML) 3.4 was triggered by injection into the crystalline basement. A potential explanation for this observation is injection directly into a highly fractured fault damage zone, as indicated by seismicity and focal mechanism analysis (13). Such damage zones likely exhibit higher α values, even within the crystalline basement, but will rarely be encountered unless specifically targeted during injection because of their hydrogeologic properties. The 1967 Rocky Mountain Arsenal sequence is arguably in a similar category given the reported extent of the seismicity, but it is not included in this analysis because of the data quality issues discussed earlier.

In addition to differences in local geology, we also observed variations in maximum event magnitudes between the two populations of induced earthquakes. Previous studies showed maximum magnitudes above M4 for earthquakes induced by hydraulic fracturing and above M5 for earthquakes induced by geothermal operations or fluid disposal (1, 2, 15, 32). Peak magnitudes for a given injection volume are generally larger at sites with steady decay, which also produced the largest event in our dataset, with a magnitude of ML 4.7 (Fig. 3B). The difference in maximum magnitudes is expected if more extensive spatial footprints at sites with steady decay increase the probability of encountering larger faults (33). Even for injection into a single well, poroelastic stresses may promote the activation of such distant, large faults that are close to failure without requiring a direct hydraulic connection. This effect is further enhanced by closely spaced, high-rate injection wells, which act as a finite source with stress decay out to ~2 to 3 source dimensions (11).

The spatial decay of induced sequences can be expressed by the product of stress perturbation and available prestressed faults that fail because of this perturbation (Fig. 4). At sites with high poroelastic coupling, we expect elastic stresses to dominate earthquake triggering in the far field of injection, showing a scale-invariant spatial decay of ~r–2, whereas direct pressure effects dominate close to injection wells. If we assume that the term for fault availability in Eq. 2 is expressed by Embedded Image, where df is the fractal dimension of the fault network and D is the geometric dimension of the density measurement (here, D = 2), we find that df = 2.2. Thus, fault density increases gradually as r0.2 with distance, so that for pressure-dominated sequences, the expected point of peak seismicity occurs at some distance from the well (Fig. 4B). The model proposed here captures large-scale, average trends in seismicity decay with distance from injection wells. Nevertheless, additional processes such as coupling between changes in pore pressure and permeability, as well as variations in regional fault network geometry, may also contribute to the functional form of the decay.

Fig. 4 The probability of inducing an earthquake at distance r is controlled by fault availability and amplitude of stress perturbation.

(A) Schematic representation of injection operation, footprint of poroelastic response (blue and red ellipses), and fault network (gray lines). (B) Earthquake probability (in events per area) as a function of distance from the injection well for pressure-dominated triggering. (C) Same as (B) in a coupled system with elastic stress dominance in the far field. [Both the x and y axes in (B) and (C) are logarithmic.]

Our global compilation of fluid injection–induced seismicity allows for a better understanding of the maximum earthquake-triggering distance from an injection site. We can describe the shape of spatial seismicity decay as either abrupt or steady. Steady decay shows power law–like behavior to distances of more than 10 km and is more gradual than spatial aftershock decay. Sites with abrupt decay are dominated by square-root space-time migration, which is consistent with pressure diffusion. Abrupt decay is limited to sites where injection is within the crystalline basement, whereas steady decay primarily occurs above the basement. The maximum magnitude is larger for sites with steady decay owing to the greater probability of activating bigger faults within the extended spatial footprint of the injection wells.

Previous strategies to mitigate induced seismicity encouraged injecting in sedimentary units instead of directly into the basement (34). However, our results suggest that injection into sedimentary rocks leads to more distant and larger earthquakes for a given volume of injection, perhaps owing to more efficient pressure and stress transmission. The larger spatial footprints of above-basement injection may be responsible for the extensive seismogenic response in some areas, such as Alberta and Oklahoma (1, 4). Nevertheless, injection into the basement also poses a source of seismic hazard. Large earthquakes, such as the 1967 Rocky Mountain Arsenal ML 5.3 (8) or the recent South Korean Mw 5.4 events (32), can occur if fluid is directly injected into a basement fault that is either specifically targeted or encountered by chance. The key result from this analysis is that injection in sedimentary units immediately overlying basement rocks is more likely to encounter a large fault by chance because of the larger spatial footprint. Far-reaching poroelastic effects complicate the assessment of distance-based induced seismic hazard and should be included in mitigation strategies.

Supplementary Materials

Materials and Methods

Supplementary Text

Figs. S1 to S12

Table S1

References (3684)

Data S1

References and Notes

  1. Materials and methods are available as supplementary materials.
Acknowledgments: We thank P. Ogwari, H. De Shon, P. Martinez-Garcon, J. Ritter, M. Bohnhoff, L. De Barros, S. Wiemer, T. Kraft, J. Albaric, V. Oye, S. White, J. Riffault, and M. Fehler for their help with accessing the data and H. Shaddox, T. Lay, and S. Siman-Tov for discussions. Funding: We gratefully acknowledge funding from the U.S. Department of Energy, grant DE-SC0015539. Author contributions: T.H.W.G. and E.E.B. contributed equally to this work. Competing interests: The authors are not aware of any competing interests. Data and materials availability: Seismicity and injection data are available in the public domain (the supplementary materials contain references for each site).
View Abstract

Stay Connected to Science

Navigate This Article