Revisiting the Timpson Induced Earthquake Sequence: A System of Two Parallel Faults

The 17 May 2012 M4.8 Timpson earthquake is the largest known earthquake in eastern Texas. It is thought to have been induced by wastewater injection from two nearby, high‐volume wells. Its cataloged aftershocks form a NW‐SE trend, which unlike other induced earthquakes sequences is unfavorably oriented for failure in the local stress field. To understand this, we enriched the catalog using PhaseNet, a deep‐learning‐based picker followed by double‐difference relocation with cross‐correlation‐based differential traveltimes. We clustered the aftershocks based on waveform similarity. Most of the seismicity falls into two clusters, which define a detailed fault structure of two parallel subfaults that are more favorably oriented than the overall trend. We inferred from waveform similarity that the sequence initiated on the northern subfault with a M3.9 foreshock and M4.8 mainshock, then extended to the southern subfault with a M4.1 aftershock, and was finally reactivated on the northern subfault with two more M4 events.


Introduction
The 17 May 2012 Timpson earthquake occurred after 5 years of wastewater injection in two nearby disposal wells. It is historically the largest recorded earthquake in east Texas, an area with otherwise low seismic hazard (Frohlich et al., 2014). The triggering mechanism for this earthquake sequence is still poorly understood. Frohlich et al. (2014) located the aftershocks and found that they occurred on a NW-SE striking fault. Lund Snee and Zoback (2016) measured the local stress field from borehole breakouts and showed that the SHmax orientation ranged from N68°E to N80°E. They also calculated the Coulomb stresses on the previously determined NW-SE striking fault and found the fault delineated by aftershocks was quite unfavorably oriented for failure. Fan et al. (2016) conducted geomechanical modeling with poroelasticity and showed that the unfavorably oriented fault could be activated under highly elevated pore pressure; however, their results relied on the strong assumption that no other faults with more favorable orientation were present. Unlike the Timpson earthquake sequence, most of the induced earthquakes in Texas occurred on favorably oriented faults (Frohlich et al., 2011;Gan & Frohlich, 2013;Justinic et al., 2013). Walsh and Zoback (2016) calculated fault slip potential on mapped faults in north-central Oklahoma by considering fault orientation uncertainty. Hennings et al. (2019) applied fault slip potential analysis to the Fort Worth Basin, Texas, to assess risk from induced seismicity in the metropolitan region. Because the assumption that induced seismicity occurs on favorably oriented faults is important to characterizing their hazard, we revisited the Timpson induced earthquake sequence to understand the discrepancy.
We first present catalog reconstruction and earthquake relocation results in section 2. Then we describe the two parallel faults structure revealed by waveform-similarity based event clustering (section 3.1). We associated the five major events in the Timpson sequence using a combination of spatiotemporal correlation and waveform similarity with the two subparallel faults and resolved the sequence initiation and ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. evolution (section 3.2). Finally, we show the two subfaults are better oriented in the local stress field, which supports the importance of accurate fault geometry in hazard assessment (section 3.3). Frohlich et al. (2014) built a catalog for the Timpson sequence and found the aftershocks coincide with a mapped fault. In this study, we enriched the catalog by reanalyzing continuous recordings from 26 May 2012 when the local stations were deployed through the end of 2013. We ran PhaseNet (Zhu & Beroza, 2019), a machine-learning-based phase picker to generate P and S phase picks and associated them using REAL (Zhang et al., 2019), a grid search associator. The associated picks were used as input for VELEST (Kissling et al., 1994) with an initial velocity model determined by Frohlich et al. (2014) and subsequently HypoDD (Waldhauser, 2001;Waldhauser & Ellsworth, 2000). We also added cross-correlation differential traveltime measurements in double difference relocation. The relocated events are shown in Figure 1a.

