Traveltime Tomography Research Articles

Illumination compensation for the adjoint-state first-arrival traveltime tomography in VTI media

The sensitivity kernel of each parameter is angle-dependent in an anisotropic traveltime multiparameter inversion. To analyze the angle illumination in a transversely isotropic media with a vertical symmetry axis (VTI) by the adjoint state traveltime tomography (AST) method, a circular observation system is designed with receivers radially distributed at different propagated angles. Then, the nonuniform angular illumination in the gradient of different parameters in the VTI AST is revealed. The nonuniform angular illumination seriously hinders the simultaneous updates of multiple parameters toward the true solution and decreases the accuracy of the anisotropic multiparameter inversion. To overcome the illumination problem, a preconditioned AST method (PAST) is proposed by utilizing the diagonal elements of the approximate Hessian matrix as the angle illumination compensation for multiparameter gradients in the VTI medium. The PAST method not only compensates for angular illumination but also balances the ray density and eliminates gradient singular values, which improves the inversion accuracy for multiple parameters. The BP synthetic model test shows that the PAST method based on angular illumination compensation can invert the VTI multiparameter models more accurately. Therefore, the PAST inversion model can predict first-arrival traveltime information more accurately. Field data application in the East China Sea indicates that the inverted velocity of the proposed PAST fits well with the logging velocity, and the inverted anisotropic parameter agrees well with the Backus average values. Consequently, the image quality based on the PAST inverted models is slightly improved. It is proven that the proposed PAST method can better realize the near-surface multiparameter inversion of the VTI medium.

Crustal velocity structure in the South Yellow Sea revealed by the joint tomographic inversion of reflected and refracted seismic waves

The crustal velocity structure in the South Yellow Sea (SYS) Basin is crucial for understanding the basin’s geological structure and evolution. OBS (ocean-bottom station) data from the OBS2013 line have been used to determine the crustal velocity structure in the SYS. The velocity model of the upper crust in the northern SYS was determined using first-arrival traveltime tomography. The model showed a higher resolution shallow crustal velocity structure but a lower resolution middle-lower crustal velocity structure. The crustal velocity structure, together with the Moho discontinuity in the SYS Basin, was also constructed using a human–computer interactive traveltime simulation, and the result was highly dependent on the prior knowledge of the operator. In this study, we reconstructed a crustal velocity model in the SYS Basin using a joint tomographic inversion of the traveltime and its gradient data of the reflected and refracted waves picked from the OBS data. The resolution of the inverted velocity structure from shallow-to-deep crust was improved. The results revealed that the massive high-velocity body below the Haiyang Sag of the Jiaolai Basin extends to the Qianliyan Uplift in the SYS; the low-velocity Cretaceous strata directly cover the pre-Sinitic metamorphic rock basement of the Sulu orogenic belt; and the thick Meso-Paleozoic marine strata are retained beneath the Meso–Cenozoic continental strata in the northern depression. The Moho depth in the SYS Basin ranges from 28 to 32 km.

Fast 3D “wave-consistent” ray tomography based on a polynomial representation of a model

In the paper, we describe an original 3D travel-time tomography approach. It is based on the new realization of the bending method which to some extent takes into account the band-limited nature of real seismic signals propagation. As a result, two-point ray tracing provides more reliable ray trajectories and travel times in complex media. Another original feature of the proposed tomography is that the model is represented using the Chebyshev polynomials. Such parameterization allows analytical calculation of travel times and their derivatives with respect to model parameters and significantly reduces the number of parameters to be recovered during inversion compared to more common grid tomography. In certain situations, the proposed approach provides significant computational advantages. Numerical examples prove its efficiency.

Supervirtual Refraction Interferometry in the Radon Domain

Accurate picking of seismic first arrivals is very important for first arrival travel time tomography, but the first arrivals appearing at far offsets are often more difficult to pick accurately due to the low signal-to-noise ratio (SNR). The conventional supervirtual refraction interferometry (SVI) method can improve the SNR of first arrivals to a certain extent; however, it is not suitable for seismic data that interfered by strong noise. In order to better process the first arrivals at far offsets with serious noise interference, we propose a modified method, in which SVI implemented in the Radon domain (RDSVI) due to the cross-correlation in the Radon domain have a better effect. According to the kinematic characteristics of first arrival refractions, SVI is performed in the linear Radon domain. Both synthetic data and field data demonstrate the proposed method can enhance the effective signal and attenuate the strong noise simultaneously, so as to significantly improve the SNR of the first arrival data. Meanwhile, the RDSVI method is tested on the first arrival data with missing traces, which proves that this method can overcome the influence of abnormal traces and is suitable for the reconstruction of sparsely sampled seismic data.

