Hypocenter relocations and tsunami simulation for the 15 November 2014 Northern Molucca Sea earthquake in Indonesia

A reverse fault earthquake (Mw 7.1) occurred in the Northern Molucca Sea, Indonesia, on 15 November 2014 at 2:31:40 UTC. The earthquake produced small tsunami waves that are recorded at Jailolo (9 cm), Tobelo (1 cm), and Menado (3 cm) tide gauges. The Indonesian Agency for Climatology, Meteorology, and Geophysics (BMKG) issued a timely (5 minutes after the earthquake) tsunami warning for the event. We used the teleseismic doubledifference seismic tomography method (teletomoDD) to relocate the hypocenters of the mainshock and the aftershocks. The relocated hypocenter of the mainshock for the 2014 Northern Molucca Sea earthquake is located at 1.923N, 126.539E, and depth of 48.87 km. In general, the relocated aftershock hypocenters are shallower than those from the BMKG catalog. The relocated hypocenters are distributed within a depth range of 6 to 64 km. The aftershock area from the relocated hypocenters is 80 km long and 55 km wide. The estimated seismic moment from the Global CMT solution (GCMT) was 4.75 1019 Nm. We simulated the tsunami from fault model of each GCMT nodal plane to find a fault model that can best explain the observed tsunami heights at Jailolo, Tobelo, and Menado tide gauges. The best single fault model for this event is dipping to the west, has fault length, width, and slip amount of 47 km, 25 km, and 1.16 m, respectively. The K value calculated using the observed and simulated tsunami heights for this best model is 1.026, suggests a very good fit to tsunami observations.


