**5. Conclusions**

In this research, we performed travel time tomographic inversion using an initial 3D velocity model. These are obtained by manual re-picking of the same number of P- and S-arrivals and simultaneous determination of the 3D structure for the Vp, Vs, and Vp/Vs ratio, as well as relocating the hypocenter location. The results show a high Vp, high Vs, and low-Vp/Vs anomaly beneath the Minahasa Peninsula, Sulawesi's north arm, with the Celebes Sea (cross-section A-A , south–north) dipping to the north at a ~40◦ angle. Beneath the Celebes Sea, Sangihe Island, Molucca Sea, the Morotai Basin (cross-section C-C , west–east), and Sulawesi's north arm, the Molucca Sea, and Halmahera Island (crosssection D-D , west–east) can be interpreted as a double subduction of the Molucca Sea Plate. The tomographic image beneath the southern part of the Halmahera volcanic arc (Bacan Island) shows low Vp, low Vs, and high Vp/Vs, which indicates a possible magma source; although the volcano has been dormant until the present. The tomographic image beneath the group of volcanoes in the northern arm of Halmahera also shows low Vp, low Vs, and high Vp/Vs, which can be interpreted as a possible magma source. Differing from

volcanoes in the south, however, some volcanoes in this area have erupted and are emitting volcanic ash. Both the volcanic groups on Bacan Island and the northern arm of Halmahera are connected by partial melting, as indicated by low Vp, low Vs, high Vp/Vs above the Molucca Sea Plate, which subducts to the east.

The upper curvature of the Molucca Sea Plate in the northern part is submerged deeper into the upper mantle layer at a depth of ~150 km compared with the southern part. It was also found that beneath the Central Ridge of the Molucca Sea, at latitude 1.10 ◦N, the low Vp, low Vs, and high Vp/Vs in the crust and mantle can be interpreted as highly fractured and water-saturated rock material; while high Vp, high Vs, and low Vp/Vs are observed in the western (beneath the ST) and eastern parts (beneath the HT), which can be interpreted as rock material that experienced micro-cracks that have since closed due to high pressure. Further north, the low Vp, low Vs, and high Vp/Vs anomalies shift eastward and extend to HT with a clear boundary between low and high velocity below the Central Ridge of the Molucca Sea.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/app122010520/s1, Figure S1: Molucca Sea Plate profile overlaid with seismic profile from 1959–1966 which forms a double subduction pattern was modified from [7]. NW: Northwest, SE: Southeast. Sign ?: the length of the plate subduction to the west and east is not yet known. The black dot denotes earthquakes.; Figure S2: Histogram of the number of phases per station; Figure S3: Horizontal and vertical ray path coverage (West-East: below and North-South: right). The brown line is the ray path connecting hypocenters locations with BMKG stations; Figure S4: Wadati diagram for all events used in this study. Tp is arrival times of the P-wave and Ts−Tp is the difference between S- and P-waves arrival times. The Gray dot is the picks of phase time from each station at a certain event, and the solid line is the slope to determine the Vp/Vs value which is 1.76; Figure S5: Trade-off curve between data variance and model variance. We chose damping parameters Vp and Vp/Vs to be 150 and 130, respectively. We use this damping value to achieve a good compromise between data variance and model variance for the inversion, both the real data and the checkerboard test; Figure S6: Vertical cross-section of Derivative Weight Sum (DWS) on profile South-North: A-A', B-B' and profile West-East: C-C', D-D'. The DWS is used to find the number of ray paths that pass through a grid node. The first column is the Vp value, and the second column is Vp/Vs value (logarithmic scale). The black color is the high-resolution zone, and the red triangle is the volcanoes; Figure S7: Vertical cross-section of Ray Hit Count (RHC) on profile South-North: A-A', B-B' and profile West-East: C-C', D-D' for making a thorough cut in a particular DWS. The first column is the Vp value, and the second column is Vp/Vs value (logarithmic scale). The black color is the high-resolution zone, and the red triangle is the volcanoes; Figure S8: Vertical cross-section of Diagonal Resolution Element (DRE) on profile South-North: A-A', B-B' and profile West-East: C-C', D-D' which concerning matrix resolution of R (damped least-squares problems), with an R-value between 0 and 1. The first column is the Vp value, and the second column is Vp/Vs value. The black color is the high-resolution zone, and the red triangle is the volcanoes; Figure S9: The distribution frequency estimated error location. (A–C) are the x, y, and z directions; Figure S10: The distribution between the number of events and the azimuthal gap; Figure S11: The distribution between P- and S-travel times and epicentral distance; Figure S12: Histogram of travel time residual (black and yellow bars color are initial and final travel time residual, respectively). The mean travel time residual decreased from −1.111 s to 0.156 s. The Binwidth histogram is 0.2 s; Figure S13: The curve between the number of iterations and weighted RMS residual. Tomography inversion became convergent after 15 iterations, and the weighted RMS residual decreased by approximately 47%; Figure S14: Histogram of RMS Residual (initial: yellow bars color and final: black bars color). The mean of RMS residual decreased by approximately 33%. The Binwidth histogram is 0.1 s; Figure S15: The horizontal section of CRT profiles from a depth of 0 to 330 km. The purple dashed line is an area that is good well-resolved. The blue color indicates positive perturbations, and the red color indicates negative perturbations relative to the initial model. The first column is Vp perturbation, and the second column is Vp/Vs perturbation; Figure S16: Horizontal tomographic inversion results for Vp, Vs, dan Vp/Vs ratio from a depth of 0 to 330 km. The bright area is an area that good well-resolved, and the blurred area is an area that was not well-resolved; Table S1: The highest value of Vp and Vp/Vs for the DWS, RHC, and DRE resolution tests.

