# Complex rupture dynamics of the extremely shallow August 2020 M5.1 Sparta, North Carolina earthquake

Miguel Neves, Lindsay Chuang, Wei Li, Zhigang Peng, Paula Figueiredo, Sidao Ni

► To cite this version:

Miguel Neves, Lindsay Chuang, Wei Li, Zhigang Peng, Paula Figueiredo, et al. Complex rupture dynamics of the extremely shallow August 2020 M5.1 Sparta, North Carolina earthquake. Communications Earth & Environment, 2024, 5 (1), pp.163. 10.1038/s43247-024-01316-8. hal-04921754

HAL Id: hal-04921754

https://hal.science/hal-04921754v1

Submitted on 31 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Complex rupture dynamics of the extremely shallow August 2020 M5.1 Sparta, North Carolina earthquake

<figure><img src="image_3.png" /></figure>

Check for updates

Miguel Neves<fnref n="1" />,⁵✉, Lindsay Y. Chuang<fnref n="1" />, Wei Li², Zhigang Peng<fnref n="1" />, Paula M. Figueiredo³ & Sidao Ni⁴

On August 9, 2020, an $M_w$ 5.1 earthquake ruptured the uppermost crust near the town of Sparta, North Carolina, creating the first co-seismic faulting surface rupture documented in the Eastern United States. Combining deep learning and matched filter earthquake detection, with differential-travel times relocation, we obtain a catalog of 1761 earthquakes, about 5.8 times the number of events listed in the standard USGS/NEIC catalog. The relocated seismicity revealed a complex fault structure with distinct planar alignments, supported by a moment tensor inversion with significant non-double-couple component. The Sparta mainshock with a centroid depth of 1.3 km is interpreted to have nucleated near the intersection of two main fault strands. The mainshock likely ruptured a blind strike-slip fault and a reverse fault associated with the identified surface rupture, both possibly part of a flower structure-like diffuse fault zone. Our observations highlight a complex behavior of extremely shallow earthquakes in stable continental regions.

On August 9, 2020, an $M_w$ 5.1 earthquake ruptured the uppermost crust near Sparta, North Carolina, generating a WNW-trending surface rupture (Fig. 1)<fnref n="1" />. This was the largest recorded earthquake in North Carolina in the past 100 years, and likely the first reported surface-rupturing event in the Eastern United States<fnref n="1" />,² and the second reported in Eastern North America³.

The mainshock occurred in a stable continental region characterized by low-strain rates resulting in episodic seismic activity, with widespread clusters (Fig. 1c) and infrequent moderate to large earthquakes with after-shock sequences that can be sustained for long time periods⁴,⁵. It ruptured in the Blue Ridge Province, a part of the larger Southern Appalachian orogen, characterized by prevailing NE-trending structures⁶. The WNW-ESE co-seismic rupture trend is anomalous, however it is consistent with WNW to E-W lineaments recognized to cross the prevailing structural fabric along the Appalachians⁷. However, their genesis and development are not entirely understood and constrained in time. Estimates of the crustal stress field for the region show that the earthquake area is located in a transition zone between strike-slip and reverse faulting, with a maximum horizontal principal stress ($S_{Hmax}$) oriented about NE-SW to ENE-WSW⁸,⁹.

Previously identified regional seismicity is located mostly in the Valley and Ridge Province that borders the Blue Ridge Province where the Sparta mainshock ruptured. Major clusters in the Valley and Ridge Province

include the Giles County Seismic Zone in Virginia¹⁰ to the north, and the Eastern Tennessee Seismic Zone¹¹ to the west. Seismicity in these clusters shows predominantly N to NE striking right-lateral strike-slip motion and also E-W striking left-lateral motion¹⁰,¹².

However, in the Sparta area no earthquakes were reported prior to the 2020 sequence in the National Earthquake Information Center (NEIC) earthquake catalog, part of the United States Geological Survey (USGS). The closest reported earthquake to the area is a magnitude 1 earthquake detected on October 8, 1989 and located about 16 km SW from the 2020 M5.1 epicenter. The $M_w$ 5.1 mainshock was preceded by at least 8 foreshocks reported by the NEIC in the 25 hr before the rupture<fnref n="1" />,². Figueiredo et al.<fnref n="1" /> documented the surface rupture associated with the M5.1 earthquake along the newly identified Little River fault. The rupture shows mostly reverse slip in a preexisting planar structure with a WNW trend and SW dip<fnref n="1" />. In contrast, Wicks and Chiu³, analyzed Interferometric Synthetic Aperture Radar (InSAR) data in the epicenter area and identified a shallow normal faulting event associated with the surface rupture and deeper reverse deformation possibly aseismic. Pollitz¹³ reanalyzed both geodetic and regional seismic data and concluded that no normal-faulting component is needed to explain observations and found evidence of thrust and left-lateral strike-slip motion.

<footnote n="1">¹School of Earth and Atmospheric Sciences, Georgia Institute of Technology, Atlanta, GA, USA. ²School of Earth and Space Sciences, University of Science and Technology of China, Hefei, China. ³Department of Marine, Earth and Atmospheric Sciences, North Carolina State University, Raleigh, NC, USA. ⁴State Key Laboratory of Geodesy and Earth's Dynamics, Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan, China. ⁵Present address: Université Côte d'Azur, IRD, CNRS, OCA, Geoazur, Valbonne, France. ✉e-mail: neves@geoazur.unice.fr</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Study area. a Map of mainshock location and stations used in the study. Blue triangles correspond to seismic stations working at the time of the mainshock rupture and red stations to the temporary stations deployed in the days after the mainshock. Mainshock location from the NEIC is denoted by the yellow star.</figcaption></figure>

<caption>b Zoom around the mainshock location and the town of Sparta (orange square).</caption>

<caption>Black dots represent the identified foreshocks in the NEIC catalog. Black lines show major US and state roads in the area. c Regional map of the Eastern United States showing Sparta mainshock location, all earthquakes listed in the NEIC earthquake catalog and historical earthquakes in the region.</caption>

At the time of the mainshock, the closest seismic station was located at more than 50 km from the epicenter, however following the mainshock, 5 temporary seismic stations were deployed by the Center for Earthquake Research and Information (CERI) of the University of Memphis in the epicenter area (Fig. 1)¹⁴. Taking advantage of this uncommon station coverage for a seismic sequence in the Eastern US, we use state-of-the-art earthquake detection and relocation techniques, to image the sequence in detail and understand the geometry of the local faults at play and the source processes of intraplate earthquakes. We show that the Sparta mainshock nucleated near the junction between a sub-vertical predominantly strike-slip fault and an SW-dipping predominantly reverse fault, and possibly propagated along these two structures.

## Results

### New earthquake catalog

