Distinct Yet Adjacent Earthquake Sequences near the Mendocino Triple Junction: 20 December 2021 $Mw$ 6.1 and 6.0 Petrolia, and 20 December 2022 $Mw$ 6.4 Ferndale Open Access

At the Mendocino Triple Junction (MTJ) in northern California, subduction of the Gorda plate beneath the North American plate at the southernmost section of the Cascadia subduction zone transitions into strike‐slip motion along the northernmost section of the San Andreas fault (Fig. 1). As one of the most seismically active areas in California, the MTJ regularly experiences moderate‐to‐large earthquakes as well as aseismic slip (Materna et al., 2018) with varied mechanisms on a diverse range of faults in the area. Notably, on 25 April 1992, the reverse‐fault $Mw$ 7.2 Cape Mendocino earthquake occurred in the accretionary wedge within the overriding North American plate (Oppenheimer et al., 1993; McCrory et al., 2012).

The MTJ recently experienced three moderate‐magnitude strike‐slip earthquakes. On 20 December 2021, an offshore strike‐slip $Mw$ 6.1 earthquake at 20:10:20 UTC on the Mendocino fracture zone (Fig. 1, blue mechanism) triggered a deeper strike‐slip $Mw$ 6.0 earthquake (Fig. 1, green mechanism), 11 s later at 20:10:31 UTC, ∼30 km away on a fault within the subducting Gorda slab. The second $Mw$ 6.0 earthquake had a distributed aftershock zone extending ∼10 km to the north (Yeck et al., 2023). These 2021 Petrolia earthquakes caused limited damage. Then a year later, on 20 December 2022 at 10:34:24 UTC, the $Mw$ 6.4 Ferndale earthquake (Fig. 1, black mechanism) nucleated ∼15 km northwest of the Petrolia aftershock zone, rupturing unilaterally in the east‐northeast (ENE) direction toward populated areas located inland. Ferndale aftershocks extended ∼40 km ENE of the hypocenter (Shelly et al., 2024). The Ferndale earthquake resulted in two deaths and extensive structural damage in several Humboldt County communities (Stein et al., 2023).

This study seeks to delineate the fault structures at depth activated by and understand the relationship between the 2021 $Mw$ 6.1 and 6.0 Petrolia earthquake sequence and the subsequent 2022 $Mw$ 6.4 Ferndale sequence. We therefore created an enhanced relocated earthquake catalog for the MTJ region, spanning both the sequences. Aided by the development of deep‐learning seismic phase pickers, automated workflows for creating high‐resolution earthquake catalogs from continuous seismic data are emerging as standard tools in observational seismology (e.g., Zhu et al., 2022) to illuminate fault structure, sequence evolution, and earthquake processes (e.g., Yoon et al., 2023).

We applied the automated workflow used in Yoon et al. (2023) (Fig. S1, available in the supplemental material to this article) to continuous seismic data between 1 December 2021 and 1 June 2023 at 66 stations (Fig. 1, inverted triangles) to create an enhanced relocated catalog for the MTJ region, which includes both the Petrolia and Ferndale earthquake sequences. These stations from the BK, CE, NC, NP, and PB networks (Data Set S0), all located within the region in Figure 1 (39.5° to 41.5° latitude, −126° to −123° longitude), included 21 three‐component broadband (HH*) seismometers, 24 three‐component (HN*) accelerometers, 8 three‐component short‐period (EH*) sensors, and 13 single‐component vertical short‐period (EHZ) sensors.

First, we applied the EQTransformer deep‐learning model for automated event detection and phase picking (Fig. S1a; Mousavi et al., 2020) at each station. EQTransformer trained on a global data set of ∼1.3 million seismograms was one of the best‐performing deep‐learning pickers in a large‐scale benchmark test (Münchmeyer et al., 2022). We used nonstandard EQTransformer input parameters (Table S1), including lower detection thresholds and greater overlap between time windows, to increase the number of detected events and picks. EQTransformer output 823,177 P and 834,530 S picks across the 66 stations.

