Systematic Search for Repeating Earthquakes Along the Haiyuan Fault System in Northeastern Tibet

Repeating earthquakes have been found at many faults around the world. The Haiyuan fault is a major left‐lateral strike‐slip fault along the northeastern boundary of the Tibetan Plateau. Two great earthquakes (1920 Haiyuan and 1927 Gulang) have occurred on and round this fault system, but the section between the ruptures of the two earthquakes, also known as the Tianzhu seismic gap, remains unbroken. Shallow creep has been observed from geodetic data at the Laohushan section of the Haiyuan fault near the eastern end of the seismic gap. However, the driving mechanism and depth extent of shallow creep are not clear. Here we conduct a systematic search for repeating earthquakes in northeastern Tibet based on seismic data recorded by permanent stations in 10 years (2009–2018). Based on waveform cross‐correlations, we find several repeating earthquake clusters at the Laohushan section. This is consistent with the shallow creep inferred from the geodetic data, indicating repeating earthquakes can be driven by nearby aseismic slip. Approximately 300 repeaters were found within clusters of intense seismicity near the rupture zones of the 1927 M8.0 Gulang and 2016 M6.4 Menyuan earthquakes. Relocation of events in the cluster near the Gulang earthquake delineates two possible unmapped faults orthogonal to the Haiyuan fault. In addition, we also identify several repeating earthquakes generated by mining activities with different waveforms and occurrence patterns. Our study suggests that repeating earthquakes around the Haiyuan fault are likely driven by long‐term postseismic relaxation process associated with the 1920 Haiyuan and 1927 Gulang earthquakes.


