Subslab ultra low velocity anomaly uncovered by and facilitating the largest deep earthquake

Deep-focus earthquakes (>300 km) remain difficult to explain given high-pressure–temperature conditions that inhibit rapid seismic slip. A prevailing model invokes transformational faulting of olivine to spinel within the metastable olivine wedge (MOW), which helps explain the cutoff of deep seismicity near the 660 km discontinuity (660-D) but struggles to account for rare M8+ events unless ruptures are highly elongated or experience strong dynamic weakening that allows propagation beyond the MOW. The two-stage hypothesis posits initial transformational faulting followed by thermal runaway instability (TRI), which requires high strain/stress accumulation within the slab. The study aims to resolve how such large stresses arise and how a small rupture within the MOW can evolve into an M8+ event. To address this, the authors apply high-frequency (up to 0.8 Hz) teleseismic waveform modeling to image small-scale P-wave velocity anomalies near 660-D at the Sea of Okhotsk and to analyze the earliest subevents of the 2013 Mw 8.3 Okhotsk earthquake. The central hypothesis is that a nearby, ultra–low-velocity, small-scale anomaly beneath the slab provides buoyancy-induced stress that helps trigger TRI and enables rupture to propagate outside the MOW.

Foundational work links deep earthquakes to metastable mantle phase transformations within subducting slabs and to plastic/thermal instabilities (e.g., MOW, TRI). Prior source studies of the 1994 Bolivia and 2013 Okhotsk deep earthquakes revealed rupture complexity and multiple subevents with sub-horizontal faulting and strong directivity. High-frequency waveform modeling has been shown to better resolve sharp velocity anomalies than travel-time tomography, which tends to smooth features and underestimates amplitude (often by a factor of 1.5–2). Global and regional tomography commonly shows low-velocity anomalies beneath slabs in the mantle transition zone (MTZ), possibly linked to plumes or entrained upper mantle material, implying localized high temperatures conducive to partial melting and volatile-rich anomalies. However, the specific small-scale structures hosting extreme stress and their role in transitioning from initial rupture to large M8+ events remained unresolved.

Data and preprocessing: The study collects global teleseismic P-wave recordings at 30–90° epicentral distances and bins stations by azimuth (e.g., China/SE Asia, Europe/Central Asia, USArray West/East, Pacific). Initial quality control and picks are refined via cross-correlation. Records are downsampled to 0.1 s. One-dimensional synthetics are computed using PREM and the gCMT focal mechanism (1 s source time function), then deconvolved from observations using Projected Landweber Deconvolution (PLD) to obtain relative source time functions (RSTFs), with a deconvolution window of 2 s before to 15 s after P. RSTFs are stacked within 2° latitude × 3° longitude boxes and bandpass filtered (0.1–0.5 Hz; 0.2–0.8 Hz). A path calibration approach uses a smaller nearby event to identify clean receiver-side paths.

Subevent selection and relative relocation: The analysis focuses on the initial subevent S1 and an early subevent S2 within E1; a smaller Sx is also identified between S1 and S2. Arrival-time differences between subevents are modeled as T = τ + Γ + LΔt + H, where τ is origin time difference, Γ the horizontal directivity parameter, Δt the vertical directivity parameter, and L and H the horizontal/vertical offsets. A grid search minimizes misfit over τ, L, H, and rupture direction Θ using observations across arrays to constrain relative timing and locations.

Multipathing detection: The USArray observations reveal in-plane multipathing for S2 in specific azimuth–distance bins (notably WA3–WA4), manifested as a ~1 s jump in peak arrivals at ~55–60° in distance profiles and absent in azimuthal profiles—characteristic of in-plane multipathing. A Multipathing Detector (MPD) quantifies pulse splitting time and amplitude ratios of split arrivals on 0.2–0.8 Hz RSTFs, mapping coherent splitting primarily for S2 toward North America.

Forward waveform modeling: To isolate the cause of multipathing, 2-D velocity models combining Slab (S), background Plume (P), and a Subslab Ultra–Low-Velocity Anomaly (U, SULVA) are tested. Only models including a small-scale SULVA reproduce the observed S2 multipathing; slab and plume features alone are insufficient. Given the dominance of in-plane effects, 2-D finite-difference simulations (CUDA/GPU-accelerated, 1 km grid; minimum wavelength ~7.25 km; numerical accuracy ~1.5 Hz) generate synthetic seismograms that are deconvolved with the same 1-D Green’s functions to produce synthetic RSTFs. Stations from paired, opposite azimuth bins (e.g., WA3 and conjugate CN bin) allow a single 2-D model to match both directions.

SULVA parameterization and fitting: The SULVA is parameterized as a “cut-donut” with inner radius r, outer radius R (measured from S1), and two sharp dipping boundaries (cut1 at the left boundary from S1; cut2 at the right boundary from S2). A uniform P-wave velocity reduction dVp is assumed. Model evaluation uses relative S1–S2 timing and waveform similarity in key bins (especially WA3). A Template Matching Cross-Correlation (TMCC) metric, computed over a 2.5 s S2 window and summed across traces, quantifies fits. Sensitivity analysis shows strong dependence on cut2 dip and dVp; changes >5° in cut2 or ~4% in dVp significantly alter the multipathing pattern. Parameter uncertainties are estimated by bootstrapping (80% station resampling, 2000 iterations) to derive 95% confidence intervals. A 2-D grid search constrains r and R. The approach is extended across multiple azimuthal bins (EA2, WA1, WA2, WA4) to map SULVA geometry.

Buoyancy stress estimation: Using the best-fit SULVA dimensions intersecting the 660-D, buoyancy-induced stress is estimated as σ = Σ Δρ_i g h_i for portions above and below 660-D. Numerical examples consider hydrous silicate melt (e.g., (Mg0.75,Fe0.25)2SiO4 with 4 wt% H2O) and a pure-volatile endmember, with and without a 10 km uplift of 660-D, to quantify stress magnitudes (Methods). Robustness tests examine alternative explanations (anisotropy, attenuation, Vp/Vs variations, out-of-plane multipathing, focal mechanism changes, and 660-D topography), none of which reproduce the observed multipathing except in extreme cases.

• The earliest rupture of the 2013 Mw 8.3 Okhotsk event comprises distinct subevents: S1 (initial), Sx (smaller, predominantly horizontally offset), and S2 (early in E1, shallower than S1 with significant vertical offset). S2’s positive distance move-out indicates larger ray parameter and shallower depth relative to S1.
• Coherent in-plane multipathing is observed for S2 in North America (notably USArray WA3–WA4), evident as a ~1 s jump in peak arrivals around 55–60° in distance profiles at 0.2–0.8 Hz; S1 shows stable waveforms without multipathing. The effect is strongest in distance profiles, not azimuthal profiles, indicating a source-side, small-scale structure as the cause.
• Existing tomography reveals only a weak (−2 to −3%) low-velocity anomaly beneath the slab, insufficient to cause the observed multipathing. Forward modeling demonstrates that only inclusion of a small-scale Subslab Ultra–Low-Velocity Anomaly (SULVA) reproduces the observations; slab or plume structures alone do not.
• Best-fit SULVA properties from 2-D finite-difference waveform modeling and TMCC fitting across multiple azimuth bins: P-wave velocity reduction dVp = 18 ± 2%; sharp right boundary dip cut2 = 40 ± 2° (strong sensitivity); horizontal dimensions ~30 × 60 km; vertical thickness ~60 km; geometry tapers by ~50% from SE to NW; the anomaly spans across the 660-D (present both above and below).
• The multipathing arises from splitting of S2 energy along the SULVA boundary; S1 likely penetrates the SULVA center and thus does not split.
• Buoyancy stress estimates from SULVA density contrast across 660-D: for hydrous silicate melt (4 wt% H2O), σ ≈ 86 MPa (flat 660-D) to ≈125 MPa (10 km uplift). For a volatile-only endmember (e.g., diamond density), σ ≈ 216–255 MPa depending on 660-D uplift.
• Comparative context: Low-velocity anomalies beneath slabs are also imaged in global/regional tomography near other M8+ deep earthquakes (1994 Bolivia; 2018 Fiji), and slab distortion near the 2018 Fiji event is consistent with buoyancy effects, suggesting a common setting.
• Interpretation: Volatile-bearing, partially molten SULVA provides significant buoyancy, adding stress that pushes the slab toward conditions for thermal runaway instability. This facilitates rupture propagation beyond the narrow MOW, enabling M8+ deep earthquakes. The thermal/velocity structure of mantle near 660-D, rather than intrinsic slab properties, may chiefly determine the occurrence of such large deep events.

The study directly addresses how a deep earthquake can grow to M8+ despite the restricted MOW by showing that a small-scale, ultra–low-velocity, volatile-bearing anomaly beneath the slab near 660-D can produce substantial buoyancy stresses. The observation of in-plane multipathing exclusively for the shallower S2, and not for S1, ties the effect to source-adjacent structure. Forward modeling constrains an ~18–20% P-wave velocity reduction within a ~30 × 60 × 60 km SULVA that spans 660-D. Estimated buoyancy stresses (≈86–125 MPa for hydrous melt; up to ≈255 MPa for volatile-rich endmembers) are of the order needed, in conjunction with realistic rheologies, to facilitate TRI. This supports a dual-mechanism rupture model: transformational faulting initiates rupture in the MOW (S1–S2), and buoyancy-enhanced stresses from the SULVA enable thermal runaway so the rupture extends beyond the MOW, producing large deep events. Tests argue against alternative explanations (e.g., anisotropy, attenuation, 660-D topography) for the multipathing, bolstering the SULVA interpretation. The presence of similar sub-slab low-velocity anomalies near other M8+ deep earthquakes implies that interactions between water-rich melts and slabs near 660-D may be a widespread prerequisite for generating such large deep events.

High-frequency teleseismic waveform analysis of the 2013 Mw 8.3 Sea of Okhotsk deep earthquake reveals a small-scale (~30 × 60 × 60 km), ultra–low-velocity (dVp ≈ 18–20%) anomaly beneath the slab that straddles the 660 km discontinuity. This SULVA produces in-plane multipathing of early subevent S2 and likely reflects volatile-bearing partial melt. Its buoyancy imparts stresses of order 10^2 MPa, bringing the slab closer to thermal runaway conditions and enabling rupture to propagate beyond the MOW—offering a mechanism for M8+ deep earthquakes. Similar sub-slab low-velocity anomalies near other M8+ events suggest the mechanism may be broadly applicable. Future work should combine higher-frequency waveform modeling with improved tomography, and pursue mineral physics and geodynamic studies to explain the large velocity reductions and constrain the composition, melt fraction, water content, and dynamics of SULVAs near 660-D.

• Model dimensionality: Forward modeling uses 2-D FD simulations emphasizing in-plane multipathing; potential out-of-plane effects are not fully captured, though tests suggest limited impact for observed features.
• Structural non-uniqueness and parameter trade-offs: The SULVA inner radius r is less well constrained; some trade-offs exist between r and R. Inferences depend on the chosen parameterization (cut-donut) and uniform dVp assumption.
• Frequency band and deconvolution: Analysis relies on RSTFs in 0.1–0.5 and 0.2–0.8 Hz bands with PLD deconvolution; precise wiggle-by-wiggle modeling at higher frequencies remains challenging.
• Background models and simplifying assumptions: Use of 1-D PREM Green’s functions and simplified velocity anomalies may omit fine 3-D heterogeneity; attenuation, anisotropy, and Vp/Vs variations were tested but cannot be ruled out in extreme scenarios.
• Interface/topography uncertainties: 660-D topography can shift arrivals uniformly but cannot reproduce multipathing; its detailed structure remains uncertain.
• Event complexity: Later, larger subevents (E2–E4) are longer/complex and not modeled at high frequency due to coda contamination; Sx may have a different focal mechanism, adding uncertainty to relative relocation.
• Composition and density: Buoyancy stress estimates depend on assumed melt composition and water content; exact SULVA composition and melt fraction are not directly measured.