TY - JOUR
AU - Keranen, Katie, M
AB - SUMMARY The orientations of faults activated relative to the local principal stress directions can provide insights into the role of pore pressure changes in induced earthquake sequences. Here, we examine the 2011 M 5.7 Prague earthquake sequence that was induced by nearby wastewater disposal. We estimate the local principal compressive stress direction near the rupture as inferred from shear wave splitting measurements at spatial resolutions as small as 750 m. We find that the dominant azimuth observed is parallel to previous estimates of the regional compressive stress with some secondary azimuths oriented subparallel to the strike of the major fault structures. From an extended catalogue, we map ten distinct fault segments activated during the sequence that exhibit a wide array of orientations. We assess whether the five near-vertical fault planes are optimally oriented to fail in the determined stress field. We find that only two of the fault planes, including the M   5.7 main shock fault, are optimally oriented. Both the M 4.8 foreshock and M   4.8 aftershock occur on fault planes that deviate 20–29° from the optimal orientation for slip. Our results confirm that induced event sequences can occur on faults not optimally oriented for failure in the local stress field. The results suggest elevated pore fluid pressures likely induced failure along several of the faults activated in the 2011 Prague sequence. Induced seismicity, Rheology and friction of fault zones, Seismic anisotropy INTRODUCTION On 6 November 2011 at 3:53:10 UTC a M  5.7 earthquake occurred near the town of Prague, Oklahoma, and was subsequently linked to nearby wastewater disposal (Keranen et al. 2013). The main shock was preceded ∼20 hr beforehand by a M  4.8 foreshock and the largest aftershock, a M 4.8, occurred ∼48 hr after the main shock. The Prague sequence was the first in a series of moderate earthquakes that occurred in central to northern Oklahoma related to wastewater disposal, including the M 5.0 Cushing, M  5.1 Fairview and M  5.8 Pawnee earthquakes (e.g. McNamara et al. 2015; Yeck et al. 2016; Chen et al. 2017; McGarr & Barbour 2017). More generally, an increase in the earthquake rates above background in Oklahoma and across the Central and Eastern United States was observed beginning in 2009 (Ellsworth 2013; Llenos & Michael 2013). Earthquake rates in Oklahoma peaked in 2014–2015 and have since declined (Norbeck & Rubinstein 2018) as regulations and economics have reduced disposal volumes. Much of the increase in earthquake rates in Oklahoma has been linked to deep disposal of saltwater resulting from oil or gas production (Keranen et al. 2014; Walsh & Zoback 2015; Weingarten et al. 2015), with a smaller portion of earthquakes attributed to hydraulic fracturing (Holland 2011; Skoumal et al. 2018). Saltwater disposal primarily targets the Arbuckle Group, a highly permeable sedimentary unit that is directly upsection from the granitic basement (Ham 1973; Morgan & Murray 2015). At its peak in 2014, the annual volume of saltwater disposal into the Arbuckle reached more than 1.05 billion barrels (Murray 2015). Pore pressure observations of tidal response in the Arbuckle indicate vertical fluid flow and have been interpreted to suggest the Arbuckle is not fully confined (Wang et al. 2018; Barbour et al. 2019). Faults and fractures that extend from the basement into sedimentary sections likely provide hydraulic pathways for vertical migration of fluids into the basement (Schwab et al. 2017; Williams 2017). Overall, the observed temporal changes in pore fluid pressures in the Arbuckle Group are complex and reflect both long-term pressure fluctuations resulting from disposal as well as short-term responses to nearby earthquakes (Kroll et al. 2017; Barbour et al. 2019). Saltwater disposal can induce nearby earthquakes by altering the stresses along faults primarily through increased pore fluid pressures (Ellsworth 2013; Norbeck & Rubinstein 2018), and to a lesser degree poroelastic stress changes (Segall & Lu 2015; Goebel & Brodsky 2018). An increase in the pore fluid pressure acts to reduce the effective normal stress along pre-existing faults and fractures in the Precambrian basement, so the Coulomb failure criterion may be exceeded at lower shear stresses allowing fault slip to occur (Handin et al. 1963). Further, increased fluid pressures may allow the failure envelope to be exceeded at lower differential stresses (Goertz-Allmann et al. 2011). Failure at lower differential stresses may be supported by observations of lower stress drops of small magnitude induced events (Boyd et al. 2017; Sumy et al. 2017; Trugman et al. 2017), although the cause of the lower stress drops and even whether the stress drops are in fact different from tectonic events is still an open question (Huang et al. 2017; Wu et al. 2018). The probability that a fault will slip as pore fluid pressures reduce the effective normal stress is higher if the fault is optimally oriented relative to the current, local tectonic stress field, for example a vertical fault striking ∼30° from the maximum horizontal principal stress direction in a strike-slip faulting regime. Indeed, multiple studies conclude that induced earthquakes occur preferentially on faults that are optimally oriented in the current stress field (Schoenball & Ellsworth 2017; Skoumal et al. 2019) and that stress state may control rupture propagation (Jones 1988). The supposition that certain fault orientations are more likely to fail is sometimes used to define the hazard (or lack of hazard) of induced seismicity near known faults depending on their orientations (Walsh & Zoback 2016; Alt & Zoback 2017; Snee & Zoback 2018). However, failure along unfavourably oriented faults (faults oriented more than ±15° from optimal orientation) can occur as increased pore fluid pressures reduce the frictional strength of pre-existing faults or fractures (Hubbert & Rubey 1959; Sibson 1985, 1990). For example, during the 1960s Denver sequence that was induced by local injection it was inferred that a range of pre-existing fracture orientations were activated as fluid pressures increased (Healy et al. 1968). Faults that are unfavourably oriented or misoriented may also be expected to slip aseismically when triggered by fluid pressure (Zoback et al. 2012). The Prague earthquake sequence was observed to rupture along three primary fault planes, defined by a M  4.8 foreshock, M  5.7 main shock and M  4.8 aftershock (Keranen et al. 2013). The focal mechanisms of the three largest events in the Prague sequence suggest variable fault strikes, which is also supported by published catalogues of the sequence (Keranen et al. 2013; Sumy et al. 2014; McNamara et al. 2015; McMahon et al. 2017). The regional stress field from inversion of focal mechanisms near the Prague earthquake show the maximum principal stress (    σ 1) is near N84E (Sumy et al. 2014; Walsh & Zoback 2016). Similarly, Alt & Zoback (2017) more recently estimated σ 1 as N84E ± 3 from focal mechanism inversion and wellbore breakouts. The foreshock fault plane strikes subparallel to previously mapped faults in the regional field (Keranen et al. 2013), and likely has low probability of slip (Walsh & Zoback 2016). The focal mechanism of the M5.7 main shock was found to be similar to the best fit regional focal mechanism (Sumy et al. 2014), and Walsh & Zoback (2016) estimated a relatively high probability of slip along this plane when considering the regional stress field. The M  4.8 aftershock occurred on a nearly east–west oriented plane but has not been evaluated to determine if it is well-oriented in the regional stress field. The relationship between the fault planes activated during the Prague sequence and the local stress field remains an open question. For hazard assessments, it is important to understand the interactions between local stress-state and the increases in pore fluid pressure for faults that are unfavourably oriented relative to the local stress field. While regional principal stress orientations have been resolved near the Prague sequence (Sumy et al. 2014; Walsh & Zoback 2016; Alt & Zoback 2017), it is not known if stresses vary more locally. It has previously been observed that stresses may rotate close to active fault structures using mapped fracture orientations (Rawnsley et al. 1992) and borehole measurements (Barton & Zoback 1994; Hickman & Zoback 2004). Spatial or temporal rotations in stress orientations have also been inferred local to faults from inversion of focal mechanisms (Provost & Houston 2001; Bohnhoff et al. 2006), although other focal mechanism inversion analyses show no such rotations even close to major faults (Provost & Houston 2003; Hardebeck & Michael 2004). In Oklahoma, Skoumal et al. (2019) identified a ∼20° rotation near the Nemaha fault. An estimate of the spatial distribution of the horizontal stress orientations around faults can also be obtained using observations of velocity anisotropy. Crustal velocity anisotropy typically results from preferential closure of microcracks normal to the maximum compressive stress, but may also reflect shear fabric aligned with fractures and faults or intrinsic anisotropy caused by aligned minerals (Aster et al. 1990; Kaneshima 1990; Zinke & Zoback 2000; Cochran et al. 2003; Peng & Ben-Zion 2004; Balfour et al. 2005; Boness & Zoback 2006; Holt et al. 2013). Shear velocity anisotropy is commonly measured using a technique known as shear wave splitting (SWS, Crampin 1985; Silver & Chan 1991; Savage et al. 2010). Shear waves travelling through an anisotropic medium become polarized parallel and perpendicular to the aligned microcracks. The shear wave propagating parallel to the strike of the cracks is essentially insensitive to the presence of the cracks (both micro- and macrocracks), while the wave propagating perpendicular to the strike of the cracks is slowed by their presence. In strike-slip faulting regimes, such as the region around the Prague sequence, σ1 is equivalent to the maximum horizontal stress, SHmax. Thus, orientation of the fast shear wave can provide an estimate of the orientation of σ1, and the delay time between the fast and slow arriving shear wave provides a measure of the density of cracks along the propagation path. With sufficient ray coverage, it is possible to identify local spatial variability in the fast direction orientations that may reflect local changes in stress orientations or the existence of a shear zone along a fault (Liu et al. 2008, 2015; Johnson et al. 2011; Cochran & Kroll 2015). Here, we use SWS to measure and infer the likely cause of the shear velocity anisotropy and to identify spatial variability of fast directions near the Prague sequence. The Prague sequence was one of the best-recorded aftershock sequences of a moderate magnitude induced earthquake, so we refine the imaging of the subfaults activated in the sequence by applying advanced event detection and fault fitting techniques. We apply template matching methods (e.g. Shelly et al. 2016; Ross et al. 2017; Skoumal et al. 2019) to extend the Prague catalogue beyond what has been previously published (Keranen et al. 2013; Sumy et al. 2014; McNamara et al. 2015; McMahon et al. 2017; Savage et al. 2017). Catalogues complete to lower magnitudes illuminate detailed fault structures and track the evolution of earthquake sequences (Shelly et al. 2016; Ross et al. 2017a,b,2019; Cochran et al. 2018; Skoumal et al. 2019). Using the extended catalogue, we can better map the subfaults activated in the sequence and their orientations. Within the context of the identified fault structures, we explore the spatiotemporal characteristics of the sequence evolution. Additionally, we examine whether the identified subfaults are optimally oriented within the local stress field determined from SWS. DATA AND METHODS A temporary seismic deployment around the 2011 Prague, Oklahoma, earthquake sequence was initiated following the M  4.8 foreshock, with a total of 31 seismometers installed in approximately one week (Keranen et al. 2013). The deployment recorded continuous seismic data for a period of 90 d (through 1 February 2012). The deployment includes ∼18 stations (network code: LC) within ∼10–15 km of the fault planes activated in the sequence (a subset of the stations is shown in Fig. 1). The supplemental information of Sumy et al. (2014) provides detailed information on the station locations, dataloggers, sensors, and data sample rates. Here, we use a catalogue of 900 relocated events developed by Sumy et al. (2017) that includes events from 7 November 2011 through the end of 2011 (Fig. 1). Manually determined P and S wave phase picks are available for each of the events, and absolute location errors are estimated to be ∼750 m horizontally and 1.3 km vertically (Sumy et al. 2017). Figure 1. Open in new tabDownload slide Map of the study area showing LC network seismic stations (blue squares) and earthquakes used in the study sized by magnitude (brown circles). The yellow stars and beach balls show the locations and USGS moment tensor, with (strike, dip, rake) given, of the three largest events in the sequence (M  4.8 foreshock, M  5.7 main shock and M  4.8 aftershock). The black lines show the Wilzetta fault system based on regional and local seismic and geological studies (Way 1983). The red inverted triangles are wastewater disposal wells sized based on the volume of fluid injection over their lifetime of operation until the end of 2010 into fault-bounded reservoirs (grey shaded regions), see Keranen et al. (2013) and references therein. An overview map is shown in the upper left with the location of the study area in the state of Oklahoma shown by the red box. Figure 1. Open in new tabDownload slide Map of the study area showing LC network seismic stations (blue squares) and earthquakes used in the study sized by magnitude (brown circles). The yellow stars and beach balls show the locations and USGS moment tensor, with (strike, dip, rake) given, of the three largest events in the sequence (M  4.8 foreshock, M  5.7 main shock and M  4.8 aftershock). The black lines show the Wilzetta fault system based on regional and local seismic and geological studies (Way 1983). The red inverted triangles are wastewater disposal wells sized based on the volume of fluid injection over their lifetime of operation until the end of 2010 into fault-bounded reservoirs (grey shaded regions), see Keranen et al. (2013) and references therein. An overview map is shown in the upper left with the location of the study area in the state of Oklahoma shown by the red box. Shear wave splitting estimation We measure SWS on all LC network stations to estimate shallow crustal velocity anisotropy. We use the automated Multiple Filter Automatic Splitting Technique (MFAST) developed by Savage et al. (2010) to determine the fast direction (ϕ) and the delay time (δt) for each source–receiver pair. MFAST implements a grid search to determine the (ϕ, δt) that minimizes the eigenvalue of the horizontal waveforms (Silver & Chan 1991). The minimum eigenvalue removes the effect of the anisotropy and restores the linearity of the shear wave. SWS parameters (ϕ, δt) can be dependent on both the frequency of the filter(s) applied and the measurement time window. MFAST uses a set of bandpass filters that can be customized for the expected frequency content of the shear waves depending on, for example, typical source–receiver distances in the data set as well as local site geology. The best filter is generally chosen to be the dominant (highest signal to noise) frequency within the measurement window. Here, we implement MFAST version 2.2 with the provided bandpass filter set appropriate for very local data. The 14 bandpass filters applied for very local data are listed in Table 1. We choose this implementation of MFAST because we have broadband data recorded at source–receiver distances of less than 20 km. Table 1. Set of filters used for very local events in the automated shear wave splitting code MFAST. Filter number . Low pass (Hz) . High pass (Hz) . 1 1 5 2 1 8 3 1 15 4 1 30 5 3 5 6 2 8 7 3 15 8 3 30 9 5 10 10 5 15 11 5 30 12 5 45 13 10 20 14 10 45 Filter number . Low pass (Hz) . High pass (Hz) . 1 1 5 2 1 8 3 1 15 4 1 30 5 3 5 6 2 8 7 3 15 8 3 30 9 5 10 10 5 15 11 5 30 12 5 45 13 10 20 14 10 45 Open in new tab Table 1. Set of filters used for very local events in the automated shear wave splitting code MFAST. Filter number . Low pass (Hz) . High pass (Hz) . 1 1 5 2 1 8 3 1 15 4 1 30 5 3 5 6 2 8 7 3 15 8 3 30 9 5 10 10 5 15 11 5 30 12 5 45 13 10 20 14 10 45 Filter number . Low pass (Hz) . High pass (Hz) . 1 1 5 2 1 8 3 1 15 4 1 30 5 3 5 6 2 8 7 3 15 8 3 30 9 5 10 10 5 15 11 5 30 12 5 45 13 10 20 14 10 45 Open in new tab MFAST applies an extensive list of criteria used to automatically grade the quality of the measurements that is expanded from the clustering analysis of Teanby et al. (2004). The quality criteria assess the stability of the (ϕ, δt) parameters across multiple measurement windows and frequencies. Quality assessments eliminate parameters for ray paths of the local events with an angle of incidence greater than 35° to ensure the measurements are within the shear wave window (Nuttli 1961; Booth & Crampin 1985). To determine the angle of incidence, the source to station ray paths are estimated from the TauP toolkit (Crotwell et al. 1999) with the local velocity model of Keranen et al. (2013). The angle between the measured fast direction and the initial source polarization is also considered; the shear wave will not be split into fast and slow shear polarizations if the initial polarization is nearly parallel or perpendicular to ϕ (Leary et al. 1990; Silver & Chan 1991). MFAST eliminates measurements if ϕ is outside of the range 20–70° from the initial polarization direction (Peng & Ben-Zion 2004; Savage et al. 2010). Here, we evaluate (ϕ, δt) evaluated to be of good quality (Quality A or B) using the default thresholds. We further select measurements with δt less than 0.20 s, although this condition only removes ∼1 per cent of the good quality measurements. Additional details on the MFAST quality measures can be found in Savage et al. (2010). The measured ϕ reflect the crustal velocity anisotropy along (at least a portion of) the path between the source and the receiver, and it can be difficult to isolate the location of the anisotropy along the path. We examine the spatial distribution of fast directions first using station-based compilations of fast directions for the case where the near-station anisotropy dominates the measurement. We also map the spatial variability of ϕ in two horizontal dimensions to represent the case where contributions of anisotropy may occur along the entire path. The Tomography Estimation and Shear-wave-splitting Spatial Average (TESSA) method developed by Johnson et al. (2011) computes the spatial average of ϕ on a grid. We use the quadtree implementation that allows us to have variable grid sizes depending on the density of ray paths in a region. We allow grid cells as small as 750 m using a gridding criterion that required between 100 and 400 ray paths in each grid cell. Decomposition of fast azimuth distributions The TESSA method (Johnson et al. 2011) generates an ensemble of rose diagrams that may contain tens of thousands of points and considerable scatter. To facilitate visual interpretation, TESSA also calculates a weighted-average azimuth in each grid cell. This approach may be adequate for many applications; however, experience shows that data often exhibit multiple modes in a single grid cell (Johnson et al. 2011). As we are interested in characterizing local deviations from the regional velocity anisotropy, we wish to highlight some of the higher-order modes in addition to the primary azimuth. One approach would be to select the highest peaks from a histogram, but adjacent peaks frequently obscure each other. To characterize multiple modes, we modify the TESSA codes to decompose the entire azimuth distribution in each grid cell into a finite set of component distributions. Our approach follows the Gaussian peak-fitting method of Brandon (1992), adapted to directional data. We envision the observed azimuth distribution as a composite probability density plot, which is the sum of the probability distributions that represent Nt individual azimuthal measurements and their uncertainties (Hurford et al. 1984). Our data, for the purpose of the decomposition, are the number of measurements per azimuthal bin. These values are given by the frequency function, $$\begin{eqnarray*} {\phi _o} \left( x \right) = \ {\rm{\Delta }}x\mathop \sum \limits_{i\ = \ 1}^{{N_t}} \frac{{{e^{\left[ {\sigma _i^{ - 2}\cos \left( {x - {\mu _i}} \right)} \right]}}}}{{2\pi {I_0}\left( {\sigma _i^{ - 2}} \right)}}{\rm{\ }}, \end{eqnarray*}$$(1) where Δx is the bin width, μ is a measured azimuth, σ2 is the variance of the azimuth measurement and I0 is the modified Bessel function of order zero. We assume that azimuth measurements are described by the Von Mises distribution, where μ is the location parameter and σ–2 is the concentration parameter (Mardia & Jupp 2009). The Von Mises distribution provides a good approximation to the wrapped normal distribution. We wish to fit the data with a model frequency function, Φp, that is constructed from the sum of Nf Von Mises distributions, where Nf << Nt, typically. Model functions have 3 Nf free parameters: |$\hat{\mu }$|⁠, |$\hat{\sigma }$| and |$\hat{n}$|⁠, the number of azimuths per component distribution. We find the best-fit parameters by minimizing the merit function, $$\begin{eqnarray*} {\chi ^2} = \sum \frac{{{{\left( {{\Phi _p}\left( x \right) - {\Phi _{o\left( x \right)}}} \right)}^2}}}{{{{\left( {s\left( x \right)} \right)}^2}}}, \end{eqnarray*}$$(2) where s is calculated from (1) by propagation of error (Brandon 1992). In order to identify an acceptable model, we iteratively fit the data with models for a range of Nf. For each iteration, we approximate the Bayesian information criterion, $$\begin{eqnarray*} 3{N_f}\ln {N_t} + {\chi ^2} \end{eqnarray*}$$(3) and select the model with the lowest value (Schwarz 1978). Although |${\chi ^2}$| only approximates the log-likelihood function when the level of noise is low, (3) performs well in our experience. For visual interpretation, we plot a grid of simplified rose diagrams showing the azimuths, |$\hat{\mu }$|⁠, and amplitudes, |$\hat{n}$|⁠, of the modeled component distributions. We restrict these plots to no more than the three highest amplitude azimuths for clarity of presentation. It is important to recall that the component distributions do not necessarily have physical meaning. A single component distribution might reflect two or more crack orientations that cannot be resolved with the available data, for example. Or, multiple component distributions might be wrongly inferred from a field of cracks with a single, noisy orientation. With this caveat in mind, the decomposition method may reveal important patterns in complex data. Template matching to enhance the catalogue To map the detailed fault structure, we extended the catalogue using a template matching method to identify smaller events that have waveforms similar to the original catalogue events. We applied the technique of Ross et al. (2017) to the continuous waveform data of the temporary LC network stations recorded between 2011/11/07 to 2012/01/31. We took as template events the relocated catalogue of Sumy et al. (2014) that was used in the shear wave splitting analysis described above. We used 2.0 s windows starting 0.5 s before the P or S phase as template waveforms. Note the P phase window was shortened to the S-P time if the S phase arrived before the end of the P phase window. Template waveforms were compared to daylong continuous data and cross correlation functions were summed across all phases and channels. We declared an event detection when at least nine channels of data exceeded a threshold of nine times the median absolute deviation. If more than one detection occurred within a 2.0 s window, the template with the largest average cross correlation was chosen. The parameters used to declare detections are comparable to previous studies that have similar station distributions (Ross et al. 2017; Cochran et al. 2018). Detections were relocated with the cluster-based relative relocation code GrowClust (Trugman & Shearer 2017). For the relocations, we use the same velocity model that was used to develop the original catalogue and relocate with HypoDD (Keranen et al. 2013; Sumy et al. 2017). Relative relocations were estimated from differential times calculated from waveform cross-correlation of pairs of events. We used the default parameters for GrowClust (Trugman & Shearer 2017), and include differential times for waveform pairs with a cross-correlation coefficient of 0.7 or larger. In the relocation scheme, the differential times are typically weighted by the cross-correlation value associated with each differential time. Here, we increased the weights of the template-template event pairs by a factor of 100. This implementation relocates the better-resolved template event first, and makes the relocations less dependent on the number and distribution of the detections (Ross et al. 2019). Location uncertainties are estimated using Growclust's bootstrap approach that is based on resampling the input P- and S-phase observation data. We implement 100 bootstrap iterations, which was found to provide stable location uncertainties (Trugman & Shearer 2017). Fault identification Using our relocated Prague catalogue, we identify seismogenic fault planes. Previously, 2-D fault traces were identified algorithmically in Oklahoma using a relative relocated template matched catalogue (Skoumal et al. 2019). This algorithm (FaultID) iteratively clusters seismicity and identifies linear features within those cluster that are assumed to be faults. Here, we expand the functionality of this FaultID algorithm to identify seismogenic fault planes (Fig. S1). First, earthquakes were spatially clustered using the density-based spatial clustering of applications with noise (DBSCAN) algorithm (Ester et al. 1996). With the application of this approach, an earthquake is classified as a core event if there are at least ten other earthquakes within the Euclidean distance Dpoint. Any earthquakes that are within distance Dpoint from a core point and had fewer than ten neighbors within distance Dpoint are considered border events. If a core event is within distance Dpoint from another core event, they are considered to be in the same cluster along with any corresponding border events. We repeat this clustering approach three times with the interevent distances (m) Dpoint in the set [1000, 500, 100]. Following each of these clustering steps, we search for planar trends within each cluster of earthquakes using the RANdom SAmple Consensus (RANSAC) algorithm (Fischler & Bolles 1981). With this approach, three earthquakes within a cluster are randomly selected and used to fit a plane. Inlier events are classified as earthquakes within a cluster that are within Euclidean distance Dplane. We repeat the process of randomly selecting three points and determining inliers 1000 times for each cluster. The model that contains the largest number of inliers for each cluster is selected, and the other models are discarded. Using the inliers in the preferred model, a least-squares regression is then performed and is used to represent the seismogenic fault plane. As this plane fitting step was performed following each of the three clustering steps, the plane fitting step was also performed three times using plane distances (m) Dplane in the set [300, 150, 100]. As a quality control step, fault planes represented by fewer than [500, 100, 20] earthquakes were discarded. Any earthquakes that have not been associated with a fault are reconsidered in the subsequent DBSCAN/RANSAC iterations. At the end of this processing, the resulting planes are then considered to represent the location and trend of seismogenic faults. RESULTS We used MFAST to measure splitting parameters for 900 events measured on the LC network stations. Note that a single waveform can yield multiple high-quality measurements because measurements are made across a range of bandpass filters. For our analysis we use only (ϕ, δt) parameters that are graded as A or B quality with δt less than 0.2 s, totaling 8569 (ϕ, δt) measurements (see Supplemental Information for SWS measurements). The majority of high-quality measurements are made at stations within 10 km of the fault (see Fig. 1) and at source–receiver distance of <15 km. Fig. 2(a) shows the distribution of ϕ versus source polarization. Most ϕ are between ∼50° and 100° from north and do not depend strongly on the initial source polarization. Similarly, we see no dependence of ϕ on angle of incidence (Fig. 2b). δt are generally scattered and do not show a dependence on source depth or hypocentral distance (Figs 2c and d), thus the location of anisotropy along the path is unknown. We focus our subsequent shear wave splitting analyses on examining the spatial variability of ϕ in the study region as we are interested in local stress orientations and leave more exhaustive exploration of δt to future work. Figure 2. Open in new tabDownload slide (a) Initial source polarizations versus ϕ are shown as grey dots. Only fast directions that are 20° to 70° from the initial source polarization are shown and used in subsequent analysis. (b) Angle of incidence versus ϕ; grey dots show individual measurements and larger black dots show 500 point moving average with a 50 per cent overlap and vertical bars show the standard error. (c) δt versus event depth; grey dots show individual measurements and larger black dots show 500 point moving average with a 50 per cent overlap and bars show the standard deviation. (d) δt versus hypocentral distance; grey dots show individual measurements and larger black dots show 500 point moving average with a 50 per cent overlap and bars show the standard deviation. Figure 2. Open in new tabDownload slide (a) Initial source polarizations versus ϕ are shown as grey dots. Only fast directions that are 20° to 70° from the initial source polarization are shown and used in subsequent analysis. (b) Angle of incidence versus ϕ; grey dots show individual measurements and larger black dots show 500 point moving average with a 50 per cent overlap and vertical bars show the standard error. (c) δt versus event depth; grey dots show individual measurements and larger black dots show 500 point moving average with a 50 per cent overlap and bars show the standard deviation. (d) δt versus hypocentral distance; grey dots show individual measurements and larger black dots show 500 point moving average with a 50 per cent overlap and bars show the standard deviation. Fig. 3 shows polar histograms of ϕ plotted on the station location where the measurements were made. Across the study area we generally observe that ϕ are largely oriented approximately east-west, or slightly east–northeast. The inset of Fig. 3 shows the polar histogram of all measurements that shows a dominant ϕ close to previously reported regional principal stress direction as N80E (Sumy et al. 2014) and N84E (Alt & Zoback 2017). However, there exists some spatial variability and secondary orientations across the study area. For example, we observe that ϕ measured at stations in the southeast portion of the study area show a subset ϕ oriented subparallel or parallel to the strike of the Wilzetta Fault (∼NE). Figure 3. Open in new tabDownload slide Polar histograms of ϕ plotted on the LC station locations. Inset shows the polar histogram of all ϕ measurements and the orange and magenta arrows show the direction of the principal compressive stress direction from Sumy et al. (2014) and Alt & Zoback (2017), respectively. Other symbols are the same as in Fig. 1. Figure 3. Open in new tabDownload slide Polar histograms of ϕ plotted on the LC station locations. Inset shows the polar histogram of all ϕ measurements and the orange and magenta arrows show the direction of the principal compressive stress direction from Sumy et al. (2014) and Alt & Zoback (2017), respectively. Other symbols are the same as in Fig. 1. The station-based compilations of ϕ are appropriate if the contribution of anisotropy near the station dominates, but the measured ϕ may also reflect the crustal velocity anisotropy along the entire path between the source and the receiver. We map the spatial variability of ϕ in two horizontal dimensions to represent the case where contributions of anisotropy may occur along the entire path. We use the quadtree gridding implementation of TESSA to divide the study region into a grid and determine the ray crossings in each grid. The resulting grid has smaller (750 m) cells close to the Prague sequence and larger (3 km) cells further from the sequence (Fig. 4). Figure 4. Open in new tabDownload slide Quadtree grid (black boxes) used to determine the local fast directions in local regions. Gridding criteria required between 100 and 400 ray paths in each box, with a minimum box size of 750 m. A total of 61 blocks were used to cover the region encompassed by the LC stations (blue squares). The red dots show regularly spaced nodes along the ray path between the station and the source. Polar histograms for the two grey shaded grid cells are shown in Fig. 5. The cyan box shows the approximate location of the study area (see Fig. 1). Figure 4. Open in new tabDownload slide Quadtree grid (black boxes) used to determine the local fast directions in local regions. Gridding criteria required between 100 and 400 ray paths in each box, with a minimum box size of 750 m. A total of 61 blocks were used to cover the region encompassed by the LC stations (blue squares). The red dots show regularly spaced nodes along the ray path between the station and the source. Polar histograms for the two grey shaded grid cells are shown in Fig. 5. The cyan box shows the approximate location of the study area (see Fig. 1). The original TESSA implementation uses the spatial average to determine the dominant azimuth in each grid cell; here we use Gaussian peak fitting to estimate the set of principal azimuths. Fig. 5 provides a comparison of the original TESSA spatial average and the principal azimuths measured by Gaussian peak fitting for two grid cells. We observe that both methods produce a similar result for cells where there is a single dominant ϕ (Fig. 5a). However, when the set of ϕ in a cell is more mixed, the principal azimuths determined by the modified Gaussian-fitting method better captures the variability. In fact, the average ϕ may do a poor job of representing the measurements; for example, the average may lie between two peaks in the distribution that has few-to-no actual ϕ observations (Fig. 5b). Note that a limitation of the Gaussian peak fitting method is that it may be more likely to fit a principal component azimuth to noise, so care should be taken not to overinterpret the output. Here we attempt to avoid including spurious peaks by only plotting up to three principal peaks and require that the plotted peaks have at least 20 measurements and standard deviations of less than 10°. Figure 5. Open in new tabDownload slide (a) Polar histogram for a grid cell that has a single dominant ϕ apparent in the measurements. This histogram is for a cell located on the west side of the study area (see west grey shaded box in Fig. 4). (b) Polar histogram for a grid cell located near the centre of the study area that has a more variable set of ϕ measurements (see east grey shaded box in Fig. 4). The dashed black bar shows the original TESSA calculation of average ϕ (Johnson et al. 2011) and the green bars show up to three principal azimuths from the Gaussian peak-fitting method, with length scaled by number of rays that were fit to determine the peak divided by the sum of the rays composing all of the fitted peaks in the grid cell. Figure 5. Open in new tabDownload slide (a) Polar histogram for a grid cell that has a single dominant ϕ apparent in the measurements. This histogram is for a cell located on the west side of the study area (see west grey shaded box in Fig. 4). (b) Polar histogram for a grid cell located near the centre of the study area that has a more variable set of ϕ measurements (see east grey shaded box in Fig. 4). The dashed black bar shows the original TESSA calculation of average ϕ (Johnson et al. 2011) and the green bars show up to three principal azimuths from the Gaussian peak-fitting method, with length scaled by number of rays that were fit to determine the peak divided by the sum of the rays composing all of the fitted peaks in the grid cell. Fig. 6(a) shows the principal azimuths determined across the study region (see Supplemental Information for principal azimuths, errors and weights). We observe that, to first order, the pattern is the same as we observed when fast directions were plotted at the station locations (Fig. 3). In 74 per cent of grid cells where principal azimuths are determined, we observe at least one principal azimuth oriented within one standard deviation of N84E (Fig. 6b). Some variation from ∼N84E is apparent between the southern portion of the main shock rupture plane and the mapped Wilzetta structure. For the set of grid cells along the Wilzetta structure, we consistently observe at least one principal azimuth that is subparallel to the previously mapped fault system (∼N45E, Fig. 6b). We do not observe a strong signature subparallel to the southern-most portion of the structure that ruptured during the M  5.7 main shock; a few grid cells show a secondary or tertiary component azimuth that is subparallel to the strike of the seismicity, but with the primary azimuth oriented ∼N84E. Supplemental Fig. S2 provides results for different sizes and distributions of grid cells to show the choice of grid does not significantly change the spatial distribution of principal azimuths. Figure 6. Open in new tabDownload slide (a) Principal azimuths (green bars) determined from the Gaussian peak-fitting method plotted at the centre of a grid cell (grid shown in Fig. 4). The length of each bar is scaled by the number of rays that were fit to determine the peak divided by the sum of the rays composing all of the fitted peaks in the grid cell. Longer bars indicate a single peak fits most of the measurements. Peaks are only shown if at least 20 measurements contribute to the peak and the standard deviation of the peak is less than 10°. The standard deviations of the azimuths are similar to, or just larger than, the width of the bars and are omitted for clarity. (b) Unweighted, principal azimuths that most closely aligned with the regional maximum compressive stress direction, N84E from Alt & Zoback ( 2017). Blue bars show the azimuths that are aligned within one standard deviation with the regional stress and red bars shows those that are significantly rotated. Other symbols are the same as in Fig. 1. Figure 6. Open in new tabDownload slide (a) Principal azimuths (green bars) determined from the Gaussian peak-fitting method plotted at the centre of a grid cell (grid shown in Fig. 4). The length of each bar is scaled by the number of rays that were fit to determine the peak divided by the sum of the rays composing all of the fitted peaks in the grid cell. Longer bars indicate a single peak fits most of the measurements. Peaks are only shown if at least 20 measurements contribute to the peak and the standard deviation of the peak is less than 10°. The standard deviations of the azimuths are similar to, or just larger than, the width of the bars and are omitted for clarity. (b) Unweighted, principal azimuths that most closely aligned with the regional maximum compressive stress direction, N84E from Alt & Zoback ( 2017). Blue bars show the azimuths that are aligned within one standard deviation with the regional stress and red bars shows those that are significantly rotated. Other symbols are the same as in Fig. 1. The 900 earthquakes in the Sumy et al. (2017) catalogue provide an image of the three major fault planes that ruptured during the Prague earthquake sequence. This catalogue also suggests there may be more complex structure, particularly along a length of ∼5 km centred on the M  5.7 main shock and extending north to M  4.8 foreshock epicentre and south to the epicentre of the M  4.8 aftershock (Fig. 1). We use template matching to identify additional events similar to those in the existing catalogue to better image the detailed fault structure of the Prague sequence. In total, the final catalogue contains 8811 events, including the 900 template events and 7911 relocated detections (see Supplemental Information for catalogue). Relocated detections have median relative uncertainties of 112 m horizontally and 133 m vertically. Fig. S3 shows an example of one catalogue template (Event     1000836, a M1.71 earthquake on 2011/11/10 05:08:26.02 UTC) and the set of 163 matching detections (magnitudes range from –0.8 to 1.86) recorded on the east channel at four LC network stations. The majority (∼94 per cent) of events occur between 1.5 and 6 km depth (Fig. S4A). We observe a moderate (∼8 per cent) number of templates and detections whose absolute depths fall within the sediments (Simpson and Arbuckle Groups) and locate above the Precambrian basement at 1.9 km depth (Keranen et al. 2013). This is a smaller percentage of events than have previously been identified as within the sediments, with McMahon et al. (2017) reporting as many as 40 per cent of the events above basement using a basement depth of 2.6 km. If we instead use a depth to basement of 2.6 km, we still only have 16 per cent of our events above the basement. The absolute locations and depths are strongly dependent on the velocity model used to locate events; here we use the velocity model of Keranen et al. (2013) where shallow velocities are constrained by well log information while McMahon et al. use a more generic Oklahoma velocity model. Almost 95 per cent of the original template events were identified in the first ∼10 days after the main shock with the remainder of the events identified before the end of 2011 [see Sumy et al. (2017) for details, Fig. S4B]. A majority of templates have fewer than ten detections, while a handful (∼10) of templates have over 100 detections each (Fig. S4c). Similar distributions of number of detections per template have been previously observed in template matching studies (Cochran et al. 2018). Across the entire spatial extent of the sequence we see template events with a moderate number of detections, but the templates with large numbers of detections appear to be clustered within ∼5 km epicentral distance from the main shock (Fig. S5). Fig. 7 shows the full, relocated catalogue of 8811 templates and detections coloured by time after the main shock. As also noted in Fig. S5, we observe large numbers of detected events within about 5 km of the main shock epicentre. Earthquakes span the full depth range of the sequence between ∼1 and 10 km depth in this location, with event depths shallower to the northeast and southwest along the main shock fault plane (Profile C-c in Fig. 7). Along cross-section B-b located at the northeastern end of the rupture, we observe what appears to be two shallow fault branches that connect to a single structure below ∼4 km depth. This cross-section is in a region of complex faulting identified by Joseph (1987) and Way (1983) and ∼2 km north of where the foreshock and main shock planes intersect. Cross-section A-a shows a single near-vertical fault plane is present for much of the southern half of the main shock rupture plane. The time evolution of the sequence suggests that the entire fault system was active immediately after the occurrence of the main shock. We do not observe any apparent growth of the aftershock zone through time either along strike or with depth. Figure 7. Open in new tabDownload slide Full catalogue, showing both templates and detections coloured by time in days after the main shock on 6 November 2011 at 3:53:10 in map and cross-section view. Green and yellow shading in cross-sections shows the approximate depth of the Hunton and Arbuckle Groups. Other symbols are as described on Fig. 1 . Figure 7. Open in new tabDownload slide Full catalogue, showing both templates and detections coloured by time in days after the main shock on 6 November 2011 at 3:53:10 in map and cross-section view. Green and yellow shading in cross-sections shows the approximate depth of the Hunton and Arbuckle Groups. Other symbols are as described on Fig. 1 . To examine the temporal evolution of the sequence and density of events in more detail we divide the main shock fault plane (Profile C-c in Fig. 7) into grid cells that are 0.5 km by 0.5 km. Of the grid cells that have at least one event, most of those cells include an event within the first day after the main shock (Fig. 8a). Over half (∼59 per cent) of grid cells with at least one earthquake are active for less than 45 d of 3-month catalogue (Fig. 8b). The grid cells that remain active over much of the study period tend to be located between the M  4.8 foreshock epicentre to the northeast and the M  4.8 aftershock epicentre to the southwest, and including the M  5.7 main shock epicentre, remain active through the end of our catalogue. This same region has the highest density of events with as many as 538 events occurring in a single grid cell (Fig. 8c). Overall, we observe the most prolific families, highest density of events, and longest active sequences are concentrated in a nest of seismicity close the main shock and between the foreshock and largest aftershock. Figure 8. Open in new tabDownload slide (a) Time in days since the main shock of the first event in each grid cell for cross-section C-c (Fig. 7) divided into 500 m by 500 m cells. (b) Same as in (a), except showing the time of the last event in each grid cell. (c) Event density (measured in number of events per 500 m by 500 m grid cell) for cross section C-c. Figure 8. Open in new tabDownload slide (a) Time in days since the main shock of the first event in each grid cell for cross-section C-c (Fig. 7) divided into 500 m by 500 m cells. (b) Same as in (a), except showing the time of the last event in each grid cell. (c) Event density (measured in number of events per 500 m by 500 m grid cell) for cross section C-c. We use the expanded catalogue and the FaultID algorithm to determine the set of faults activated in the Prague sequence. The fitted faults highlight the three primary rupture planes of the M  4.8 foreshock, M 5.7 main shock and M4.8 aftershock and 7 minor planes (Fig. 9, Supplementary Movie M1). Table 2 and Fig. 9 show the primary planes agree well with the moment tensor strikes and dips reported in the Advanced National Seismic System Comprehensive catalogue (U.S. Geological Survey 2019) and Global Centroid Moment Tensor catalogue (Dziewonski et al. 1981; Ekström et al. 2012). The minor fault planes show a range of orientations, with the majority having strikes that are subparallel to one of the three primary planes. The exception is a small fault striking almost perpendicular to the foreshock fault plane towards the northern end of the sequence that has a shallower dip (071°) than the remaining fault planes. The minor fault planes are all within the range of focal mechanisms previously reported (Sumy et al. 2014; McNamara et al. 2015; Walsh & Zoback 2016). Figure 9. Open in new tabDownload slide Fault planes (thick coloured lines) determined by the FaultID algorithm using the full catalogue are coloured by their deviation from the optimal orientation (°). The optimal fault plane orientation is assumed to be a fault striking 30° from the orientation of the local maximum horizontal stress. Local stress orientations are taken to be the azimuth closest to each 1 km segment of the fault that is most closely aligned with the regional stress direction (green bars are the primary azimuth with their standard deviations shown by dotted black lines). Note that faults with dips more than ±10° from vertical are not assessed (grey bars). Other symbols are as described in Fig. 1, except here seismicity is shown by open circles here. Figure 9. Open in new tabDownload slide Fault planes (thick coloured lines) determined by the FaultID algorithm using the full catalogue are coloured by their deviation from the optimal orientation (°). The optimal fault plane orientation is assumed to be a fault striking 30° from the orientation of the local maximum horizontal stress. Local stress orientations are taken to be the azimuth closest to each 1 km segment of the fault that is most closely aligned with the regional stress direction (green bars are the primary azimuth with their standard deviations shown by dotted black lines). Note that faults with dips more than ±10° from vertical are not assessed (grey bars). Other symbols are as described in Fig. 1, except here seismicity is shown by open circles here. Table 2. Fault orientations identified in this study compared with moment tensors for the M4.8 foreshock, M5.7 main shock and M4.8 aftershock. . Fitted plane . USGS moment tensor . Global centroid moment tensors . Foreshock N31°E, 091° N32°E, 098° N27°E, 107° Main shock N52°E, 083° N56°E, 095° N54°E, 088° Aftershock N84°E, 087° N88°E, 097° N91°E, 074° Fault 1 N48°E, 103° – – Fault 2 N84°E, 106° – – Fault 3 N38°E, 077° – – Fault 4 N25°E, 089° – – Fault 5 N44°E, 074° – – Fault 6 N126°E, 071° – – Fault 7 N41°E, 096° – – . Fitted plane . USGS moment tensor . Global centroid moment tensors . Foreshock N31°E, 091° N32°E, 098° N27°E, 107° Main shock N52°E, 083° N56°E, 095° N54°E, 088° Aftershock N84°E, 087° N88°E, 097° N91°E, 074° Fault 1 N48°E, 103° – – Fault 2 N84°E, 106° – – Fault 3 N38°E, 077° – – Fault 4 N25°E, 089° – – Fault 5 N44°E, 074° – – Fault 6 N126°E, 071° – – Fault 7 N41°E, 096° – – Open in new tab Table 2. Fault orientations identified in this study compared with moment tensors for the M4.8 foreshock, M5.7 main shock and M4.8 aftershock. . Fitted plane . USGS moment tensor . Global centroid moment tensors . Foreshock N31°E, 091° N32°E, 098° N27°E, 107° Main shock N52°E, 083° N56°E, 095° N54°E, 088° Aftershock N84°E, 087° N88°E, 097° N91°E, 074° Fault 1 N48°E, 103° – – Fault 2 N84°E, 106° – – Fault 3 N38°E, 077° – – Fault 4 N25°E, 089° – – Fault 5 N44°E, 074° – – Fault 6 N126°E, 071° – – Fault 7 N41°E, 096° – – . Fitted plane . USGS moment tensor . Global centroid moment tensors . Foreshock N31°E, 091° N32°E, 098° N27°E, 107° Main shock N52°E, 083° N56°E, 095° N54°E, 088° Aftershock N84°E, 087° N88°E, 097° N91°E, 074° Fault 1 N48°E, 103° – – Fault 2 N84°E, 106° – – Fault 3 N38°E, 077° – – Fault 4 N25°E, 089° – – Fault 5 N44°E, 074° – – Fault 6 N126°E, 071° – – Fault 7 N41°E, 096° – – Open in new tab Given the local stress orientations determined from shear wave splitting, we can determine whether each of the determined fault planes is optimally oriented for slip. We assume the stress field remains constant over the duration of the sequence, as suggested by the consistency between the local fast directions and the regional stress direction (Fig. 6b). It is possible that smaller scale changes in the stress occurred that we are unable to resolve with shear wave splitting, however most of the faults we evaluate are similar or larger in scale to the resolution of the shear wave splitting estimates. Since most of the faults are near-vertically oriented, we assume a 90° dip and a coefficient of friction of 0.6 such that faults striking 30° from the maximum (horizontal) principal stress are optimally orientated to slip, similar to the method of Skoumal et al. (2019). Fig. 9 shows the deviation of individual (100 m long) fault segments from the optimal orientation, considering the nearest estimated primary azimuth with orientation most similar to the regional stress direction (see Fig. 6). We only consider faults with dips within 10° from vertical. Faults that have deviations close to zero are optimally oriented for slip, while those with larger deviations are expected to have a lower probability of failure. We find that the main shock plane is well oriented for failure across its entire length with deviations from optimal orientation that range from 1° to 4°. In contrast, the M4.8 foreshock and M4.8 aftershock fault planes are both unfavourably oriented in the local stress field. The M4.8 foreshock fault plane deviates between 20° and 25° from the optimal orientation. And, the M4.8 aftershock is more significantly deviated from optimal orientation, with deviations between 26° and 29°. DISCUSSION Disposal-induced earthquakes, including the 2011 Prague sequence in Oklahoma, have in some cases caused damage to structures in the epicentral region (Keranen et al. 2013) as well as concerns about the safety of critical energy infrastructure (McNamara et al. 2015). Further, there may be less immediately visible impacts to local populations, including decreased housing prices (Cheung et al. 2018) as well as positive correlations of earthquake rates with vehicular collisions (Casey et al. 2019) and Google searches related to anxiety (Casey et al. 2018). Thus, it is critical to understand the occurrence of moderate induced earthquake sequences with the hope of better forecasting future events. Several current forecasting methods use regional stress field information and orientations of mapped faults in order to predict which structures are more (or less) hazardous in the presence of local fluid disposal (Alt & Zoback 2017; Snee & Zoback 2018). And, a number of studies show many of the active structures in Oklahoma are generally optimally oriented in the regional stress field (McNamara et al. 2015; Walsh & Zoback 2016; Schoenball & Ellsworth 2017; Skoumal et al. 2019). However, some studies have suggested slip may sometimes occur on faults that are unfavourably oriented, at least in the regional stress field (Sumy et al. 2014; Skoumal et al. 2019). Here, we determine the local stress field and detailed fault sequence of the 2011 M5.7 Prague earthquake sequence to determine whether the set of activated fault structures are well oriented for failure. We made 8569 high-quality shear wave splitting measurements (ϕ, δt) on the temporary stations deployed around the Prague sequence. We find that ϕ measurements are predominantly oriented N84°E, parallel to the regional principal stress orientation estimated from focal mechanism inversions (Sumy et al. 2014; Walsh & Zoback 2016). Fracture orientations striking predominantly N75°E were mapped in the upper 60 m in northcentral Oklahoma (Queen & Rizer 1990), which E. Liu et al. (1991) suggest are oriented close to the maximum compressive stress direction and may result in anisotropic permeability. Thus, we infer that the majority of determined ϕ are sensitive to, and aligned with, the maximum horizontal stress. We are unable to constrain stress orientations to a specific depth range and the variability in measured δt does not allow us to constrain the depth extent of the anisotropy. The consistency between ϕ and the principal stress orientation determined from focal mechanism inversions, that reflect stresses at source depths, allows us to infer that the stress orientations do not change significantly with depth in this region. The dominant azimuths measured from the Gaussian peak fitting procedure in most of the grid cells are nearly uniformly, consistent with the regional stress direction (∼N84°E) across the entire study area. We observe no significant change in the principal horizontal stress direction for grid cells directly adjacent to the M5.7 main shock rupture. The existence or lack of rotations of the principal stress axes following moderate or large ruptures can tell us about the strength of the crust relative to the stress drop of the earthquake (Hardebeck & Okada 2018). If the stress drop is of the same order as, or a large fraction of, the crustal strength then one would expect to see a rotation in fast directions. However, if the stress drop represents only a small fraction of the total strength of the crust, then no such rotation should occur. A few studies have suggested a rotation of the fast directions following moderate earthquakes measured from focal mechanism inversions and shear wave splitting (Crampin et al. 1990; Hauksson 1994). However, as noted by Peng & Ben-Zion (2004), extreme care should be taken to ensure that changing source–receiver paths are not the cause of an apparent temporal change in ϕ for shear wave splitting studies. Further, several studies have looked for, but failed to observe, evidence of local changes in fast directions close to recent moderate earthquakes (Cochran et al. 2006; Liu et al. 2008). Here, we do not have sufficient pre-event data to determine whether the orientations change due to the Prague earthquake sequence, but the observed consistency between local and regional stress orientations suggest no significant changes to the stress field occurred. We also see no obvious influence of local injection on the stress orientations at the northeastern end of the sequence that has been implicated in triggering (Keranen et al. 2013). This is in contrast to a study in Kansas that suggested temporal rotations in ϕ resulting from injection (Nolte et al. 2017), although the study does not adequately control for differing source–receiver paths and their observations may instead reflect spatially heterogeneous ϕ. Further, a high-resolution study of anisotropy at a hydraulic fracturing site showed ϕ remained constant during and after the injection stages, but observed a variation in fracture densities (δt) (Farghal 2018). In addition to the dominant ϕ oriented at N84°E, we observe a subset of fast directions that are parallel to mapped faults or seismicity lineaments. There is a stronger fault-parallel signature for grid cells close to the previously mapped Wilzetta fault system (Way 1983; Joseph 1987). This may indicate that the Wilzetta fault has a more well-developed shear fabric, and thus may be a more mature structure than the southwest portion of the main shock fault. Alternately, the shear wave splitting may be more sensitive to structures in the sedimentary units, which may suggest the southwestern portion of main shock fault does not extend updip into sedimentary groups as also suggested by a lack of mapped structures in the region. The extended catalogue of events provides a high-resolution view of the seismicity over the 90 d following the main shock rupture. We note that the original catalogue of template events was biased, in that significant attention was given to identifying events in the first 10 d or so after the main shock (see Sumy et al. (2017) for details). This initial catalogue includes events distributed across the three primary rupture planes, so we assume here that the full catalogue adequately captures events across these major features. The greatest density of earthquakes is observed to occur where three primary rupture planes of the M4.8 foreshock, M5.7 main shock and M4.8 aftershock intersect. In this same active region, the seismicity extends across the full depth range of the sequence and includes both the shallowest (∼1 km, in the Arbuckle Group) and deepest (∼10 km depth) earthquakes. The area of vigorous seismicity is located where the mapped Wilzetta system is comprised of a complex set of faults that have been interpreted to act as boundaries to fluid flow (Keranen et al. 2013). Way (1983) reports that the Wilzetta fault system is a ridge structure comprised of a series of en echelon transverse faults rather than a single fracture. Plane fitting to the seismicity also suggests a small fault segment subparallel to the main shock fault plane in this same region. And, to the north, the seismicity appears to become sparser past the location of a cross fault that bisects the foreshock rupture plane. While the M5.7 main shock fault plane is optimally oriented for failure in the background stress field, we find that neither the M4.8 foreshock nor the M4.8 aftershock planes are optimally oriented to fail in the local stress field with deviations from optimal between 20° and 29°. Faults that deviate more than ±15° from optimal orientation are considered unfavourably oriented (Sibson 1990). Keranen et al. (2013) suggested the sequence was likely induced due to injection into closed fault compartments that caused increased fluid pressures on the foreshock fault. Our finding that the foreshock fault is poorly oriented may further suggest that elevated fluid pressures may have induced failure of this structure. The foreshock is subparallel to, and thought to align with, faults along the mapped Wilzetta fault that extends updip into the Arbuckle and Hunton Groups (Way 1983). The observation of fault parallel-aligned fast shear wave orientations is suggestive of shear fabric and may suggest a fault damage zone exists along the Wilzetta fault. Fault damage zones have been shown to have higher permeabilities and may act as conduits for fluid flow parallel to faults (Zhang & Sanderson 1995; Wibberley & Shimamoto 2003). Thus, we surmise that the Wilzetta fault may have provided a conduit for fluids, leading to increasing fluid pressure in the Precambrian basement. Our results also suggest that the aftershock fault plane is not optimally oriented for slip in the local stress field. Additionally, the static stress changes from the main shock imparted on M4.8 aftershock epicentre discouraged failure (Sumy et al. 2014). Norbeck & Horne (2016) invoke transient fluid flow resulting in increased pore fluid pressures to explain the occurrence of the largest aftershock. The penetration of fluid along the eventual plane of the M4.8 aftershock may be supported by observations of preceding activity within the damage zone of this structure (Savage et al. 2017). The Prague sequence defines ten primary fault planes that have a wide range of strikes and dips, showing that a heterogeneous set of structures were activated in the sequence. Of the five near-vertical structures that we evaluate within the context of the regional stress field, we find that only the main shock and an isolated fault segment located in the SE portion of the study area are well-oriented for failure. The remaining three structures that host events as large as M4.8 are not optimally oriented to slip. Our results confirm that even structures that are not optimally oriented in the regional, or local, stress field may host moderate-sized induced events (M∼5). Some recent studies have presented assessments of the likelihood of inducing events based on the orientation of previously identified faults in the regional stress field, often described as the fault slip potential (Walsh & Zoback 2016; Alt & Zoback 2017; Snee & Zoback 2018). The activation of unfavourably oriented faults during the Prague sequence suggests the hazard potential of a given structure may be difficult to quantify without additional information on pre-existing fault stress and hydromechanical considerations. CONCLUSIONS We investigate the 2011 Prague earthquake sequence to determine the local principal stress directions, activated fault structures and other characteristics of the sequence. We find that fast orientations measured via shear wave splitting primarily reflect the maximum horizontal stress direction, with secondary orientations suggesting the presence of shear fabric aligned subparallel to the Wilzetta fault. We find the principal horizontal stress directions are uniform across spatial scales as small as 750 m and are consistent with regional estimates from focal mechanism inversions (Sumy et al. 2014; Walsh & Zoback 2016; Alt & Zoback 2017). We find a majority of aftershocks occur within a 5-km-long portion of the sequence near the M5.7 main shock epicentre that is bounded to the southwest by the epicentre of the largest aftershock (M4.8) and to the northeast by the epicentre of a M4.8 foreshock. We identify ten distinct fault segments that are activated during the sequence. The faults have a wide array of orientations, suggesting the activation of a complex fault network. We find that the M5.7 main shock fault is optimally oriented, but both the M4.8 foreshock and M4.8 aftershock occur on fault planes that are deviated from the optimal orientation for slip. Our results suggest that elevated pore fluid pressures were likely required for activation of these unfavourably oriented planes. This sequence shows that slip can occur on faults not well-oriented for failure in the local stress field, and care should be taken when forecasting hazard from estimates of fault slip potential. SUPPORTING INFORMATION Figure S1. Algorithm flowchart describing the FaultID method. N represents the minimum number of events to define a cluster and D represents the maximum interevent distance (m) for two points to be considered in the same cluster and P represents the minimum number of inlier events for a plane to be accepted. DBSCAN, density-based spatial clustering of applications with noise; EQ, earthquake; QC, quality control; RANSAC, RANdom SAmple Consensus. Figure S2. Spatial distribution of primary azimuths using different grid configurations. Same as Fig. 6 in the main text except for grids configured as: (left-hand panel) quadtree grid with minimum grid size of 500 m, (middle panel) quadtree grid with minimum grid size of 1500 m and (right-hand panel) uniform grid with 1000 m grid spacing. Figure S3. Plot of matching detections associated with Template 1 000 836 as recorded on CHE channel for stations: (A) LC02, (B) LC03, (C) LC04 and (D) LC05. Each horizontal line in the plot is one detection, and red/blue indicates waveform polarity. Traces are aligned on the Pwave but the main energy that is observed in these figures is from the Swave. Figure S4. (a) Stacked histogram by depth of detections (blue) and templates (grey). Green and yellow bars show approximate location of the Simpson and Arbuckle Groups, respectively. (b) Stacked histogram of number of detections (blue) and templates (grey) versus time (days) after the main shock. (c) Histogram showing the number of detections per template event. Figure S5. Template events coloured by the number of matching detections identified with the template matching procedure shown in map view (lower left) and in three cross-sections A-a (perpendicular to the mid-section of the main shock fault), B-b (perpendicular to the foreshock fault and the northern end of the main shock fault), and C-c (parallel to the main shock fault). The cross-sections include all events that occur within 2 km of the profile. Mapview symbols are the same as Fig. 1 of the main text. Cross sections also show the approximate location of the Simpson (green band) and Arbuckle (yellow band) groups. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. ACKNOWLEDGEMENTS We thank Daniel Trugman for making the GrowClust codes available and providing advice on an appropriate weighting scheme for relocation of the template and detected events. We acknowledge D. Sumy and C. Neighbors for useful discussions of early SWS findings. We thank two anonymous journal reviewers and USGS reviewers N. Farghal and J. Hardebeck for constructive reviews. The Generic Mapping Tools software was used to generate several of the figures in this paper (Wessel & Smith 1991; Wessel et al. 2013). Any use of trade, firm or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. REFERENCES Alt R.C. , Zoback M.D., 2017 . In situ stress and active faulting in Oklahoma , Bull. seism. Soc. Am. , 107 , 216 – 228 . 10.1785/0120160156 Google Scholar Crossref Search ADS WorldCat Crossref Aster R.C. , Shearer P.M., Berger J., 1990 . Quantitative measurements of shear wave polarizations at the Anza seismic network, Southern California: implications for shear wave splitting and earthquake prediction , J. geophys. Res. , 95 , 12 449 – 12 473 . Google Scholar Crossref Search ADS WorldCat Balfour N.J. , Savage M.K., Townend J., 2005 . Stress and crustal anisotropy in Marlborough, New Zealand: evidence for low fault strength and structure-controlled anisotropy , Geophys. J. Int , 163 , 1073 – 1086 . 10.1111/j.1365-246X.2005.02783.x Google Scholar Crossref Search ADS WorldCat Crossref Barbour A.J. , Xue L., Roeloffs E., Rubinstein J.L., 2019 . Leakage and increasing fluid pressure detected in Oklahoma's wastewater disposal reservoir , J. geophys. Res. , 124 , 2896 – 2919 . 10.1029/2019JB017327 Google Scholar Crossref Search ADS WorldCat Crossref Barton C.A. , Zoback M.D., 1994 . Stress perturbations associated with active faults penetrated by boreholes: possible evidence for near-complete stress drop and a new technique for stress magnitude measurement , J. geophys. Res. , 99 , 9373 – 9390 . 10.1029/93JB03359 Google Scholar Crossref Search ADS WorldCat Crossref Bohnhoff M. , Grosser H., Dresen G., 2006 . Strain partitioning and stress rotation at the North Anatolian fault zone from aftershock focal mechanisms of the 1999 Izmit Mw = 7.4 earthquake , Geophys. J. Int. , 166 , 373 – 385 . 10.1111/j.1365-246X.2006.03027.x Google Scholar Crossref Search ADS WorldCat Crossref Boness N.L. , Zoback M.D., 2006 . Mapping stress and structurally controlled crustal shear velocity anisotropy in California , Geology , 34 , 825 – 828 . 10.1130/G22309.1 Google Scholar Crossref Search ADS WorldCat Crossref Booth D.C. , Crampin S., 1985 . Shear-wave polarizations on a curved wavefront at an isotropic free surface , Geophys. J. Int. , 83 , 31 – 45 . 10.1111/j.1365-246X.1985.tb05154.x Google Scholar Crossref Search ADS WorldCat Crossref Boyd O.S. , McNamara D.E., Hartzell S., Choy G., 2017 . Influence of lithostatic stress on earthquake stress drops in North America , Bull. seism. Soc. Am. , 107 , 856 – 868 . 10.1785/0120160219 Google Scholar Crossref Search ADS WorldCat Crossref Brandon M.T. , 1992 . Decomposition of fission-track grain-age distributions , Am. J. Sci. , 292 , 535 – 564 . 10.2475/ajs.292.8.535 Google Scholar Crossref Search ADS WorldCat Crossref Casey J.A. , Elser H., Goldman-Mellor S., Catalano R., 2019 . Increased motor vehicle crashes following induced earthquakes in Oklahoma, USA , Sci. Total Environ. , 650 , 2974 – 2979 . 10.1016/j.scitotenv.2018.10.043 Google Scholar Crossref Search ADS PubMed WorldCat Crossref Casey J.A. , Goldman-Mellor S., Catalano R., 2018 . Association between Oklahoma earthquakes and anxiety-related Google search episodes , Environ. Epidemiol. , 2 , e016 , doi:10.1097/EE9.0000000000000016 . 10.1097/EE9.0000000000000016 Google Scholar Crossref Search ADS WorldCat Crossref Chen X. et al. ., 2017 . The Pawnee earthquake as a result of the interplay among injection, faults and foreshocks , Sci. Rep. , 7 , 4945 , doi:10.1038/s41598-017-04992-z . 10.1038/s41598-017-04992-z Google Scholar Crossref Search ADS PubMed WorldCat Crossref Cheung R. , Wetherell D., Whitaker S., 2018 . Induced earthquakes and housing markets: Evidence from Oklahoma , Reg. Sci. Urban. Econ. , 69 , 153 – 166 . 10.1016/j.regsciurbeco.2018.01.004 Google Scholar Crossref Search ADS WorldCat Crossref Cochran E.S. , Kroll K.A., 2015 . Stress- and structure-controlled anisotropy in a region of complex faulting—Yuha Desert, California , Geophys. J. Int. , 202 , 1109 – 1121 . 10.1093/gji/ggv191 Google Scholar Crossref Search ADS WorldCat Crossref Cochran E.S. , Li Y.-G., Vidale J.E., 2006 . Anisotropy in the shallow crust observed around the San Andreas Fault before and after the 2004 M 6.0 Parkfield earthquake , Bull. seism. Soc. Am. , 96 , S364 – S375 . 10.1785/0120050804 Google Scholar Crossref Search ADS WorldCat Crossref Cochran E.S. , Ross Z.E., Harrington R.M., Dougherty S.L., Rubinstein J.L., 2018 . Induced earthquake families reveal distinctive evolutionary patterns near disposal wells , J. geophys. Res. , 123 , 8045 – 8055 . 10.1029/2018JB016270 Google Scholar Crossref Search ADS WorldCat Crossref Cochran E.S. , Vidale J.E., Li Y.-G., 2003 . Near-fault anisotropy following the Hector Mine earthquake , J. geophys. Res. , 108 , doi:10.1029/2002JB002352 . 10.1029/2002JB002352 OpenURL Placeholder Text WorldCat Crossref Crampin S. , 1985 . Evaluation of anisotropy by shear-wave splitting , Geophysics , 50 , 142 – 152 . 10.1190/1.1441824 Google Scholar Crossref Search ADS WorldCat Crossref Crampin S. , Booth D.C., Evans R., Peacock S., Fletcher J.B., 1990 . Changes in shear wave splitting at Anza near the time of the North Palm Springs Earthquake , J. geophys. Res. , 95 , 11 197 – 11 212 . 10.1029/JB095iB07p11197 Google Scholar Crossref Search ADS WorldCat Crossref Crotwell H.P. , Owens T.J., Ritsema J., 1999 . The TauP Toolkit: flexible seismic travel-time and ray-path Utilities , Seismol. Res. Lett. , 70 , 154 – 160 . 10.1785/gssrl.70.2.154 Google Scholar Crossref Search ADS WorldCat Crossref Dziewonski A.M. , Chou T.-A., Woodhouse J.H., 1981 . Determination of earthquake source parameters from waveform data for studies of global and regional seismicity , J. geophys. Res. , 86 , 2825 – 2852 . 10.1029/JB086iB04p02825 Google Scholar Crossref Search ADS WorldCat Crossref Ekström G. , Nettles M., Dziewoński A.M., 2012 . The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes , Phys. Earth planet. Inter. , 200–201 , 1 – 9 . 10.1016/j.pepi.2012.04.002 Google Scholar Crossref Search ADS WorldCat Crossref Ellsworth W.L. , 2013 . Injection-induced earthquakes , Science , 341 , 1 225 942 – 1 225 942 . 10.1126/science.1225942 Google Scholar Crossref Search ADS WorldCat Crossref Ester M. , Kriegel H.-P., Sander J., Xu X., 1996 . A density-based algorithm for discovering clusters a density-based algorithm for discovering clusters in large spatial databases with noise , in Proceedings of the Second International Conference on Knowledge Discovery and Data Mining KDD’96 , pp. 226 – 231 ., AAAI Press . Retrieved from http://dl.acm.org/citation.cfm?id=3001460.3001507 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Farghal N.S. , 2018 . Fault and Fracture Identification and Characterization in 3D Seismic Data From Unconventional Reservoirs , Stanford University . Retrieved from http://purl.stanford.edu/yg354fd2885 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Fischler M.A. , Bolles R.C., 1981 . Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography , Commun. ACM , 24 , 381 – 395 . 10.1145/358669.358692 Google Scholar Crossref Search ADS WorldCat Crossref Goebel T.H.W. , Brodsky E.E., 2018 . The spatial footprint of injection wells in a global compilation of induced earthquake sequences , Science , 361 , 899 – 904 . 10.1126/science.aat5449 Google Scholar Crossref Search ADS PubMed WorldCat Crossref Goertz-Allmann B.P. , Goertz A., Wiemer S., 2011 . Stress drop variations of induced earthquakes at the Basel geothermal site , Geophys. Res. Lett. , 38 , doi:10.1029/2011GL047498 . 10.1029/2011GL047498 OpenURL Placeholder Text WorldCat Crossref Ham W.E. , 1973 . Regional Geology of the Arbuckle Mountains, Oklahoma. Oklahoma Geological Survey Field Guidebook , pp. 1 – 56 ., Norman . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Handin J. , Hager R.V. Jr, Friedman M., Feather J.N., 1963 . Experimental deformation of sedimentary rocks under confining pressure: pore pressure tests , AAPG Bull. , 47 , 717 – 755 . OpenURL Placeholder Text WorldCat Hardebeck J.L. , Michael A.J., 2004 . Stress orientations at intermediate angles to the San Andreas Fault, California , J. geophys. Res. , 109 , doi:10.1029/2004JB003239 . 10.1029/2004JB003239 OpenURL Placeholder Text WorldCat Crossref Hardebeck J.L. , Okada T., 2018 . Temporal stress changes caused by earthquakes: a review , J. geophys. Res. , 123 , 1350 – 1365 . 10.1002/2017JB014617 Google Scholar Crossref Search ADS WorldCat Crossref Hauksson E. , 1994 . State of stress from focal mechanisms before and after the 1992 Landers earthquake sequence , Bull. seism. Soc. Am. , 84 , 917 – 934 . OpenURL Placeholder Text WorldCat Healy J.H. , Rubey W.W., Griggs D.T., Raleigh C.B., 1968 . The Denver earthquakes , Science , 161 , 1301 – 1310 . Google Scholar Crossref Search ADS PubMed WorldCat Hickman S. , Zoback M., 2004 . Stress orientations and magnitudes in the SAFOD pilot hole , Geophys. Res. Lett. , 31 , doi:10.1029/2004GL020043 . 10.1029/2004GL020043 OpenURL Placeholder Text WorldCat Crossref Holland A. , 2011 . Examination of Possibly Induced Seismicity from Hydraulic Fracturing in the Eola Field, Garvin County, Oklahoma (No. OF1-2011) , pp. 1 – 28 ., Oklahoma Geological Survey . Retrieved from http://theweeks.org/tmp/FILES/AustinHollandsEarthquakePaperAndFrackingOF1_2011.pdf . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Holt R.A. , Savage M.K., Townend J., Syracuse E.M., Thurber C.H., 2013 . Crustal stress and fault strength in the Canterbury Plains, New Zealand , Earth planet. Sci. Lett. , 383 , 173 – 181 . 10.1016/j.epsl.2013.09.041 Google Scholar Crossref Search ADS WorldCat Crossref Huang Y. , Ellsworth W.L., Beroza G.C., 2017 . Stress drops of induced and tectonic earthquakes in the central United States are indistinguishable , Sci. Adv. , 3 , e1700772 , doi:10.1126/sciadv.1700772 . 10.1126/sciadv.1700772 Google Scholar Crossref Search ADS PubMed WorldCat Crossref Hubbert M.K. , Rubey W.W., 1959 . Role of fluid pressure in mechanics of overthrust faulting: I. Mechanics of fluid-filled porous solids and its application to overthrust faulting , Bull. geol. Soc. Am. , 70 , 115 – 166 . Google Scholar Crossref Search ADS WorldCat Hurford A.J. , Fitch F.J., Clarke A., 1984 . Resolution of the age structure of the detrital zircon populations of two Lower Cretaceous sandstones from the Weald of England by fission track dating , Geol. Mag. , 121 , 269 – 277 . Google Scholar Crossref Search ADS WorldCat Johnson J.H. , Savage M.K., Townend J., 2011 . Distinguishing between stress-induced and structural anisotropy at Mount Ruapehu volcano, New Zealand , J. geophys. Res. , 116 , doi:10.1029/2011JB008308 . 10.1029/2011JB008308 OpenURL Placeholder Text WorldCat Crossref Jones L.M. , 1988 . Focal mechanisms and the state of stress on the San Andreas Fault in southern California , J. geophys. Res. , 93 , 8869 – 8891 . Google Scholar Crossref Search ADS WorldCat Joseph L. , 1987 . Subsurface analysis, “Cherokee” Group (Des Monesian), portions of Lincoln, Pottawatomie, Seminole, and Okfuskee counties, Oklahoma , Shale Shaker , 12 , 44 – 69 . OpenURL Placeholder Text WorldCat Kaneshima S. , 1990 . Origin of crustal anisotropy: Shear Wave splitting studies in Japan , J. geophys. Res. , 95 , 11 121 – 11 133 . Google Scholar Crossref Search ADS WorldCat Keranen K.M. , Savage H.M., Abers G.A., Cochran E.S., 2013 . Potentially induced earthquakes in Oklahoma, USA: Links between wastewater injection and the 2011 Mw 5.7 earthquake sequence , Geology , 41 , 699 – 702 . Google Scholar Crossref Search ADS WorldCat Keranen K.M. , Weingarten M., Abers G.A., Bekins B.A., Ge S., 2014 . Sharp increase in central Oklahoma seismicity since 2008 induced by massive wastewater injection , Science , 345 , 448 – 451 . Google Scholar Crossref Search ADS PubMed WorldCat Kroll K.A. , Cochran E.S., Murray K.E., 2017 . Poroelastic properties of the Arbuckle Group in Oklahoma derived from well fluid level response to the 3 September 2016 Mw 5.8 Pawnee and 7 November 2016 Mw 5.0 Cushing earthquakes , Seismol. Res. Lett. , 88 , 963 – 970 . Google Scholar Crossref Search ADS WorldCat Leary P.C. , Crampin S., McEvilly T.V., 1990 . Seismic fracture anisotropy in the Earth's crust: An overview , J. geophys. Res. , 95 , 11105 – 11114 . Google Scholar Crossref Search ADS WorldCat Liu E. , Crampin S., Queen J.H., 1991 . Fracture detection using crosshole surveys and reverse vertical seismic profiles at the Conoco Borehole Test Facility, Oklahoma , Geophys. J. Int. , 107 , 449 – 463 . Google Scholar Crossref Search ADS WorldCat Liu Y. , Zhang H., Thurber C.H., Roecker S., 2008 . Shear wave anisotropy in the crust around the San Andreas fault near Parkfield: spatial and temporal analysis , Geophys. J. Int. , 172 , 957 – 970 . Google Scholar Crossref Search ADS WorldCat Liu Y. , Zhang H., Zhang X., Pei S., An M., Dong S., 2015 . Anisotropic upper crust above the aftershock zone of the 2013 Ms 7.0 Lushan earthquake from the shear wave splitting analysis , Geochem. Geophys. Geosyst. , 16 , 3679 – 3696 . Google Scholar Crossref Search ADS WorldCat Llenos A.L. , Michael A.J., 2013 . Modeling earthquake rate changes in Oklahoma and Arkansas: possible signatures of induced seismicity , Bull. seism. Soc. Am. , 103 , 2850 – 2861 . Google Scholar Crossref Search ADS WorldCat Mardia K.V. , Jupp P.E., 2009 . Directional Statistics , John Wiley & Sons . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC McGarr A. , Barbour A.J., 2017 . Wastewater disposal and the earthquake sequences during 2016 Near Fairview, Pawnee, and Cushing, Oklahoma: induced earthquake sequences during 2016 , Geophys. Res. Lett. , 44 , 9330 – 9336 . Google Scholar Crossref Search ADS WorldCat McMahon N.D. , Aster R.C., Yeck W.L., McNamara D.E., Benz H.M., 2017 . Spatiotemporal evolution of the 2011 Prague, Oklahoma, aftershock sequence revealed using subspace detection and relocation: Prague, OK, Aftershock Sequence , Geophys. Res. Lett. , 44 , 7149 – 7158 . Google Scholar Crossref Search ADS WorldCat McNamara D.E. et al. ., 2015 . Earthquake hypocenters and focal mechanisms in central Oklahoma reveal a complex system of reactivated subsurface strike-slip faulting: Earthquake Source Parameters in Oklahoma , Geophys. Res. Lett. , 42 , 2742 – 2749 . Google Scholar Crossref Search ADS WorldCat McNamara D.E. et al. ., 2015 . Reactivated faulting near Cushing, Oklahoma: increased potential for a triggered earthquake in an area of United States strategic infrastructure , Geophys. Res. Lett. , 42 , 8328 – 8332 . Google Scholar Crossref Search ADS WorldCat Morgan B.C. , Murray K.E., 2015 . Characterizing small-scale permeability of the Arbuckle Group, Oklahoma (Oklahoma Geological Survey Open-File Report No. OF2-2015) , pp. 1 – 32 . Murray K.E. , 2015 . Class II saltwater disposal for 2009–2014 at the annual-, state-, and county-scales by geologic zones of completion, Oklahoma (Oklahoma Geological Survey Open-File Report No. OF5-2015) , pp. 1 – 12 . Nolte K.A. , Tsoflias G.P., Bidgoli T.S., Watney W.L., 2017 . Shear-wave anisotropy reveals pore fluid pressure–induced seismicity in the U.S. midcontinent , Sci. Adv. , 3 , e1700443 , doi:10.1126/sciadv.1700443 . Google Scholar Crossref Search ADS PubMed WorldCat Norbeck J.H. , Horne R.N., 2016 . Evidence for a transient hydromechanical and frictional faulting response during the 2011 M w 5.6 Prague, Oklahoma earthquake sequence , J. geophys. Res. , 121 , 8688 – 8705 . Google Scholar Crossref Search ADS WorldCat Norbeck J.H. , Rubinstein J.L., 2018 . Hydromechanical earthquake nucleation model forecasts Onset, Peak, and falling rates of induced seismicity in Oklahoma and Kansas , Geophys. Res. Lett. , 45 , 2963 – 2975 . Google Scholar Crossref Search ADS WorldCat Nuttli O. , 1961 . The effect of the earth's surface on the S wave particle motion , Bull. seism. Soc. Am. , 51 , 237 – 246 . OpenURL Placeholder Text WorldCat Peng Z. , Ben-Zion Y., 2004 . Systematic analysis of crustal anisotropy along the Karadere—Düzce branch of the North Anatolian fault , Geophys. J. Int. , 159 , 253 – 274 . Google Scholar Crossref Search ADS WorldCat Provost A.-S. , Houston H., 2001 . Orientation of the stress field surrounding the creeping section of the San Andreas Fault: Evidence for a narrow mechanically weak fault zone , J. geophys. Res. , 106 , 11373 – 11386 . Google Scholar Crossref Search ADS WorldCat Provost A.-S. , Houston H., 2003 . Stress orientations in northern and central California: Evidence for the evolution of frictional strength along the San Andreas plate boundary system , J. geophys. Res. , 108 , doi:10.1029/2001JB001123 . 10.1029/2001JB001123 OpenURL Placeholder Text WorldCat Crossref Queen J.H. , Rizer W.D., 1990 . An integrated study of seismic anisotropy and the natural fracture system at the Conoco Borehole Test Facility, Kay County, Oklahoma , J. geophys. Res. , 95 , (B7) , 11 255 – 11 273 . Google Scholar Crossref Search ADS WorldCat Rawnsley K.D. , Rives T., Petti J.-P., Hencher S.R., Lumsden A.C., 1992 . Joint development in perturbed stress fields near faults , J. Struct. Geol. Mech. Instabil. Rock. Tecton. , 14 , 939 – 951 . OpenURL Placeholder Text WorldCat Ross Z.E. , Hauksson E., Ben-Zion Y., 2017 . Abundant off-fault seismicity and orthogonal structures in the San Jacinto fault zone , Sci. Adv. , 3 , e1601946 , doi:10.1126/sciadv.1601946 . 10.1126/sciadv.1601946 Google Scholar Crossref Search ADS PubMed WorldCat Crossref Ross Z.E. , Rollins C., Cochran E.S., Hauksson E., Avouac J.-P., Ben-Zion Y., 2017 . Aftershocks driven by afterslip and fluid pressure sweeping through a fault-fracture mesh: aftershocks from afterslip and fluids , Geophys. Res. Lett. , 44 , 8260 – 8267 . Google Scholar Crossref Search ADS WorldCat Ross Z.E. , Trugman D.T., Hauksson E., Shearer P.M., 2019 . Searching for hidden earthquakes in Southern California , Science , 364 , 767 – 771 . Google Scholar Crossref Search ADS PubMed WorldCat Savage H.M. , Keranen K.M., Schaff D.P., Dieck C., 2017 . Possible precursory signals in damage zone foreshocks: damage zone earthquakes , Geophys. Res. Lett. , 44 , 5411 – 5417 . Google Scholar Crossref Search ADS WorldCat Savage M.K. , Wessel A., Teanby N.A., Hurst A.W., 2010 . Automatic measurement of shear wave splitting and applications to time varying anisotropy at Mount Ruapehu volcano, New Zealand , J. geophys. Res. , 115 , doi:10.1029/2010JB007722 . 10.1029/2010JB007722 OpenURL Placeholder Text WorldCat Crossref Schoenball M. , Ellsworth W.L., 2017 . A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and Southern Kansas: induced seismicity evolution in Oklahoma , J. geophys. Res. , 122 , 10 189 – 10 206 . Google Scholar Crossref Search ADS WorldCat Schwab D.R. , Bidgoli T.S., Taylor M.H., 2017 . Characterizing the potential for injection-induced fault reactivation through subsurface structural mapping and stress field analysis, Wellington Field, Sumner County, Kansas: injection-induced fault slip, Kansas , J. geophys. Res. , 122 , 10132 – 10154 . Google Scholar Crossref Search ADS WorldCat Schwarz G. , 1978 . Estimating the dimension of a model , Ann. Statist. , 6 , 461 – 464 . Google Scholar Crossref Search ADS WorldCat Segall P. , Lu S., 2015 . Injection-induced seismicity: poroelastic and earthquake nucleation effects , J. geophys. Res. , 120 , 5082 – 5103 . Google Scholar Crossref Search ADS WorldCat Shelly D.R. , Ellsworth W.L., Hill D.P., 2016 . Fluid-faulting evolution in high definition: connecting fault structure and frequency-magnitude variations during the 2014 Long Valley Caldera, California, earthquake swarm , J. geophys. Res. , 121 , 1776 – 1795 . Google Scholar Crossref Search ADS WorldCat Sibson R.H. , 1985 . A note on fault reactivation , J. Struct. Geol. , 7 , 751 – 754 . Google Scholar Crossref Search ADS WorldCat Sibson R.H. , 1990 . Rupture nucleation on unfavourably oriented faults , Bull. seism. Soc. Am. , 80 , 1580 – 1604 . OpenURL Placeholder Text WorldCat Silver P.G. , Chan W.W., 1991 . Shear wave splitting and subcontinental mantle deformation , J. geophys. Res. , 96 , 429 – 416 . Google Scholar Crossref Search ADS WorldCat Skoumal R.J. , Kaven J.O., Walter J.I., 2019 . Characterizing seismogenic fault structures in Oklahoma using a relocated template-matched catalog , Seismol. Res. Lett. , 90 , 1535−1543 . 10.1785/0220190045 OpenURL Placeholder Text WorldCat Crossref Skoumal R.J. , Ries R., Brudzinski M.R., Barbour A.J., Currie B.S., 2018 . Earthquakes induced by hydraulic fracturing are pervasive in Oklahoma , J. geophys. Res. , 123 , 10 918 – 10 935 . Google Scholar Crossref Search ADS WorldCat Snee J.-E.L. , Zoback M.D., 2018 . State of stress in the Permian Basin, Texas and New Mexico: implications for induced seismicity , Leading Edge , 37 , 127 – 134 . Google Scholar Crossref Search ADS WorldCat Sumy D.F. , Cochran E.S., Keranen K.M., Wei M., Abers G.A., 2014 . Observations of static Coulomb stress triggering of the November 2011 M 5.7 Oklahoma earthquake sequence , J. geophys. Res. , 119 , 1904 – 1923 . Google Scholar Crossref Search ADS WorldCat Sumy D.F. , Neighbors C.J., Cochran E.S., Keranen K.M., 2017 . Low stress drops observed for aftershocks of the 2011 M w 5.7 Prague, Oklahoma, earthquake: low stress drops for Prague aftershocks , J. geophys. Res. , 122 , 3813 – 3834 . Google Scholar Crossref Search ADS WorldCat Teanby N. , Kendall J.-M., Jones R.H., Barkved O., 2004 . Stress-induced temporal variations in seismic anisotropy observed in microseismic data , Geophys. J. Int. , 156 , 459 – 466 . Google Scholar Crossref Search ADS WorldCat Trugman D.T. , Dougherty S.L., Cochran E.S., Shearer P.M., 2017 . Source spectral properties of small to moderate earthquakes in Southern Kansas , J. geophys. Res. , 122 , 8021 – 8034 . Google Scholar Crossref Search ADS WorldCat Trugman D.T. , Shearer P.M., 2017 . GrowClust: a hierarchical clustering algorithm for relative earthquake relocation, with application to the Spanish Springs and Sheldon, Nevada, earthquake sequences , Seismol. Res. Lett. , 88 , 379 – 391 . Google Scholar Crossref Search ADS WorldCat U.S. Geological Survey , 2019 . ANSS comprehensive earthquake catalog , Advanced National Seismic System Comprehensive Earthquake Catalog. Accessed from https://earthquake.usgs.gov/earthquakes, Retrieved July 17, 2019 . OpenURL Placeholder Text WorldCat Walsh F.R. , Zoback M.D., 2015 . Oklahoma's recent earthquakes and saltwater disposal , Sci. Adv. , 1 , e1500195 , doi:10.1126/sciadv.1500195 . 10.1126/sciadv.1500195 Google Scholar Crossref Search ADS PubMed WorldCat Crossref Walsh F.R. , Zoback M.D., 2016 . Probabilistic assessment of potential fault slip related to injection-induced earthquakes: application to north-central Oklahoma, USA , Geology , 44 , 991 – 994 . Google Scholar Crossref Search ADS WorldCat Wang C.-Y. , Doan M.-L., Xue L., Barbour A.J., 2018 . Tidal response of groundwater in a leaky aquifer—application to Oklahoma , Water Resour. Res. , 54 , 8019 – 8033 . Google Scholar Crossref Search ADS WorldCat Way H.S.K. , 1983 . Structural Study of the Hunton Lime of the Wilzetta Field T12-13 N, R5E, Lincoln County, Oklahoma, Pertaining to the Exploration for Hydrocarbons. Accessed from https://shareok.org/handle/11244/16448 . Weingarten M. , Ge S., Godt J.W., Bekins B.A., Rubinstein J.L., 2015 . High-rate injection is associated with the increase in U.S. mid-continent seismicity , Science , 348 , 1336 – 1340 . Google Scholar Crossref Search ADS PubMed WorldCat Wessel P. , Smith W.H.F., 1991 . Free software helps map and display data , EOS, Trans. Geophys. Am. Un. , 72 , 441 – 446 . Google Scholar Crossref Search ADS WorldCat Wessel P. , Smith W.H.F., Scharroo R., Luis J., Wobbe F., 2013 . Generic mapping tools: improved version released , EOS, Trans. Geophys. Am. Un. , 94 , 409 – 410 . Google Scholar Crossref Search ADS WorldCat Wibberley C.A.J. , Shimamoto T., 2003 . Internal structure and permeability of major strike-slip fault zones: the Median Tectonic Line in Mie Prefecture, Southwest Japan , J. Struct. Geol. , 25 , 59 – 78 . Google Scholar Crossref Search ADS WorldCat Williams J.A. , 2017 , Geologic, permeability, and fracture characterization of the Arbuckle Group in the Cherokee Platform, Oklahoma , Master's thesis, Emporia State University. Retrieved from http://hdl.handle.net/123456789/3562 . Wu Q. , Chapman M., Chen X., 2018 . Stress-drop variations of induced earthquakes in Oklahoma , Bull. seism. Soc. Am. , 108 , 1107 – 1123 . Google Scholar Crossref Search ADS WorldCat Yeck W.L. , Weingarten M., Benz H.M., McNamara D.E., Bergman E.A., Herrmann R.B., Rubinstein J.L., Earle P.S., 2016 . Far-field pressurization likely caused one of the largest injection induced earthquakes by reactivating a large preexisting basement fault structure , Geophys. Res. Lett. , 43 , 10 198 – 10 207 . Google Scholar Crossref Search ADS WorldCat Zhang X. , Sanderson D.J., 1995 . Anisotropic features of geometry and permeability in fractured rock masses , Eng. Geol. , 40 , 65 – 75 . Google Scholar Crossref Search ADS WorldCat Zinke J.C. , Zoback M.D., 2000 . Structure-related and stress-induced shear-wave velocity anisotropy: observations from microearthquakes near the Calaveras fault in central California , Bull. seism. Soc. Am. , 90 , 1305 – 1312 . Google Scholar Crossref Search ADS WorldCat Zoback M.D. , Kohli A., Das I., Mcclure M.W., 2012 . The importance of slow slip on faults during hydraulic fracturing stimulation of shale gas reservoirs , in Presented at the SPE Americas Unconventional Resources Conference, Society of Petroleum Engineers , doi:10.2118/155476-MS . 10.2118/155476-MS Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Crossref Published by Oxford University Press on behalf of The Royal Astronomical Society 2020. This work is written by (a) US Government employee(s) and is in the public domain in the US.
TI - Activation of optimally and unfavourably oriented faults in a uniform local stress field during the 2011 Prague, Oklahoma, sequence
JF - Geophysical Journal International
DO - 10.1093/gji/ggaa153
DA - 2020-07-01
UR - https://www.deepdyve.com/lp/oxford-university-press/activation-of-optimally-and-unfavourably-oriented-faults-in-a-uniform-i2FNiQ7vDp
SP - 153
VL - 222
IS - 1
DP - DeepDyve
ER -