Local seismotectonic analysis of the July 2019 Molucca Sea earthquake sequence based on moment tensor solutions

We investigate the local seismotectonic of the Molucca Sea area using moment tensor calculations for the earthquakes that occurred in July 2019 at a depth of 10–55 km. The mainshock of Mw 6.8 occurred on July 7, followed by aftershocks until July 18, with magnitudes ranging from Mw 4.6 to Mw 5.8. Moment tensor solutions are calculated by applying Isolated Asperities (ISOLA) software using the full waveform data recorded at regional seismic stations. The analyzed frequency bands used in this study are 0.01–0.03 Hz and 0.04–0.05 Hz for the event with Mw ≥ 5 and Mw < 5, respectively. We provide validations of new moment tensor solutions for Mw < 5 events in the Molucca Sea region for the period during the earthquake sequence. The results show that thrust and oblique faults are dominant during this event, which indicate a compressional stress of divergent double subduction (DDS) of the Sangihe and Halmahera arcs. Only one full moment tensor solution reveals the normal fault mechanism, which may indicate the manifestation of strain release of compressional stress in the surrounding area. Furthermore, these results also support the previous studies suggesting that the Talaud-Mayu Ridge located in the middle of the Molucca Sea has developed as a consequence of the transpressional tectonic activity.