Using our detection and location framework combining deep learning earthquake detection, matched-filter technique (MFT) detection, and differential travel-time relocation (Fig. 2), we improved the earthquake catalog of the sequence. We performed the detection from June 1 to November 18, 2020, about two months prior to the mainshock and two and a half months after, when most of the temporary stations were recovered. Our improved earthquake catalog includes 1761 earthquakes with high-precision locations (see Supplementary Figs. 1–7), about 5.8 times the number of earthquakes in the initial catalog by the NEIC with 301 events (Fig. 3). Moreover, our catalog includes 18 events before the mainshock, which initiated 47 hours prior to the mainshock. Additional detected signals prior to the first foreshock were either relocated outside the earthquake area or did not pass the quality thresholds that we imposed. In terms of the magnitude-frequency distribution (Supplementary Figs. 6–8), our improved catalog has an estimated magnitude of completeness of −1.5 compared to the magnitude of completeness of −0.5 of the NEIC catalog. The largest aftershock in our catalog has a magnitude of 2.87 lower than the expected 3.9 according to the Bath's law¹⁵, and the frequency magnitude distribution further shows a low

number of detections in the M 0.5–3 interval. Similar observations of lower-than-expected largest aftershock magnitude and a low number of aftershock detections have also been reported in a moderate shallow rupturing event in southeast France, the 2019 Mₗ 4.9 Le Teil earthquake¹⁶.

