On the potential for megathrust earthquakes and tsunamis off the southern coast of West Java and southeast Sumatra, Indonesia

High seismicity rates in and around West Java and Sumatra 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 seismic data sourced from the Agency for Meteorology, Climatology, and Geophysics (BMKG) of Indonesia and the International Seismological Centre (ISC) for the period April 2009 through to July 2020, in order to conduct earthquake hypocenter relocation using a teleseismic double-difference method. Our results reveal a large seismic gap to the south of West Java and southeast Sumatra, which is in agreement with a previous GPS study that finds the region to be a potential future source of megathrust earthquakes. To investigate this further, tsunami modeling was conducted in the region for two scenarios based on the estimated seismicity gaps and the existence of a backthrust fault. We show that the maximum tsunami height could be up to 34 m along the west coast of southernmost Sumatra and along the south coast of Java near the Ujung Kulon Peninsula. This estimate is comparable with the maximum tsunami height predicted by a previous study of southern Java in which earthquake sources were derived from the inversion of GPS data. However, the present study extends the analysis to southeast Sumatra and demonstrates that estimating rupture from seismic gaps can lead to reliable tsunami hazard assessment in the absence of GPS data.


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 earthquake prone Sumatra segment. Historical earthquakes in western Java in the early twentieth century have been studied by Newcomb and McCann (1987), who suggested that this region experienced two major earthquake events with magnitudes in excess of 7.5 in 1903 and 1921. Over the last 2 decades, a further two major earthquakes occurred in this area, namely the Mw 7.8 megathrust event on July 17, 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 intraslab earthquake on September 2, 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); another interpretation is that the Java subduction zone cannot accommodate large megathrust events (e.g., Okal 2012). Indeed, over the last 100 years, the regions along the southern coast of western Java, e.g., Pelabuhan Ratu, Pangandaran, and south of Banten (Fig. 1), have only been subjected to moderately large earthquakes (Mw < 8) and their associated tsunamis. However, recent studies of tsunamigenic deposits along the south coast of Java (e.g., Harris et al. 2019) suggest that large megathrust events do occur in this region and have a return period of ~ 500 years. Consequently, it is important that seismic gaps to the south of western Java are studied in more detail, similar to what has been done in the 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 presence of seismic gaps south of western Java and southeast Sumatra through relocation of earthquake data from April 2009 to July 2020; such gaps are a possible indicator of strain accumulation along a subduction interface that may ultimately be released seismically via a large earthquake. After confirming that large seismic gaps exist, we derive the fault geometry of what we assume is the locked region, based on the identified gap segments and backthrust fault, and investigate its tsunamigenic hazard potential.