Introduction
The Molucca Sea, an area well known for its tectonic complexity, is a region with high seismicity, located between the eastern arm of Sulawesi Island and Halmahera Island in the eastern part of Indonesia (see Fig. 1a). This region lies within a complex triple junction between three major plates: the southeastern section of the Eurasian Plate, the Philippine Sea Plate, and the Australian Plate (Rangin et al. 1999;Bird 2003). The Molucca Sea Plate, one of the region's microplates, is bounded by the convergence of the Eurasian Plate in the western part and the Philippine Sea Plate in the eastern part. The Eurasian Plate and the Philippine Sea Plate move to the east and west direction, respectively (e.g., Hall and Smyth 2008;Waltham et al. 2008). On the other hand, the Molucca Sea Plate is subducted under the arc-arc collision between Sangihe in the west and Halmahera in the east. This region is the only unique example of an active arcarc collision on earth that consumes an oceanic basin via subduction in a different direction. It is formed by the convergence of the Sangihe and Halmahera arcs into the Molucca Sea Plate in eastward and westward directions, respectively (see Fig. 1b). This complexity is known as the Molucca Sea Collision Zone. This double subduction system causes the Molucca Sea Plate to be consumed by the Eurasian Plate in the west and the Philippines Sea Plate in the east, and the upliftment of the Talaud-Mayu Ridge at the center of the Molucca Sea. Furthermore, the accretional complex has been thrust toward the forearc Open Access *Correspondence: aditya.dwi.prasetio@lipi.go.id; adityadprasetio@gmail.com 1 Research Center for Physics -LIPI, Kawasan PUSPIPTEK Serpong, South Tangerang, Indonesia Full list of author information is available at the end of the article  Supendi et al. (2020) shown by colored dots. The red rectangle in the inset presents the coverage area of the seismicity map. Black lines represent the plate boundaries. The yellow rectangle represents the location of this study. The blue star represents the largest recent earthquake in the Molucca Sea region (Mw 7.1), which occurred on November 14, 2014. b Illustration of A-B cross section of Molucca Sea that describes the structure of arc-arc collision zone (modified from Hall and Smyth 2008;Waltham et al. 2008)) sedimentary basin (Silver and Moore 1978;Hamilton 1979;McCaffrey et al. 1980;Hall 1987Hall , 1997Hall and Blundell 1996;Hinschberger et al. 2005;Hall and Smyth 2008). Such condition makes the tectonic setting of Molucca unique compared to other regions.
Several studies have successfully identified the geological structure of the subsurface surrounding the Molucca Sea, as evidenced by Bouguer gravity anomaly, seismicity, and geological studies (Bemmelen 1949;Setyanta and Setiadi 2011). The geological studies show that the tectonic setting of the Molucca Sea from west to east consists of a volcanic arc ridge in northern Sulawesi (the Sangihe Ridge), a non-volcanic arc ridge in Talaud (the Talaud-Mayu Ridge), a volcanic arc ridge in northern Halmahera, and a non-volcanic arc ridge in eastern Halmahera. Several geophysical investigations have also been conducted in the Molucca Sea to get a better understanding of the area's tectonic setting. For example, the double subduction of the Molucca Sea Plate has been modeled in 3D by Zhang et al. (2017). Their results suggest that the subduction process occurs asymmetrically on the east and west sides of the slab. In addition, their model agrees well with the P-S wave travel-time tomography study of Zenonos et al. (2019), suggesting that the Sangihe slab penetrates to the lower mantle while the Halmahera slab only reaches a depth of 400-500 km below sea level. Widiwijayanti et al. (2003Widiwijayanti et al. ( , 2004 used satellite gravity data to determine the subsurface structure of the northern Molucca Sea. They identified the central ridge as an active deformed region based on the large density variation found in their research. Furthermore, the geodetic investigation suggests that the plate convergence rate between the Sangihe and Halmahera subduction systems in the Molucca Sea is estimated to be about 80 mm/year (Rangin et al. 1999). The latest finding, based on GPS observation of the Mw 7.1 Molucca Sea earthquake on November 14, 2014, Gunawan et al. (2016) suggests that the earthquake caused a coseismic slip of ± 36 cm around the hypocenter located in the central ridge. In other previous seismic studies conducted in the Molucca Sea, hypocenter relocation was investigated to enhance the location accuracy of earthquake events Shiddiqi et al. 2015;Yulianto et al. 2017;Supendi et al. 2020). These studies apply the double-difference method, where the distance of two events is smaller than the distances to measuring stations, implying that they have similar ray paths (Waldhauser and Ellsworth 2000). This method can be applied to large numbers of earthquakes over large distances.
Local seismotectonic conditions in the surrounding area of the Molucca Sea are poorly understood. To explain the local seismotectonic setting in this area, we analyze an earthquake sequence that occurred in the area using regional seismic data. Analyzing the waveform inversion of mainshock and aftershock events provides a more comprehensive identification of the source mechanism (Saputra et al. 2021). An analysis of moment tensor inversion is very important to map the variety, complexity of tectonic deformation and also to understand the geodynamic process in a region. In this study, we examine 14 earthquakes with magnitude Mw > 4, the largest of which at Mw 6.8 occurred on July 7, 2019 and was also the largest earthquake in the area during 2019.
Examination of earthquake sequences based on moment tensor in this study is used to determine the seismotectonic trend of the Molucca Sea. Focal mechanism obtained from the moment tensor inversion provides a manifestation of the activity of geological structures (Stein and Wysesion 2003). The analysis of the mainshock-aftershock sequence provides information on the process of relative stress release associated with major and minor faults around the study area. The results of full-waveform inversion produce accuracy regarding the type of source mechanisms and may indicate the character of local fault structures. The moment tensor solution of full-waveform inversion can be effectively computed using ISOLA software, which has been broadly used in other seismotectonic studies (Agurto et al. 2012;Choi et al. 2012;Benetatos et al. 2013;Cambaz and Mutlu 2016;González et al. 2017;Supendi et al. 2018). This software has been developed to automatically calculate centroid moment tensor solutions (Triantafyllis et al. 2016;Vackáař et al. 2017). Thus, we employ the ISOLA software to retrieve the centroid moment tensor solutions. This study may provide important information for seismotectonic pattern in the Molucca Sea area.

Seismic stations and dataset
In this study, we analyze seismic waveforms obtained from July 2019 Molucca Sea earthquakes sequence. The waveforms were recorded by the Indonesian Meteorological, Climatological, and Geophysical Agency (BMKG) seismic network (IA-Network). We selected the earthquake data related to July 2019 Molucca Sea earthquake using the GFZ catalog. The initial dataset used in this study is obtained from 18 broadband IA stations. Each station is equipped with a three-component broadband seismometer (NS, EW, and Z). The closest (GLMI) and farthest (BKSI) of the seismic stations are about 230 km and 1000 km away from the earthquake epicenter, respectively. The distributions of the seismic stations and sequence of earthquakes are shown in Fig. 2.
In this study, we analyze 14 earthquakes that occurred in the Molucca Sea from July 7 to July 18, 2019 (see Table 1). The obtained datasets are in Standard for the  Exchange of Earthquake Data (SEED) format. We then convert these waveforms into SAC format. We check the data manually by plotting the waveforms to ensure that we process the waveforms with a good signal-to-noise ratio. The main event occurred on July 7, 2019, with a magnitude of Mw 6.8 at a depth of 29 km and the epicenter location about 110 km west of Ternate Island. This largest event was followed by 13 smaller earthquakes (magnitudes ranging from Mw 4.6 to Mw 5.8) that took place within a 40 km radius.

