Strategies to prevent damage to critical infrastructure due to induced seismicity

The hazard framework is focused on a potential induced seismicity source due to hydraulic fracturing operations in horizontal wellbores and to a lesser extent wastewater disposal. Since about 2008, tight oil and gas reservoirs in western Canada have increasingly been developed by drilling wells horizontally through the reservoir rock, at depths of a few km. The horizontal segment of the well typically extends ∼1.8–3 km, over which distances there is a controlled, high-pressure injection of fluid and proppant to fracture the target formation. Hydraulic fracturing in western Canada has, in some instances, been associated with significant induced seismicity. In the last few years, several sequences of seismicity including events at the M > 4 level have been induced by hydraulic fracturing in western Alberta and northeastern British Columbia (Atkinson et al. 2016; Bao and Eaton 2016; Schultz et al. 2017).

Figure 1 shows the association between oil and gas activity and seismicity in the Western Canada Sedimentary Basin (WCSB) on a regional basis, as determined by Ghofrani and Atkinson (2016) in a statistical study. They subdivided the region into cells of 10 km radius and evaluated for each cell whether there was a relationship between seismicity and hydraulic fracturing in horizontal wells (referred to as HF wells). The rationale for dividing the region into cells was to be able to evaluate susceptibility to seismicity in spatially discrete areas, rather than looking at individual wells. Such an approach is more useful for hazard forecasting in specific areas. Every cell that contains at least one HF well is classed as active, whereas the remainder of the cells are classed as null. A possible temporal correlation between M ≥ 3 seismicity and HF wells is identified if an event (M ≥ 3) occurred in a cell within a window beginning with the commencement of hydraulic fracturing and ending 3 months after the completion of treatment. If an active cell contains at least one event that meets the temporal criterion for correlation, then its grid point is a hit for potential induced seismicity. The total number of hit points (N _Hit) and the total number of active cells (N _Active) can be used to calculate the ratio R_{h 10} = N _Hit/N _Active, which is an estimate of the regionally averaged probability that operations in a 10 km radius cell may be associated with seismicity. We consider that the quantity R _hA is a rough estimate of the prior probability of inducing seismicity by commencing HF operations in a small area in the WCSB in a Bayesian sense; note that the area in this case is that associated with a circle of the given radius (i.e., an area of π(10 km)² = 314 km²). Based on this logic, and using a Monte Carlo approach to assess the statistical significance of the observed hit rates, Ghofrani and Atkinson (2016) concluded that the regionally averaged likelihood that HF well operations in a 10 km cell will trigger M ≥ 3 seismicity is ∼0.01–0.03. We can normalize this to obtain a per annum rate of M ≥ 3 events in a 10 km radius cell by considering the following conditions: (i) the average number of M ≥ 3 per hit cell = 2.13; (ii) a total of 71 hit cells; and (iii) the hit cells represent a total of 136.5 well years of operation, for the temporal association window of 3 months following each operation. Thus, the cellular activity rate can be normalized to a per annum rate by multiplying it by the factor (2.13)(71)/(136.5) = 1.1, which is near unity. We conclude that the per annum rate of induced M ≥ 3 events is ∼0.01–0.03 for a 10 km cell (i.e., area of 314 km²) in which HF well operations are being performed in the year. This net rate for induced events for HF wells (which includes the activation probability) is a key source parameter in hazard assessments. By expressing it per annum on a per area basis, it is straightforward to define an areal source zone for hazard analysis and specify its recurrence parameters.