Next, the Rapid Earthquake Association and Location (REAL) grid‐search event associator (Fig. S1b; Zhang et al., 2019) merged P and S picks from different stations into their originating events. To focus on the recent MTJ seismicity, REAL was setup to associate only events within 0.7° (∼78 km) of the region centered around (latitude 40.5° and longitude −124.3°) and depths 0–40 km. To maximize the number of associated events, we set low association thresholds for REAL, requiring the minimum of three P and three S picks with a total of six picks on at least three stations (Table S2). To calculate travel times, we used a 1D velocity model (Fig. S2, Table S3) from the Northern California Seismic System (NCSS) for the Mendocino area (U.S. Geological Survey [USGS] Menlo Park, 1966). REAL associated 18,940 total events.

Third, the Hypocenter Inversion with Stein Variational Inference (HypoSVI) algorithm computed absolute event locations (Fig. S1c; Smith et al., 2021). HypoSVI, a probabilistic Bayesian sampling method for event location, overcomes limitations of earthquake location methods based on linearized inversion. Initialized with randomly selected locations, HypoSVI undergoes an iterative Stein Variational Gradient Descent optimization process, converging to estimate the posterior distribution of the 3D hypocenter and its uncertainties (Fig. S3). HypoSVI inputs include parameters in Table S4, P and S arrival times for associated events, and a trained EikoNet model (Smith et al., 2020) for the MTJ region. Given a known 1D MTJ velocity model (Fig. S2, Table S3), EikoNet trains a neural network to solve the eikonal (high‐frequency) wave equation within a predefined 3D volume bounded by latitude 39.5° to 41.5°, longitude −126° to −123°, and depth −2 to 41 km (Fig. 1), which includes all 66 stations and all earthquakes from the Petrolia and Ferndale sequences. EikoNet can then rapidly compute grid‐free travel times between any two locations in this volume. HypoSVI located 18,911 total events.

We then removed false detections occurring in the coda of larger (M > 4) earthquakes (Fig. S1d) using the same empirical algorithm from Yoon et al. (2023). These false detections usually were from 15 to 45 s after the origin time of the preceding earthquake, had relatively few (6–12) picks with low EQTransformer probabilities, and large (M >4) calculated magnitudes. Some false detections may remain in, or a few real earthquakes may be mistakenly missing from, the HypoSVI catalog. After eliminating 1150 such false detections comprising 6.1% of the 18,911 HypoSVI‐located events, 17,761 HypoSVI events remained. HypoSVI locations and $ML$ magnitudes were also computed for 131 missed events in the ComCat catalog (U.S. Geological Survey [USGS] Earthquake Hazards Program, 2017) that EQTransformer and REAL failed to identify for a total of 17,892 events in the HypoSVI catalog (Data Set S1).

Finally, HypoSVI locations were further refined with the GrowClust algorithm for precise relative earthquake location (Fig. S1e; Trugman and Shearer, 2017), keeping the same $ML$ magnitudes. Along with input parameters (Table S5), GrowClust needs P and S cross‐correlation differential travel times computed between pairs of nearby events at a common station, which are more precise than EQTransformer pick times used in HypoSVI locations. We computed ∼51 million differential times in the time domain, with longer time windows for S than for P picks (parameters in Table S6). 12,373 events (∼70% of HypoSVI events) were relocated by GrowClust (Data Set S2), though several M > 3 earthquakes and all three M > 6 earthquakes were not relocated. To include these larger earthquakes, our final enhanced relocated MTJ catalog (Fig. 1) combines 12,373 GrowClust locations with 29 HypoSVI locations for all M > 3 events not belonging to a GrowClust cluster (Data Set S3).

We first compare the HypoSVI deep‐learning‐enhanced catalog against the ComCat catalog (USGS Earthquake Hazards Program, 2017) contributed by NCSS (Data Set S4; USGS Menlo Park, 1966) to evaluate what earthquake information was added by the automated deep‐learning workflow (Fig. 2). Counting only events within the blue box in Figure 1, a match between HypoSVI and ComCat events (Fig. 2a,b, blue) was declared if their origin times matched within 5 s and their epicenters matched within 25 km. We use the HypoSVI catalog as a basis for comparison, since all earthquakes were located in a consistent way, instead of the GrowClust relocated catalog, where some isolated events not belonging to a cluster of nearby events, such as larger events with complex waveforms, were not relocated. HypoSVI locations are more uncertain farther west offshore due to poor station coverage (Fig. S4a, second row). HypoSVI 95% confidence interval uncertainties were 5.14 ± 2.60 km in the north–south (latitude) direction, 9.75 ± 3.32 km in the east–west (longitude) direction, and 3.44 ± 1.76 km for depth, and 0.36 ± 0.41 s for origin time (Fig. S5), with lower uncertainties at higher magnitudes (Fig. S5, first row and Fig. S6). Calculated $ML$ values (equation 1) reasonably matched the ComCat magnitudes $Mcat$⁠, which were mostly coda duration magnitudes for $Mcat<3$⁠, local magnitudes for $3<Mcat<3.5$⁠, and moment magnitudes for $Mcat>3.5$⁠, with mean and standard deviation of −0.25 ± 0.26 units for their residual value $ML−Mcat$ (Fig. S7a,b).