Unlike the original NEIC catalog where no clear seismicity clustering or linear structures can be easily identified (Fig. 3a), our catalog shows two clear lineaments, which are interpreted to be associated with two faults that intersect at a low angle, being that intersection point the locus of the relocated mainshock (Fig. 4a). To obtain a more detailed geometry of the structures highlighted by the seismicity we apply a density based clustering algorithm¹⁷ (Fig. 4a and Supplementary Fig. 9). With this clustering approach, we are able to identify five clusters (see Methods). Three of the clusters form the two distinct lineaments that meet near the epicenter. The seismicity in the lineament eastward of the epicenter (cluster shown on AA' line on Fig. 4 and AA' cross-section in blue on Fig. 5) concentrates at shallower depths between 0.5 and 1.5 km (Fig. 5a) with a spatial distribution roughly following the co-seismic surface faulting trend¹. At a westward lineament (shown on BB' line on Fig. 4 and cross-section on Fig. 5) we identify two clusters, one westernmost shallower spreading from the near-surface to about 1 km deep and a second deeper and closer to the mainshock hypocenter at 2 to 3 km deep. Both clusters strike similarly and show a near vertical spreading almost aligned, suggesting that the clusters are possibly part of the same structure, although between them there is an area with a lower density of earthquakes. This lower density of earthquakes at depths 0.5 to 1 km can be due to the limited time period analyzed or our limited detections in the early aftershock period. But we also note that these depths correspond to the same depth of the high density of earthquakes in the main cluster in the eastward lineament. Additionally, we identify a fourth and fifth minor clusters striking sub-parallel to the eastern cluster and north of the two main lineaments. Most of the largest earthquakes (M 1.5–3.3), including the mainshock and foreshocks in our catalog, are not included in any of the retrieved clusters. These are mostly located in the epicentral area between the identified eastern and western clusters, an area showing a diffuse group
of earthquakes that results in inconsistent density values and is therefore not selected by the clustering algorithm.

We produce a simple fault model by fitting planes to each of the clusters (Fig. 4b). The eastern main cluster is best explained by a structure with a N118.7° ± 1.4' strike and 45.4° ± 2.1' SW dip, with a reduced chi-squared of 16.2. The clusters west of the epicenter are best explained by near vertical structures with an 79.9° ± 0.1' NE best-fitting dip and a 142.6° ± 0.1' strike for the westernmost and shallower cluster and an 88.6° ± 0.1' NE dip and a 138.5° ± 0.1' for the deeper cluster. The smaller fourth and fifth clusters are best explained by structures of similar strike of N114.8° ± 5.2' and N114.6° ± 63.45' and dips of 59.2° ± 9.89' SW and 14.64° ± 35.49' SW, respectively. The fits to these clusters show a reduced chi-squared of 43.4 for the fourth cluster and 65.34 in the fifth cluster. We highlight the poor fit of

<figure><img src="image_5.png" /></figure>

**Fig. 2** | The workflow followed to compile the earthquake catalog. We employed a multistep approach to compile the earthquake catalog (see Methods). First, we obtained a deep learning detection catalog using the Earthquake Transformer that was merged with the NEIC catalog. This combined catalog was then used for matched filter detection (MFT), followed by relative relocation using differential travel times and relative magnitude estimation.

the fifth cluster that is due to the low number of events associated with it, leading to high uncertainties and a non-robust result.

## Focal mechanisms

We perform a moment tensor inversion of the mainshock and find that it has an oblique-thrust focal mechanism with a moment magnitude $M_w$ of 5.14 and best focal planes 350°/70°/115° and 116°/31°/40° (strike/dip/rake) and an optimal centroid depth of 1.3 km. These results are somewhat consistent with a plane with a strike of 118.7° and a dip of 45.4° obtained for the shallow main eastern cluster. The moment tensor inversion also returns a significant non-double-couple component of 26%, indicating that the mainshock rupture likely occurred along a fault plane with complex geometries. Regarding the centroid depth, we obtain a shallower depth than the relocation hypocentral depth of 3 km. This suggests the rupture initiated in a deeper section and progressed to the surface, where it released most of the energy at a shallower structure.

Analyzing the first motion polarities together with the S and P phase amplitude ratios of smaller earthquakes in the catalog, we obtain two types of focal mechanisms: mostly strike-slip focal mechanisms, and mostly reverse focal mechanisms (Fig. 6c). Most of these events are not located within any of the clusters we retrieved, but are instead located in the diffuse area where the main lineaments meet. Still, we can associate one oblique-reverse focal mechanism with the eastern cluster (Fig. 6c, blue focal mechanism), similar to the obtained mainshock moment tensor inversion. We also associate a strike-slip event with the deeper western cluster (Fig. 6c, pink focal mechanism). Considering these observations and the inferred planes' geometry, we refer to the eastern cluster as a mostly reverse cluster, mostly strike-slip cluster 1 to the shallow western cluster, and mostly strike-slip cluster 2 to the deeper western cluster from this point on. We are also able to constrain a focal mechanism showing an oblique-reverse motion associated with the fourth cluster. It is also intriguing that in our first motion polarities analysis of the mainshock, we invert a strike-slip focal mechanism with a preferred plane of 113°/88°/1°. Altogether, from our polarity analysis, we find no seismic evidence of normal faulting suggested by Wicks and Chiu² to be responsible for the surface rupture. These results are consistent with the crustal stress estimates that show a transition zone of strike-slip to reverse faulting³.

## Stress drop

Released shear stress estimates (i.e., stress drop) have long been an important source parameter, due to its impact on the high-frequency ground motions¹⁸,¹⁹. By using the smaller earthquakes as Empirical Green's Functions (EGF)²⁰,²¹, we calculate spectral ratios and determine corner

<figure><img src="image_6.png" /></figure>

**Fig. 3** | Comparison of the earthquake catalog for the different steps of our detection and location framework. a Earthquakes listed in the NEIC catalog. b NEIC earthquake catalog combined with detections performed using the

EQTransformer that were used as templates for a matched filter detection. c Final relocated earthquake catalog. Color denotes depth.
<figure><img src="image_7.png" /><figcaption>**Fig. 4** | Detected earthquake clusters using an unsupervised approach. **a** Map view of the identified clusters, with the mainshock denoted by the yellow diamond. **b** 3D plot of the best fitting planes to each of the identified three clusters shown in (a). **A**</figcaption></figure>

fully interactive 3D plot with the fault planes and the scattered seismicity is included in Supplementary Figs. 10 and 11.

frequencies for the Sparta mainshock, the largest foreshock and two of the largest aftershocks (Supplementary Table 1). We attempted to extend the spectral ratio analysis to smaller events, however, we could not obtain stable spectral ratios for them, most likely due to the large source-station distances in our study which led to low signal-to-noise ratios (SNR). Our results show a corner frequency for the mainshock of 1.07 Hz (Fig. 7a). This results in stress drop estimates of 16.3 MPa considering a Brune source model²², 89.3 MPa considering a Madariaga model²³ and 47.1 MPa considering the Kaneko and Sheare²⁴ model for the mainshock.

The stress drop estimates for the mainshock, the largest foreshock, and two of the largest aftershocks show a clear relationship between magnitude and stress drop, with stress drop increasing with moment magnitude (Fig. 8a). Our estimated stress drops considering the Brune source model are compared with published estimates for other earthquakes in Eastern North America in Fig. 8²⁰,²⁵⁻³¹. We note that our stress drop estimates, including for the mainshock, are lower than the median of stress drop estimates of earthquakes of similar magnitudes in Eastern North America (Fig. 8a).

## Rupture directivity

Using the obtained spectral ratios of the mainshock, we also investigate possible effects caused by a rupture directivity (Fig. 7b). We analyzed the spectral ratios in both along-strike and across-strike directions separately²¹,³², and we observe small to null variations in the spectral ratios and a weak azimuthal variability of corner frequencies. These results suggest that the mainshock ruptured bilaterally.

## Discussion

The analysis of our improved earthquake catalog reveals that at least two structures played a role in the 2020 $M_w$ 5.14 Sparta earthquake sequence and that it may have been a complex event (Fig. 4). A structure is highlighted by the mainly reverse cluster in the seismic sequence, and is inferred as a shallow ( < 2 km) N118', SW dipping fault corresponding to the Little River reverse fault that co-seismically ruptured to the surface¹. The shallow centroid of 1.3 km inferred from the moment tensor inversion further suggests that most of the seismic energy was released in this structure and explains the observation of a surface rupture. The second structure is possibly highlighted by the two identified clusters to the west, clusters mostly strike-slip 1 and mostly strike-slip 2. This structure is interpreted as a sub-vertical

fault with an inferred strike of N143' to N139', and left-lateral strike-slip kinematics that had a blind rupture during the mainshock. This inferred strike-slip structure has a strike consistent with a dominant joint set striking N130-N150' that locally at outcrops shows evidence of displacements even though their timing is not constrained, and that likely favors the development of linear valleys occupied by third to fifth-order streams³³. The hypocenter is relocated near the intersection of these structures, suggesting that this convergence is likely to correspond to an area of stress concentration that eventually led to the mainshock rupture. Strike-slip and thrust faults have previously been identified to interact seismically and rupture together³⁴.

Complex ruptures associated with moderate earthquakes in stable continental regions have been reported for other significant earthquakes. The nearby 2011 $M_w$ 5.8 Mineral, Virginia earthquake was observed to be composed of three subevents from the mainshock teleseismic waveforms³⁵. Analysis of geologic mapping and airborne geomagnetic and gravity surveys also suggest that the 2011 mainshock occurred near a juncture between geologic contacts and near a fault bend or tip³⁶. The 2019 $M_w$ 4.9 Le Teil earthquake, a surface rupturing event in southeast France also broke two distinct fault planes at shallow depths with the nucleation point suggested being in-between the two planes³⁷. Analysis of waveforms for the 1989 Ms 6.3 Ungava earthquake³⁸, another intraplate surface rupturing event in eastern Canada, also shows that the mainshock was composed of two subevents along a thrust plane and a strike-slip plane with different trends. In the SE United States, the historical 1886 M 7 Charleston, South Carolina earthquake has been reported to possibly have two epicenters based on intensity and damage information suggesting a complex rupture³⁹, and modern seismicity reveals that there are two active structures with reverse and strike-slip behavior⁴⁰.

Aftershock locations of the Sparta sequence suggest that the shallower predominantly reverse fault zone may root into a sub-vertical and strike-slip structure, suggesting a complex fault structure. In addition, the small variation in strike between the strike-slip (N143' to N139') and the mostly reverse cluster (N118'), and the gradual variation in the inferred strikes of the strike-slip events in the area of diffuse seismicity (as highlighted by the retrieved focal mechanism for the mostly strike-slip cluster 2 with a strike of N118') (Fig. 6c), further suggest a network of discontinuities or structures that may have been reactivated or even deformed and rotated by the
<figure><img src="image_8.png" /><figcaption>**Fig. 5** | Cross-sections of the seismicity shown in Fig. 4. a Cross-section along and across-strike of AA'. b Cross-section along and across-strike of BB'. Mainshock location is denoted by the yellow star. We include both depth estimates from phase arrival location (deep yellow star) and optimal centroid location (shallow yellow</figcaption></figure>

star). The focal mechanism shows the moment tensor inversion solution plotted in the corresponding side view. The different clusters are plotted according to their colors in Fig. 4. Dashed gray lines correspond to rectangles shown in Fig. 3.

combined rupture. These observations can be interpreted as 1) a change of the strike-slip fault zone geometry, 2) a strike-slip fault zone with several strands with small angles between them 3) a possible additional fault structure, or 4) both clusters are part of a complex structure similar to a positive half "flower structure" corresponding to a diffuse and broader fault zone where a sub-vertical and a SW-dipping strands were activated during the 2020 sequence and merge near the epicentral location. A fourth and minor cluster sub-parallel to the mostly reverse cluster with a dip of 69', and oblique reverse kinematic of one of its earthquakes is further evidence of additional reverse faults rooting at the strike-slip fault, corroborating possibly a positive flower structure ⁴¹. The strike and dip obtained from the aftershock distribution of the mostly reverse cluster are consistent with the parameters of the moment tensor inversion and the location of the surface rupture. Still, the mostly reverse cluster locations are diffuse, which is highlighted by the fitted plane reduced chi-square coefficient. If we projected upwards the fault plane to the surface, it would cross the surface slightly North of the mapped surface rupture (Fig. 4b). This can be explained by a broader fault zone with several strands that is not explained by a simple fault plane, or a possible change in the dip of the structure at shallower depths.