The rate of occurrence of earthquakes of different sizes is expressed by the Gutenberg–Richter relation (Richter 1958), truncated at a maximum magnitude (Mx). Mx is the largest event that is physically possible in the zone. The Gutenberg–Richter relation states that the cumulative number of events of magnitude ≥M, considering events below Mx, is given by log N(≥M) = a − bM, where a is the rate parameter and b is the slope. An example is shown in Fig. 2 for the Fox Creek area of western Alberta (the middle purple box in Fig. 1), in which the average yearly recurrence statistics are shown for 3 years starting in mid-2013, along with the average rates per annum over the 3-year period. The earthquake data for this plot are derived from the Alberta Composite Seismicity Catalog (inducedseismicity.ca). A Gutenberg–Richter relation with the slope b = 1 is plotted for three different rates of activity. The lower two levels of activity correspond to the average cellular per annum rates (cell of 10 km radius) of 0.01 and 0.03 for M ≥ 3, which have been normalized to account for the area of the rectangular source defined for the Fox Creek region (e.g., by multiplying by a factor of (195 km × 165 km)/(π 10²)). Note that the defined box containing the Fox Creek region is more active than average for induced seismicity sources in western Canada. The Fox Creek activity rates in the last 3 years are greater than average rates for natural seismicity sources in the North American craton and Canadian shield, as determined by Atkinson and Martens (2007), by nearly a factor of 100. For the induced seismicity hazard analysis, we define a large rectangular box (such as the purple box in Fig. 1) that acts as an induced seismicity source zone and assign it the rate parameters shown in Fig. 2, assuming the middle rate (the cellular rate of 0.03 p.a. for a cell with a 10 km radius).

Among the event clusters in the WCSB, there have been several events of M > 4. The largest event attributed to hydraulic fracturing to date was M = 4.6 (Atkinson et al. 2016). The maximum magnitude of the events that could be triggered is not yet known. It has been postulated that for events triggered by fluid disposal, the maximum magnitude may be limited by the cumulative volume of fluid injected into the area (McGarr 2014). However, the hypothesis that magnitude is limited by the injected volume does not appear to hold for earthquakes triggered by hydraulic fracturing in western Canada. This is illustrated in Fig. 3, which shows that several of the larger events triggered by hydraulic fracturing in the WCSB have clearly exceeded the upper bounds established based on injected volume, even if we consider the maximum reasonable uncertainty in both the moment magnitude and the volume of fluid involved.

As a general principle, the maximum magnitude is limited by the size of the fault that is re-activated, whereas the expected frequency of large events (relative to small events) can be calculated based on Gutenberg–Richter scaling (Petersen et al. 2015; van der Elst et al. 2016). This is illustrated for seismicity in the WCSB in Fig. 4. The figure shows the number of events of M ≥ 3.5 in 5-year windows, ending at the end of 2016, as listed in the Alberta Composite Catalog (inducedseismicity.ca); the area covered is that from 48° to 53° N and 110° to 121° W, plus the area from 53° to 60° N and 110° to 124° W. In most 5-year windows since 1970, the total number (N) of M ≥ 3.5 events is ∼10, and thus for Gutenberg–Richter scaling (Richter 1958) with a slope (b-value) of 1, we would expect the maximum observed event to be near M = 4.5; in 1 out of 10 such 5-year windows, we would expect the maximum observed to be ∼1 unit higher, or M ∼5.5. Overall this is what we observe. Specifically, the maximum observed magnitude in a time window tends to be larger in windows with more events. The most active window is the most recent (ending at the end of 2016), for which the observed maximum magnitude of 4.6 is just slightly lower than would be expected based on the rate of events.

It is important to recognize that the maximum observed values in a time window are not the maximum possible values (Mx). We cannot determine Mx directly from observations because we cannot adequately sample events with low recurrence rates. Maximum possible events may have annual likelihoods of only ∼1/1000, for example, and thus be much greater than the 90th percentile expected values for any given time window (i.e., greater than the 90th percentile values plotted in Fig. 4). Considering that the maximum possible magnitude for induced events is poorly constrained and may not be significantly lower than that for natural events in the region, the hazard analysis will consider a range of possible values for maximum magnitude (Mx), ranging from 5.0 to 6.5. This range is slightly lower than that considered for natural earthquakes in the region (Halchuk et al. 2014). The analysis is not particularly sensitive to reasonable variations in the assumed values of maximum magnitude (Atkinson et al. 2015a).

