The Raton Basin had the highest number of earthquakes in Colorado and New Mexico from 2008 to 2010. The rate of both wastewater injection and earthquakes in the basin increased dramatically starting in 1999 and 2001, respectively. We compare seismicity ( M L 0.0–4.3) in the Raton Basin from 2008 to 2010 with the location of modeled pore pressure increases, estimated from cumulative wastewater injection volume beginning at the onset of well injection to present for all 28 injection wells in the basin. We find that modeled pore pressures in the New Mexico portion of the basin (above 0.08 MPa) reached that necessary to induce seismicity (0.01–0.1 MPa). We simulate a fault plane, 20 km long, inferred from seismicity in Vermejo Park (1355 of 1881 total earthquakes), in our model. We find that the relatively permeable fault allows pressures to migrate deeper into the basin at the onset of our study in 2008, providing an explanation for the observed seismicity in the basement. The Tercio lineament of earthquakes is similar to Vermejo Park fault in strike, but has fewer earthquakes (129) and is shorter in length (9 km). Seismicity in Vermejo Park occurs continuously, but earthquakes occur episodically in the remainder of the basin. The number of earthquakes we observe in seven seismic clusters increases as the cumulative injected volume from wells within 5 km increases. The modeled pore pressures, earthquake locations, and relationship between cumulative volume and number of earthquakes indicate that seismicity in the Raton Basin is likely induced.

