The 2001–Present Induced Earthquake Sequence in the Raton Basin of Northern New Mexico and Southern Colorado

We investigate the ongoing seismicity in the Raton Basin and find that the deep injection of wastewater from the coal-bed methane field is responsible for inducing the majority of the seismicity since 2001. Many lines of evidence indicate that this earthquake sequence was induced by wastewater injection. First, there was a marked increase in seismicity shortly after major fluid injection began in the Raton Basin in 1999. From 1972 through July 2001, there was one M ≥4 earthquake in the Raton Basin, whereas 12 occurred between August 2001 and 2013. The statistical likelihood that such a rate change would occur if earthquakes behaved randomly in time is 3.0%. Moreover, this rate change is limited to the area of industrial activity. Earthquake rates remain low in the surrounding area. Second, the vast majority of the seismicity is within 5 km of active disposal wells and is shallow, ranging between 2 and 8 km depth. The two most carefully studied earthquake sequences in 2001 and 2011 have earthquakes within 2 km of high-volume, high-injection-rate wells. Third, injection wells in the area are commonly very high volume and high rate. Two wells adjacent to the August 2011M 5.3 earthquake injected about 4.9 million cubic meters of wastewater before the earthquake, more than seven times the amount injected at the Rocky Mountain Arsenal well that caused damaging earthquakes near Denver, Colorado, in the 1960s. The August 2011 M 5.3 event is the second-largest earthquake to date for which there is clear evidence that the earthquake sequence was induced by fluid injection. Online Material: Gutenberg–Richter plots for varying decade-long catalogs.


Introduction
Earthquakes induced by deep underground injection of fluids were first recognized in the 1960s during the Rocky Mountain Arsenal earthquake sequence that was induced near Denver, Colorado (Evans, 1966;Healy et al., 1968). As proposed by Hubbert and Rubey (1959), fluid injection can raise pore pressure within nearby fault zones, thus lowering the effective stress and frictional resistance on faults. The lowered frictional resistance makes earthquake slip more likely. Earthquakes can be induced in this way if there is hydraulic communication between the injection interval and a seismogenic fault zone. The injection of fluids causes a pore pressure increase, which is transmitted into a seismogenic fault zone to induce earthquakes, even though the injected fluid itself may not reach the fault. Study of injection-induced earthquakes has been extensive, including a field experiment in earthquake control (Raleigh et al., 1976). This experiment demonstrated that fluid pressure controlled the rate of earthquake occurrence in part of the Rangely Oil Field in northwestern Colorado. Raleigh et al. (1976) found that when they increased the pressure within a portion of the field, earthquake rates increased; and, when it was decreased, the earthquake rate dropped. This was the first demonstration of controlling earthquakes by adjusting pore pressure at depth.
Since these studies in the 1960s and 1970s, many other earthquakes have been identified as having been induced by fluid injection, including earthquakes in Ashtabula, Ohio (Seeber and Armbruster, 1993;Seeber et al., 2004), Dallas-Fort Worth, Texas (Frohlich et al., 2011), and El Dorado, Arkansas (Cox, 1991). There was a recent spate of earthquakes believed to be induced, including the 2011 M 5.7 Prague, Oklahoma, earthquake (Keranen et al., 2013;Llenos and Michael, 2013;Sumy et al., 2014), the 2011 M 4.0 Youngstown, Ohio, earthquake (Kim, 2013), the Paradox Valley, Colorado, earthquakes (M max 4.4; Block et al., 2014), and the 2011 Guy-Greenbrier, Arkansas, earthquake sequence (M max 4.7; Horton, 2012). These are part of a larger trend of increased seismicity in the central and eastern United States since 2001, much of which is believed to be induced by wastewater injection (Ellsworth, 2013). All of these earthquakes occurred in areas with little or no previous seismicity, and, assuming that they are induced, this indicates that earthquake hazard needs to be reassessed in these areas as well as other areas where wastewater is injected.
It is often difficult to determine with certainty whether a specific earthquake was induced or not. Many kinds of data are needed to describe the physical state of the focal volume; for example, initial stress state, location and static strength of faults, hydrogeologic structure, fluid injection volumes, injection rates, injection pressures, and precise earthquake locations. Without these kinds of data, it is very difficult to determine whether earthquakes are induced or natural. Seismologists are often left looking for correlations in time and space between injection activities and earthquakes. To deal with this problem, Davis and Frohlich (1993) offered a set of seven criteria that help to determine whether or not earthquakes are induced. They were intended as guidance to help determine whether earthquakes were induced and not as a decision tool. Although these criteria provide a good framework, the lack of data may preclude answering one or more of the questions, thus limiting the utility of the criteria. Additionally, studies of induced earthquakes in the past 20 years have shown that these criteria can be too restrictive, such that some clearly induced earthquakes would show attributes inconsistent with an affirmative answer to some of the Davis and Frohlich (1993) criteria.
In this work, we focus on an extended earthquake sequence in the Raton Basin of southern Colorado and northern New Mexico (Fig. 1). Historical and instrumental data show that the Raton Basin had a low level of seismic activity until August 2001. This changed with an earthquake sequence that started in August 2001, followed by increased seismicity in the vicinity of the initial sequence that has continued to the present. This includes an M 5.3 earthquake that occurred on 23 August 2011.
Here, we address the question of whether this earthquake sequence is related to wastewater injection and find multiple lines of evidence supporting this hypothesis. To evaluate this, we first examine the industrial activities and seismicity of the broader region and explore the relationship between them. Statistical analysis of the earthquake rate change that occurred in 2001 indicates that it is highly unlikely the rate change could arise from random fluctuations of the ambient seismicity, given a constant background rate. These earthquakes are limited to the area of fluid injection, they occur shortly after major fluid injection activities began, and the earthquake rates track the fluid injection rates in the Raton Basin. We also examine earthquake sequences in 2001 and 2011 in detail. These sequences lie very close (≤ 2 km) to high-volume, high-injection-rate wells, which shows that these sequences, specifically, were induced by nearby wastewater injection.

Regional Tectonics and Deformation
The Raton Basin is a coal-bearing sedimentary basin situated along the Colorado-New Mexico border, between the western edge of the Great Plains to the east and the Sangre de Cristo Range and Rio Grande rift to the west. It is approximately 150 km long (north-south) and 75 km wide at its maximum. On a regional scale, the Raton Basin lies within a broad crustal zone of east-west extension astride the Rio Grande rift (Heidbach et al., 2008;Berglund et al., 2012).
Geologic mapping within the Raton Basin reveals little evidence for faulting. Maps commonly show that the basin is bounded by thrust faults on the western edge of the basin and a west-dipping normal fault striking northwest that runs the length of the eastern side of the basin (Fig. 1). In addition to the basin-bounding faults identified by many maps, one map also identifies faults lying within the basin (Robson and Banta, 1987). This map includes two normal faults striking to the northeast within the southeastern portion of the Colorado side of the basin. These faults are believed to be buried in the Precambrian basement and are not seen to outcrop within the Raton formation, which lies at the surface above these mapped faults (upper Cretaceous/lower Paleocene). We could not identify any seismicity associated with these faults. They are also not found within the U.S. Geological Survey (USGS) Quaternary Fault and Fold Database (see Data and Resources), so we do not believe them to be active. The closest faults that are known to be active within the Quaternary Period are west of the basin in the Rio Grande Rift (Fig. 1).