A key issue in the assessment of hazard from induced seismicity is the specification of ground motion amplitudes associated with induced events, particularly at close distances (within tens of km of the source). These are commonly estimated using ground motion prediction equations (GMPEs) that express peak ground amplitudes and response spectra (5% damped pseudo-spectral acceleration (PSA)) as a function of moment magnitude, distance, and other variables. There is a paucity of well-calibrated and robust GMPEs derived strictly from induced events, although this situation is rapidly changing as new data sets become available. Therefore, induced seismicity hazard applications have tended, of necessity, to draw on models of ground motion amplitudes developed from natural events. This was the rationale for the Atkinson (2015) (denoted A15) GMPE, which examined the magnitude and distance scaling for small-to-moderate events (M 3–6) at hypocentral distances (R _hypo) < 40 km. The A15 GMPE was derived by empirical analysis of natural events in California. However, it explicitly addressed the point-source scaling attributes and near-distance saturation behavior expected for induced events, thus providing a middle ground between GMPEs developed entirely for natural events, and GMPEs developed entirely for induced events. A follow-up study (Atkinson and Assatourians 2017) demonstrated that at least three published GMPEs have attributes that make them good proxies for induced earthquakes, as follows: the Atkinson (2015; A15), Yenier and Atkinson (2015; YA15), and Abrahamson et al. (2014; ASK14) GMPEs. (Note that an unspecified depth to top of rupture should be prescribed with ASK14 if it is to be used for this purpose; see Atkinson and Assatourians 2017 for details.)

Figure 5 shows that these GMPEs are appropriate in functional form and overall amplitude level, as they match ground motion observations from induced events in Oklahoma and Alberta. We may, therefore, use these GMPEs to define a representative suite that models amplitudes and their scaling for use in hazard analysis. A two-GMPE suite is used to characterize ground motion amplitudes, as shown in Fig. 5; the upper GMPE is the A15 model, whereas an alternative model has been defined by applying a distance-dependent scaling factor that reduces the motions, especially at very short distances.

Of note in Fig. 5 is the high variability (scatter) of amplitudes about the median GMPEs. This is an important contributor to seismic hazard, as it means that ground motions are often significantly stronger than the median level, increasing the probability that a moderate event at close distance can cause strong ground motion. Another point of note is the level of observed ground motions for the magnitude range shown, M 4.0–4.5, relative to the damage threshold, which is taken as the motion that would typically result in a Modified Mercalli Intensity (MMI) of VI or greater. The ground motions corresponding to MMI = VI are estimated using the empirical relationships of Worden et al. (2012) between MMI and ground motion amplitudes. For high frequency ground motion parameters such as peak ground acceleration (PGA), the amplitudes for events of M ∼4 often exceed the damage threshold, whereas for lower frequency measures such exceedances will be less likely. This implies that the induced seismicity hazard is primarily a high frequency hazard, and thus, structures that are not sensitive to strong high frequency vibrations (i.e., long-period structures) would be less affected. The damage threshold is discussed in more detail later in the paper.

The frequency characteristics of typical induced events can be better appreciated by looking at plots of response spectra, which show the amplitude level versus vibration frequency for selected records. Figure 6 shows a sample of such records for induced events in Oklahoma of M 4.5–5.0 recorded in the hypocentral distance range from 5 to 10 km (corrected for site conditions, from Class C to B/C boundary). At high frequencies, the sample records tend to lie within the bounds that would be expected based on the Atkinson (2015) GMPE; they are slightly weaker than predicted at low frequencies. The ground motions that correspond to an MMI of VI, according to the relations of Worden et al. (2012), are also shown for this magnitude and distance range. This illustrates that these motions are near the damage threshold at high frequencies and below it for low frequencies.

There is no single definitive measure of what makes ground-motion damaging. In general, the damage potential of ground motion depends on a combination of factors including amplitude, frequency content, and duration. Small-to-moderate events tend to be of short duration and are rich in high frequency content (at close distances), whereas larger events are longer in duration and have significant energy at longer periods. Moreover, the damage potential of motions will also depend greatly on the structural characteristics of the facilities in question, such as their overall robustness and frequency range of response. We can use MMI as a simple generic measure of damage potential that has a sound and well-calibrated empirical basis. In general, MMI = VI is the threshold for minor damage such as cracking to walls and chimneys. By comparison, MMI = III is the felt threshold, whereas MMI = VIII is the threshold for significant damage to well-engineered structures.

