Brittle failure of a rock under triaxial stress is governed by the Coulomb failure criterion (Sibson, 1990): τ = C + μ σ ′ = C + μ ( σ − P ) , (1) where τ is the shear stress acting on the plane, C the rock cohesive strength, μ the rock internal friction, σ′ the effective normal stress, and P the fluid pressure in rock. Shear failure occurs when the shear stress (τ) on the failure plane exceeds the failure resistance, i.e., rock cohesive strength (C) and rock internal friction (μ) multiplied by the effective normal stress (σ′) (Oppenheimer et al., 1988; Sibson, 1990). Thus the Coulomb failure stress function is written as Eq. 2a (Oppenheimer et al., 1988): | τ | ≥ μ ( σ − P ) + C , (2) CFS =   | τ | − μ ( σ − P ) − C , (2a) CFS =   | τ | + μ ( σ ″ + P ) − C . (2b)

Failure will occur when CFS ≥ 0. Eq. 2a has a positive sign for compression. This is changed to negative for compression in Eq. 2b (σ″ = −σ). Coulomb failure stress function basically calculates the changes in stress. This is due to the lack of original stress state knowledge (Harris, 2000). Hence, Eq. 2b becomes Eq. 3. Cohesion vanished as it is assumed to be constant over time (Harris, 2000). Δ CFS =   Δ | τ | + μ ( Δ σ + Δ P ) . (3)

Positive ΔCFS is associated with the subsequent event, while negative ΔCFS is associated with stress shadow effect (King et al., 1994). This can be clearly seen from the ΔCFS basic equations, Eqs. 2, 2a. In order for an earthquake (failure) to occur, shear stress acting on a plane has to be greater than the failure resistance; hence, CFS needs to be positive (vice versa for negative ΔCFS).

When the fault is in a reactivation condition, its cohesive or cementation strength is rather low (C ≈ 0) (Sibson, 1990). Under this condition, the Coulomb failure criterion (with the Mohr circle) is plotted in Figure 2. Vavryčuk (2014) formulates fault instability in the above condition (C ≈ 0). If the linear line touches the Mohr circle, i.e., the red dot in Figure 2, failure on a plane will occur. The red dot marks the most unstable and most susceptible fault to failure, which is called the principle fault (Vavryčuk, 2014). Vavrycuk scaled the most unstable fault to 1 and the most stable fault to 0 and defined fault instability using the following equation: I = τ − μ ( σ − σ 1 ) τ C − μ ( σ C − σ 1 ) , (4) where τ _c and σ _c are the shear and effective normal stress along the principle fault (red dot in Figure 2), while τ and σ are the shear and effective normal stress along the observed fault (black dot in Figure 2). σ ₁, σ ₂, and σ ₃ (in Eq. 4; Figure 2) are the effective regional stresses (Vavryčuk, 2014). Fault instability (I) has a range between 0 and 1, with 1 signifying the most unstable fault plane orientation in a certain regional stress field.

Vavryčuk (2014) modified Eq. 4 so that it can be evaluated from friction coefficient (μ), shape ratio (R), and directional cosines (n). Directional cosines (n) define the fault plane inclination from the regional stress axes. The modified equation is as follows: I = τ − μ ( σ − 1 ) μ + 1 + μ 2 , (5) where σ and τ are σ = n 1 2 + ( 1 − 2 R ) n 2 2 − n 3 2 , (6) τ = n 1 2 + ( 1 − 2 R ) 2 n 2 2 + n 3 2 − σ 2 . (7)

Aside from μ, R, and n, a regional stress (maximum horizontal stress) orientation is also important in fault instability calculations, specifically in calculating n. The slip direction (rake) of the observed fault plane is assumed to be parallel to the maximum horizontal stress (S _H) (Kinoshita et al., 2019). The S _H orientation is incorporated in n, as shown in Eq.8 (n for strike-slip regime S _H > S _v > S _h): n = ( n 1 , n 2 , n 3 ) = ( − sin ⁡ δ ⁡ sin ⁡ α , cos ⁡ δ , sin ⁡ δ ⁡ cos ⁡ α ) T , (8) where δ is the observed fault plane’s dip and α is the observed fault plane’s strike subtracted by S _H. We applied Eq. 5 in this study to analyze the instability of faults related to the Pidie Jaya earthquake.

Regional stress information is needed to calculate fault instability. Thus, appropriate regional stress data should be calculated or inferred carefully. The stress inversion of focal mechanism data could be applied to obtain regional stress data, typically the orientation (Vavryčuk, 2014). In our study area (the North Sumatra Basin, see Figure 3A), it was not possible to infer the in situ stress using focal mechanism; this was mainly due to the small number of available earthquake focal mechanisms within a wide time span. Hence, in this study, we inferred the stress orientation and magnitude from existing studies.

