Increasing seismicity in the U . S . midcontinent : Implications for earthquake hazard

Earthquake activity in parts of the central United States has increased dramatically in recent years. The space-time distribution of the increased seismicity, as well as numerous published case studies, indicates that the increase is of anthropogenic origin, principally driven by injection of wastewater coproduced with oil and gas from tight formations. Enhanced oil recovery and long-term production also contribute to seismicity at a few locations. Preliminary hazard models indicate that areas experiencing the highest rate of earthquakes in 2014 have a short-term (one-year) hazard comparable to or higher than the hazard in the source region of tectonic earthquakes in the New Madrid and Charleston seismic zones.


Introduction
In July 2014, the U. S. Geological Survey (USGS) released an update of the 2008 National Seismic Hazard Map (NSHM) (Petersen et al., 2014) for the coterminous United States (Figure 1).The NSHM provides guidance for the seismic provisions of building codes and portrays ground motions with a 2% chance of being exceeded in an exposure time of 50 years.In the tectonically active western United States, the hazard model is derived from a combination of geologic studies of active faults, geodetic measurement of crustal deformation, and seismological observations of earthquakes and the shaking they produce.Over most of the central and eastern United States (CEUS), the hazard model is derived by projecting the spatially smoothed historical earthquake rate, with aftershocks removed, forward in time under the assumption that seismicity is described by a timeindependent Poisson process.Parts of the midcontinent, however, have experienced increased seismicity levels since 2009 -locally by two orders of magnitude -that are incompatible with the underlying assumption of a constant rate of earthquake occurrence (Ellsworth, 2013).
The 2014 NSHM acknowledged this problem but disregarded the increased activity.This was deemed appropriate for the specific purpose of underpinning seismic design standards for new construction.The 2014 NSHM excluded selected earthquakes in 14 William L. Ellsworth 1 , Andrea L. Llenos 1 , Arthur F. McGarr 1 , Andrew J. Michael 1 , Justin L. Rubinstein 1 , Charles S. Mueller 1 , Mark D. Petersen 1 , and Eric Calais 2 areas, of which eight were identified as sources of induced seismicity in the 2008 edition of the NSHM (Petersen et al., 2014).Both the developers of the NSHM and its users acknowledge that the rise of seismicity in the midcontinent must be understood if we are to capture fully the earthquake hazard in both space and time.
The space-time distribution of the post-2009 increased seismicity, as well as numerous published case studies, strongly implies that much of the increase is of anthropogenic origin and primarily is associated with wastewater disposal in UIC Class II wells (Frohlich, 2012;Horton, 2012;Keranen et al., 2013;Keranen et al., 2014;Kim, 2013;Andrews and Holland, 2015).If so, the assumptions and procedures used to forecast natural earthquake rates from past rates might not be appropriate.Here we examine the changing earthquake activity rates and discuss key issues that must be resolved to quantify the hazard posed by increased seismicity in the midcontinent.2011 M 5.3 Trinidad, Colorado, and2011 M 5.6 Prague, Oklahoma, earthquakes (Figure 2).To study how seismicity rates have changed, it is necessary to establish the magnitude completeness of the catalog.This can be done by examining the cumulative magnitude-frequency distribution (MFD, or Gutenberg-Richter relation).The MFD for earthquakes invariably follows an exponential distribution that can be written as log 10 (N(M)) = a -bM, where N(M) is the number of earthquakes in a sample with magnitudes ≥ M.Over the magnitude range where a catalog is complete, the cumulative MFD follows a straight line with slope -b on a linear-log plot.