First-Break Picking of Large-Offset Seismic Data Based on CNNs with Weighted Data

Deep reflection seismic data are usually accompanied by large-offset data, and the accurate and rapid identification of the first arrivals of seismic records plays an important role in eliminating the effects of topography and other factors that increase with the increasing offsets. In this paper, we propose a method based on convolutional neural networks (CNNs) that can accurately identify the first arrivals of large-offset seismic data. A time window for linear dynamic correction was established to convert the raw seismic data into rectangular images so as to reduce the amount of invalid sample data and improve the training efficiency. In order to enhance the prediction effect of the far-offset first arrivals, we propose the strategy of adjusting the weight of the far-offset data to increase the weight of the far-offset data in the training dataset and, thus, to improve the first arrival accuracy. The manually picked first arrivals are used as labels and the input to the CNNs for training, and the full-offset first arrivals are the output. The travel time tomography velocity is modeled and compared based on the first arrivals obtained through manual picking, industrial software automatic picking, and CNN prediction. The results show that the application of CNNs to large-offset seismic datasets can help researchers to obtain the first arrivals at different offsets, while the inclusion of far-offset weights can effectively improve the modeling depth of the tomography inversion, and the accuracy of the results is high.

Numerical solution of the 3-D travel time tomography problem

The first numerical solution of the 3-D travel time tomography problem is presented. The globally convergent convexification numerical method is applied.

Phase-preserving theory and its linearization approximation for forward scattering field of scalar acoustic wave equation

The conventional wave-equation linearization methods, such as the first-order Born or Rytov approximation, always implicitly imply a weak-scattering assumption, making it valid only for weak perturbation models. To extend the wave-equation linearization theory to strong perturbation models, we consider a scenario that the reference model is smooth within the scale of the incident wave length, and propose a phase-preserving method which can predict the phase perturbation of forward scattering wave field. First, we introduce the WKBJ approximation to the scattered- and incident wave fields so that the integral of the unknown solution (i.e. the scattered field) in the nonlinear Ricatti integral equation can be replaced by the integral of scattering-angle and model perturbation, yielding an explicit expression of the scattered field. Theoretical derivation shows that the proposed phase-preserving method can accurately predict the phase-perturbation of forward scattered wave field regardless of the strength of velocity perturbations for one-dimensional wave propagation problem. To apply the phase-preserving approximation to the inverse problem, we further consider a scenario of small-angle forward propagation. In this case, the phase-preserving approximation can be linearized by neglecting the influence of scattering angles, leading to a linear relation between the scattered field and the model perturbation, which we refer to as the generalized Rytov approximation. Numerical experiments demonstrate that the generalized Rytov approximation can predict the phase perturbation of the scattered field with higher accuracy for small-angle forward propagation, and is suitable for strong model perturbations. The generalized Rytov approximation extends the validity and the scope of application of the traditional Rytov approximation. In specific application fields such as the seismic traveltime tomography or medical ultrasonic transmission imaging, a new traveltime/phase sensitivity kernel can be derived by replacing the conventional Rytov approximation with the proposed method, which can increase the inversion accuracy and speed up the convergence.

Bayesian seismic travel-time cross-hole tomography in vertically transversely isotropic media

We apply a transdimensional, hierarchical Markov chain Monte Carlo sampling algorithm (McMC) for 2-D cross-hole travel-time tomography in transversely isotropic media with vertical symmetry axis. The McMC approach has several advantages compared to classical inversion approaches: It is a global search, the high number of tested models allows the statistical analysis including the calculation of a reference model as well as uncertainty estimation, no initial models or regularization parameters are needed, the amount of data noise is automatically determined, and the model parametrization is data dependent and self-adjusting. For the forward solution a FD Fast Marching method utilizing second-order Godunov schemes is used. The performance of the approach is first tested on synthetic datasets to evaluate the potential and possible limitation to recover anisotropic models. We have shown that the recovery of models described by 2 anisotropic parameters (Thomsen parameters) and the vertical velocity is possible for observation scenarios with good distribution of sources and receivers. For more realistic observational geometries (i.e. cross-hole experiments), the recovery of the 3 parameters is limited, but still possible for example for the elliptical anisotropic case (ε = δ) or regarding the horizontal velocity. Finally we applied the McMC approach to a well-studied real cross-hole data set from the MALLIK 2002 research program and compared the results with previous conventional inversions.