Industrial Activities in the Raton Basin
Coal mining began in the Raton Basin in 1862 and continued through 2002, although production greatly decreased in the early 1960s. In 1994, energy companies began producing coal-bed methane from the same formations. Methane production was initially in the Colorado portion of the basin and expanded to New Mexico in 1999. The production formations are the Raton, Vermejo, and Trinidad formations, with production wells typically drilled to the top of the underlying Pierre formation, which ranges from 200 to 800 m depth.
Along with methane, there is considerable formation water produced at the same time. Some of this formation water is injected deep underground for disposal. Wastewater injection in Colorado began in 1994, expanding into New Mexico in 1999 and is primarily in the Dakota formation, a buff, conglomeratic sandstone (Johnson, 1969), with injection intervals ranging between 1250 and 2100 m below the ground surface, depending on location in the Raton Basin (Colorado Oil and Gas Information System [COGIS]; see Data and Resources).
The Dakota formation has a large lateral extent. The closest outcrops of the Dakota formation to the Raton Basin are 40 km north of Trinidad and 5 km south of the southwestern limit of the Raton Basin (Johnson, 1969;New Mexico Bureau of Mines and Mineral Resources, 2003). East-west hydraulic continuity in the Dakota formation is believed to be on the order of 80 km in the Raton Basin, with the exception of the northwestern portion of the basin, which appears to be hydrologically isolated (Nelson et al., 2013). There are no disposal wells in the northwestern portion of the basin because the produced formation water there meets water-quality standards for surface discharge.
Many of the hydrologic units within the Raton Basin are underpressured, including the Dakota sandstone (Johnson and Finn, 2001;Nelson et al., 2013). On average, the hydraulic head within the Dakota unit lies approximately 500 m below the surface, which is approximately 4.9 MPa underpressured (Nelson et al., 2013). Because the Dakota unit is underpressured, injection throughout much of the Colorado portion of the Raton Basin can be done with gravity feed (D. Onyskiw, personal comm., 2012). Of 21 injection wells in the Colorado portion, only 5 have ever injected under pressure, 2 of which have been operating under gravity feed since 2005. In the case of gravity feed, even when there is no injection pressure, the weight of the water column in the well bore still causes a stress change at depth. This is not the only source of stress changes due to fluid injection. Other sources of stress change come from the redistribution of fluids at depth upon injection and poroelastic effects whereby the medium is forced to accommodate the injected fluids. Information on injection pressures on the New Mexico side of the basin indicates that the wells have always injected under gravity feed.
The Dakota sandstone lies anywhere between 800 and 1400 m below the bottom of the Trinidad formation (the lower production unit; Johnson and Finn, 2001). Given that both the production and injection formations are believed to be hydrologically isolated from adjacent geologic units (Nelson et al., 2013), we expect little or no hydrologic commu-nication between them. Likewise, without a fluid pathway (e.g., a fault), we would expect little communication between the injection formation (Dakota) and underlying geology. This includes the Precambrian basement, where the majority of the earthquakes have occurred. Between the Dakota hydrologic unit and the basement, there is a Jurassic hydrologic group that is primarily sandstones, below which is a Permian-Pennsylvanian hydrologic unit composed of limestones, sandstones, and shales (Geldon, 1989). This is underlain by Miocene metamorphics and volcanics, which lie on top of the Precambrian basement. Given the numerous seals between these hydrologic units, communicating hydraulic pressures to depth would require a connecting fluid pathway, like a fault.

Wastewater Injection in Colorado
The Colorado Oil and Gas Conservation Commission (COGCC) regulates wastewater injection wells in Colorado and maintains the COGIS online database with all the records for these wells. Wastewater injection began in the Raton Basin in November 1994 at the Cottontail Pass well in Colorado (Table 1). Two months later, in January 1995, two additional injection wells came online: Apache Canyon 10 and Apache Canyon 19. The field rapidly expanded from 1997 to 2001, with eight additional injection wells coming online. Six of these wells are located in the eastern portion of the production field. Since 2001, 10 more injection wells have opened, giving a total of 21 presently active injection wells in the Colorado portion of the Raton Basin.
Prior to the rapid expansion of production and injection activities in the Colorado portion of the Raton Basin (in mid- 2000), injection rates remained under 600;000 barrels=month ( Fig. 2). These rates began to rise with the expansion of the field; and, since mid-2000, monthly injection rates have ranged between 600,000 and 1:9 million barrels=month. With increasing production, there has been a corresponding increase in the number of disposal wells in the area. Because barrels are the standard measure of volume used within the oil and gas industry, we have elected to use barrels as the measure of volume in this article instead of the metric measure of cubic meters. For reference, there are approximately 6:29 barrels=m 3 .

Wastewater Injection in New Mexico
In New Mexico, the Oil Conservation Division of the Energy Minerals and Natural Resources Department regulates wastewater injection wells and maintains an online database with records pertaining to wastewater injection wells.
Injection data prior to June 2006 are unavailable for New Mexico. We are able to get a field-wide sense of injection prior to June 2006 by using produced-water data as a proxy for injection, but associating injection rates and total volumes with individual wells is not possible. The production records for the Raton Basin in New Mexico extend back to the beginning of production in the Raton Basin in 1999, and these data include produced-water rates. We believe that water production is an appropriate proxy for injection for two reasons: (1) for the period when both injection and production values are available, the totals match each other fairly well (Fig. 3a); (2) essentially 100% of all produced water is injected in this area by requirement of the landowner (personal comm. with a local operator, 2012). Based on well permitting and drilling information, there was an apparent lag between the initiation of production and the initiation of injection in the New Mexico portion of the basin. The first injection well that was completed was VPRA042, with a completion date of May 2000. It is unclear what happened to the water prior to the completion of this well, but certainly production rates prior to this are not an appropriate proxy for injection rates in the New Mexico portion of the Raton Basin. Given that these rates are quite small, we do not expect it to strongly influence later analysis.
Injection rates in New Mexico, inferred from producedwater rates, slowly climbed from the start of injection in 1999 until early 2004 when rates began to level off (Fig. 3). Since Mexico. Monthly wastewater injection rates and produced-water rates follow each other reasonably well for much of this time period, but differences can be as large as 580,000 barrels in an individual month. It is unclear why at some points the rates match each other very well and other times there are large differences. (b) Number of earthquakes per six months (bars/stems) and monthly injection rates in the New Mexico portion of the Raton Basin (line). The earthquake catalog is complete for earthquakes M ≥ 3:8. One can see that the earthquakes start shortly after injection rates increase. No information is available on injection in New Mexico before June 2006, so produced-water data is used as a proxy. The color version of this figure is available only in the electronic edition.
2004, the median injection rates have been approximately 1:2 million barrels=month.

Seismicity of the Raton Basin
Here, we discuss the seismicity of the Raton Basin in a broader context. A more detailed discussion of the earthquake sequences that began in August 2001 and August 2011 is provided below in the Case Study: The August-September 2001 Earthquake Sequence and Case Study: The August-September 2011 Earthquake Sequence sections.

