Hydraulic fracturing volume is associated with induced earthquake productivity in the Duvernay play

See allHide authors and affiliations

Science  19 Jan 2018:
Vol. 359, Issue 6373, pp. 304-308
DOI: 10.1126/science.aao0159

Seismicity curbed by lowering volume

Determining why hydraulic fracturing (also known as fracking) triggered earthquakes in the Duvernay Formation in Canada is important for future hazard mitigation. Schultz et al. found that injection volume was the key operational parameter correlated with induced earthquakes in the Duvernay. However, geological factors also played a considerable role in determining whether a large injection volume would trigger earthquakes. These findings provide a framework that may lead to better forecasting of induced seismicity.

Science, this issue p. 304


A sharp increase in the frequency of earthquakes near Fox Creek, Alberta, began in December 2013 in response to hydraulic fracturing. Using a hydraulic fracturing database, we explore relationships between injection parameters and seismicity response. We show that induced earthquakes are associated with completions that used larger injection volumes (104 to 105 cubic meters) and that seismic productivity scales linearly with injection volume. Injection pressure and rate have an insignificant association with seismic response. Further findings suggest that geological factors play a prominent role in seismic productivity, as evidenced by spatial correlations. Together, volume and geological factors account for ~96% of the variability in the induced earthquake rate near Fox Creek. This result is quantified by a seismogenic index–modified frequency-magnitude distribution, providing a framework to forecast induced seismicity.

Subsurface injection of fluid may induce earthquakes (1) through anthropogenic alteration of crustal stresses (2). In the case of hydraulic fracturing (HF), high-pressure injection of fluid intended to increase the permeability of tight shales has been known to trigger earthquakes (3), some of which are large enough to be recorded or felt regionally (48). Within the Western Canada Sedimentary Basin, the recent increase in seismicity has been largely attributed to HF (9). Moreover, earthquakes in the Duvernay Formation (10) have been among the largest-magnitude events caused by HF completions globally (11, 12). These events have appreciably increased the seismic hazard in the area, and felt ground motions have resulted in the implementation of a traffic light protocol (TLP) (13).

Despite recent progress in characterizing the Duvernay-related earthquakes (14, 15), many scientifically critical questions have remained unresolved. For example, it is not clear why there was a large time delay (~3 years) between the first Duvernay play completion and the initiation of HF-related earthquakes (14), after which time such earthquakes became frequent occurrences. Moreover, it is not well understood why only a small subset of operations appear to be seismogenic (9). Within the Duvernay play specifically, only those completions located in the Kaybob region are seismogenic, whereas all HF completions in the Willesden Green and Edson regions have been seismically quiescent (fig. S1). Geological factors have been suggested to contribute to the spatial distribution of these earthquakes (16, 17). However, to date, there has been little understanding or documentation of the relative contributions of surface injection pressure, rate, and volume to the seismogenic process.

To address these questions, we first examined the timing and location of earthquakes in relation to HF completions in the seismogenic Kaybob region of the Duvernay play (Fig. 1). We compiled a database of all (~300) horizontal HF well completions in the Kaybob Duvernay up to February 2016 (because of the ~1-year period of confidentiality) from public, hard-copy regulator records. Individual wells are aggregated into ~180 well “pads” on the basis of the proximity of multiple wells oriented along similar trajectories (fig. S2). A cursory examination of the time-averaged evolution of injected volumes per pad (Fig. 1B) indicates increasing pad design complexity and HF completion volumes, typical of maturing development in shale plays (18). The Kaybob Duvernay has been injected with more than 8.5 × 106 m3 of fracturing fluid (as of February 2016) to stimulate well productivity. The observation of a relative increase in pad volumes before the first recorded earthquakes suggests that injected volume may be a controlling factor in the Kaybob Duvernay earthquake activity.

Fig. 1 HF-induced seismicity within the Duvernay play up to February 2016.

