Introduction

On January 14, 2021, a destructive earthquake (Mw 6.2) between Mamuju and Majene, West Sulawesi, Indonesia, occurred at 18:28 UTC (red star in Fig. 1). This event was preceded by a Mw 5.9 foreshock on the same day at 06:35 UTC (Fig. 2). Based on the shakemap from the Agency for Meteorology, Climatology, and Geophysics (BMKG), the earthquake produced strong ground motion with a V–VI MMI scale in West Sulawesi, caused severe damage to more than 279 houses and resulted in 84 fatalities (National Disaster Management Authority/NDMA). Peak Ground Acceleration (PGA) values of the Mw 6.2 mainshock from two of the nearest stations in Mamuju and Majene are 95.9 and 92.8 Gals, respectively, equivalent to VI on the MMI scale (Additional file 1: Figure S1). Over the past five decades, there have been two deadly earthquakes around Mamuju and Majene (see red stars in Fig. 1). The M 7.0 (23 February 1969) earthquake and tsunami resulted in 61 casualties (Prasetya et al. 2001), and the M 7.0 (8 January 1984) earthquake resulted in two casualties (Lander et al. 2003). The purpose of this study is to analyze the 2021 (Mw 6.2) Mamuju–Majene earthquake and its foreshocks and aftershocks through hypocenter relocation, focal mechanism solutions, and analysis of fault rupture process and stress changes.

Fig. 1
figure 1

Map of the Sulawesi region, with the study area denoted by a blue box. Geological features are taken from Watkinson (2011). The three red stars depict the 15 January 2021 (Mw 6.2) mainshock from BMKG and the historical destructive earthquakes around Mamuju are taken from the USGS earthquake catalog. The historical seismicity (M ≥ 6.0) and Centroid Moment Tensor (CMT) solutions for the entire Sulawesi region in the time period 1977 to 2020 are taken from the global CMT catalog, except for the 1969 (M 7.0) event taken from Fitch (1972) and the 2021 (M 6.2) event from the BMKG catalog. The blue inverted triangles are the BMKG seismic stations used in earthquake relocation process; red traces represent major crustal faults in the region extracted from Irsyam et al. (2017). The inset map shows the Indonesian region, together with Indo-Australian and Pacific plate motions relative to the Eurasian plate

Full size image
Fig. 2
figure 2

Magnitude vs time for foreshocks, mainshock, and aftershocks of Mamuju–Majene earthquake. The detection threshold of our analysis is M 2.3

Full size image

The island of Sulawesi is located within a complex tectonic region at the confluence of the Eurasian, Indo-Australian and Philippine plates (Hall et al. 2011; Watkinson 2011). The geology of the area records the complex history of subduction, extension, obduction, and collision of continental fragments during the Mesozoic–Cenozoic Era (e.g., Katili 1978; Hamilton 1979; Hall 2002). The island forms a unique “k” shape with four elongated arms (known as the north, east, southeast, and south arms) consisting of distinct lithological assemblages (Fig. 1). The E–W oriented north arm consists of a Neogene island arc rock assemblage with a mix of oceanic continental crustal fragments (Elburg et al. 2002; Leeuwen et al. 2007). This arm connects to the main island by a narrow mountainous range composed of metamorphic and ultramafic basement (Sukamto 1973). The east and southeast arms are characterized by highly deformed ophiolite imbricated onto Mesozoic and Cenozoic sediment units (e.g., Hamilton 1979; Simandjuntak 1986; Parkinson 1991). Western Sulawesi is characterized by Miocene–Pliocene magmatic arc rocks, dominated by granitic intrusions through a basement of Cenozoic volcaniclastic and Mesozoic metamorphic rocks (Elburg et al. 2002; Hamilton 1979; Hall and Wilson 2000). All of these rock units have been accreted onto the southeastern margin of Sundaland, the continental core of southeast Asia, since the Middle Cretaceous (Hamilton 1979; Parkinson et al. 1998; Hall et al. 2009), with rapid uplift occurring during the Pliocene and continuing to present day (Bellier et al. 2006). All of these rock units are dissected by many deep-seated major fault systems, such as the Palu-Koro, Matano, Gorontalo, Lawanopo, and Walanea faults (Fig. 1).

