Potential megathrust earthquakes and tsunamis off the southern coast of West Java, Indonesia

High seismicity rates in and around West Java occur as a result of the Indo-Australian plate converging with and subducting beneath the Sunda plate. Large megathrust events associated with this process likely pose a major earthquake and tsunami hazard to the surrounding community, but further effort is required to help understand both the likelihood and frequency of such events. With this in mind, we exploit catalog and phase data sourced from the Agency for Meteorology, Climatology, and Geophysics (BMKG) of Indonesia and the International Seismological Centre (ISC) for the period of April 2009 through July 2020, in order to conduct earthquake hypocenter relocation using a teleseismic double-difference method. Our results reveal a seismic gap to the south of West Java, which is in agreement with a previous GPS study that nds the same region to be a potential future source of megathrust earthquakes. Tsunami modeling was conducted in the region based on two scenarios using sources from the Indonesian National Center for Earthquake Studies, and show that the maximum tsunami height may be up to ~ 8 m. This estimate is approximately half the maximum tsunami height predicted by a previous study in which earthquake sources were derived from GPS data inversion.


Introduction
West Java is located at an active plate boundary between oblique subduction of the Australian plate beneath Sumatra and orthogonal subduction along Java (DeMets et al. 2010). Despite both regions being located along the same active subduction margin, historical data suggest that the portion of western Java in the Sunda arc is relatively aseismic compared to the highly seismically active Sumatra segment. Historical earthquakes in western Java in the early 19th century have been studied by Newcomb and McCann (1987), who suggested that this particular region experienced two major earthquake events with magnitudes in excess of 7.5, i.e. in 1903 and 1921. Over the last two decades, there were also two major earthquakes that occurred in this area, i.e. the Mw 7.8 thrust fault slow slip event on 17 July 2006, which generated a devastating tsunami in Pangandaran (Fujii and Satake 2006;Mori et al. 2007;Hanifa et al. 2014;Gunawan et al. 2016); and the Mw 6.8 reverse fault earthquake on 2 September 2009, which struck south of western Java (Suardi et al. 2014;Gunawan et al. 2019).
Great earthquakes close to the Java trench are typically interplate faulting events along the slab interface between the Australia and Sunda plates; these earthquakes generally have high tsunamigenic potential due to their shallow depths (Jones et al. 2014). The absence of recent great earthquakes may indicate that even more powerful tsunamigenic events along the south coast of western Java are a potential threat (Widiyantoro et al. 2020). The regions along the southern coast of western Java, e.g., Pelabuhan Ratu, Pangandaran, and south of Banten ( Fig. 1) are prone to large earthquakes and their associated tsunamis. However, seismic gaps to the south of western Java are less studied compared to southwest of Sumatra , which is known to generate large megathrust events (Natawidjaja et al. 2004;Chlieh et al. 2008;Konca et al. 2008).
In this study, we investigate the potential for megathrust earthquakes and consequent tsunamis south of western Java through relocation of earthquake data from April 2009 to July 2020. We also conduct tsunami simulations using megathrust sources provided by the Indonesian National Center for Earthquake Studies (PuSGeN; Irsyam et al. 2020) and compare our results with the interplate coupling sources from GPS data inversions undertaken in previous studies (Hanifa et al. 2014, Widiyantoro et al. 2020. Two different rupture scenarios are tested, with a focus on maximum wave height along the south coast of West Java. Based on tsunamigenic deposit studies by Harris and Major (2017) and Harris et al. (2019), the predicted recurrence interval of large megathrust events is around 500 years.

Data And Method a. Seismicity
We used new data taken from the Indonesian Tsunami Early Warning System (InaTEWS) catalog reported by the Indonesian Agency for Meteorology, Climatology, and Geophysics (BMKG) combined with International Seismological Centre (ISC) catalog data for the period April 2009 through July 2020. Prior to 2009, BMKGs arrival time catalogue is quite limited due to sparse station coverage; this greatly improved with the establishment of InaTEWS. We combined P-and S-waves arrival-time data from 436 seismic stations at local, regional, and teleseismic stations (Fig. S1a). The initial number of earthquakes available in our catalogue is 1,302. However, prior to relocation, we implement selection criteria to eliminate poorly constrained events. The selection criteria are: (1) each event must produce at least 10 P-and S-phase arrival times, and (2) the azimuthal gaps must be less than 210 o for local/regional stations in the Indonesian network. A total of 1,165 events passed the selection process and were subsequently relocated. The epicentral shifts of relocated events are in general perpendicular (north-south direction) to the Java Trench (Fig. S1b), which may be caused by the stations being mainly distributed to the north on Java.
We used a teleseismic double-difference (DD) relocation algorithm (teletomoDD) that is an extension of the DD tomography method (Zhang and Thurber 2003) to teleseismic distances (Pesicek et al. 2010(Pesicek et al. , 2014. However, in this study we hold the initial seismic velocity model xed and used the relocation capabilities only, following (Pesicek et al. 2010). Travel times were calculated using a 3D regional seismic velocity model of the Indonesian region with a grid size of 1°×1° van der Hilst 1996, 1997;Widiyantoro et al. 2011), and the global 1D model ak135 (Kennett et al. 1995) for regions outside Indonesia, using pseudo-bending ray tracing (Um and Thurber 1987) adapted for use in a spherical coordinate system by Koketsu and Sekine (1998). This ray tracing method allows the use of local, regional, and teleseismic P-and S-wave arrival-time data, thereby helping to mitigate location errors due to insu cient data and large azimuthal gaps in the local seismic network. A total of 20 iterations of the double difference algorithm were applied, which yielded 1,082 well located events (see Fig. 2) and 83 poorly located events that were discarded from the nal dataset. The average horizontal change in location from before to after DD relocation is 9.38 km (Fig. S1b). Comparison of cross sections of hypocenters from the BMKG catalog and the teletomoDD relocation result are show in Fig S2. b. Tsunami We carried out tsunami modeling with two different megathrust scenarios from PuSGeN (Fig. S3): (1) West-Central Java segments (Mw 8.7), and (2) Combination of Sunda Strait, West-Central Java, and East Java segments (Mw 9.1). The second scenario is based on the fact that the 2004 Sumatra and the 2011 Tōhoku giant earthquake and tsunami events were generated by more than one megathrust segment rupturing simultaneously (Maercklin et al. 2012;Lay 2015; Bilek and Lay 2018).
We modeled the tsunami height using the TUNAMI modelling code of Imamura (1996) and Imamura et al. (1996). To do so, we projected the General Bathymetric Chart of the Oceans (GEBCO) 15 arc-second interval datasets (https://www.gebco.net) onto a rectangular grid of points in latitude (5 o to 14 o S) and longitude (103 o to 110 o E). The model domain, which is projected into Cartesian coordinates, has 3120 and 2160 nodes in the x (longitude) and y (latitude) directions, respectively, with a uniform grid spacing of 463 m. The vertical wall boundary condition, which does not allow for the transfer of momentum ux, is used at the coast. This is a common boundary condition for tsunami modelling that prevents on-shore inundation (e.g. Imamura, 1996). We model tsunami propagation using two different source scenarios with various slip distributions by solving the long wave equation of Imamura (1996) with a nite difference method. The rst and second scenarios are shown in Fig. S6. Scenario 1 is based on a large event with an estimated magnitude of Mw 8.7, whereas Scenario 2 combines the Sunda Strait, West-Central Java, and East Java segments to produce an even larger event of Mw 9.1. The various slip distribution sources were modeled using an elastic deformation model (Okada 1985) for each segment. It was assumed that all segments ruptured simultaneously and no account was taken of rupture speed or nite rise time (Fig. 4). The model was run for 5 hours to obtain maximum tsunami height along the south coast of Java.

Results And Discussion
Our earthquake relocation results for Mw ≥ 4.0 illuminate regions located between the coast of West Java and the Java Trench that lack seismicity and hence can be identi ed as seismic gaps (Fig. 2a). The strongly curved arc of seismicity highlighted in Fig. 2c shows an excellent agreement with the slip de cit map of Hanifa et al. (2014). Since this slip de cit map was derived using GPS data from 2008-2010, our new result indicate that the strain continuously accumulates, at least until June 2020, in this particular region. Further east, the Mw 6.8 earthquake occurred in 2009 (Fig. S4). The seismicity gaps identi ed in Fig. 2 (a) are still present in a much larger dataset (1963 to 2008) taken from the ISC-EHB catalog (Fig.  S5). Fig. 3 clearly show that events occur along the megathrust. However, we also detected a near-vertical cluster of events that may be related to a backthrust, as has been found in other locations, e.g. in Sumatra (Pesicek et al. 2010) and in Central Java (Ramdhan et al. 2019). A backthrust occurs as a result of layer-parallel shortening in a late stage of thrust sequences in which it is displaced in an opposite direction to that of the main thrust direction (Xu et al. 2015). However, Sirait et al. (2020) has shown that such a feature can also be caused by upward uid migration.

Vertical cross-sections A and B in
Ground motion generated by a great earthquake off Java may cause severe damage to megacities such as Jakarta (Nguyen et al. 2015), which is located on a deep basin lled with unconsolidated sediment that ampli es the shaking (Ridwan et al. 2017;Saygin et al. 2017). For instance, the signi cant earthquakes of magnitude Mw 7.8 (17 July 2006), Mw 6.8 (02 September 2009), Mw 6.5 (15 December 2017), Mw 6.1 (23 January 2018), and Mw 6.9 (02 August 2019) (see the yellow stars in Fig. 2b) shook Jakarta quite severely (MMI up to V; http://shakemap.bmkg.go.id/). Hence, information on seismic gaps that point to the possibility of future megathrust earthquakes along the Java trench is of particular importance for mitigation of seismic and tsunami hazards in the region.
The resulting tsunami heights along the south coast of western Java from Scenario 1 and 2 show that a maximum tsunami height of up to ~ 8 m is likely to occur in the region with an average of 1.8 m for Scenario 1 and 3.3 m for Scenario 2. In comparison, a tsunami height of up to ~ 20 m with an average of 4.5 m was estimated using slip de cit determined from GPS data (Widiyantoro et al. 2020). Thus, our estimated tsunami height is less than half that of the previous study that assumes the same earthquake magnitude. According to the pioneering work of Aki (1966), the seismic moment is equal to the fault area multiplied by the rigidity and slip of the earthquake source. Assuming the same magnitude and rigidity, the slip of the PuSGeN source is smaller than the source inferred from GPS because the fault area of the PuSGeN source is much larger than that of the GPS source. This explains the difference in maximum estimated wave heights.

Concluding Remarks
We conducted hypocenter relocation by applying the teletomoDD method to BMKG and ISC catalogs. A total of 1,082 relocated earthquakes (Mw ≥ 4.0) in West Java for the period April 2009 through July 2020 reveal a seismic gap to the south of the island of Java that is in agreement with results from a previous GPS study. This gap is clearly de ned by a strongly curved band of seismicity to the south of Banten. Our relocated events show a near-vertical cluster to the south of West Java, which is likely related to backthrust faulting. The tsunami modeled for two plausible megathrust segments that rupture shows that tsunami heights can reach ~ 8 m on the south coast of West Java.