(A) Spatial distribution of induced earthquakes (red circles) associated with HF wells (dark gray “tadpoles”; associated wells are bolded and black) in the Duvernay play (purple shaded area) near Fox Creek (dark gray) and Crooked Lake (blue). The inset map shows the location of the Duvernay Formation in North America. (B) Timing and magnitude of induced earthquakes (red circles) alongside stage volumes (gray) and a 30-pad running average of pad volumes (black line). Colored areas indicate Alberta Energy Regulator (AER) traffic light protocol (TLP) cut-offs (13). Initial wells stimulated during March 2010 were included in the analysis but are not depicted.

Disposal is not likely to be a major factor in the induced seismicity in this area. The closest water-injection well to Crooked Lake that was actively injecting at depths similar to those of the Duvernay during the seismogenic period is ~35 km away. Disposal wells within a 50-km radius of Crooked Lake injected only ~1.2 × 106 m3 of fluid during the Duvernay’s seismogenic period (much less than the total volume involved in the HF operations).

To examine the role of operational factors more closely, we associated clusters of seismicity with seismogenic pads on the basis of a spatiotemporal association filter (SAF), which identifies the causally closest pad(s) in time and then space (supplementary materials). The validity of these associations can be demonstrated by comparison with prior case studies (14, 15) and direct communication with the regulator and the anonymous companies responsible. Based on the SAF, ~10% of the pads and ~15% of the wells in the Kaybob Duvernay are associated with seismicity. Of the completions associated with induced earthquakes, ~50% are single-well pads. Using this subset of associated pads, we then contrasted operational parameters at seismogenic pads, wells, and stages with those of their parent distributions, as derived from the entire Kaybob region of the Duvernay play (Fig. 2). Statistical distributions of operational parameters (fig. S3) were analyzed using the Kolmogorov-Smirnov (KS) test (19) to discern whether it is likely that the seismogenic distributions are sampled randomly from their parent distribution. In these tests, we used the standard significance level of 0.05. The computed P values (Fig. 2) allowed a rejection of the hypothesis that the seismogenic distributions are subsampled from the parent Kaybob Duvernay distribution. On the basis of finding nonrandomness, the Mann-Whitney (MW) U test (20) was applied next to determine whether the seismogenic subsets have significantly larger median values of key operational parameters than do their parent distributions. We performed both one- and two-tailed MW tests to first determine whether subset median values are different from the parent distribution (two-tailed) and then further assess whether the subset median values are larger than their parent distributions (one-tailed). This analysis (P < 0.05) demonstrated that seismogenic pads, wells, and stages are associated with larger injected volumes at a statistically significant level (Fig. 2). Analogously, an analysis of HF in the Horn River Basin found a relationship between induced earthquake productivity and volume (21). Examination of pad, well, and stage pressures suggests no significantly compelling association (fig. S3). Although injection rate has been suggested as a driving factor for disposal-induced seismicity in the central United States (22), we did not observe a meaningful association with injection rates. Potentially, this discrepancy could be due to differences in injection operations for disposal and HF; for example, the slowest Duvernay injection rates are nearly an order of magnitude faster than the critical rate threshold identified for seismogenic disposal (22). Thus, the effects of rate on HF-induced seismicity in the Kaybob region are either secondary, indiscernible, or negligible. The robustness of these findings was established using bootstrap (23) resampling sensitivity tests (supplementary materials and figs. S4, S5, and S6).

Fig. 2 Distributions of HF volumes in the Kaybob Duvernay.

Volumes are plotted as histograms on a per-stage (A), -well (B), and -pad (C) basis for seismogenic (red) and all (gray) pads. P values from comparisons of complete volume distributions with seismogenic subsets are inset [KS, Kolmogorov-Smirnov; MW, Mann Whitney; 1(2)T, one(two)-tailed]. The dotted line depicts the required detection threshold volume, Vc.