We can estimate ground motions that would typically result in MMI VI using empirical relations between ground motion amplitude and MMI. There are several such relations in the literature. For this exercise, the relations of Worden et al. (2012) were adopted, because they are based on a rich empirical database (from California) that explicitly includes the effects of amplitude, magnitude, and distance on MMI. It is important to realize that an inherent limitation of MMI is that it is a single number, whereas the underlying ground motion contains a range of amplitudes over an entire spectrum of frequencies. Thus, the relationship between MMI and ground motion varies with frequency; consequently, a single ground-motion record could be associated with a range of MMI values, with the most appropriate being either some average measure or possibly being dependent on the application (i.e., depending on whether high frequency or low frequency is of primary interest). Numerous authors have noted that, in general, PGA tends to be well correlated with felt effects (and lower MMI levels), whereas peak ground velocity (PGV) tends to be well correlated with damage (and higher MMI levels) (e.g., Worden et al. 2012). Worden et al. (2012) concluded that an optimal approach for estimating MMI from ground motion is to use the average value as predicted from PGA and PGV. Atkinson et al. (2014) followed this approach in their development of GMPEs for intensity in North America. This explains why MMI tends to be higher in the east than the west for the same magnitude of earthquake, even at close distances. The higher frequency content of eastern earthquakes and slower attenuation of high frequency motions within the competent eastern crust result in greater PGA compared with western events (whereas the PGV is similar) (e.g., Atkinson and Boore 2011).

The relations of Worden et al. (2012) predict that the median PGA associated with MMI = VI is 170 cm/s² for moderate events at close distances (M ∼4.5 at ∼10 km). This intensity level is also associated with median response spectral amplitudes of 69 cm/s² for PSA at 1 Hz and 363 cm/s² for PSA at 3.33 Hz. These values are plotted as the damage threshold in Figs. 3–5. Note that these are median values, and there is significant variability in the amplitudes associated with MMI = VI (standard deviation of about a factor of two).

For PGV, the relations of Worden et al. (2012) predict that MMI = VI is associated with a median PGV = 9 cm/s. This is in good agreement with the range of values of PGV that are associated with damage potential based on blasting criteria and used to limit vibrations from quarry blasting and similar operations. For example, Westaway and Younger (2014) suggested a value of PGV = 5 cm/s for plaster cracking, with PGV = 6 cm/s being the damage threshold, whereas Siskind et al. (1980) suggested that damage to mortar or block foundations is associated with PGV = 12.7 cm/s. Considering this range of estimates, a value of PGV = 10 cm/s was adopted as the damage threshold for induced earthquakes.

Finally, a guide to the overall damage potential of a recorded ground motion can be obtained using the relations of Worden et al. (2012) to estimate the MMI associated with the PGA and PGV values (at the median level), then taking the average of these two estimates. If this average MMI value exceeds VI, we consider that the motion has damage potential. Whether damage would actually result depends on the characteristics of the structures at risk.

A conceptually simple approach to prevent damaging ground motions is to preclude events that could cause damage by specifying setbacks of operations from critical facilities—an exclusion zone within which operations will not occur. This ensures that earthquakes are not generated at very close distances. Such an approach requires the definition of a scenario event upon which to base the setback distance, which is a subjective exercise. If we consider hydraulic fracturing operations, Atkinson et al. (2016) estimated that ∼0.3% of HF wells will be associated with earthquakes of M ≥ 3. If we assume Gutenberg–Richter scaling with a typical slope of b = 1, it follows that the 1/10 000 event associated with completion of a HF well has a magnitude of M ∼4.5. We can calculate the PGA and PGV associated with an M = 4.5 event as a function of distance using GMPEs. Based on the discussion in the preceding sections, we take the geometric mean of the values predicted by the three GMPEs shown to be good proxies for induced events (Abrahamson et al. 2014; Atkinson 2015; Yenier and Atkinson 2015). The motions are evaluated at the median plus one standard deviation level; the use of the median plus sigma is a common practice in deterministic hazard assessments, to account for variability (e.g., Reiter 1990). Finally, the values of PGA and PGV are converted to MMI using the median relations of Worden et al. (2012); the averages of the PGA- and PGV-based MMI estimates are used to assess the expected MMI. This intensity is plotted as a function of distance from the hypocenter in Fig. 7.

From Fig. 7, we conclude that events of M ∼4.5 are potentially damaging (MMI = VI) if they occur within a hypocentral distance of 6 km. If we assume that any triggered events would be initiated very close to the well, at a focal depth of ∼3 km we would require an exclusion distance for HF wells of 5 km laterally (i.e., in map view) for critical infrastructure to preclude potentially damaging motions. The appeal of a deterministic approach is its conceptual simplicity. However, it is limited by considering only a single magnitude level at a fixed sigma level, and that we have considered likelihood in only a very general way. The likelihood of the considered scenario is ∼1/10 000 (assuming an activation probability of 0.3%), but the total likelihood would, in reality, include other scenarios that we have not considered. It is thus difficult to assess whether or not the likelihood of potentially damaging motions is acceptable.