The HypoSVI catalog provides a more complete view of seismicity in both the sequences, finding 14,084 previously unknown, mostly $0<ML<2$⁠, earthquakes (Fig. 2a,b, cyan), ∼16 times greater than the 860 ComCat events during this time. Thanks to permissive EQTransformer parameter settings (Table S1; Mousavi et al., 2020), the HypoSVI catalog is missing only 8 of 860 (Fig. 2a,b, red) ComCat events, mostly from incorrect association: five events were within <6 s of another event, two offshore and deeper (>30 km depth) events were poorly localized, and one event (origin time 31 March 2023 22:57:19.55 UTC) was missing continuous seismic data. The workflow separately associated the $Mw$ 6.1 and 6.0 Petrolia earthquakes, which happened only 11 s apart (Fig. S8), but errors in pick labeling and association led to mislocating the offshore $Mw$ 6.1 event by >20 km from its human‐reviewed ComCat location (Figs. 3d–f, 4c, and 5a, blue stars). Notably, NCSS uses stricter criteria than our three‐station minimum (Table S2) to release event origins to ComCat.

We also compare the HypoSVI (Smith et al., 2021) deep‐learning‐enhanced catalog against template‐matching earthquake catalogs (Fig. 3). These template‐matching catalogs (Fig. S6d) available for the Petrolia sequence from 20 December 2021 to 11 January 2022 (Yeck et al., 2023; Fig. 3b) and the Ferndale sequence from 19 December 2022 to 21 January 2023 (Shelly et al., 2024; Fig. 3c) used ComCat (USGS Earthquake Hazards Program, 2017) earthquake waveforms as templates, which were cross correlated with continuous seismic data to detect smaller unknown earthquakes with similar waveforms. Magnitudes $Mtemp$ in the template‐matching catalogs were calculated in a pairwise correlation manner, so they differ slightly from $ML$ values (equation 1) in the HypoSVI catalog, with mean and standard deviation of 0.18 ± 0.31 units for their residual value $ML−Mtemp$ and systematically higher residuals at lower magnitudes (Fig. S7c); these magnitude estimates are nevertheless close enough for comparison.

Although the HypoSVI and template‐matching catalogs overlap by 34%, with 4,120 of 12,245 common events (Fig. 3a–c, blue, Fig. 3e, and Fig. S4c), we found events that belong only to the template‐matching catalog or the HypoSVI catalog, but not to both, consistent with other studies comparing deep‐learning and template‐matching catalogs (Scotto di Uccio et al., 2022). 4,724 of 12,245 (∼38%) events, only in the template‐matching catalog (Fig. 3a–c, magenta), have the lowest magnitudes between −0.5 and 1, concentrate closer in time to larger earthquakes, and locate in tight clusters near known larger events (Fig. 3d). In contrast, 3,401 of 12,245 (∼28%) events, only in the HypoSVI catalog (Fig. 3a–c, cyan), have slightly higher magnitudes from 0.5 to 2 and occur more evenly distributed in time with more scattered locations throughout the region (Fig. 3f).

These results suggest that deep‐learning and template‐matching earthquake catalogs contain complementary information and could be combined to create a more complete catalog than either approach alone. Deep learning can identify background seismicity at lower magnitudes missed by template matching, whereas template matching might better fill in spatial and temporal details of seismicity near already known larger earthquakes. One could identify new earthquakes with a deep‐learning approach and then use them as template waveforms in another round of template‐matching event detection.