The present-day maximum horizontal stress (S _H) orientation is generally assumed to be aligned with the plate convergence vector. For instance, Kinoshita et al. (2019) used this assumption in their study of slip tendency at the forearc area of the Nankai Zone. In Sumatra, the relative plate convergence vector defined from Global Positioning System (GPS) data is N14°E (Sieh and Natawidjaja, 2000). However, studies show that, on a local scale, Sumatra’s SH does not strictly align with N14°E. From the focal mechanism inversion, Sahara et al. (2018) found a variation of the SH orientation at −1° to 2.2°N of the GSF, i.e., N12°E ± 12°, N32°E ± 10°, and N10°E ± 10° from −1° to 2.2°N, respectively. Mount and Suppe (1992) observed an arc-normal SH orientation at the back-arc basin. They found that borehole breakouts in central and southern Sumatra oil field basins show consistent elongation of the maximum horizontal stress at N39°E ± 3.9° and N50°E ± 4.1°, respectively. Though not strictly parallel with the plate convergence vector (N14°E), Mount and Suppe (1992) suggested that strong coupling between the subducting Indo-Australian Plate and overriding Sunda Plate is the underlying force in these basins. Tingay et al. (2010) suggested the orientation in these areas is perpendicular to the adjacent subduction zone, which is slightly oriented NNE-SSW in northern Sumatra and primarily oriented NE-SW in the central to the southern part.

The Pidie Jaya earthquake occurred in the back-arc basin of the northern part of Sumatra Island. Due to lack of local SH information about this area, we inferred the S_H orientation from the nearest basin, that is, from the central Sumatra Basin, to be N39°E ± 3.9° (Mount and Suppe, 1992). In this calculation, we considered the uncertainty range (±3.9°) to cover the possibility of SH which may be less than N39°E in our study area. We inferred the stress magnitude from the in situ stress measurement in the Suban field, South Sumatra, i.e., SH 340 bar/km, Sv 240 bar/km, and Sh 180 bar/km (Hennings et al., 2012).

For the first objective in this study, we aim to conduct a fault instability analysis of the Pidie Jaya earthquake to locate the possible causative fault. We took two faults related to the possible causative fault into consideration: the Samalanga-Sipopok and the Panteraja Faults. The strike of the Samalanga-Sipopok Fault was approximately ∼180° (Barber et al., 2005). The strike and dip of the Panteraja Fault were inferred from Pijay-Net aftershock distribution done by Muzli et al. (2018). The uncertainty was approximated from the width of the aftershock distribution, as shown in Figure 4. We found that the strike is 43.6° ± 4.7° and the dip is 61.2° ± 1.6°. Due to the lack of knowledge about the fault dip for the Samalanga-Sipopok Fault, we calculated fault instability of all possible dips (0°–90°) for both faults.

For the fault’s friction coefficient (µ), we estimated the maximum friction coefficient from the Mohr diagram of the stress magnitude data (SH 340 bar/km, Sv 240 bar/km, Sh 180 bar/km). We found that the failure line is tangent to this Mohr circle at µ of 0.32. Though this value might represent the principle fault's friction coefficient at the place where this data is taken, we still used it as a constraint for estimating the maximum value that we should consider in this study. In addition to the maximum friction coefficient, we also used several smaller friction coefficient values: 0.30 and 0.28.

In the second objective, we aimed to compare fault instability and ΔCFS results to the Pidie Jaya aftershocks. We calculated the fault instability of the focal mechanisms of the Pidie Jaya aftershocks. For the ΔCFS, we calculated the stress change imparted by the Pidie Jaya mainshock to the aftershock focal mechanisms. We used Coulomb3.3 software (Lin and Stein, 2004; Toda et al., 2005; Toda et al., 2011) to calculate the ΔCFS. We assessed both the nodal planes of the focal mechanism in the fault instability as well as the ΔCFS calculation. The Pidie Jaya mainshock slip model and focal mechanisms we used (Figure 5) refer to the results of Muzli et al. (2018) and Rusli (2017).

The focal mechanisms of the Pidie Jaya aftershocks might provide additional focal mechanism data for this area. Revisiting the possibility of inverting stress orientation from focal mechanism data, we plotted the P/T axes of the focal mechanism of the Pidie Jaya aftershocks in Figure 6, using STRESSINVERSE code (Vavryčuk, 2014). The P/T axes of the focal mechanisms could show the direction of maximum and minimum principal stress direction acting on the focal mechanism. From Figure 6, we can see that most of the aftershocks create a convergent cluster of principal stress axes; however, some points deviate, which is shown by red open circles in the lower half-circle and blue positive signs in the left half-circle. This deviation might show slight variation in the direction of local stress acting on the focal mechanisms. The azimuth of (average) maximum principal stress or sigma 1 is oriented approximately NNW, as shown by the green dot in Figure 6. The orientation of maximum principal stress could represent the SH direction. However, in this case, inverting stress orientation from aftershock data might represent the local stress orientation as the result of the mainshock stress perturbation. Thus, we did not invert the regional stress orientation from the focal mechanism of the Pidie Jaya aftershocks.