**Author Contributions:** Conceptualization, G.R., B.J.S., A.D.N. and S.R. (Supriyanto Rohadi); methodology, G.R., A.D.N., B.J.S. and S.R. (Supriyanto Rohadi); software, G.R., S.R. (Shindy Rosalia), P.S., M.R., A.A., F.M. and H.A.; validation, G.R., A.D.N., S.R. (Shindy Rosalia) and F.M.; formal analysis, G.R., A.D.N., B.J.S. and S.R. (Supriyanto Rohadi); investigation, A.D.N., B.J.S., S.R. (Supriyanto Rohadi) and S.R. (Shindy Rosalia); resources, G.R., A.D.N., B.J.S., S.R. (Supriyanto Rohadi), S.S. and S.R. (Shindy Rosalia); data curation, G.R., A.D.N. and S.R. (Supriyanto Rohadi); writing—original draft preparation, G.R. and F.M.; writing—review and editing, A.D.N., S.R. (Shindy Rosalia), Z.Z., D.P.S. and S.S.; visualization, G.R., F.M., A.A. and H.A.; supervision, A.D.N., S.R. (Shindy Rosalia), Z.Z., D.P.S. and S.S.; project administration, G.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by BUDI-DN Ministry of Research, Technology and Higher Education (Kemristekdikti) in collaboration with the Indonesian Endowment Fund for Education (LPDP) Number: PRJ-5440/LPDP.3/2016; Addendum Number: PRJ-193/LPDP.4/2019.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The waveform data used in this research were obtained from BMKG. The data supporting the findings of this study are available from the corresponding author upon request.

**Acknowledgments:** Our grateful appreciation goes to the BUDI-DN Ministry of Research, Technology, and Higher Education (Kemristekdikti) in collaboration with the Indonesian Endowment Fund for Education (LPDP), which provides scholarship and research funds. We would also like to thank the Volcanology and Geothermal Laboratory (Geophysical Engineering, Faculty of Mining and Petroleum Engineering, Institut Teknologi Bandung) for providing facilities for conducting research. The Agency for Meteorology, Climatology, and Geophysics (BMKG) provided waveform data. Anthony Lomax created and developed the NLLoc and Seisgram2Kv.7.0 software used in this study. All 2D images are plotted using the GMT tool, and 3D images using Matlab 2018b.

**Conflicts of Interest:** The authors declare no conflict of interest.
