Low-frequency earthquakes track the motion of a captured slab fragment
Links
- More information: Publisher Index Page (via DOI)
- Document: XML
- Open Access Version: USGS Accepted Manuscript
- Download citation as: RIS | Dublin Core
Main Text
With the advent of plate tectonics theory came the recognition that the behavior of triple junctions, where three tectonic plates meet, is central to understanding the evolution of the present-day plate configuration (1, 2). Along the west coast of northern California, the Mendocino triple junction marks the convergence of the Pacific, North American, and Gorda (southern Juan de Fuca) plates (3, 4) and forms the boundary between the strike-slip San Andreas Fault system and the Cascadia subduction zone (Fig. 1). Although often idealized as a single point, in three dimensions the triple junction is a broad, complex region. It encompasses a network of major and minor faults that exhibit high rates of seismicity and accommodate both seismic and aseismic deformation across the three interacting plate boundaries, with high rates of deformation extending well into the Gorda plate (5). The triple junction at the northern end of Cascadia is also seismically active and may be equally complex (6–8).
The nature of the transition between the Cascadia subduction zone and the San Andreas fault system has long been debated (4, 9). Despite extensive study, uncertainty persists regarding fundamental tectonic features, particularly the nature of the southern edge of Gorda subduction. One model holds that a slab gap (window) forms at the southern edge of the northward migrating subducting Gorda plate, leading to shallow asthenospheric upwelling (9–14). A contrasting model proposes that this region is underlain by a subducted Farallon slab fragment captured by the Pacific plate and sharing its northward motion (15). In the latter scenario, this captured slab, termed the Pioneer fragment (16, 17) would lie beneath North American plate crust and track the northward motion of the southern edge of Gorda subduction at depth, preventing formation of slab window in this location. However, no evidence has directly constrained the modern-day kinematics.
Tremor and low-frequency earthquakes (LFEs) might provide such evidence. Tremor is a low-amplitude, long-duration seismic signal first identified in the early 2000s (18), which is associated with slow slip events downdip of the seismogenic zone in subduction zones (19, 20). Meanwhile, LFEs are very small individual earthquakes with rapid recurrence that form the longer duration tremor (21–23). LFEs are often detected using matched-filter analysis, with events identified by a given waveform template considered a “family” of events. While tremor frequently has large location and depth uncertainties, LFE source locations can sometimes be determined much more accurately by identification of P and S phases, often on stacked records combining many LFE detections (24). LFE and tremor activity worldwide dominantly occurs on major plate bounding faults (25).
An anomalous zone of tremor, offset 50‒100 km westward of the main band of Cascadia tremor, was identified by Wech (26). Shelly et al. (27) identified and located 27 families of constituent LFEs in this zone, showing that they aligned along a zone at 22‒29 km depth dipping ~45º toward the north‒northeast, extending eastward and downdip from a neighboring zone of persistent microseismicity (M<3) (Fig. 1). Given the association with major plate bounding faults and the short (~2 day) recurrence times, Shelly et al. (27) argued that these Mendocino LFEs likely reflected slip on a high-rate plate bounding fault, perhaps at the southern edge of Cascadia subduction. However, that study was unable to constrain the sense of slip occurring in the LFEs, leaving unanswered fundamental questions as to the nature of this boundary and its associated tectonics.
Plausible arguments could be made for various orientations of LFE slip, each with distinct tectonic implications (Fig. 2). LFEs might reflect slip between the Gorda slab and North American plate crust at the edge of a slab window, leading to normal faulting on a structure defined by LFE hypocentral alignment (here determined as strike=293º, dip=48º – see Methods). In contrast, if the Pioneer captured slab fragment hypothesis is correct, then we might observe relative motion between the Pacific and Gorda plates, for which we would expect right-lateral strike-slip motion on the LFE-defined fault. Finally, it is conceivable that the Pioneer fragment has been captured by the Pacific plate but that North American crust lies to the north of the LFEs, leading to oblique reverse and right-lateral strike-slip faulting to accommodate relative motion between the Pacific and North American plates (Fig. 2). The actual faulting orientation, if it could be determined, would distinguish between these distinctly different tectonic scenarios.
To address this uncertainty, we performed two complementary analyses aimed at constraining the orientation of fault slip and tectonic significance of Mendocino LFE activity.
Focal mechanism analysis
We computed a composite first-motion focal mechanism for the LFEs families. Because low signal-to-noise for these LFEs usually precludes the unambiguous determination of P-wave first motions, even on stacked records, we developed a strategy based on comparison of P-wave arrivals with nearby earthquakes of known polarity (28, 29) (see Methods). These nearby earthquakes are located within a zone of persistent seismicity at ~22 km depth (3, 30, 31) with diverse focal mechanisms, which has previously been referred to as the “Garberville swarm” (32) (Fig. 1).
The resulting focal mechanism is shown in Fig. 2. The northeast-dipping nodal plane (strike=292º, dip=34º) aligns remarkably well with the orientation determined by the hypocentral alignment of the families, (strike=293º, dip=48º), which is also similar to that previously estimated visually (strike=290º, dip=45º) (27). The focal mechanism’s rake angle of 177º on this plane indicates strike-slip motion. Notably, despite the variety of mechanisms observed in the reference earthquakes (Fig. 1), none closely matched the LFE mechanism.
Tidal analysis
Strong tidal modulation of tremor and LFEs has been well documented in multiple fault zones, including the Nankai subduction zone (33, 34), Cascadia subduction zone (35, 36), and the central San Andreas fault (37, 38). Because the LFE families documented by Shelly et al. (27) occur with short recurrence intervals, as opposed to being driven by large transient slow slip events such as those that occur farther north on the Cascadia subduction zone, we expect the LFE rate to be proportional to the tidal stress amplitude, as has been shown in other settings that host nearly continuous LFEs (e.g., the San Andreas fault near Parkfield). The stresses imposed by tides vary according to fault and slip orientation, so we can estimate the optimal fault geometry by determining which fault orientation and sense of slip maximizes the occurrence of LFEs during positive (slip encouraging) shear stress changes imparted by the tides. We modeled these effects to determine whether the Mendocino LFEs are also tidally modulated and the range of slip orientations consistent with this modulation (see Methods).
For all families combined, we find strong tidal modulation that is broadly consistent with either normal faulting or strike-slip faulting, meaning that an excess number of LFEs occur during tidally encouraging shear stress, but inconsistent with reverse faulting (Fig. 2C,D). For a strike of 293º and dips of 35‒55º, we find that the strongest and most systematic modulation (i.e. the maximum slope) occurs for dip=40º, rake=180º (Fig. 2D), closely matching one nodal plane from the independently determined focal mechanism (strike=292º, dip=34º, rake=177º) (Fig. 2). A secondary peak corresponding to left-lateral strike-slip motion for dip=35º is tectonically implausible and inconsistent with the P-wave focal mechanism.
Tectonic Implications
Fig. 2A shows the predicted focal mechanisms for multiple tectonic scenarios, where the fault strike and dip are set by the LFE hypocentral alignment, with the rake varied according to that required to accommodate the MORVEL 2010 model plate motions (39). We find that the mechanism theoretically predicted by Pacific-Gorda relative plate motions closely matches the computed LFE composite focal mechanism, while other tectonic scenarios are outside the range of acceptable solutions (Fig. 2B). Furthermore, the frequent, ~2-day recurrence times observed for this LFE activity would likely be difficult to produce with a low-slip-rate fault, suggesting they occur on a high-slip-rate fault (27) similar to other worldwide observations of plate-bounding tremor and LFEs (25). Therefore, we argue that the Mendocino LFEs likely reflect relative motion between the Gorda and Pacific plates.
The determined strike-slip mechanism is inconsistent with the hypothesis that a slab window forms directly south of the Gorda slab in this area. Instead, it argues that Pacific-like plate motion is actively occurring southeast of the surface expression of the triple junction, beneath the North American plate, matching the hypothesis of capture of a partially subducted former Farallon slab fragment by the Pacific plate (15, 17). We therefore infer that the LFEs reflect slip on an eastward extension of the Mendocino transform fault (Fig. 1).
A puzzling feature noted by Shelly et al. (2025) is that the LFEs occur on the southern edge of an aseismic zone that extends ~15 km south from the southern limit of Gorda slab seismicity (Fig. 3a,b), with distinctly lower P-wave velocities at ~23-30 km depth compared with corresponding depths to the north and (especially) south (17, 31). Although this low velocity zone is not apparent in seismic refraction data (11), perhaps due to limited resolution (Figs. S1-S2), it is prominent in tomography studies and has been previously interpreted as thickened North American crust (17, 31). However, the fault motion of the LFEs implies that rather than being a thickened portion of present-day North American crust, this low-velocity zone is instead accreted to the southern edge of the subducting Gorda slab, shifting the subduction boundary south from the edge of the petrologic slab.
A plausible source of this low-velocity material is the southernmost portion of North American accretionary prism, which formed as sediments were scraped off the subducting Gorda plate. Fig. 3C diagrams our proposed model. Because the western boundary of the accretionary prism (the deformation front) is west of the northward projection of the San Andreas fault (Fig. 1), the intervening zone is subjected to high rates of north‒south compression, more so than the underlying Gorda slab, which is subducting obliquely northeastward relative to North America (Fig. 2A). The Gorda slab has a distinct southward dip at its southern edge, as previously recognized from seismic velocity structure and seismicity (17, 31, 40) (Fig. 3c). Although the Gorda plate has been previously interpreted from seismic velocity structure as purely east-dipping beneath the coast (41, 42), a southward dip at its southernmost edge (Fig. 3c) does not conflict with these data (Figs. S1, S3). To escape the compression imparted by the Pacific plate, a portion of accretionary prism above the Gorda slab may have been subducted, effectively transferring these low-velocity sediments back to the Gorda plate from which they originated before accretion onto the continent. This model bears some resemblance to a previously proposed model in which the southernmost portion of the North American accretionary prism was subducted due to eastward migration of the San Andreas fault and Cascadia subduction front (43, 44). However, we envision this accretionary prism subduction to have occurred progressively as the Mendocino fault moved northward with Pacific plate motion, rather than limited to the time of an eastward jump in the San Andreas fault. Subduction of this low-velocity zone may continue until it reaches the eastern edge of the Pioneer fragment and escapes the associated north‒south compression, east of the observed LFEs (Fig. 3b).
Multiple lines of evidence support the interpretation of a wedge of material subducting above the petrologic Gorda slab at its southern edge. For example, the McCrory et al. (40) model suggests the slab is >10 km deep at the southern end of the Cascadia deformation front, implying subduction of a thick layer above the slab (Fig. 1). In addition, repeating earthquakes observed on the Mendocino fault as shallow as ~10 km depth suggest that right-lateral strike-slip motion occurs at slip rates similar to those reflected by repeating earthquakes occurring much deeper (>20 km depth), adjacent to the inferred subducting Gorda slab (45).
Further support for this model is provided by the 1992 M 7.2 Cape Mendocino earthquake, which occurred at ~10 km depth, ~10 km above the inferred subducting slab, with low-angle thrust slip consistent with subduction motion (Figs. 1,3). Although this event was originally interpreted to represent rupture of the subduction interface (46, 47), more recent studies, which showed the rupture was much shallower than the Gorda slab (as expressed by the slab’s seismic velocity structure and seismicity), instead interpreted this event as thrusting on a minor fault within the upper plate (40, 42). We infer that the 1992 M 7.2 event occurred on the primary subduction interface fault, and we propose that this interface is ~10 km above the petrologic slab at its southern end, with the depth perhaps influenced by a stress concentration imparted by the upper boundary of thickened lower crust of the colliding Pacific plate (48). This interpretation is consistent with the eastward curve of the southernmost offshore deformation front (Fig. 1), which does not project to the slab itself but aligns with the proposed shallower interface and the 1992 M 7.2 rupture (Fig. 3d). This is also consistent with thrust-type mechanisms observed for the 1991 M 6.0 Honeydew earthquake (49). The Northern California Seismic System focal mechanism solutions for this event all indicate reverse slip with northeast-dipping (22-46º) nodal planes (Fig. 1). A shallow dip is most consistent with potential subduction interface slip. Furthermore, ongoing nearby thrust-type microseismicity extends downdip toward the east (Figs. 1,3). This hypothesis effectively reconciles disparate prior interpretations.
Accretionary prisms are well recognized for their considerable fluid content, initially forming from a collection of saturated sediments scraped off the downgoing plate (50). As the former accretionary prism subducts, it would be increasingly heated and compressed, likely leading to high fluid pressures. High fluid pressures are strongly associated with tremor and LFE activity, as they may enable brittle failure under conditions where distributed ductile deformation would otherwise dominate (51–55). Tremor has been previously observed offshore Japan within the accretionary prism (e.g., 56), though at somewhat shallower depths than the 22-29 km determined for Mendocino LFEs. Alternatively, fluids enabling Mendocino tremor and LFEs could originate from dehydration of the Gorda or Pioneer slabs. Fluids might also contribute to the generation of persistent, sometimes swarm-like microearthquakes (32) adjacent to LFE activity. These microearthquakes exhibit diverse focal mechanisms reflecting a combination of strike-slip and normal faulting (Figure 1), which seems to preclude their occurrence on the same fault as the LFEs. However, their consistent, long-lived activity suggests they may be associated with tectonic deformation; they might reflect fluid-assisted tear faulting immediately above the Pioneer/Gorda boundary. Hybrid shear-dilational faulting, associated with high fluid pressure and fracture mesh faulting (29, 57) is also possible.
We explore the proposed plate geometry in Fig. 4. The moderately dipping LFE zone suggests that the LFEs may occur along the upper surface of a northward-dipping Pioneer fragment. If so, despite its northward dip, the Pioneer fragment is not actively subducting to the north as the southeasterly strike accommodates Pacific‒Gorda relative motion (Fig. 2). The proposed geometry is attractive because it is similar to that associated with the main band of tremor and LFEs on the subduction interface in Cascadia (Fig. 1) and elsewhere, downdip of the primary seismogenic zone, except that the relative motion in this case is strike-slip rather than thrust.
In this captured slab model, the San Andreas fault would sole into a shallowly eastward-dipping, midcrustal detachment on the surface of the Pioneer fragment, as it translates northward beneath westernmost North America (15, 17) (Fig. 4). In fact, active-source seismic imaging has resolved a pronounced highly reflective, eastward dipping structure beneath the northern San Andreas system, though it has been debated whether or not the San Andreas cuts through this structure (12, 58). The Maacama and possibly Bartlett Springs faults would also be expected to sole into this detachment, depending on its eastward extent. Although no associated seismicity on the detachment itself has been reported, the existence of such a structure might impart shear tractions that could help explain the partitioning of slip across multiple similarly oriented faults of the northern San Andreas system (17). Although the detachment might slip aseismically, it alternatively could be capable of generating infrequent large earthquakes, perhaps extending megathrust rupture from the Cascadia subduction zone itself. Similar lower-crustal structures have been observed farther south in California, leading to a proposal that a midcrustal detachment fault on the surface of a fossil slab may also underlie other portions of the San Andreas fault system (59, 60) and indeed much of California (61). This idea has been revisited more recently (62), yet debate continues regarding the prevalence of fossil slabs versus slab windows and asthenospheric upwelling (13).
Our results have important implications for subduction zone modeling and seismic hazards. The potentially large discrepancy between petrologic slab boundaries and the faults that accommodate subduction contrasts with usual assumptions that these are practically interchangeable (40, 63). The shallower interface location proposed here, although a small portion of the Cascadia subduction zone, could bring slip closer to the surface and alter slab coupling or dynamic rupture models, which have relied on the geometry of slab models (64, 65). Perhaps more critically, the proposed detachment fault along the surface of Pioneer fragment does not appear in current fault models (66, 67) and is not considered in the USGS National Seismic Hazard Model (68). If this fault is seismogenic, it could represent a substantial unaccounted seismic hazard, given the high rate (~51mm/yr) of relative motion between the Pacific and North American plates, exceeding the rate of Gorda subduction (Fig. 2).
The extended length of interaction between the Gorda and Pacific/Pioneer plates at depth onshore might increase the likelihood of interactions between large earthquakes on the Cascadia subduction interface and those along the San Andreas fault (69, 70). Furthermore, high seismicity rates near the Mendocino triple junction, including three earthquakes of M 6+ since 2021, contrast starkly with the relative seismic quiescence observed throughout most of Cascadia. This may increase the likelihood that the next Cascadia megathrust earthquake would nucleate within the triple junction, perhaps contributing to an apparently shorter average recurrence interval for southern Cascadia compared to northern Cascadia (71). Continued monitoring of tremor and LFEs in this area might help to ascertain any future variations in deep fault slip rate, particularly in the aftermath of neighboring large earthquakes. Such slip variations may be subtle and not detectable by other means.
References and Notes
29. D. R. Shelly, J. L. Hardebeck, W. L. Ellsworth, D. P. Hill, A new strategy for earthquake focal mechanisms using waveform-correlation-derived relative polarities and cluster analysis: Application to the 2014 Long Valley Caldera earthquake swarm. J. Geophys. Res. Solid Earth 121, 8622–8641 (2016).
47. D. Oppenheimer, J. Eaton, A. Jayko, M. Lisowski, G. Marshall, M. Murray, R. Simpson, R. Stein, G. Beroza, M. Magee, G. Carver, L. Dengler, R. McPherson, L. Gee, B. Romanowicz, F. Gonzalez, W. H. Li, K. Satake, P. Somerville, D. Valentine, The Cape Mendocino, California, Earthquakes of April 1992: Subduction at the Triple Junction. Science 261, 433–438 (1993).
66. A. Plesch, S. Marshall, J. Shaw, SCEC Community Fault Model (CFM), Statewide California Earthquake Center (2024); https://zenodo.org/records/13685611.
68. M. D. Petersen, A. M. Shumway, P. M. Powers, E. H. Field, M. P. Moschetti, K. S. Jaiswal, K. R. Milner, S. Rezaeian, A. D. Frankel, A. L. Llenos, A. J. Michael, J. M. Altekruse, S. K. Ahdi, K. B. Withers, C. S. Mueller, Y. Zeng, R. E. Chase, L. M. Salditch, N. Luco, K. S. Rukstales, J. A. Herrick, D. L. Girot, B. T. Aagaard, A. M. Bender, M. L. Blanpied, R. W. Briggs, O. S. Boyd, B. S. Clayton, C. B. DuRoss, E. L. Evans, P. J. Haeussler, A. E. Hatem, K. L. Haynie, E. H. Hearn, K. M. Johnson, Z. A. Kortum, N. S. Kwong, A. J. Makdisi, H. B. Mason, D. E. McNamara, D. F. McPhillips, P. G. Okubo, M. T. Page, F. F. Pollitz, J. L. Rubinstein, B. E. Shaw, Z.-K. Shen, B. R. Shiro, J. A. Smith, W. J. Stephenson, E. M. Thompson, J. A. Thompson Jobe, E. A. Wirth, R. C. Witter, The 2023 US 50-State National Seismic Hazard Model: Overview and implications. Earthq. Spectra 40, 5–88 (2024).
69. C. Goldfinger, K. Grijalva, R. Bürgmann, A. E. Morey, J. E. Johnson, C. H. Nelson, J. Gutiérrez-Pastor, A. Ericsson, E. Karabanov, J. D. Chaytor, J. Patton, E. Gràcia, Late Holocene Rupture of the Northern San Andreas Fault and Possible Stress Linkage to the Cascadia Subduction Zone. Bull. Seismol. Soc. Am. 98, 861–889 (2008).
70. C. Goldfinger, J. Beeson, B. Black, A. Vizcaino, C. H. Nelson, A. Morey, J. R. Patton, J. Gutiérrez-Pastor, C. Romsos, M. D. Walzcak, Unravelling the dance of earthquakes: Evidence of partial synchronization of the northern San Andreas fault and Cascadia megathrust. Geosphere, doi: 10.1130/GES02857.1 (2025).
71. C. Goldfinger, C. H. Nelson, A. E. Morey, J. E. Johnson, J. Patton, E. Karabanov, J. Gutiérrez-Pastor, A. Ericsson, E. Gràcia, G. Dunhill, R. J. Enkin, A. Dallimore, T. Vallier, Turbidite Event History-Methods and Implications for Holocene Paleoseismicity of the Cascadia Subduction Zone. Prof. Pap. No 1661-F US Geol. Surv., doi: 10.3133/pp1661F (2012).
73. NCEDC, Northern California Earthquake Data Center, UC Berkeley Seismological Laboratory (2014); https://doi.org/10.7932/NCEDC.
74. D. R. Shelly, A Catalog of Low-Frequency Earthquakes at the Southern Edge of Cascadia Subduction, U.S. Geological Survey data release (2025); https://doi.org/10.5066/P1TCKK7G.
Additional Information
Funding: This work was supported in part by National Science Foundation Award #1848302 (AMT).
Author Contributions: Conceptualization: DRS, AMT, KZM; Methodology: DRS, AMT, RJS; Investigation: DRS, AMT, RJS; Visualization: DRS, AMT, RJS; Writing – original draft: DRS; Writing – review & editing: DRS, AMT, KZM, RJS
Competing interests: The authors declare that they have no competing interests
Data and materials availability: Waveform, phase, and event catalog data were obtained from the Northern California Earthquake Data Center (NCEDC) (73) and the EarthScope Seismological Facility for the Advancement of Geoscience (SAGE) Data Management Center (https://ds.iris.edu/ds/nodes/dmc/, last accessed April 15, 2025). The catalog of LFE detections is available from Shelly (74). Updated version of Waldhauser and Schaff (75) catalog obtained from https://nocaldd.ldeo.columbia.edu/, last accessed May 1, 2025. Fault information is from the U.S. Geological Survey and California Geological Survey, Quaternary fault and fold database for the United States, accessed May 15, 2025, at: https://www.usgs.gov/natural-hazards/earthquake-hazards/faults. This study involved no new materials preparation. All other data needed to evaluate the conclusions in the paper are present in the paper or the supplementary materials.
Acknowledgements
Stations used in this study are operated by the U.S. Geological Survey (NC network, doi: 10.7914/SN/NC), UC Berkeley (BK network, doi:10.7932/BDSN), and the Earthscope Consortium (PB network, https://www.fdsn.org/networks/detail/PB/, last accessed May 12, 2025). Figure 1 was created using GMT6 (72). We thank Hao Guo for providing a copy of the velocity model from Guo et al. (42). We greatly appreciate reviews from Michael Bostock, Jeff McGuire, Tom Parsons, and Robert McPherson, and discussions with Kevin Furlong and Janet Watt. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
Supplementary Materials
Figures
Mendocino triple junction seismicity and focal mechanisms. (A) Cascadia margin. White lines show subducting plate contours from McCrory et al. (40), in 5-km increments. Black lines show faults; line with triangles indicates the Cascadia subduction deformation front (western edge of the accretionary prism). Blue dots mark tremor epicenters from pnsn.org/tremor, (last accessed June 13, 2025) from 2015-01-01 through 2024-07-31. Red box indicates the region shown in (B). The anomalous tremor zone studied here is highlighted. (B) Mendocino triple junction region. Low-frequency earthquake (LFE) family locations are shown as white-outlined diamonds, color-coded by depth (27). Seismicity 1984-2021, updated from Waldhauser and Schaff (75), is shown as small dots, also color-coded by depth. Focal mechanisms are similarly color-coded and show the 1992 M 7.2 Cape Mendocino earthquake (47), the 1991 M 6.0 Honeydew earthquake (ID=nc228027), and examples of other nearby thrust events including a 2022 M 2.7 (focal mechanism, ID=nc73661611) and a 2021 M 3.6 (moment tensor, ID=nc73629686) earthquake. Also shown are the LFE composite mechanism determined in this study and a representative sample of focal mechanisms of microseismicity (M<3) events used for LFE focal mechanism analysis. White triangles show the locations of seismic stations used for LFE analysis. The zone of high crustal compression (Figure 3c) is also indicated (white dashed box). The inferred buried eastward extension of the Mendocino fault is shown as a black line near LFEs, with a dashed line connecting it to the offshore Mendocino fault. The shaded region just to the north is the nearly aseismic region of observed low P-wave velocity (Low-V) at depth (17, 31). Motions for Gorda and Pacific plates (See Methods) relative to the North American plate are also shown. Black lines show faults from the Quaternary Fault and Fold Database (major fault zones labeled), which are generally shallow faults not directly related to deeper structures. Topography/bathymetry is from Ryan et al. (76).
Comparison of theoretical and determined slip orientations. (A) Relative plate motion from the MORVEL 2010 model (39) is shown with vectors (see Methods). Focal mechanisms corresponding with these relative motions resolved onto the low-frequency earthquake (LFE) hypocentral plane fit (strike=293º, dip=48º) are plotted. Strike, dip, and rake of these mechanisms are indicated above each. (B) Composite focal mechanism for all LFE families computed LFE P-wave polarity analysis. Symbols indicated estimated compressional (plus) or dilatational (circle) polarity observation projected onto the focalsphere, labeled with corresponding station name (corresponding colors). Light lines indicate potential acceptable solutions assuming hypocentral uncertainties of 0.5 km and a 10% polarity error (see Methods). (C) Tidal modulation slope for a strike of 293º, showing a maximum tidal modulation slope for dip=40º, rake=180º. Gray shading shows the range of rake angles corresponding with each faulting type (labeled), with unshaded regions denoting oblique faulting. (D) Bar plot showing the ratio of actual to expected event numbers as a function of tidal shear stress for the maximum slope orientation. Lower inset shows fit and associated root mean squared error (RMSE). Tidal modulation increases systematically with increasing stress. Corresponding mechanism for this orientation is plotted in inset. Note the close correspondence between the theoretical mechanism representing Pacific-Gorda plate convergence in (A) (light blue), and those determined from P-wave analysis and tidal modulation.
Velocity structure, seismicity, and proposed formation of low-velocity zone. Relocated earthquakes and low-frequency earthquakes (LFEs) are from Shelly et al. (27), colored by depth. White stars show 1991 M 6.0 (Honeydew) and 1992 M 7.2 (Cape Mendocino) thrust-faulting earthquakes, with locations from (75). (A) Contours of P-wave velocity in km/s at 27.25 km depth (near mean depth of LFEs) from Furlong et al. (17) are shown in blue; note the alignment of LFEs with contours separating a high-velocity zone to the south with the low-velocity zone to the north. Locations of transects shown in parts C-E are indicated, with dashed boxes showing regions of plotted seismicity for each. (B) Interpreted plate structure at 27 km depth, corresponding with velocity structure. (C) Cartoon illustration showing a south-north cross section a-a’ near the coastline, as indicated in (A), dashed boxes indicate zone), with proposed evolution of downwarping at the southern edge of the Gorda slab, and progressive subduction of North American accretionary prism to present-day low-V corridor at the southern edge of the Gorda slab. Circle with central dot indicates subduction motion of the Gorda plate out of the page. Third panel (at right) overlays actual relocated seismicity. (D) West-east cross-section b-b’, as indicated in part (A), showing shallow seismicity lineation containing 1991 M 6.0 and 1992 M 7.2 thrust earthquakes, projecting to the deformation front, as defined by bathymetry (Fig. 1). We hypothesize that a portion of the former North American (N.A.) accretionary prism subducts above the south-dipping southern edge of the Gorda slab. Dashed black line indicates the approximate upper boundary of Gorda intraslab seismicity, while the white dashed line indicates the same at the Gorda’s southern edge. (E) Southwest-northeast cross section c-c’ (oriented at 20º), across the LFE zone, showing the interpreted structure and the alignment of LFE hypocenters. Arrows indicate approximate motion relative to North America. The circle with cross indicates motion into the page. Pacific-Gorda relative motion is estimated to be nearly pure strike-slip along the LFE strike.
Proposed plate geometry at the Mendocino triple junction. Low-frequency earthquakes (LFEs) are interpreted to occur on the boundary between the captured Pioneer fragment and a low-velocity zone accreted to the subducting Gorda slab, adjacent to a zone of microearthquakes (MicroEQ). Circles with a central dot denote motion out of the page, while the circle with cross indicates motion into the page. We propose that the low-velocity (low-V) zone is accreted to the Gorda slab, and that LFEs occur along the southern boundary of Gorda subduction. A detachment surface with unknown seismic potential likely separates the Pioneer fragment from the North American (N.A.) crust. LFE motion is nearly pure strike-slip, with the southeasterly strike serving to accommodate north-south compression between the Pacific and Gorda plates (Fig. 2). Diagram view inspired by Hole et al. (58).
Supplementary Materials for Low-frequency earthquakes track the evolution of a captured slab fragment
Materials and Methods
Plane fit for LFE alignment
We used a principal component regression to determine the optimal planer fit for the alignment of LFE hypocenters determined by Shelly et al. (27). This regression, which minimizes the sum of squared distances from the hypocenters to the plane, is implemented using an eigenvalue decomposition in MATLAB. This produces an optimal orientation of strike=293º, dip=48º. This is very similar to the visually estimated alignment of strike=290º, dip=45º reported by Shelly et al (27).
Relative plate motions
Relative plate motion from the MORVEL 2010 model (39) were determined using latitude/longitude positions 40.5º/-125º (Gorda) and 40º/-124.5º (Pacific), with the plate motion calculator here: https://www.unavco.org/software/geodetic-utilities/plate-motion-calculator/plate-motion-calculator.html (last accessed May 29, 2025).
Focal mechanisms
To determine the LFE focal mechanism, we utilized stacked seismograms developed by Shelly et al. (27) for each LFE family, derived by stacking the 100 strongest LFE detections for each family from 2018-2024. P and S phases were manually picked from these stacks and used by Shelly et al. (27) to determine source locations.
Although it is often possible to pick the P-wave arrival times on these stacked waveforms, in most cases, due to low signal-to-noise, it is not possible to unambiguously determine the P polarity (the direction of first motion of the P-wave arrival). To overcome this challenge, we used waveforms from 22 nearby earthquakes located near the western edge of the LFE zone which exhibit a variety of focal mechanisms (Figure 1). We then determined relative polarities by the sign of the maximum cross-correlation function, following the approach of Ide et al. (28) and Shelly et al. (29). A positive peak correlation implies the two events likely have the same polarity, with negative correlation implying opposite polarities. We cut earthquake waveforms to 2.5 seconds, beginning 0.3 s prior to the NCSS-listed phase arrival. Earthquake waveforms were filtered 2-8 Hz. We initially aligned earthquake and LFE waveforms on their respective P arrivals, for data at any station with both earthquake and LFE phase arrivals. We allowed a shift of up to +/- 1.0 seconds, recording the time and sign of the peak absolute value correlation coefficient. Observations were weighted based on the difference in amplitude between the absolute values of the larger and second-largest correlation peaks, with the secondary peak separated by at least 0.05s from the primary peak. We required a minimum correlation coefficient of 0.4. Finally, we converted each relative polarity to an absolute polarity based on the NCSS catalog polarity determination and all weighted observations for a given station are added together (multiple channels for a given station and multiple earthquakes correlated with a given LFE family). We repeated this process for each data channel (73 channels from stacked LFE analysis), earthquake (22) and LFE family (27). Vertical and horizontal seismic data channels were used for relative polarity determination.
An example of LFE-earthquake waveform correlations is shown in Fig. S4. Although individual measurements using this approach are often uncertain, increased reliability is attained with redundancy. In total, we used 22 reference earthquakes. Given the 27 LFE families, this provides 594 potential measurements contributing to each single-component station polarity, and 1782 for each 3-component station polarity. Because P-wave arrivals are only visible on a subset of stations for each LFE family, we combined estimated polarities from all families to compute a composite focal mechanism, using the source location of each family. Migration observed among LFE families (27) suggests that they occur within a connected fault zone, supporting the use of a composite mechanism approach.
To compute the final focal mechanism, we used the sign of the weighted summed polarity measurements from these 22 earthquakes to compute a composite focal mechanism solution with SKHASH (77). We determine the source-receiver ray paths needed for the mechanism inversion using the Furlong et al. (17) velocity model assuming a horizontal and vertical hypocentral uncertainty of 0.5 km. To estimate uncertainties, SKHASH computes suites of focal mechanism solutions by varying the input parameters across multiple iterations and finding the solutions that match the observations. For each iteration, the takeoff and source-receiver azimuths are varied by assuming horizontal and vertical hypocentral uncertainties of 0.5 km and the polarities are varied assuming 10% of the first-motion polarities are incorrect.
Earthquake included in this analysis have the following NC network identification numbers: 75009596, 74042641, 74038841, 74038106, 74037811, 74021581, 74016661, 74013766, 73990031, 73985951, 73951815, 73753416, 73709971, 73611556, 73579996, 73299615, 73084215, 72706106, 72570136, 72507356, 72503886, 72499536.
Tidal analysis
We compute time-dependent stresses resulting from the application of Earth’s body and load tides at the locations of LFE hypocenters using SPOTL (78). SPOTL calculates load tides by utilizing the OSU TPXO 7.2 Atlas model (79) which provides the amplitude and phases of these ocean tides at a grid of points globally for the k1, k2, m2, n2, o1, p1, q1, and s2 constituents. For greater accuracy, we also employ the regional osu.usawest.2010 regional model (79) across the west coast of the USA. 3D strains at depth are computed using the Boussinesq solution for a point load on a half space from Malvern (80). The full 3D strain tensor resulting from the load tides is computed at depth. In addition, we use Love numbers and a spherically symmetric Earth model to compute the strains at a specified geographic location. Horizontal strains resulting from the body tides, eθθ and eλλ are computed at the surface and the vertical strain, err, can be expressed as a function of the two horizontal strains as err=-ν/(1-ν)*(eθθ +eλλ) where ν=0.25 is the Poisson’s ratio. We assume this surface strain tensor is a good approximation of the strains resulting from the body tides at depth. We combine the strain tensors for the load and body tides. Assuming a rigidity of 30 GPa and linear elasticity we can convert the time-dependent strain tensor to a time dependent stress tensor.
We first explored the general influence of the tides on the Mendocino LFEs. We computed the time-dependent shear and normal stresses for each orientation by projecting the stress tensor onto the plane defined by its strike, dip, and rake. We found that, like tremor and LFEs elsewhere, their timing is strongly influenced by shear stresses imparted by the solid Earth and ocean tides. To further constrain the slip orientation, we examined variations in the LFE rate as a function of tidally induced shear stress compared to the expected LFE rate if LFEs occurred randomly in time. For each orientation, we computed the ratio of the actual versus the expected number of LFEs as a function of shear stress, assume that distribution is approximately linear and compute the slope, and we also compute the root-mean-square error (Fig. S5).
To narrow the parameter space, we fixed the strike at 293º, following the alignment of LFE hypocenters, and allowed the dip to vary from 35-55º (in increments of 5º) to accommodate greater uncertainty in dip. Testing Coulomb stress with a coefficient of friction 0.1 (38) produced similar results.
Supplementary Figures:
Comparison with active source imaging lines. (A,B) Same as Fig. 3A,B, but with active source imaging lines of Beaudoin et al. (41) (Fig. S3) and Beaudoin et al. (11) (Fig. S2) shown. Their shot points are labeled (black stars) as in the original studies. Dashed black line denotes the slab edge proposed by Beaudoin et al. (11).
Comparison with Beaudoin et al. (11) seismic refraction data (Line 9). Shotpoint locations are shown in Fig. S1. Cross-section c-c’ (Figure 3E) is approximately overlain, shifted ~3 km deeper to roughly account for slab dip and the more easterly location of line 9. Upper panel shows the velocity model of Beaudoin et al. (11). Although the low-velocity zone observed in tomography studies (17, 31) is not apparent here, the velocity contrast seen in tomography is greatly reduced at this more easterly location (Fig. S1). Furthermore, the refraction ray density (lower panel) is low in this area, suggesting limited resolution in that model. The overall velocities reported in Beaudoin et al. (11) are generally lower than those reported from tomography by Furlong et al. (17) (Fig. 3) at 25-30 km depth; however, a north dipping high velocity zone can be in the refraction model near shot point 905 (see dipping 7 km/s contour), which is qualitatively similar to the tomography model.
Comparison of cross-section a-a’ (Figure 3) with seismic velocity results. (A) Seismic velocity model from Beaudoin et al. (41). See Fig. S1 for location of shot points. Ray paths for shot point 4 are shown. Cross-section a-a’ from Fig. 3C is approximately overlain on these results. Note the 2:1 vertical exaggeration of the Beaudoin et al. (41) cross-section. The a-a’ cross section lies on the very edge of this model. The simple velocity model on the southern edge is likely reflective of limited resolution in the Beaudoin et al. (41) model here. (B) Ratio of P-wave to S-wave velocities (Vp/Vs) along line a-a’ (Fig. 3) as determined by Guo et al. (42). Model is partially masked in areas where Guo et al. (42) reported limited resolution. Plotted seismicity same as Fig. 3C. Although Guo et al. (42) inferred a horizontal upper boundary of the Gorda along a similar cross-section, the dipping boundary proposed here is also consistent with high Vp/Vs inferred by Guo et al. (42) to reflect the Gorda crust.
Example of relative polarity analysis, showing correlation of P waveforms for one earthquake (ID=nc74042641, 2.5 s duration, red or blue) with longer duration stacked waveforms for one LFE family (family #8, black). Station and channel names are given at right. At left is signed cross-correlation (cc) coefficient, and the difference between the peak and secondary absolute value correlation coefficient (ccdiff). Positive correlation indicates that the earthquake and LFE family likely have the same polarity at that station. Negative correlation indicates that they are likely opposite. ccdiff is a measure of confidence in this relative polarity. Here it is signed by estimated polarity (earthquake polarity times relative polarity), where blue text and waveform represents likely compressional LFE motion and red represents dilatational. Orange and purple bars indicate P and S picks, respectively for the LFE family from Shelly et al. (27).
Tidal modulation as a function of dip (separate panels), strike, and rake. Color indicates slope, with blue indicating positive slope, while red indicates negative slope. Circle size scales with linearity of the actual/expected numbers of events versus tidal shear stress. The orientations most consistent with tidal stresses are generally those with greatest (positive) tidal modulation slope and at least moderate linearity. The most plausible parameter space is explored in more detail in Fig. 2.
Authors
Disclaimers
Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
Although this information product, for the most part, is in the public domain, it also may contain copyrighted materials as noted in the text. Permission to reproduce copyrighted items must be secured from the copyright owner.
Suggested Citation
Shelly, D.R., Thomas, A.M., Materna, K.Z., and Skoumal, R.J., 2026, Low-frequency earthquakes track the motion of a captured slab fragment: Science, v. 391, no. 6782, p. 294-299, https://doi.org/10.1126/science.aeb2407.
| Publication type | Article |
|---|---|
| Publication Subtype | Journal Article |
| Title | Low-frequency earthquakes track the motion of a captured slab fragment |
| Series title | Science |
| DOI | 10.1126/science.aeb2407 |
| Volume | 391 |
| Issue | 6782 |
| Publication Date | January 15, 2026 |
| Year Published | 2026 |
| Language | English |
| Publisher | American Association for the Advancement of Science |
| Contributing office(s) | Geologic Hazards Science Center - Seismology / Geomagnetism |
| Description | 6 p. |
| First page | 294 |
| Last page | 299 |