Historical Seismicity
To provide historical context, we review the historical seismicity for a broad region centered on the Raton Basin. We define this region as between 103°and 106°W longitude and 36°and 38°N latitude (Fig. 1).
Our review of the earthquake history of the region relies on sources published in the scientific literature (Northrop and Sanford, 1972;Stover and Coffman, 1994;Kirkham and Rogers, 2000;Sanford et al., 2002) and annual summaries (see Data and Resources) of U.S. seismicity (U.S. earthquakes) published by the Department of Commerce and covering the years from 1928 through 1975. Studies of earthquakes in the preinstrumental period commonly rely heavily on newspapers and other written accounts. Settlements in the Colorado portion ( Fig. 1) of the basin include the towns of Weston (founded in the late nineteenth century) and Cokedale (founded in 1906). There are no towns in the New Mexico portion of the basin where coal-bed methane is being produced today.
Our search of historical reports through 1970 did not identify any earthquakes that could be associated with the Raton Basin. Given that The Chronicle-News (based in Trinidad, Colorado, 10 km east of the coal-bed methane field) has been in continuous publication since 1877 and that the 10 August 2005 M 5.0 earthquake was widely felt at modified Mercalli intensity (MMI) IV in Trinidad (USGS "Did You Feel It?" system; see Data and Resources), we expect historic accounts to be complete throughout the Raton Basin at M ≥ 5:0 since 1877.
The earliest reported earthquake in this general region occurred on 14 June 1878 and was reported to have broken windows at Cimarron, New Mexico (Stover and Coffman, 1994). An earthquake with a maximum reported MMI of V occurred on 13 August 1924, about 20 km east of Wagon Mound, New Mexico (Northrop and Sanford, 1972). On 12 July 1936, an earthquake with a magnitude of 3.4 occurred near the New Mexico-Oklahoma border (Neuman, 1938). An MMI V earthquake was felt near Cimarron, New Mexico, on 4 August 1952 (Murphy and Cloud, 1952). Later that year (7 October 1952), an earthquake with MMI V occurred, with strongest intensities observed at Antonito, Colorado (Murphy and Cloud, 1952;Kirkham and Rogers, 2000). Locations of these towns are shown in Figure 1.

Spatiotemporal Evolution of Instrumentally Recorded Seismicity
For the purposes of examining the instrumentally recorded seismicity, we narrow our study region so it contains little more than the Raton Basin. We nominally select a box with the following bounds: 36.7°N to the south, 37.6°N to the north, 105.2°W to the west, and 104.4°W to the east. This covers the boundaries of the Raton Basin fairly well, only omitting small portions of the basin at its edges (Fig. 4). This region extends 20 km or more beyond every injection well in the area, so it is an appropriate size to examine the relationship between earthquakes and wastewater injection wells.
The earthquake catalog is primarily composed of earthquakes listed in the Advanced National Seismic System (ANSS) Comprehensive Catalog (ComCat; see Data and Resources). ComCat contains earthquake hypocenters, magnitudes, phase picks, moment tensor solutions, macroseismic information, tectonic summaries, and maps, which are produced by contributing seismic networks. ComCat replaces the ANSS composite catalog. This catalog is supplemented by catalogs from the International Seismological Centre (ISC) Bulletin, the Los Alamos Scientific Laboratory seismic array (1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984), a catalog produced by Meremonte et al. (2002) to study the August 2001 earthquake sequence, and a catalog we prepared for this study that examines the August 2011 earthquake sequence. Earthquake locations in the ComCat, ISC Bulletin, and the Meremonte et al. (2002) catalogs have been refined using an updated velocity model. Further details on the location methodologies can be found in Appendix A. Analysis of the complete catalog indicates a magnitude of completeness of M c 3.8 from 1970 to present (Appendix B).
Before August 2001, seismicity in the Raton Basin was widely distributed and infrequent (Fig. 4a). Seismicity was primarily found near the northeastern margin of the basin. Additional seismicity could also be seen outside the basin to the northeast and to the west.
Since 2001, the behavior of seismicity in the Raton Basin has changed significantly. Earthquakes are far more frequent than they were in the preceding 30 years. Sixteen M ≥ 3:8 earthquakes occurred in the period August 2001-2013, as opposed to one in the preceding 31 years and 7 months. This represents a 40-fold increase in earthquake rate, from 0:032 earthquakes=year to 1:27 earthquakes=year (Fig. 5). The spatial distribution of the seismicity also changed in 2001. Most of the recent seismicity is concentrated near the center of the Raton Basin, whereas earlier seismicity was located on the periphery of or outside the basin (Fig. 4).

Notable Earthquakes and Earthquake Sequences
Three notable earthquake sequences have occurred since August 2001: the August-September 2001 earthquake sequence; the August-September 2005 earthquake sequence that included an M 5.0 event; and the August-September 2011 earthquake sequence that included the largest recorded earthquake in the area, the 23 August 2011 M 5.3 earthquake. All of these sequences are located near the center of the Raton Basin (Figs. 4, 6, and 7).
The 2001 earthquake sequence occurred in the Colorado portion of the basin and included an M bLg 4.0 event that occurred on 4 September 2001 and an M bLg 4.5 event that occurred on 5 September 2001. This earthquake sequence marks the beginning of a higher seismicity rate in the Raton Basin. To better understand these events, the USGS responded by deploying a 12-station temporary seismic network. From precise hypocenter locations, Meremonte et al. (2002) identified a plane of seismicity striking to the northeast and dipping steeply to the southeast (approximately N37°E). Our relocations are consistent with this finding (Fig. 6). Focal mechanisms of the two largest earthquakes in this sequence indicate east-west extension, although the inferred fault planes do not dip as steeply as the plane imaged by the well-located hypocenters (Fig. 8). Generally, this evidence is consistent with the regional extensional tectonics (Heidbach et al., 2008;Berglund et al., 2012). More details about this earthquake sequence can be found in the Case Study: The August-September 2001 Earthquake Sequence section.
The 2005 earthquake sequence includes the second-largest instrumentally recorded earthquake in the Raton Basin, the 10 August 2005 M 5.0 earthquake (Fig. 4b). This earthquake sequence was located in the New Mexico portion of the basin, about 10 km south of the state line. No temporary instruments were deployed at this time, limiting station coverage to permanent stations (of which there were nine stations within 200 km) and a maximum azimuthal gap of 146°, using stations up to 350 km from the event. There is no apparent structure to the seismicity in this earthquake sequence, probably owing to earthquake locations with uncertainties of about 15 km or more. An M b 4.1 foreshock occurred 6 s before the mainshock, and an M L 3.0 aftershock occurred 16 min later. There were no other M ≥ 3 earthquakes within one month of these earthquakes. The magnitude of completeness is M c 3:0 (Appendix B). The moment tensor for the mainshock indicates normal faulting on a nodal plane striking north-south (Fig. 8). As already noted, this is consistent with the regional east-west extension (Heidbach et al., 2008;Berglund et al., 2012).
The 2011 earthquake sequence began on August 21 and included an M 4.6 earthquake on August 22, which was followed by the M 5.3 mainshock 6 hr later on 23 August. The largest aftershock was an M 4.0 earthquake that occurred later in the day on 23 August. The 2011 earthquake sequence immediately abuts the 2001 earthquake sequence, extending to the southwest, with virtually no spatial overlap between the epicenters of the 2001 and 2011 sequences (Fig. 7a). Like the 2001 sequence, the epicenters form a steeply dipping tabular structure striking northeastward. As with the August 2001 and August 2005 sequences, mainshock focal mechanisms are consistent with normal faulting on northeast-striking structures (Fig. 8). The USGS responded to the earthquake sequence by deploying a four-station temporary seismic network to record the aftershocks. An extended discussion of this earthquake sequence can be found in the Case Study: The August-September 2011 Earthquake Sequence section.