The results of the KS and MW tests show that volume is a controlling factor for HF-induced earthquakes in the Kaybob Duvernay. Potentially, this observation may be the result of greater injection volumes allowing for larger stimulated reservoir volumes and thus greater likelihood of intersecting a critically stressed fault (24). However, this does not necessarily imply that it controls the maximum possible magnitude. This finding echoes similar conclusions from other induced cases (25). The relationship between volume and induced earthquakes has important implications for the management of HF-related seismic hazard. For example, controls on seismicity related to volume would be affected by the time-dependent nature of increasing pad completion volumes within a maturing play (Fig. 1B and fig. S7). On the other hand, the use of lower pad completion volumes since the first TLP red light (Fig. 1B) could have resulted in an overall reduction in earthquake response from the Duvernay play. Whether this decrease in volume after the first TLP red light resulted in a net decrease in seismic hazard is more complicated, however, because hazard is dependent on multiple factors, some of which are independent of injection volume (26, 27). These points would require a physical model to validate the statistical earthquake-volume association.

To validate this relationship between injected volume and earthquakes (28), we first considered the Gutenberg-Richter frequency-magnitude distribution (GR-FMD) (29, 30)NM = 10a10bM(1)In this formulation, NM is the number of events greater than magnitude M, the a-value governs the rate at which earthquakes occur, and the b-value measures the proportion of relatively smaller events to larger ones. To consider the time-varying rate of induced earthquakes, modifications have been suggested on the basis of solutions to the diffusion equation that incorporate the nonstationary effects of injection volume V(t)—e.g., a = Σ + log10[V(t)] (25, 3133). In this equation, Σ is the seismogenic index, an injection-invariant parameter that represents the seismotectonic response to increasing fluid pressure within the study area (25). Incorporation of Σ in the GR-FMD modifies Eq. 1 to (25, 3133)

NM = V(t) · 10Σ10bM(2)

We investigated the implications of Eq. 2 for the Kaybob Duvernay HF earthquakes by first compiling the time-dependent histories of completion volumes and the number of earthquakes above the detection threshold (fig. S8). We observed a strong correlation between these variables, especially if only the SAF-associated pads were considered (Fig. 3). The best-fit parameters of Eq. 2 are b = 0.90 ± 0.03 and Σ = –1.8 ± 0.2 for the detection threshold of local magnitude (ML) 1.3 (supplementary materials and fig. S8). This regionally averaged Σ value near Fox Creek is very similar to Σ values computed for earthquakes in the Brazeau Cluster, which are related to wastewater disposal in the Cordel Field (–2.1 ± 0.2) (34), and to Σ values for numerous case studies worldwide (35). This is an important finding because it suggests that the regional seismic response to fluid injection in central Alberta is similar for both HF and disposal wells.

Fig. 3 Number of earthquakes above the detection threshold (ML 1.3) versus cumulative injection volume.

(A) Volumes from all HF operations within the Kaybob Duvernay play (red circles) are compared with the best fit of the data points (black line). (B) Analogous to (A), except only seismogenic HF pads are considered. Both panels are during a period where the detection threshold remains constant (July 2014 to February 2016). Well flow-back was not considered in computing volumes. In both panels, the goodness-of-fit of the data to the expected line is displayed with the R2 (coefficient of determination) values.

We repeated the fitting process for individual clusters associated with seismogenic wells and found Σ and b-values that vary from –2.5 to –0.5 and 0.7 to 1.7, respectively (fig. S9). The relatively high variability of these parameters likely reflects the difficulty in making robust parameter determinations from small data subsets. More complete earthquake catalogs would enable a more robust analysis of seismic parameters for the individual clusters, improving confidence and possibly reducing variability. Overall, we found seismic values for individual clusters that are roughly similar to Σ values previously determined for clusters on a local array (15). The determination of these values provides a powerful tool with which to forecast the expected number of future earthquakes at seismogenic wells and their likely magnitude distribution.