Segmentation of continental Indian plate by the Narmada‐Son diffuse plate boundary

The Narmada‐Son lineament is a major tectonic feature in Central India. Here we employ data from geodetic measurements of crustal deformation and upper mantle P‐wave travel time tomography and propose that the Narmada‐Son palaeo‐rift is a present‐day diffuse plate boundary that segments the continental part of the Indian “composite” plate into two component plates, the North and South‐Indian plate. Tomographic images reveal that the collision of an indenter, associated with the relatively faster moving South‐Indian plate, with the Narmada‐Son diffuse plate boundary acted as a pivot point exerting a torque on one of the component plates along with mechanical rotation. This new proposal along with a new identity of the Narmada‐Son diffuse plate boundary is in better agreement with the estimated location of the relative Euler pole across two respective component plates. Our model is also able to better explain the Indian “composite” plate hypothesis, diversity in the topographic architecture along the rift, and contrasting deformation styles on either side of the diffuse plate boundary.

An efficient hybrid mesh method for seismic traveltime computation with irregular surfaces

Accurate and efficient traveltime computation with irregular surfaces is a prerequisite for traveltime tomography, seismic location, and seismic exploration in mountainous areas. The eikonal equation method is computationally favourable and has gained increasing popularity for obtaining traveltime in recent years. Nowadays most algorithms that can solve eikonal equation with irregular surfaces are based on fully unstructured triangular mesh or surface-fitting mesh. However, these methods have poor performance in terms of computational cost and accuracy. To cope with this problem, we develop a hybrid mesh method based on fast-sweeping-based eikonal solver with rectangular and unstructured triangular mesh for obtaining traveltime with irregular surfaces. In other words, we propose to divide the model into several different regions, which are discretized by rectangular mesh and unstructured triangular mesh respectively. Near the irregular surfaces of model, we use fast sweeping method based on unstructured triangular mesh which releases the full potential of unstructured triangular mesh’s flexibility in accurately describing the irregular geometry of computation domain. For the rest part of the model, we adopt rectangular mesh to make full use of high efficiency of the corresponding high-order fast sweeping algorithm. At the interface of two types of mesh, we need interpolate the known nodes’s traveltimes to get the boundary condition of another mesh, so as to start the fast sweeping method therein. It is found that the accuracy of hybrid mesh method can be increased by 83%, the computation time is reduced by 72%, and the storage requirement is reduced by 72%, compared with the fast sweeping method with full triangular mesh, when the rectangular mesh proportion is 80% in our model. Therefore, this proposed hybrid mesh method can greatly enhance the computational efficacy and serve as an potential forward engine in seismic inversions.

Time-Lapse Cross-Well Monitoring of CO2 Sequestration Using Coda Wave Interferometry

In this study, we explored the capability of coda wave interferometry (CWI) for monitoring CO2 storage by estimating the seismic velocity changes caused by CO2 injection. Given that the CWI method is highly efficient, the primary aim of this study was to provide a quick detection tool for the long-term monitoring of CO2 storage safety. In particular, we looked at monitoring with a cross-well geometry. We also expected that CWI could help to reduce the inversion errors of existing methods. Time-lapse upgoing waves and downgoing waves from two-component datasets were utilized to efficiently monitor the area between the wells and provide a quick indication of possible CO2 leakage. The resulting mean velocity changes versus the depth indicated the depth where velocity changes occurred. Combining the upgoing and downgoing wavefields provided a more specific indication of the depth range for changes. The calculated velocity changes were determined using the time shift between the time-lapse wavefields caused by CO2 injection/leakage. Hence, the resulting velocity changes were closely related to the ratio of propagation path length through the CO2 injection/leakage layer over the length of the entire travel path. The results indicated that the noise level and repeatability of the time-lapse datasets significantly influenced the results generated using CWI. Therefore, denoising and time-lapse processing were very important for improving the detectability of any change. Applying CWI to time-lapse cross-well surveys can be an effective tool for monitoring CO2 in the subsurface at a relatively low computational cost. As a highly efficient monitoring method, it is sensitive to changes in the seismic response caused by velocity changes in the subsurface and provides additional constraints on the inversion results from conventional travel time tomography and full waveform inversion.