The Relationship between Fluid Production and Seismicity in the Raton Basin
We do not believe that the production of gas or water is directly related to the earthquakes in the Raton Basin. Given that the oil production in the area is minimal, and water is far heavier than natural gas, the extraction of water will have the largest stress effect. Examining the Colorado portion of the basin through June 2012, we find that the maximum amount of produced water in 2012 in a 5 km × 5 km square (25 km 2 ) within 20 km of the 2011 earthquake sequence was approximately 66 million barrels, or approximately 10 million cubic meters. Assuming the production is uniform across the 25 km 2 , this would give a withdrawal of approximately 40 cm of fluid across the area, which is equivalent to a 4 kPa stress change. Because fluid withdrawal is likely to be uneven, we consider a factor of 5 uncertainty in the stress change, giving a 20 kPa stress change, which is of the same order of magnitude as the minimum threshold for natural earthquake triggering.
Studies of static stress triggered earthquakes typically only see earthquakes triggered at static stresses of 10 kPa or more (Hardebeck et al., 1998), although Ziv and Rubin (2000) saw triggering at smaller stresses given specific criteria. Earthquakes dynamically triggered by the passage of teleseismic waves are typically triggered at 30 kPa or greater (Gomberg and Johnson, 2005), but triggering stresses have been seen to be as small as 5 kPa in areas particularly susceptible to earthquake triggering (Brodsky and Prejean, 2005). Because the stress changes from the production of water are at the lower end of stresses where earthquakes are triggered by natural processes (i.e., earthquakes are not triggered in the vast majority of locations experiencing these stress changes), we find it unlikely that fluid production contributes to the occurrence of the earthquakes in the Raton Basin. Additionally, the maximum stress change that we observe is located 15-20 km to the northwest of the seismicity in the 2011 sequence.  Davis and Frohlich (1993) to this earthquake sequence but were equivocal as to whether the earthquake sequence was related to the wastewater injection in the region. It is worth noting, however, that Meremonte et al. (2002) stated that if the earthquake sequence were natural, they would expect the seismicity to tail off and return to the lower background seismicity rate. Instead, the seismicity rate in the Raton Basin has continued at a much higher rate than the pre-2001 period. In this section, we examine the seismicity in the Raton Basin in its entirety with respect to wastewater injection activities across the Raton Basin. We find strong evidence that this earthquake sequence is induced by fluid injection in the area.  Table 2). In this section, we examine the likelihood that the change in earthquake rate that we observe could occur as a random fluctuation in earthquake rate given a constant, long-term background rate of seismicity (i.e., earthquake occurrence follows a uniform Poisson process).
Although we are confident that the catalog is complete for earthquakes M ≥ 3:8 (Appendix B), we use a catalog of the M ≥ 4 earthquakes in the Raton Basin to be even more certain that no events are missing. In this catalog, we find 1 M ≥ 4 earthquake from 1970-July 2001 and 12 M ≥ 4 earthquakes from August 2001-2013 (Table 2).
With this catalog, we apply the binomial test to determine how likely it is that the earthquake rate change that we observe could be produced by random variation in earthquake rate, given a constant, background seismicity rate. The binomial test takes the background rate of seismicity and, given this rate, determines how likely it is to get X number of earthquakes within a time period Y. One must be careful with the binomial test, because it assumes that earthquake behavior is random and not clustered as earthquakes commonly are (e.g., aftershocks). This can be addressed by declustering the catalog with respect to time and space. Declustered catalogs have been shown to be Poissonian (Gardner and Knopoff, 1974), and thus the binomial test is appropriate. Given that there are clear foreshock-aftershock sequences in our catalog (Table 2), we decluster our catalog using the Gardner and Knopoff (1974) method and identify four different sequences of potentially related earthquakes (Table 2). In the declustered catalog, we count each cluster as only one earthquake. This yields a catalog of seven earthquake sequences.    Filled circles with thick borders are earthquakes that occurred before the temporary network was installed and have location uncertainties of 2:3 km. Thin, unfilled circles are earthquakes that occurred after the temporary network was installed and have relative location uncertainties of 300 m. Thick, unfilled circles are the M 5.3 mainshock, foreshocks, and early aftershocks located using the polarization method. Relative uncertainty for polarization locations is 2 km. The ×'s are early foreshocks and aftershocks for which no magnitude has been computed. (a) A broader picture of seismicity during this period. Squares are injection wells sized according to the cumulative volume of fluid injected through August. The lines show mapped faults in or near the Raton Basin (Johnson, 1969;Robson and Banta, 1987;Scott and Pillmore, 1993;USGS Quaternary Fault and Fold Database [see Data and Resources]). More detailed mapping has been conducted in the Colorado portion of the basin than on the New Mexico side, hence there are more mapped faults. We do not expect there to be a significant difference in the faulting on either side of the border, we just have more detailed knowledge of faults on the Colorado side. The large box in (a) is the study area shown in (b), and the smaller box is the study area for the 2001 earthquake sequence in Figure 6 With the declustered catalog, we treat the background rate as seven earthquakes occurring in the 44-year study period. Assuming earthquakes occur randomly in time, we can determine the likelihood that an individual earthquake would occur in a period of time as simply the length of that period divided by length of the entire study (1970-2013, or 44 years). Here, we are interested in whether the earthquake rate change that occurred at the time fluid injection commenced could be produced by random fluctuations in the background rate. Thus, the likelihood that any individual earthquake would occur in the coinjection period (November 1994-December 2013 is approximately 44%. There were six independent earthquakes within this time period and one in the preinjection period. We use combinatorics to determine the likelihood that six or seven earthquakes would occur in the coinjection period (Appendix C). This analysis yields that it is unlikely (3% probability) that the observed seismicity rate change can be the result of random fluctuations in earthquake rate, given a constant, long-term background seismicity rate.
Although large earthquake rate changes have been seen to be correlated to industrial activities (e.g., Healy et al., 1968;Horton, 2012), natural causes for earthquake rate changes have also been observed. The most common causes of prolonged, large earthquake rate changes are fluids and fluid movement in volcanic, geothermal, and hydrothermal systems (e.g. Hill, 2006). Thus, this earthquake rate change alone is not sufficient evidence to demonstrate that these earthquakes are caused by fluid injection.

Spatial Relationship between Seismicity and Wastewater Injection Wells
Seismicity before 2001 was concentrated at the edges of the Raton Basin (Fig. 4a). From 2001 forward, most of the seismicity is found in the center of the basin, in the area of the injection wells (Fig. 4b,c). In fact, the earthquake rate change appears to be solely coming from the area of the wells. No earthquakes larger than the magnitude of completeness (M c 3:8) are found outside the area of the wells before or after injection began, thus the increased rate occurs in the center of the basin near the wells. The broadscale coincidence of the location of increased seismicity and location of the increased injection activity suggests that the injection activities are related to the increased seismicity rate.
Since the start of injection, one can see that most of the seismicity is located quite close to active injection wells, whereas earlier seismicity is not (Fig. 4a-c). The majority of this seismicity was located without the benefit of dense seismometer deployments, giving high epicentral uncertainties (2σ of approximately 15 km). Thus, the absence of a precise spatial correlation for these earthquakes with wells does not mean that they are not related to the wells. Additionally, there is evidence that induced earthquakes can occur at distances of greater than 10 km from causative injection wells (Herrmann et al., 1981;Horton, 2012;Ellsworth, 2013;Keranen et al., 2014).
Although for much of this extended earthquake sequence, location uncertainties are large, when uncertainties are lower (as with the 2001 and 2011 sequences), we see that their epicenters lie within 2 km of high-volume, high-injection-rate wells (Figs. 6 and 7). The well-located seismicity is situated at depths extending downward from approximately 1 km below the bottom of the injection wells.