Alternatively, it can be related to uncertainties in the relative earthquake locations, biases in the absolute earthquake locations (Supplementary Method 1) or the inferred cluster plane. The apparent shift of the mostly reverse cluster to the east of the identified surface rupture can be explained by a bias on the absolute locations due to the geometry of the temporary stations (Supplementary Figs. 12–16).

Of the 18 events identified prior to the mainshock in our new catalog, no foreshock is associated with the retrieved clusters and most are located in the epicentral area between the two structures. The first motion polarities of the mainshock also indicate a strike-slip focal mechanism and the hypo-center location based on phase arrivals is deeper than the mostly reverse cluster (Figs. 4b and 5). These observations suggest that the rupture likely initiated in the central area of diffuse seismicity possibly in a strike-slip structure at a deeper region, and then moved upwards to the reverse structure where most of the energy was released. This is similar to other recent observations of moderate-to-large earthquake ruptures that typically started on a fault either conjugate to the main fault rupture (sometimes as foreshocks) ⁴²,⁴³, or with different focal mechanisms⁴⁴,⁴⁵. This hypothesis is consistent with the reported co-seismic surface rupture documents only
<figure><img src="image_9.png" /><figcaption>**Fig. 6** | Focal mechanisms of the *M*<sub>w</sub> 5.1 Sparta mainshock earthquake sequence. a Moment tensor inversion results of the mainshock using gCAP method. b Optimal centroid depth derived from the moment tensor inversion. c Focal mechanisms inverted from first motion polarities. Focal mechanism colors indicate the cluster</figcaption></figure>

they are associated with, and black corresponds to events that were not included in a cluster. For the mainshock, we show both the moment tensor inversion and the first motion polarities result in yellow.

<figure><img src="image_10.png" /><figcaption>**Fig. 7** | Spectral ratio analysis. a The median spectral ratio of the 2020 Sparta mainshock with different EGFs is shown in dark red. The black lines show the different spectral ratios with each EGF at different stations. The dashed green line shows the best fit of the spectral ratio model through Bayesian inference, and the</figcaption></figure>

dashed black lines the inferred corner frequencies. b Median spectral ratios at different along and across-strike quadrants. Stations in each quadrant are shown in the inset map.
<figure><img src="image_11.png" /><figcaption>**Fig. 8** | Comparison of Brune stress drop estimates of earthquakes in the Eastern United States. **a** Stress drop estimates with magnitude. Background solid red line shows the median of stress drops in a moving window of magnitude 1 at each 0.1</figcaption></figure>

magnitude, dashed lines mark the median absolute deviation. **b** Stress drop estimates with depth. Background solid red line shows the median of stress drops in a moving window of 5 km at each 0.5 km, dashed lines mark the median absolute deviation.

reverse faulting¹, even though the authors using InSAR and geodetic surveys inferred a left-lateral strike-slip component, which has been corroborated by Pollitz¹³, with all observations supporting that the strike-slip component occurred mostly at larger depths. The weak azimuthal variation of spectral ratios also suggests that the rupture propagated bilaterally, possibly along both the identified structures.

Our stress drop estimates for four of the largest events in our catalog show a magnitude dependence (Fig. 8a). The stress drop estimates for the area are lower than the median stress drops estimated for other earthquake sequences in Eastern North America, particularly the analyzed foreshocks and aftershocks show stress drops below the estimated median absolute deviation of our compiled stress drop estimates in Eastern North America. Looking at a possible relation with depth, which has been suggested to impact stress drops for earthquakes in Eastern North America¹⁶, our four stress-drop measures do not show a clear depth relationship. Other explanations for the lower stress drops of earthquakes at Sparta can include the spatial dependence observed by Wu and Chapman²⁰ in the 2011 M<sub>w</sub> 5.8 Mineral, Virginia earthquake sequence or the complex fault structure at Sparta. The Sparta mainshock occurred in the Blue Ridge Province whereas other significant earthquakes in the region occurred in different geologic units, potentially leading to differences in the rupture process and stress drop. This is consistent with the recent observation that rupture velocities can potentially impact stress-drop estimates of earthquakes⁴⁷. Our results also clearly show that the Sparta mainshock occurred in a complex fault system, potentially rupturing two faults with distinct tectonic kinematics, which also impacts stress drop estimates⁴⁸. We note that our stress drop estimates are obtained using a simple model that considers an instantaneous circular crack with a simple pulse and constant rupture velocity²², which may not be adequate given our observations at Sparta. Additional analysis is necessary to answer this question and better understand the expected ground motions from moderate Eastern North America earthquakes, possibly using higher-sampling rate seismometers to extend our measurements to smaller earthquakes in the area or using a more complex source model to estimate stress drops⁴⁸.

We show that detailed earthquake catalogs are essential tools to better understand earthquake sequences in stable continental regions. Our analysis reveals that the 2020 Sparta earthquake sequence was a complex sequence that involved multiple structures with different tectonic faulting styles, possibly a diffuse fault zone similar to a positive half flower structure, and the mainshock nucleated near the intersection point of two faults. Such complicated fault ruptures are thought to only exist for large earthquakes⁴⁹, although recent studies have shown that moderate-size earthquakes (M3-5) can also have complex ruptures⁴⁵. Our stress drop estimates for earthquakes

in the sequence also show a low-stress drop estimate compared with other Eastern North America earthquakes, suggesting that more complex models might be necessary to fully comprehend stress drops and as consequence the ground motions of earthquakes in stable continental settings with complex ruptures. The observations of this shallow earthquake sequence provide relevant information that can be used to improve earthquake hazard and risk assessment in stable continental areas, and estimate the impacts of future similar events rupturing in areas with increased hazard as populated areas.

## Methods

### Earthquake detection

To compile the earthquake catalog for this sequence, we used different detection and location techniques that are summarized in a workflow in Fig. 2. We first compiled a template catalog, combining the NEIC catalog with additional detections from a deep learning model. We then used the template catalog to perform additional event detections using MFT. The resulting MFT catalog was then relocated using cross-correlation derived differential travel times.

We first performed an earthquake and phase detection using the Earthquake Transformer (EQTransformer) deep learning model⁵⁰. We used the default EQTransformer model that is trained on a globally distributed dataset of waveforms⁵¹ (Supplementary Method 2 and Supplementary Figs. 17 and 18). The model simultaneously predicts the probability of P and S phase arrivals within a time window. We set the probabilistic thresholds for earthquake detection as 0.3 and 0.1 for P and S wave arrivals. The picked phase arrivals were then associated using the simple association method provided in the EQTransformer package⁵⁰ that searches for phase picks in a moving time window of 15 seconds. The associated events were then located by using the Hypoinverse²² program with a 1-D velocity gradient model derived by Chapman³⁵ for the 2011 M<sub>w</sub> 5.8 Virginia earthquake (Supplementary Table 2).

The resulting deep learning catalog was then combined with the NEIC earthquake catalog in the same region. When an event from the deep learning catalog was detected within 2 seconds of an event in the NEIC catalog we considered it to be the same event. In this case, only the NEIC location and origin time were kept in the template catalog. The combined catalog includes 501 events from August 1, 2020 to November 30, 2020 (Fig. 3).

