The Kalaotoa Fault: A Newly Identified Fault that Generated the $Mw$ 7.3 Flores Sea Earthquake Open Access

On 14 December 2021, an $Mw$ 7.3 earthquake occurred beneath the Flores Sea (Fig. 1). According to Badan Meteorologi, Klimatologi, dan Geofisika (BMKG), otherwise known as the Indonesian Agency for Meteorology, Climatology, and Geophysics, the mainshock occurred at 03:20:23 UTC and was located at 7.59° S, 122.24° E, and 10 km depth (green star in Fig. 1). BMKG used both regional and teleseismic data with the maximum azimuthal gap of 14° to locate this event. From about one minute following the mainshock, a sequence of aftershocks was recorded that continued until at least 20 March 2022 (Fig. S1, available in the supplemental material to this article), the date on which the most recent data used in this study was downloaded. BMKG issued tsunami warnings after the earthquake, although it was later confirmed to be a strike‐slip mechanism. Subsequently, tide gauge data from Flores Island recorded a 7 cm high tsunami. Based on the BMKG shakemap, the earthquake produced ground motion measured at between III and VI on the modified Mercalli intensity scale in Flores Island caused severe damage to more than 346 houses and injured at least 11 people.

Flores Island is located in the Sunda–Banda arc transition at the junction between the eastern limit of the Australian plate subduction margin, defined by the Sunda arc, and the western section of the collision zone between the Australian plate and the Banda arc (Hall, 2002; Audley‐Charles, 2004). In the past 50 yr, there have been three earthquakes that were followed by a tsunami (yellow stars in Fig. 1), that is, $Mw$ 8.3 in 1977, $Mw$ 7.5 in 1992, and $Mw$ 6.5 in 1995 (Supendi et al., 2020). Earthquakes in the region are a consequence of subduction, collision, back‐arc thrust faults, and terrestrial faults (Irsyam et al., 2017). Since 9–10 Ma, the Australian continent has been actively colliding with the western section of the Banda arc (Keep and Haig, 2010) beneath the islands of Timor, Sumba, and Flores (Hall, 1997). Based on seismicity data from 1963 through to 2021 (Fig. S2a), there have not been any shallow earthquakes (depth <100 km) with a magnitude $Mw$ ≥ 7.0 in the region, and very few moderate earthquakes (⁠$Mw$ ≥ 5.0) at depths <100 km (see Figure S2b). Because of the structural complexity of the western Banda arc and lack of data, not all faults have been properly mapped, including the one that ruptured to produce the 14 December 2021 $Mw$ 7.3 earthquake.

To identify the fault responsible for the event, we relocate the mainshock and aftershocks of the earthquake sequence to determine more precise hypocenter locations compared to standard BMKG locations, undertake focal mechanism analysis, and investigate changes in fault stress. The goal is to characterize the fault in terms of its type, its kinematics, and the static stress changes caused by the $Mw$ 7.3 main event.

We use earthquake arrival times from the BMKG catalog between 14 December 2021 and 20 March 2022, which were extracted from continuous data recorded by the broadband seismic stations shown in Figure 1 (red inverted triangles). A total of 1456 events with a magnitude range 1.9–7.3 were identified from 19,151 P‐wave and 5239 S‐wave picks obtained from the mainshock and aftershocks (Fig. S1). BMKG uses LocSAT (Bratt and Nagy, 1991)—a linearized inversion routine that is part of the SeisComP3 program (Hanka et al., 2010) to help construct their earthquake catalog from the arrival time picks they make. As part of this process, hypocenter locations were determined using the IASP91 reference velocity model (Kennett and Engdahl, 1991).

The double‐difference method (Waldhauser and Ellsworth, 2000), implemented via the HypoDD program (Waldhauser, 2001), is used to relocate the mainshock and associated aftershocks based on data recorded by nearby stations (see Fig. 1). This is achieved by minimizing the residuals of the differences between theoretical and observed travel times for nearby earthquake pairs recorded at each station, under the assumption that nearby sources have almost identical ray paths and travel times if the receiver is distant. Therefore, event pairs are identified by locating nearby earthquakes with similar P and S travel times (based on locations, origin times, and arrival time picks from the BMKG catalog). The hypocentral separation maximum was set to 200 km, the maximum number of neighbors per event was set to 100, and the minimum number of links needed to define neighbors was set to 10. We imposed the maximum distance between cluster centroid and station of 200 km. These input parameter choices were made after extensive testing to determine the combination of values that achieved the best results, noting that the final distribution of relocated hypocenters does not change significantly across a reasonable range of input values. We used a 1D velocity model based on CRUST 1.0 (Laske et al., 2013) in the region of interest, combined with AK135 (Kennett et al., 1995), which provides velocities in the uppermost mantle.