Recent changes in seismicity
Magnitude-frequency plots for the eastern and western areas are divided into time periods corresponding to the early instrumental era  and four subdivisions of the modern instrumental era (1967-2014) (Figure 3).The observed MFD distributions in the eastern area for all time periods (Figure 3b) have consistent aand b-values, indicating that a time-independent Poisson process should have predictive power for future earthquake occurrence over a 50-year period appropriate for building codes.The catalog also appears to be complete for M ≥ 3.5 over the entire time period as well as complete for M ≥ 3.0 after 1967.
The situation is quite different in the western area (Figure 3a).Although the magnitude of completeness agrees with the eastern area for the same time periods, the earthquake rate (a-value) increases systematically with time after 2000.We can be confident that the increased seismicity is not a reporting artifact because the same procedures were used to analyze the earthquakes throughout the CEUS, and hence any procedural change would show up in both regions.The constant rate of earthquake activity (M ≥ 3.5) from 1930 through 2000 increased slightly between 2001 and 2008, and it underwent a major increase in 2009.The excess earthquakes in 2001-2008 appear to be located exclusively in the Raton Basin on the Colorado-New Mexico border (Rubinstein et al., 2014).From 2009 to the present, most of the excess activity has occurred in Oklahoma, with contributions from Arkansas, Colorado, Kansas, New Mexico, and Texas.
The evolution of the earthquake rate changes can be visualized by plotting the running count of earthquakes (Figure 4).The cumulative curve in the eastern area behaves as would be expected for a constant random process that combines a time-independent Poisson process augmented by an aftershock process.The west region has intervals that might be satisfied by such a model, but there are two clear violations of the constant-rate requirement.The first occurs between 1962 and 1970, and it can be attributed entirely to the Denver earthquakes induced by fluid injection in a deep well on the Rocky Mountain Arsenal (Healy et al., 1968).After removing these events from the running count, it agrees with a constant-rate random process until 2009, when the rate and newspapers.The first documented reports of earthquakes with epicenters in the midcontinent west of the Mississippi River, for example, date from the 1840s.Felt reports of earthquakes provided the primary data used to estimate location and magnitude until the 1930s, when instrumental recordings became the primary data source.
Using a newly compiled earthquake catalog for the CEUS (Petersen et al., 2014), we examine the earthquake rate in two equal-sized areas: an eastern zone that includes the M 7+ 1811-1812 New Madrid earthquakes, the M 6.9 1886 Charleston, South Carolina, earthquake, and the M 5.6 2011 Virginia earthquake; and a western zone that includes the principal areas where seismicity has increased in the past decade, including the  increases manyfold.We also can model the seismicity rate in each region using an inhomogeneous (time-varying) random process.A simple model for the rate parameter derived by fitting a smooth function to the cumulative curves for the eastern and western areas (without the Denver earthquakes) shows that the eastern rate is nearly constant, although the western rate increases slightly in 2001 and then accelerates after 2009 (Figure 4).
We can identify quantitatively where the excess seismicity is occurring in both space and time by comparing the observed seismicity with the expected number of earthquakes from a prior model.By comparing the expected and observed number of earthquakes, the probability of observing at least as many events as in the actual catalog can be determined.Areas with improbable numbers of earthquakes correspond to low p values ( p ≤ 0.05) in Figure 5.Note that the anomalous seismicity extends beyond the zones defined in the 2014 NSHM in northern Oklahoma and southern Kansas (Figure 1), which used earthquake data only through 2012.