Preconditioned transmission + reflection joint traveltime tomography with adjoint‐state method for subsurface velocity model building

AbstractAdjoint‐state‐based transmission and reflection traveltime inversion are promising approaches for reconstructing a macro velocity field in both shallow and deep depths. However, the contaminated gradient by the uneven illumination impedes the recovery of an optimal velocity model and results in slow convergence rate in inversion. To mitigate this issue, we propose a novel preconditioned transmission + reflection joint traveltime tomography method with transmission and reflection illumination compensation. The proposed illumination compensation operator is theoretically demonstrated as the diagonal elements of the approximated Hessian, which can be efficiently computed in a way similar to gradient calculation. The proposed method can construct a more accurate velocity model compared to the uncompensated counterpart. Convergence is also accelerated by the relatively uniform velocity updates from the shallow to the deep parts. Furthermore, velocity updating and reflector imaging are simultaneously iterated during the inversion. Synthetic numerical example based on the inclusion model indicates that the proposed preconditioned method can result in uniform gradient values of both transmitted and shot/receiver‐side reflected traveltime. The complex foothill model test shows the improved inversion accuracy of the proposed method with an increased convergence rate and improved reflection traveltime predictions using the inverted velocity model. Experiment on the small foothill model indicates that the reflector and velocity can be updated simultaneously, and the inversion accuracy can be significantly improved by selecting multi reflectors. Finally, the effectiveness and practicality of the proposed method are demonstrated by a field data application in East China Sea.

Crustal Deformation in the Sierra Nevada and Walker Lane Region Inferred From P‐Wave Azimuthal Anisotropy

AbstractCrustal deformation in the Sierra Nevada and Walker Lane is mainly investigated by geologic and geodetic constraints which suffer from a lack of depth information. In this study, we construct a new depth‐dependent azimuthally anisotropic P‐wave velocity (Vp) model of this area to investigate regional deformation regimes in the upper and middle crust. The model is built based on adjoint‐state traveltime tomography of ∼650,000 local direct P arrivals collected from 1967 to 2021. The average Vp structure agrees well with previous tomographic models with refined velocity features benefiting from the ray‐free adjoint‐state tomography and more arrival time data. The azimuthal anisotropy of Vp reveals distinct differences between the rigid Sierra Nevada block and the Walker Lane shear zone: (a) The western Sierra Nevada shows weak anisotropy (<2%) and NW‐SE oriented fast velocity directions that are parallel to the Pacific‐North American plate boundary, mainly reflecting preserved paleofabrics developed in past subduction processes. (b) In the shallow Walker Lane (<4 km), the orientation of fast‐velocity directions is NNE controlled by regional compressive stress, and changes to NNW at greater depths in the north under shear deformation regimes. (c) In the central Walker Lane, strong anisotropy with ENE‐oriented P fast axes is imaged at 6–16 km from the area south of Lake Tahoe (2%–4%) to the Long Valley volcanic area (4%–8%) because of clockwise crustal block rotations and crustal fluids. The proposed anisotropic Vp model provides new insight into how past and present‐day deformation is accommodated in this tectonic complex plate boundary region.

Neural Eikonal solver: Improving accuracy of physics-informed neural networks for solving eikonal equation in case of caustics

The concept of physics-informed neural networks has become a useful tool for solving differential equations due to its flexibility. A few approaches use this concept to solve the eikonal equation that describes the first-arrival traveltimes of waves propagating in smooth heterogeneous velocity models. However, the challenge of the eikonal is exacerbated by the velocity models producing caustics, resulting in instabilities and deterioration of accuracy due to the non-smooth solution behavior. In this paper, we revisit the problem of solving the eikonal equation using neural networks to tackle caustic pathologies. We introduce the novel Neural Eikonal Solver (NES) for solving the isotropic eikonal equation in two formulations: the one-point problem is for a fixed source location; the two-point problem is for an arbitrary source-receiver pair. We present several techniques which provide stability in the case of caustics: improved factorization; non-symmetric loss function based on Hamiltonian; gaussian activation; symmetrization. In our tests, NES showed the relative mean-absolute error of 0.2-0.4% from the second-order factored Fast Marching Method with a similar inference time, and outperformed existing neural-network solvers giving 10-60 times lower errors and 2-30 times faster training. The one-point NES provides the most accurate solution, while the two-point NES gives slightly lower accuracy but an extremely compact representation with all spatial derivatives. It can be useful in many seismic problems: massive computations of traveltimes for millions of source-receiver pairs in Kirchhoff migration; modeling of ray amplitudes using spatial derivatives; traveltime tomography; earthquake localization; ray multipathing analysis.