In the north, the present day tectonic configuration is dominated by the North Sulawesi subduction zone with a convergence rate of 42–50 mm/year (Socquet et al. 2006). In the center, the structural configuration is mainly controlled by the Palu-Koro fault. This left-lateral strike–slip fault has a slip-rate of 42 mm/year (Socquet et al. 2006) and was the source of the 2018 (Mw 7.5) destructive earthquake and tsunami (Gusman et al. 2019; Supendi et al. 2020). The western region of Sulawesi is characterized by the east-dipping Majene–Mamuju fold and thrust zone. The fold and thrust propagates to the west, where it is bounded by the Makassar Strait Thrust (MST) fault, which is composed of four segments (Irsyam et al. 2017). This fold and thrust zone accommodates most of the westerly accretion of Sulawesi towards the eastern limit of the Sundaland continent, and is the main source of seismicity in this region (Irsyam et al. 2017). Our study area is located in the vicinity of the Makassar Strait Thrust zone (Fig. 1).

Data and method

For this research, we use data from the BMKG earthquake catalog between 14 and 20 January 2021, which is based on recordings from BMKG permanent broadband stations (see the blue inverted triangles in Fig. 1). A total of 34 events of magnitude 2.3–6.2 were recorded by our stations during this period, including foreshocks, the mainshock, and aftershocks (Fig. 2) with 1316 and 187 identified P- and S-wave arrival times, respectively. The initial hypocenters from the BMKG catalog were determined by using the linearized inversion scheme LocSAT (Bratt and Nagy 1991), which is part of the SeisComP3 program (Helmholtz-Centre Potsdam-GFZ German Research Centre For Geosciences and GEMPA GmbH 2008), in the presence of the IASP91 reference velocity model (Kennett and Engdahl 1991). For the relocation of foreshock, mainshock, and aftershock hypocenters, we use the HypoDD program (Waldhauser 2001) to perform the double-difference method (Waldhauser and Ellsworth 2000). HypoDD minimizes residuals between the observed and calculated travel-time differences via an iterative procedure, with the locations and partial derivatives updated after each iteration. We first applied a coupled velocity-hypocenter inversion using the program Velest (Kissling et al. 1994; Kissling 1995) to update a 1-D P-wave velocity model (Vp) derived from CRUST 1.0 (Laske et al. 2013) by exploiting data recorded by 14 of the BMKG seismic stations from 34 earthquakes. Due to the limited number of S-wave arrival times, we did not try and invert for a 1-D S-wave velocity model (Vs). Relative location uncertainties are estimated using a bootstrap method (Efron 1982; Billings 1994) for all relocated events. The synthetic data were calculated by adding Gaussian noise with a standard deviation of 0.1 s to all residuals. We relocated all events using the synthetic dataset and determined the perturbation from the original double-difference location. The process was repeated 1000 times, which allowed error ellipsoids at the 95% confidence level to be estimated for each event. In order to further assess the influence of the dataset on the final event locations, we also randomly perturbed initial locations by adding Gaussian noise with a standard deviation of 10 km prior to relocation.

BMKG’s focal mechanism solutions were determined by moment tensor inversion using the SCMTV program that is included in SeisComP3. This program applies the stable moment tensor inversion algorithm of Minson and Dreger (2008). The Green’s functions were calculated using GEMINI (Friederich and Dalkolmo 1995) based on the IASP91 velocity model (Kennett and Engdahl 1991).

We obtain a rupture model of the mainshock by performing a finite fault inversion. The inversion is constrained by P-wave teleseismic displacement data with good azimuthal coverage (Additional file 1: Figure S2a). The method is an inversion in the wavelet domain using a simulated annealing approach to search for the best estimates of the source model (Ji et al. 2002; Shao et al. 2011; Hao et al. 2013). The fault plane is constructed as a single rectangular east-dipping plane with a strike of 334° and a dip of 32°, taken from BMKG’s focal mechanism solution (Additional file 1: Figure S2b). We divided the fault plane into 77 2 km × 2 km sub-faults. The size of this plane is first estimated by scaling according to Wells and Coppersmith (1994); however, we adjust the size of the plane during the inversion. In this model, the rupture was initiated from the relocated hypocenter (at a depth of 20 km). Five unknowns are constrained for each sub-fault, i.e., the slip amplitude, rake angle, rupture initiation time, and two parameters describing the asymmetric source time function (Ji et al. 2003). A comparison between the observed and synthetic seismograms is available in Additional file 1: Figure S3. The finite fault inversion depends on many factors. Because we use only teleseismic displacement seismograms in this study (distance 30°–90°), which have more reliable temporal resolution compared to of spatial and moment resolution, we can obtain the general picture of rupture propagation and a reliable slip amplitude. They can be represented by the source time functions we provide. Based on the finite fault model, we determined the Coulomb stress change of the co-seismic slip using Coulomb (Toda et al. 2011). We calculated stress changes for the optimal thrust fault (east-dipping plane with a strike of 334° and a dip of 32°) at ~ 20 km depth, assuming a friction coefficient of 0.4 and a Poisson ratio of 0.25.