Deep‐learning‐based automated workflows for earthquake catalog generation have key advantages: They can easily be applied directly to continuous seismic data without needing a preexisting catalog with event waveforms and phases, and they need minimal human intervention and quality control with advance parameter tuning. Such workflows could help situational awareness of an evolving earthquake sequence in a response situation (Yoon et al., 2023), and are being explored by regional seismic networks for improving completeness and quality of routine earthquake catalogs (e.g., Walter et al., 2021).

The enhanced relocated MTJ earthquake catalog, which combines GrowClust (Trugman and Shearer, 2017) relocations with HypoSVI (Smith et al., 2021) locations for all nonclustered M >3 events, illuminates the relationship between the 2021 Petrolia and 2022 Ferndale earthquake sequences (Fig. 4), and their detailed fault structure at depth (Fig. 5). GrowClust locations for M > 1 events are overall comparable to ComCat (USGS Earthquake Hazards Program, 2017) locations (Fig. S6a,c, second row). GrowClust locations are more tightly clustered compared to those from the HypoSVI catalog (Fig. S6b), where errors from automatic EQTransformer (Mousavi et al., 2020) picks resulted in more scattered locations, but their hypocentral parameters show good overall agreement (Fig. S4b). Using a 1D velocity model could bias our absolute event locations, especially offshore where station coverage is poor. GrowClust relative locations are less affected, although an inaccurate velocity model may bias the dip angle of relocated seismicity. Going forward, enhanced catalog locations could be improved with 3D velocity models, and EQTransformer‐derived picks could in turn also constrain 3D velocity models.

The 2021 $Mw$ 6.1 and 6.0 Petrolia sequence (Fig. 4a) was a complex, multiple‐fault, dynamically triggered rupture (Yeck et al., 2023). It started offshore on 20 December 2021 20:10:20 UTC as an $Mw$ 6.1 right‐lateral strike‐slip earthquake (Figs. 4c and 5a, blue stars) on the east–west trending Mendocino fracture zone (MFZ), which forms the boundary between the Pacific and Gorda plates. Eleven seconds later at 20:10:31 UTC, a deeper $Mw$ 6.0 earthquake happened on a southwest–northeast‐trending fault located ∼30 km eastward (Figs. 4c and 5a,d, green stars) onshore in the subducting Gorda slab, with a left‐lateral strike‐slip mechanism (Yeck et al., 2023). Aftershocks in a distributed zone extended ∼10 km to the north of the onshore $Mw$ 6.0, and offshore east of the $Mw$ 6.1 along the MFZ (Fig. 4a), with mostly strike‐slip moment tensors. However, an ∼10 km gap remained between the easternmost offshore MFZ aftershocks of the $Mw$ 6.1 and aftershocks of the onshore $Mw$ 6.0.

Then a year later, on 20 December 2022 10:34:24 UTC, the $Mw$ 6.4 Ferndale earthquake (Fig. 4b) started ∼15 km northwest of the Petrolia aftershock zone (Figs. 4c and 5a,b, black stars) and ruptured unilaterally ENE on a left‐lateral strike‐slip fault within the subducting Gorda slab (Shelly et al., 2024), activating an aftershock zone extending ∼40 km ENE of the hypocenter (Figs. 4b and 5a,b). This intraslab aftershock zone at depth is oriented differently from the northwest–southeast trend of mapped surface faults, which Gulick et al. (2001) interpreted as evidence for partial decoupling between the Gorda and North American plates. Aftershock locations along this zone were more concentrated in space and persistent in time to the west near the hypocenter, while more distributed in space inland to the east (Fig. 5a,b). Most Ferndale aftershock mechanisms in this zone were strike slip like the $Mw$ 6.4 mainshock, with some normal‐fault mechanisms farther east (Fig. 4b). One exception was a spatially isolated, deeper $Mw$ 5.4 aftershock on 1 January 2023; the trend of its immediate aftershocks suggests right‐lateral strike‐slip motion on a northwest–southeast‐trending fault (Figs. 4b and 5a,c; Stein et al., 2023). This $Mw$ 5.4 aftershock located ∼20 km away from the primary ENE‐trending Ferndale aftershock zone may have been a triggered earthquake in the slab mantle. The Ferndale sequence occurred where historical seismicity levels from 1974 to December 2021, according to ComCat (Fig. S9; USGS Earthquake Hazards Program, 2017), were relatively lower compared to the surrounding region.