The resulting template catalog was then used to perform MFT detection, which detects additional earthquakes similar to those in a template catalog through cross-correlation⁵³. We followed a similar procedure to Neves et al.⁵⁴.
Each of the templates is a set of truncated waveforms within 0.5 s before and 3.5 s after the P and S wave arrivals. The MFT detections were made by first performing moving window cross-correlation of the template waveforms with daily continuous data, and then by identifying the periods when high cross-correlation coefficients are present. In order to avoid biases at individual stations, we stack and take the average of the cross-correlation functions of each channel available. Detections are then declared when the coefficients are higher than nine times the median absolute deviation. Detections were performed using a 2–32 Hz two-way fourth-order band-pass filtered data.

## Relocation

Once the smaller events were detected, we constrained their relative locations using the XCORLOC package⁵⁵. Using the same 1-D velocity gradient model from Chapman³⁵ (Supplementary Method 3 and Supplementary Fig. 19), we first derived station term corrections using the template catalog. We then constrained the absolute and relative locations performing phase arrival location and cross-correlation derived differential travel-time relocation⁵⁵. To calculate the cross-correlation derived differential travel-times, we cross-correlated each event with all possible pairs in a 5 km radius. We used −0.5 to 1.5 seconds around the P-phase arrival for the vertical channels, truncating if the S phase was included in the time window, and −0.5 to 2.0 seconds around the S-phase arrival for the horizontal channels. To ensure a millisecond precision, we used the highest sampling rate available and used a spline interpolation to 1 ms in the samples around peak cross-correlation. We used only differential times from measurements of cross-correlation coefficients that are greater than 0.5. Each event pair must have at least 5 differential time measurements and at least 4 with cross-correlation coefficients greater than 0.6 to be relocated. An event was only relocated if it paired with at least two other earthquakes.

The relocated events were then subjected to relative magnitude measurements comparing the peak amplitude ratios of the new events with well paired NEIC templates using the same cross-correlation measurements described before. The magnitudes estimated by the NEIC were kept in the final catalog. The magnitude of completeness ($M_c$), the magnitude above which we consider the catalog to be complete, of the new catalog was estimated using the new magnitudes and the maximum curvature method with a +0.2 correction⁵⁶, returning a magnitude of completeness of −1.5 (Supplementary Method 4 and Supplementary Fig. 20).

## Clustering and fault planes

To identify clusters in the new relocated catalog, we used an unsupervised density-based spatial clustering (DBSCAN) technique¹⁷,⁵⁷. The DBSCAN algorithm takes as parameters the maximum distance for a point to be considered in a neighborhood (EPS), and the minimum number of samples in a neighborhood for an event to be considered a core event. Since we do not have a ground truth, we iteratively searched for the best clustering parameters guided by the R-squared index $R^2$⁵⁸,⁵⁹. $R^2$ compares the variance in the data with the within-cluster variance to search for a solution with lower within-cluster variance and larger variance between different clusters. Its value varies between negative infinity to 1, with higher values indicating generally a better clustering solution. We tested EPS values of 0.1 to 0.5 km with a step of 0.01 km and a minimum number of samples between 25 and 100 with a step of 5 samples (Supplementary Table 3). We used the 3-dimensional Euclidean distance as the metric. Other than the highest $R^2$ index, we also required that the clustering returns more than 1 cluster and assigns clusters to at least 40% of the earthquakes in our relocated catalog. We obtained an EPS value of 0.28 km and a minimum of 25 samples as the optimum parameters for the cluster search.

To obtain a simple geometric model of the seismogenic structures involved in the sequence, we fitted a plane to each of the identified clusters by DBSCAN. We fitted the plane using only events above magnitude of completeness and used an outlier robust Huber regressor⁶⁰.

## Focal mechanisms

To better understand faulting and deformation in the area, we analyzed the radiation patterns of earthquakes in our catalog using two methods.

For the $M_w$ 5.1 mainshock we inverted the moment tensor using the generalized “Cut-and-Paste method”⁶¹,⁶². We calculated the Green’s functions for a depth range of 0.5 to 2.0 km with an interval of 0.1 km using a frequency-wave number method⁶³. We used a 30 s long Pnl wave segment filtered between 0.02 to 0.15 Hz and a 60 s long surface wave segment filtered between 0.02 to 0.1 Hz for the inversion. The search range of strike, dip, and rake angle are 0°–360°, 0°–90° and −180°–180°, respectively, at an interval of 5°, and we set the source duration as 1s⁶⁴. We considered stations at distances less than 3' from the epicenter. The inversion results show that the mainshock is a $M_w$ 5.14 reverse event with some strike-slip component. The two double-couple nodal planes are 350°/70°/115° and 116°/31°/40° at an optimal centroid depth of 1.3 km (Fig. 6a, b). The percentage of the non-DC components is about 26%.

For smaller earthquakes in the catalog, we determined the focal mechanisms using the P-wave first motion polarities constrained by S/P amplitude ratios⁶⁵. P-wave polarities were manually identified for all events with $M > 1.5$ after correcting for the instrument response. To measure S/P amplitude ratios we followed Yang et al.⁶⁶ methodology. We integrated the seismograms to displacement and selected time windows of −0.5 to 1.5 s relative to the phase arrivals as the signal windows and −2.5 to −0.5 s relative to the P-phase arrival as the noise window. We considered the signal amplitude as the difference between the maximum and minimum amplitude values and sum over the 3-component seismograms before computing the SNR. We required at least 15 polarity measurements and a minimum P-to-noise amplitude ratio of 3 to attempt the inversion. We performed the grid search angle with a 5° step, perform 100 trials, and set the fraction of assumed bad polarities measurements to 0.2, the assumed noise in amplitude ratios to 2, and the angle for computing the mechanism probability to 45°. Using these parameters, we were able to estimate the first-motion polarity focal mechanisms for 17 events, including the mainshock (Fig. 6c).

## Spectral ratios, directivity effects, and stress-drop estimates

To estimate the stress-drop of the mainshock and assess if there are rupture directivity effects, we performed a deconvolution of the mainshock seismograms using the smaller aftershock and foreshocks events as EGFs. We followed Ross and Ben-Zion²¹ methodology to perform deconvolution using multiple events as EGF to minimize source effects from the smaller events. We performed this analysis using the direct S-phases to retain directivity information. Due to the large magnitude difference between the mainshock and the largest foreshocks and aftershocks, we selected as EGFs events with a magnitude difference of 2.5 to the mainshock. This resulted in 4 EGF events. We estimated the mainshock corner-frequency ($f_c$) using the Ross and Ben-Zion²¹ proposed two-step framework considering an initial time window of 7 s starting 0.15 s prior to the S-arrival time ($T_s$) based on an estimate from the seismic moment, followed by a more accurate time window of 5.1 s from the first step corner frequency estimate. We calculated spectrums for the considered time windows using a multitaper technique⁶⁷ using only channels with sampling rates equal or higher than 100 Hz and then compute the spectral ratio of the mainshock and EGF event. We then resampled the spectral ratios to 100 points on a log-space between 0.1 and 40 Hz, which is 80% of the lowest Nyquist frequency²⁰. We stacked the spectral ratios of all events at each station to obtain a station spectral ratio, and then stacked all station spectral ratios to obtain the final smoothed spectral ratio ($r$) that we used to estimate the mainshock corner frequency. Only segments of the spectral ratios with SNR larger than 3 were used in the stacking process. The noise window corresponds to a window with the same length as the signal window that ends 5 seconds prior to the P-phase arrival. To estimate the corner-
frequencies we fitted to the spectral ratios by Bayesian inference a source spectral model following<sup>20,22</sup>:

$$r(f) = \frac{\Omega^1}{\Omega^2} \left( \frac{(1 + f/f_c^2)^{ny}}{(1 + f/f_c^1)^{ny}} \right)^{1/\gamma}, \quad (1)$$

where 1 and 2 refer to the mainshock and EGF, respectively, Ω corresponds to the low-frequency asymptote, γ and n are constants. We fixed γ to 2 in our model but allow n that corresponds to the fall-off rate which is usually fixed at 2 to be inferred. We assessed the uncertainty of the corner frequency estimate by considering the 95% confidence interval. The fit was performed using the PyMC3<sup>68</sup> library that uses a Markov chain Monte Carlo algorithm to sample the posterior distributions. As described by Wu and Chapman<sup>20</sup>, for the prior distributions we considered log-normal distributions of the corner frequencies and the moment ratio, and a normal distribution of the likelihood function.

With the estimated corner-frequency, we calculated the stress drop assuming a circular crack model<sup>69</sup>:

$$\Delta\sigma = \frac{7M_0}{16} \left(\frac{f_c}{\kappa\beta}\right), \quad (2)$$

where Δσ is the stress-drop, M₀ the moment, β the shear-wave velocity and κ is a constant that depends on the source model considered<sup>22-24</sup>.

To assess the directivity effects, we split the analyzed stations into four quadrants in relation to the moment tensor fault plane of 116°, respectively two along-strike and two across-strike quadrants. We then stacked the previously obtained station spectral ratios of all the stations in each quadrant to obtain a “quadrant” spectral ratio. Similarly to Ross and Ben-Zion<sup>21</sup> and Zhou et al.<sup>32</sup>, we compared the quadrant spectral ratios to identify possible differences and determine if these are significant differences using t-statistics.

We additionally estimated the stress drop of the largest foreshock and two of the largest aftershocks in our catalog using the multi-window coda spectral ratio method<sup>20</sup>. The spectral ratio was obtained by comparing the coda window of these events with the mainshock. This method returns more stable spectral ratio estimates but cancels out both sources’ directivity effects. We only compared events with similar first motion polarities as the mainshock. Similarly to Wu and Chapman<sup>20</sup>, to measure the coda spectral ratios we used 5 subwindows of 5.1 s with a 50% overlap starting at 1.5T_s for stations at distances less than 100 km and 2T_s for more distant stations. For this, we required again a SNR larger than 3. Finally, the corner-frequencies are once again estimated using Bayesian inference.

## Data availability

