Geophysical Journal International

SummaryIn this work, we study the induced seismicity recorded during an injection operation at the Muara Laboh geothermal plant (Indonesia). The swarm, consisting of three bursts activating a normal fault zone, is characterized by rapid earthquake (km/day) migration. We use a 2D rate-and-state asperity model to better understand the physical mechanisms controlling the evolution of this induced swarm. The model suggests that the observed rapid seismic migration can be explained by the interaction among asperities through the expansion of slow postseismic slip fronts. Also, it shows that the amount of seismicity generated by the fluid injection is strongly controlled by the background seismicity of the system, that is by the seismicity determined by the tectonic load charging the fault. This close correlation between natural and induced seismicity suggests that the injection in Muara Laboh principally stimulates critically stressed faults, which release the seismicity determined by their natural seismic cycle.

SummaryThe Mw 6.8 Murghob earthquake is the third earthquake in an Mw 6.4+ sequence occurring in the Pamir initiated by the 2015 Sarez Mw 7.2 earthquake. It is of great significance to investigate their interactions and to assess future seismic hazards in the region. In this paper, we use Sentinel-1 radar interferometric data to retrieve coseismic deformation, invert for the slip distributions of the four events, and then investigate their interactions. The cumulative Coulomb failure stress changes (ΔCFS) suggest that the 2023 Murghob earthquake was promoted by the three prior earthquakes in the sequence. Pre-stress from historical earthquakes is a key factor in explaining the triggering mechanism of the two 2016 Mw 6.4+ earthquakes. Stress loading and unloading effects on major faults in the region indicate that future attention should be paid in 1) the segment of the Sarez-Karakul fault north of the Kokuibel Valley, 2) the segment of the Sarez-Murghab thrust fault west of the Sarez-Karakul fault, and 3) the east segments of the Pamir thrust fault system, all with a large positive ΔCFS.

SummaryThe Lijiang–Xiaojinhe fault (LXF) and its vicinity are located in the transition zone among the Tibetan Plateau (TP), the South China block and the Indochina block. Researchers believe that this area has acted as a key tectonic zone during the evolution of the TP. Owing to the continuous growth and SE-ward expansion of the TP, the LXF and its vicinity have experienced intense deformation. Although different models, such as the rigid block extrusion and mid-lower crustal flow models, have been proposed to explain this intense deformation, a consensus has not yet been achieved. To better understand the deformation of the LXF and its vicinity, a high-resolution image of the subsurface structure must be constructed. In this study, we construct images of wave velocity and azimuthal anisotropy structures by using an eikonal equation-based traveltime tomography method. We collect high-quality seismic data from 276 broadband seismic stations and manually pick a total of 48,037 first arrivals for the tomography study. Our tomographic results reveal a strong low-velocity body below the LXF and its vicinity. In addition, a strong azimuthal anisotropy structure with a N–S-oriented fast velocity direction (FVDs) is distributed along the low-velocity body. These features indicate the occurrence of mid-lower crustal flow, that penetrates across the LXF and extends to the Dianzhong block (DZB). In addition, we find obvious low-velocity perturbations in the mid-lower crust and uppermost mantle beneath the DZB. The low-velocities may be attributed to the upwelling of hot materials from the upper mantle. We consider the limited distribution of mid-lower crustal flow on the margin of the SE TP, and mid-lower crustal flow may not play a significant role in the expansion of the TP.

SummaryThe recent rapid improvement of machine learning techniques had a large impact on the way seismological data can be processed. During the last years several machine learning algorithms determining seismic onset times have been published facilitating the automatic picking of large data sets. Here we apply the deep neural network PhaseNet to a network of over 900 permanent and temporal broad band stations that were deployed as part of the AlpArray research initiative in the Greater Alpine Region (GAR) during 2016-2020. We selected 384 well distributed earthquakes with ML ≥ 2.5 for our study and developed a purely data-driven pre-inversion pick selection method to consistently remove outliers from the automatic pick catalog. This allows us to include observations throughout the crustal triplication zone resulting in 39,599 P and 13,188 S observations. Using the established VELEST and the recently developed McMC codes we invert for the 1D P- and S-wave velocity structure including station correction terms while simultaneously relocating the events. As a result we present two separate models differing in the maximum included observation distance and therefore their suggested usage. The model AlpsLocPS is based on arrivals from ≤ 130 km and therefore should be used to consistently (re)-locate seismicity based on P & S observations. The model GAR1D_PS includes the entire observable distance range of up to 1000 km and for the first time provides consistent P- & S-phase synthetic travel times for the entire Alpine orogen. Comparing our relocated seismicity with hypocentral parameters from other studies in the area we quantify the absolute horizontal and vertical accuracy of event locations as ≈ 2.0 km and ≈ 6.0 km, respectively.