In light of these findings, we consider the hypothesis that the ~3-year-delayed response in earthquake productivity was simply the result of the minimum injection volume (Vc) required to raise the seismicity rate to a sufficient level for observation (i.e., so that it produces an earthquake larger than the regional detection limit). Conservative upper-bound estimates using a detection limit of ML 2.0 (the TLP requires operators to report all events of ML 2.0 and greater) and 95% probability of exceedance suggest a Vc of (8.0 ± 0.2) × 103 m3 (supplementary materials). This number is comparable to the average total injected volume at pads at the time of their first corresponding earthquakes, V1 = (1.0 ± 1.0) × 104 m3. Although Vc and V1 roughly agree, it is likely that V1 is systematically larger than Vc owing to insufficient resolution with which to discern aseismic stages or wells within a seismogenic pad. Corroborating these results, studies in the Horn River Basin found that seismic response appeared to “turn on” during months where HF injection volumes were greater than 2.0 × 104 m3 (21). Because fewer than 10 pads in the Kaybob Duvernay have a volume per pad less than Vc, we argue that it is unlikely that the spatial distribution of our catalog has been seriously biased by Vc. We justify this claim through comparison of our regional catalog (14, 36) with catalogs supplemented by local operator networks (15). In both cases, we observe the same spatial distributions of earthquakes and the same associations with seismogenic wells and regions.

Although pad completion volumes have been increasing with time, it is interesting that more than 20 pads (~50%) that were completed before the first induced earthquake have volumes greater than 1.5 × 104 m3. This finding indicates that although volume appears to be a controlling factor for Kaybob Duvernay induced seismicity, other factors are also playing a nontrivial role. Prior work has suggested that the spatial distribution of all induced seismicity within central Alberta has been influenced by geological factors, with earthquakes preferentially occurring along underlying reef margins (16, 37) or within regions of relatively higher formation overpressure (17).

Following this rationale, we accommodate spatial factors using a modification Embedded Image that explicitly introduces a spatial variable Embedded Image so that Eq. 2 becomesEmbedded Image(3)In our formulation, we define Embedded Image to have binary values of 1 or 0 to indicate regions that do or do not experience earthquakes, respectively. In this limiting case, Σ′ becomes equivalent to Σ [i.e., when Embedded Image]. Introducing this spatial term is a useful refinement, because numerous plays have been stimulated using similar per-well volumes (38, 39) with limited or undocumented seismic response. Even considering only the Duvernay, induced earthquakes have been restricted to one area, the Kaybob region (fig. S1). Furthermore, numerous high-volume Kaybob pads, wells, and stages are not associated with earthquakes (Fig. 2). In fact, the first seismogenic pads observed in the Duvernay (10) were also the first pads to complete in the most seismically susceptible region, ~10 km from Crooked Lake (fig. S7). The efficacy of introducing this spatial term for the Kaybob Duvernay is quantified by the significant enhancement of goodness-of-fit to the data when using the SAF (Fig. 3). In this sense, the SAF represents a rudimentary and empirical estimation of Embedded Image. Taken in the context of the original formulation, the spatial variability of Σ can be cast as Σ = –1.8 + log10(SAF). The incorporation of the SAF as an empirical estimate of Embedded Image in a Σ-modified GR-FMD results in a model that accounts for ~96% of the seismic response variability within the Kaybob data set (Fig. 3B). This means that complexities related to other potential factors account for only the remaining 4% of variability or constitute part of the proposed Embedded Image. This is an important finding because these additional factors are numerous; they include well flow-back, effects of staged stimulation pressures, well shut-in response, lagged seismicity relative to seismogenic pads, petroleum production, earthquake aftershock sequences, operator seismic mitigation strategies, varieties of fracturing fluids and proppants, SAF resolution limitations, interpad communication, poroelastic triggering effects, magnitude uncertainties, and spatiotemporal variability in b-values or Σ′. This observation is further confirmed by applying additional filtering to the data (figs. S10 and S11).

