Next Article in Journal
A Hybrid Linear Iterative Clustering and Bayes Classification-Based GrabCut Segmentation Scheme for Dynamic Detection of Cervical Cancer
Next Article in Special Issue
Research on Initial Model Construction of Seismic Inversion Based on Velocity Spectrum and Siamese Network
Previous Article in Journal
Documenting Archaeological Petroglyph Sites with the Use of 3D Terrestrial Laser Scanners—A Case Study of Petroglyphs in Kyrgyzstan
Previous Article in Special Issue
Reverse Time Migration Imaging Using SH Shear Wave Data
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Applied Sciences
Volume 12
Issue 20
10.3390/app122010520
applsci-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Gazali Rachman
1,2,*,
Bagus Jaya Santosa
1,
Andri Dian Nugraha
3,4,
Supriyanto Rohadi
5,
Shindy Rosalia
3,4,
Zulfakriza Zulfakriza
3,4,
Sungkono Sungkono
1,
David P. Sahara
3,4,
Faiz Muttaqy
6,
Pepen Supendi
5,
Mohamad Ramdhan
6,
Ardianto Ardianto
7 and
Haunan Afif
8
1
Physics Department, Institut Teknologi Sepuluh Nopember (ITS), Surabaya 60111, Indonesia
2
Physics Education-FKIP, Pattimura University, Kota Ambon 97233, Indonesia
3
Global Geophysics Research Group, Faculty of Mining and Petroleum Engineering, Institut Teknologi Bandung, Bandung 40132, Indonesia
4
Center for Earthquake Science and Technology (CEST), Research Center for Disaster Mitigation, Institut Teknologi Bandung, Bandung 40132, Indonesia
5
Agency for Meteorology, Climatology and Geophysics (BMKG), Jakarta 83124, Indonesia
6
Research Center for Geological Disaster, National Research and Innovation Agency of Indonesia (BRIN), Bandung 40135, Indonesia
7
Geophysical Engineering Department, Faculty of Mining and Petroleum Engineering, Institut Teknologi Bandung, Bandung 40132, Indonesia
8
Center for Volcanology and Geological Hazards Mitigation, Geological Agency, Ministry of Energy and Mineral Resources, Bandung 40122, Indonesia
*
Author to whom correspondence should be addressed.
Appl. Sci. 2022, 12(20), 10520; https://0-doi-org.brum.beds.ac.uk/10.3390/app122010520
Submission received: 20 September 2022 / Revised: 10 October 2022 / Accepted: 14 October 2022 / Published: 18 October 2022
(This article belongs to the Special Issue Technological Advances in Seismic Data Processing and Imaging)
Browse Figures
Review Reports Versions Notes

Abstract

:
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.

1. Introduction

The Molucca Sea region is a very complex area where three main plates, the Eurasian, Australian, and Pacific plates, interact and create a collision zone known as the Molucca Sea Collision Zone [1,2,3,4,5]. The activity of these plates and the subduction of the Molucca Sea Plate are the main factors causing the high level of seismicity in this area. Tectonic activity is also increased by the Molucca Sea Plate which completely submerges the Halmahera microplate [6] and the Sangihe microplate [7], as shown in Figure S1.
The Molucca Sea Collision Zone is located between the Halmahera microplate in the east and the Sangihe microplate in the west with the Molucca Sea in the center. One collision resulted in the formation of a high ridge in the middle of the two microplates, known as the Mayu-Talaud Ridge, which is characterized by intense, shallow earthquakes and an low gravity anomaly [3,5,7,8]. The collision is ~250 km apart [4] and convergently collides with the Molucca Sea Plate [4,9] at a movement of 1.5 cm/yr [9]. The Molucca Sea Plate has an inverted U- or V-shape [2,8,10,11,12,13] and has slab subducting to the west under the Sangihe Arc and subducts to the east under the Halmahera Arc. The subduction of the Molucca Sea Plate in both eastern and western directions makes for arc volcanisms in the Halmahera and Sangihe Arcs that run parallel to the subduction zone [8,13,14] with its two Wadati–Benioff Zone active of mantle earthquakes [4,7,10]. The collision of the Halmahera and Sangihe Arcs with the subduction of the Molucca Sea Plate remains active today [15,16,17,18].
The Molucca Sea is one of the most seismically active areas in the world [19], as shown in Figure 1. There are many seismic zones in this area, such as the seismic zone concentrated in the central Molucca Sea at a depth of <60 km beneath the Mayu-Talaud Ridge, the seismic zone beneath the Sangihe and Halmahera Arcs, the seismic zone around the Gorontalo Basin, along the Minahasa peninsula, around Morotai Island, the Sulawesi Sea, the Molucca-Sorong Fault, the Bacan-Sorong Fault, the Sula-Sorong Fault, the Sulawesi slab in the northern part of Sulawesi’s north arm, and the area around the Philippine slab in the eastern part of Halmahera Island. The red dots in Figure 1 indicate shallow earthquakes; these are mostly concentrated in the central Molucca Sea area beneath the Mayu-Talaud Ridge and the Gorontalo Basin; the brown dots indicate intermediate earthquakes that are mostly concentrated beneath the Sangihe and the Halmahera Arcs; the yellowish blue dots indicate deep earthquakes that are concentrated in the Sulawesi Sea area.
Numerous studies on seismicity that relate to the structure and geometry of the Molucca Sea Plate have been carried out by previous researchers. One such study covers the two Benioff zones that lie in opposite directions and form a double subduction pattern in the Molucca Sea Plate [5,7,8,10,13]. The subducted slab under the Sangihe Arc can be identified as far down as the mantle transition zone layer at a depth of ~650 km; the subducted slab under the Halmahera Arc reaches the solid upper asthenosphere layer at a depth of ~250 km [8]. According to Cardwell et al. [10], the subducting slab to the east under the Halmahera Arc subducts at an angle of ~45° and reaches a depth of ~230 km; Hall and Spakman [20] show that it reaches depths of ~400 km, which is in line with recent research [21] that used 3D earthquake data distribution plots from ISC-EHB to identify slab configurations that were not modeled by slab2. Hutchings and Mooney [12] showed that the subducted slab in the Sangihe Arc at depths of ~250 to 350 km has a seismic gap, while the tomography image of the slab appears to be continuous.
Tomographic studies have been carried out by previous researchers, including Puspito et al. [22] who used P-phase travel time data from ISC bulletins; their work shows a high-velocity zone that forms a double subduction of the Molucca Sea Plate. The slab is subducted to the west under the Sangihe Arc and reaches to the lower mantle. Meanwhile, the slab that subducts to the east under the Halmahera Arc reaches a depth of 400 km, which is conformable with the results of the study by Widiyantoro and Hilst [23]. This research was later updated by Widiyantoro [24] who shows that the subducted slab to the west reaches the lower mantle with a folded-looking slab. Huang et al. [25] also detected a high-velocity zone in the form of bipolar subduction. Meanwhile, Zhang et al. [26] showed that the Molucca Sea Plate has a positive velocity anomaly with unique asymmetry. Fan and Zhao [27,28] showed that the subduction of the Molucca Sea Plate ~5 Ma ago propagated from south to north; they were able to detect the Molucca Sea Plate using anisotropic tomography. Zenonos et al. [29], using data from ISC-EHB, show the high-velocity Vp and Vs of the double subduction of the Molucca Sea Plate in which the high Vp reached a depth of ~700 km beneath the Sangihe Arc and ~400 km beneath the Halmahera Arc. The high Vs reached a depth of ~500 km, which is in positive agreement with a study performed by Chen et al. [21] who used only P-wave data. Zenonos et al. [30] showed that negative Vp/Vs anomalies can reach ~160 km, and positive Vp/Vs anomalies can reach depth of ~200 km. Another study by Cao et al. [31] observed the mantle flow using seismic anisotropy. Their study shows that beneath the Halmahera Arc the sub-slab mantle flow is oriented oblique to the Halmahera trench, while the sub-slab mantle flow beneath the Sangihe Arc runs in a parallel direction to the Sangihe trench.
According to previous research, tomography is generally based on data from global teleseismic earthquakes within a broad region, producing a tomographic image with a broad resolution that displays only relatively large objects. Tomography around the Molucca Sea region usually only shows the main features of the Molucca Sea Plate [22,23,26,29]. Several local researchers in the Molucca Sea area have used different methods [2,17,18]. However, their results extend only up to the layer of the Earth’s crust; the tomographic results have not shown anomalies of 3D velocity structures due to volcanic activity in the subduction zone. The research undertaken in this study made use of local and regional data from BMKG (Figure 1) to obtain more detailed features of the Molucca Sea region, especially the estimated subduction angle of the Molucca Sea Plate from south to north beneath Sulawesi’s north arm; the relation between the subduction zone of the Molucca Sea Plate subducts to the east and its volcanic activity in the Halmahera Volcanics Arc, and the depth of submerged of the double subduction of the Molucca Sea Plate in the south and north.