1 Introduction The Raton Basin along the Colorado‐New Mexico border hosts coal‐bed methane deposits, extraction of which has taken place in the basin since 1999 (Hoffman & Brister, 2003). Underground Injection Control (UIC) class II wells (also referred to as wastewater disposal wells) were permitted beginning in 1988 to dispose of formation wastewater (Colorado Oil and Gas Conservation Commission, 2017; New Mexico Oil and Gas Commission (NMOCD), 2016), although the injection start date of the wells varies from 1994 to 2006, depending on the well. There are currently 28 wastewater disposal wells in the Raton Basin (Figure 1). Historical reports of earthquakes indicate that the Raton Basin region has been seismically active since at least 1966. The earthquake activity increased dramatically starting in the year 2001 making this region an area of interest for possible induced seismicity. Figure 1 Open in figure viewer PowerPoint 2014 2007 2014 2014 M W 5.3 earthquake. Black stars are the earthquakes after 1973 discussed in section M Wc 5.0 earthquake in Vermejo Park (Herrmann, 2016 Study area and seismicity of the Raton Basin. Seismicity from Rubinstein et al. () relocated catalog (1963–2014; green) and our catalog from 2008 to 2010 (blue). Dikes and sills (dark red) and faults (black lines) are shown. Axes of anticlines and synclines are dashed gray lines (Papadopulos and Associates, Inc.,). Black dashed line is a modeled fault surface projection from Barnhart et al. (). White diamonds are the saltwater disposal injection wells and are numbered. The Raton Basin is filled in with a tan color. High‐volume wells discussed in Rubinstein et al. () are purple‐filled diamonds. Red‐filled diamonds are the highest cumulative volume wells (summed as point source, since injection began in 1994) in the basin, 22 and 23. Red‐filled stars are the earthquakes over magnitude 4.0 from 2001 and 2011 sequences. The white star is the 20115.3 earthquake. Black stars are the earthquakes after 1973 discussed in section 1 . The 1966 and 1968 earthquakes are outside the boundaries of the map, but near the basin. Bright green star in New Mexico is the5.0 earthquake in Vermejo Park (Herrmann,). Earthquake clusters described in text are (a) Gulnare, (b) Spanish Peaks Foothills, (c) Boncarbo, (d) Weston, (e) Tercio, (f) Vermejo Park, and (g) Valle Vidal. Inset map shows the location of Raton Basin in Colorado and New Mexico. Section A‐A′ is shown in Figure 2 Because of geographic boundaries and limited seismic coverage, answering questions about the likelihood of induced seismicity has been difficult. The basin is nearly equally split between Colorado and New Mexico, which have different state geological surveys and different historical injection reporting requirements. Reporting of wastewater volumes was not recorded for New Mexico UIC wells from the beginning of injection (late‐1990s) to 2006. Previous seismicity studies have focused on the Colorado portion of the Raton Basin because it is the location of a 2001 earthquake swarm and the 2011 M W 5.3 earthquake, both leading to limited augmentation of seismic monitoring in that area. With the availability of the EarthScope Transportable Array (TA) stations to provide additional seismic station coverage throughout the basin, improved earthquake catalogs are now available (Nakai, Sheehan, & Bilek, 2017), allowing a more comprehensive view of seismicity in the entire Raton Basin in order to understand the impact of over a decade of wastewater injection in the basin. The Raton Basin is positioned to the east of the Rio Grande Rift, at the boundary of the Sangre de Cristo Mountains in the Great Plains physiographic province (Figure 1). In the basin, the Pierre shale lies atop the Dakota, Morrison, Entrada, and Glorieta Mesozoic sediments (Pillmore, 1969). These formations overlie Permian deposits and finally, Precambrian basement at varying depths from 1.2 to 2.8 km (Hemborg, 1996; Suleiman & Keller, 1985) (Figure 2). Coal deposits, mined in the last century and now exploited as coal‐bed methane, occur throughout the basin, but are highly altered in the southern basin from contact with widespread intruding sills (Bankey et al., 2002; Pehlivan, 2015; Pillmore, 1969). Figure 2 Open in figure viewer PowerPoint 1985 2015 1992 1992 Cross section of the Raton Basin (A‐A′ in Figure 1 ). Data are constrained by drill holes for formations above the Dakota and by a basement map of New Mexico (Suleiman & Keller,). Vertical exaggeration is 10.4:1. Names of formations: Tpc = Poison Canyon formation, TKr = Raton formation, Kv = Vermejo formation, Kt = Trinidad sandstone, Kp = Pierre shale, Kdp = Dakota sandstone and Purgatoire formation, Jmra = Morrison formation, TRig = Johnson Gap formation, PRsc = Sangre de Cristo formation, PCgr = Precambrian granite. Wastewater disposal wells inject primarily into the Dakota sandstone (Kdp) and the Morrison formation (above the Entrada and Glorieta formations) (Jmra). Basement contour (bottommost dashed line) is from an interpolated map in Weingarten (). Figure adapted from Stevens et al. (). Stevens et al. () do not differentiate between the Cretaceous Niobrara and Benton formations and the Pierre shale and Dakota sandstone, although the formations lie between the Pierre and Dakota. Wells 22, 23, 25, and 28 are projected approximately on the section. The Raton Basin first formed in the Pennsylvanian was deepened later in the Cretaceous during Laramide compression and is characterized by volcanic intrusions and Laramide faulting (Baltz, 1965). Coal‐bed methane, extracted from the Vermejo and Raton formations (Figure 2), is produced with hydraulic fracturing techniques (horizontal perforations in discontinuous coal seams) (Carlton, 2006; Hoffman & Brister, 2003; Environmental Protection Agency, 2015; Papadopulos et al., 2008). Intruding dikes, sills, Laramide faulting, and close proximity to the active Rio Grande Rift characterize the Raton Basin, whose geologic history is a unique combination of volcanism and Tertiary faulting not seen in Paradox Basin (Laramide faulting) and Greeley, CO (basement, but no surface, faulting), or basins in Oklahoma (basement faulting), three previously documented regions of induced seismicity (Keranen et al., 2014; King et al., 2014; Yeck et al., 2016). The western portion of the Raton Basin is marked by thrust faults dipping westward beneath the sedimentary cover, exposing the overturned Dakota formation (Bachman & Dane, 1962; Baltz & Bachman, 1956; Bolyard, 1959; Johnson, 1958; Northrop et al., 1946). The Dakota formation, the primary wastewater injection interval for disposal of wastewater associated with coal‐bed methane production (Colorado Oil and Gas Conservation Commission, 2017), is 9–18 m (30–60 feet) to 30–46 m (100–150 feet) in thickness from the southern to the northern Raton Basin and is up to 2.2 km below surface elevation (Figure 2). At several wells, the injection interval extends into the Entrada and Glorieta sandstones, intervals just below the Dakota. There is an interval of Pennsylvanian and Early Permian Sangre de Cristo sediments beneath the Glorieta (Baltz, 1965), separating the wastewater disposal interval from Precambrian basement. Volcanic intrusive activity in the Raton Basin occurred from ~30 to 22 Ma. The Spanish Peaks intrusive complex dominates the far northern basin, but dikes and sills abound across the basin. The radial dikes around the Spanish Peaks are generally oriented with the local state of stress created by the intrusions, while the east‐west 25 Ma dikes and the northwest‐southeast 22 Ma dikes reflect the regional stress orientation at the time of emplacement (Miggins, 2002). Currently, the state of stress is approximate east‐west extension as evidenced by 25 earthquake moment tensors (Herrmann, 2016). Laramide faulting is profuse on the eastern edge of the Sangre de Cristo Mountains, and features such as the Tres Valles fault in Colorado along the western basin and the Tercio anticline are evidence of thrust faulting at depth (Baltz, 1965; Stevens et al., 1992) (Figure 2). The Raton Basin, along with the Rocky Mountains and Rio Grande Rift, has high heat flow, with temperatures in the basin ranging from 200°C to more than 350°C (Blackwell et al., 2011). The Sangre de Cristo formation, suspected to have high permeability, is a target for geothermal circulation given the elevated geotherm in the formation (Bohlen, 2012, 2013). The Yellowstone Caldera (Wyoming) and the Valles Caldera (New Mexico) both have high heat flow (Blackwell et al., 2011). Yellowstone is seismically active while the Valles Caldera is nearly aseismic (Nakai et al., 2017; U.S. Geological Survey, 2017). The Raton Basin has been seismically active historically. Shaking from an m b 3.8–4.1 earthquake was felt near Cimarron, NM (south edge of the Raton Basin), and Weston, CO (within the Raton Basin) (Stover, Reagor, & Algermissen, 1983; von Hake & Cloud, 1968) on 24 September 1966. Felt reports are documented near Aguilar, Segundo, and Trinidad, CO (within the Raton Basin), from a M L 4.6 earthquake on 3 October 1966 (Stover, Reagor, & Algermissen, 1984). The earliest event in the U.S. Geological Survey (USGS) Comprehensive Earthquake Catalog (ComCat) is an M L 4.2 on 23 September 1973 near Valdez, CO (within the Raton Basin) (U.S. Geological Survey, 2017), and four other earthquakes were reported but not located instrumentally from 19 to 23 September 1973 (Coffmann et al., 1975; U.S. Geological Survey, 2017). The New Mexico Institute of Mining and Technology's seismic catalog from 1962 to 1998 contains four events >M d 1.3 in our study area (Figure 1) (Sanford et al., 2002). Catalogs developed by Sanford et al. (2006) and Pursley, Bilek, and Ruhl (2013) for 1999 to 2008 (prior to the start of this catalog) contain 78 earthquakes in the study area. The closest station used in the New Mexico Institute of Mining and Technology locations is SDCO, ~80 km northeast of the earthquakes, and there was a large azimuthal gap in station coverage to the north, east, and south. The location uncertainties in the Raton Basin are high for these catalogs (~10 km) (Sanford et al., 2006). There are 129 earthquakes >M L 3 in the USGS ComCat catalog from 1999 to 2016, and 15 earthquakes >M L 4, with eight in Colorado and seven in New Mexico (U.S. Geological Survey, 2017). A swarm of 12 earthquakes in August and September 2001, ranging in magnitude from 2.8 to 4.6, prompted the USGS to deploy a temporary seismometer array west of Trinidad, CO (Meremonte et al., 2002). That deployment recorded 39 earthquakes from 3 to 6 km depth in a 6 km long, northeast striking, steeply dipping (70–80°) fault plane. There were several fluid disposal wells within 5 km of the earthquake sequence. The 2001 earthquake sequence satisfied some but not all of the criteria popular at the time for distinguishing induced earthquake sequences (Davis & Frohlich, 1993). The sequence was not uniformly identified as induced. This was based on the prior history of earthquakes in the area, the 1 year time lag between the start of nearby well injection and the earthquakes, and the underpressured nature of the Dakota formation, which Meremonte et al. (2002) (incorrectly) suggested would preclude high pore pressure increases in the Dakota. On 23 August 2011, an M W 5.3 earthquake in the Raton Basin caused shaking and damage in nearby Trinidad, CO. The USGS deployed a temporary seismometer array from August to December 2011 to study the aftershocks of this earthquake. Rubinstein et al. (2014) recognized two fault strands, one striking north‐south and one striking northeast‐southwest, and concluded that the earthquakes occurred along a fault ~8 km deep, and were the result of wastewater injection in the VPRC 14/39 wells and the Wild Boar and PCW wells farther north. The existence of an 8 to 10 km long fault plane was corroborated by Barnhart et al. (2014) using interferometric synthetic aperture radar observations (Figure 1). The majority of analysis has focused on the Colorado portion of the Raton Basin because it is the location of the 2001 earthquake swarm and the 2011 M W 5.3 earthquake, though the number of moderate magnitude (>magnitude 3) earthquakes have increased as well in the New Mexico portion of the basin since 2001 (Figure S1 in the supporting information). In this work, we compare the seismicity patterns observed during the USArray and CREST seismometer deployments (from 2008 to 2010) with hydrogeologic modeling to show that the pressure generated from the cumulative injection in New Mexico could far exceed the pressure necessary to induce earthquakes in the southern basin, and a permeable fault we simulated in the model allows pressures to migrate deeper into the basement. We explore the relationship between cumulative injection volume and the number of earthquakes for the observed time period, and note differences in seismicity in Vermejo Park in comparison to the earthquake clusters in Tercio and the remainder of the basin.

2 Methods 2.1 Location, 1‐D Inversion, and Relocation of Earthquakes The earthquake location process is described by Nakai et al. (2017), who used the Antelope software package (Pavlis et al., 2004) to detect, associate, determine initial single‐event locations, and calculate local magnitudes. A total of 1881 earthquakes in the Raton Basin were located manually with P and S waves in Nakai et al. (2017). We also consider the earthquakes with no calculated magnitude in this work (135 earthquakes), which are part of the same catalog but unpublished. The lack of magnitude is due to difficulty measuring amplitude of low signal‐to‐noise ratio waveform recordings at individual stations. The earthquakes were initially located with a global standard International Association of Seismology and Physics 1991 (Kennett, 1991) velocity model. We used a subset of the 249 largest (>M L 1.8) earthquakes from 36.42° to 37.8° latitude and −105.5° to −104.1° longitude (the Raton Basin) (Figure 1) to invert for a 1‐D velocity model using VELEST, a simultaneous velocity model, earthquake location, and station correction least squares inversion (Kissling, Kradolfer, & Maurer, 1995). VELEST requires an a priori velocity model to calculate initial travel times. We tested a range of models based on two initial models, the first local to the Raton Basin used by Pioneer Natural Resources (P. Friberg, personal communications, 2015) and the second based on a western U.S. velocity model (Herrmann, 2016). Models began with five to eight layers that we eventually collapsed to four layers because we could not resolve near‐surface velocity layers close (~20 km) to the source. We used only arrivals with an epicentral distance of less than 170 km in the inversion to avoid fitting scatter in the travel time data due to P g to P n crossovers. We constructed Wadati diagrams (Wadati, 1933) to determine the V p /V s ratio for the starting S wave velocity models, but VELEST independently inverts for the S wave velocity model. Two of our final velocity models were very similar, the “slow” and “medium” models. The slow model had a slower upper crust and a faster lower crust in both P and S wave velocity in comparison to the medium model. The medium model was the model used for the final locations. The “fast” model was fast throughout the lower and upper crust for P and S wave velocities (all models given in Table S1 in the supporting information). The VELEST average RMS residual and the data variance for these models were similar (~0.14 s RMS residual), and the final model was chosen based on visual fit of the data to the reduced predicted travel time curves. We fixed the upper mantle P wave speed to 8.00 km/s as observed on the reduced travel time curves and S wave upper mantle speeds to 4.3 km/s for S wave velocities (Shen, Ritzwoller, & Schulte‐Pelkum, 2013). Earthquakes were relocated with Bayesloc (Myers, Johannesson, & Hanley, 2007, 2009), and only arrivals within 1.7° were used in the relocations to avoid the P n phase. The epicentral errors for the catalog are ±2 km for the Colorado and New Mexico catalog based on mine blast relocations (Nakai et al., 2017). The shallow structure (upper 5 km) is not well resolved due to distant station spacing, and we do not have enough arrivals at the close distance range needed to resolve shallow layers. These shallow structures are important for resolving the depths of near‐surface earthquakes, and we adopt the Bayesloc depth errors with the understanding that there is an inherent bias error in earthquake depths. Bayesloc depths are with respect to the average elevation in the Raton Basin. Relocated epicenters were nearly identical for the three models determined from the 1‐D inversion, although the depths varied. Regardless of which of the three models we used, the majority of earthquake depths were between 4 and 6 km, at basement depths. This agrees with the depths Rubinstein et al. (2014) determined in the Colorado part of the Raton Basin. Depths for four larger earthquakes were within 2 to 5 km of the regional moment tensors determined by Herrmann (2016). The depth of the regional moment tensors is chosen as the maximum value of best fit as a function of depth. The uncertainty varies with the depth sensitivity for waveform mechanisms. 2.2 Injection Well Data Monthly injection volumes by injection well were obtained online from the Colorado Oil and Gas Conservation Commission (2017) and NMOCD (2016). Data from October 1999 to May 2006 are not available for New Mexico wells, which led Rubinstein et al. (2014) to use produced water data from coal‐bed methane wells in New Mexico as a proxy for injection volumes. The produced water is a reasonable proxy for injection data because nearly 100% of water produced is injected in the area, and post‐2006, when production and injection data are both available, the two volumes correlate well, with the exception of 580,000 barrel in 1 month of 2009 (Rubinstein et al., 2014). We use the produced water data to estimate October 1999 to May 2006 injection volumes from six wells in New Mexico. The June 2006 to May 2008 average monthly volumes were used as a baseline of well injection. Well 25 was not operating until June 2008 and all wastewater disposal data are available online, so it was not included in the calculations. The average percent contribution from each well for 2006–2008 volumes for each well was calculated from NMOCD (2016) data. Each well began injecting on different months and years, so if only two wells had begun operating by May 2000, their relative percent contribution for May 2000 is estimated using their 2006–2008 relative percents. This allows us to reconstruct an approximate history of wastewater disposal volumes in the New Mexico Raton Basin. 2.3 Hydrogeologic Model 2000 (1) We developed a three‐dimensional hydrogeologic model of pore pressure diffusion from injection wells operating in the Raton Basin over the period from November 1994 to December 2010. Injection at well 8 (other wells began injecting subsequently) began in 1994, and the pore pressure in the basin is cumulative since the beginning of injection. Our seismicity catalog begins in 2008 and ends in 2010 (coincident with the time of USArray in Colorado and New Mexico), so we consider this 2 year time period in the pore pressure modeling. We simulate pressure diffusion in a layered, heterogeneous, and anisotropic medium using MODFLOW, a modular finite difference code developed at the USGS (Harbaugh et al.,). MODFLOW solves the groundwater flow equation in three dimensions for a fluid of constant density and dynamic viscosity in a heterogeneous and anisotropic aquifer with sources and sinks: where h is the hydraulic head (L); K xx , K yy , and K zz are the principal components of the hydraulic conductivity tensor (L/T); S s is the specific storage coefficient (L−1); and x, y, z, and t are the spatial and temporal coordinates. Q i is the fluid source or sink (T−1). Hydraulic conductivity is relatable to permeability by accounting for the specific weight and dynamic viscosity of water (Figure 3). Figure 3 Open in figure viewer PowerPoint Parameterization, boundary conditions, and vertical discretization of the base hydrogeologic model. Parameterization of the fault zone sensitivity analysis and plan view discretization are shown in Table S4 and Figure S6 1989 1988 2006 1974 2003 (2) k is the permeability (m2) and z is the depth (km). This relationship generally gives permeabilities that are about an order of magnitude lower for crystalline basement at similar depths than values from other notable k‐z relationships such as Manning and Ingebritsen ( 1999 The model inputs are based on the reported hydrostratigraphy of the sedimentary layering from Geldon () and use the reported ranges in bulk permeability from field studies and calibrated numerical models of similar strata in other basins (Table S4 ; Belitz & Bredehoeft,; Garavito, Kooi, & Neuzil,; Miller & Rahn,). Permeability of the intact crystalline basement relies on the permeability‐depth relationship derived by Shmonov et al. ():whereis the permeability (m) andis the depth (km). This relationship generally gives permeabilities that are about an order of magnitude lower for crystalline basement at similar depths than values from other notablerelationships such as Manning and Ingebritsen (). The layered heterogeneity of the sedimentary system is reflected in the modeled permeabilities with depth. A description of the formations above and below the injection horizon as well as heterogeneous layering included in the model is detailed in Figure S2 and Tables S4 and S5 (Bredehoeft, Neuzil, & Milly, 1983; Harbaugh & Davie, 1964; Jensen et al., 1954; Topper, Scott, & Watterson, 2011; Watts, 2006). The pressures calculated in the model runs are well below the reported underpressures in the basin, meaning that wastewater disposal can occur without the need for wellhead pressure (Rubinstein & Mahani, 2015). In this study, we adapt the numerical model developed by Weingarten et al. (2015) to simulate fluid pressure changes in the Raton Basin that includes the Vermejo Park fault zone. The adapted model specifically considers vertical flow in a complex of fault damage zones (Vermejo Park fault zone) nested inside an otherwise lower permeability sediment and crystalline basement (Figure 3 and Figures S2 and S3). Figure 3 shows the permeability structure and boundary conditions of the hydrostatigraphic units for the base case model (Table S4). While damage zones are highly anisotropic, several case studies of flow along fault damage zones have quantified the typical range of bulk diffusivities, the ratio between hydraulic conductivity and specific storage, to be 0.01 to 0.5 m2/s (Goebel et al., 2017; Xue et al., 2013). We simulate a range of permeabilities that correspond to bulk diffusivities in this range for the Vermejo Park fault zone and to assess the sensitivity of simulated pressure changes at depth to fault zone permeability (Table S5). We include a sensitivity run in which the Vermejo Park fault zone is not a conduit to fluid pressure with permeabilities reflecting the intact crystalline basement in the model.

3 Results Magnitudes for earthquakes in the basin in our catalog range from M L 0.0 to 4.3, and we find a b value of 0.90 using the maximum curvature method (Wiemer & Wyss, 2000). The magnitude completeness is M L 1.3 in the Raton Basin, as in the Nakai et al. (2017) catalog. Estimated b values for induced earthquakes span a wide range, from 0.6 to 2.1 (van der Elst et al., 2016). Earthquakes in the Raton Basin account for 29.4% of the cumulative seismic moment in the USGS ComCat catalog from 1973 to 2016 in Colorado and New Mexico. In the context of the regional earthquake catalog in Colorado and New Mexico from February 2008 to February 2010, the Raton Basin was by far the most seismically active area (Figure S4). Nearly 62% (1881/2764) of the earthquakes in the Rio Grande Rift catalog (Nakai et al., 2017) locate within the Raton Basin. Within the Raton Basin itself, seismicity varies both spatially and temporally. The extent of concentrated earthquake activity in the basin is in an area of ~80 km × 40 km with wells interspersed ~10 km apart. The catalog in the Raton Basin consists of 1881 earthquakes from February 2008 to February 2010 (Figure 1). The majority of these events are located inside the Raton Basin, although there are isolated earthquakes outside the geological basin that are included in the maps for local context of seismicity. The emphasis on earthquake activity thus far in the Raton Basin has been on the Colorado section of the Raton Basin (Barnhart et al., 2014; Rubinstein et al., 2014). Figure S1 shows that the New Mexico part of the Raton Basin has a similar number of M W 3.0 and greater earthquakes. The 2001 earthquake swarm (up to magnitude 4.6) and the 2011 M W 5.3 both occurred in the Colorado part of the Raton Basin. Moment release within the Raton Basin in New Mexico from 1973 to 2016 is 28% of the total compared to 72% in Colorado. When comparing estimated cumulative volume versus number of earthquakes (Figure 4; data in Figure S5 and Table S2), we found that cumulative injected volume is correlated with the number of earthquakes in nearby clusters. Well volumes are the sum since the time the wells began injection, which varies for individual wells, to February 2010. The cumulative volume is a sum of the well volumes from wells within 5 km (from the closest point, not the centroid) of each defined cluster, although we also considered distances of 10 km and 15 km. The clusters are chosen as continuous groups of earthquakes which have a well within 5 km. Figure 4 shows cumulative volume injected versus number of earthquakes on a log‐log plot. For wells within 5 km of our identified earthquake clusters, the cumulative volume correlates with earthquake frequency (Pearson's correlation coefficient of 0.77) (Figure 4). Considering wells within 10 km and 15 km of clusters, the correlation coefficients decrease to 0.68 and 0.72, respectively. As the cumulative volume from nearby wells increases, the number of earthquakes observed increases. We show the number of earthquakes and cumulative injection volume in Figures S6a and S6b. Figure 4 Open in figure viewer PowerPoint y = 214466x0.89. Wells within 5 km of the groups are included in the cumulative volume calculation. Cumulative volume of well injection and number of earthquakes is shown in Figure Cumulative injection volumes versus number of earthquakes for seven geographic groupings of earthquakes on a log‐log scale. Curve is= 214466. Wells within 5 km of the groups are included in the cumulative volume calculation. Cumulative volume of well injection and number of earthquakes is shown in Figure S5 The Vermejo Park fault is continuously seismically active throughout the 2008 to 2010 observation period, and the Tercio fault and the remaining clusters have periods of quiescence in time (Figure 5). In the 2008–2010 time period, we observe only one small earthquake within 5 km of the epicenter of the 2011 M W 5.3 earthquake (Figure 1). One ML 1.33 earthquake was located 3 km southwest of Well 9 on the NNE‐SSW fault identified by Rubinstein et al. (2014) and Barnhart et al. (2014), and 9 km north of the fault we find that a cluster of earthquakes occurs near Boncarbo in 2008–2010 (Figure 1). Our catalog (Nakai et al., 2017) is not a result of a network deployed specifically to study aftershocks of a large earthquake (e.g., 2001 and 2011). It seems unlikely that the earthquakes in this study are a long‐lived aftershock sequence, given that the number of magnitude 3 and above earthquakes returned to background levels the year following the 2011 M W 5.3 earthquake in Colorado, and the next most recent earthquake >M W 5.0 was a M Wc 5.0 in 2005. The Vermejo Park seismicity may very well be a continuation of a swarm that began before 2008; in fact, several earthquakes have been observed in the area beforehand. The pore pressure has potentially been high enough (Figure S8) such that these events may have been induced. Figure 5 Open in figure viewer PowerPoint (a) Cumulative number of earthquakes as a function of time for Vermejo Park, Tercio, and other clusters. Vermejo Park (blue) is continuously active throughout 2008–2010, and the remainder of the earthquake clusters and Tercio are more episodic (long stretches of time with no or very few earthquakes). (b) Ten‐day moving average of derivative of cumulative number of earthquakes per day in Vermejo Park (black), Tercio (blue), and the remaining clusters (red). There are five additional clusters of earthquakes that we observed in the 2008 to 2010 time period: from the north, Gulnare (91 earthquakes), Spanish Peaks Foothills (18), Boncarbo (60), Weston (37), and Valle Vidal (57) (Figure 1). There are numerous east west striking dikes near Gulnare, Boncarbo, Weston, and Spanish Peaks Foothills, and they are bounded to the north by the west and east Spanish Peak intrusive mountains. An 11 km long southeast striking mafic sill strikes divide well 26 from the Vermejo Park fault. On the southwest side of this sill, at a northwest‐southeast strike, are the Valle Vidal earthquakes. All of these clusters had wells within 5 km (Table S2). We took the derivative of the cumulative number of earthquakes and plotted a 10 day average derivative for three earthquake groups (Figure 5b). In 668 days, the average derivative for Vermejo Park is 2.0 (median of 1.5); for the grouped clusters of earthquakes the average derivative is 0.4 (median of 0) and 0.2 for Tercio (median of 0). A derivative of 0.5 or less means that 0 to 1 earthquakes occurred the day of or after a given day. The derivative is clearly much higher for Vermejo Park seismicity than Tercio and the rest of the earthquake clusters combined together (Figure 5b). We examine seismicity patterns in the two areas with the highest number of earthquakes in the studied period (2008–2010): Vermejo Park and Tercio (Figure 6). The majority of the 2008–2010 events (1355) are located in a distinct lineament (that we infer is a fault) surrounded by six wells in New Mexico. These events relocate into a 20 km long lineation that strikes north‐south which we define as the Vermejo Park fault (Figures 6a and 6c). Except for an absence of earthquakes several days prior to a M W 4.0 on 29 July 2009, earthquakes in Vermejo Park occur continuously. Some earthquakes located in 2011 were in the vicinity of Vermejo Park (Rubinstein et al., 2014), but the Vermejo Park seismic lineament has not been previously recognized as a highly active area of small‐magnitude seismicity with a coherent linear pattern. We do not observe a spatial migration of these earthquakes through time or depth for the short time period from 2008 to 2010. The depths range from 1 to 15 km, measured from a surface elevation of ~2 km with the majority of earthquakes in the 4–5 km depth range, with uncertainty of 2–5 km (see section 2). The depth range is corroborated by seismic moment tensors (Herrmann, 2016), which place basin‐wide earthquakes from 1 to 9 km (majority at 3–4 km depth) and New Mexico earthquakes from 3 to 9 km (majority at 4 km depth). There is a seismic gap along strike of the fault 6 km from the north end (Figure 6a). Allowing for error in the location of the 2005 M Wc 5.0 earthquake, it is possible that the seismic gap is due to this event. The fault appears to be steeply dipping (Figure 6) but given hypocenter uncertainty should be interpreted with caution. Figure 6 Open in figure viewer PowerPoint 2015 2015 Seismicity at Vermejo Park and Tercio (groups e and f in Figure 1 ). (a) Vermejo Park fault. The lineament contains 1355 earthquakes (circles, color coded by depth). Injection wells are shown with unfilled diamonds. Wells 22 and 23 make up more than half the projected cumulative volume in the entire basin. A southeast striking sill is dark red, while the basin syncline axis, coincident with the alignment of earthquakes, is a dashed black line. Earthquakes are colored by time from 1 May 2008. The Tercio anticline is shown by a dashed black line on the far left. Line with double‐sided arrows is the location of the section in Figure 6 c. (b) Tercio lineament in the western basin, near the Sangre De Cristo mountain front. The heavy black line is the edge of the basin, the gray lines are the mapped surface faults, and the dashed black lines are the Cuatro syncline (left) and the Tercio anticline (right). Dark red lines are mapped dikes/sills. Line with double‐sided arrows is the location of the section in Figure 6 d. (c) Vermejo Park fault depth cross section looking north. All earthquakes and moment tensors in Figure 6 a are projected onto Figure 6 c. Normal faulting is the dominant mode of faulting in the lineament. Earthquakes are colored by time (day) from 1 May 2008. Wells are unfilled diamonds above average elevation (2 km). Dashed black line is the interpolated basement depth at 36.9° from Weingarten () interpolated map. (d) Tercio fault depth cross section looking north. All earthquakes in Figure 6 b are projected onto Figure 6 d. Wells are unfilled diamonds above average elevation (2 km). Dashed black line is the interpolated basement depth at 36.9° from Weingarten () map. The feature that we refer to as the Tercio fault (Figures 6b and 6d), inferred from seismicity, also strikes north‐south, similar to the Vermejo Park fault, but is only 9 km in length with 129 earthquakes in the 2008 to 2010 period. The fault dip as estimated from the seismicity is near vertical and the earthquakes are shallower than for the Vermejo Park fault (subject to depth uncertainties). The southern half of seismicity is aligned with the Tercio anticline (Stevens et al., 1992), which is along the western edge of the Raton Basin, adjacent to the exposed Dakota formation and mapped surface thrust faults (Figures 1 and 6). We modeled pore pressures across the Raton Basin, incorporating the derived 1999–2006 injection volumes in New Mexico into the cumulative volume history. As described in section 2, we used the produced water data to estimate 1999–2006 injection volumes from six wells in New Mexico. All wells operating in the basin contribute to the pore pressure changes in the basin. By February 2010, the pressure perturbation from injection has spread throughout most of the dominant injection reservoir, the Dakota formation (Figure 7). We focus on the pore pressure changes in the New Mexico portion of the basin, specifically along the Vermejo Park fault zone above 0.08 MPa pressure (green contour in Figure 7), which is 8 times the threshold at which earthquakes may be triggered from injection (Hornbach et al., 2015; Keranen et al., 2014), includes most of the Vermejo Park earthquake cluster at the onset of this study in 2008 (Figure S8). When the Vermejo Park seismic lineament is modeled as a permeable fault, the model predicts pressure perturbations reaching depths of 12–13 km by February 2010 (Figure 7b). Pressure changes above 0.2 MPa can be observed at basement depths below ~2.5 km. The highest‐volume wells in the basin were operating in the region of the highest modeled pressures (red in Figure 7). We identify high cumulative volume wells over the time period 1994 to 2008 (the onset of this study) as those with greater than 15 million barrels of water injected. Those wells are 8 and 19 in Colorado and 22, 23, 27, and 28 in New Mexico (Figure 1). The wells 22 and 23 are collocated and as a point source, have injected more than double the water of any other well in the basin. Figure 7 Open in figure viewer PowerPoint (a) Plan view of pore pressure change in February 2010 at the depth of Dakota formation in the Raton Basin. Earthquake epicenters from February 2008 to February 2010 are plotted as white circles. The black north‐south line is a no‐flux boundary condition representing the westernmost extent of the Raton Basin. The brown line depicts the extent of the Raton formation at the surface, a demarcation of the Raton Basin proper. The black star represents the location of calculated pressure histories shown in Figure 7 c. (b) Cross section B‐B′ of pore pressure change running north‐south through the Raton Basin along strike of the Vermejo Park fault zone. The black star denotes the 4.5 km depth of calculated pressure histories shown in Figure 7 c. Solid black lines denote the hydrostratigraphic units described in Figure 3 and Tables S3 and S4 . (c) Pressure versus time for four modeled Vermejo Park fault zone permeabilities, including a run where no fault is included. The highlighted yellow region of the pressure history denotes the time frame of interest in the study from 2008 to 2010. Sensitivity analysis of the modeled pressure perturbation in the Vermejo Park fault zone shows that pressure changes through time are sensitive to fault zone permeability (Figure 7c). Varying the fault zone permeability from that of intact crystalline basement (i.e., no fault) to as high as 3 × 10−5 m2 yields pressure changes at 4.5 km depth from 0.004 to 0.1 MPa. All models show pressure changes increasing with time, including throughout the time period of interest of this study from 2008 to 2010. The pressure history from the lower bound of the modeled permeable fault zones (10−16 m2) shows that pressure changes at depth are significant and likely large enough to trigger seismicity at seismogenic depths.

4 Discussion The Pierre shale above the Dakota serves as a fluid flow barrier, which effectively prevents fluids from migrating upward; thus, modeled fluid pressure in the Raton Basin propagates primarily laterally in the Dakota, Entrada, and Sangre de Cristo formations and diffuses downward into basement rock (Figure 7b). We did not observe any seismicity in 2008–2010 near the M W 5.3 earthquake in 2011. The active seismicity in New Mexico has not been documented to the extent that the Colorado seismicity has, but the number of small‐magnitude earthquakes during the 2008 to 2010 time period is much higher in the New Mexico basin than in Colorado (Figure 1), and slightly higher in New Mexico for earthquakes >M W 3 (Figure S1). In terms of moment release during 2008 to 2010, 26% of the moment release occurred in Colorado while 74% occurred in New Mexico, while 72% of long‐term moment release (1973 to 2016) is in Colorado and 28% is in New Mexico (magnitude 2.5 and above). We note that hydraulic fracturing, associated with earthquakes in other areas (e.g., Holland, 2013; Schultz et al., 2017; Skoumal, Brudzinski, & Currie, 2015), does occur in the Raton Basin in the Raton and Vermejo formations (Figure 2) (Carlton, 2006). We have not assessed the impact hydraulic fracturing may have on the potential for earthquakes in the Raton Basin. Although there is strong evidence that earthquakes in the Raton Basin are induced (Rubinstein et al., 2014), no quantified causative mechanism has previously been shown to drive these earthquakes. Our pore pressure modeling shows that pore pressure in the basement rock can exceed the threshold at which earthquakes are induced (>0.01–0.1 MPa). The pore pressure for the majority of the basin populated by wells exceeds triggering thresholds observed in Oklahoma (0.07 MPa) (Keranen et al., 2014), the Ellenburger formation in north Texas (0.09 MPa) (Hornbach et al., 2016), and in Azle, TX (0.01 to 0.2 MPa) (Hornbach et al., 2015). The highest pore pressure in the basin is in the region where we observe the highest number of earthquakes in 2008–2010, and the high‐volume wells (22 and 23) are within 10 km of six earthquakes >M W 3. We simulated a permeable (relative to the surrounding crystalline basement) fault that allows pore pressures to migrate below the sedimentary disposal interval into the basement. The pattern of seismicity in Vermejo Park is consistent with high pore pressures near high‐volume wells and a basement fault. Figure 7b demonstrates the effect of adding a permeable fault to the model, which is the increase in pore pressure at seismogenic depths. This is one explanation for the spatial pattern of seismicity in Vermejo Park we observe: the modeled fault has allowed pressures to migrate far along the fault to the south of the high‐volume wells. Poroelastic stress effects may play a role in earthquake triggering in the Raton Basin, and a geomechanical model (Fan, Eichhubl, & Gale, 2016; Goebel et al., 2017) of the basin would be potentially beneficial. The regional horizontal least compressive stress orientation is east‐west in and around the Rio Grande Rift (Berglund et al., 2012), although the orientation of stress changes in eastern Colorado, according to recent earthquake moment tensors (Herrmann, 2016) and historical stress data (Zoback & Zoback, 1980). Lund Snee and Zoback (2016) inverted moment tensors in the Raton Basin and determined a maximum horizontal stress direction of 167.1°. Three moment tensors in the Vermejo Park seismic lineament are consistent in strike with our interpretation of a north‐south striking normal fault (Figure 6c). If the Vermejo Park fault is one continuous fault, the length is greater than the Pawnee fault in Oklahoma on which a M W 5.8 earthquake occurred (Yeck et al., 2017) and similar in length to the Fairview fault, also in Oklahoma, on which a M W 5.1 earthquake occurred (Yeck et al., 2016). In 2005, a M Wc 5.0 earthquake occurred near the north of the Vermejo Park fault, the largest earthquake in New Mexico Raton Basin in the USGS ComCat (U.S. Geological Survey, 2017). The strike and dip of the 2005 earthquake was 160° and 40° (or 5° and 53°) with a depth of 4 km. The moment tensor (Herrmann, 2016) for this earthquake located 2.5 km to the west of the USGS location, hinting at the epicentral uncertainty of the location, although the uncertainty in regional network locations is ~15 km (Rubinstein et al., 2014). The moment tensor location is coincident with the Vermejo Park lineament and it is possible that the 2005 earthquake took place on the fault. The dip of the earthquakes is difficult to reconcile with the seismicity we observe, which is may be taking place on en echelon splays off the larger fault and appears as a vertical cloud in section. The Tercio fault is similar, although shorter in length (9 km), less seismically active (129 earthquakes), and the seismicity is shallower (1–11 km). The lineation of seismicity and the presence of mapped surface faults 10–15 km west of Tercio are strong evidence that this is also a fault striking north‐south that has been activated by high pore pressures. Displacement of basin sediments along the Tres Valles fault ends north of the Tercio anticline, but likely displaces the Precambrian granite in the anticline (Stevens et al., 1992). A M W 3.5 earthquake occurred 4 km east of the Tercio fault (4 December 2012) with a moment tensor indicating east‐west normal faulting, while a moment tensor for an M W 3.53 earthquake on 29 September 2009 at the northeast terminus of the fault indicates slight southwest‐northeast directed normal faulting (Figure 6b). We note that our earthquake depths extend from shallow (1–2 km depending on the velocity model used) to 15–20 km. The shallow depths are not well determined, but the majority of earthquakes occur at 4–6 km in all three models. Although we do not have good constraints on depth due to station spacing, we are confident that many of the earthquakes lie in the basement. In the Paradox Basin of western Colorado, where continuous injection of salt water at a single well began in 1996, the vast majority of earthquakes lie in the sedimentary layers above the basement near the depth of injection (4.88 km), although earthquakes occur up to 6.5 km in depth (Block, 2010). In Fairview, OK, earthquakes occur up to 6.5 km below the sedimentary‐basement interface, and the majority of earthquakes occur from 6 to 9 km depth (Yeck et al., 2016). In Greeley, CO, the majority of the earthquakes occur in the basement, at less than 6 km depth (Yeck et al., 2016). McNamara et al. (2015) found aftershocks of the M W 5.6 Prague, OK, earthquake at depths of 1–10 km. Induced earthquakes in north Texas are <8 km deep, ranging from the top of the Ellenburger disposal formation (1 km to 2.8 km) into the basement (Hornbach et al., 2016). The depth extent of the Raton earthquakes is unusually deep, and although we focus on the observation that the majority occur at 4–6 km depth, this is an interesting topic for future study. These Vermejo Park and Tercio faults are likely suitably oriented normal faults in the east‐west stress regime if two of these seismic lineaments are observed. These high‐density linear patterns of seismicity with a scarcity of earthquakes surrounding the lineament are similar between Vermejo Park, Tercio, and the Fairview, OK, faults, on which a M W 5.1 occurred at 8–10 km depth (Yeck et al., 2016). This is similar to lineaments in Arkansas and Kansas (Choy et al., 2016; Horton, 2012). If we take a subsurface rupture length of 10–20 km, the Vermejo Park and Tercio faults could potentially produce a magnitude 6 to 6.5 earthquake (Wells & Coppersmith, 1994). As cumulative volume injected increases, the number of earthquakes tends to increase (Frohlich et al., 2011). For wells within 5 km of our identified earthquake clusters, the cumulative volume correlates with earthquake frequency (Pearson's correlation coefficient of 0.77) (Figure 4). Considering wells within 10 km and 15 km of clusters, the correlation coefficients decrease to 0.68 and 0.72. A similar relationship between cumulative volume and earthquake frequency has been reported in the Carthage Cotton Valley Gas Field, TX (Rutledge, Phillips, & Mayerhofer, 2004). The Raton Basin has a cumulative injected volume of 45 million barrels of water since 1994, compared to 1.7 billion barrels of water in the Ellenburger formation, TX, since 2006. We see two types of seismicity in the basin, continuous and episodic (Figure 5). We define continuous earthquakes as those with an average derivative (earthquakes/day) of 1.0 or higher (Figure 5b), although this may vary based on the magnitude completeness of a catalog and detection threshold of the seismic network. Continuous seismicity refers to relatively steady increase in number of earthquakes over the catalog period, whereas discrete pulses of seismic activity with gaps of time between earthquakes characterize the episodic seismicity (Figure 5). Vermejo Park seismicity is continuous and is near the locus of the highest changes in pore pressure in the basin (when a fault is simulated), whereas the episodic seismicity is in regions of lower changes in pore pressure (~0.08 MPa in Figure 7). Episodic seismicity has been observed in areas of induced seismicity, for example, in the 2008–2009 Dallas‐Fort Worth sequence and the 2009–2010 Cleburne, TX, sequence (Frohlich et al., 2011; Justinic et al., 2013). In northern Texas, induced seismicity does not follow a main shock‐aftershock pattern but is characterized by swarms of small earthquakes (Hornbach et al., 2016). Episodic seismicity from 1975 to 1977 with a magnitude detection threshold of 2.0 is also observed in the Permian Basin near the Texas‐New Mexico border and is spatially clustered (Rogers & Malkiel, 1979). Given the fact that a M Wc 5.0 has occurred in New Mexico, a M W 5.3 has occurred in Colorado, and the wells continue to dispose of wastewater, continually raising the basin pore pressure, we expect that the earthquake hazard is substantial in the Raton Basin.

5 Conclusions We show that elevated pore pressures can exist in the Raton Basin, especially in New Mexico, at sufficient levels (>0.03 MPa) to induce earthquakes as early as May 2005, provided that the Vermejo Park fault is simulated. The two most seismically active faults during this time period are the Vermejo Park and Tercio lineaments. We model the Vermejo Park fault as a fault plane in the hydrogeologic model, which results in pressures of ~0.08 MPa at 9 km depth (well below basement). In Vermejo Park, this is indicated by three north‐south striking normal faulting moment tensors within the cluster and north‐south seismicity alignment. Tercio is 10 to 15 km east of the faulted and exposed Sangre de Cristo formation and Precambrian granite of the Sangre de Cristo Mountains, is a 9 km north‐south alignment of seismicity at basement depths, with two normal faulting moment tensors within 5 km of the lineament. Although the depths are not well constrained due to station spacing, the majority are likely at basement depths, which is common in other cases of induced seismicity. The spatial patterns of seismicity we observe are reflected in distribution of wastewater injection and modeled pore pressure change when a fault is simulated. Modeled pore pressures at seismogenic depths are well above the earthquake‐triggering threshold under a variety of possible fault permeability scenarios. There is a power law relationship between cumulative wastewater disposal volume and number earthquakes. There is a similarity between of patterns of seismicity (basement earthquakes, seismic lineaments, and spatial clustering) in the Raton Basin and other regions with induced seismicity. The number of earthquakes of the observed clusters is directly correlated to the cumulative injection volume from wells within 5 km of the seismicity. As the cumulative volume of wastewater injection increases, the pore pressure nearby increases, increasing the likelihood that earthquakes will be triggered on existing faults. If this relationship holds, we would expect to see an overall increase of small‐magnitude seismicity as the pore pressure continually increases in the basin. We observe two types of seismicity in the Raton Basin, continuous and episodic. Vermejo Park is continuously active, and the remainder of the clusters and the Tercio fault are episodically active.Four main lines of evidence support the contention that seismicity in the Raton Basin is induced from long‐term wastewater injection:

Acknowledgments The authors would like to thank Justin Rubinstein for generously sharing the produced water data in New Mexico and the relocated USGS earthquake catalog in the Raton Basin from 1963 to 2011, Stephen Myers from LLNL for the updated Bayesloc location code, Craig Jones for the helpful comments and suggestions, and Paul Friberg and Pioneer Natural Resources for the Raton Basin starting velocity model. The authors would also like to thank Heather DeShon and an anonymous reviewer for thorough comments which have improved the manuscript. Many figures were created with GMT (Wessel et al., 2013). Thanks to Josh Stachnik and the PIs of the CREST experiment for supplying the CREST data. Funding came from an AGEP and NSF GRF to J.N., NSF EAR 1053596, and NSF EAR 1053597. MW was partially supported by the Stanford Center for Induced and Triggered Seismicity. Data from the TA network were made freely available as part of the EarthScope USArray facility, operated by Incorporated Research Institutions for Seismology (IRIS) and supported by the National Science Foundation, under cooperative agreement EAR‐1261681. The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to waveforms, related metadata, and/or derived products used in this study. IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope (SAGE) Proposal of the National Science Foundation under cooperative agreement EAR‐1261681.

Supporting Information Filename Description jgrb52328-sup-0001-AGU_Supporting_Information.docxWord 2007 document , 1.7 MB Supporting Information S1 Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.