Fluid Injection Volumes, Injection Rates, and Seismicity Rates
In this section, we examine the fluid injection history on both sides of the Colorado-New Mexico border and compare this history to the temporal variations in seismicity. Overall, we find that both earthquake rates and fluid injection rates have remained high over the last 10 years. Injection rates are known in Colorado extending back to the beginning of methane production in the field (1994) Raton Basin in Colorado. From the start of injection in November 1994 through July 2001, no earthquakes were detected in the Colorado portion of the Raton Basin. During much of this period, fluid injection rates in this area were fairly low. From 1994 through February 2001, the maximum rate was 725;000 barrels=month, and the median rate was 300;000 barrels=month ( Fig. 2). From the start of 2001 until the first earthquakes in this sequence were felt (August 2001), injection rates in the Colorado portion of the field dramatically increased, rising from a median rate of 500;000 barrels=month in the year 2000 to 1:2 million barrels=month in August 2001 (Fig. 2). The earliest earthquakes were located in the eastern portion of the gas field, shortly after six wastewater injection wells were put into operation. The level of seismicity, although quite variable, has been fairly high in the Colorado portion of the field since injection increased in March 2001 (Fig. 2). Examining the cumulative injection volumes and number of earthquakes in the Colorado portion of the field, we can see, from a broad perspective, that total injection volumes and the number of earthquakes roughly track each other (Fig. 9a). There is a fairly constant injection rate on the Colorado side of the border. Likewise, the earthquake rate is roughly constant with considerable temporal clustering.
Raton Basin in New Mexico. We use water production in New Mexico as a proxy for injection from 1999 to May 2006. Production rates slowly increase from the start of production in 1999 through mid-2005. At the time of the first M ≥ 3:0 earthquakes detected in the New Mexico portion of the Raton Basin (January 2002), injection rates were low (∼250;000 barrels=month). By the time earthquakes became more frequent (early 2004), injection rates had risen to 750;000 barrels=month ( Fig. 3). Since mid-2005, injection rates have remained largely constant (median: 1:2 million barrels=month). Similarly, earthquake rates have been generally constant since this time and closely track the injection rates (Fig. 9b).
The total volumes injected in New Mexico and Colorado since 1999 are comparable, with ∼110 million barrels on both sides between June 2006 and September 2013. During the time period of the increased seismicity (2001-present), Colorado and New Mexico also had similar numbers of M ≥ 3 earthquakes (54 and 45, respectively).