AbstractCombining elastic full waveform inversion (FWI) with rock physics holds promise for expanding the application of FWI beyond seismic imaging to predicting and monitoring reservoir properties. Distributed Acoustic Sensing (DAS), a rapidly developing seismic acquisition technology, is being explored for its potential in supporting FWI applications. In this study, we implement a sequential inversion scheme that integrates elastic FWI and Bayesian rock physics inversion, using a vertical seismic profile (VSP) dataset acquired with accelerometer and collocated DAS fiber at the Carbon Management Canada’s Newell County Facility. Our aim is to establish a baseline model of porosity and lithology parameters to support later monitoring of CO2 storage. Key strategies include an effective source approach for addressing near-surface complications, a modeling strategy to simulate DAS data comparable to field data, and a Gaussian mixture approach to capture the bimodality of rock properties. We conduct FWI tests on accelerometer, DAS, and combined accelerometer-DAS data. While our inversion results accurately reproduce either dataset, the elastic models inverted from accelerometer data outperform the other two in matching well logs and identifying the target reservoir. We attribute this outcome to the limited complementarity of DAS data with accelerometer data in our experiment, along with the limitations imposed by single-component measurements on DAS. The porosity and lithology models predicted from accelerometer-derived elastic models are reasonably accurate at the well location and exhibit geologically meaningful spatial distribution.

SummaryFractured hydrate-bearing reservoirs are extensively discovered worldwide and show notable anisotropic geophysical properties. Hydrate distribution in fractures significantly affects the anisotropic properties, and hence plays an important role in the accurate assessment of hydrate resources. However, the knowledge about how the hydrate distribution impacts the anisotropic geophysical properties of fractured reservoirs, which is the premise for the identification and quantification of hydrate in fractured reservoirs, is still poorly understood. To obtain such knowledge, we forward study the effects of various hydrate distribution (including floating, bridging and evolving hydrate distribution) in aligned fractures on the anisotropic elastic, electrical and joint elastic-electrical properties of a digital core using validated numerical methods. We show that for all the hydrate distribution, the anisotropic velocities increase while the conductivities decrease with increasing hydrate saturation, with the effects of the floating and bridging distribution being the least and greatest, respectively. We also show that the anisotropic velocities and conductivities for the floating and bridging distribution vary approximately linearly with hydrate saturation, leading to linear correlations between the elastic and electrical properties. Further investigation illustrates that the difference in the slopes of the linear joint correlations between the two distribution is significantly greater than that of the individual properties as a function of hydrate saturation. The results have revealed the distinct effects of hydrate distribution on the anisotropic elastic and electrical properties of fractured reservoirs, and have confirmed the superiority of the joint elastic-electrical properties for the distinguishment of hydrate distribution in fractures over individual physical properties.

SummaryPrevious studies have shown that it is difficult to determine whether the 2015 Pishan earthquake occurred on a uniform fault or a ramp-flat fault with variable dip angles due to the similar goodness of data fit to co-seismic and afterslip models on these two fault models. Here, we first present the InSAR deformation obtained from both ascending and descending orbits, covering the coseismic period and cumulative five-year period after the 2015 Pishan earthquake. We then determine the preferred fault geometry by the spatial distributions between the positive Coulomb failure stress change triggered by mainshock and the afterslip. Based on the preferred fault model, we finally use a combined model to determine the contributions of elastic and viscoelastic deformation in the post-seismic deformation. We find that the Pishan earthquake prefers to occur on a ramp-flat fault, and the coseismic slip is mainly distributed at a depth of 9-13 km, with a maximum slip of about 1.3 m. The post-seismic deformation is primarily governed by afterslip, as the poroelastic rebound-induced deformation fails to account for the observed post-seismic deformation and the contributions from the viscoelastic relaxation mechanism can be considered negligible in the combined model. Moreover, the modelled stress-driven afterslip and observed kinematic afterslip have good consistency, and the difference between the root mean square error (RMSE) of the two afterslip models is only 4.3 mm. The results from the afterslip model indicate that both of the up- and down-dip directions distribute the afterslip, and slip in the up-dip direction is greater than that of the down-dip direction. Meanwhile, the maximum cumulative afterslip after five years is approximately 0.26 m which is equivalent to a released seismic moment of a Mw 6.47.