Results and discussion

We relocated seven foreshocks, a mainshock, and 24 aftershocks of the 2021 Mamuju–Majene earthquake from 14 to 20 January 2021. We compared the relocated aftershocks with the initial locations (from the BMKG catalog) in map view and the vertical cross-section (Additional file 1: Figure S4). The events that had previously been held fixed at 10 km depth from the BMKG catalog have now been relocated using a double-difference method. Relative location errors are shown in Additional file 1: Figure S5. Based on the bootstrap analysis method described in “Data and method” section, average horizontal and vertical mislocations are generally less than 0.8 km, and the corresponding maximum mislocations are less than 1.5 km (Additional file 1: Table S1, Figure S6). The test involving a random perturbation of the initial event locations prior to application of HypoDD reveals an average shift in relocation of 3.5 km, which suggests that the data more strongly influences the final locations than the choice of initial locations. Additional file 1: Figure S7 illustrates the results from this test, which shows that the pattern of relocations is not significantly different.

The sequence associated with the Mw 6.2 Mamuju–Majene earthquake started with a foreshock of Mw 5.9 at 06:35 UTC; this was followed by six foreshocks with magnitudes less than M 5.0 in 12 h before the mainshock. Note that the detection threshold of the array is around M 2.0, so it is likely that there were many undetectable small foreshocks and aftershocks. Based on BMKG data analysis, both the Mw 5.9 foreshock and the Mw 6.2 mainshock occur on thrust faults (Fig. 3). We use the distribution of the relocated earthquakes and their focal mechanism solutions to assess candidates for the causative fault of this earthquake sequence (Fig. 3a). Figure 3b shows the profile of the fault rupture and our interpretation of the fault geometry based on the focal mechanism solution. The up-dip projection of the earthquake distribution in Fig. 3b suggests that the fault structure is located at MST Mamuju; the mainshock fault plane and aftershock hypocenter pattern appears to converge to MST Mamuju at the surface. The vertical cross-section shows most of the aftershocks are distributed across a depth range of 9–12 km, and within close proximity of the foreshock locations. There is only one aftershock recorded near the mainshock hypocenter. Based on the depth distribution of the foreshock and the mainshock–aftershock sequence, we interpret that they occurred on two different fault planes. It is likely that the foreshock modified the stress state of a nearby fault segment, causing it to rupture. The relative locations and behavior of these two fault planes are depicted in Fig. 3b.

Fig. 3
figure 3

a Map view of relocated events and focal mechanism solutions for a Mw 5.9 foreshock and the Mw 6.2 mainshock from the BMKG catalog; b cross-section A–B in West–East direction showing relocated events, focal mechanisms, and an interpretation of how multiple faults cause the earthquake sequence. The solid and dashed red lines suggest that the Mw 6.2 mainshock fault plane rupture is probably associated with MST Mamuju, while the solid blue line highlights the fault plane rupture based on the Mw 5.9 foreshock focal mechanism and aftershock hypocenter pattern

Full size image

The development of fold and thrust belts in a highly compressive tectonic regime such as in Sulawesi is commonly known (Bergman et al. 1996; Puspita 2005; Morley et al. 2011; Wu et al. 2020). This fold and thrust belt may form a series of imbricated synthetic thrust faults (roof) seated on larger decollement thrust fault (floor thrust) (Yan et al. 2016). The 3D seismic images from Brackenridge et al. (2020) reveal the presence of a fold and thrust belt that cuts Pliocene-aged rock units off the west coast of Mamuju. This fold and thrust belt, is known as the Majene Fold and Thrust Belt. The Mw 5.9 foreshock and its associated foreshock events likely ruptured one of the roof thrusts. Meanwhile, the Mw 6.2 mainshock likely occurred on the floor thrust. The down-dip section of these two fault planes might be connected to one another via a decollement (hidden fault) of unknown depth.

Unlike other destructive crustal earthquakes of similar scale (~ M6) in Indonesia, i.e., the 2006 (Mw 6.4) Yogyakarta earthquake (Walter et al. 2008; Budiman et al. 2019), the 2016 (Mw 6.6) Poso earthquake (Supendi et al. 2018; Wang et al. 2019), the 2016 (Mw 6.5) Pidie Jaya earthquake (Muzli et al. 2018), all of which had many aftershocks, the 2021 (Mw 6.2) Mamuju–Majene earthquake had fewer aftershocks.