Natural or induced?
This is the most challenging and perhaps the most contentious question being asked about the recent change in seismicity.Although some people have argued that more research is needed or that the scientific case for a connection between industrial activities and earthquakes has yet to be made, an extensive and solid body of scientific investigations, theory, and even experiments documents how earthquakes can be induced (Nicholson and Wesson, 1990).As early as 1968, the potential for inducing earthquakes by wastewater injection was acknowledged as a risk factor for underground waste management (Galley, 1968).That also was the year the link between injection at the Rocky Mountain Arsenal and the Denver earthquakes was established firmly (Healy et al., 1968).The Denver earthquakes occurred as fasr as 10 km from the injection point and persisted for more than a decade after injection ceased.Until 2011, the Denver earthquakes were the largest (M W 4.9) known to have been induced by injection.
It was hypothesized at the time that the increase in the porefluid pressure from injection was the mechanism that activated an ancient fault.The pressure required to overcome the frictional resistance to slip is given by the effective-stress relation τ crit = μ(σ n -P) + τ o .At failure, the critical shear stress τ crit equals the product of the coefficient of friction μ and the effective normal stress given by the difference between the applied normal stress σ n and the pore pressure P plus the cohesive strength of the fault τ o .The fault will remain locked as long as the applied shear stress is less than the strength of the contact.Increasing the shear stress, reducing the normal stress, and/or elevating the pore pressure can bring the fault to failure, triggering the nucleation of the earthquake.
This theory was tested in the field by the USGS in the Rangely, Colorado, oil field between 1969 and 1973 (Raleigh et al., 1976).The shear stress and initial pore pressure were measured in situ, and the coefficient of friction of the rocks was measured in the laboratory, from which the critical pore pressure needed to destabilize the faults was predicted.The experiment not only established the validity of the effective-stress model but also demonstrated that the pore-pressure perturbation extended kilometers from the injection point and responded rapidly to pressure changes at the wellhead.
It is rare to have the extensive knowledge of the stress, pore pressure, and hydrologic conditions that was available for the Rangely experiment.Consequently, most studies must rely on spatial and temporal correlations when assessing the possible involvement of industrial activity in the seismicity.Many of the areas with significantly elevated seismicity identified in  geologic data.In many cases, the evidence points to wastewater disposal by injection as the cause (Frohlich, 2012;Horton, 2012;Keranen et al., 2013;Kim, 2013;Keranen et al., 2014;Rubinstein et al., 2014;Andrews and Holland, 2015;Hornbach et al., 2015).
Although hydraulic-fracturing treatments are known to induce earthquakes occasionally (British Columbia Oil and Gas Commission, 2014), there was no detectable seismic response (M ≥ 2.5) in Pennsylvania during the development of the Marcellus shale (Ellsworth, 2013).It is unknown to what degree earthquakes induced during fracking operations have contributed to the rise in seismicity in the western area because the times and locations of fracking jobs are not presently available in most states.Earthquakes also accompany enhanced oilrecovery operations in some fields, such as at Rangely, Colorado, and in the Cogdell field in Texas, but not others, even in the same structural province (Gan and Frohlich, 2013).Productioninduced earthquakes are common in geothermal fields and have been observed in gas fields in southern Texas (Nicholson and Wesson, 1990) as well as in the Groningen field in the Netherlands (Muntendam-Bos et al., 2015;van Thienen-Visser and Breunese, 2015), among others.
Several alternatives to the fluid-injection hypothesis, including drought, recharge of reservoirs as drought eases, and natural variations in earthquake activity, have been put forward as explanations for the increased seismicity in Oklahoma and Texas, although none that we are aware of appears in the published literature.Because drought had been particularly severe in north Texas, Hornbach et al. (2015), in a study of the Azle-Reno area earthquakes of 2013-2014, consider the effects of both a depressed water table and changes in lake level in a nearby reservoir on the source volume of the earthquakes.They show that the maximum stress changes from either mechanism were one to three orders of magnitude smaller than stress changes associated with stress triggering, whereas the pressure increase at the fault from nearby injection wells was sufficient to induce earthquakes.
The magnitude of the earthquake rate change in the western area after 2008 (Figure 4) precludes a random fluctuation about the historic mean as a viable explanation.Could we be witnessing a tectonic event?To test this possibility, we examined the crustaldeformation rate in Oklahoma by using continuous GPS data from 2002 to the end of 2014, following the procedures of Calais et al. (2006).Across Oklahoma, GPS data contain no evidence for contemporary crustal deformation (Figure 6).Only one site, OKDT, has a displacement rate that is nominally different from zero at the 95% confidence level; however, it is not confirmed by neighboring sites OKAO, PRCO, and OKAD.The apparent differential motion between OKDT and OKTU, located on opposite sides of the seismicity, corresponds to an extensional strain of 1 × 10 -8 /year between the sites.This, however, is in the direction of maximum horizontal compression.If real, this strain would reduce the stress driving the earthquakes.We found no evidence for deformation anywhere in the state that would increase the stress driving the contemporary seismicity and concluded that there is no evidence for a natural tectonic origin of the increased seismicity.
If disposal of wastewater by injection is the principal cause of the excess seismicity, as now appears almost certain (Keranen et al., 2014;Andrews and Holland, 2015), it nonetheless needs to be stated clearly that disposal of wastewater by injection in UIC Class II wells more often than not results in no detectable seismic response.Consequently, the existence of a well has low predictive power for seismicity by itself.Frohlich et al. (2015) study the relation between earthquakes and disposal wells in the Bakken Shale of North Dakota and Montana using data from the USArray Transportable array from 2008 through 2011.Only three earthquakes were found that could be associated with disposal wells.Using the same methods and data, Frohlich (2012) finds dozens of earthquakes in the Barnett Shale of north Texas that clustered near several high-volume injection wells, but none associated with other injection wells.This suggests that for increased fluid pressure to induce earthquakes, three conditions must be met: (1) a preexisting fault must be present; (2) the fault must be oriented suitably in the tectonic stress field to slip; and (3) the pore-pressure perturbation must be sufficient to overcome the frictional strength of the fault.