Similar to the study by Supendi et al. (2022), we assess location uncertainty through application of a statistical resampling scheme based on the “bootstrap” method (Billings, 1994) for all relocated events. All residuals had Gaussian noise (standard deviation of 0.1 s) added to simulate realistic levels of noise in the test dataset. Using this dataset, we relocated all events and calculated their perturbation with respect to the initial double‐difference location. This entire process was repeated 200 times, thus allowing for the estimation of 95% confidence ellipsoids for each event.

By constraining fault geometry and slip using the aftershock distribution and focal mechanisms computed earlier, we determined the Coulomb static stress change (King et al., 1994) caused by the coseismic slip through application of Coulomb 3.3 (Toda et al., 2011). The subsurface ruptures (length = 98.5 km, width = 22.5 km, and slip = 1.5 m for cluster 1; length = 38.9 km, width = 13.2 km, and slip = 0.8 m for cluster 2; and length = 53.1 km, width = 15.7 km, and slip = 0.9 m for cluster 3) were defined using the empirical equation from Wells and Coppersmith (Wells and Coppersmith, 1994). We calculated static stress changes assuming an elastic half‐space, as implemented in Okada (1992), for the three strike‐slip source faults (average of fault parameters for each cluster are shown in Table 1) at 12.2 km depth, in which we assume a Poisson ratio of 0.25 and a coefficient of friction of 0.4 (Song et al., 2018). Four separate scenarios were tested, each corresponding to a different receiver fault on which future ruptures may occur. These include the central back‐arc thrust near Flores, segments 1 and 2, and the Selayar fault near the southern tip of Sulawesi.

We relocated the mainshock and a total of 1403 aftershocks from the 2021 Flores Sea earthquake from 14 December 2021 to 20 March 2022 (Fig. 2); the remaining 53 events (less than 4% of the total) were relocated above the surface and, therefore, discarded. This reduction in the total number of events is typical, and likely due to data noise and use of an approximate reference model conspiring to introduce errors into the relocation. The BMKG catalog, from which initial locations and arrival time’s picks are sourced, typically uses fixed depths for shallow events, which is likely why this inconsistency was not picked up on earlier. The difference between the initial locations (from BMKG) and relocated aftershocks can be compared in cross section and map view (Fig. S3). Following HypoDD relocation, the magnitudes of the travel‐time residuals are clearly reduced (see Fig. S4), which is suggestive of more robust hypocenters. A significant proportion of events in the BMKG catalog are fixed at 10 km depth, but these are now distributed in depth (Fig. S3) following relative relocation. Relative errors in location based on the bootstrap analysis method are shown in Figure 3. Mean horizontal and vertical location errors are typically less than 3 km (represented by the longest radius of the ellipsoid), and the corresponding maximum mislocations are less than 7 km (Fig. 3 and Fig. S5). However, due to limited station coverage to the east–northeast of the earthquakes we study (Fig. 1), horizontal location uncertainty is larger in the east–west direction compared to the north–south direction, although depth uncertainty is still greater.

The distribution of the relocated aftershocks exhibit an interesting pattern, with at least three clusters (Fig. 2) appearing to be associated with the $Mw$ 7.3 earthquake, which we interpret as a sign of slip along what we call the Kalaotoa fault system. Cluster 1 is caused by slip on the ∼100 km long east–west main fault that also accommodates the $Mw$ 7.3 earthquake, which appears to trigger the ∼40 km long southeast‐oriented segment in the east (cluster 2) and ∼50 km long west‐northwest‐oriented segment in the west (cluster 3). The width of these three segments, based on visual analysis of the aftershock distribution, is ∼20 km. Such faults appear to be quite rare in Indonesia, but can be found in other regions, for example, the $Mw$ 6.1 Ludian earthquake, China (Niu et al., 2019), the 2016 $Mw$ 7.1 Kufeastmamoto earthquake, Japan (Lin and Chiba, 2017), and on the 2019 Ridgecrest earthquake sequence, California (Fialko and Jin, 2021). In addition, there is another cluster to the west of cluster 3 that we interpret to occur on an extension of the Selayar fault (Fig. 4), which is a normal fault (Irsyam et al., 2017). Furthermore, another cluster to the south of cluster 2 may be associated with the central back‐arc thrust fault that was triggered by the $Mw$ 7.3 event (Fig. 4). In this study, it is not easy to differentiate between what could be regarded as an aftershock or a triggered event, which is why these terms are used somewhat interchangeably. However, given that the mainshock occurs on the fault segment that defines cluster 1, it would be reasonable to define events in the other clusters as triggered events, particularly as finite fault predictions from teleseismic data estimate the rupture length of the $Mw$ 7.3 event to be less than 100 km.