Mechanics of Injection-Induced Seismicity
The mechanics of injection-induced earthquakes have been long established (e.g., Healy et al., 1968;Raleigh et al., 1976). Injection increases the pore fluid pressures in the fault, reducing the normal stress and hence frictional resistance, and in turn facilitating slip under the ambient tectonic shear stress (Hubbert and Rubey, 1959   Even in an underpressured, extensional system like the Raton Basin, where fluids can typically be injected with no wellhead pressure, earthquakes can still be induced. In these cases, pore pressure still increases at the injection interval as a result of the excess hydraulic head in the well bore. In fact, the fluid pressures from hydraulic head alone have been observed to both hydraulically fracture and cause slip on faults in wells in underpressured, extensional environments (Stock et al., 1985;Moos and Morin, 1991). When the increased fluid pressure from the water column alone is large enough to induce failure, fluid injection need not require wellhead pressure to cause slip. The critical stress hypothesis (Townend and Zoback, 2000), which states that the crust is everywhere near failure, also supports the notion that small stress changes can be enough to induce an earthquake.
In addition to stresses coming directly from the wellhead, there are other, more broadscale stress changes that come from fluid injection. With the injection of fluids, there will be significant fluid movement to accommodate the increase in fluid volume. This fluid redistribution will change the stress field. Poroelastic effects are also important, as pore spaces will get filled and additional pores may be opened or closed due to the increased fluid pressures.

Case Study: The August-September 2001
Earthquake Sequence Description of the Earthquake Sequence In response, members of the USGS Geologic Hazards Team from Golden, Colorado, deployed a temporary array of 12 seismometers to record the ongoing earthquake sequence. The stations were deployed beginning 8 September 2001 (Meremonte et al., 2002), after the two largest earthquakes in the sequence. A total of 39 earthquakes, that occured between 10 September 2001 and 15 October 2001, were located by Meremonte et al. (2002). The hypocenters of these earthquakes define a northeast-striking plane of seismicity spanning approximately 7 km in length and dipping steeply to the southeast.
We relocate this seismicity using an updated velocity model. Our locations for the 2001 hypocenters are much the same as those of Meremonte et al. (2002). Details on the location methodology can be found in Appendix A. Our locations also define a northeast-striking, steeply dipping fault between 1.5 and 6 km depth, centered below the Wild Boar UIC Class II disposal well (Fig. 6b). The plane imaged by the seismicity runs parallel to a fault mapped by Robson and Banta (1987) but lies approximately 5 km to the northwest of it (Fig. 6a).
Uncertainties for the relative earthquake locations when the temporary network was running are approximately 200 m (2σ). Absolute uncertainties (2σ) are approximately 2 km. Earthquakes that occurred before 10 September 2001 were located using the regional seismic network, which has interstation spacing of hundreds of kilometers. Consequently, the 2σ location uncertainty of these earlier earthquakes is approximately 15 km. The depth to basement beneath the Wild Boar well is estimated to be between 2 and 3 km, based on seismic reflection sections (personal comm. with a local oil and gas operator, 2012). Thus, it appears that the earthquakes during the late stages of the sequence were in the lower sedimentary section and uppermost basement (Fig. 6b). Regional moment tensors were computed for the two largest earthquakes in this sequence by Saint Louis University (SLU; Herrmann et al., 2011). Both earthquakes are normal-faulting events striking approximately northeast (Fig. 8a,b). This is consistent with the regional tectonics of the area, which are east-west extension, dominated by the Rio Grande rift system.

Evidence that the 2001 Earthquake Sequence was Induced by Wastewater Injection
The August 2001 earthquake sequence is effectively centered below the Wild Boar injection well (Fig. 6) with many epicenters within hundreds of meters of the injection well. Most events in the 6 km lineation of seismicity lie between 1 and 3 km below the injection depth (Fig. 6b). At the time of the 2001 sequence, there were four other active injection wells within 10 km of the earthquake sequence: PCW (2 km), Sawtooth (4.5 km), Long Canyon (8 km), and La Garita (8.3 km).
The injection rates at Wild Boar in the months leading up to the earthquake sequence were higher than any of the other wells listed above. Injection rates at Wild Boar routinely exceeded 200;000 barrels=month in 2001, whereas rates at the PCW well averaged 150;000 barrels=month and rates at Sawtooth, Long Canyon, and La Garita never exceeded 93,000, 52,000, and 36;000 barrels=month, respectively. Wild Boar also began injecting in August 2000, one year before felt earthquakes began. PCW, the other high-injection-rate, high-volume well began injection three years before Wild Boar without any detected earthquakes. Sawtooth, Long Canyon, and La Garita also began injection shortly before the first felt earthquakes (April 2000, August 2001, and April 2001. All of these wells have always injected under gravity feed only.
Given that Wild Boar is the highest-rate well, began injecting shortly before earthquakes started in the area, and is the closest well to the earthquakes, it is the most likely well to have induced the earthquakes. Because PCW is a highinjection-rate well and the other wells began injecting shortly before the earthquakes started, we cannot rule out these wells as contributing to inducing the earthquakes.

Case Study: The August-September 2011
Earthquake Sequence The largest earthquake documented in the Raton Basin is an M 5.3 earthquake that occurred on 23 August 2011. This earthquake was preceded by several foreshocks, including an M 4.7 earthquake approximately 6 hr before the mainshock. The aftershock sequence decayed quickly, with most aftershocks ending within approximately one month of the mainshock. The USGS Geologic Hazards Science Center in Golden, Colorado, installed a four-station, temporary seismic network the day after the mainshock to record and locate the aftershocks. It would have been desirable to install addi-tional instruments, but this was not possible due to the need to respond to the M 5.8 central Virginia earthquake that occurred later on the same day.
Earthquake location quality varies for this earthquake sequence. Because the temporary seismic network was not installed at the time of foreshocks, mainshock, and early aftershocks, we applied standard location techniques with the permanent regional stations. Location uncertainties were similar to those of the preceding 10 years (∼15 km). Fortunately, we were able to reduce the uncertainty in epicenters of these early temblors to approximately 0:9 km. We do this by measuring the P-wave polarization and S-P time at station T25A and finding the intersection of this ray with two fault segments defined by the well-located aftershocks of the 2011 earthquake. Smaller, early aftershocks were located using a hypocentral decomposition approach and have 2σ epicentral uncertainties of 2:3 km. Aftershocks recorded by the temporary network have 2σ relative epicentral uncertainties of 300 m, relative depth uncertainties of 250 m, and absolute uncertainties of 2 km. Further details on the methods used can be found in Appendix A.

Foreshocks
The 23 August 2011 M 5.3 Trinidad earthquake was preceded by a number of foreshocks. Within ComCat, three earthquakes were identified in the 24 hr preceding the mainshock: the previously mentioned M 4.7 earthquake and events of M L 2.9 and 3.0. To identify additional foreshocks, we manually scanned the continuous seismic records at station T25A for July and August 2001. This analysis revealed a foreshock sequence with 36 events that began with an M ∼ 1:1 event at 10:05 UTC on 21 August, which is the smallest foreshock that we detected. Although there are not enough earthquakes to conduct a formal analysis, Gutenberg-Richter statistics suggest that the magnitude of completeness is 2.1 or greater for the period during which we manually scanned the waveforms. For 2001-2013, the overall magnitude of completeness is approximately 3.0 (Appendix B).
The early foreshock sequence (prior to the M 4.7) is largely centered on the first foreshock (Fig. 7d). The earliest foreshocks were in the north. Following the M 4.7 earthquake, foreshocks were largely concentrated to the south. The M 4.7 is very close (< 2 km) to two high injection-rate wells, VPRC 14 and VPRC 39. The later foreshocks extend further south than the majority of the aftershocks (Fig. 7d).
Many of the largest events produced prominent phases on the seismograms that are diagnostic of earthquakes occurring at shallow depth. We were able to confirm this shallow depth of the foreshocks by comparing seismograms recorded at T25A (Fig. 10a) with synthetic seismograms computed for varying depths. We find the best match for surface waves (which are the waves most diagnostic of earthquake depth) and the P and S arrival times (which are also diagnostic of depth) for the M L 2.9 foreshock at 13:52 UTC on 22 August 2011 is with the synthetic seismograms for hypocentral depths of 3-4 km. We find hypocentral depths of 3-4 km are also appropriate for the M L 3.0 earthquake foreshock that occurred on 23 August 2011 at 02:48 UTC.

The Mainshock
The M 5.3 mainshock occurred on 23 August 2011 at 05:46 UTC. It is located near the center of the foreshock sequence, approximately 5 km north of the southern end of the sequence and 10 km south of the northern end (Fig. 7). Regional waveform modeling indicates that it is a normalfaulting earthquake, with fault planes striking to the northnortheast and dips ranging between 40°and 50° . The lineation of seismicity of both the foreshocks and the aftershocks is consistent with the strikes of the two nodal planes defined by the moment tensor. The lineation is also roughly consistent with the strike of a normal fault identified by Robson and Banta (1987) (Fig. 7a,b). Considering that the closest part of the seismicity is approximately 4.5 km from the mapped fault, we do not expect that this similarity reflects a direct physical connection between the two structures. Both structures, though, are consistent with the ambient state of stress. The focal mechanism is also consistent with the regional tectonic setting of east-west extension (Heidbach et al., 2008;Berglund et al., 2012) Standard arrival-time location techniques place the mainshock hypocenter at 4.3 km depth, with an uncertainty of 15 km but we assume that the earthquake lies within the aftershock zone, extending over 2.5-8 km depth. The USGS/ SLU regional moment tensor places the centroid depth at 3 km.
To add an additional constraint on the depth of the mainshock, we compute synthetic seismograms for two smaller earthquakes that were close to the mainshock and well recorded at T25A. We examine the seismograms of these earthquakes instead of the mainshock because the shorter source duration of these events results in simpler seismograms that more clearly show the effects of depth on wave propagation. Computing synthetic seismograms for the M L 2.9 foreshock discussed in the Foreshocks section and the M 4.0 aftershock at 14:11 UTC on 23 August 2011 confirms the above estimate of 3 km, whereby the surface waves and body-wave arrival time in the synthetics best match the observed waveforms at T25A for depths of 3 and 4 km (Fig. 10). This places both earthquakes in the uppermost crystalline basement.
Analysis of Interferometric Synthetic Aperture Radar (InSAR) image pairs that span the earthquake confirms these shallow hypocentral depths (W. Barnhart, personal comm., 2013). W. Barnhart identified deformation associated with the earthquake sequence that indicates a normal-faulting earthquake occurred with a similar strike and dip as seen in the seismic moment tensor. Results of analyzing the In-SAR images indicate that slip is concentrated within a zone of approximately 2.5-6 km in depth.

Postmainshock Seismicity: 23 August-15 December 2011
We analyzed the seismicity in the Raton Basin following the M 5.3 earthquake primarily using the temporary seismic network that was deployed by the USGS Geologic Hazards Team in Golden, Colorado. This network was fully operational by 25 August 2011. We examine the data through 15 December 2011, because we are most interested in the aftershocks of the M 5.3 mainshock, which had largely abated by then. The surface waves match best with the synthetics for earthquakes having depths of 3 and 4 km. Surface waves are the waves that are most diagnostic of earthquake depth, so this match is more important that a good bodywave match. Synthetics computed for depths greater than 5 km do not generate strong surface waves, which are evident in both earthquakes. Thus, we infer the earthquakes to lie between 3 and 4 km depth. Additionally, we find that the best P-and S-arrival time match is for depths of 3 and 4 km earthquake depth. Similar results are found for other early aftershocks and foreshocks. The color version of this figure is available only in the electronic edition.
Aftershocks. From the time of the mainshock through 2013, there have been 10 M ≥ 3:0 earthquakes in the vicinity. Seven of these events occurred within the first 48 hr, including the largest aftershock, an M 4.0 earthquake that occurred 8.5 hr after the mainshock. The three other M ≥ 3 earthquakes occurred in October 2011, December 2011, and May 2012. These three later events all occurred at least 4 km to the northeast of the mainshock hypocenter.
Examining the distribution of the seismicity, we find two intersecting lineations (Fig. 7b). The southern lineation strikes nearly north-south, and the northern lineation strikes approximately northeast-southwest. We infer that these semiplanar zones of seismicity represent faults. Although we were unable to find any discussion of these faults in the scientific literature, the strong alignment of the seismicity indicates that they are faults. The faults cover a region approximately 14 km in length. Although these lineations describe the general trend of seismicity, precisely located hypocenters from the temporary network do not define a simple plane for either trend. Rather, the lineations appear to be composed of en echelon subfaults 2-3 km in length that strike anywhere from slightly west of north to northeast (Fig. 7b). This contrasts with the seismicity in the 2001 earthquake sequence, which forms a single 7 km long, steeply dipping plane. It is also possible that the segmentation of the earthquakes may arise from an insufficient number of instruments that recorded these earthquakes, such that earthquakes appear to be preferentially moved toward the station closest to the earthquake instead of their true location.
The majority of the seismicity in the 2011 earthquake sequence lies at depths between 4 and 8 km below ground level, with the deepest seismicity appearing at the southern end of the aftershock zone. This increase in depth may be an artifact of bias in the location process, as there are no stations close to these events.
A composite focal mechanism (Snoke et al., 1984) from the aftershock sequence indicates normal faulting on a northsouth-trending fault plane (Fig. 8). This is consistent with the focal mechanisms for the larger earthquakes in the 2001 and 2011 sequences, the tectonic regime, and the lineations of seismicity in the area. As with the 2001 sequence, the eastdipping nodal plane implied by the composite mechanism appears to have a shallower dip than the plane defined by the seismicity.
Seismicity Outside the Aftershock Zone. In addition to the aftershocks of the M 5.3 mainshock, two smaller sequences occurred elsewhere in the Raton Basin between 25 August 2011 and 25 October 2011. One is about 20 km to the west of the August 2011 earthquake sequence, and the other is approximately 20 km to the southwest (Fig. 7a). The sequence to the southwest is notable in that it includes an M 3.9 earthquake that occurred 16 September 2011. We believe that neither of these additional earthquake sequences are related to the seismicity near the M 5.3 mainshock, given that its aftershock sequence diminished greatly after three weeks and because these earthquake sequences are ∼20 km from the obvious aftershocks of the 23 August mainshock.

Evidence that the 2011 Earthquake Sequence Was Induced by Wastewater Injection
The August-September 2011 earthquake sequence lies within 10 km of five injection wells in the Raton Basin: VPRC 14, VPRC 39, Hill Ranch, PCW, and Wild Boar (Fig. 7a). With the exception of Hill Ranch, these are highinjection-rate, high-volume wells. Cumulative injection volumes ranged between 15.1 million and 22.5 million barrels of wastewater through the end of August 2011 for all but Hill Ranch. Monthly injection rates at the same wells between January and August 2011 range from 220,000 to 262,000, 177,000 to 217,000, 78,000 to 91,000, and 88,000 to 124,000 barrels per month at VPRC 14, VPRC 39, PCW, and Wild Boar, respectively. The low-volume, low-rate Hill Ranch well injected 5000-46,000 barrels per month in the previous eight months.
Because the two VPRC wells are nearly colocated, we treat them as a single well for our analysis. At the time of the earthquake sequence, these wells were individually injecting at higher rates than the other wells in the area, and when summed they were injecting at far higher rates than any of the other wells. In the eight months preceding the earthquake sequence, the summed injection rate of the VPRC wells ranged between 408,000 and 479,000 barrels/month, whereas the well with the next highest injection rate in the same time period was Wild Boar, which injected 124; 000 barrels=month at most. Since 2005, the VPRC wells have been operating under gravity feed.
The proximity of the VPRC wells to the 2011 earthquake sequence also suggests that they are the wells most likely to have induced the earthquake sequence. The VPRC wells are the closest wells to the mainshock (3.7 km) and foreshock sequence (2.6 km from the first identified foreshock). The other wells (PCW, Wild Boar, and Hill Ranch) lie significantly further from the mainshock and foreshocks than the VPRC wells (9.9, 10.6, and 4.9 km from the mainshock, respectively). It is worth noting that Wild Boar and PCW lie closer than the other wells to the northernmost aftershocks of the earthquake sequence, so they may have played a role in inducing these earthquakes. The northernmost events could also be typical aftershocks, which commonly extend well beyond the initial mainshock rupture zone (Mendoza and Hartzell, 1988). Although we prefer the model whereby the VPRC wells are primarily responsible for this earthquake sequence, the cumulative effect of all or a subset of the nearby wells may have induced these earthquakes.

Summary
The Raton Basin of southern Colorado and northern New Mexico has seen a dramatic rise in the seismicity rates, beginning in 2001. In the past 13 years, there have been 16 M ≥ 3:8 earthquakes. In the previous 30 years only one M ≥ 3:8 earthquake occurred. The seismicity in the last 12 years includes M 5.3 and 5.0 earthquakes; the former caused damage to some unreinforced masonry buildings in the Colorado towns of Segundo and Valdez and minor damage in Cokedale and Trinidad (Morgan and Morgan, 2011). Based on the three broad lines of evidence below, we argue that the majority of these earthquakes have been induced by wastewater injection activities in the area.
1. We observe a large increase in earthquake rate in the Raton Basin shortly after major wastewater injection began in the area. This rate increase is limited to the area of industrial activity, with no significant change to earthquake behaviors outside the basin. Statistical analysis of the seismicity shows that the activity from 1970 to 2013 is unlikely to result from a constant-rate Poisson process.
Comparing the rate of earthquakes before and after the initiation of wastewater injection into the basin, there is a 3.0% probability that the observed activity could be produced by random variations in the background seismicity rate. 2. The vast majority of the seismicity lies within 5 km of wastewater injection wells. Careful investigation of two locally recorded earthquake sequences places them in close proximity to high-volume, high-injection-rate wells. The 2001 earthquake sequence epicenters surround the Wild Boar injection well, and the precisely located aftershocks are primarily less than 3 km below the injection interval. The 2011 earthquake sequence initiated within 2.6 km of two high-volume, high-injection-rate wells: VPRC 14 and VPRC 39. Combined, these wells were injecting more than 400,000 barrels of wastewater per month in the 16 months of injection prior to the M 5.3 mainshock. The injection intervals of these wells lie within 2 km depth of the shallowest hypocenters of the 2011 earthquake sequence. 3. The wells in the area are high volume and high injection rate. Individual wells greatly exceed the injection rates and cumulative volumes at the classic case of injectioninduced seismicity at the Rocky Mountain Arsenal. On a field-wide basis, injection in the Raton Basin is orders of magnitude larger still. Across the field, wastewater injection rates and earthquake rates show similar time histories. In particular, the seismicity starts shortly after injection rates increased significantly. Within the Colorado portion of the basin, both injection rates and rates of earthquake occurrence have remained fairly consistent since 2001. Beginning in 2006, when injection data becomes available, injection rates and earthquake rates in the New Mexico side of the basin are also constant. We also note that the total injection volumes since 2006 on both sides of the border are similar, and the numbers of detected earthquakes that have occurred on both sides are similar. This suggests that the volume of injected fluid or injection rate has some control on the seismicity that occurs.
Although there are many lines of evidence showing that the seismicity in the Raton Basin has been induced by wastewater injection activities in the area, it is very difficult to say whether an individual earthquake was caused by injection because natural seismicity has also been recorded there. Unfortunately, earthquakes that occurred when there was no local network deployed cannot be located with much accuracy, so the probability of fully describing these is low. For future research, a longer-term study with dense network coverage on both sides of the border would be especially useful in understanding the inducing relationship between the earthquakes and fluid injection in the Raton Basin.
It is also difficult to disentangle whether injection rate or cumulative volume controls if and when earthquakes are induced in the area. We expect that both cumulative injection volume and injection rates affect the rate and maximum magnitude of induced earthquakes, as has been suggested previously (McGarr, 2014;McGarr and Rubinstein, 2014).

Data and Resources
Los Alamos Seismic Laboratory catalog data is found in The facilities of the Incorporated Research Institutions for Seismology Data Management System (IRIS-DMS), and specifically the IRIS-DMS, were used for access to waveform, metadata, or products required in this study. The IRIS-DMS is funded through the National Science Foundation, specifically the GEO Directorate through the Instrumentation and Facilities Program of the National Science Foundation under Cooperative Agreement EAR-1063471. Some activities are supported by the National Science Foundation EarthScope Program under Cooperative Agreement EAR-0733069. Data from the USArray Transportable Array were made freely available as part of the EarthScope USArray facility, operated by IRIS and supported by the National Science Foundation under Cooperative Agreements EAR-0323309, EAR-0323311, and EAR-0733069. Waveform data from the 2011/2012 USGS temporary deployment can be accessed at IRIS. Additional waveform data used in the location of the 2011 earthquake sequence came from the Advanced National Seismic System backbone network and the USGS Albuquerque Seismological Laboratory Network (IU). IU stations are the part of the Global Seismic Network that are installed, maintained, and operated by the USGS Albuquerque Seismological Laboratory (http://earthquake.usgs.gov/regional/ asl/ and http://www.liss.org; last accessed January 2014). The latest information on U.S. stations and data availability may be viewed at http://earthquake.usgs.gov/monitoring/anss/ (last accessed January 2014).
The Upon publication, catalog data for the earthquakes located in this paper and pick data for the 2001 temporary seismic deployment will be inserted into ComCat. USGS/Saint Louis University regional moment tensor data comes from http://www.eas.slu.edu/eqc/eqc_mt/MECH.NA/index.html (last accessed January 2014).
For the time periods outside the two temporary, dense seismic deployments in 2001(i.e., 1963-August 2001and November 2001-July 2011, we use pick data from the International Seismological Centre Online Bulletin, supplemented by readings provided by Jim Dewey (written comm., 2012). We then use VELEST to compute the earthquake locations. Bootstrap resampling methods give horizontal uncertainties of 15 km. Varying the velocity model and earthquake starting locations suggests that location uncertainties may be even larger.
For those earthquakes that occurred during the temporary deployment in 2001, we use pick data from these stations. Bootstrap resampling methods give relative horizontal uncertainties of 200 m. Varying starting locations and velocity models gives absolute uncertainties of 2 km.
For those events that occurred while the temporary network was deployed in 2011, we use VELEST to locate these events. We compute a 1000 realization bootstrap analysis (Efron, 1979) of the VELEST locations and find a 1σ uncertainty across the entire region, located 300 m horizontally and 250 m in depth. We explored using the double-difference earthquake relocation algorithm hypoDD (Waldhauser and Ellsworth, 2000) to refine our locations, but we found that it did not significantly reduce uncertainties in earthquake locations, nor did the locations reduce to more planar features, so we chose to use the VELEST locations as our final locations.
We use locations from the Comprehensive Catalog (Com-Cat) for earthquakes that occurred after the dense seismometer deployment in 2011 (i.e., events between 15 December 2011 and 31 December 2013). Location uncertainties are likely on the order of 10-15 km.

Polarization Location Method: August 2011 Mainshock and Foreshocks
We use a method based upon P-wave polarization and S-P time to determine the location of the 2011 mainshock and its foreshocks. We use this method because we only have one station (T25A) that is within 30 km of the earthquake sequence at the time it initiated. T25A is located approximately 30 km to the east-northeast of the earthquake sequence (Fig. 10).
In this method, we first rotate the seismograms to identify the direction of P-wave polarization to establish the direction of wave propagation. For the mainshock, the initial Pwave motion is rectilinear on the horizontal components at a back azimuth of S73°W, implying that the earthquake lies to the south-southwest of T25A. It should be noted that the orientation of the components at T25A are N1°E and S89°E, instead of north and east (G. Ekström, personal comm., 2013). We then determine where a ray with that back azimuth intersects a two-segment model of the fault defined by the well-located aftershocks. Under the assumption that these earthquakes lie within the area defined by the aftershocks, the intersection between the ray described by the P-wave polarization and the simple fault model is deemed to be the epicenter of the earthquake. Observed S-P times at T25A are consistent with this assumption, with most events lying within 3 km of the fault segments.
In this analysis, epicentral uncertainty is controlled by (1) uncertainties in the polarity analysis, (2) uncertainty in the orientation of the horizontal components of T25A, and (3) uncertainty in the location of the fault segments. We are also making the assumption that these earthquakes lie on the same fault planes as the later aftershocks that were recorded by the temporary network. As noted above, S-P times indicate that these earthquakes are close to these fault planes. A conservative estimate of the uncertainty in the polarization analysis is 3°. Given an average source-receiver distance of 30 km, this gives an epicentral uncertainty of 1:5 km. After correcting for the slight misorientation of the horizontal components of T25A, the uncertainty in their orientation was determined to be 1.2° (Ekström and Busby, 2008). This gives an epicentral uncertainty associated with the orientation of the sensor of 0:6 km. The maximum width of the aftershock zone is approximately 2.8 km; therefore, given that the faults are projected to the center of the seismicity, this gives an epicentral uncertainty of 1:4 km. It should also be noted that the uncertainty in the aftershock locations, and thus the location of the faults that are used to locate the earthquakes, is approximately 300 m. In a worst-case scenario, in which all of the possible errors are additive, epicentral errors for the events located using the polarization method should not exceed 3:8 km. We expect that the errors are much smaller than this.
To assess the polarization method, we compare its locations with those from VELEST. We test it on the six earthquakes in the USGS Preliminary Determination of Epicenters (PDE) catalog that are M ≥ 2:5 during the period of the temporary deployment for which we computed VELEST locations (25 August 2011-15 December 2011). The VELEST locations have uncertainties on the order of 300 m in epicentral location. For the earthquakes examined, we find the mean distance between the VELEST locations and the polarization locations is 0.9 km, and the maximum is 1.4 km. Thus, a conservative estimate of the epicentral uncertainty of the polarization locations would be 2 km. This is likely a more accurate estimate of the uncertainty in the P-wave polarization locations than just summing all the possible sources of error.
Early Aftershocks of the 2011 Earthquake Sequence and Earthquakes in 2011 Prior to the August 2011 Earthquake Sequence We use a hypocentral decomposition location method (Jordan and Sverdrup, 1981;Hayes et al., 2013;McNamara et al., 2014) for two populations of earthquakes: (1) earthquakes that occurred in 2011 prior to the foreshock sequence that began 21 August 2011 and (2) early aftershocks that were recorded prior to the installation of the temporary seismic network. The uncertainty in these locations is 2:3 km, determined by the mean difference between these locations and the polarization locations used for the foreshocks and mainshock. Although the polarization method that we used to locate the foreshocks has lower uncertainties than the hypocentral decomposition event method, we do not use it because it became too time intensive, given the dramatic increase in earthquake rate following the mainshock.

Magnitude of Completeness
To statistically analyze the seismicity in the Raton Basin, it is important to be certain that the earthquake catalog is complete. Effectively, we need to compute the magnitude threshold for which all earthquakes are found within the earthquake catalog. We use the entire-magnitude-range (EMR) method (Woessner and Wiemer, 2005) and the maximum curvature method (Wiemer and Wyss, 2000) to compute the magnitude of completeness for our earthquake catalogs. We utilized the ZMAP software package to compute these values (Wiemer, 2001). For each computation, we treat the larger minimum magnitude from the EMR and maximum curvature methods as the magnitude of completeness.
The earthquake catalog that we use is a combination of the PDE, the U.S. National Hazard Map (Petersen et al., 2008), ComCat, and Los Alamos Scientific Laboratory seismic array catalogs. We also add in supplementary data from the temporary seismometer deployments in 2001 and 2011. Because we are concerned with the variation of the magnitude of completeness with time (i.e., if the magnitude of completeness is higher at earlier points in the catalog), we divide the catalog into 10-year increments and compute magnitudes of completeness. We compute this for each 10-year period beginning in 1970, with the last period extending through the end of 2013. For example, for a catalog beginning in 1970, we would examine the periods 1970-1979, 1980-1989, 1990-1999, and 2000-2013. We also increment the time periods by years, such that we examine every possible starting year, for example, 1971-1980, 1972-1981, 1973-1982, etc. We examine these time periods using catalogs with radii of 100, 200, 300, and 500 km surrounding a point in the center of the Raton Basin. We use these large radii to ensure that there is enough seismicity to compute a magnitude of completeness for the early years of the catalog when there was little seismicity in the region. Ⓔ Gutenberg-Richter plots are shown for each of the catalog combinations in Figures S1-S10 (available in the electronic supplement to this article).
For the different decade and distance-defined catalogs, we find that the magnitude of completeness varies from 2.7 to 3.8. Thus, a conservative approach for examining the seismicity would be to only consider earthquakes of M 3.8 and larger in our statistics. We also examine the catalog com-pleteness for 2000-present. We find that the magnitude of completeness for this most recent period is 3.0.