Seismic hazard
Quantitative estimates of the likelihood and severity of future earthquake shaking in the NSHM are based on the method of probabilistic seismic hazard analysis (PSHA) Cornell, 1968).In PSHA, the likelihood of shaking at a site (not the occurrence of an earthquake) derives from the distribution of potential future earthquakes around the site (some could be hundreds of kilometers distant), the rate of occurrence for each earthquake source (as a function of earthquake magnitude), and the expected shaking at the site for each possible earthquake.
The increased earthquake activity in the midcontinent presents some particular challenges for probabilistic seismic hazard analysis.Conventionally, the goal of PSHA is to provide quantitative guidance for seismic provisions of building codes corresponding to the serviceable life of buildings, bridges, and other structures.How can this be done when the earthquake process is nonstationary in both space and time?We might be able to make viable forecasts of tomorrow's earthquakes based on the past week, month, or year.However, it is difficult to compute hazard at the exceedance levels of the NSHM for the next 50 years.Are there alternatives for quantifying the hazard over shorter time periods, and for whom would they be of value?
These questions were explored at the National Seismic Hazard Workshop on Induced Seismicity, cohosted by the Oklahoma Geological Survey and USGS in Midwest City, Oklahoma, on 18 November 2014.The workshop was attended by more than 150 participants representing petroleum producers; geophysical service providers; university, private-industry, and government researchers; state geological surveys; state and federal regulators; the reinsurance industry; and users of hazard models, including local government, state departments of transportation, and state architects.There was broad agreement that a one-year forecast would have value, and this is the present focus of the USGS effort for the NSHM (Petersen et al., 2015).The workshop also highlighted some key challenges for making one-year hazard forecasts, including (1) rapidly evolving temporal and spatial patterns of seismicity, (2) uncertainty in the MFD for induced earthquakes, and (3) potential differences in ground-motion prediction equations between natural and induced earthquakes.
Epistemic uncertainty in the components of the hazard model can be incorporated through the use of a logic tree that considers alternative models.Petersen et al. (2015) explore in detail the sensitivity of the hazard to alternative choices for the parameters in the logic tree.A representative example of one path through the logic tree shows the hazard to be elevated broadly across Oklahoma and surrounding states (Figure 7).In north-central Oklahoma, hazard is comparable to or higher than the hazard in the source region of tectonic earthquakes in the New Madrid and Charleston seismic zones.

Discussion
As of this writing in early 2015, the seismicity rate in north-central Oklahoma and south-central Kansas continues at unprecedented levels.In just the month of January 2015, residents of Oklahoma reported feeling more than 130 earthquakes, including 59 with magnitudes between 3.0 and 4.3.In southern Kansas, 21 earthquakes strong enough to be felt occurred, including nine with magnitudes between 3.0 and 3.9.What risk do future earthquakes pose to those living in the affected areas?
The current seismicity level in the midcontinent also needs to be kept in perspective.Since the upsurge in earthquake activity in 2009, a total of 10 events with M ≥ 4.5 has occurred; seven of these in 2011 and one each in the following three years.The largest (the M 5.6 2011 Prague, Oklahoma, earthquake) resulted in insured losses in excess of $10 million and sent some injured people to the hospital for treatment.Fortunately, no one was killed in that earthquake or any of the other events.Damage to unreinforced masonry buildings occurred in several of the other events, including the M 5.3 2011 Raton Basin, Colorado, earthquake and the M 4.8 Timpson, Texas, earthquake.Although reports of damage have been received for others, the main effect of most of the hundreds of felt earthquakes has been the nuisance of unexpected shaking, fraying nerves, and anxiety over the unknown potential for stronger shaking in the future.
Similar concerns have been raised in many other countries.For example, residents of the Netherlands have been experiencing shaking from earthquakes induced by long-term production of gas from the Groningen field (Van Eck et al., 2006).The Dutch Petroleum Company and the government are working together to manage the hazard by adjusting production to reduce the risk by strengthening structures and to compensate those who have suffered losses.
Going forward, the most probable risks in areas of increased seismicity include life-threatening injuries caused by falling objects and economic loss from damage to structures with low capacity to absorb moderate earthquake shaking.Rational decisions to improve life safety and to reduce losses require sound scientific input.Although hazard-model estimates such as the example shown in Figure 7 can be improved, there is no question that increased hazard accompanies higher levels of earthquake activity.
The potential for catastrophic loss in a larger-than-experienced earthquake must not be neglected, particularly in areas where buildings in service were constructed before modern seismic-design standards were established.Although the 2007 M 5.6 earthquake in San Jose, California, caused little more than cosmetic damage in a heavily urbanized area with strong earthquake codes and decades of preparation, the 2011 M 5.6 earthquake in largely rural central Oklahoma resulted in losses in the millions to unreinforced masonry structures, and the 1986 M 5.7 San Salvador, El Salvador, earthquake took more than 1000 lives, injured another 10,000 people, and made 100,000 homeless.The lesson to be drawn from these examples is that earthquake safety depends not only on the strength of shaking but also on the capacity of structures and society to survive it.
Learning how to prepare for, behave during, and respond after earthquakes is a proven means of improving life safety and reducing losses.Many resources are freely available, including Dart et al. (2011), which provides clear guidance on how to prepare for, survive, and recover from earthquakes.Participation in annual earthquake drills such as the Great Central U. S. ShakeOut also provides an excellent means for increasing earthquake safety awareness at school and in the workplace.