A more rigorous way to assess the ground motion hazard from induced events is to use a probabilistic approach, in which we can more readily consider the range of scenario events and their likelihoods and address their uncertainties. This is the aim of a classic PSHA. As described previously, the key inputs for a PSHA are the seismic source model, the recurrence statistics within seismic sources, and the GMPEs—along with the uncertainties and variabilities in these parameters. The use of the PSHA approach for induced seismicity applications was described by Atkinson et al. (2015a, 2015b). In this application, the model elements developed in the preceding sections are used to assess the induced seismicity hazard from hydraulic fracture operations in the area surrounding a site of interest. This facilitates the development of strategies that would limit the likelihood of potentially damaging ground motion to an acceptable level. For critical facilities, it is assumed that the maximum acceptable probability of potentially damaging motions (defined here as MMI = VI) is 1/10 000 p.a. Note that this considers only the induced seismicity hazard, whereas the total hazard is the sum of the induced and natural hazards. However, because the calculations carry a high level of uncertainty, this is a factor that is neglected for this exercise.

The PSHA is conducted using the Cornell–McGuire method for a source area that surrounds a site of interest, assuming that HF operations will be occurring at average regional rates for the WCSB (i.e., such as those over the last few years). The cellular rates of Ghofrani and Atkinson (2016), normalized for the area under consideration, are used to define the rate parameter (a value) of the Gutenberg–Richter relation for the source zone, as described previously. Specifically, a square source zone is defined that is 50 km × 50 km (2500 km²), with the site in the center. The rate of M ≥ 3 events for this square area is ∼0.075–0.25 p.a. for the assumed typical density of operations. The upper range of this rate is similar to that observed in the Fox Creek area, as noted previously. This rate (i.e., 0.25 M ≥ 3 p.a. in a box 50 km × 50 km) is adopted as the default for the calculations, assuming a b-value of 1.

The PSHA is conducted using the EQHaz software of Assatourians and Atkinson (2013). EQHaz conducts a Cornell–McGuire analysis based on a Monte Carlo approach, in which an earthquake catalog is simulated for the specified geometry and rate parameters of the source zones. Ground motions at the site are calculated for each event in the simulated catalog, drawing from a suite of GMPEs and their associated variability. The frequency statistics of the simulated ground-motion catalog define the hazard curves and can be used to explore their variability and uncertainty.

A minimum magnitude of 4.0 is assumed for the hazard calculations (see Atkinson et al. 2015b), meaning that only events above this magnitude will contribute to the calculated hazard. The maximum magnitude is assigned a distribution within the range of 5.0 to 6.5, using weights of 0.2, 0.3, and 0.5 for Mx = 5.0, 6.0, and 6.5, respectively. These values are slightly lower than those used in the national building code calculations for natural earthquakes in this region (Halchuk et al. 2014). The events are assumed to be distributed at depths from 3 to 6 km, using depth distributions of events of 2.0, 4.0, and 6.0 km with weights of 0.2, 0.5, and 0.3, respectively, and allowing a random 1 km perturbation in all depth values. The reader is referred to Assatourians and Atkinson (2013) for details of the hazard software and calculations; more details on hazard calculations for induced seismicity sources are given in Atkinson et al. (2015a, 2015b).

For the GMPE suite, the three candidate GMPEs discussed earlier (A15; YA15; ASK14) are used to define an uncertainty band around the A15 GMPE, as shown in Fig. 5; this is based on the representative suite approach (Atkinson and Adams 2013). The central GMPE is the A15 median GMPE, and a lower alternative is defined to represent the values suggested by alternative GMPEs. The lower alternative is defined by subtracting the value Δ, in log units, from the A15 GMPE, whereThis is a function that reduces the median amplitudes relative to A15 by 0.2 log units at very close distances, transitioning to a reduction of 0.1 log units at hypocentral distances (R _hypo) of 10 km and greater. The A15 GMPE was convenient to use as the central GMPE because it is defined in terms of hypocentral distance. If larger events (>6.5) are considered, adjustments could be needed to accommodate potential saturation problems with this GMPE at large magnitudes. It is assumed that the standard deviation of ground motions about the GMPEs is 0.25 log units at low frequencies (and PGV), increasing to 0.3 log units at high frequencies (and PGA).