P-, SV- and SH-wave eikonal approximations for HTI media based on perturbation theory

SUMMARY The eikonal equation for transversely isotropic media is a complex quartic partial differential equation that causes instability in traveltime inversion. The phenomenon of shear wave splitting in a transversely isotropic model with a horizontal symmetry axis (HTI) provides information for inversion. However, the involvement of shear waves makes the computation of forward modelling and inversion more time-consuming than when only using Pwaves. Most studies of traveltime approximation based on the eikonal equation are based on an acoustic assumption, which does not work for shear waves. Therefore, based on perturbation theory, we derive approximate solutions of the eikonal equation in a heterogeneous elastic HTI model for P, SV and SH waves to improve the traveltime computation efficiency. The approximate solutions of linear, quadratic and Shanks forms are based on the linearization of the traveltime function in terms of fracture weaknesses in an isotropic background, which are suitable for joint P-, SV- and SH-wave traveltime inversion and tomography for fracture parameters. When applied for a homogenous model, the proposed formulation results in a simple and accurate traveltime approximation. The first- and second-order perturbation coefficients are used to analyse the sensitivity of the traveltime to fracture weakness. To improve the accuracy of our approximation, we also expand the approximation of the P-wave traveltime with respect to the anellipticity parameter η in an elliptic background, which is evaluated by performing an error test. Based on the accuracy and stability of the perturbation, we select a set of linear equations with respect to the fracture weaknesses and a set of more accurate nonlinear equations to serve different needs.

Cross-Well Reflection Imaging at the Verkhnekamskoe Potash Deposit

Abstract —We propose the technique of extended processing of cross-well seismic data for reflection imaging. Based on the analysis of the wave field and synthetic modeling, a graph for processing cross-well data is developed to separate the reflections in the presence of a boundary of a sharp change in the velocity of elastic waves. The migration of the reflected waves from the time scale to the depth is performed by solving a forward problem for each source–receiver pair. As a result, the position of the reflection points is calculated, after which all the traces are stacked based on the binning grid. The input information for performing migration and solving a forward problem is the velocity characteristic of the massif, which is calculated using traveltime tomography on direct body waves obtained from the same data set of cross-well survey and vertical seismic profiling. The resulting depth seismic section has a much higher resolution than that of vertical seismic profiling and ground-based shallow seismic studies. This opens up the opportunity of identifying different small natural or technogenic objects in the interwell space.

Rapid magma ascent beneath La Palma revealed by seismic tomography

For the first time, we obtained high-resolution images of Earth's interior of the La Palma volcanic eruption that occurred in 2021 derived during the eruptive process. We present evidence of a rapid magmatic rise from the base of the oceanic crust under the island to produce an eruption that was active for 85 days. This eruption is interpreted as a very accelerated and energetic process. We used data from 11,349 earthquakes to perform travel-time seismic tomography. We present high-precision earthquake relocations and 3D distributions of P and S-wave velocities highlighting the geometry of magma sources. We identified three distinct structures: (1) a shallow localised region (< 3 km) of hydrothermal alteration; (2) spatially extensive, consolidated, oceanic crust extending to 10 km depth and; (3) a large sub-crustal magma-filled rock volume intrusion extending from 7 to 25 km depth. Our results suggest that this large magma reservoir feeds the La Palma eruption continuously. Prior to eruption onset, magma ascended from 10 km depth to the surface in less than 7 days. In the upper 3 km, melt migration is along the western contact between consolidated oceanic crust and altered hydrothermal material.

Advances in seismic imaging of magma and crystal mush