Vavryčuk (2011) suggests that the major earthquake could impart large stress perturbation which might rotate the orientation of the regional stress. Therefore, it is necessary to analyze the possible stress rotation in the study region. Stress rotation can be observed indirectly from ΔCFS distribution. When the earthquake's stress drop is much larger than the regional deviatoric stress, the optimum slip planes near the fault might be rotated (King et al., 1994).

For that purpose, we conducted a ΔCFS analysis of the major events preceding the Pidie Jaya earthquake. The nearest stress perturbation source in our study area comes from GSF activity in northern Sumatra. Therefore, we conducted static ΔCFS preanalysis using these events. GSF focal mechanisms data during the period of June 1976–December 2016 were compiled from the Global Centroid Moment Tensor Catalog (Figure 3B). We estimated the ΔCFS and found that ΔCFS stress change is small in the target study area (less than 1 bar), as shown in Figure 3C. Thus, the SH rotation due to GSF events is likely insignificant in our study area.

Another possible stress perturbation source arises from subduction zone activity, especially the 2004 Sumatra-Andaman M 9.2 earthquake. Hardebeck (2012) analyzed the stress rotation due to the Sumatra-Andaman earthquake in its rupture zone area, excluding inland Sumatra. Hardebeck found that the earthquake generates stress rotation in the rupture zone which rotated back after several months in the southern rupture zone area (≈2°–5°N). In addition, Rafie et al. (2019) analyzed stress rotation in inland Sumatra. The moderate rotation was found in the northern part of GSF (≈3°–6°N), yet Rafie et al. noted that this result should be cautiously interpreted because of limited available input data (earthquake data). Our study area (4.5°–5.5°N) is in the Northern Sumatra Basin and is off the GSF; therefore, the impact of the megathrust event might be even lower. Moreover, though the area might be exposed to stress rotation, there is a possibility of return rotation to occur.

We calculated the fault instability of the Panteraja and Samalanga-Sipopok Faults using approximated strikes of 43.6° and 180° for each fault. The maximum horizontal stress orientation we used is N39° ± 3.9°E and the magnitudes are SH 340 bar/km, Sv 240 bar/km, and Sh 180 bar/km. We varied the friction coefficient into three values: 0.32, 0.30, and 0.28. The results based on these input parameters are plotted in Figure 7. The right pane shows instability vs. dip values, while the left pane shows the Mohr diagram for strike and dip combinations in the right pane. Though the Samalanga-Sipopok Fault has higher instability in these figures, we could see that the Panteraja Fault also shows increasing instability values at dips of 40°–70° regardless of the friction coefficient. Optimum dip values for the Panteraja Fault to undergo a slip in this case might lie in dips of ≈40°–70°. Moreover, the Pidie Jaya aftershock distribution shows a dip of 61.2° ± 1.6°, which lies inside the optimum dip values. The higher instability values of the Samalanga-Sipopok Fault might be related to the slip of preceding events, e.g., the mb 6.1 1967.

Instability, as shown in Figure 7, is plotted with varying friction coefficients but constant SH values. Hence, the Mohr circle is intersected at two points with the failure line. In real conditions, the Mohr circle will decrease and become tangent only to the failure line at one point. Therefore, we again calculated the instability of the Panteraja Fault with a smaller SH magnitude, utilizing Eq. 9: sin ⁡ ψ =   ( 1 / 2 ) ( SH − Sh ) ( 1 / 2 ) ( SH + Sh ) , (9) where ѱ is the arcus tangent of the friction coefficient, which we varied into 0.32, 0.30, and 0.28. We found that decreases in SH magnitude will cause a slight increase in the peak of instability curve of the Panteraja Fault, as shown in Figure 8.

Fault instability was calculated on both nodal planes of the Pidie Jaya aftershocks in Figure 5. We grouped the nodal planes into NE-SW and NW-SE based on their strike orientation. Using the nodal plane parameter, we calculated instability for each event and plotted these in the Mohr diagram in Figure 9A. The majority of the aftershocks have higher instability for NE-SW nodal planes, except for events 9 and 12 (Figure 9B)). In the Mohr diagram, these NE-SW nodal planes are clustered mostly near the left inner circle where the mainshock would be, as seen in Figure 7. We also calculated ΔCFS as imparted by the mainshock to the NE-SW and NW-SE aftershock nodal planes. Most of the aftershocks experienced stress change in the range of ±1 bar. We plotted the ΔCFS results in Figure 9C. The friction coefficient used to generate results in Figure 9 is 0.32. We only used µ 0.32, because we found that smaller friction coefficients (0.30 and 0.28) in this case did not flip the result. The list of nodal planes and the results are in Supplementary Table S1.