The limited aftershocks caused by this earthquake are likely related to the kinematics of its source, which involves a relatively high moment release over a small source area and with a short duration (Fig. 4). We find that the primary moment release occurred within 4 s after the rupture initiated, with the peak slip of up to 180 cm at a depth of about 17–21 km. The asperity (with slip > 50 cm) has an average slip of about 106 cm, across an effective length of only ~ 8 km. The rupture appears to have propagated monotonically updip. The static stress drop is estimated from the Kanamori and Anderson (1975) formula for circular rupture, by assuming a rigidity of 42.9 GPa as derived from the CRUST 1.0 model (Laske et al. 2013), with an average slip of 106 cm, and effective length of 8 km, which produces a stress drop at the asperity of about 78 bars, much higher than the average stress drop of ~ 40 bars typically reported at active deformation zones (Ye et al. 2016). Based on source studies of aftershocks of the 2010 (Mw 8.8) Maule earthquake, Sen et al. (2015) suggested that depth-dependent variations in rigidity, rather than frictional conditional stability at the plate interface, is responsible for variations in source durations at subduction zones. This may also be the case in our study, even though it is not at a plate margin. Similarly, Sallarès and Ranero (2019) suggest that shallower events (usually near the trench) have a longer duration, with lower radiation frequencies. By contrast, deeper events like the 2021 Mamuju earthquake (depth ~ 20 km) tend to be of shorter duration and exhibit higher radiation frequencies.

Fig. 4
figure 4

a Moment rate function and b slip distribution, of the Mw 6.2 Mamuju–Majene mainshock. The gray arrows indicate the slip (rake) direction while the numbers (e.g., 3) indicate the rupture propagation time in seconds. The color represents the co-seismic slip in cm

Full size image

Seismic stations near the earthquake epicenter recorded high-frequency ground motion (Additional file 1: Figure S8). The primary moment release duration of 4 s for the Majene earthquake is short and is certainly much faster than the 13 s of a typical rupture of a comparable magnitude 6.2 estimated from an empirical relationship (Ekström and Engdahl 1989). As observed in other studies, the production of aftershocks is largely controlled by the rupture’s aspect ratio and stress drop, and to a lesser degree, by the down-dip width (Dascher-Cousineau et al. 2020). A recent study from Miller (2020) also shows that a lack of aftershocks can be due to the absence of high-pressure fluid sources at depth.

The Coulomb stress change relative to an Mw 6.2 mainshock (Fig. 5) shows areas of increasing stress (depicted in red), towards the northwest of the mainshock, while the areas that have a stress drop (depicted in blue) are towards the southwest relative to the mainshock. The Coulomb stress change may explain the aftershock distributions and the approximate distribution of future earthquakes (Stein and Lisowski 1983; King et al. 1994). It is likely that the mainshock energy release immediately reduced stress at the rupture location; this stress then transferred to the northwest and caused aftershocks, since this area exhibits a high Coulomb stress change (close to 1 Bar).

Fig. 5
figure 5

Modeled Coulomb stress change caused by the Mw 6.2 mainshock of the Mamuju–Majene earthquake. We calculated stress changes for an optimal thrust fault at ~ 20 km depth, with a coefficient of friction of 0.4 and Poisson’s ratio of 0.25. The blue and red colors depict negative and positive Coulomb stress changes, respectively. The white circles denote foreshocks, the mainshock, and aftershocks

Full size image

Based on the PGA map at bedrock (with Vs between 750 and 1500 m/s) for 2% probability of exceedance in 50 years (Irsyam et al. 2020), West Sulawesi has a modeled PGA value of less than 1 g. However, many houses and buildings were severely damaged, and even some government buildings collapsed in Mamuju and Majene. The geological map from Hall and Wilson (2000) shows that the damaged areas are located in quaternary sediments, which likely explains the increased ground motion which peaks at around 96 Gal. Destructive earthquakes have struck in this area in the past (Fig. 1 and Additional file 1: Figure S9). Therefore, spatial and regional planning should use an updated earthquake hazard map that accounts for this earthquake sequence and the ground motion it produced.

Concluding remarks

The Mamuju–Majene earthquakes were likely caused by movement on faults within the Majene Fold and Thrust belt, in which the foreshock and mainshock occurred on two different fault planes. The unusually small number of aftershocks that occurred after this earthquake is likely related to the geometry of the fault rupture, and also to the short duration of the rupture. The mainshock increased stresses in the northwest area, as indicated by the high Coulomb stress change, which coincides with the distribution of aftershocks in this area.