Conclusions
Earthquake activity has undergone a manifold increase in the U. S. midcontinent since 2009, principally in Oklahoma but also in Arkansas, Colorado, Kansas, New Mexico, and Texas.The nature of the space-time distribution of the increased seismicity, as well as numerous published case studies, strongly indicates that the increase is of anthropogenic origin, principally driven by injection of wastewater coproduced with oil and gas from tight formations.Enhanced oil recovery and long-term production also contribute to the rise in seismicity at a few locations.In contrast, the earthquake rate in the regions of highest long-term earthquake hazard in the central and eastern United States has remained steady for at least the last 85 years.
Incorporating rapidly changing rates of earthquakes into hazard models challenges the standard methodologies that are predicated on a time-independent rate of earthquake occurrence.Computing hazard for the next year appears feasible and has value for decisions aimed at managing the hazard or reducing vulnerability.Currently, such forecasts project the immediate past rate of earthquakes forward.
Forecasting more than "last year's earthquakes" will require a deeper understanding of the physical processes and conditions that link perturbations to the earth system to its response in seismic events.A key challenge is to develop a hazard-modeling capability that not only accounts for temporally varying activity rates but also anticipates where induced earthquakes could start or stop in response to changing industrial drivers.Despite these many challenges, there are reasons to be optimistic that the hazard can be managed.
Earthquakes have been known to be a natural part of the landscape in the CEUS since colonial times.The historical record of earthquakes generally begins after settlement by literate people who recorded their observations in journals, letters, 1 U. S. Geological Survey. 2 École Normale Supérieure, PSL Research University.http://dx.doi.org/10.1190/tle34060618.1.

Figure 1 .
Figure 1.2014 National Seismic Hazard Map (NSHM) for the coterminous United States showing peak acceleration with a 2% probability of being exceeded in 50 years.Polygons identify areas where recent seismicity was removed in preparation of the map, corresponding to areas where earthquakes are known to be induced or where seismic activity increased after 2006.After Petersen et al. (2015), Figure 1.Courtsey of U. S. Geological Survey.

Figure 3 .
Figure 3. Annual frequency of occurrence of earthquakes with magnitude greater than or equal to magnitude value for (a) the western area and (b) eastern area.Earthquakes are from the catalog used to produce the 2014 NSHM, divided into time intervals indicated by color key.
Figure 4. Cumulative count of M ≥ 3.5 earthquakes from 1930 through 2014 for the eastern (green) and western areas (blue) defined in Figure 2. The dashed-dotted blue line corresponds to cumulative count in the western area after removal of earthquakes induced by fluid injection at the Rocky Mountain Arsenal, Colorado, between 1962 and 1969.Annual rate of earthquake occurrence derived from a model fit to the cumulative curves shown by the orange and red lines for the eastern and western areas, respectively.

Figure 5 .
Figure 5. Map showing probability of observing at least as many

Figure 6 .
Figure 6.Crustal-movement rate in Oklahoma and surroundings determined from continuous Global Positioning System (GPS) measurements.Uncertainties have a 95% confidence level and include white plus time-correlated noise.A rigid rotation has been removed so that the residuals are deviation from rigid behavior, i.e., strain.Inset compares the cumulative distribution of observed velocities (black bins) with the theoretical χ 2 distribution expected if residuals were distributed normally in two dimensions with a unit variance (red line).The residuals are well described by a random, nontectonic process.Earthquakes of M ≥ 3 in 2013-2014 are shown by orange circles.

Figure 7 .
Figure 7. Uniform hazard maps for (a) 1% and (b) 0.04% probability of exceedance in one year.This example calculation uses a 2014 nondeclustered catalog with M ≥ 2.5, b = 1.5, 5-km smoothing, eight NSHM ground-motion models, and USGS craton M max model (mean M 7).The 5-Hz (0.2-s) spectral accelerations are in units of g (acceleration of gravity) and correspond approximately to the resonance frequency of two-story structures.After Petersen et al. (2015), Figure 16.Courtesy of U. S. Geological Survey.