Relocation of the Timpson Events
To study the detailed fault structure and evolution of the sequence, we applied hierarchical clustering to waveform similarity using Ward linkage as a clustering metric (Akuhara & Mochizuki, 2014). This is a classical method to group the events by minimizing the variance within clusters. Compared with using simple linkage, it better divides the cloud with connecting events between two clusters, which is important in this case as the aftershocks locate close to one another. We filtered the waveforms from 1 to 20 Hz and calculated pairwise cross-correlation coefficients in 1 s windows around both the P and S wave arrivals. By windowing, we avoided having both the P and S wave in the same window at the cross-correlation step. This retains the sensitivity of cross-correlation coefficients to event waveforms rather than to S-P times, or event location. The distance in clustering space between each event pair was defined as one minus the maximum of cross-correlation coefficients for all stations and phases. We took a maximum to keep the correlation coefficient of the waveform with the highest quality because the station coverage and recording quality changes with time. As similar events tend to have high cross-correlation coefficients, they would have short distances in the clustering space by doing one minus the cross-correlation coefficients. We determined the number of clusters by balancing a trade-off between distance within clusters and cluster size. One can always divide a cluster into two and make them more tightly connected inside, but splitting into too many small clusters would render it not useful for fault structure interpretation. In balancing this trade-off, we decided to split the cloud into five clusters (supporting information Figure S1), which resulted in two major clusters that contain most events: Cluster 1 on the north and Cluster 2 on the south (Figure 1a). Events in the two main clusters are marked by dark gray (Cluster 1) and light gray (Cluster 2). Events from the same cluster are closely located in space even if we only provide waveform similarity information in clustering without using traveltimes. The two September 2013 M4 events are located in Cluster 1 while the January 2013 event is near Cluster 2. These larger events occurrence time coincided with high seismicity rate in time as well. As shown in Figure 1b, the January 2013 event occurred when Cluster 2 was most active and the two September 2013 events occurred during the peak activity of Cluster 1. In Figure 1a, the 10 May 2012 foreshock and 17 May 2012 mainshock are plotted with original catalog locations from Frohlich et al. (2014). They suggest these early event locations have large uncertainty due to a lack of local stations. We relocated the 10 May 2012 foreshock and 17 May 2012 mainshock using their traveltime differences at the more distant stations relative to other large events. The relocated early foreshock and mainshock are shown in Figure 1a. Although they seem to locate closer to the northern cluster, considering their location uncertainty, it is still not clear how the sequence initialized only from event locations and occurrence time. Below we use early event waveform similarity with the aftershocks as additional information to infer how this sequence initiated and evolved.

The Case for Two Parallel Faults
We investigated fault structure by first separately relocating Cluster 1 and Cluster 2. Events in each cluster locate in a unique area and approximate a plane, suggesting that they occurred on distinct fault structures. To verify the geometry of the two clusters, we ran HypoDD (Waldhauser, 2001;Waldhauser & Ellsworth, 2000) in SVD mode with all events from both clusters (Figure 2). Compared with the LSQR mode, the SVD mode removes damping and gives unbiased estimate of the location errors. Within uncertainties, the locations clearly delineate two subparallel faults with about a 1 km right-step offset. The local maximum horizontal stress orientation is N68°E to N80°E (Lund Snee & Zoback, 2016). In this stress field, the two faults we delineate are more favorably oriented for slip than the overall trend. The Coulomb stress analysis quantifies this statement in section 3.3 with updated information on the stress orientation. Frohlich et al. (2014) also proposed two groups based on the observation that events S-P time at ETX02 (station marked in Figure 1a) formed two groups. Figure 3b visualizes characteristics of the seismograms in the two clusters. We find systematic differences between the two clusters in the S wave first motion and P coda. Events in the same cluster are similar with one another. We can use this waveform similarity to infer early event locations during the initiation of the sequence and before local stations were deployed.

Initiation and Evolution of the Sequence
The Timpson sequence initiated with a M3.9 foreshock on 10 May 2012 and the M4.8 mainshock followed on 17 May 2012. The first temporary local station started operation on 26 May 2012 and the nearest permanent station, NATX, is located 25 km to the southwest. As a consequence, the locations of early events, including the M3.9 foreshock and the M4.8 mainshock, are uncertain; however, we can infer their locations from their waveform similarity with the well-located aftershocks at NATX, which is the only regional station available at all times. We computed cross-correlation coefficients between early events and all the aftershocks recorded by NATX in the two well-located clusters. The waveform analysis of the M4.8 mainshock is shown in Figure 3a as an example. The events in the northern cluster and the southern cluster are shaded by their waveform correlation with the M4.8 mainshock. We can see from Figure 3a that the waveform of the M4.8 mainshock is more similar to the northern cluster than the southern cluster. A similar pattern can be seen for the M3.9 foreshock and other early events in 2011 as well. This suggests that the Timpson sequence likely initiated on the northern subfault. Together with the aftershock spatial and temporal pattern show in Figure 1, we infer the sequence evolution history as follows. The Timpson sequence initiated on the