2. Materials and Methods

2.1. Data

The earthquake data used in this research was obtained from the BMKG network of broadband stations in the form of waveform data from three-component seismograms (BHZ, BHN, and BHE) from 1 January 2010 to 31 December 2017. A total of 32 stations were used, covering the area 4.8° S–4.25° N and 120.5° E–134.5° E (Figure 1 and Figure 2).
Considering the distribution and locations of the regional stations, earthquakes with a magnitude ≥ 4 were chosen to facilitate the identification of P- and S-wave arrival times. A total of 1647 earthquake events have been recorded in the study area. Picking was performed manually to obtain P- and S-wave arrival times using SeisGram2K version 7.0 software [32]. P- and S-wave arrival times were determined in pairs at each station; the arrival times are clearly identified and recorded by at least four stations (Figure 3). The arrival times of a total of ~17,628 P-waves and ~17,628 S-waves were identified. The highest number of P- and S- phases were recorded in the southern part of the earthquake epicenter, while a smaller number of phases were recorded in the eastern part of Seram Island. The arrival times distribution for each station are shown in Figure 2, and its histogram in Figure S2.

2.2. Methods

In NLLoc version 7.00 [33], the location algorithm adopts the inversion approach from Tarantola and Valette [34], and the event location procedure from Tarantola and Valette [34], Moser et al. [35], and Wittlinger et al. [36]. This algorithm generated the misfit function for each sample cell, estimated the probability density function (PDF), and determined the optimal hypocenter location. The location search process began by defining the search area in 3D by setting the boundaries of the x, y, z regions and then applying the 1-D AK135 velocity model [37,38]. The optimal hypocenter location for each event was determined by estimating the probability density function (PDF) using an oct-tree importance sampling.
The result of the hypocenter location from the NonLinLoc program with an equal amount of P- and S-waves used to process the travel time tomography inversion was performed using SIMULPS12 [39] to obtain the 3D structure of both Vp and Vp/Vs ratio and simultaneously hypocenters relocation by using pseudo-bending ray tracing [40], as shown in Figure S3. This program performs an inversion of Vp and Vp/Vs models from P- arrival times and S-P times directly, not Vp and Vs separately, due to differences in data quality between P- and S- arrival times [41,42]. The initial velocity model used in the inversion process is a 3D regional-global seismic velocity model (Vp) beneath Indonesia obtained from high-resolution tomographic imaging [23]. Meanwhile, the value of the Vp/Vs ratio was obtained from the Wadati Diagram (Figure S4), which is 1.76 and is consistently used for each layer.
We performed the inversion through three stages: Firstly, a damping analysis was performed to obtain the optimum damping parameter input by making variations in the data distribution, number, size, and grid nodes distance [41]. The damping value was obtained by comparing the variance model (km2/s2) and the variance data (s2) for both Vp and Vp/Vs values. The optimum value for Vp and Vp/Vs damping obtained was 150 and 130, respectively (Figure S5). This value was used to obtain a stable solution when performing the tomographic inversion and a checkerboard resolution test (CRT). Secondly, we ran a checkerboard resolution test using several grid spacings to define the best model parameterization used in this study. We decided to use a 60 × 60 km grid and a vertical grid size of 30 km after considering the results and dataset distribution. Thirdly, we performed the tomographic inversion using the damped least square where the norm of the perturbations model is weighted and combined with the misfit data squared. We set high residual data < 10 s at the early inversion. It is automatically discarded if a value greater than 10 s is found. We set the maximum number of iterations to 20 times and stopped at 15 iterations. The weighted Root Mean Square (RMS) residual value decreases each iteration until convergence. The stopping condition provided by an F-test indicates that variance reduction is insignificant after more iterations [39].

2.3. Resolution

Several methods were used to determine whether the inversion results are reliable for interpretation or not. Firstly, the Checkerboard Resolution Test (CRT) was used to test the resolution of the tomography model by making a synthetic velocity model. The CRT results showed a pattern of high and low-velocity anomalies relative to the 1D initial model, which is ±10% relative to the 1D initial model. Secondly, the Derivative Weight Sum (DWS) calculated the density of ray paths that pass through the grid nodes by giving them a certain weight. A high DWS value represents a high-resolution zone (Figure S6). Thirdly, the Ray Hit Count (RHC) was used to locate the number of ray paths that pass through the grid nodes. An increased number of ray paths represents high resolution zones that can be well-interpreted (Figure S7). Fourthly, the Diagonal Resolution Element (DRE) was used, as shown in Figure S8, which relates to the resolution matrix R (damped least-squares problems); these have an R-value between 0 and 1 [39,43,44].
The focus of this research on the vertical cross-section for the resolution test and tomography inversion consisted of four areas in the north–south cross-section: profile A-A′ (longitude 123.80°) and profile B-B′ (longitude 127.58°), and the west–east cross-section, profile C-C′ (latitude 3.26°) and profile D-D′ (latitude 1.10°). The horizontal grid size used in this study was 60 km × 60 km; the vertical grid size ranged from 30 km to a depth of 330 km, as shown in Figure 4.

3. Results

Based on the re-picking phase and distribution of available BMKG stations, we determined the location of the hypocenter using the NonLinLoc Global mode. The estimated mean error in x, y, and z directions using the covariances values are 10,198 km, 14,081 km, and 19,035 km, respectively. The distribution frequency estimated error location in the x, y, and z directions is shown in Figure S9. The distribution between the number of events and the azimuthal gap is shown in Figure S10, and the epicentral distance distribution is shown in Figure S11. We limit the value of travel time residuals in the range of −10 ≤ initial travel time residuals (itr) ≤ 10 to control the quality of the input data for the tomographic inversion process. The bad picks can be traced in unit 20 of simulps output file, as shown in Figure S12. Before the inversion, the initial travel time residual values varied in the range of −9.993 s to 9.858 s, with a mean value of −1.111 s. After inversion, the final travel time residuals varied in the range of −7.148 to 6.197 s, with a mean value of 0.156 s, forming a Gaussian distribution centered around 0 s. Tomographic inversion achieved the convergence after the 15th iteration with the weighted Root Mean Square (RMS) residual reduced by 47%, initially from 0.82312 s down to 0.43871 s (Figure S13). The RMS residual before the inversion ranged from 0.06 to 1.26 s with an average value of 0.54 s. After inversion, the range of values was down to 0.04 to 0.83 s with a mean value of 0.36 s (Figure S14) or reduced by about 33%.
The Interpretation of the tomogram becomes valid and has a high confidence level by performing a resolution test of the seismic tomography model. Therefore, we present the CRT results shown in Figure 5 for Vp and Vp/Vs of profiles A-A′, B-B′, C-C′, and D-D′ and Figure S15 horizontal section for Vp and Vp/Vs with the recovered area is delineated by dashed purple lines. The CRT positive recovery for A-A′ is beneath the Gorontalo Basin, Sulawesi’s north arm, and the Celebes Sea up to a depth of 210 km. The CRT positive recovery for the B-B′ profile is beneath the Halmahera Volcanic Arc up to a depth of 240 km; for profile C-C′ the CRT positive recovery is shown beneath the Celebes and Molucca Seas, and the Morotai Basin up to a depth of 210 km. For profile D-D′, the CRT positive recovery is beneath Sulawesi’s north arm, the Molucca Sea, and Halmahera Island up to a depth of 300 km. These results are consistent with the density distribution of P- and S-rays as revealed by the DWS, RHC, and DRE results (Shown in Figures S6–S8); the highest Vp and Vp/Vs values are listed in Table S1. The primary resolution test we used to interpret the tomographic results was the CRT, supported by the DWS, RHC, and DRE results. Hence, this research focused on the travel time tomography results in four areas in the subduction zone beneath Sulawesi’s north arm, the subduction zone beneath the Halmahera Volcano Arc, and the double subduction zone of the Molucca Sea.
Tomographic inversion results of cross-section A-A′, as shown in Figure 6 as follows:
Tomographic inversion results of cross-section B-B’, as shown in Figure 7 as follows:
Tomographic inversion results of cross-section C-C′, as shown in Figure 8 as follows:
Tomographic inversion results of cross-section D-D′, as shown in Figure 9 as follows:
Three-dimensional tomographic inversion results combine cross-section profile A-A′ and profile D-D′, as shown in Figure 10 as follows:

4. Discussion

4.1. Subduction Zone Beneath Sulawesi’s North Arm

Subducting slabs are usually indicated by high-velocity anomalies, which are caused by denser material that has a lower temperature than its surroundings [25,46,47], along with the presence of intermediate and deep earthquakes [48,49]. Based on the results of the tomographic inversion, the section A-A′ (south–north) shows a velocity anomaly of high Vp, high Vs, and low Vp/Vs at the A1 area, down to a depth of ~210 km (Figure 6). These features are interpreted as double subduction of the Molucca Sea Plate, which subducts westward beneath Sulawesi’s north arm and the Sangihe Islands with an angle of subduction from south to north, an approximate angle of ~40°. These results are consistent with previous research performed in different locations [29,30,46,50], which show that high Vp, high Vs, and low Vp/Vs can be identified as subducted slabs. It also has positive agreement the slab model of Hayes et al. [45]. In addition, most earthquakes are concentrated in a high-velocity area, these continued in the subduction zone as far down as a depth of ~500 km due to the very active movement of the subduction plate.
Three major earthquakes with magnitude ≥ 6 are observed in this area between January 2010 and June 2021. Two of these earthquakes (23 June 2020 and 19 January 2020 events) had a depth of ~100 km and were located in the South Bolaang Mongondow Regency around the Minahasa Peninsula of Sulawesi. These earthquakes are associated with the double subduction activity of the Molucca Sea Plate, which subducts to the west. Meanwhile, one shallow earthquake with a magnitude of 6.1 is located near the Gorontalo Basin and is associated with an active fault in the sea. The results of tomography inversion and hypocenters horizontally for each depth section are shown in Figure S16.
The A2 area in Figure 6 and the D2 area in Figure 9 are located in the northern part of Sulawesi’s north arm, precisely in the North Sulawesi Trench (NST) Zone. This area reveals an anomalous pattern of high Vp, high Vs, and low Vp/Vs that can be interpreted as a slab of Sulawesi that subducts to the south with the structure above being the Molucca Sea Plate. This result is consistent with previous studies conducted by Chou et al. [50], Zhao [46], and Zenonos et al. [29,30], suggesting that high Vp, high Vs, and low Vp/Vs are associated with the presence of subducted slabs. It is also consistent with Fan and Zhao [27] which identified the slab of Sulawesi from P-wave tomography and is similar to the subduction zone geometry model (yellow line) described by Hayes et al. [45]. The existence of this slab is also followed by earthquake events, even though the level of seismic activity in this cross-section is less than the earthquake activity that occurs in the A1 area of the Molucca Sea Plate zone.

4.2. Subduction Zone Beneath Halmahera Volcano Arc

The subduction zone beneath the Halmahera Volcano Arc is shown in the vertical cross-section of B-B′ (south–north) in Figure 7. The B1 area is located under three large islands in the southern part of Halmahera: Bacan, Kasiruta, and Mandioli Islands. This area is also crossed by the Bacan–Sorong Fault, which is a branch, or splay, of the Sorong Fault that passes to the east [51], and is responsible for the seismicity in this region. Based on historical records, there was a strong and destructive earthquake, magnitude 7.1, on 16 April 1963 with an intensity scale of VIII Modified Mercalli Intensity (MMI). The most recent earthquake occurred on 26 February 2020, magnitude 5.2. The earthquake shocks were felt to be quite strong and were estimated on the V MMI intensity scale. There are several volcanoes on Bacan Island, consisting of the Bibinoi Volcano in the south, and the Batusibela, Amasing, and Buku Rica (Meng) Volcanoes in the north; all of these are stratovolcanoes. The tomographic inversion results show low-Vp, low-Vs, and high-Vp/Vs anomalies, which can be interpreted as the possible source of this magmatic arc [25,52,53] as also represented by Huang et al. [25] in the same region. These volcanoes have remained dormant until the present and without any historical eruptions. However, these volcanoes may possibly erupt in the future. This possibility is supported by the flow of fluids released from partial melting (B3 area) at a depth of ~100 km above the Molucca Sea Plate that subducts to the east under Halmahera Island [54].
The B2 area shown in Figure 7 is located beneath the northern arm of Halmahera Island and belongs to the West Halmahera Regency. The vertical cross-section passes through several active volcanoes: Jailolo, Todoku-Ranu, Gamkonora, Ibu, and Tobaru. Jailolo is a stratovolcano [54] that has experienced several earthquake swarms [55,56], but whose last eruption is unknown. Young lava flows have been found on the east side of the volcano, indicating that there has been an eruption. There is also no information about the last eruption of the Todoko–Ranu Volcano, which is a caldera volcano. Gamkonora, a stratovolcano, experienced its last eruption in 2013. The Ibu Volcano, also a stratovolcano, had its last eruption in 2012. The latest update was on 16 June 2021, in the form of a volcanic ash advisory. The Tobaru Volcano, also known as the Lolodai Volcano, is a stratovolcano and its last eruption is unknown [54].
The tomographic inversion results show low-Vp, low-Vs, and high-Vp/Vs anomalies, which can be interpreted as possible sources of a magmatic arc [25,52,53]. This interpretation is also supported by the ascending flow of fluid that originated from partial melting in the B3 area at depths of ~50–70 km. The partial melting and the subducted Molucca Sea Plate (the green line in Figure 7) are shallower in this area in comparison with the B1 area. We suggest that the difference in fluid depth is probably caused by the higher activity of volcanoes in the north of Halmahera compared with volcanoes in the south. However, further and more detailed research is necessary to ascertain this.
The B3 area is located beneath the Halmahera Volcanic Arc at a depth of ~100 km and above the Molucca Sea Plate that subducts to the east. The tomographic inversion results show low-Vp, low-Vs, and high-Vp/Vs anomalies, indicating the presence of fluids that promotes partial melting. Its good agreement with the previous studies by Nakajima et al. [57] and Zhao [53], the low Vp, low Vs, and high Vp/Vs are interpreted as partial melting materials. It is observed that the upwelling process is connected with zone B1 beneath the volcanic group located on Bacan Island, and zone B2 located beneath the volcanic group on the northern arm of Halmahera. High Vp/Vs above the subduction zone can be interpreted as slab dehydration, as shown by Chou et al. [50] who stated that Vp/Vs is an effective indicator of hydration, or melting, on the mantle wedge. In studies by Lin et al. [58] and Maruyama et al. [59], hydration is indicated by low Vp, Vs, and high Vp/Vs.