Although in reality, a fault is either reactivated or not, incomplete information about the subsurface often prevents definitive assessments a priori. Our spatial parameter is a useful concept in this regard. Statistically, we can extend Embedded Image to represent the seismogenic activation potential, defined as the likelihood of a well inducing a detectable earthquake at a given location. We interpret this parameter as the probabilistic intersection of all geological conditions required to cause an induced earthquake at a given location—i.e., the spatial variability of the geological susceptibility to induced seismicity. This interpretation is intuitive, owing to the derivation of Embedded Image from Σ, which constrains the seismotectonic state of the point of injection. For a demonstration of this interpretation, we consider the effects of distance to underlying fossil reef margins (16, 37) and formation overpressure (17) as regional proxies for faulting and stress, respectively. The fractions of seismogenically associated pads as a function of these proxies are plotted (Fig. 4), confirming that regions of development that are closer to the reef margin or more highly overpressured have been more likely to induce earthquakes (fig. S12). For western Canada, an entire basin-wide average of this activation probability in regions that are also coincident with viable HF plays appears to be less than 0.3% (9).

Fig. 4 Statistics of HF operations in relation to the Swan Hills Formation and Duvernay overpressure.

(A) Distance of HF operations to the Swan Hills Upper Bank edge for seismogenic (red) and all (gray) pads, shown as a histogram. (B) Similar to (A), except each bin has been normalized (red circles), and error bars have been added. (C) Pressure gradient for seismogenic (red) and all (gray) pads. (D) Fraction of seismogenic pads as a function of pressure gradient (red circles) with error bars from the standard deviation.

Solely on the basis of the Σ modifications, it is possible to improve hazard estimates for future injections in real time during stimulation. Knowledge of fault size, hydraulic connectivity to proximal stimulation stages, and injected volumes that may be directed to reactivated faults can serve as input to calibrated models to estimate induced earthquake rates and magnitude distributions. Conversely, microseismic monitoring of HF completions may be used to scrutinize the number and rate of induced events resulting from individual stage stimulations as an indicator of hydraulic connectivity to nearby faults. For example, HF pads that are oriented subparallel to the north-south–oriented fault planes (12, 14, 15) are likely to be in extended hydraulic communication with seismogenic faults as compared with northeast-southwest–oriented pads. This rationale may explain why the three largest-magnitude clusters occurred at north-south–oriented pads: Greater volume pumped into these faults would allow for more numerous induced events and thus an increased likelihood of a larger event. Extended to the regional scale, this study can be used to better inform earthquake rate models in induced seismic hazard forecasts. Coupled with an estimation of the seismogenic activation potential, our proposed framework (Eq. 3) would allow for the quantification of both induced earthquake rate and location models. Rate and location models are some of the most critical parameters in forecasting hazard related to induced earthquakes (26, 27). Although this study has focused on the Fox Creek area, the proposed framework can be applied to other jurisdictions to improve the management of induced seismic hazard.

We find that the most important operational parameter controlling induced earthquakes in the Duvernay play near Fox Creek is injected volume, which scales linearly with the total number of earthquakes in a Σ-modified Gutenberg-Richter formulation. Conversely, injection pressure and rate appear unrelated to induced seismicity response near Fox Creek. Furthermore, wells that exhibit seismicity appear to display a strong spatial bias related to geological factors, which has a pronounced effect on the resultant seismic response. Last, this study provides a framework with which to incorporate the seismogenic activation potential into seismic hazard analysis for induced earthquakes.

Supplementary Materials

Materials and Methods

Figs. S1 to S12

References (40, 41)

References and Notes

Acknowledgments: Waveform data are available through the Incorporated Research Institutions for Seismology, earthquake catalog data are available through the Alberta Geological Survey (14, 36), and tour reports of operator hydraulic fracturing parameters are available through the Alberta Energy Regulator. We thank K. Haug, G. Jean, S. Kolenosky, M. Speta, M. Araneta, and J. Chadi for their help in compiling the hydraulic fracturing database. We also thank H. Corlett and three anonymous reviewers for their comments. All authors declare no conflicts of interest.

Stay Connected to Science

Navigate This Article