Changes of Coulomb Stress With Fault Structure
We calculated the corresponding Coulomb stresses for the two parallel faults versus the single-fault interpretation. To estimate the fault orientation with uncertainties, we took the event relocations with errors and simulated 1,000 set of earthquake locations using Monte Carlo method and used them for fault plane fitting (Browaeys, 2020). The fitted fault plane orientations calculated from the mean value of the 1,000 slopes are plotted with 2σ errors in Figure 4. We first assumed the entire sequence occurred on a single fault with a best-fit fault orientation of N139.8°E, which agrees with N138°E trend reported by Frohlich et al. (2014). Then we estimated the best fit fault plane orientations from the events in the northern cluster and the southern cluster, respectively. We also checked whether the fault planes outlined by the aftershocks matched the  strikes given by focal mechanisms of large events in the sequence. For the four large events with reported moment tensors by the U.S. Geological Survey (USGS) catalog (USGS National Earthquake Information Center, 2016), the reported strikes ranged from 141°to 171°, which is more uncertain than the fault plane orientation. The large uncertainty in the fault strike likely resulted from limited station coverage and velocity model uncertainties. Therefore, we cannot distinguish the two parallel faults scenario from the single-fault scenario simply based on reported focal mechanisms. Next we calculated the stress conditions of the faults for each of the two scenarios. Here we applied the updated stress measurements from Lund Snee and Zoback (2020) shown in Figure 4. In their previous analysis, Lund Snee and Zoback (2016) calculated an average stress orientation of N74°E based on field measurements in two local wells. Their updated measurements (Lund Snee & Zoback, 2020), however, further refine the local stress field. As shown in Figure 4 inset, the local maximum horizontal stress orientation with the new stress data gradually rotates from N85°E in the southwest to N68°E in the northeast. Given the new stress measurements, the local maximum horizontal stress orientation would likely to be around N80°E, which we chose for the Coulomb stress calculation.
We visualize the stress conditions in a Mohr circle diagram ( Figure 4). As we can see from Figure 4, the stresses acting on both faults, yellow for the northern Cluster 1 and blue for the southern Cluster 2, are closer to the failure line compared with the single-fault scenario (red arrow). The two parallel faults are better oriented in the local stress field. By considering the two-fault structure, we require only half of the pore-pressure perturbation to activate the fault as opposed to the single-fault scenario. We also calculated the required pore-pressure on the conjugate fault plane of the M4.8 mainshock shown with the green arrow in Figure 4. Lund Snee and Zoback (2016) suggested this conjugate plane was almost ideally oriented in the observed stress state. However, this conjugate plane is not consistent with observed seismicity. Our results show that the two parallel subfaults are close to failure while matching the seismicity distribution. Our calculated pore-pressure change also agrees with observation of Shirzaei et al. (2016) of little surface uplift, which argues against a large pore pressure increase. According to the geomechanical modeling results of Fan et al. (2016), the two parallel subfaults are within the range where fault slip would occur before hydraulic fracturing. Our findings emphasize that stability of a fault can be very sensitive to its orientation in the stress field. Here, a 10°change in the fault orientation or the local stress field orientation turned an unfavorably oriented fault into one plausibly activated by fluid injection. Therefore, precisely determined fault geometry can be a key factor for assessing fault slip potential.

Conclusions
We revisited the 17 May 2012 Timpson earthquake sequence using machine-learning-based phase detection and uncovered complex fault structure by a combination of precise hypocenter relocation and waveformsimilarity-based clustering. The Timpson earthquakes ruptured two nearly parallel subfaults that are separated by a 1 km right step. This structure was verified by double-difference relocation with cross-correlation differential traveltime measurements. We associated the five major events in the Timpson sequence with its aftershocks based on spatiotemporal correlation and waveform similarity and inferred the sequence evolution pattern as follows. The Timpson earthquakes started on the northern subfault which hosted the M3.9 foreshock and the M4.8 mainshock in May 2012. The southern subfault was then activated with increasing seismicity that culminated in a M4.1 event in January 2013. Later in September 2013, the northern fault was reactivated with a peak in seismicity and another two M4 aftershocks. We investigated how the two-subfault structure affected fault slip potential by calculating resolved Coulomb stress with and without this structure and the corresponding pore-pressure perturbation required in each scenario. The calculation suggests that the two subfaults are closer to failure than the single-fault scenario and thus better oriented in the local stress field. Our findings provide a possible explanation to the discrepancy between stress and aftershock distribution in the Timpson induced sequence. A key implication for hazard assessment is that detailed fault structural knowledge may significantly change the fault slip potential and the hazard level of a fault.

Data Availability Statement
The waveforms can be downloaded from Incorporated Research Institutions for Seismology (IRIS) Data Services (https://ds.iris.edu/ds/nodes/dmc/data/types/waveform-data/). The relocated earthquake catalog is available in the supporting information (Table S1).