Introduction
Repeating earthquakes (also known as repeaters) are families of seismic events generated by repeated loading and failure of a single fault patch (e.g., Nadeau et al., 1995;Vidale et al., 1994). Because they are typically driven by aseismic slow slip surrounding them (Beeler et al., 2001), repeaters provide new insight into diverse fault slip behavior at depth, which is usually difficult to characterize from surface observations alone. These include postseismic afterslip, triggered creep and episodic slow slip events, and steady fault creep during the interseismic period (Uchida, 2019;Uchida & Bürgmann, 2019, and references therein). In addition, repeating earthquakes can be used to quantify temporal changes of seismic velocities (e.g., Peng & Ben-Zion, 2006;Poupinet et al., 1984;Rubinstein & Beroza, 2004;Schaff & Beroza, 2004) or seismic scatterers at depth (Niu et al., 2003;Taira et al., 2009).
Repeating earthquakes were first identified in Central California along the Calaveras Fault (Vidale et al., 1994) and the Parkfield section of the San Andreas Fault (Nadeau et al., 1995). Since then, they have been found on major plate boundary faults elsewhere around the world, such as the Japan trench (e.g., Igarashi et al., 2003;Uchida et al., 2003), Taiwan , Tonga (e.g., Yu, 2013), Costa Rica , and Turkey (Peng & Ben-Zion, 2005. There are also increasing reports on repeating earthquakes in intraplate settings. For example, repeaters have been found in seismic zones in the central United States (Bisrat et al., 2012). Based on cross-correlating of regional seismic waveforms, Li et al. (2007) found repeaters in eastern China along the aftershock zone of the 1976 M7.6 Tangshan earthquake, and Li et al. (2011) Zhang et al. (1988) obtained a slip rate of 8 ± 2 mm/year along the eastern section of the HYF. Based on offset geomorphic features and age constraints, Li et al. (2009) inferred a slip rate of 4.5 ± 1.0 mm/year on the same section. However, Lasserre et al. (1999) estimated a high slip rate (12 ± 4 mm/year) in the middle section from the geomorphic mapping and dating of offset terrace risers. This was recently reevaluated and updated to be 5-9 mm/year based on the new dating, high-resolution airborne Light Detection and Ranging topography, and both upper-and lower-terrace reconstructions of geomorphic offsets (Yao et al., 2019).
At the east end of the middle Laohushan section, shallow creep has been observed from geodetic data and is estimated to be 5 ± 1 mm/year (Jolivet et al., 2012(Jolivet et al., , 2013, close to or slightly smaller than the geologic rate. However, creep in this section is poorly understood as compared to other known creeping faults around the world (e.g., Harris, 2017). In addition, the depth extent of the creep is not well constrained. It is still not clear whether creep on this section of the HYF is a transient phenomenon following the 1920 Haiyuan mainshock or reflects a long-term slip behavior (Chen et al., 2018).
In order to better understand the seismicity pattern and fault slip behaviors along the HYF, we conduct a systematic search for repeating earthquakes in this region based on 10 years of microseismic data. Specifically, we identify repeating earthquake pairs using waveform cross-correlations (CCs), and then we group the pair into clusters. Finally, we use these repeating clusters to better understand the aseismic process and slip rates along the HYF and other faults in NE Tibet.

Data and Methods
Within the following geographic boundaries (longitudes between 96°and 107°E, latitudes between 36°and 41°N), there are 25,888 M ≥ 1 events from 2009 to 2018 based on the regional catalog archived by the China Earthquake Networks Center (CENC) (Figure 1). The earthquakes are mostly distributed along the HYF and other major faults/boundaries in this region, such as near the eastern end of the Altyn Tagh fault and the western boundary of the Ordos basin. We obtain raw seismogram data (100 Hz sampling rate) recorded by 51 permanent stations from the Data Management Center of the China National Seismic Network (Zheng et al., 2010). These include 44 stations deployed before 2009 and seven stations deployed at the end of 2014 (supporting information Figure S1).
Hypocenter location and waveform similarity are two main methods to identify repeating earthquakes (Uchida, 2019). Here we use the later one because it is the most frequently used method and we can relocate the events afterward. Our analysis procedure includes the following six steps. First, we compute the P and S arrival time using the Taup program (Crotwell et al., 1999) according to a local velocity model (Deng et al., 2018; Table S1) modified from IASPEI91 model (Kennett & Engdahl, 1991). Because there are some handpicked P and S phases in the catalog, we use them when available and use the computed arrival times for the rest traces. We then apply a 2-16 Hz band-pass filter on the data, which is the predominant frequency range for microearthquakes in the magnitude range of 1-4 (Uchida, 2019). Next, we perform a quality control on the raw data, based on the signal-noise ratio (SNR). We use a 20 s window (5 s before the P wave) as the signal window and the same 20 s window 25 s before the P wave as the noise window. We choose the data with the SNR higher than 4 and remove the rest of the seismograms with low SNRs.  Zheng et al. (2017). Active faults are adopted from Tapponnier et al. (2001). The "Tianzhu seismic gap" (Gaudemer et al., 1995) is highlighted in yellow; the surface ruptures of the 1920 Haiyuan and 1927 Gulang earthquakes are highlighted with red lines; stars denote the epicenters of the Haiyuan, Gulang, and Menyuan earthquakes.
We then compute CCs among all M ≥ 1 events beneath all stations and identify event pairs with locations less than 200 km apart with CC value above 0.80. Because the time windows are usually set to contain both P and S phases to ensure the same P-S time and thus the same hypocentral distance (Uchida et al., 2003), here we choose 2 s before and 18 s after the P arrivals. The predicted P arrival may be not accurate due to the complicated velocity structure and inaccurate epicentral location, so we allow a maximum time shift of 5 s to obtain the highest CC between each waveform pair. Next, we group repeating pairs into clusters using an Equivalency Class algorithm (Peng & Ben-Zion, 2005;Press et al., 1986). When two pairs have a common event and meet a certain threshold (median CC value ≥ 0.9 with number of at least two stations), we group them into the same cluster. Because the final number of repeating clusters depending strongly on the threshold values (e.g., Peng & Ben-Zion, 2005), we also group repeating clusters based on a threshold CC value from 0.91 to 0.95. As will be shown in the following sections, although the total number of repeating clusters are different, their spatial distributions and the resulting slip rate are similar. Finally, we relocate the repeaters in each cluster using the Growclust program (Trugman & Shearer, 2017) based on the differential arrival times from waveform CCs.
Based on the identified repeaters, we can estimate their cumulative slip according to their local magnitudes M L . The individual slip d can be estimated by assuming a standard crack model where r is the rupture size, M 0 is the seismic moment, and μ is shear modulus that is set to be 3 × 10 10 N/m 2 . The rupture size r can be obtained from (Kanamori & Anderson, 1975) where Δσ is the static stress drop. The seismic moment M 0 can be estimated from the local magnitude M L with the following equation (Abercrombie, 1996) Finally, the average slip rate is computed from dividing the cumulative slip with the total duration of the repeating cluster (Li et al., 2011). While the local magnitudes M L can be obtained from the CENC catalog, the stress drops Δσ for these events are not available. Hence, we use nominal stress drop values of 1, 5, and 10 MPa (Li et al., 2011).

Results
Based on the median CC threshold of 0.9, we identify 929 clusters (~2,500 events) in NE Tibet with at least two repeaters, which accounts for 9.7% of earthquakes in this region. Figure 2 shows the spatial distribution of all repeating clusters, together with the background seismicity. Even though repeaters are widely distributed, clusters with more than eight events are only found at certain regions along and outside of the HYF. In the following sections, we discuss them in more details. If we change the threshold to be 0.92 and 0.95, the resulting numbers of clusters are 485 (1,178 events) and 72 (154 events), corresponding to 4.6% and 0.6% of seismicity in this region (Table 1), respectively. The corresponding spatial distributions are shown in Figure S2, similar to those shown in Figure 2. In the subsequent sections, unless otherwise notified, the results are based on the threshold of 0.9. Figure 3 shows the average magnitude within clusters ranges mostly from 1 to 2.5, and the standard deviation of magnitude is mostly less than 0.3. The average distance between event pairs within clusters is less than 10 km, and the standard deviation is mostly less than 0.5 km. Figure S3 shows the histograms of recurrence intervals for all repeating clusters. Only a few pairs have very short recurrence intervals, likely representing a burst or short-term mainshock-aftershock pair (e.g., Lengliné & Marsan, 2009), the majority of them have intervals on the range of a few tens to hundreds of days, likely representing those driven by tectonic creep (e.g., Uchida & Bürgmann, 2019).

10.1029/2020JB019583
Journal of Geophysical Research: Solid Earth

The Creeping Section of the HYF (Middle Laohushan Section) (Region 1)
We find the 31 repeating earthquakes within 12 clusters along the LHSF, which is in the central section of the HYF (Figure 4). From the map view, the repeaters are mostly located around 103.7°E, to the west side of the peak creep region as reported by Jolivet et al. (2013). The 87% and 94% of repeaters have the depth shallower than 10 km before and after relocation, respectively. Relative location error is less than 0.8 km in the horizontal direction, less than 1.4 km in the vertical direction ( Figure S4).       (Li et al., 2016). Based on the matched filter detection and relocation of early aftershocks, Liu et al. (2019) inferred that the 2016 Menyuan mainshock occurred on a steeply dipping secondary fault rather than the major LLLF. This interpretation is consistent with the geodetic and geological observations (Li et al., 2016). Liu et al. (2019) also found 26 repeating clusters (~172 events) in the aftershock zone of the Menyuan mainshock.
From 10 years of seismic data, we identify 81 events within 37 clusters in this region ( Figure 6). Most of them occurred following the 2016 Ms6.4 mainshock. A subtle increase of repeating events is also observed after the 2016 M4.7 and 2017 M4.1 events (likely aftershocks of the Ms6.4 mainshock) but not following the 2013 M5.3 event. Because repeating aftershocks are mostly driven by afterslip Schaff & Beroza, 2004), their recurrence times increase following the Menyuan mainshock . Hence, we do not compute their cumulative slip and average slip rate in this region.

Gulang Seismic Zone, North of the HYF (Region 3)
We found 210 repeating earthquakes (75 clusters) that occurred in the Gulang seismic zone near the epicenter of the 1927 Gulang earthquake (Figure 7). From the map view, these repeating earthquakes are located between the strike-slip HYF in the south and the south-dipping Huangcheng-Shuangta fault in the north. The repeating earthquakes are within the clusters of intensive background seismicity that delineate two nearly N-S striking features. But there are no mapped faults on the surface. These repeating earthquakes occurred every year from 2009 to 2018, and the occurrence time is almost uniformly distributed during local  In comparison, some repeating clusters are not related with any moderate-sized earthquakes. For example, four or five repeating events are in the four clusters C1-C4 (Table S3), and their recurrence intervals are mostly within days. Since these sequences do not follow a clear moderate-sized mainshock (earthquakes with M ≥ 4; Figure 7), they may be considered as earthquake swarms (Vidale & Shearer, 2006), but we could not rule out the triggering by nearby fault patches because there are several other events occurred at the similar period of time.

Mining-Related Repeaters
In the CENC catalog, there are 332 marked mining explosions at the boundary between Neimenggu and Ningxia Provinces (39.0-39.2°N, 105.95-106.15°E; marked as Region 5 in Figure 2). Because mining explosions almost occur in the same position, the corresponding seismic waveforms have very high similarities and can be detected as repeaters in this study. In our study region we find 136 possible repeating events within 25 clusters that are likely mining related.
We use the following evidence to confirm that these events are likely related to mining activities. First, some events listed in the repeating clusters are marked as mining explosion in the CENC catalog. Second, the local times of these events have peaks in the afternoon instead of a uniform distribution like Regions 1 and 3 (Figure 8), which is consistent with the expected mining explosion schedule (e.g., Ruan et al., 2017). Third, we can identify possible mining-related images with Google Map in this region ( Figure S5). Finally, the event magnitudes are tightly clustered (Figure 9), and there are no M > 4 events in this region.

Journal of Geophysical Research: Solid Earth
In addition to this region, we also find three other regions (Regions 6-8 in Figure 2) with possible mining explosions. Compared with clear P and S arrivals for regular earthquakes, explosion-generated waveforms have different phases (Figures 5 and 9). We note that the first few cycles of the P waves are not quite the same among these regions (Figure 9), which may indicate different styles of explosions (e.g., delayed vs. single fire) (Stump et al., 2002), as well as possible structural difference.
There are also some repeating events in other places, such as Madongshan that is marked as Region 4 in Figure 2. In other regions, if there are no M > 4 events and no marked surface fault, we generally consider them as possible mining explosions.

The Characteristics of Repeating Earthquakes at the Laohushan Section
In this study, we identified 12 clusters with 31 repeating events along and near the Laohushan section of the HYF, where surface creep has been identified (Jolivet et al., 2013). As mentioned before, it is still not clear what is the cause of creep in this region. There is some evidence for the surface-breaking rupture of paleoearthquakes in this section. In addition, the cumulative slip offset in this region is similar to noncreeping brittle faults (Chen et al., 2018). Based on these observations, Chen et al. (2018) inferred that either the fault is capable of switching between creeping and brittle faulting over time or the fault is partially creeping. In the second possibility, the creep on the Laohushan section is a shallow phenomenon, with the fault remaining locked at depth, similar to the Ismetpasa segment of the North Anatolian fault in Turkey (Kaneko et al., 2013;Karabacak et al., 2011;Ozener et al., 2010) and the Hayward fault in northern California (Schmidt et al., 2005;Simpson et al., 2001).
Based on an assumed stress drop of 1, 5, and 10 MPa (Li et al., 2011), we computed the cumulative slip and slip rate of the repeating earthquakes in this region. We found that the largest slip rate is 6.6 mm/year and the average slip rate is 2.25 ± 2.24 mm/year. Although with large uncertainties, our estimated slip rate is less than the shallow creep rate of~5 mm/year from InSAR observations (Jolivet et al., 2013) and geological rate of 5-9 mm/year (Yao et al., 2019).
As mentioned before, the final numbers (and percentage) of repeating earthquakes depend strongly on the choice of CC threshold. If we increase the CC threshold from 0.90 to 0.95 with a 0.01 step, the total numbers of repeaters are 2,502, 1,809, 1,178, 712, 386, and 154 (Table 1). The corresponding percentage range is from 0.6% to 9.7%. The latter is roughly consistent with the finding that~10% seismic events around China are repeating events (Schaff & Richards, 2004). However, it is worth noting that Schaff and Richards (2004) used primarily Lg waves recorded at regional distances, and their CC threshold was 0.8. Hence, it is possible that at least some of them were burst-type or short-term mainshock-aftershock type (Lengliné & Marsan, 2009), rather than driven by sustained fault creep. In our study, we only found a very small percentage of event pairs with very short recurrence intervals (e.g., within 1 day; Figure S3). Hence, we believe that most repeaters identified in this study reflect long-term creep around the locked fault patches. Finally, although the total numbers of repeaters are quite different with different CC threshold, the resulting average slip rates at the Laohushan section are around 2 mm/year (Table 1), again indicating that they are less sensitive to the threshold parameters.
Our estimation of slip rate based on repeating clusters also has several limitations. First, we only take the M ≥ 1 events into account. Hence, seismic slip released by events with smaller magnitudes is not included. Second, the M ≥ 1 events catalog may be incomplete, which could result in a biased estimation of slip rate. Finally, it is possible that some asperities creep at interseismic time periods (Beeler et al., 2001). Hence, what we estimate here can be considered as the lower bound of the actual slip rate, although we do not expect to observe several folds increase in the slip rate.
After the relocation, the repeaters are mostly located between 4 and 8 km depth (Figure 10), below the inferred depth of shallow creep (Jolivet et al., 2013). In addition, these repeaters occurred to the west of the peak creep observed from InSAR observation (Jolivet et al., 2013). Our observation suggests that the geodetically observed creep is constrained at shallow depth, while faults at deeper depth remain locked. Hence, these repeating earthquakes likely mark the boundary between the creep and locked region, similar to the Parkfield section of the San Andreas Fault (Lengliné & Marsan, 2009;Nadeau & McEvilly, 1999,

Journal of Geophysical Research: Solid Earth
the Morgan Hill section of the Calaveras Fault Rubin, 2002;Schaff et al., 2002), and the Hayward Fault (Bürgmann et al., 2000;Shirzaei et al., 2013). This is also compatible with the observation of earthquake swarms driven by shallow aseismic slip in Salton Trough, California (Lohman & McGuire, 2007), and anticorrections between afterslip and aftershocks (including repeating earthquakes) following the 2012 M7.6 Nicoya, Costa Rica, earthquake Yao et al., 2017). These studies suggest that while microseismicity (including earthquake swarms and repeating earthquakes) can be driven by nearby aseismic slip, they likely occur in slightly different regions, indicating varying frictional behavior along dip and strike directions.
Considering that the creeping Laohushan section on the HYF is located immediately to the west of 230 km surface rupture of the 1920 Haiyuan earthquake (Liu-Zeng et al., 2015), Chen et al. (2018) speculated that long-term postseismic deformation following the 1920 Haiyuan mainshock might have a heightened effect on the recently observed creep rate. The creeping section is also located at the eastern end of the over-

10.1029/2020JB019583
Journal of Geophysical Research: Solid Earth 200 km long "quiet" seismic gap, known as the Tianzhu gap (Gaudemer et al., 1995). Hence, the Laohushan section is somewhat similar to the creeping section of the San Andreas Fault from Parkfield to Hollister in Central California, which was sandwiched between the 1857 Fort Tejon and 1906 San Francisco earthquakes. However, the creeping section of the HYF is only 35 km (Jolivet et al., 2012), much shorter than the~150 km length for the San Andreas Fault. While the creeping section is generally considered as "barrier" to seismic ruptures, Noda and Lapusta (2013) demonstrated that dynamic earthquake ruptures (combined with coseismic weakening) can break through long portions of creeping faults, indicating the possibility a total rupture on long strike-slip faults such as the San Andreas or HYFs.

The Implication of Intense Seismicity at the Gulang Seismic Zone
Another region (Region 3 in Figure 2) with intense background seismicity and repeating earthquakes is located near the epicenter of the 1927 Gulang earthquake, north of the HYF (Figure 7). Yang (2017) performed repeating earthquake detections in this region and found similar patterns of repeating clusters. While there are no corresponding faults mapped on the surface, relocated seismicity is concentrated along two linear zones trending NNW to nearly NS, indicating two possible hidden faults at depth (Figure 7). Such interpretation is consistent with the available focal mechanisms of moderate-sized events in this region, which are predominately NNW-SEE trending right-lateral strike slip events. Wang (2018) named the intensive seismic activity in this region as "Gulang seismic window," which are stress sensitive regions in the aftershock zones following large earthquakes. Because this region is spatially close to the rupture zone of the 1927 M~8 Gulang earthquake, it is possible that they are extended aftershocks of the Gulang earthquake. This interpretation is consistent with general observations of long-tailed aftershock activity, especially in intraplate regions around the world (Ebel et al., 2000;Stein & Liu, 2009). Even at plate boundary regions, large earthquakes can potentially affect the deformation pattern and cycles of small-to moderate-sized earthquakes at nearby distances. For example, Ben-Zion et al. (1993) employed 3-D finite-element modeling to infer that the M6-type Parkfield earthquakes are driven by time-dependent loading from the 1857 M8 Fort Tejon earthquake. Hence, we argue that intensive earthquake swarms and repeating earthquakes in this region are possibly driven by long-term relaxation process induced by the 1927 M~8 Gulang earthquake.
Based on seismic velocity inversions, Deng et al. (2018) found that this region is mechanically weaker as compared with the North China Craton in the north and the central Qilian in the south. This is comparable with the earthquake swarms at Belo Jardim, NE Brazil (Lopes et al., 2010), and the Salton Trough South California (Lohman & McGuire, 2007), which are likely occurred at the preexisting weak zones.

Conclusions
Based on waveform CCs and relocations, we systematically search for repeating earthquakes in NE Tibet with 10 years of seismic data. The repeating earthquakes are found in certain regions. Laohushan section of HYF has several repeating clusters, whose epicenters are located to the west of the peak creep region indicated by InSAR observations. The slip rate estimated from the repeating earthquakes is slightly smaller than the geodetic and geological observation. However, our estimation is based on assumed constant stress drops and has several limitations and, hence, could be considered as a lower bound for the slip rate. In addition, the relocated microearthquakes are mostly located deeper than 4 km, which is the depth extent of geodetically observed shallow creep. Hence, we infer that these repeating earthquakes are located at the transition between the shallow creeping and deep locked regions. We also find many repeating earthquakes following the 2016 M6.4 Menyuan mainshock. The intense seismicity at the Gulang seismic zone is aligned along two linear zones, indicating two possible hidden faults. This region may be intrinsically weak, and hence, these repeating events could be driven by the long-term relaxation process of the 1927 Gulang earthquake. In comparison, repeating events related to mining explosions in this region have different waveforms from natural earthquakes and can be readily identified based on several diagnostic features.

Data Availability Statement
Waveform data for this study are provided by Data Management Center of China National Seismic Network at Institute of Geophysics' China Earthquake Administration (doi:10.11998/SeisDmc/SN, http://www.

10.1029/2020JB019583
Journal of Geophysical Research: Solid Earth seisdmc.ac.cn). Our final repeating catalogs obtained in this research are included in its supporting information files and through an online website (https://figshare.com/articles/Dataset_for_Deng_s_2020_JGR_ paper/12470057).