The seismic waveforms used in this study were retrieved from the IRIS Web Services (https://service.iris.edu/), including the following seismic networks: ET (CERI Southern Appalachian Seismic Network), CO (South Carolina Seismic Network, https://doi.org/10.7914/SN/CO), KY (Kentucky Seismic and Strong Motion Network, https://doi.org/10.7914/SN/KY), N4 (Central and Eastern US Network, https://doi.org/10.7914/SN/N4), SE (Southeastern Appalachian Cooperative Seismic Network), US (United States National Seismic Network, https://doi.org/10.7914/SN/US). The compiled high-resolution catalog (Supplementary Data 1), the estimated first motion polarities focal mechanism catalog (Supplementary Data 2) and the detection catalog (Supplementary Data 3), together with 3-Dimensional plots of the earthquake catalog with best fitting planes (Supplementary Fig. 10) and the inferred fault planes (Supplementary Fig. 11) can be accessed at https://doi.org/10.6084/m9.fishshare.22689250.

## Code availability

Analyses for this study were performed using established codes referenced in the methods section. We additionally used Obspy<sup>70</sup>, Numpy<sup>71</sup>, Scipy<sup>72</sup>, Matplotlib<sup>73</sup>, Plotly and Generic Mapping Tools<sup>74</sup> for data handling and

figure creation. Scripts applying the codes are available from M.N. (email: neves@geoazur.unice.fr) on request.

Received: 2 May 2023; Accepted: 12 March 2024;

Published online: 29 March 2024

## References

1. Figueiredo, P. M. et al. The Mw 5.1, 9 August 2020, Sparta Earthquake, North Carolina: The first documented seismic surface rupture in the Eastern United States. GSA Today 32, 4–11 (2022)

2. Wicks Jr, C. W. & Chiu, J.-M. Surface Rupture on a Secondary Fault Associated with the 8 August 2020 Mw 5.1 Sparta North Carolina Earthquake. Seismic Record 2, 59 (2022).

3. Adams, J., Wetmiller, R. J., Hasegawa, H. S. & Drysdale, J. The first surface faulting from a historical intraplate earthquake in North America. Nature 352, 617 (1991).

4. Stein, S. & Liu, M. Long aftershock sequences within continents and implications for earthquake hazard assessment. Nature 462, 87 (2009).

5. Liu, M. & Stein, S. Mid-continental earthquakes: Spatiotemporal occurrences, causes, and hazards. Earth-Sci. Rev. 162, 364 (2016).

6. Levine, J. S. F. et al. Kinematic, deformational, and thermochronologic conditions along the gossan lead and fries shear zones: constraining the western-eastern blue ridge boundary in Northwestern North Carolina. Tectonics 37, 3500 (2018).

7. Thomas, W. A. The lapetan rifted margin of southern Laurentia. Geosphere 7, 97 (2011).

8. Levandowski, W., Herrmann, R. B., Briggs, R., Boyd, O. & Gold, R. An updated stress map of the continental United States reveals heterogeneous intraplate stress. Nat. Geosci. 11, 433 (2018).

9. Lund Snee, J.-E. & Zoback, M. D. Multiscale variations of the crustal stress field throughout North America. Nat. Commun. 11, 1951 (2020).

10. Bollinger, G. A. & Wheeler, R. L. The Giles County, Virginia, Seismic Zone. Science 219, 1063 (1983).

11. Daniels, C. & Peng, Z. A 15-year-long catalog of seismicity in the Eastern Tennessee Seismic Zone (ETSZ) using matched filter detection. Earthquake Res. Adv. 3, 100198 (2023).

12. Powell, C. A. et al. A seismotectonic model for the 300-kilometer-long Eastern Tennessee Seismic Zone. Science 264, 686 (1994).

13. Pollitz, F. F. Shallow Fault Slip of the 2020 M 5.1 Sparta, North Carolina, Earthquake. Seismolog. Res. Lett. 94, 2831–2839 (2023) https://doi.org/10.1785/0220230127.

14. Horton, S., Withers, M. M., Cramer, C. H. & Withers, H. Source parameters of the M5.1 Sparta, NC earthquake of 9 August 2020. Seismolog. Res. Lett. 92, 1439 (2021).

15. Richter, C. F. Elementary Seismology Series of books in Geology (W. H. Freeman, 1958)

16. van den Ende, M. et al. Distributed Acoustic Sensing for aftershock monitoring: The case of the 2019 Mw 4.9 Le Teil earthquake. EarthArXiv6646 [Preprint], https://doi.org/10.31223/X5DX26 (2024).

17. Ester, M., et al. A density-based algorithm for discovering clusters in large spatial databases with noise. In kdd (1996) pp. 226–231.

18. Boore, D. M. Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra. Bull. Seismolog. Soc. Am. 73, 1865 (1983).

19. Hough, S. E. Shaking from injection-induced earthquakes in the central and Eastern United States. Bull. Seismolog. Soc. Am. 104, 2619 (2014).

20. Wu, Q. & Chapman, M. Stress-Drop Estimates and Source Scaling of the 2011 Mineral, Virginia, Mainshock and Aftershocks. Bull. Seismolog. Soc. Am. 107, 2703 (2017).

21. Ross, Z. E. & Ben-Zion, Y. Toward reliable automated estimates of earthquake source properties from body wave spectra. J. Geophys. Res.: Solid Earth 121, 4390 (2016).
22. Brune, J. N. Tectonic stress and the spectra of seismic shear waves from earthquakes. *J. Geophys. Res.* (1896-1977) **75**, 4997 (1970).

23. Madariaga, R. Dynamics of an expanding circular fault. *Bull. Seismolog. Soc. Am.* **66**, 639 (1976).

24. Kaneko, Y. & Shearer, P. M. Variability of seismic source spectra, estimated stress drop, and radiated energy, derived from cohesive-zone models of symmetrical and asymmetrical circular and elliptical ruptures. *J. Geophys. Res.: Solid Earth* **120**, 1053 (2015).

25. Shi, J., Kim, W.-Y. & Richards, P. G. The 2002 M5 Au Sable Forks, NY, earthquake sequence: Source scaling relationships and energy budget. *Bull. Seismolog. Soc. Am.* **88**, 531 (1998).

26. Viegas, G., Abercrombie, R. E., & Kim, W.-Y. The 2002 M5 Au Sable Forks, NY, earthquake sequence: Source scaling relationships and energy budget. *J. Geophys. Res.: Solid Earth* **115**, B07310, https://doi.org/10.1029/2009JB006799 (2010).

27. Boatwright, J. & Seekins, L. Regional spectral analysis of three moderate earthquakes in Northeastern North America. *Bull. Seismolog. Soc. Am.* **101**, 1769 (2011).

28. Viegas, G. Source parameters of the 16 July 2010 Mw 3.4 Germantown, Maryland, Earthquake. *Seismolog. Res. Lett.* **83**, 933 (2012).

29. Boyd, O. S., McNamara, D. E., Hartzell, S. & Choy, G. Influence of lithostatic stress on earthquake stress drops in North America. *Bull. Seismolog. Soc. Am.* **107**, 856 (2017).

30. Huang, Y., Ellsworth, W. L. & Beroza, G. C. Stress drops of induced and tectonic earthquakes in the central United States are indistinguishable. *Sci. Adv.* **3**, e1700772 (2017).

31. Daniels, C. et al. The 15 February 2014 Mw 4.1 South Carolina Earthquake Sequence: Aftershock Productivity, Hypocentral Depths, and Stress Drops. *Seismolog. Res. Lett.* **91**, 452 (2019).

32. Zhou, Y. et al. Seismological Characterization of the 2021 Yangbi Foreshock-Mainshock Sequence, Yunnan, China: More than a Triggered Cascade. *J. Geophys. Res.: Solid Earth* **127**, e2022JB024534 (2022).

33. Merschat, A. et al. Long-lived crustal deformation and reactivation in the Eastern Blue Ridge near Sparta, NC. In *Geological Society of America 2022 Meeting Abstracts with Programs*, Vol. 54–5, https://doi.org/10.1130/abs/2022AM-383089 (2022).

34. Rivero, C., Shaw, J. H. & Mueller, K. Oceanside and Thirtymile Bank blind thrusts: Implications for earthquake hazards in coastal southern California. *Geology* **28**, 891 (2000).

35. Chapman, M. C. On the Rupture Process of the 23 August 2011 Virginia Earthquake. *Bull. Seismolog. Soc. Am.* **103**, 613 (2013).

36. Shah, A. K. et al. Subsurface geologic features of the 2011 central Virginia earthquakes revealed by airborne geophysics. In *The 2011 Mineral, Virginia, Earthquake, and Its Significance for Seismic Hazards in Eastern North America* Horton, J. W., Jr., Chapman, M. C., and Green, R. A. (Geological Society of America, 2015).

37. De Novellis, V. et al. Coincident locations of rupture nucleation during the 2019 Le Teil earthquake, France and maximum stress change from local cement quarrying. *Commun. Earth Environ.* **1**, 1 (2020).

38. Bent, A. L. The 1989 (MS 6.3) Ungava, Quebec, earthquake: A complex intraplate event. *Bull. Seismolog. Soc. Am.* **84**, 1075 (1994).

39. Dutton, C. E., The Charleston earthquake of August 31, 1886, U.S. Geological Survey, Ninth Annual Report 1887–88 203 (1889).

40. Chapman, M. C., Beale, J. N., Hardy, A. C. & Wu, Q. Modern seismicity and the fault responsible for the 1886 Charleston, South Carolina, Earthquake. *Bull. Seismolog. Soc. Am.* **106**, 364 (2016).

41. Huang, L. & Liu, C.-y Three types of flower structures in a divergent-wrench fault zone. *J. Geophys. Res.: Solid Earth* **122**, 10,478 (2017).

42. Yao, D., Huang, Y., Peng, Z. & Castro, R. R. Detailed Investigation of the Foreshock Sequence of the 2010 Mw 7.2 El Mayor-Cucapah Earthquake. *J. Geophys. Res.: Solid Earth* **125**, e2019JB019076 (2020).

43. Chuang, L. Y. et al. Foreshocks of the 2010 Mw 6.7 Yushu, China earthquake occurred near an extensional step-over. *J. Geophys. Res.: Solid Earth* **128**, e2022JB025176 (2023).

44. Ozacar, A. A. & Beck, S. L. The 2002 Denali Fault and 2001 Kunlun Fault earthquakes: complex rupture processes of two large strike-slip events. *Bull. Seismolog. Soc. Am.* **94**, S278 (2004).

45. Kato, A. and Ueda, T. Source fault model of the 2018 Mw 5.6 northern Osaka earthquake, Japan, inferred from the aftershock sequence. *Earth Planets Space* **71**, 11, https://doi.org/10.1186/s40623-019-0995-9 (2019).

46. Long, L. T. The mechanics of natural and induced shallow seismicity: a review and speculation based on studies of Eastern U.S. Earthquakes. *Bull. Seismolog. Soc. Am.* **109**, 336 (2019).

47. Chounet, A., Vallée, M., Causse, M. & Courboulex, F. Global catalog of earthquake rupture velocities shows anticorrelation between stress drop and rupture velocity. *Tectonophys. Phys. Earthquake Rupture Prop.* **733**, 148 (2018).

48. Liu, M., Huang, Y., & Ritsema, J., Characterizing Multisubevent Earthquakes Using the Brune Source Model. *Bull. Seismolog. Soc. Am.* **113**, 577–591 (2023).

49. Hamling, I. J. et al. Complex multifault rupture during the 2016 Mw 7.8 Kaikoura earthquake, New Zealand. *Science* **356**, eaam7194 (2017).

50. Mousavi, S. M., Ellsworth, W. L., Zhu, W., Chuang, L. Y. & Beroza, G. C. Earthquake transformer—an attentive deep-learning model for simultaneous earthquake detection and phase picking. *Nat. Commun.* **11**, 3952 (2020).

51. Mousavi, S. M., Sheng, Y., Zhu, W. & Beroza, G. C. STanford Earthquake Dataset (STEAD): a global data set of seismic signals for AI. *IEEE Access* **7**, 179464 (2019).

52. Klein, F. W., User's Guide to HYPOLINVERSE-2000, a Fortran Program to Solve for Earthquake Locations and Magnitudes, USGS Numbered Series 2002-171 (U.S. Geological Survey, 2002).

53. Peng, Z. & Zhao, P. Migration of early aftershocks following the 2004 Parkfield earthquake. *Nat. Geosci.* **2**, 877 (2009).

54. Neves, M., Peng, Z., & Lin, G., A high-resolution earthquake catalog for the 2004 Mw 6 Parkfield earthquake sequence using a matched filter technique. *Seismolog. Res. Lett.* **94**, 507 (2022).

55. Lin, G. The source-specific station term and waveform cross-correlation earthquake location package and its applications to California and New Zealand. *Seismolog. Res. Lett.* **89**, 1877 (2018).

56. Woessner, J. & Wiemer, S. Assessing the quality of earthquake catalogues: estimating the magnitude of completeness and its uncertainty. *Bull. Seismolog. Soc. Am.* **95**, 684 (2005).

57. Pedregosa, F. et al. Scikit-learn: machine learning in python. *J. Mach. Learn. Res.* **12**, 2825 (2011).

58. Edwards, A. W. F. & Cavalli-Sforza, L. L. A method for cluster analysis. *Biometrics* **21**, 362 (1965).

59. Scott, A. J. & Symons, M. J. Clustering methods based on likelihood ratio Criteria. *Biometrics* **27**, 387 (1971).

60. Huber, P. J. Robust estimation of a location parameter. *Annals Math. Stat.* **35**, 73 (1964).

61. Zhu, L. & Helmberger, D. V. Advancement in source estimation techniques using broadband regional seismograms. *Bull. Seismolog. Soc. Am.* **86**, 1634 (1996).

62. Zhu, L. & Ben-Zion, Y. Parametrization of general seismic potency and moment tensors for source inversion of seismic waveform data. *Geophys. J. Int.* **194**, 839 (2013).

63. Zhu, L. & Rivera, L. A. A note on the dynamic and static displacements from a point source in multilayered media. *Geophys. J. Int.* **148**, 619 (2002).

64. Somerville, P., Collins, N., Abrahamson, N., Graves, R., and Saikia, C., Earthquake source scaling and ground motion attenuation relations for the central and Eastern United States, Final Report to the US Geological Survey, Contract No. 99HQGR0098 (2001).
65. Hardebeck, J. L. & Shearer, P. M. Using S/P amplitude ratios to constrain the focal mechanisms of small earthquakes. *Bull. Seismolog. Soc. Am.* **93**, 2434 (2003).

66. Yang, W., Hauksson, E. & Shearer, P. M. Computing a large refined catalog of focal mechanisms for Southern California (1981–2010): Temporal stability of the style of faulting. *Bull. Seismolog. Soc. Am.* **102**, 1179 (2012).

67. Prieto, G. A. The multitaper spectrum analysis package in python. *Seismolog. Res. Lett.* **93**, 1922 (2022).

68. Salvatier, J., Wiecki, T. V. & Fonnesbeck, C. Probabilistic programming in Python using PyMC3. *PeerJ Comput. Sci.* **2**, e55 (2016).

69. Eshelby, J. D. & Peierls, R. E. The determination of the elastic field of an ellipsoidal inclusion, and related problems. *Proc. R. Soc. London. Ser. A. Math. Phys. Sci.* **241**, 376 (1957).

70. Beyreuther, M. et al. ObsPy: a python toolbox for seismology. *Seismolog. Res. Lett.* **81**, 530 (2010).

71. Harris, C. R. et al. Array programming with NumPy. *Nature* **585**, 357 (2020).

72. Virtanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in python. *Nat. Methods* **17**, 261 (2020).

73. Hunter, J. D. Matplotlib: A 2d graphics environment. *Comput. Sci. Eng.* **9**, 90 (2007).

74. Wessel, P. et al. The Generic Mapping Tools Version 6. *Geochem. Geophys. Geosyst.* **20**, 5556 (2019).

## Acknowledgements

This work was partially supported by a USGS NEHRP grant G20AP00039. M. N. is supported by a Ph.D. fellowship from the Portuguese research foundation Fundação para a Ciência e Tecnologia (FCT), Grant SFRH/BD/ 139033/2018. This work is also supported by the National Natural Science Foundation of China (NSFC, 42030311). We thank Dr. Steve Horton from CERI/University of Memphis for sharing the temporary aftershock deployment data through IRIS DMC. Finally, we thank Luigi Passarelli, Bing Li and one anonymous reviewer and the editor Luca Dal Zilio for their insightful comments.

## Author contributions

M.N., L.C. and Z.P. conceived the study. M.N. performed the matched filter detection, earthquake relocation, and first motion polarities, clustering, stress drop and directivity analysis. L.C. performed the Earthquake

Transformer detection and deep learning detection tests. W.L. and S.N. performed mainshock moment tensor inversion. All the authors contributed to interpret the results. M.N., Z.P. and P.F. wrote the article. All the authors contributed to review the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-024-01316-8.

**Correspondence** and requests for materials should be addressed to Miguel Neves.

**Peer review information** *Communications Earth & Environment* thanks Yijian Zhou, Luigi Passarelli and Bing Qiuyi Li for their contribution to the peer review of this work. Primary Handling Editors: Luca Dal Zilio, Joe Aslin and Aliénor Lavergne. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024