Method
To provide valuable information on the deformation process, the determination of the source mechanisms requires full waveform analysis. Our determination of the earthquake mechanism is carried out using the full waveform inversion method executed by ISOLA software, which was developed by Sokos and Zahradnik (2008). Wen et al. (2015) showed that the waveform inversion could be used to obtain precise fault-plane solutions. This enables us to better understand the seismotectonic process around a region. The ISOLA software packages consist of two processing parts, which are ISOLA Fortran code and ISOLA-GUI written in MATLAB. ISOLA Fortran code is used for moment tensor calculation. On the other hand, the ISOLA-GUI interface is used for the data pre-processing and plotting. We use ISOLA in the data processing to retrieve the moment tensor parameters and identify the source mechanism of the earthquake sequences. The moment tensors are computed by minimizing the difference between synthetic and observed seismograms using the least-squares method. Then, the source mechanisms are used to interpret the local seismotectonic setting in the surrounding area of the earthquake.
The synthetic seismogram is determined by Green's functions in the Fortran code of ISOLA, which uses the discrete wavenumber method (Bouchon 1981(Bouchon , 2003Coutant 1989). This method is carried out in the time-domain with the lowest-possible frequency, which reduces the effect of subsurface inaccuracy of the adopted crustal model. The calculation of Green's function using this method is applicable to the regional and local events (Sokos and Zahradnik 2008). The method optimizes the moment tensor solutions through spatiotemporal grid search conducted in two steps to find the best-fitting source position and time. First, we create a model consisting of trial point sources distributed along the vertical line below the epicenter with a depth interval of 1 km and a time interval of 5 s around the origin time with 0.1-s increments. In the second step, we design trial sources distributed along the lateral plane grid at the preferred source depth previously obtained from the first step with a grit interval of 2 to 4 km. Since there is no appropriate velocity model around the study area, thus we use the AK135 earth model (Kennett et al. 1995) for generating synthetics that allows the inversion routine to be more reliable for this study (see Fig. 3). Consequently, this restricts the inversion to low-frequency bands that are less affected by structural heterogeneities as well as by inaccurate velocity model (Larson and Ekström 2001). Furthermore, the waveforms used for this study are also limited for low-frequency waves since we invert waveforms from low magnitude events providing signals higher than the noise level only at low-frequency waves. These two restrictions allow us to use low-frequency bands for the inversion process. The optimum moment tensor solution for the largest earthquake (Mw 6.8), which occurred on July 7, 2019 (event #1), is within the range of frequency band of 0.01 to 0.03 Hz used in our analyses. On the other hand, the optimum moment tensor solution for the smallest earthquakes (Mw 4.6), which occurred on July 8, 2019 (events #4 and #10), is within the range of frequency band of 0.04 to 0.05 Hz. In this study, the optimum result shows the range of frequency band of analysis is narrower for the smaller earthquakes (see Table 2). Fig. 3 The AK135 velocity model which is applied for the calculation of the Green's functions in this study and adopted from Kennett et al. (1995) The double-couple (DC) percentage of the moment tensor solution indicates how good the inversion result is. The selection of the best-fit solution considers the highest correlation with the optimum DC parameter. A previous study shows that the double-couple percentage of less than 0.5 (DC < 0.5) provides significant uncertainty relating to the focal mechanism (Carvalho et al. 2016). The quality of the moment tensor solution is eventually inspected by the variance reduction (VR), which measures the residual between the observed and synthetic waveforms. The VR value is good if it is higher than 0.4 (VR > 0.4) (Sokos and Zahradnik 2008). We also calculate condition number (CN), which shows the resolvability of the moment tensor solution. A smaller CN value indicates better resolvability (for the exact value, see Křížová et al. 2013). Thus, the best solution is selected as that with the highest VR and DC but the smallest CN, which represents the best fit between observed and synthetic waveforms. The best solutions will be interpreted to understand the local seismotectonic setting of the July 2019 Molucca Sea earthquake sequence.

Moment tensor solutions
An example of observed and synthetic waveforms from the earthquake sequence moment tensor analysis is shown in Fig. 4. The best-fitting solution between observed and synthetic waveforms is indicated by the parameters of VR, DC, and CN (see Table 3). We only examine moment tensor solutions with high VR and DC values and low CN values. In this study, most of the results have VR and DC values higher than 0.4 and 0.6, respectively, and CN values below 6. This suggests that our inversion results are well fit for the observed waveforms.
Our fault-plane solutions of the July 2019 Molucca Sea earthquake sequence are shown in Fig. 5. The results of moment tensors show that the epicenter depths can be classified as shallow earthquakes based on the earthquake depth classification of Spence et al. (1989) with a depth range of 10-63 km (see Table 4). We observe that some events show differences in focal depth between our inversion result and the earthquake catalog (Table 4). This is also observed by many previous studies (e.g., Semmane et al. 2005;Bai et al. 2017;González et al. 2017;Jian et al. 2018). We suggest that the differences may be due to lack of station coverage and complex velocity structure in the region. In this study, all tectonic activities occur in the upper lithosphere and the Molucca Sea Plate, according to the cross-section of tectonic settings used by Silver and Moore (1978) and Zhang et al. (2017). McCaffrey et al. (1980) suggested that the tectonic regime in the upper lithosphere at this region has a thrust fault pattern with general azimuth of fault plane in the NNE-SSW direction.
Our results suggest that a transpressional structure pattern with the NE-SW strike direction is the major pattern in this area, such as shown by the #1, #5, #7, #12, #13, and #14 events. The strike-slip structure pattern for the major fault plane is N-S direction, as indicated by events #2, #3, #6, and #8. The compressional structure pattern with the major fault plane of NW-SE azimuth is depicted by events #4, #9, and #10. On the other hand, event #11 has a tensional structure pattern in the form of a normal fault, which may correlate to the strain release structure from the compressional force in this area. Therefore, we suggest that this earthquake cluster maps the existence of a major fault plane dominated by thrust movements. Figure 5a depicts the interpretation of the fault plane solutions of all earthquakes in this study. It reveals the main structure pattern of thrust fault orienting NE-SW. This thrust fault elongates in the middle of the Molucca Sea and may contribute to the uplift of the Mayu Island, which is located in the northern side of the study area. A cross section of this major structure. as shown in Fig. 5b may indicate the presence of a west dipping thrust fault structure delineated under the Molucca Sea. This interpreted thrust fault system is represented by the ten compressional structures obtained from our inversion result. Different orientation of the thrust faults is suggested as the segmentation of the structure in the same fault zone under the same compressional regime. This finding is consistent with the geological interpretation from the studies conducted by Hall and Smyth (2008) and Waltham et al. (2008) as shown in the geological tectonic cross-section across the Molucca Sea region (see Fig. 1b). They suggest that the presence of trust faulting zone around the area may be associated with the Talaud-Mayu ridge above the crest of the Molucca Sea Plate.
In this study, we observe several distinct fault planes during the July 2019 Molucca Sea earthquakes sequences. This condition may be interpreted as an event on July 7, 2019, that is considered as a mainshock. The event may be associated with a transpressional movement that resulted from the stress release process on the major fault. The mainshock may be triggered by the movement of the fault structures in the surrounding area based on the static stress change. Assuming that aftershocks are triggered by a change in static stress, then there is a shift on the mainshock fault that causes shear stress change in the surrounding area. If the shear stress change exceeds the coefficient of static friction, this will trigger a compressional movement. In contrast, if the shear stress is significantly reduced, a tensional (normal fault) movement will be triggered.
The obtained 13-moment tensor solutions represent the compressional force and shear stress on the sedimentary basin between the Sangihe and Halmahera arcs. The thrust fault mechanism within the north-south azimuth is the main structure combined with the strike-slip fault mechanism. On the other hand, event #11 has tensional stress as a normal fault mechanism, which has different principal stress to the other events. This event may be determined as the strain release of compressional and shear-stress mechanisms in this region. We infer that tectonic activity with the compressional regime in this area may control the uplift of the Talaud-Mayu Ridge, which is elongated on the north-south azimuth in the middle of the Molucca Sea caused by the over-thrusting of the sedimentary basin. This result is in good agreement with the previous studies (Silver and Moore 1978;McCaffrey et al. 1980;Hall and Smyth 2008), suggesting that the thrust faulting in the middle Molucca Sea is associated with the development of the Talaud-Mayu Ridge.

Comparison of moment tensor between our result and global data
In this study, the GFZ catalog provides a moment tensor solution of five events (events #1, #5, #6, #7, and #14) based on the calculation by Saul et al. (2011). Meanwhile, the GCMT presents that of seven events (events #1, #5, #6, #7, #12, #13, and #14) (Dziewonski et al. 1981;Ekström et al. 2012). We compute the percentage of double-couple (DC) components of GCMT and GFZ catalogs using Vavryčuk's (2001Vavryčuk's ( , 2015 calculation (see Eq. 1). The percentage of DC data is expressed as:   where M is the total scalar seismic moment, and M DC is the scalar seismic moment of DC component computed as the function of the eigenvalue of a seismic event (see Vavryčuk 2015). We then examine the accuracy of our inversion solutions by comparing them with those published in global catalogs (GCMT and GFZ), as shown by Table 3. We found that the double-couple percentages of the events obtained in this study are generally comparable to those obtained from the global catalogs. However, the double-couple percentages for events #1 and #7 differ from GCMT and GFZ. These differences may explain the differences in the orientation of the source mechanism as indicated by the differences of eigenvalue between the results (Dahm and Krüger 2014;Vavryčuk 2015).
Furthermore, our study provides the completeness of moment tensor solution for lower magnitude portion during the July 2019 Molucca Sea earthquake sequence with reliable moment tensor solution condition. However, it is difficult to precisely determine the source parameter of small earthquakes, such as moment magnitude, due to the lack of data quality. Hence, we calculate local magnitude (ML) to compare with our new moment magnitude (Mw) from the inversion result. ML is commonly used to estimate the size of an earthquake precisely due to less affected by attenuation. In this analysis, we compare our newly obtained moment magnitude (Mw) from moment tensor solution with local magnitude (ML) to clarify the scaling relationship, especially for the lower magnitude portion. We clearly explain this analysis and its result in Additional file 1.

Conclusion
The moment tensor solutions from the Mw 6.8 July 7, 2019, Molucca Sea earthquake and its 13 aftershocks with Mw > 4 occurring along the Mayu Ridge have been computed using regional full-waveform inversion. We present 14 focal mechanisms to give a better comprehension of source mechanisms from the earthquake sequence occurring in the July 2019. Our results are complementing the GCMT and GFZ moment tensor solutions i.e., events #2, #3, #4, #8, #9, #10, and #11. In general, we obtain good and reliable moment tensor solutions in this analysis, as shown by high DC and VR percentages, and small CN values. The obtained focal mechanisms for the July 2019 Molucca Sea earthquake sequence indicate domination of thrust and strike-slip faults mechanism. There is only one normal fault event which possibly associated with strain release in the study area. The earthquake sequence occurs at shallow depth (< 70 km) and lies below the sedimentary basin (1) DC (%) = 1 M (M DC ) × 100% of the Molucca Sea Plate. Therefore, we interpret that the moment tensor solutions obtained from this study indicate the presence of a major thrust fault system combined with the strike-slip mechanism, which may be responsible for the development of the Talaud-Mayu Ridge. In addition, the normal faulting occurring in this area may act as a minor mechanism representing the strain release of the compressional and strike-slip stresses. We may assume that the tectonic setting did not change significantly since the sequence of our analyses occurred within one month. Thus, this study provides updated information about the seismotectonic pattern during the July 2019 Molucca Sea earthquake sequence.
Additional file 1. Relationship between Mw and ML.