Unsupervised machine learning reveals slab hydration variations from deep earthquake distributions beneath the northwest Pacific

Deep-focus earthquakes (>300 km depth) occur under high-pressure and high-temperature conditions where brittle failure should be inhibited, making their mechanisms enigmatic. Three main hypotheses have been proposed: thermal shear instability, dehydration embrittlement, and transformational faulting associated with the olivine (α) to wadsleyite (β-spinel) phase transition. Recent observations and experiments increasingly support transformational faulting within a metastable olivine wedge (MOW), particularly near its rim, but direct seismic constraints on the rim’s geometry are lacking due to detection limits for small events. The Gutenberg–Richter b-value (slope of frequency–magnitude distribution) offers insight into rupture processes and has been linked to slab thermal state and mechanisms for deep earthquakes. Prior work suggested nucleation within the MOW and possible propagation outside via thermal runaway for larger events, yet small-magnitude completeness has limited constraints on nucleation. This study investigates whether variations in b-values across magnitude ranges and regions can reveal the presence and thickness of an unstable MOW rim and assess the role of slab hydration in modulating deep earthquake mechanisms in the northwestern Pacific.

Past studies identified metastable olivine wedges in subducting slabs and proposed transformational faulting as a mechanism for deep earthquakes, with slab thermal parameters controlling MOW geometry. While dehydration embrittlement has been invoked for intermediate-depth seismicity, its applicability to deep-focus events is debated, particularly regarding mantle earthquakes. Hydrous defects can catalyze phase transformations and potentially reduce olivine metastability, implying slab hydration also shapes MOW structure. Observed b-value variations correlate with slab thermal state, with proposed changes at larger magnitudes indicating shifts from transformational faulting to thermal shear instability and/or fault width saturation. However, global catalogs have high completeness thresholds, limiting inferences at small magnitudes critical for nucleation mechanisms.

Data: The study analyzed 10,660 earthquakes at depths ≥300 km from the Japan Meteorological Agency (JMA) catalog (2002–2016) across Kuril, Japan, and Izu–Bonin subduction zones. Additional analyses revisited the global CMT catalog (with three more years than prior work) and the IRIS catalog to examine larger magnitude behavior. Depth range selection (300–700 km) is detailed in Supplementary Note 3.

Clustering: Unsupervised K-means clustering was applied to earthquake hypocenter locations to objectively partition deep events. The optimal number of clusters (four) was determined via silhouette analysis, producing spatial groupings corresponding to Bonin, Izu, Honshu, and Kuril slab segments. Robustness checks using spectral clustering and Gaussian mixture models yielded similar partitions (Supplementary Notes 5–7).

Magnitude processing and completeness: JMA magnitudes (M_J) were converted to moment magnitude (Mw) using an empirical relation M = a + b M_J + c with a = 0.053 ± 0.003, b = 0.33 ± 0.02, c = 1.68 ± 0.03. Magnitude of completeness (Mc) was estimated from non-cumulative (binned) distributions (Supplementary Note 8), yielding Mc ≈ 3.3 for Bonin (sparser stations) and ≈3.0 for the other clusters.

b-value estimation: For each cluster, b-values were computed above Mc using three estimators: (1) maximum likelihood (b0), (2) modified maximum likelihood correcting for binning (b1), and (3) Tinti–Mulargia estimator less sensitive to magnitude errors (b2). Corresponding standard deviations were computed for uncertainty characterization. To test robustness, Monte Carlo bootstrapping resampled events within clusters to assess the stability and potential bimodality of b-value distributions. Sensitivity tests confirmed that observed kinks were not artifacts of specific large events (e.g., Tohoku-Oki) affecting stress fields.

Cross-catalog comparison at large magnitudes: Gutenberg–Richter distributions were also computed for deep earthquakes in global “warm” subduction zones (Japan–Kuril, Izu–Bonin–Mariana, South America, Philippine) and Tonga (cold) using CMT and IRIS catalogs to evaluate kinks at larger Mw and relate to MOW thickness.

• Four deep-earthquake clusters (Bonin, Izu, Honshu, Kuril) were identified via K-means, each corresponding to distinct slab segments and inferred hydration states. Event counts: Kuril 1643, Honshu 5338, Izu 2108, Bonin 1571.
• Honshu and Izu clusters exhibit a distinct Gutenberg–Richter kink at Mw ≈ 3.7–3.8 for depths >300 km. Below the kink, b ≈ 1.4 (Honshu) and 1.7 (Izu); above the kink, b ≈ 0.6 (Honshu) and 0.7 (Izu). The change is far larger than standard deviations (<0.03), and bootstrapping reveals bimodal b-value distributions supporting the kink.
• Kuril and Bonin clusters show nearly linear GR relations (no kink) with b ≈ 1.0 across analyzed magnitudes.
• The high b-values at small magnitudes (Mw < 3.7–3.8) in Honshu and Izu indicate enhanced small-event nucleation attributed to transformational faulting in a narrow, unstable rim of the MOW, catalyzed by hydrous defects.
• The characteristic size of Mw 3.7–3.8 events implies a seismogenic thickness ≤1 km, constraining the unstable MOW rim thickness to ~1 km (between ≈1 and ≈4 km). Larger events with Mw > 3.7–3.8 can rupture inward from the rim into the more homogeneous MOW interior with characteristic thickness of 10–20 km.
• Regional tectonic context explains hydration contrasts: Honshu and Izu feature obliquely oriented pre-existing fabrics and outer-rise bending faults that enhance deep hydration; Kuril’s plate fabric orientation limits new bending faults and hydration; Bonin’s Cretaceous–Paleogene magmatic overprint (Ogasawara Plateau) implies a warmer, drier lithosphere and reduced deep hydration.
• Global “warm” subduction zones (CMT) show a kink at Mw ≈ 6.7, with b approximately doubling from ~0.6 below to ~1.1 above; Tonga (cold) shows no significant large-magnitude kink and maintains b ≈ 1.1, consistent with a thicker MOW accommodating large ruptures without mechanism change.
• Table 1 b-values (three estimators) corroborate: Kuril ~1.0; Bonin ~1.0; Izu and Honshu high b at Mw <3.8 (Izu: 1.7–2.1 depending on estimator; Honshu: 1.4–1.7) and low b at Mw >3.8 (Izu: ~0.73–0.80; Honshu: ~0.60–0.65).

The observed b-value kink at Mw 3.7–3.8 in Honshu and Izu directly addresses the study’s aim to constrain rupture domains and mechanisms at small magnitudes in deep-focus seismicity. High b for small events implies strong heterogeneity and dense nucleation facilitated by hydrous defects within a narrow unstable MOW rim, rather than high pore pressure. The abrupt drop in b above the kink signifies a transition from primary transformational faulting in the rim to ruptures propagating into the more homogeneous interior of the MOW via transformation-driven stress transfer (TDST), yielding lower fractal dimension and b-values. The absence of a small-magnitude kink in Kuril and Bonin aligns with drier slabs and reduced hydrous defects, limiting small-event nucleation in the rim. At larger magnitudes (global catalogs), the kink near Mw ~6.7 and increased b-values suggest additional rupture propagation beyond the MOW into warmer wadsleyite via thermal runaway, introducing new heterogeneity and changing controlling mechanisms. Together, the results support a three-domain framework: (1) unstable MOW rim (transformational faulting, high b at small Mw), (2) MOW interior (TDST, lower b), and (3) surrounding wadsleyite (thermal runaway for largest events, increased b again). Spatial variations in slab hydration and thermal state modulate these domains and b-value behaviors across subduction zones. The findings provide seismic constraints on the MOW rim thickness (~1 km) and MOW interior thickness (10–20 km), linking deep earthquake size distributions to slab structure and hydration history.

Using an unsupervised clustering of a highly complete regional catalog (JMA, 2002–2016), the study identifies four deep-earthquake clusters in the northwestern Pacific and reveals a robust b-value kink at Mw 3.7–3.8 in the Honshu and Izu segments. Below this threshold, b is anomalously high (1.4–1.7), transitioning to low b (0.6–0.7) above, indicating a shift from nucleation within a narrow unstable rim of the metastable olivine wedge (MOW) to rupture propagation into the MOW interior. This behavior is absent in the drier Kuril and Bonin segments (b ~1), tying the anomaly to slab hydration. The results provide a seismic constraint on the MOW rim thickness (~1 km) and, together with global catalog analyses, support a three-stage mechanism with thermal runaway enabling rupture beyond the MOW at large magnitudes (kink near Mw ~6.7). These insights link deep earthquake statistics to slab hydration and thermal structure. Future work should: (1) extend small-magnitude completeness to other subduction zones (e.g., Tonga) to test for low-Mw kinks; (2) integrate high-resolution slab structure imaging to refine MOW and rim geometry; (3) improve magnitude calibration and completeness assessments; and (4) develop models coupling hydrous defect distributions with transformational kinetics and rupture dynamics.

• Detection limits and completeness: Despite low Mc (~3.0–3.3) in JMA, smaller events remain undetected, potentially biasing the low-Mw GR slope. Bonin has fewer nearby stations (higher Mc), affecting comparisons.
• Temporal coverage: JMA analysis spans 2002–2016, limiting the number of large deep events and precluding observation of high-Mw kinks within this dataset; large-Mw inferences rely on global catalogs.
• Magnitude conversion and uncertainties: Mw derived from M_J via empirical relations introduces magnitude errors that can affect b-value estimates; multiple estimators were used to mitigate but not eliminate this issue.
• Regional generalizability: Findings for low-Mw kinks are robust for Honshu and Izu but require more complete catalogs to verify in other slabs (e.g., cold Tonga). Structural interpretations (hydration, MOW thickness) are inferred from statistical patterns and tectonic context rather than direct imaging.
• Clustering assumptions: K-means uses spatial hypocenter distributions and may be sensitive to cluster number selection, though validated with silhouette scores and alternative algorithms.