SummaryThe Earth’s magnetic field is generated by a dynamo in the outer core and is crucial for shielding our planet from harmful radiation. Despite the established importance of the core-mantle boundary heat flux as driver for the dynamo, open questions remain about how heat flux heterogeneities affect the magnetic field. Here, we explore the distribution of core-mantle boundary heat flux on Earth and its changes over time using compressible global 3-D mantle convection models in the geodynamic modeling software ASPECT. We discuss the use of the consistent boundary flux method as a tool to more accurately compute boundary heat fluxes in finite element simulations and the workflow to provide the computed heat flux patterns as boundary conditions in geodynamo simulations. Our models use a plate reconstruction throughout the last 1 billion years—encompassing the complete supercontinent cycle—to determine the location and sinking speed of subducted plates. The results show how mantle upwellings and downwellings create localized heat flux anomalies at the core-mantle boundary that can vary drastically over Earth’s history and depend on the properties and evolution of the lowermost mantle as well as the surface subduction zone configuration. The distribution of hot and cold structures at the core-mantle boundary changes throughout the supercontinent cycle in terms of location, shape and number, indicating that these structures fluctuate and might have looked very differently in Earth’s past. We estimate the resulting amplitude of spatial heat flux variations, expressed by the ratio of peak-to-peak amplitude to average heat flux, q*, to be at least 2. However, depending on the material properties and the adiabatic heat flux out of the core, q* can easily reach values >30. For a given set of material properties, q* generally varies by 30-50% over time. Our results have implications for understanding the Earth’s thermal evolution and the stability of its magnetic field over geological timescales. They provide insights into the potential effects of the mantle on the magnetic field and pave the way for further exploring questions about the nucleation of the inner core and the past state of the lowermost mantle.

SummaryThe wide-swath altimeter satellite Surface Water and Ocean Topography (SWOT) will provide high spatiotemporal resolution sea surface heights (SSHs), which is crucial for studying the impact of observation errors on marine gravity recovery. This study uses simulated SWOT data to derive deflection of the vertical (DOV) and gravity anomalies in the northern South China Sea. We quantified the impact of SWOT errors on DOV and gravity anomalies, and analyzed the contributions from different directions of geoid gradient. The results show that the geoid gradient in the cross-track direction significantly improves gravity field recovery by enhancing the precision of east component of DOV. For one-cycle SWOT observations, phase errors emerge as the most impactful error affecting both DOV and gravity anomalies, followed by random errors. Two-dimensional Gaussian filtering and the tilt correction proposed in this study could effectively mitigate their impact. Using the corrected data for DOV computation, the precision in the east and north components improves by 75.32% and 46.80%, respectively, while enhancing the accuracy of the gravity field by 70.23%. For 17-cycle data, phase errors and random errors remain the predominant factors affecting DOV and gravity anomalies, but their impact diminishes with an increase in SWOT observations. Our results indicate that marine gravity accuracy improves by approximately 70% compared to a single cycle. Whether for single-cycle or multi-cycle data, the impact of phase errors is roughly twice that of random errors. These data processing strategies can serve as valuable references for wide-swath altimeter data processing, aiming to advance the precision and resolution of marine gravity field recovery.