We found that the 2021 $Mw$ 6.1 and 6.0 Petrolia and 2022 $Mw$ 6.4 Ferndale earthquake sequences were distinct in space and time. The sequences were adjacent to each other in space but more separated in time (Fig. 2c); most Petrolia aftershocks occurred in the first few hours following the $Mw$ 6.1 and 6.0 events (Fig. 2d,e, magenta) before dropping to nearly background seismicity levels for several months, until December 2022 when the $Mw$ 6.4 Ferndale mainshock started its own aftershock sequence. Neither sequence had detected foreshocks in the 20 days before their respective mainshocks. Both the sequences exhibited similar rapid decay rates and low aftershock productivities (Fig. 2d) typical of the MTJ region, which has unusually short‐lived and unproductive aftershock sequences (Hardebeck et al., 2019; Gomberg and Bodin, 2021), possibly due to the young age (<10 Ma) of the Gorda plate (Dascher‐Cousineau et al., 2020; Wu et al., 2023). The two sequences had little overlap in space (Figs. 4a,b and 5a,d), as they happened on different faults at different depths, where subsurface conditions in nearby areas within the MTJ can differ considerably (Chen and McGuire, 2016). The onshore Petrolia sequence mostly occurred deeper (25–35 km depth) in the Gorda slab mantle (Fig. 5a,d, red), whereas most Ferndale aftershocks were confined to a shallower and narrower 17–20 km depth range in the uppermost Gorda slab (Fig. 5a,b, purple‐green; Shelly et al., 2024), near the subduction interface where fluid pressure is elevated (Guo et al., 2021). A small <5 km gap (Fig. 5d, red arrow) separates the northernmost end of the 2021 $Mw$ 6.0 onshore Petrolia aftershock zone (Fig. 5a,d, red) from the 2022 Ferndale aftershock zone (Fig. 5a,b, purple‐green), where a cluster of seismicity extends slightly southward (Fig. 5a). Inside this gap, one small shallower (<20 km depth) Petrolia aftershock cluster adjoins the southward‐extending Ferndale aftershock cluster (Fig. 5d). But it is unknown why the spatially adjacent Ferndale sequence happened one year later, whereas the onshore $Mw$ 6.0 Petrolia earthquake happened only 11 s after the offshore $Mw$ 6.1 with a larger (∼10 km) spatial gap between their aftershocks and ∼30 km gap between areas of peak slip (Yeck et al., 2023). Fault complexity was also observed after the 1992 $Mw$ 7.2 Cape Mendocino earthquake, where the reverse‐fault mainshock occurred in the overriding North American plate (Fig. 5a,b,d, gray), and, on the next day, strike‐slip $Mw$ 6.5 and 6.6 aftershocks happened offshore ∼30 km to the west on different faults (Oppenheimer et al., 1993). Complex ruptures and sequences on multiple faults with variable timing, including those not on the main plate boundary interfaces, are important to consider within the range of possible earthquakes in a complicated region like the MTJ (Hagerty and Schwartz, 1996).

We created a deep‐learning‐enhanced relocated catalog encompassing two recent MTJ earthquake sequences: (1) 20 December 2021 $Mw$ 6.1 and 6.0 Petrolia, and (2) 20 December 2022 $Mw$ 6.4 Ferndale. This new catalog identified over 14,000 earthquakes not in the ComCat catalog, mostly magnitude 0–2, while missing only 8 of 860 ComCat events; it also separated the $Mw$ 6.1 and 6.0 Petrolia earthquakes that occurred 11 s apart, though with a >20 km location error for the offshore $Mw$ 6.1. Deep‐learning and template‐matching approaches improved earthquake catalog completeness in complementary ways; deep learning found more background seismicity, whereas template‐matching catalogs from Shelly et al. (2024) and Yeck et al. (2023) found the smallest M < 0 earthquakes. We found that the Petrolia sequence was distinct in space and time from, but spatially adjacent to, the Ferndale sequence. Both sequences with similar fast decay rates and low aftershock productivities typical of the MTJ region occurred in the downgoing Gorda slab. But the shallower Ferndale sequence was confined to the uppermost slab near the subduction interface, while the onshore Petrolia sequence happened deeper in the slab mantle. While the Petrolia and Ferndale sequences were close in space (<5 km gap) but farther apart in time (1 yr), the $Mw$ 6.1 and 6.0 Petrolia earthquakes happened farther apart in space (∼10–30 km gap; Yeck et al., 2023) but closer in time (11 s). Deep‐learning‐enhanced earthquake catalogs could help situational awareness during an evolving earthquake sequence, especially in a complex plate boundary region such as the MTJ where earthquake behavior spans a wide range of faults, ruptures, and sequences with variable location and timing.