EQHaz is used to simulate catalogs of events and associated ground motions that would be expected if there are no exclusions on HF operations around the site. We can then examine this catalog to find what distance range of events needs to be precluded or controlled to ensure that the ground motions at the 1/10 000 p.a. level do not exceed MMI = VI. For this purpose, a synthetic catalog that is 1 000 000 years in length is generated. The ground motion having a probability of exceedance of 1/10 000 p.a. is the 100th largest amplitude in this catalog (i.e., 100 exceedances in 1 000 000 years are equivalent to an average exceedance rate of 1 in 10 000 years). Note that the catalog length is not intended to represent an actual duration of 1 000 000 years; rather, it can be thought of as 1 000 000 random realizations of a 1-year catalog of events. The large majority of these catalogs (>99%) will not result in a potentially damaging ground motion at the site, but rarely (<1% chance per annum), a large event will occur close to the site resulting in very strong ground motions.

Figure 8 shows the expected intensity (from the average values from PGV and PGA) at the site for the 1 000 000 year catalogue. It is observed that exceedances of MMI = VI come from moderate events at close distances, plus the rare possibility of larger events (M > 6) that would typically be at larger distances. The larger events are very rare (from the Gutenberg–Richter relation), so the odds of getting one very close to the site are very low. The total number of exceedances of MMI = VI (for the assumed cellular rate parameter of 0.03), with no exclusion zones, is >200, or more than double the target allowable number of 100 for 1/10 000 p.a. Figure 8 also shows the effect of an exclusion zone on the number of exceedances of MMI = VI. As we remove events within a certain distance range, by prescribing an exclusion zone about the site, the number of exceedances of MMI = VI will decrease, but the effect only becomes significant for exclusion distances of >5 km. The data in Fig. 8 suggest that if we are relying entirely on an exclusion zone to reduce the hazard, the exclusion distance would need to be ∼10 km to reduce the likelihood of potentially damaging motions to <1/10 000 p.a. This is greater than the deterministically derived exclusion distance of 5 km, because the probabilistic analysis indicates that there are significant hazard contributions from events of M > 4.5 at distances beyond 5 km.

An alternative strategy to reduce the hazard is to combine an exclusion zone with a plan to ensure that the rate of induced events in the region, to a distance of ∼25 km, remains low—thus also reducing the likelihood of strong ground motions coming from more distant events. Figure 8 shows the effect of combining an exclusion zone with a stipulation that the cellular rate in the region, to 25 km, shall not exceed 0.01 p.a.; this has the effect of dividing the number of exceedances by a factor of 3 (from a rate of 0.03 to 0.01). From Fig. 8, we infer that the combination of an exclusion zone of 5 km laterally (6 km in hypocentral distance), plus a provision that the maximum cellular rate (for 10 km radius cell) be kept to <0.02 p.a. in the area extending out to 25 km, will achieve the objective of limiting the likelihood of potentially damaging motions to <1/10 000 p.a. The equivalent rate for the 25 km radius is a factor of 6.25 times greater (scaled for the larger radius), giving a maximum allowable rate parameter of 0.13 for M ≥ 3 or ∼1 event per annum for M ≥ 2 in the 25 km radius.

Finally, Fig. 8 shows that an overall control on the seismicity rate within 25 km could achieve the target probability without any exclusion zone. It may, therefore, be tempting to conclude that no exclusion zones for operations are required. However, such a conclusion would imply that there is no risk distinction between operations directly beneath a critical facility and those 20 km away. By contrast, common sense and an understanding of the scaling of ground motion with distance dictates that measures to avoid the potential for significant earthquakes should be focused preferentially on the near-distance range—especially considering the uncertainties in assessing the likelihood of such events. It is primarily the uncertainty in assessing the likelihood of triggering earthquakes, coupled with the high consequences of an error, which drive the cautionary logic behind an exclusion zone (i.e., the need to couple mitigation with avoidance).