The misfit between the observed and synthetic traces following inversion for the mainshock and aftershock focal mechanisms has a value of ∼0.45, indicating a satisfactory fit according to the criterion of Cesca et al. (2010). We produce double‐couple parameters that correspond to a right‐lateral fault, based on the aftershock distribution for the earthquake, with an average strike of N82°W and dip 74° SW (segment 1), an average strike of N42°W and dip 71° SW (segment 2), and an average strike of N76°W and dip 78° SW (segment 3), as shown in Figure 4 and Table 1. Interestingly, both the aftershock distribution and focal mechanisms suggest that segment 3 is offset to the north of segment 1 by ∼10–15 km; this “gap” appears to be accompanied by a decrease in aftershock activity and is worthy of further investigation. Uncertainties associated with the focal mechanisms determined for the mainshock, and two representative aftershocks for segment 2 and segment 3 are shown in Figures S9–S11, and indicate that they are well constrained by the data.

Insight into the possible locations of future earthquakes (Stein and Lisowski, 1983; King et al., 1994) can be gained by comparing the Coulomb stress change and aftershock distribution. The stress changes caused by the $Mw$ 7.3 mainshock (Fig. 5) reveal several interesting patterns. In the case of Figure 5a, it seems clear that there is a stress increase at the endpoints of the source faults, which is unsurprising because that would be, in which we might expect future failure on the same fault surface to occur; Figure 5b is fairly similar to Figure 5a, as one might expect; Figure 5c is interesting, because it shows that there is a stress increase on the back‐arc thrust fault, in which we see aftershocks, which may reflect rupture on that fault. However, there are regions of stress increase on either side of these aftershocks, which may suggest that failure is now more likely in these regions of the back‐arc thrust fault. Similarly, Figure 5d is interesting because there is a clear zone of stress increase at the southern end of the Selayar fault, but not a lot of aftershocks, except on its western segment. This suggests that risk of failure on this fault has increased. A more detailed analysis of the implications of relatively high static stress transfer to the northwest and southeast with respect to the mainshock is urgently needed to investigate possible future fault rupture scenarios.

The complex fault system in this region accommodates the varying orientations of the ongoing plate convergence in the eastern Sunda–Banda arc, which forms a curvilinear feature bounded to the south by the northern Australian continental margin and to the east by the Arafura Sea margin (Norvick, 1979; Hinschberger et al., 2005; see inset map in Fig. 1). Further complexity is added by the presence of the transition from subduction to arc‐continent collision to the south of Flores, which involves a slab tear and underthrusting fore‐arc lithosphere (Zhang and Miller, 2021). It is, therefore, unsurprising that a variety of fault types across a range of orientations are present to accommodate these changes. Identification and understanding of these fault systems is a major challenge, but essential if we are to understand seismic hazard in southeast Indonesia.

The $Mw$ 7.3 Flores Sea earthquake and associated aftershocks were caused by slip along the newly identified dextral strike‐slip Kalaotoa fault system, which appears to be composed of three segments, with lengths of ∼100, ∼50, and ∼40 km. The segment to the west is offset to the north of the central segment by ∼10–15 km. Our Coulomb stress change predictions based on the three source segments were undertaken for a number of realistic scenarios that include strike‐slip, normal, and thrust faults, all of which occur in the study region as a result of its complex tectonic setting. We find that there are a lack of aftershocks in regions of high Coulomb stress change that may eventually yield a large earthquake on adjoining fault systems, including the central back‐arc thrust adjacent to Flores and the Selayar fault near Sulawesi. A more detailed study of the implications of relatively high static stress transfer to these two areas with respect to the main shock is urgently needed to properly assess the region’s seismic hazard potential.

The earthquake dataset exploited in this study is sourced from Badan Meteorologi, Klimatologi, dan Geofisika (BMKG). All figures in the main article and supplemental material were generated using The Generic Mapping Tools (Wessel and Smith, 1998). Topography and bathymetry data were obtained from the Digital Elevation Model Nasional (https://tanahair.indonesia.go.id/demnas/#/demnas) and Batimetri Nasional (https://tanahair.indonesia.go.id/demnas/#/batnas). The earthquake relocation catalog data are available at https://doi.org/10.5281/zenodo.6592435. The supplemental material includes Figures S1–S11. All websites were last accessed in March 2022.

The authors acknowledge that there are no conflicts of interest recorded.

Badan Meteorologi, Klimatologi, dan Geofisika (BMKG) is thanked for providing data from their earthquake archive for use in this study. The University of Cambridge funded this research through a Herchel Smith Research Fellowship awarded to Pepen Supendi Konsorsium Gempabumi dan Tsunami. BMKG is also acknowledged for their support.