Data and method
We use 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 catalog is quite limited due to sparse station coverage; this greatly improved with the establishment of InaTEWS. We combined P-and S-wave arrival-time data from 436 seismic stations at local, regional, and teleseismic distances (Fig. S1a). The initial number of earthquakes 1 3 available in our catalog 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° 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 largely 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 fixed and use 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  Altamimi et al. (2017). The topography and bathymetry data were taken from Batimetri Nasional data (BATNAS) (https:// tanah air. indon esia. go. id/ demna s/#/ batnas). Inset shows the location of the study area (blue rectangle) with respect to southeast Asia. Red lines correspond to the major crustal faults in the region extracted from Irsyam et al. (2020) to mitigate location errors due to insufficient data and large azimuthal gaps in the local seismic network.
We also conducted a jackknife resampling analysis (Tichelaar and Ruff 1989;Waldhauser and Ellsworth 2000) to estimate the uncertainties in the earthquake relocations. We created 100 datasets by randomly removing 10% of the travel-time observations and rerunning the hypocenter relocations for all the data, following Halpaap et al. (2019). The standard deviation of these relocations was then used to compute the estimated uncertainties. The location uncertainties are likely to be underestimated due to uncertainties in the velocity model, which are not considered in this analysis. Nevertheless, the jackknife test is still a useful tool for evaluating the robustness of hypocenter estimates, at least in a relative sense.
We modeled the tsunami height using the TUNAMI modeling code of Imamura (1995), which solves the long wave equation with finite differences. To do so, we projected the Batimetri Nasional data (BATNAS) 6 arc-second interval datasets provided by Geospatial Information Agency of Indonesia (https:// tanah air. indon esia. go. id/ demna s/#/ batnas) onto a rectangular grid of points in latitude (3-12 o S) and longitude (100-110 o E). The model domain, which is projected into Cartesian coordinates, has 1200 and 1080 nodes in the x (longitude) and y (latitude) directions, respectively, with a uniform grid spacing of 936 m. The vertical wall boundary condition, which does not allow for the transfer of momentum flux, is used at the coast. This is a common boundary condition for tsunami modeling that prevents on-shore inundation (e.g., Imamura 1995), which means that our synthetics are relevant to deep ocean tsunamis, not coastal run-ups.

Results of seismicity analysis
For the earthquake relocation component of our study, a total of 20 iterations of the doubledifference algorithm were applied, which yielded 1,082 well-located events (see Fig. 2 and Table 1 in the supplementary material) and 83 poorly located events that were discarded from the final dataset. The average horizontal change in location from before to after DD relocation is 9.4 km (Fig. S1b). The temporal and magnitude distribution of the relocated seismicity is show in Fig. 3, and a comparison of cross sections that show hypocenters from the BMKG catalog and the teletomoDD relocation result is provided in Fig S2. The average horizontal and depth uncertainties from the jackknife test for all 1,082 well-located events are found to be ~ 3.7 and ~ 6.2 km, respectively (Fig. 4), which suggest that the patterns of seismicity we obtain south of West Java and southeast Sumatra are robust.
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 identified as seismic gaps (Fig. 2a). The strongly curved arc of seismicity highlighted in Fig. 2c shows good agreement with the slip-deficit map of Hanifa et al. (2014). This slip deficit was derived by inverting 21 baseline length changes and the absolute rate of vertical displacements at 10 continuously recording GPS stations located in the western part of Java for the period 2008-2010, using the geodetic inversion method of Yabuki and Matsu'ura (1992). The fault plane is assumed to be on the plate interface and extends from the shallow subsurface near the trench to a depth of 60 km. The slip-deficit modeling undertaken by Hanifa et al. (2014) only accounts for the effects of afterslip following the 2006 Mw 7.8 earthquake, and therefore, the estimates produced should be regarded as a lower limit (see Widiyantoro et al. 2020). Since this slip-deficit map was derived using GPS data from 2008 to 2010, our new results suggest that the strain continued to accumulate-assuming no aseismic slip-until at least June 2020, which represents the limit of our data. Further east, an Mw 6.8 earthquake occurred in 2009 (Fig. S3), in the region adjacent to the deeper part of the zone of slip deficit. The 2009 earthquake was an intraslab earthquake dipping westward with a steep dip angle (Gunawan et al. 2019), thus implying that the slip deficit from the megathrust has not been released. The seismicity gaps identified in Fig. 2(a) are still present in a much larger dataset  taken from the ISC-EHB catalog (Fig. S4). To investigate the slip deficit from the seismicity, we calculated the cumulative moment release in the same area that was studied by Hanifa et al. (2014) for all relocated earthquakes between 0 and 50 km depth in the time period from April 2009 to July 2020 as ~ 2.4 × 10 25 dyne cm (Fig. S5); according to Hanifa et al. (2014), the available moment estimated from the expected accumulated slip in the same period is ~ 1.12 × 10 28 dyne cm. This significant difference points to the presence of a sizable slip deficit, under the assumption that aseismic slip is not significant. Although an absence of aseismic slip is difficult to prove, it is noteworthy that Feng et al. (2015) were unable to detect any slow slip events from 10 years of data from the Sumatran GPS array. Nevertheless, since this only pertains to one mode of aseismic slip, we should regard our estimate of moment release deficit as being at the upper limit of what is possible.
Vertical cross sections A, B, and C in Fig. 5 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). Several focal mechanisms in Fig. 5 show thrust A backthrust occurs as a result of layer-parallel shortening in a latestage thrust sequence that has a vergence opposite to the main thrust (Xu et al. 2015). However, Sirait et al. (2020) have shown that such a feature can also be caused by upward fluid migration.
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 filled with unconsolidated sediment that amplifies the shaking (Ridwan et al. 2017;Saygin et al. 2016Saygin et al. , 2017Cipta et al. 2018). For instance, the significant 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:// shake map. bmkg. go. id/). Griffin et al. (2019) showed that Jakarta experienced an MMI VIII when a megathrust scenario occurred on January 5, 1699 (Mw 7.4). This value, however, is the same as when an intraslab earthquake was used in the modeling. If a megathrust earthquake of Mw 8.9 occurs (as in the scenario used in this study), then the ground motion in, for example Jakarta, should be greater compared to the 1699 event, since it has a larger seismic moment. 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. As an aside, we acknowledge that moderate to large earthquakes on the shallow crustal faults will produce more severe shaking than the larger offshore megathrusts. Griffin et al. (2019) showed that a smaller magnitude shallow crustal fault (Mw 5.5 on October 10, 1834) generated an identical MMI VIII in Jakarta to the megathrust event (Mw 7.4 on January 5, 1699).

Fault models for tsunami simulation
We carried out tsunami modeling by identifying gaps in the distribution of seismicity (Fig. 2a), which were then used to define two possible rupture segments, each with a maximum moment magnitude of Mw 8.9 (Fig. 2d), which is consistent with a return period of ~ 400 years (Okal 2012;Harris et al. 2019;Widiyantoro et al. 2020). The western segment has a trench-parallel extent of 325 km, width of 120 km, and a homogenous slip of 24 m, while the eastern segment is 442 km long, width of 109 km, with a homogeneous slip of 20 m (Fig. S6); in both cases, we assume a shear rigidity of 30 Gpa (Davies and Griffin 2020). We also take into account the possible backthrust fault in the south of West Java based on the distribution of seismicity (Fig. 3) and a previous study (Sirait et al. 2020), which may amplify the potential tsunami height along the coast (Heidarzadeh 2011). To include the backthrust fault in the tsunami simulation, it was defined as having a length and width of 312 km and 55 km, respectively, and a homogeneous slip of 16 m (Fig. S6).
We modeled tsunami propagation using two scenarios: Scenario 1 uses a slip of 24 m and 20 m in the western and eastern segments, respectively, without a backthrust segment, while scenario 2 uses the western and eastern segments plus a backthrust segment (slip of 16 m). Thus, the only difference between scenario 1 and scenario 2 is the use of slip along the backthrust segment; these scenarios are shown in Figs. S6 and S7, respectively. We assumed a homogenous slip distribution in the tsunami modeling, which implies that all segments ruptured simultaneously, with no account taken of rupture speed or finite rise time. It has recently been demonstrated (Melgar et al. 2019) that assuming a homogeneous slip distribution may result in lower tsunami potential energies and an underprediction of peak tsunami amplitudes. As such, our results may lie at the lower end of what is possible. In this study, the numerical simulation was run for 6 h to obtain maximum tsunami height along the south coast of Java.

Result of tsunami simulation
The tsunami predictions undertaken in this study show that a maximum wave height could be as large as ~ 34 m on the west coast of the southernmost part of Sumatra and the south coast of Java near the Ujung Kulon Peninsula (Fig. 6). Based on field surveys, the observed tsunami height for the Mw 9.1 2004 Sumatra earthquake is 20-30 m in northern Sumatra (Borrero 2005). This is broadly consistent with our predictions for the Mw 8.9 event in our study, bearing in mind that magnitude is only one factor in determining maximum tsunami height. The average tsunami height along the Sumatran coast and Java coast is 11.8 m and 10.6 m, respectively, a result which incorporates the effects of the backthrust. In comparison, a tsunami height of up to ~ 20 m with an average of 4.5 m was estimated using slip deficit determined from GPS data in Java (Widiyantoro et al. 2020). Thus, our estimated tsunami height is higher than 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 estimated slip used in this study is higher than the average slip of the source inferred from GPS. Homogenous slip implies similar initial wave height for the whole extent of the coast facing the tsunami generation area. This will not allow for discrepancies in the initial height of the tsunami above the rupture that would otherwise contribute to the range of tsunami heights along the coast that is also influenced by the coastal morphology and bathymetric effects during its propagation from deep ocean to shallow shelf regions. The different results are also likely caused by the different megathrust segmentation that was defined; in particular, our study area includes more of Sumatra than Java compared to the Widiyantoro et al. (2020) study. Furthermore, our analysis also used the BATNAS 6 arc-second interval bathymetric dataset, while Widiyantoro et al. (2020) used the General Bathymetric Chart of the Oceans (GEBCO) 15 arc-second interval dataset. However, overall, the pattern of tsunami heights in the south of West Java is relatively similar (Fig. S12).
The existence of the backthrust fault in scenario 2 elevates the negative wave during initial sea surface displacement, thus decreasing the resulting tsunami height along the coast, particularly to the south of Java (Pelabuhanratu), by ~ 3.8 m. This is caused by the dip of the backthrust, which is different to that of the megathrust. The new results for SW Sumatra indicate that a large megathrust event will produce coastal wave heights at least as large as those predicted for West Java, with Enggano Island near the edge of the shelf region at particular risk. It is important, therefore, for further work to be undertaken using complementary datasets and methods to help quantify this risk.  (Hayes et al. 2018) and red lines depict the megathrust zone. The focal mechanisms are taken from the Global Centroid Moment Tensor catalog (https:// www. globa lcmt. org) 1 3

Concluding remarks
Through application of teletomoDD to BMKG and ISC data, a total of 1,082 earthquakes (Mw ≥ 4.0) in West Java and southeast Sumatra for the period April 2009 to July 2020 were relocated. We use a three-dimensional velocity model for the earthquake relocations, which has performed well compared to the one-dimensional velocity model for the initial locations, with estimated depths commensurate with the tectonic setting (subducted slab) in the Fig. 6 a Maximum tsunami height throughout the model region over the duration of the simulation with megathrust and backthrust sources. b Maximum tsunami height along the south coast of West Java, Sumatra, and the Sunda Strait study area. The distribution of hypocenters reveals a seismic gap to the south of the island of Java that is in agreement with results from a previous GPS study and an adjacent gap off the coast of southeast Sumatra. These gaps are clearly delineated by strongly curved bands of seismicity above the subducting Sunda slab. Our relocated events show a near-vertical cluster to the south of West Java, which is likely related to backthrust faulting. The tsunami that is modeled for two plausible megathrust segments and a backthrust that rupture simultaneously shows that tsunami heights can reach up to ~ 34 m on the south coast of southernmost Sumatra and West Java, with average wave heights of around 11 m. In West Java, the results are broadly comparable with an equivalent study that used slip deficits derived from GPS data, although the average wave height is somewhat larger due to our use of a higher resolution bathymetric dataset, greater estimated fault slip, and different megathrust segments. This consistency in results suggests that our estimate from southern Sumatra is likely robust, which has important implications for seismic and tsunami hazard in this region.