4.3. Double Subduction of the Molucca Sea Plate

The double subduction zone of the Molucca Sea Plate is shown in the vertical cross-section of C-C′ (west–east), crossing the Celebes Sea, Sangihe Island, Molucca Sea, and the Morotai Basin (Figure 8); and the vertical cross-section D-D′ (west–east), crossing Sulawesi’s north arm, the Molucca Sea, and Halmahera Island (Figure 9). Based on tomographic inversion results in the C1 area, the high-Vp, high-Vs anomalies subduct to the west and reach the asthenosphere layer at a depth of ~300 km. Based on the resolution test (Figures S6–S8), the subduction pattern is shown to reach a depth of ~400 km (CRT result: not good resolved areas), while a low Vp/Vs is seen in the partially molten asthenosphere layer at a depth of ~210 km. High-Vp, high-Vs, and low-Vp/Vs anomalies subducting eastward to a depth of ~200 km are also observed. This feature can be interpreted as a double subduction of the Molucca Sea Plate, which subducts westward beneath the Sangihe Islands and eastward beneath the Morotai Basin. This is also followed by the presence of intermediate and deep earthquakes in the upper mantle layer that forms a double subduction slab pattern, as indicated by high seismic velocity [53]. These results are consistent with those of Chou et al. [50], Zhao [46], and Zenonos et al. [29,30] who suggest that high Vp, high Vs, and low Vp/Vs can be associated with the feature of subducted slab. Previous studies of P-wave tomography also identified a double subduction of the Molucca Sea Plate [9,22,23,24,27,60] in which the plate subducting to the west can reach the mantle transition layer at a depth of ~600 km. According to Fan and Zhao [27], the Molucca Sea Plate is submerged deeper in the north than in the south. These studies also indicate that the arch of double subduction of the Molucca Sea Plate has submerged into the upper mantle at a depth of ~150 km. Based on the geometry model of slab2 (green line) from Hayes et al. [45], the arch of the double subduction is closer to the surface in the uppermost mantle layer.
The C2 area is located beneath the Sangihe Trench (ST), which is a collision zone of the Sangihe continental crust that moves eastward with an overriding plate structure [26] and the Molucca Sea Plate that subducts to the west. The tomographic inversion results show that high-Vp, high-Vs, and low-Vp/Vs anomalies in the mantle crust can be interpreted as rock material that experienced micro-cracks that have since closed due to high pressure [30,61,62,63]. The C3 area is located beneath the Halmahera Trench (HT), which is a collision zone of the Halmahera continental crust that moves westward, having an overriding plate structure [26] with the Molucca Sea Plate subducting to the east. The tomographic inversion results show that low-Vp, low-Vs, and high-Vp/Vs anomalies in the crust and mantle can be interpreted as highly fractured, water-saturated rock material [63]. Our results obtained beneath the Central Ridge of Molucca Sea, indicate a clear presence of boundary between ST which has high-velocity anomalies and low Vp/Vs and HT, which has low-velocity anomalies and high Vp/Vs. This boundary also has a signature of high seismicity area. An earthquake of magnitude 6.2 occurred on 10 July 2021 on this boundary that can be interpreted as a collision zone between fault planes or fractures [10]. The boundary may represents material with different densities between the western and eastern parts of the Central Ridge of the Molucca Sea which have upraised ophiolite bodies [17,18].
In the vertical cross-section D-D′ of area D1 (Figure 9), the results of the inversion tomography show the presence of a double subduction of the Molucca Sea Plate indicated by high-Vp, high-Vs anomalies. It subducts westward and eastward, reaching a solid asthenosphere layer at a depth of ~300 and ~250 km, respectively. The Molucca Sea Plate is not well-imaged as far as the Vp/Vs ratio, and there is a discontinuity in the slab that subducts westward at depths of ~100–150 km. On the other hand, the slab that subducts to the east is well-imaged, as indicated by low Vp/Vs. The seismic pattern observed at intermediate and deep depths had a high seismic velocity, which can be associated with a double subduction slab [53]. This finding is consistent with Chou et al. [50], Zhao [46] and Zenonos et al. [29,30]; i.e., high Vp, high Vs, and low Vp/Vs are interpreted as a slab subduction. Previous studies of P-wave tomography have also identified a double subduction of the Molucca Sea Plate to depths of ~600–650 km subducting to the west, and depths of ~250–400 km subducting to the east [9,22,23,24,27,60]. The arch of double subduction in this area is in the uppermost mantle layer near the surface, which is consistent with the geometry of the slab2 model from Hayes et al. [45]. Fan and Zhao [27] showed that the arch of the Molucca Sea Plate submerged shallower in the south than in the north.
The D3 and D5 areas are located in ST and HT which is a collision zone between the continental crust and the denser, heavier oceanic crust. The ST area is the collision zone between the continental crust of Sulawesi’s north arm located above the Molucca Sea Plate that subducts to the west. The HT area is the collision zone between the continental crust of Halmahera above the Molucca Sea Plate and subducting to the east [26]. According to Puspito et al. [22], low velocity dominates the Molucca Sea collision zones and volcanic zones. The results of tomography inversion in these two areas show high-Vp, high-Vs, and low-Vp/Vs anomalies in the crust and mantle. This feature can be interpreted as rock material that experienced micro-cracks that have since closed due to high pressure [30,61,62,63]. The D4 area, located beneath the island of Tifore in the Central Molucca Sea may correspond to upraised ophiolite body [17,18] which has an active thrust fault [11,12]. The tectonic activity in this area causes many earthquakes of both small and large magnitudes. Tomographic inversion results in the Central Molucca Sea show low-Vp, low-Vs anomalies, and high Vp/Vs in the crust and mantle which indicate the presence of highly fractured brittle material and water-saturated rock [63]. In the west and east (ST and HT), we suggest that the rock material is probably not water-saturated because the micro-cracks have been covered by high pressure. However, further north (Figure 8, cross-section C-C’), the cracked and water-saturated rock materials shift to the east (there is a clear boundary in the Central Molucca Sea) and extend to HT.
Figure 10 shows a 3D visualization of the double subduction of the Molucca Sea Plate beneath Sulawesi’s north arm, the Molucca Sea, and Halmahera Island based on P-wave tomography inversion results. These results show that the Molucca Sea Plate subducts to the west beneath Sulawesi’s north arm. In perspective from south to north, the state of the slab field also dips north from the Minahasa Peninsula to the Celebes Sea, having numerous earthquake events that follow the high-Vp anomaly pattern and is interpreted as a double subduction of the Molucca Sea Plate. This is consistent with the green three-dimensional surface plot of the slab2 models of the Molucca Sea Plate [45] in Figure 6.

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 (cross-section 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://0-www-mdpi-com.brum.beds.ac.uk/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.

References

  1. Hall, R.; Nichols, G.; Ballantyne, P.; Charlton, T.; Ali, J. The Character and Significance of Basement Rocks of the Southern Molucca Sea Region. J. Southeast Asian Earth Sci. 1991, 6, 249–258. [Google Scholar] [CrossRef]
  2. McCaffrey, R.; Silver, E.A.; Raitt, R.W. Crustal Structure of the Molucca Sea Collision Zone, Indonesia. In The Tectonic and Geologic Evolution of Southeast Asian Seas and Islands; American Geophysical Union (AGU): Washington, DC, USA, 1980; pp. 161–177. [Google Scholar]
  3. Moore, G.F.; Kadarisman, D.; Evans, C.A.; Hawkins, J.W. Geology of the Talaud Islands, Molucca Sea Collision Zone, Northeast Indonesia. J. Struct. Geol. 1981, 3, 467–475. [Google Scholar] [CrossRef]
  4. Morris, J.D.; Jezek, P.A.; Hart, S.R.; Hill, J.B. The Halmahera Island Arc, Molucca Sea Collision Zone, Indonesia: A Geochemical Survey. In The Tectonic and Geologic Evolution of Southeast Asian Seas and Islands: Part 2; Geophysical Monograph Series; American Geophysical Union (AGU): Washington, DC, USA, 1983; Volume 27, pp. 373–387. [Google Scholar]
  5. Silver, E.A.; Moore, J.C. The Molucca Sea Collision Zone, Indonesia. J. Geophys. Res. Solid Earth 1978, 83, 1681–1691. [Google Scholar] [CrossRef]
  6. Hall, R.; Audley-Charles, M.G.; Banner, F.T.; Hidayat, S.; Tobing, S.L. Basement Rocks of the Halmahera Region, Eastern Indonesia: A Late Cretaceous–Early Tertiary Arc and Fore-Arc. J. Geol. Soc. London. 1988, 145, 65–84. [Google Scholar] [CrossRef]
  7. Hamilton, W.B. Tectonics of the Indonesian Region; U.S. Govt. Print. Off.: Washington, DC, USA, 1979. [Google Scholar]
  8. Hatherton, T.; Dickinson, W.R. The Relationship between Andesitic Volcanism and Seismicity in Indonesia, the Lesser Antilles, and Other Island Arcs. J. Geophys. Res. 1969, 74, 5301–5310. [Google Scholar] [CrossRef]
  9. Hafkenscheid, E.; Buiter, S.J.; Wortel, M.J.; Spakman, W.; Bijwaard, H. Modelling the Seismic Velocity Structure beneath Indonesia: A Comparison with Tomography. Tectonophysics 2001, 333, 35–46. [Google Scholar] [CrossRef]
  10. Cardwell, R.K.; Isaacks, B.L.; Karig, D.E. The Spatial Distribution of Earthquakes, Focal Mechanism Solutions, and Subducted Lithosphere in the Philippine and Northeastern Indonesian Islands. In The Tectonic and Geologic Evolution of Southeast Asian Seas and Islands; Geophys. Monogr., Ser.; Hayes, D.E., Ed.; AGU: Washington, DC, USA, 1980; Volume 23, pp. 1–35. [Google Scholar]
  11. McCaffrey, R. Lithospheric Deformation within the Molucca Sea Arc-Arc Collision: Evidence from Shallow and Intermediate Earthquake Activity. J. Geophys. Res. 1982, 87, 3663. [Google Scholar] [CrossRef]
  12. Hutchings, S.J.; Mooney, W.D. The Seismicity of Indonesia and Tectonic Implications. Geochem. Geophys. Geosyst. 2021, 22, e2021GC009812. [Google Scholar] [CrossRef]
  13. Hall, R. Plate Boundary Evolution in the Halmahera Region, Indonesia. Tectonophysics 1987, 144, 337–352. [Google Scholar] [CrossRef]
  14. Milsom, J.; Masson, D.; Nicols, G. Three Trench Endings in Eastern Indonesia. Mar. Geol. 1992, 104, 227–241. [Google Scholar] [CrossRef]
  15. Sufni Hakim, A.; Hall, R. Tertiary Volcanic Rocks from the Halmahera Arc, Eastern Indonesia. J. Southeast Asian Earth Sci. 1991, 6, 271–287. [Google Scholar] [CrossRef]
  16. Moore, G.F.; Silver, E.A. Collision Processes in the Northern Molucca Sea. In The Tectonic and Geologic Evolution of Southeast Asian Seas and Islands: Part 2; Geophys. Monogr., Ser.; Hayes, D.E., Ed.; AGU: Washington, DC, USA, 1983; Volume 27, pp. 360–372. [Google Scholar]
  17. Widiwijayanti, C.; Tiberi, C.; Deplus, C.; Diament, M.; Mikhailov, V.; Louat, R. Geodynamic Evolution of the Northern Molucca Sea Area (Eastern Indonesia) Constrained by 3-D Gravity Field Inversion. Tectonophysics 2004, 386, 203–222. [Google Scholar] [CrossRef]
  18. Widiwijayanti, C.; Mikhailov, V.; Diament, M.; Deplus, C.; Louat, R.; Tikhotsky, S.; Gvishiani, A. Structure and Evolution of the Molucca Sea Area: Constraints Based on Interpretation of a Combined Sea-Surface and Satellite Gravity Dataset. Earth Planet. Sci. Lett. 2003, 215, 135–150. [Google Scholar] [CrossRef]
  19. Hédervári, P.; Papp, Z. Seismicity Maps of the Indonesian Region. Tectonophysics 1981, 76, 131–148. [Google Scholar] [CrossRef]
  20. Hall, R.; Spakman, W. Mantle Structure and Tectonic History of SE Asia. Tectonophysics 2015, 658, 14–45. [Google Scholar] [CrossRef] [Green Version]
  21. Chen, P.F.; Chien, M.; Bina, C.R.; Yen, H.Y.; Antonio Olavere, E. Evidence of an East-Dipping Slab beneath the Southern End of the Philippine Trench (1°N–6°N) as Revealed by ISC-EHB. J. Asian Earth Sci. X 2020, 4, 100034. [Google Scholar] [CrossRef]
  22. Puspito, N.T.; Yamanaka, Y.; Miyatake, T.; Shimazaki, K.; Hirahara, K. Three-Dimensional P-Wave Velocity Structure beneath the Indonesian Region. Tectonophysics 1993, 220, 175–192. [Google Scholar] [CrossRef]
  23. Widiyantoro, S.; Hilst, R. Mantle Structure beneath Indonesia Inferred from High-Resolution Tomographic Imaging. Geophys. J. Int. 1997, 130, 167–182. [Google Scholar] [CrossRef] [Green Version]
  24. Widiyantoro, S. Complex Morphology of Subducted Lithosphere in the Mantle below the Molucca Collision Zone from Non-Linear Seismic Tomography. ITB J. Eng. Sci. 2003, 35, 1–10. [Google Scholar] [CrossRef] [Green Version]
  25. Huang, Z.; Zhao, D.; Wang, L. P Wave Tomography and Anisotropy beneath Southeast Asia: Insight into Mantle Dynamics. J. Geophys. Res. Solid Earth 2015, 120, 5154–5174. [Google Scholar] [CrossRef]
  26. Zhang, Q.; Guo, F.; Zhao, L.; Wu, Y. Geodynamics of Divergent Double Subduction: 3-D Numerical Modeling of a Cenozoic Example in the Molucca Sea Region, Indonesia. J. Geophys. Res. Solid Earth 2017, 122, 3977–3998. [Google Scholar] [CrossRef]
  27. Fan, J.; Zhao, D. Evolution of the Southern Segment of the Philippine Trench: Constraints from Seismic Tomography. Geochem. Geophys. Geosyst. 2018, 19, 4612–4627. [Google Scholar] [CrossRef]
  28. Fan, J.; Zhao, D. P-Wave Anisotropic Tomography of the Central and Southern Philippines. Phys. Earth Planet. Inter. 2019, 286, 154–164. [Google Scholar] [CrossRef]
  29. Zenonos, A.; De Siena, L.; Widiyantoro, S.; Rawlinson, N. P and S Wave Travel Time Tomography of the SE Asia-Australia Collision Zone. Phys. Earth Planet. Inter. 2019, 293, 106267. [Google Scholar] [CrossRef]
  30. Zenonos, A.; De Siena, L.; Widiyantoro, S.; Rawlinson, N. Direct Inversion of S-P Differential Arrival Times for Ratio in SE Asia. J. Geophys. Res. Solid Earth 2020, 125, e2019JB019152. [Google Scholar] [CrossRef]
  31. Cao, L.; He, X.; Zhao, L.; Lü, C.; Hao, T.; Zhao, M.; Qiu, X. Mantle Flow Patterns Beneath the Junction of Multiple Subduction Systems Between the Pacific and Tethys Domains, SE Asia: Constraints From SKS-Wave Splitting Measurements. Geochem. Geophys. Geosyst. 2021, 22, e2021GC009700. [Google Scholar] [CrossRef]
  32. Lomax, A.; Michelini, A. Mwpd: A Duration-Amplitude Procedure for Rapid Determination of Earthquake Magnitude and Tsunamigenic Potential from P Waveforms. Geophys. J. Int. 2009, 176, 200–214. [Google Scholar] [CrossRef] [Green Version]
  33. Lomax, A.; Virieux, J.; Volant, P.; Berge-Thierry, C. Probabilistic Earthquake Location in 3D and Layered Models; Springer: Dordrecht, The Netherlands, 2000; pp. 101–134. [Google Scholar]
  34. Tarantola, A.; Valette, B. Inverse Problems = Quest for Information. J. Geophys. 1982, 50, 159–170. [Google Scholar]
  35. Moser, T.J.; van Eck, T.; Nolet, G. Hypocenter Determination in Strongly Heterogeneous Earth Models Using the Shortest Path Method. J. Geophys. Res. 1992, 97, 6563. [Google Scholar] [CrossRef]
  36. Wittlinger, G.; Herquel, G.; Nakache, T. Earthquake Location in Strongly Heterogeneous Media. Geophys. J. Int. 1993, 115, 759–777. [Google Scholar] [CrossRef] [Green Version]
  37. Kennett, B.L.N.; Engdahl, E.R. Traveltimes for Global Earthquake Location and Phase Identification. Geophys. J. Int. 1991, 105, 429–465. [Google Scholar] [CrossRef]
  38. Kennett, B.L.N.; Engdahl, E.R.; Buland, R. Constraints on Seismic Velocities in the Earth from Traveltimes. Geophys. J. Int. 1995, 122, 108–124. [Google Scholar] [CrossRef] [Green Version]
  39. Evans, J.R.; Eberhart-Phillips, D.; Thurber, C.H. User’s Manual for SIMULPS12 for Imaging vp and vp/vs; a Derivative of the “Thurber” Tomographic Inversion SIMUL3 for Local Earthquakes and Explosions; Open-File Rep. Open-File Report 94-431; U.S. Geological Survey: Reston, VI, USA, 1994; p. 101. [Google Scholar] [CrossRef] [Green Version]
  40. Um, J.; Thurber, C. A Fast Algorithm for Two-Point Seismic Ray Tracing. Bull. Seismol. Soc. Am. 1987, 77, 972–986. [Google Scholar] [CrossRef]
  41. Eberhart-Phillips, D. Three-Dimensional Velocity Structure in Northern California Coast Ranges from Inversion of Local Earthquake Arrival Times. Bull. Seismol. Soc. Am. 1986, 76, 1025–1052. [Google Scholar] [CrossRef]
  42. Thurber, C.H. Local Earthquake Tomography: Velocities and Vp/Vs-Theory; Chapman and Hall: London, UK, 1993; ISBN 0412371901. [Google Scholar]
  43. Nugraha, A. Tomografi Seismik, 1st ed.; ITB Press: Bandung, Indonesia, 2017; ISBN 978-602-5417-48-1. [Google Scholar]
  44. Rawlinson, N.; Sambridge, M. Seismic Traveltime Tomography of The Crust and Lithosphere; Elsevier: Amsterdam, The Netherlands, 2003; pp. 81–198. [Google Scholar]
  45. Hayes, G.P.; Moore, G.L.; Portner, D.E.; Hearne, M.; Flamme, H.; Furtney, M.; Smoczyk, G.M. Slab2, a Comprehensive Subduction Zone Geometry Model. Science 2018, 362, 58–61. [Google Scholar] [CrossRef] [PubMed]
  46. Zhao, D. Multiscale Seismic Tomography; Springer: Japan, Tokyo, 2015; ISBN 978-4-431-55359-5. [Google Scholar]
  47. Zhao, D.; Yamashita, K.; Toyokuni, G. Tomography of the 2016 Kumamoto Earthquake Area and the Beppu-Shimabara Graben. Sci. Rep. 2018, 8, 15488. [Google Scholar] [CrossRef] [Green Version]
  48. Hasegawa, A.; Umino, N.; Takagi, A. Double-Planed Structure of the Deep Seismic Zone in the Northeastern Japan Arc. Tectonophysics 1978, 47, 43–58. [Google Scholar] [CrossRef]
  49. Hasegawa, A.; Umino, N.; Takagi, A. Double-Planed Deep Seismic Zone and Upper-Mantle Structure in the Northeastern Japan Arc. Geophys. J. Int. 1978, 54, 281–296. [Google Scholar] [CrossRef] [Green Version]
  50. Chou, H.-C.; Kuo, B.-Y.; Chiao, L.-Y.; Zhao, D.; Hung, S.-H. Tomography of the Westernmost Ryukyu Subduction Zone and the Serpentinization of the Fore-Arc Mantle. J. Geophys. Res. Solid Earth 2009, 114, 12301. [Google Scholar] [CrossRef] [Green Version]
  51. Hall, R.; Ali, J.R.; Anderson, C.D.; Hall, R.; Ali, J.R.; Anderson, C.D. Cenozoic Motion of the Philippine Sea Plate: Palaeomagnetic Evidence from Eastern Indonesia. Tecto 1995, 14, 1117–1132. [Google Scholar] [CrossRef]
  52. Kim, K.-H.; Chiu, J.-M.; Pujol, J.; Chen, K.-C.; Huang, B.-S.; Yeh, Y.-H.; Shen, P. Three-Dimensional VP and VS Structural Models Associated with the Active Subduction and Collision Tectonics in the Taiwan Region. Geophys. J. Int. 2005, 162, 204–220. [Google Scholar] [CrossRef] [Green Version]
  53. Zhao, D. Multiscale Seismic Tomography and Mantle Dynamics. Gondwana Res. 2009, 15, 297–323. [Google Scholar] [CrossRef]
  54. Siebert, L.; Simkin, T.; Kimberly, P. Volcanoes of the World; Smithsonian Institution: Washington, DC, USA, 2010; ISBN 0520268776. [Google Scholar]
  55. Nugraha, A.D.; Shiddiqi, H.A.; Widiyantoro, S.; Puspito, N.T.; Triyoso, W.; Wiyono, S.; Daryono; Wandono; Rosalia, S. Hypocenter Relocation of Earthquake Swarm in West Halmahera, North Molucca Region, Indonesia by Using Double-Difference Method and 3D Seismic Velocity Structure. IOP Conf. Ser. Earth Environ. Sci. 2017, 62, 012053. [Google Scholar] [CrossRef] [Green Version]
  56. Nugraha, A.D.; Supendi, P.; Widiyantoro, S.; Daryono; Wiyono, S. Hypocenter Relocation of Earthquake Swarm around Jailolo Volcano, North Molucca, Indonesia Using the BMKG Network Data: Time Periods of 27 September–10 October 2017. AIP Conf. Proc. 2018, 1987, 020093. [Google Scholar] [CrossRef]
  57. Nakajima, J.; Matsuzawa, T.; Hasegawa, A.; Zhao, D. Three-Dimensional Structure of Vp, Vs, and Vp/Vs beneath Northeastern Japan: Implications for Arc Magmatism and Fluids. J. Geophys. Res. Solid Earth 2001, 106, 21843–21857. [Google Scholar] [CrossRef]
  58. Lin, J.Y.; Hsu, S.K.; Sibuet, J.C. Melting Features along the Ryukyu Slab Tear, beneath the Southwestern Okinawa Trough. Geophys. Res. Lett. 2004, 31, 20862. [Google Scholar] [CrossRef] [Green Version]
  59. Maruyama, S.; Hasegawa, A.; Santosh, M.; Kogiso, T.; Omori, S.; Nakamura, H.; Kawai, K.; Zhao, D. The Dynamics of Big Mantle Wedge, Magma Factory, and Metamorphic-Metasomatic Factory in Subduction Zones. Gondwana Res. 2009, 16, 414–430. [Google Scholar] [CrossRef]
  60. Amaru, M.L. Global Travel Time Tomography with 3-D Reference Models. Geol. Ultraiectina 2007, 274, 81–88. [Google Scholar]
  61. Wagner, L.S.; Anderson, M.L.; Jackson, J.M.; Beck, S.L.; Zandt, G. Seismic Evidence for Orthopyroxene Enrichment in the Continental Lithosphere. Geology 2008, 36, 935–938. [Google Scholar] [CrossRef] [Green Version]
  62. Zheng, Y.; Lay, T. Low Vp/Vs Ratios in the Crust and Upper Mantle beneath the Sea of Okhotsk Inferred from Teleseismic PMP, SMP, and SMS Underside Reflections from the Moho. J. Geophys. Res. Solid Earth 2006, 111, 1305. [Google Scholar] [CrossRef] [Green Version]
  63. Bariş, S.; Nakajima, J.; Hasegawa, A.; Honkura, Y.; Ito, A.; Balamir, S. Three-Dimensional Structure of Vp, Vs and Vp/Vs in the Upper Crust of the Marmara Region, NW Turkey. Earth Planets Sp. 2005, 57, 1019–1038. [Google Scholar] [CrossRef]
Applsci 12 10520 g001 550
Figure 1. The horizontal distribution of seismic map of the research location around the Molucca Collision Zone. SLW: Sulawesi island; ML: Molucca island; PA: Papua island; HLM: Halmahera island; SNA: Sulawesi’s north arm; SS: Sulawesi slab; PS: Philippine slab; HT and ST: Halmahera and Sangihe Trench; HA and SA: Halmahera and Sangihe volcanoes arc; CS: Celebes sea; MP: Minahassa peninsula; GB: Gorontalo Basin; MS: Molucca Sea; MTR: Mayu-Talaud Ridge; BSF, MSF, and SSF: Bacan, Molucca, and Sula Sorong Faults, respectively.
Figure 1. The horizontal distribution of seismic map of the research location around the Molucca Collision Zone. SLW: Sulawesi island; ML: Molucca island; PA: Papua island; HLM: Halmahera island; SNA: Sulawesi’s north arm; SS: Sulawesi slab; PS: Philippine slab; HT and ST: Halmahera and Sangihe Trench; HA and SA: Halmahera and Sangihe volcanoes arc; CS: Celebes sea; MP: Minahassa peninsula; GB: Gorontalo Basin; MS: Molucca Sea; MTR: Mayu-Talaud Ridge; BSF, MSF, and SSF: Bacan, Molucca, and Sula Sorong Faults, respectively.
Applsci 12 10520 g001
Applsci 12 10520 g002 550
Figure 2. Arrival times distribution of P-and S-wave for all stations. Colored inverted triangles denote the number of phases per station.
Figure 2. Arrival times distribution of P-and S-wave for all stations. Colored inverted triangles denote the number of phases per station.
Applsci 12 10520 g002
Applsci 12 10520 g003 550
Figure 3. Example of the 22 November 2017 event waveform from a three-component seismogram (BHZ, BHN, and BHE) which was picked using SeisGram2K version 7.0 software [32]. The blue, green, and red colors from each station represent the vertical component (BHZ), the north–south component (BHN), and the east–west component (BHE). From top to bottom, the recorded seismograms are shown at the KMSI, TMSI, GTOI, SMSI, SANI, TNTI, LBMI, and TOLI2 stations. The vertical P and S lines show the picked arrival times for P- and S-waves, while Ts-Tp is the time difference between the arrival times of the S- and P-waves.
Figure 3. Example of the 22 November 2017 event waveform from a three-component seismogram (BHZ, BHN, and BHE) which was picked using SeisGram2K version 7.0 software [32]. The blue, green, and red colors from each station represent the vertical component (BHZ), the north–south component (BHN), and the east–west component (BHE). From top to bottom, the recorded seismograms are shown at the KMSI, TMSI, GTOI, SMSI, SANI, TNTI, LBMI, and TOLI2 stations. The vertical P and S lines show the picked arrival times for P- and S-waves, while Ts-Tp is the time difference between the arrival times of the S- and P-waves.
Applsci 12 10520 g003
Applsci 12 10520 g004 550
Figure 4. The horizontal and vertical seismic map of the research location in the Molucca Collision Zone. The bold black lines A-A′ (south–north), B-B′ (south–north), C-C′ (west–east), D-D′ (west–east) denote the areas of the vertical cross-section tomography. The black plus sign represents the distribution of the grid nodes for tomographic inversion; the colored dots indicate the distribution of earthquakes; the blue triangles are BMKG stations, and the black triangles signify volcanoes.
Figure 4. The horizontal and vertical seismic map of the research location in the Molucca Collision Zone. The bold black lines A-A′ (south–north), B-B′ (south–north), C-C′ (west–east), D-D′ (west–east) denote the areas of the vertical cross-section tomography. The black plus sign represents the distribution of the grid nodes for tomographic inversion; the colored dots indicate the distribution of earthquakes; the blue triangles are BMKG stations, and the black triangles signify volcanoes.
Applsci 12 10520 g004
Applsci 12 10520 g005a 550Applsci 12 10520 g005b 550
Figure 5. The vertical cross-section of CRT profiles. Profile A-A′ (south–north), B-B′ (south–north), C-C′ (west–east), and D-D′ (west–east) are vertical cross-sections of the CRT recovery model from the inversion process. The purple dashed line is the positive CRT recovery area for tomographic interpretation. The first column is Vp perturbation, the second column is Vp/Vs perturbation, and the red triangles signify volcanoes.
Figure 5. The vertical cross-section of CRT profiles. Profile A-A′ (south–north), B-B′ (south–north), C-C′ (west–east), and D-D′ (west–east) are vertical cross-sections of the CRT recovery model from the inversion process. The purple dashed line is the positive CRT recovery area for tomographic interpretation. The first column is Vp perturbation, the second column is Vp/Vs perturbation, and the red triangles signify volcanoes.
Applsci 12 10520 g005aApplsci 12 10520 g005b
Applsci 12 10520 g006 550
Figure 6. Vertical cross-section of A-A′ (south–north) of the Molucca Sea Plate, crossing the Minahasa Peninsula, Sulawesi’s north arm, and the Celebes Sea. The red color indicates low velocity and low Vp/Vs ratio, and the blue color indicates high velocity and high Vp/Vs ratio. The A1 and A2 are areas for the interpretation of tomographic inversion results. The black dots signify earthquakes within a radius of 30 km ≤ A-A′ ≤ 30 km, and the yellow stars are significant earthquakes. The green line and the green three-dimensional surface plot are the slab2 models of the Molucca Sea Plate. The yellow line is the slab2 model of Sulawesi [45]. NST is the North Sulawesi Trench. The red triangles signify volcanoes. The blurred area is an area that is not well-resolved.
Figure 6. Vertical cross-section of A-A′ (south–north) of the Molucca Sea Plate, crossing the Minahasa Peninsula, Sulawesi’s north arm, and the Celebes Sea. The red color indicates low velocity and low Vp/Vs ratio, and the blue color indicates high velocity and high Vp/Vs ratio. The A1 and A2 are areas for the interpretation of tomographic inversion results. The black dots signify earthquakes within a radius of 30 km ≤ A-A′ ≤ 30 km, and the yellow stars are significant earthquakes. The green line and the green three-dimensional surface plot are the slab2 models of the Molucca Sea Plate. The yellow line is the slab2 model of Sulawesi [45]. NST is the North Sulawesi Trench. The red triangles signify volcanoes. The blurred area is an area that is not well-resolved.
Applsci 12 10520 g006
Applsci 12 10520 g007 550
Figure 7. Vertical cross-section of B-B’ (south–north) crossing beneath the Halmahera Volcanic Arc. The red color indicates low velocity and low Vp/Vs ratio. The blue color indicates high velocity and high Vp/Vs ratio. The B1, B2, and B3 are areas for the interpretation of tomographic inversion results. The green line and the green three-dimensional surface plot are the slab2 models of the Molucca Sea Plate [45]. The black dots signify earthquakes within a radius of 30 km ≤ B-B′ ≤ 30 km. SeT: Seram Trench; HT: Halmahera Trench; Bi: Bibinoi Volcano; Bs: Batusibela Volcano; Am: Amasing Volcano; Br: Buku-Rica Volcano; Rr: Rogi-Rogi Volcano; Ja: Jailolo Volcano; To: Todoko-Ranu Volcano; Ga: Gamkonora Volcano; Ib: Ibu Volcano; Tb: Tobaru Volcano. The blurred area is an area that is not well-resolved.
Figure 7. Vertical cross-section of B-B’ (south–north) crossing beneath the Halmahera Volcanic Arc. The red color indicates low velocity and low Vp/Vs ratio. The blue color indicates high velocity and high Vp/Vs ratio. The B1, B2, and B3 are areas for the interpretation of tomographic inversion results. The green line and the green three-dimensional surface plot are the slab2 models of the Molucca Sea Plate [45]. The black dots signify earthquakes within a radius of 30 km ≤ B-B′ ≤ 30 km. SeT: Seram Trench; HT: Halmahera Trench; Bi: Bibinoi Volcano; Bs: Batusibela Volcano; Am: Amasing Volcano; Br: Buku-Rica Volcano; Rr: Rogi-Rogi Volcano; Ja: Jailolo Volcano; To: Todoko-Ranu Volcano; Ga: Gamkonora Volcano; Ib: Ibu Volcano; Tb: Tobaru Volcano. The blurred area is an area that is not well-resolved.
Applsci 12 10520 g007
Applsci 12 10520 g008 550
Figure 8. Vertical cross-section of C-C′ (west–east) crossing the Celebes Sea, Sangihe Island, Molucca Sea, and Morotai Basin. The red color indicates low velocity and low Vp/Vs ratio. The blue color indicates high velocity and high Vp/Vs ratio. The C1, C2, and C3 are areas for the interpretation of tomographic inversion results. The black dots signify earthquakes within a radius of 30 km ≤ C-C′ ≤ 30 km, and the yellow star is a significant earthquake. The green line and the green three-dimensional surface plot are the slab2 models of the Molucca Sea Plate, and the brown line is the slab2 model of the Philippines [45]. ST: Sangihe Trench; HT: Halmahera Trench; the blurred area is an area that is not well-resolved.
Figure 8. Vertical cross-section of C-C′ (west–east) crossing the Celebes Sea, Sangihe Island, Molucca Sea, and Morotai Basin. The red color indicates low velocity and low Vp/Vs ratio. The blue color indicates high velocity and high Vp/Vs ratio. The C1, C2, and C3 are areas for the interpretation of tomographic inversion results. The black dots signify earthquakes within a radius of 30 km ≤ C-C′ ≤ 30 km, and the yellow star is a significant earthquake. The green line and the green three-dimensional surface plot are the slab2 models of the Molucca Sea Plate, and the brown line is the slab2 model of the Philippines [45]. ST: Sangihe Trench; HT: Halmahera Trench; the blurred area is an area that is not well-resolved.
Applsci 12 10520 g008
Applsci 12 10520 g009 550
Figure 9. Vertical cross-section of D-D′ (west–east) crossing Sulawesi’s north arm, the Molucca Sea, and Halmahera Island. The red color indicates low-velocity areas and low Vp/Vs ratio, and the blue color indicates high-velocity areas and high Vp/Vs ratio. The D1, D2, D3, D4, and D5 are areas for the interpretation of tomographic inversion results. The black dots signify earthquakes within a radius of 30 km ≤ D-D′ ≤ 30 km, and the yellow star is the magnitude of a significant earthquake. The green line and the green three-dimensional surface plot are the slab2 models of the Molucca Sea Plate; the yellow line is the slab2 model of Sulawesi [45]. Thus, So: Soputan Volcano; ST: Sangihe Trench; the reserve arrow is Tifore island; HT: Halmahera Trench; Ja: Jailolo Volcano; To: Todoko-Ranu Volcano; Wa: Wato-Wato Volcano; the blurred area is an area that is not well-resolved.
Figure 9. Vertical cross-section of D-D′ (west–east) crossing Sulawesi’s north arm, the Molucca Sea, and Halmahera Island. The red color indicates low-velocity areas and low Vp/Vs ratio, and the blue color indicates high-velocity areas and high Vp/Vs ratio. The D1, D2, D3, D4, and D5 are areas for the interpretation of tomographic inversion results. The black dots signify earthquakes within a radius of 30 km ≤ D-D′ ≤ 30 km, and the yellow star is the magnitude of a significant earthquake. The green line and the green three-dimensional surface plot are the slab2 models of the Molucca Sea Plate; the yellow line is the slab2 model of Sulawesi [45]. Thus, So: Soputan Volcano; ST: Sangihe Trench; the reserve arrow is Tifore island; HT: Halmahera Trench; Ja: Jailolo Volcano; To: Todoko-Ranu Volcano; Wa: Wato-Wato Volcano; the blurred area is an area that is not well-resolved.
Applsci 12 10520 g009
Applsci 12 10520 g010 550
Figure 10. The 3D visualization of the Molucca Sea Plate, combining vertical cross-section profile A-A′ (south–north) and profile D-D’ (west–east) from the P-wave tomography inversion results. The earthquake events follow the high Vp anomaly pattern. The red color indicates low-velocity areas and low Vp/Vs ratio, and the blue color indicates high velocity and high Vp/Vs ratio. The white dots signify earthquakes; the yellow stars signify earthquakes with significant magnitudes. The red triangle is a volcanic arc located in Sangihe and Halmahera.
Figure 10. The 3D visualization of the Molucca Sea Plate, combining vertical cross-section profile A-A′ (south–north) and profile D-D’ (west–east) from the P-wave tomography inversion results. The earthquake events follow the high Vp anomaly pattern. The red color indicates low-velocity areas and low Vp/Vs ratio, and the blue color indicates high velocity and high Vp/Vs ratio. The white dots signify earthquakes; the yellow stars signify earthquakes with significant magnitudes. The red triangle is a volcanic arc located in Sangihe and Halmahera.
Applsci 12 10520 g010
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Rachman, G.; Santosa, B.J.; Nugraha, A.D.; Rohadi, S.; Rosalia, S.; Zulfakriza, Z.; Sungkono, S.; Sahara, D.P.; Muttaqy, F.; Supendi, P.; et al. Seismic Structure Beneath the Molucca Sea Collision Zone from Travel Time Tomography Based on Local and Regional BMKG Networks. Appl. Sci. 2022, 12, 10520. https://0-doi-org.brum.beds.ac.uk/10.3390/app122010520

AMA Style

Rachman G, Santosa BJ, Nugraha AD, Rohadi S, Rosalia S, Zulfakriza Z, Sungkono S, Sahara DP, Muttaqy F, Supendi P, et al. Seismic Structure Beneath the Molucca Sea Collision Zone from Travel Time Tomography Based on Local and Regional BMKG Networks. Applied Sciences. 2022; 12(20):10520. https://0-doi-org.brum.beds.ac.uk/10.3390/app122010520

Chicago/Turabian Style

Rachman, Gazali, Bagus Jaya Santosa, Andri Dian Nugraha, Supriyanto Rohadi, Shindy Rosalia, Zulfakriza Zulfakriza, Sungkono Sungkono, David P. Sahara, Faiz Muttaqy, Pepen Supendi, and et al. 2022. "Seismic Structure Beneath the Molucca Sea Collision Zone from Travel Time Tomography Based on Local and Regional BMKG Networks" Applied Sciences 12, no. 20: 10520. https://0-doi-org.brum.beds.ac.uk/10.3390/app122010520

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop