Mechanism for seismic supershear dynamic rupture based on in-situ stress: a case study of the Palu earthquake in 2018

Abstract The Palu earthquake, with a magnitude of 7.5, severely devastated the region in Central Sulawesi, Indonesia, but the mechanism of dynamic rupture considering in-situ conditions is not well understood. The in-situ stress field and fault geometry play important roles in the evolution of the supershear rupture. This study investigated the in-situ stress field, critical factors including the inversion of principal stress orientation by the focal mechanism solution, and constraint magnitude using the improved lateral pressure coefficient polygon. Furthermore, the dynamic rupture process is simulated using the finite element method (FEM) considering the strike and in-situ stress magnitudes of the seismogenic fault. The result shows that the principal stress orientation rotates counterclockwise about 15° from north to south, and the northern part of the fault zone is less straight than the fault from Palu Bay to the south, which are two decisive factors for supershear in the north. The damage zones considering the peak acceleration and co-seismic displacement are enlarged in Palu City, which is attributed to the rupture mode. It provides plausible explanations for this supershear event and sheds light on different dynamic rupture processes and seismic hazards that can be predicted by fully understanding the regional in-situ stress field.


Introduction
On 28 September 2018 at 10:02 (UTC), the M W 7.5 earthquake struck the Indonesian island of Sulawesi and triggered a tsunami in Palu Bay. This disaster caused the death of two thousand individuals and destroyed nearly seventy-thousand buildings. After the earthquake, a growing number of studies using synthetic aperture radar (SAR) or geodetic surveys have provided robust evidence of the severe damage caused by this event (Bao et al. 2019;Socquet et al. 2019). Dynamic models of supershear earthquakes show that the rupture velocity of this event significantly exceeded the S wave velocity of the region's crust, and was the dominant cause for the event. The longer the duration of the supershear rupture, the larger is the number of disasters that can occur. The 1999 M W 7.6 event in Izmit (Bouchon et al. 2001), the 2002 M W 7.9 event in Alaska (Walker and Shearer 2009), and the 2010 Mw6.9 Yushu earthquake (Zhu and Yuan 2020) are examples of past supershear earthquakes. These earthquakes are mainly strike-slip, and the Palu earthquake belongs to this model of supershear earthquakes. However, according to the results of space geodesy research, both the significant supershear rupture and the largest surface deformation occurred in the Palu Bay area, located 80 km south of the epicentre, instead of the region at the onset of seismic rupture (Socquet et al. 2019). Thus, it is necessary to study the mechanism of the seismic dynamic rupture process of the Palu earthquake. Figure 1 shows that the Palu-Koro fault zone is a partially straight strike-slip structure located between the Makassar and Sula blocks. Wu et al. (2021) reported new fault structures in the northern Palu Basin near Palu City, marking the previously unmapped Palu-Koro fault, which accurately reflects the actual fault structure and rupture range of the event. Supershear earthquakes that occur on straight faults have ubiquitously been validated by geodetic surveys (Valkaniotis et al. 2018;Bao et al. 2019) and numerical simulation (Liu and Shi 2020). Bao et al. (2019) used teleseismic Figure 1. The geological structure of the Palu region. The epicentre (yellow stars) is located south of Palu Bay. The earthquake induced a new surface fault with co-seismic rupture (red line and dot; Wu et al. 2021). The earthquake ruptured through Palu Bay and continued to propagate on the identified faults (orange line and dot), and many previously mapped faults that did not rupture (black line) belong to the Palu-Koro fault zone. backprojection to reveal that the rupture rapidly reached a steady velocity of approximately 4.1 km/s, exceeding the local S-wave velocity. Fang et al. (2019) used Advanced Land Observing Satellite-2 (ALOS-2) interferometry synthetic aperture radar (InSAR) data along with broadband regional seismograms to demonstrate that the earthquake was a supershear event. Socquet et al. (2019) used geodesy to conclude that the Palu earthquake probably ruptured an extremely linear segment of the 30 km section of the rupture south of Palu city at supershear velocities of 4.3-5.2 km/s. Ulrich et al. (2019) revealed that a physics-based earthquake and coupled tsunami 3 D model can reasonably provide information on the mechanisms and competing hypotheses, but the calculation of the stress magnitude depends on two stress shape ratio assumptions of 0.5 and 0.7, without further constraints. Therefore, a systematic study was conducted, including but not limited to the Palu event, indicating that the shear and normal stresses on seismic faults control the rupture velocity transition on strike-slip faults (Burridge 1973;Andrews 1976;Burridge et al. 1979;Andrews 1985;Dunham 2007). Stress heterogeneities and fault geometry are two important factors affecting the rupture mode.
Some researchers have summarised the regularity of rupture under non-uniform stress distribution on faults (Hu et al. 2021), and others have also proved that supershear rupture could propagate in a stepwise manner on curved faults (Hu et al. 2016;Tadapansawut et al. 2021). Unfortunately, few studies have analysed natural seismic events while simultaneously considering the fault geometry and stress of in-situ conditions. One important reason is the barrier to determining in-situ stresses at the hypocentre depth, which is much shallower than the focal depth; therefore, the initial stress is one of the more difficult parameters to estimate in a dynamic rupture. Generally, empirical estimates and theoretical predictions are utilised to define the stress state in deep rock masses. The former primarily relies on empirical fitting formulae (Hoek and Brown 1980), physical experimental models (Zoback et al. 1987;Starr, 2011), and numerical simulations (Parsons 2006;Matsuki et al. 2009;Zhong et al. 2018).
Regarding the theoretical calculation of in-situ stresses, Zoback et al. (1987) introduced the fault frictional strength theory (i.e. stress polygon theory) and constrained the possible maximum and minimum effective principal stresses with the frictional strength. The results are typically inaccurate (the scale is too large). The results are usually too vast; therefore, a secondary constraint using borehole damage data is necessary. In the Kontinentales Tiefborh program der Bundesrepublik Deutschland (KTB), Brudy et al. (1997) established the relative relationships between two horizontal principal stresses by examining the distribution angle of drilling-induced fractures and rock tensile strengths on the borehole wall at depths of approximately 8 km, as well as estimated the total stress tensors at depths of 7.7 and 8.6 km using the polygon theory. Using the same method, Chang et al. (2010) constrained the deep stresses in two of the four boreholes in the Nankai wedge based on wellbore ruptures and drilling-induced tension cracks in the four boreholes. The images of boreholes several hundred meters below the seafloor captured at different times after completion were compared. Moore et al. (2011) discovered that borehole wall damage estimations from resistivity imaging are broader, and the horizontal principal stress estimations using stress polygon methods are higher. Secondary constraints using borehole damage data are typically required. However, the stress state of the hypocentre is not well understood.
Driven by these results, we first used a method of constraining the in-situ stress using a lateral pressure coefficient polygon and a focal mechanism solution (FMS; Wang et al. 2019). To improve the accuracy of the constraint results, a normal distribution analysis was conducted on the shape ratio of the FMS. Second, the threesigma rule is used to further constrain the deep lateral pressure coefficient in the Palu area, and the deep stress state of the fault plane is calculated according to the occurrence of the fault. Finally, we used the open-source software, Pylith (Aagaard et al. 2013) to establish a 2 D Palu fault model (mainly considering the strike feature) and simulate the earthquake rupture process. The results show that stable supershear rupture occurred approximately 7 s after the onset of rupture when the rupture front reached Palu Bay, by considering both the in-situ stress field and fault structure. Simultaneously, there is an apparent Mach cone phenomenon causing a much greater surface displacement than in other regions, which shows a high degree of similarity to the SAR analysis and geodetic data. The numerical simulation of earthquake rupture, based on the in-situ stress and fault structure, also aids in predicting the level of seismic damage and provides implications for the mitigation of seismic risks.

Method background
In light of Anderson's theory and the Mohr-Coulomb criterion, Zoback et al. (1987) constrained the in-situ stresses in the crust. Based on these theories, we can derive the relationship between the maximum and minimum effective principal stresses versus frictional strength for a fault region without considering the cohesive force, as follows: where r 1 and r 3 are the maximum and minimum effective principal stresses, respectively; S 1 and S 3 are the maximum and minimum principal stresses, respectively; P p is the pore water pressure; and l is the fault frictional coefficient. According to Anderson's fault classification theory, Zoback et al. (1987) introduced the concept of a stress polygon and constrained the possible maximum and minimum effective principal stresses with the frictional strength, but the use of this method is complicated. Therefore, Wang et al. (2019) modified Equations (1) and (2) after considering the corresponding relationships among the principal stresses for the three basic fault types: where k max and k min are the maximum and minimum lateral pressure coefficients for effective stresses, expressed as k max ¼ r H/ r v and k min ¼ r h/ r v , respectively; r H and r h are the maximum and minimum horizontal stresses, respectively; r v is the vertical stress in the strike-slip fault stress state; r H is r 1 ; r h is r 3 ; and the other symbols are the same as defined above. The shape ratio R is an essential parameter in the inversion of FMS and describes the relative level of the three principal stresses. Additionally, it is an important tool for investigating the interactions between FMS and the in-situ stress field at the focus. It is expressed as follows (Gephart and Forsyth 1984): where r 2 is the intermediate principal stress and the other symbols are the same as those defined above. Notably, R obtained by the inversion of multiple FMS satisfies the normal distribution; thus, we can use the three-sigma rule to further constrain the range of R (denoted as R'). According to Equation (3), by introducing the corresponding relationships among the principal stresses for a particular stress state into Equation (3), we can obtain the functional expression between the horizontal principal stresses and R' as the coefficient: The symbols in Equation (4) are the same as those defined above. The derivation process from Equations (2) to (4) is shown in Supplement I.
From Figure 2, the constraint lines of R 0 for different faulting stress states always pass point (1, 1), and when R 0 is an extreme value (0 or 1), it coincides with the boundaries of the polygons for different faulting states. This further verifies that the stress factor can constrain the stresses for a second time based on constraints by a lateral pressure coefficient or a stress polygon. Figure 2 illustrates the constraint principle of lateral pressure coefficients according to this method.

Seismogenic structure and earthquake catalogue
The Mw7.5 Palu earthquake occurred at the intersection of the Australian, Sunda, and Philippine plates with strong tectonic activity. The geological structure of Sulawesi Island in Indonesia is complex. The basin where Palu City is located is a near-south-north valley formed by the long-term left-lateral strike-slip activity of the Palu-Kolo fault, and active faults have developed in the eastern and western margins of the basin (Bellier et al. 2006). The rupture with left-lateral strike-slip is compatible with normal fault characteristics. The M W 7.5 earthquake occurred on 28 September 2018 at 10:02 (UTC); the epicentre was located at 0.178 S, 119.84 E and the hypocentre depth was 13.5 km, as reported by USGS (2020), all FMSs are listed in Supplement II. Figure 3A shows that 82 small FMSs were collected from 1 January 1976 to 28 September 2018 from the Global Centroid Moment Tensor (GCMT) catalogue (Dziewonski et al. 1981;Ekstr€ om et al. 2012), which are all listed in Supplement I. These events are located in the area of 119 E-121 E, 1 N-2 S. For the sake of simplicity, all FMSs are divided into the three regions comprising grids with a size of 2 Â1 , as shown in Figure 3B, which is a graphical visualisation used to display the FMS properties of groups of earthquakes (Frohlich and Apperson 1992). The results show that events in regions I and II mainly belong to the strike-slip type, and events in region III are mainly normal strikes. As described in section 2.1, an FMS similar to the standard Anderson's faulting model is convenient for constraining the deep stress field. Below, we constrain the lateral pressure coefficients and invert the stress orientation for the three subregions using this method.

Inversion of the in-situ stress field in the study area
The constraint of the regional stress field embodies two fundamental parts: the inversion of stress direction and the constraint of stress magnitudes. We used the joint iterative inversion method introduced by Vavrycuk (2014). The inversion results of the three regions, including the orientations of the three principal stresses and distribution of the shape ratio, are shown in Figure 4A. In parallel, the normal distribution curve (red line in the histogram) and the three-sigma range of the stress shape ratio (purple shadow in the histogram) were also obtained by statistical analysis. The type of stress field in the three regions is compatible with the classified results of FMS; notably, the stress fields of regions I and II are approximately strike-slip faulting in Anderson's model, and that of region III is strike-slip normal. Thus, the principal stress orientations of the three regions could be approximated as 279.85 , 295.23 , and 325.75 . Additionally, the ranges of the stress shape factors and confidence intervals up to 99.73% were 0.161-0.379, 0.462-0.654, and 0.384-0.625, respectively.
The magnitudes of the principal stresses can be constrained by combining them with the lateral pressure coefficient polygon. First, we assume that the stress state is limited by Coulomb frictional sliding, and the ratio between the two lateral pressure coefficients are above certain values (boundary of the lateral pressure coefficient polygon) defined by the coefficient of friction (ls) of 0.6 (Byerlee 1978). Subsequently, the shape ratio from the FMS inversion was used to further constrain the lateral pressure coefficient. Figure 4B and 4C show the ultimate constraint results for each region. The constrained ranges of the lateral pressure coefficient, which are shown below, are good guidelines for the in-situ stress magnitudes on seismogenic faults. According to the lateral pressure coefficient calculated by the stress polygon and rock density obtained from CRUST1.0 (Laske et al. 2013), completed stress profiles at a depth of 20 km in the Palu fault, including regions I and II, are shown in Figure 5. Moreover, the in-situ stress magnitude of the two regions at the hypocentre depth can be determined by these stress profiles. Table 1 presents the in-situ stress magnitude and orientation of the hypocentre depth based on the lateral pressure coefficient.

Model configuration and parameter setup
Given that the seismogenic fault was nearly strike-slip in the Palu earthquake, the real 3 D rupture problem can be idealised as a 2 D plane strain with a multi-segment  curved fault. Model geometric features mainly refer to the study of surface ruptures and preexisting geological structures (Wu et al. 2021). Figure 6B shows the integral model and cell grid. We observed that the finite element model (FEM) has a square geometry of 20 Â 100 km to include the entire seismogenic fault, the details of which are listed in Table 2. Moreover, we set two virtual seismic stations near the epicentre (S1) and Palu City (S2) to record the seismic wave.
The mechanical parameters mainly include the initial stress, fault frictional constitutive, and boundary conditions. With a good understanding of the regional stress field, the magnitude of initial shear stress and normal stress on the fault plane can be calculated by the orientation of principal stress and fault occurrence (Jaeger et al. 2007): where l, m, and n are the direction vectors of the fault plane; s and r are the shear stress and normal stress, respectively; and the other symbols are the same as those defined above. Given that the Palu fault is mainly located in regions I and II, where the stress fields are approximately strike-slip faults, we considered that r 1 ¼ r H , r 2 ¼  r v, and r 3 ¼ r h . The density of the overlying strata can calculate the vertical stress. Table 2 presents the strata parameters queried by the CRUST1.0 (Laske et al. 2013). Subsequently, it is easy to determine rH and rh if the lateral pressure coefficient and vertical stress are known. The lateral pressure coefficients were obtained by considering the boundary of the results shown in Figure 4B. Therefore, considering a hypocentre depth of 13.5 km, as reported by the USGS (2020), Figure 6A shows the shear stress and normal stress on multi-segment seismogenic faults at a depth of 13.5 km. Furthermore, the widely used friction laws mainly include the slip-weakening friction law and the rate-state friction law. The slip-weakening friction law was proposed by Ida (1972), which is based on the fault dislocation theory and fracture mechanics and is mainly used to describe the dynamic rupture process of an earthquake (e.g. Day 1982;Harris and Day 1993). The rate-state friction law (Dieterich 1981;Ruina 1983), which includes the fault slip rate or other state quantities, can completely describe the co-seismic and inter-seismic stress accumulation processes. To simulate the dynamic rupture process of the Palu earthquake, we assume that the fault is governed by a form of slip-weakening friction law. This was validated to simulate dynamic rupture using the open-source software Pylith in the community benchmark problem TPV14 (Harris et al., 2018); the parameters used are listed in Table 2. This uniform assumption is necessary for analysing the relationship between the in-situ stress field and rupture mode in the Palu earthquake. The model borders are enclosed by an absorbing boundary that can efficiently remove wave reflections.
where l s is the static friction coefficient, l d is the dynamic friction coefficient, jdj is the fault slip, and d c is the characteristic slip distance. In addition to the above setup, it is necessary to define a nucleation zone on the fault for spontaneous rupture in the hypocentre area (red stars in Figure 6B). Generally, there are two general nucleation approaches for dynamic rupture models, including the time-weakening (TW; Andrews 1985) and overstressed patch (Galiset al., 2015). We used the latter strategy, one of the most widely used methods of nucleation (Bizzarri 2010;Liu et al. 2014), to set the nucleation zone at the hypocentre of the Palu earthquake. The size of the nucleation zone was 2 km, and the  (Laske et al. 2013) initial shear stress of nucleation is set to be slightly above the static fault strength ($0.5% of the strength;$1.5 Mpa).

Simulated fault rupture propagation
To analyse the mechanism of the supershear rupture of the Palu Earthquake in 2018, we built a model configuration and initial parameters according to the structural features and in-situ stress magnitudes of the seismogenic fault. The dynamic rupture process of the Palu earthquake was simulated using the open-source software, Pylith (Aagaard et al. 2013). The results demonstrated the process of dynamic rupture propagation of this supershear event. Figure 7 shows that the spontaneous rupture propagates from the nucleation of the seismogenic fault with low velocities. After 5.26 s, the rupture front travelled to the south of the epicentre (about 30 km), and the Mach cone, characteristic of supershear rupture, appeared as the particle velocities near the rupture front increased. However, the particle velocities of rupture decrease because the stress on the next segment is not prone to supershear. The stable supershear rupture accompanies the Mach cone occurring at Palu Bay and continues to propagate for the next 4 s. The complete rupture propagation from slow to supershear in our simulation is in good agreement with the study results from the SAR and geodetic recordings (Valkaniotis et al. 2018;Socquet et al. 2019). Figures 8 and 9 show the velocities and accelerations recorded by the two seismic stations. Notably, the seismic wave travels to the S1 station approximately 0.3 s after nucleation. There are similar accelerations and velocities in the x and y directions, and the peak velocity and peak acceleration are 0.67 m/s and 17.39 m/s 2 , respectively. In contrast, velocity and acceleration recorded by S2 in the y direction are greater than those in the x direction because the Palu earthquake was mainly a left-lateral strike-slip along the N11 S. As a standard seismic wave with low noise, the velocity increased sharply to a peak value and gradually decreased to approximately zero. Owing to the impact of the Mach wave, the peak velocity of Palu City was up to 3.52 m/s, which is approximately five times that of S1. Similarly, the peak acceleration of S2 was twice that of S1. Additionally, owing to the seismic reflections caused by the continuous perturbation of the nucleation zone near the epicentre, the S1 station recorded the wave signal again after 7.5 s, which remained for approximately 5 s.
Compared to non-curved faults, curved faults cause different physical processes of rupture. For example, different angles affect the propagation and arrest of ruptures . For a multi-curved fault system, the interaction between multicurved segments may also cause a supershear rupture on one segment of the fault under certain conditions (Bhat et al. 2004;Kame et al. 2003;Aochi et al. 2000). Similarly, as shown in Figures 6-9, our simulation results show that the fault structure and in-situ stress background play crucial roles in the evolution of the fault dynamic rupture model. On the one hand, the northern part of the rupture, from the epicentre to its intersection with the Palu Bay coast, is less straight than the rupture from Palu City to the south, which provides a prerequisite for rupture acceleration and stable propagation. Conversely, the stress fields in the two regions, wherein the Palu fault was located in both, belong to strike-slip faulting states, but the principal stress direction has a rotation about 15 between regions I and II. The strikes of  multiple segments for all faults are roughly consistent. The angle between the fault strike and principal stress direction of region II is less than that of region I owing to rotation.
Moreover, the larger the angle enhancement, the larger is the normal stress and the lower is the shear stress. For planar faults with homogeneous stress, the occurrence of supershear ruptures depends on the stress ratio S (Dunham 2007), S ¼ (s p Às 0 )/(s 0 Às r ), where s 0 is the background shear stress, and s p and s r are the peak and residual frictional strengths, respectively. Supershear ruptures occur if S < 1.77 in two dimensions (Andrews 1976;Das and Aki 1977). As we assume that the fault frictional strength is homogeneous based on Byleer's law (Byerlee 1978) and that the stress ratio S is determined by the shear stress and normal stress, S depends only on the in-situ stress and fault occurrence according to Equation (5). For the sake of simplicity, Figure 10 displays the relationship of the angle h (between the fault strike and principal stress azimuth) with the stress ratio S under the in-stress in regions I and II. Figure 10A shows S values varying on the fault along-strike, which indicates that most segments near Palu Bay prefer to undergo supershear than sub-Rayleigh rupture (S > 1.77), and the opposite is true in the region near the nucleation. As shown in Figure 10B, the angle h between the fault strike and principal stress azimuth is all converted within 90 . Considering the in-situ stress in the two regions, the dependence of h and S reveals that the appropriate interval for supershear in region I is significantly larger than that in region II. Moreover, almost all h frequencies of fault segments near Palu Bay tend to undergo supershear rupture, including all in region I and Palu Bay in region II. This could explain why segments of faults in Palu City and Palu Bay are more prone to supershear rupture than others based on the in-situ stress field. Additionally, the high shear stress and extreme inhomogeneity of the stress field between regions III (strike-slip compatible normal) and II (strike-slip) may explain the nucleation location in the northern part of the seismogenic fault.

Discussion
In this study, the lateral pressure coefficient polygon was improved by the stress polygon, and the stress shape ratio was used to constrain the in-situ stress at the hypocentre depth of the 2018 M w 7.5 Palu earthquake. Simultaneously, an equal proportion fault model is constructed to simulate the dynamic rupture process of this event by considering the geometric structure. The results show that the rupture mode of this event is controlled by the stress state of the fault plane, and only a stable super-shear rupture is formed after the rupture front travels to the Palu Bay area, which is consistent with the results of space geodesy (Socquet et al. 2019). According to simulation results, the dynamic mechanism of this event is revealed by a sufficient understanding of the in-situ stress field and fault geometry. The key parameters of earthquake disasters, including the peak acceleration and surface displacement, are discussed below.

Peak acceleration of the co-seismic area
Seismic damage can be enlarged when supershear ruptures engender Mach waves that are directly related to the ground motion. Through the simulation results, we extracted the maximum amplitude of the resultant acceleration (peak acceleration) from all the steps of the dynamic rupture. To study seismic damage, we analysed the peak acceleration in the study area. Figure 11 shows that standardised acceleration was used to compare the relative damage in different regions. The 2 D model can only simulate the acceleration at the depth of the hypocentre rather than the peak ground acceleration (PGA) of shallow ground. We selected four different measuring lines perpendicular to the fault strike to obtain the acceleration profile in the epicentre north of Palu Bay, Palu City, and south of the Palu fault.
Notably, the peak acceleration contour range increases significantly from north to south, which represents the zone of co-seismic damage. In particular, significant expansion of the seismic damage zone with a peak acceleration above 5 m/s 2 , located between lines II and III (where the stable supershear rupture occurs as the simulation result), and consistent with the point of supershear rupture, can enlarge the damage in most cases (Das 2010;Vallee and Dunham 2012;Liu et al. 2014;Xu et al. 2021). In reality, the supershear rupture caused by the heterogeneous stress on the fault segments is the decisive factor affecting the distribution of the damage zone in this event. Additionally, our dynamic rupture simulation could not describe the damage near the shoreline area in detail (for example, Qa Pagerawu and Bali Asih), which is due to the homogeneous material used in our model, and did not account for field effects by heterogeneous material properties, such as sand. This problem should be addressed and optimised in future work.

Co-seismic displacement
Co-seismic displacements are vital for characterising earthquake-induced faults and understanding earthquake dynamics. Bao et al. (2019) reported along-track displacements from the ALOS-2 SAR offsets and bathymetry after the Palu earthquake. As shown in Figure 12A, the along-track displacements are almost parallel to the fault strike and show that the maximum slip is located near the city of Palu and a clear difference in slip magnitude exists between the northern and southern segments. Many subsequent studies have used geodetic data, teleseismic data, and other satellite optical images to show similar results (Valkaniotis et al. 2018;Bao et al. 2019).
Similarly, simulated final displacements in the y direction (approximate alongstrike) after the termination of the rupture were extracted to analyse the co-seismic displacement. As shown in Figure 12B, the maximum displacement in the y-direction was located near Palu City. As a result of the supershear rupture, the average displacement along the fault in the northern part of the Palu fault is significantly higher than that in the southern part. The relative co-seismic displacement based on numerical simulations shows a high degree of similarity with the SAR and geodetic data. The co-seismic displacement on the fault near Palu City is significantly larger than that at the epicentre, and the maximum displacement near Palu City is up to 7 m, which is also similar to the field measurement reported by Valkaniotis et al. (2018). Therefore, the simulation of the Palu earthquake validated that the level of seismic hazards, including the peak acceleration and displacement, could be predicted by dynamic simulation with a full understanding of the regional insitu stress field. In future studies, a 3 D simulation can be considered to fully and quantitatively evaluate the actual surface displacement and peak ground acceleration instead of a relative qualitative evaluation. Simultaneously, it is essential to establish a complete stress profile and fine 3 D fault structure for accurate simulation studies.

Conclusion
1. The lateral pressure coefficient polygon combined with the shape ratio R 0 based on the three-sigma rule, referred to as the stress polygon, is a novel and convenient way to repeatedly constrain the lateral pressure coefficient. 2. The orientation of the regional principal stress in the Pallu fault zone rotates counterclockwise about 15 from north (region II) to south (region I). Moreover, the northern part of fault zone is less straight than the rupture from Palu Bay to the south, which is a joint influence on the stable supershear in the north. 3. The supershear rupture of the Palu earthquake enlarged the zone of peak acceleration and surface displacement to 7 m near Palu City. The severity of damage in Palu City is associated with the significance of the peak acceleration and surface displacement attributed to the supershear.
The simulation results of the Palu earthquake, considering the in-situ stress condition and the fault structure, provide an explanation for initiating the supershear propagation of this event. This further provides the possibility for locating the seismic failure zone before the large earthquake and the prediction of the fault dynamic rupture process, and also provides a reference for disaster prevention and mitigation.