Seismic imaging methods have provided detailed three-dimensional constraints on the physical properties of magmatic systems leading to invaluable insight into the storage, differentiation and dynamics of magma. These constraints have been crucial to the development of our modern understanding of magmatic systems. However, there are still outstanding knowledge gaps resulting from the challenges inherent in seismic imaging of volcanoes. These challenges stem from the complex physics of wave propagation across highly heterogeneous low-velocity anomalies associated with magma reservoirs. Ray-based seismic imaging methods such as travel-time and surface-wave tomography lead to under-recovery of such velocity anomalies and to under-estimation of melt fractions. This review aims to help the volcanologist to fully utilize the insights gained from seismic imaging and account for the resolution limits. We summarize the advantages and limitations of the most common imaging methods and propose best practices for their implementation and the quantitative interpretation of low-velocity anomalies. We constructed and analysed a database of 277 seismic imaging studies at 78 arc, hotspot and continental rift volcanoes. Each study is accompanied by information about the seismic source, part of the wavefield used, imaging method, any detected low-velocity zones, and estimated melt fraction. Thirty nine studies attempted to estimate melt fractions at 22 different volcanoes. Only five studies have found evidence of melt storage at melt fractions above the critical porosity that separates crystal mush from mobile magma. The median reported melt fraction is 13% suggesting that magma storage is dominated by low-melt fraction crystal mush. However, due to the limits of seismic resolution, the seismological evidence does not rule out the presence of small (&amp;lt;10 km3) and medium-sized (&amp;lt;100 km3) high-melt fraction magma chambers at many of the studied volcanoes. The combination of multiple tomographic imaging methods and the wider adoption of methods that use more of the seismic wavefield than the first arriving travel-times, promise to overcome some of the limitations of seismic tomography and provide more reliable constraints on melt fractions. Wider adoption of these new methods and advances in data collection are needed to enable a revolution in imaging magma reservoirs.

An Envelope Travel-Time Objective Function for Reducing Source–Velocity Trade-Offs in Wave-Equation Tomography

In conventional cross-correlation (CC)-based wave-equation travel-time tomography, wrong source wavelets can result in inaccurate velocity inversion results, which is known as the source–velocity trade-off. In this study, an envelope travel-time objective function is developed for wave-equation tomography to alleviate the non-uniqueness and uncertainty due to wrong source wavelets. The envelope of a seismic signal helps reduce the waveform fluctuations/distortions caused by variations of the source time function. We show that for two seismic signals generated with different source wavelets, the travel-time shift calculated by cross-correlation of their envelopes is more accurate compared to that obtained by directly cross-correlating their waveforms. Then, the CC-based envelope travel-time (ET) objective function is introduced for wave-equation tomography. A new adjoint source has also been derived to calculate the sensitivity kernels. In the numerical inversion experiments, a synthetic example with cross-well survey is first given to show that compared to the traditional CC travel-time objective function, the ET objective function is relatively insensitive to source wavelet variations and can reconstruct the elastic velocity structures more reliably. Finally, the effectiveness and advantages of our method are verified by inversion of early arrivals in a practical seismic survey for recovering near-surface velocity structures.

Seismic Structure Beneath the Molucca Sea Collision Zone from Travel Time Tomography Based on Local and Regional BMKG Networks

The Molucca Sea Plate, and Sangihe and Halmahera plates have a complex tectonic setting and interact to create the Molucca Sea Collision Zone. We re-picked 1647 events recorded from 2010 to 2017 from 32 of The Agency for Meteorology, Climatology, and Geophysics (BMKG) stations and obtained P- and S-arrivals of ~17,628 phases. Hypocenter locations were determined using the software NonLinLoc. Then, we performed a travel time tomography in order to image the subsurface and approximate the Molucca Sea Plate subduction angle beneath Sulawesi’s north arm, the relationship subduction zone and volcanic activity in Halmahera, and depth comparison of the Molucca Sea Plate. Our results show that (i) high Vp, high Vs, and low Vp/Vs are associated with the Molucca Sea Plate beneath Sulawesi’s north arm, and the approximate subduction angle is ~40°. (ii) Low Vp, low Vs, and high Vp/Vs beneath the northern and southern Halmahera Volcanic Arc are associated with a possible magma source. Volcanoes in the north have experienced eruptions and are dormant in the south. This group of volcanoes is connected by partial melting above the Molucca Sea Plate subducts to the east. (iii) High Vp, high Vs, and low Vp/Vs are interpreted as double subduction of the Molucca Sea Plate. It is submerged deeper in the north compared with the south, which is nearer to the surface.