ComCat catalog (doi: 10.5066/F7MS3QZH), contributed by Northern California Seismic System (NCSS; https://earthquake.usgs.gov/fdsnws/event/1/, last accessed July 2023). Continuous seismic data were accessed from Northern California Earthquake Data Center (NCEDC) (https://ncedc.org/, last accessed July 2023), with data contributed by BK (doi: 10.7932/BDSN), CE (doi: 10.7914/b34q-bb70), NC (doi: 10.7914/SN/NC), NP (doi: 10.7914/SN/NP), and PB (https://www.fdsn.org/networks/detail/PB/, last accessed July 2023) networks. Seismology software: EQTransformer (https://github.com/smousavi05/EQTransformer, last accessed May 2023), REAL (https://github.com/Dal-mzhang/REAL, last accessed May 2023), EikoNet (https://github.com/Ulvetanna/EikoNet, last accessed May 2023), HypoSVI (https://github.com/Ulvetanna/HypoSVI, last accessed May 2023), GrowClust (https://github.com/dttrugman/GrowClust, last accessed May 2023). Scripts for creating enhanced catalog from seismic data, plotting (https://gitlab.com/cyoon1/ferndale2022, last accessed March 2024). ObsPy 1.4.0 (doi: 10.1088/1749-4699/8/1/014003; https://www.obspy.org/, last accessed May 2023) was used for seismological processing and visualization, and Generic Mapping Tools (GMT) 6.4.0 (doi: 10.1029/2019GC008515; Wessel et al., 2019; https://www.generic-mapping-tools.org, last accessed May 2023) for maps. U.S. Geological Survey (USGS) Quaternary Faults database (Figs. 1 and 3–5; https://www.usgs.gov/natural-hazards/earthquake-hazards/faults; last accessed May 2023). Major plate boundaries (Fig. 1; doi: 10.1029/2001GC000252). The supplemental material includes Tables S1–S6 with input parameters; Figures S1–S9 with analysis details; Data Sets S0–S4 with station data and earthquake catalogs in text format (Yoon and Shelly, 2024; doi: 10.5281/zenodo.10621116); and References.

The authors acknowledge that there are no conflicts of interest recorded.

Waveform data, metadata, or data products for this study were accessed through the Northern California Earthquake Data Center (NCEDC; doi: 10.7932/NCEDC). The authors acknowledge the contributions of the technical and monitoring staff at the Northern California Seismic System, jointly operated by the U.S. Geological Survey (USGS) and Berkeley Seismo Lab for seismic data collection and the ComCat earthquake catalog used in this study. The authors thank all developers and maintainers of open‐source seismology software listed in Data and Resources. Jeff McGuire, Daniel Trugman, and an anonymous reviewer provided helpful article reviews, and Editor‐in‐Chief Keith Koper facilitated the submission process. Luke Blair, Jeanne Hardebeck, and Shane Detweiler helped us with the data release process for the supplemental data sets. Jeff McGuire sent us the velocity model, Gabrielle Tepp shared the most recent dML station correction file (Uhrhammer et al., 2011) to calculate local magnitudes, Sarah Minson helped with Generic Mapping Tools (GMT) plots, and Rob Graves reviewed conference abstracts about this work. Jack Wilding and Jonny Smith helped with setting up and running EikoNet and HypoSVI. The authors had helpful discussions about Mendocino Triple Junction (MTJ) earthquakes with Andy Barbour, Max Liu, Kathryn Materna, Saeed Mohanna, Jay Patton, and the USGS Operational Aftershock Forecasting team. The authors thank all presenters and participants in the special Ferndale earthquake session at the virtual 2023 Northern California Earthquake Hazards Workshop for sharing information. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.