INTRODUCTION
The Molucca Sea collision zone lies in the complex junction between large plates of Sunda, Australian, Pacific and Philippine Sea plates (Bird, 2003).Both Sangihe volcanic arc on the west and Halmahera arc on the east are active and convex toward the Molucca Sea (Silver and Moore, 1978).The Molucca Sea collision zone that is surrounded by Sulawesi, Sangihe, Molucca, and Mindanao Islands has produced tsunamis from tsunamigenic earthquakes and volcanic eruptions as described in previous studies (Iida et al., 1967;Berninghausen, 1969).The relocated hypocenters distribution show reversed V-shape Molucca sea slab beneath the area subduct to the west and east direction (Yulianto et al., 2017;Nugraha et al., 2015;Shiddiqi et al., 2014).The seismicity rate is high in this region and also some anomalous of earthquake swarm activities occured in West Halmahera, North Molucca region 1 i C 2017, Himpunan Ahli Geofisika Indonesia which may be related to deep magmatism (Nugraha et al., 2017).
On 15 November 2014 at 2:31:40 UTC, a magnitude Mw 7.1 reverse fault earthquake with a north-south plane orientation and a small oblique component.The Global Centroid Moment Tensor (GCMT) solution (http://www.globalcmt.org) for the earthquake gives an eastward dipping nodal plane with strike/dip/rake of 197/58/77 and a westward dipping nodal plane strike/dip/rake of 42/35/110.The aftershock area of the earthquake was approximately 80 km long and 55 km wide.The steep thrusting caused sea floor uplift that generated small tsunamis along the coastlines in the surrounding islands.A previous study of Shiddiqi et al. (2016) concluded that the mainshock was the reverse type with strike/dip of 1920/550, while the dominant mechanism of aftershocks were consistent with the mainshock.The seismic moment estimated by the Global CMT solution is 4.75 1019 Nm.A tide gauge in a nearby coastal town of Jailolo recorded a tsunami with maximum height of 9 cm.The Indonesian Agency for Meterological, Climatology and Geophysics (BMKG -Badan Meteorologi, Klimatologi dan Geofisika) issued a tsunami warning for the region.The initial tsunami warning was issued at 2:36:01 UTC or approximately 5 minutes after the earthquakes origin time, and the tsunami warning was terminated four hours later.
In this study, we summarize the tsunami early warnings given by the BMKG for the event.We use teleseismic doubledifference seismic tomography method (teletomoDD) to relocate the mainshock and aftershock hypocenters from the BMKG catalog.Single fault models based on the two nodal planes of GCMT solution are evaluated using tsunami simulations.Then we estimate a single fault model for the earthquake that can best explain the observed tsunami data.

TSUNAMI EARLY WARNING FOR THE 2014 MOLUCCA SEA EARTHQUAKE
Tsunami early warning system in Indonesia (InaTEWS) is run by the BMKG, and the agency has the authority to issue tsunami warnings based on earthquake parameters and tsunami simulations in the InaTEWS (Setiyono et al., 2017).For each tsunami event, every location (tsunami warning segment) in Indonesia is given with one of the following warning levels: Warning, Advisory, and No Threat.Each location would also be given an estimated time of arrival (ETA) of the tsunami.However, the current system does not make available the estimation of tsunami height at any location to the public through the InaTEWS website (https://inatews.bmkg.go.id/new/).After an initial tsunami warning bulletin, BMKG would issue updates for each event.
For the 2014 Molucca Sea earthquake and tsunami event, BMKG issued tsunami early warning 5 minutes after the earthquake origin time or at 2:36:01 UTC.The initial tsunami warning bulletin was updated for four times at 2:44:27, 4:02:28, 4:25:01, 5:15:15 UTC.For the initial tsunami warning bulletin, only warning levels are given.The estimated time of arrival only available after the second tsunami warning bulletin that was issued 13 minutes after the earthquake origin time.Table 1.shows the tsunami warning levels and ETAs for locations which are available from the second tsunami warning bulletin provided by BMKG to the public.The 3rd, 4th, 5th tsunami bulletins were issued with observed tsunami heights as they became available.The observed tsunami heights at Jailolo (West Halmahera regency, North Maluku province), Tobelo (North Halmahera regency, North Maluku province), Manado (North Sulawesi province) were 9, 1, and 3 cm, respectively.Summary of the observed tsunami heights are available in Table 2.

Double-Difference earthquake relocation
The Double-Difference hypocenter relocation method (DD) minimized residual between observed and calculated travel time for each pair of earthquakes at each station.(Waldhauser and Ellsworth, 2000.The DD method assumes that if two earthquakes that are separated by relatively small distance compare to the distance to the station, then the ray path of those two earthquakes would be very similar.Therefore, the travel time difference of those two earthquakes recorded at a station would be a function of distance between the hypocenters.We used teleseismic doubledifference seismic tomography method (teletomoDD), an extended version of double-difference tomography (Zhang and Thurber 2003), which allows us to use stations at teleseismic disctances (Pesicek et al. 2008;Pesicek et al. 2010).Travel-time calculation is performed using a 3D ray tracing conducted through nested regionalglobal velocity models and applied the method to relocate the seismicity for in the SumatraAndaman region (Pesicek et al. 2010).
In this study, we use the teletomoDD method to relocate the mainshock and aftershocks of the 2014 Northern Molucca Sea earthquake.Because this method can use teleseismic data to improve the relocation accuracy, the method is ideal for locations with sparse local seismic network such as what we have around the Northern Molucca Sea.As described in (Pesicek et al., 2010), to determine raypaths and calculate travel times in a spherical earth, we utilized the pseudobending (PB) method of Um and Thurber (1987) and extended to a spherical Earth with discontinuities by Koketsu and Sekine (1998).The regional model of 3D velocity structure for Indonesia (Widiyantoro and van der Hilst, 1998) and the global one of ak135 (Kennet et al., 1998) are used.The seismic data recorded at stations around Indonesia and surrounding regions, and the preliminary earthquake hypocenters data provided by BMKG were used as initial locations.We only relocated events that have minimum eight arrival times and maximum azimuthal gap of 210o to ensure the quality of relocated hypocenters.

Tsunami numerical model
The linear shallow water equations can be solved by finite difference method using the staggered grid scheme in the spherical coordinate system to simulate a tsunami propagation (Gusman et al., 2010).The initial sea surface displacement is calculated by equations in Okada (1985) using earthquake source parameters of location, strike, dip, rake, depth, fault dimension, and slip amount.The strike, dip, rake, and depth of the fault are available from the GCMT solution.While the fault dimension and the slip amount will be estimated in this study.We use the GEBCO bathymetry dataset (Weatherall et al., 2015) with grid size of 30 arc-sec for the tsunami simulation.

Single fault model estimation using tsunami data
We conduct multiple runs of tsunami forward modeling using single fault models with different slip amounts, fault sizes, and fault orientations.The purpose of these simulations is to find a single fault model that produce tsunami heights (peak to trough) most similar to those observed at a tide gauge in Jailolo (9 cm), Tobelo (1 cm), and Manado  (3 cm).We evaluate the two fault orientations with eastward dipping and westward dipping planes provided from the GCMT solution.We also evaluate scaling relations of seismic moment or moment magnitude to fault dimension from Blaser et al. (2010) and Murotani et al. (2013).
A fault area can be estimated from seismic moment using a scaling relation of Murotani et al. (2013), and then the respective fault area is used to estimate the slip amount by assuming a rigidity of 4 1010 Nm-2.Fault length and width are calculated from the fault area by assuming length = 2 width.Blaser et al. (2010) provides scaling relations to calculate fault length and fault width from moment magnitude.The relocated mainshock hypocenters is used as the center of the single fault model.We use a trial and error approach using the observed tsunami height (Gusman et al., 2017a;2017b) to obtain the final estimates of fault orientation, dimension and slip amount.

Relocated hypocenter
A total of 148 aftershocks with magnitude 4 were relocated.The relocated hypocenter of the mainshock for the 2014 Northern Molucca Sea earthquake is at 1.9014N, 126.5147E, and depth of 48.87 km.It was shifted from a deeper location of 1.95E, 126.490E, and depth of 60 km based on the BMKG catalog.In general, the relocated aftershock hypocenters are shallower than those from the BMKG catalog.There is no hypocenter from the BMKG catalog that has a depth shallower than 10 km.Whereas the relocated hypocenters are located as shallow as 4 km.Since there is no station within 100 km from the hypocenters, the depth of shallow events are not well constraint.The aftershock area from the relocated hypocenters is very similar to that from the BMKG catalog, which is 80 km long and 55 km wide.From the aftershock distribution, it is still not so obvious which one from the two GCMT nodal planes is the fault plane.
This event was not occurred on the Halmahera Thrust, a west-dipping thrust located at the west of Halmahera arc.Gunawan et al. (2016) conducted a coseismic slip study of the 2014 Molucca sea event using GPS data and suggested that this event was occurred on a west-dipping splay-fault about 50 km west of the Halmahera thrust.

Single fault model for the 2014 Northern Molucca Sea earthquake
We used a scaling relation of Murotani et al. (2013) to calculate the fault area from the seismic moment of the earthquake.The fault length of 60 km and width of 30 km are calculated from the respective area and a simple relation of length (L) = 2 width (W).Slip amount of 0.68 m is obtained from the seismic moment of the earthquake.The scaling relations of Blaser et al. (2010) give smaller fault dimension (47 km 25 km) compared to that from scaling relation of Murotani (2013).The calculated slip amount for this smaller fault dimension is 1.16 m.Hence, we have 4 fault models which are the westward dipping fault model Murotani (WFM-M), westward dipping fault model Blaser (WFM-B), The vertical uplift calculated from the westward dipping fault plane is slightly larger than the one from the eastward dipping fault model.This is because the dip angle for the eastward dipping fault model ( 58) is slightly steeper than the dip angle for the westward dipping fault model ( 35).As a result, the average simulated tsunami from the westward dipping fault model is slightly larger than the one from the one from the eastward dipping fault model.
To evaluate the accuracy of tsunami prediction we use the K value of Aida (1978), which is calculated using the simulated and observed tsunami heights.Table 3. shows that the westward fault model Blaser (WFM-B) gives the closest K value (1.026) to 1, compared with other fault models.The tsunami heights of the first cycle of simulated tsunami waveforms at Jailolo, Tobelo, and Manado from the eastward dipping fault model (WFM-B) are 10.67, 1.20, and 1.95 cm, respectively.The K value can be used to adjust the initial slip amount for each fault model to make the simulated tsunami heights closer the observations.The adjusted slip amounts using the K values (Table 3.) for WFM-M, WFM-B, EFM-M, and EFM-B are 0.82, 1.19, 0.82, and 1.20 m, respectively.
Tsunami simulations using both nodal planes and scaling relations of Murotani et al. (2013) or Blaser et al. (2010) can predict the tsunami rather well with K value smaller than 1.3.Our result suggests that the technique to predict the tsunami using a single fault model is appropriate for a magnitude 7 class earthquakes such as the 2014 Northern Molucca Sea earthquake (Mw 7.1).The application of the technique can be more challenging for larger or more complex earthquake events, such as the 2004 Sumatra-Andaman and the 2010 Tohoku earthquakes.

CONCLUSION
Tsunami early warning was issued by BMKG to the public in 5 minutes after the earthquake occurred.This is a timely tsunami warning and would give enough time for coastal communities to evacuate during the tsunami event.The first tsunami warning bulletin only give warning level for each of warning segment.The updated tsunami warning bulletin include warning level and estimated time of arrival.The time for evacuation until the tsunami arrive was very short because according to BMKG estimation, the tsunami will arrive along the Halmahera segment (the closest segment) at 2:38:48 UTC or in 7 minutes after the earthquake.But the second tsunami warning bulletin was issued a little late at 2:44:27 UTC or after the tsunami reached the shore of Halmahera segment.This estimation is for the time when water level beginning to rise, however, the tsunami peak amplitude will arrive latter.For the case of Jailolo, the tsunami simulation shows that the tsunami initial rise and peak amplitude were arrived in 18 and 25 minutes after the earthquake, respectively.Improvements for InaTEWS might be needed to make the initial tsunami warning bulletin includes the estimated time of arrival (ETA) and issued in 5 minutes after the earthquake.
We have relocated the mainshock and 148 aftershocks using teletomoDD method.The relocated hypocenters distributed along 80 km north-northeast plane.There is no hypocenter from the BMKG catalog that has a depth shallower than 10 km.Whereas the relocated hypocenters are located as shallow as 4 km.Since there is no station within 100 km from the hypocenters, the depth of shallow events is not well constrained.While the preferred fault plane is not well distinguished from the aftershocks distribution, previous study suggested that the earthquake occurred on a west-dipping fault plane.
We simulated the tsunami heights at Jailolo, Tobelo and Manado using four fault models.Each fault orientation available from the GCMT solution is used to make two fault models using scaling relations of Murotani et al. (2013) and Blaser et al. (2010).The scaling relations from the two publications give different fault dimensions and slip amounts.Then the simulated tsunami heights are compared with the observations to fine a best fault model.The results show that the fault model with westward dipping direction (strike/dip/rake of 197/58/77), fault length of 47 km, fault width of 25 km, and slip amount of 1.16 m, is the best fault model.

Fig. 1 .
Fig. 1.Map of Northern Molucca Sea showing focal mechanism of the 2014 earthquake (orange beach ball) and other tsunamigenic earthquakes (blue beach balls) with moment magnitudes lager or equal than 7.0.Red star and gray circles represent the relocated mainshock and aftershocks, respectively.

Fig. 2 .
Fig. 2. Aftershock distribution of the 2014 Northern Molucca Sea earthquake from the BMKG catalog a) and the relocated one using the teletomoDD method b).Red triangles are active volcanoes, two major active faults in this region, i.e.Sangihe Thrust (ST) and Halmahera Thrust (HT) are depicted with black lines.

Fig. 3 .
Fig. 3. a) Distribution of relocated hypocenter along with cross-section lines (dashed lines), b) Cross-section A , c) Cross-section B along with GCMT focal mechanism solution.Interpreted fault plane and top of subducted Molucca Sea Plate are depicted by red and black dahsed lines, respectively.

Fig. 4 .
Fig. 4. Vertical displacements from a) the westward dipping fault model estimated by scaling relation of Murotani et al., (2013) (WFM-M), b) the westward dipping fault model estimated by scaling relation of Blaser et al., (2013) (WFM-B), c) the eastward dipping fault model estimated by scaling relation of Murotani et al., (2013) (EFM-M), d) the eastward dipping fault model estimated by scaling relation of Blaser et al., (2013) (EFM-B).The contour interval for seafloor vertical displacement is 0.02 m.Black rectangular is the fault plane, thick black line represents the shallowest edge of the fault, red star represents the relocated earthquake epicenter, gray circles represent the aftershock epicenters.

Fig. 5 .
Fig. 5. Maximum tsunami amplitude distribution of the best fault model (WFM-B).Triangles represent coastal points at which tsunami waveforms are simulated.

Table II .
Observed tsunami heights at tide gauge stations for the 2014 Northern Molucca Sea earthquake.

Table III .
Simulated tsunami heights (in cm) for the four single fault models.The corresponding accuracy of each fault model is represented by K and.