SummaryAround the world the Earth’s crust is blanketed to various extents by sediment. For continental regions, knowledge of the distribution and thickness of sediments is crucial for a wide range of applications including seismic hazard, resource potential, and our ability to constrain the deeper crustal geology. Excellent constraints on the sediment thickness can be obtained from borehole drilling or active seismic surveys. However, these approaches are expensive and impractical in remote continental interiors such as central Australia. Recently, a method for estimating the sediment thickness using passive seismic data, the collection of which is relatively simple and low-cost, was developed and applied to seismic stations in South Australia. This method uses receiver functions, specifically the time delay of the P-to-S converted phase generated at the sediment-basement interface, relative to the direct-P arrival, to generate a first order estimate of the thickness of sediments. In this work we expand the analysis to the vast array of over 1500 seismic stations across Australia, covering an entire continent and numerous sedimentary basins that span the entire range from Precambrian to present-day. We compare with an established yet separate method to estimate the sediment thickness, which utilises the autocorrelation of the radial receiver functions to ascertain the two-way travel-time of shear waves reverberating in a sedimentary layer. Across the Australian continent the new results match the broad pattern of expected sedimentary features based on the various geological provinces. We are able to delineate the boundaries of many sedimentary basins, such as the Eucla and Murray Basins, which are Cenozoic, and the boundary between the Karumba Basin and the mineral rich Mount Isa Province. Contrasts in seismic delay time across these boundaries are upwards of 0.4 s. The delay signal is found to diminish to <0.1 s for older Proterozoic basins, likely due to compaction and metamorphism of the sediments over time. As an application of the method, a comparison with measurements of sediment thickness from local boreholes allows for a straightforward predictive relationship between the delay time and the cover thickness to be defined. This offers future widespread potential, providing a simple and cheap way to characterise the sediment thickness in under-explored areas from passive seismic data.

SummaryEarthquake forecasting poses significant challenges, especially due to the elusive nature of stress states in fault systems. To tackle this problem, we employ features derived from seismic catalogs obtained from acoustic emission (AE) signals recorded during triaxial stick-slip experiments on natural fractures in three Westerly granite samples. We extracted 47 physically explainable features from AE data that described spatio-temporal evolution of stress and damage in the vicinity of the fault surface. These features are then subjected to unsupervised clustering using the K-means method, revealing three distinct stages with a proper agreement with the temporal evolution of stress. The recovered stages correspond to the mechanical behavior of the rock, characterized as initial stable (elastic) deformation, followed by a transitional stage leading to an unstable deformation prior to failure. Notably, AE rate, clustering-localization features, fractal dimension, b-value, interevent time distribution, and correlation integral are identified as significant features for the unsupervised clustering. The systematically evolving stages can provide valuable insights for characterizing preparatory processes preceding earthquake events associated with geothermal activities and waste-water injections. In order to address the upscaling issue, we propose to use the most important features and, in case of normalization challenge, removing non-universal features, such as AE rate. Our findings hold promise for advancing earthquake prediction methodologies based on laboratory experiments and catalog-driven features.

SummarySubduction earthquakes show complex spatial and temporal rupture patterns, exhibiting events of varied sizes, which rupture distinct or overlapping fault segments. Elucidating first-order controlling conditions of rupture segmentation and return periods of large earthquakes is therefore critical for seismic and tsunami hazard estimations. The Chilean subduction zone frequently hosts several Mw > 8 earthquakes, with heterogeneous recurrence rates and locations. Here, we implement 3D quasi-dynamic rate and state frictional models to investigate the role of plate interface geometry on the distribution of interseismic coupling and coseismic ruptures in Central Chile. First, we develop synthetic-parametric models that show how dip and strike variations may increase the probabilities to produce partial seismic barriers, which tend to avoid the production of large earthquake ruptures and modulate rupture lengths. Then, we simulate the subduction seismic cycle processes on Central Chile (25ºS-38°S), imposing depth-dependent frictional properties on a realistic non-planar 3D subduction interface geometry. Similar to results obtained for synthetic-parametric models, after 5000 years of simulation, regions with abrupt dip or strike changes increase the probabilities of stopping coseismic propagation of simulated Mw8.0-9.0 earthquakes. Our simulated earthquake sequences on the Central Chile subduction zone delimit rupture areas that match geometrical interface features and historical earthquakes, results that point to the crucial role of fault interface geometry on seismic cycle segmentation along strike.

SummaryTrain traffic is a powerful source of seismic waves, with many applications for passive seismic imaging. Seismic signals were recorded a few meters away from the railway track. These records display harmonious waveforms below 15 Hz for trains driving at speeds of around 100 km/h. The sensors record an apparent wavelet emitted by the interaction of the axle on a few of the closest sleepers. From this, we build a simple modeling tool using shifted wavelets to simulate a train signal. The relationship involves the varying train speed, the distances between each axle, and the wavelet emitted by each axle. We propose a non-linear deconvolution method to invert this relationship. We use a local minimisation algorithm with gradients derived by the adjoint state method, and use a frequency continuation technique. A linearised picking-based inversion initializes the non-linear inversion. On real data, we apply this automatic workflow to 300 train passages, with an excellent match between the best simulation and the data. We identify the trains decelerating as they enter a train station. We also identify the train type based on inverted wheel spacing with centimetric accuracy. The inverted wavelets are consistent with the assumption that trains emit seismic waves by bending the rail above sleepers, although the theory does not explain why the inverted wavelet is not zero-phase. This automated kinematic inversion algorithm may allow for contact-less railway monitoring, and be used for source characterisation for subsurface monitoring below railway tracks.

SummaryThis study presents a new robust statistical framework, in which to measure relative differences, or deviations from a hypothetical reference value, of Gutenberg-Richter b-value. Moreover, it applies this method to recent seismicity in Italy, to find possible changes of earthquake magnitude distribution in time and space. The method uses bootstrap techniques, which have no prior assumptions about the distribution of data, keeping their basic features. Excluding Central Italy, no significative b-value variation is found, revealing that the frequency-magnitude distribution exponent is substantially stable or that data are not able to reveal hidden variations. Considering the small size of examined magnitude samples, we cannot definitively decide if the higher b-values in Central Italy, consistently founded by all applied tests, have a physical origin or result from a statistical bias. In any case, they indicate short-lived excursions which have a temporary nature and, therefore, cannot be associated solely to spatial variations in tectonic framework. Both the methodological issues and the results of the application to seismicity in Italy show that a correct assessing of b-value changes requests appropriate statistics, that accurately quantify the low accuracy and precision of b-value estimation for small magnitude samples.

SummaryInterferometric synthetic aperture radar (InSAR) decorrelation that creates great challenges to phase unwrapping has been a critical issue for mapping large earthquake deformation. Some studies have proposed a “remove-and-return model” solution to tackle this problem, but it has not been fully validated yet, and therefore has rarely been applied to real earthquake scenarios. In this study, we use the 2023 Mw 7.8 and 7.6 earthquake doublet in Turkey and Syria as a case example to develop an iterative modeling method for InSAR-based coseismic mapping. We first derive surface deformation fields using Sentinel-1 offset tracking and Sentinel-2 optical image correlation, and invert them for an initial coseismic slip model, based on which we simulate InSAR coseismic phase measurements. We then remove the simulated phase from the actual Sentinel-1 phase and conduct unwrapping. The simulated phase is added back to the unwrapped phase to produce the final phase measurements. Comparing to the commonly-used unwrapping method, our proposed approach can significantly improve coherence and reduce phase gradients, enabling accurate InSAR measurements. Combining InSAR, offset tracking and optical image correlation, we implement a joint inversion to obtain an optimal coseismic slip model. Our model shows that slip on the Çardak Fault is concentrated on a ∼100 km segment; to both ends, slip suddenly diminished. On the contrary, rupture on the East Anatolian Fault Zone propagated much longer as its geometry is fairly smooth. The iterative coseismic modeling method is proven efficient and can be easily applied to other continental earthquakes.

SummarySite amplification factors (SAFs) of seismic ground motions are essential in evaluating and estimating seismic hazards. In our previous study, the authors proposed a simple and cost-effective method to estimate a SAF based on a deep neural network (DNN) model and microtremor horizontal-to-vertical spectral ratio (MHVR). Since the previous DNN model was based on the observed SAFs and MHVRs within a limited district in Japan, the applicability of the previous model to non-source regions with different site conditions was limited. This study explored the application of a transfer learning (TL) technique to adapt an existing (pre-trained) DNN model to new regions and a different database. The SAFs obtained through generalized spectral inversion technique (GIT) at the seismic observation stations (K-NET and KiK-net) in Japan were collated as the ground truth for site effects. MHVRs recorded at the stations in several districts of Japan were collected to construct a dataset for the development of the TL model. Subsequently, a TL model was constructed, leveraging the neural network layers and their weights from the pre-trained model while incorporating additional neural network layers to enhance the performance. During the training process, a total data set of 112 sites was divided into training set, validation set, and external test set by 1:1:5. Utilizing a cross-validation strategy, the residuals between pSAFs (pseudo-SAFs) estimated by the TL model and the observed SAFs were analyzed for the external test set containing 80 sites. The results showed that the TL model outperformed the pre-trained DNN model. The cross-validation results demonstrated that almost consistent prediction results were obtained regardless of any combination of 16 sites selected as the training set. Furthermore, by contrasting the influence of varying training set sizes on the performance of the TL model and comparing the TL model to a DNN model with an extended training set, the effectiveness of constructing the model with the limited number of data (16 sites) was ascertained. Finally, the effectiveness and limitations of the TL model were evaluated using MHVRs with peak frequencies falling outside the training set's range.

SummaryAs the high-frequency analog to field-scale earthquakes, acoustic emissions (AEs) provide a valuable complement to study rock deformation mechanisms. During the load-stepping creep experiments with CO2-saturated water injection into a basaltic sample from Carbfix site in Iceland, 8791 AE events are detected by at least one of the seven piezoelectric sensors. Here, we apply a cross-correlation-based source imaging method, called geometric-mean reverse-time migration (GmRTM) to locate those AE events. Besides the attractive picking-free feature shared with other waveform-based methods (e.g., time-reversal imaging), GmRTM is advantageous in generating high-resolution source images with reduced imaging artifacts, especially for experiments with relatively sparse receivers. In general, the imaged AE locations are found to be scattered across the sample, suggesting a complicated fracture network rather than a well-defined major shear fracture plane, in agreement with X-ray computed tomography imaging results after retrieval of samples from the deformation apparatus. Clustering the events in space and time using the nearest-neighbor approach revealed a group of “repeaters”, which are spatially co-located over an elongated period of time and likely indicate crack, or shear band growth. Furthermore, we select 2196 AE events with high signal-noise-ratio (SNR) and conduct moment tensor estimation using the adjoint (back-propagated) strain tensor fields at the locations of AE sources. The resulting AE locations and focal mechanisms support our previously assertion that creep of basalt at the experimental conditions is accommodated dominantly by distributed micro-cracking.

SummaryMuography is an imaging technique that relies on the attenuation of the muon flux traversing geological or anthropogenic structures. Several simulation frameworks help to perform muography studies by combining specialised codes: for muon generation through muon transport to muon detector performance. This methodology is precise but requires significant computational resources and time. We present an end-to-end python-based MUographY Simulation Code, which implements a muography simulation framework capable of rapidly estimating muograms of any geological structure worldwide. This framework considers the generated muon flux as the observation point; the energy loss of muons passing through the geological target; the integrated muon flux detected by the telescope and estimates the 3-dimensional density distribution of the target using Algebraic Reconstruction Techniques. The simulations ignore the relatively small muon flux variance caused by geomagnetic effects, solar modulation, and atmospheric conditions. We validate the code performance by comparing our simulation results with data from other frameworks.

SummaryThe diffusion field theory has been widely used to interpret ambient noise wave fields. Based on this theory, 1D subsurface velocity structure inversion method is developed. However, few studies have referred to the Noise Horizontal to Vertical spectral ratio (NHV) inversion of 3D subsurface velocity structures, and almost no effective 3D NHV inversion tools have been developed. To develop a useful tool for obtaining 3D soil layer velocity structures, we combined the NHV forward calculation formula derived from diffusion field theory with the guided Monte Carlo algorithm and then extended the single-point NHV inversion to multipoint joint inversion through a joint objective function. Subsequently, a new three-dimensional soil layer velocity structure inversion method was proposed. Subsequently, a synthetic 2D case was used to verify the proposed method. Finally, the proposed method was applied to the Xiangtang Array in Tangshan, China, to identify the 3D velocity structures of the site based on noise observations. The results show that the proposed multipoint joint 3D inversion method is effective for identifying 3D underground velocity structures.

SummaryThe physical mechanisms controlling deep-focus earthquakes, or those observed at depths greater than 300 km, remain enigmatic. The leading processes by which deep-focus earthquakes are thought to occur include transformational faulting, thermal runaway, and dehydration embrittlement, but distinguishing observations in support of one or more prevailing mechanisms are needed. In this study, we use a modified back-projection method, data recorded by the Hi-net array in Japan, and a three-dimensional velocity model to produce source images of 19 deep-focus earthquakes within the Izu-Bonin subduction zone. We find that the rupture properties and fault plane orientations of imaged events separate according to reported moment magnitude, indicating the distinct operation of two moment-dependent causal mechanisms of deep-focus earthquakes in this region. We discuss these results in the context of previous observational, laboratory, and numerical studies and emphasize the importance of continued research to validate the dual-mechanism hypothesis both in and outside Izu-Bonin. Such work may not only improve our understanding of the nucleation and propagation of deep-focus earthquakes, but also help clarify slab structure and subduction zone dynamics.