Fault Instability and Its Relation to Static Coulomb Failure Stress Change in the 2016 Mw 6.5 Pidie Jaya Earthquake, Aceh, Indonesia

Herein, we applied the fault instability criterion and integrated it with the static Coulomb stress change (ΔCFS) to infer the mechanism of the 2016 Mw 6.5 Pidie Jaya earthquake and its aftershock distribution. Several possible causative faults have been proposed; however, the existence of a nearby occurrence, the 1967 mb 6.1 event, created obscurity. Hence, we applied the fault instability analysis to the Pidie Jaya earthquake 1) to corroborate the Pidie Jaya causative fault analysis and 2) to analyze the correlation between ΔCFS distribution imparted by the mainshock and the fault instability of the reactivated fault planes derived from the focal solution of the Pidie Jaya aftershocks. We performed the fault instability analysis for two possible source faults: the Samalanga-Sipopok Fault and the newly inferred Panteraja Fault. Although the maximum instability value of the Samalanga-Sipopok Fault is higher, the dip value of the Panteraja Fault coincides with its optimum instability. Therefore, we concluded that Panteraja was the causative fault plane. Furthermore, a link between the 1967 mb 6.1 event and the 2016 Mw 6.5 earthquake is discussed. To analyze the correlation between the fault instability and the ΔCFS, we resolved the ΔCFS of the Pidie Jaya mainshock on its aftershock planes and compared the ΔCFS results with the fault instability calculation on each aftershock plane. We discussed the possibility of conjugate failure as shown by the aftershock fault instability. Related to the ΔCFS and fault instability comparison, we found that not all the aftershocks have positive ΔCFSs, but their instability value is high. Thus, we suggest that the fault plane instability plays a role in events that do not occur in positive ΔCFS areas. Apart from these, we also showed that the off-Great Sumatran Fault (Panteraja and Samalanga-Sipopok Faults) are unstable in the Sumatra regional stress setting, thereby making it more susceptible to slip movement.


INTRODUCTION
According to the Coulomb failure criterion (Sibson, 1985;Oppenheimer et al., 1988;Zoback, 2007), failure on a plane of rock can happen if shear stress acting on the plane exceeds the failure resistance. Failure-resistance components include the rock cohesive strength and internal friction multiplied by normal stress acting on the plane. This basic concept is widely used by researchers to obtain two main methods, which relate to 1) fault reactivation potential or fault instability (Vavryčuk, 2011;Leclère and Fabbri, 2013) and 2) earthquake interaction (King et al., 1994). Both methods are important for gaining knowledge about earthquake occurrences. Fault reactivation potential estimates the ratio in which a fault plane is optimally oriented, based on the Mohr-Coulomb diagram (Sibson, 1985;Leclère and Fabbri, 2013). The ratio is high when a fault plane is optimally oriented and low when the fault plane is disoriented. The potential of fault reactivation was introduced for the first time by Sibson (1985) and further developed by other researchers, each using a unique name for this phenomena: slip tendency (Morris et al., 1996), reactivation-tendency (Tong and Yin, 2011), fault instability (Vavryčuk, 2011), and 3D fault reactivation (Leclère and Fabbri, 2013). Earthquake interaction is explained using a Coulomb failure assumption in a well-known method, namely, the Coulomb Failure Stress Change (ΔCFS) (King et al., 1994). Based on calculating stress changes before and after an earthquake, an increase in ΔCFS is associated with subsequent events, while a decrease in ΔCFS is related to the stress shadow effect (King et al., 1994).
The Pidie Jaya earthquake struck the Pidie Jaya district in Aceh, Indonesia, on December 7, 2016, at 05:03:33 local time. It was followed by aftershocks which were detected up to one month after the event. Regional and local networks in Indonesia recorded these aftershocks (Supendi et al., 2017;Muzli et al., 2018). The regional stations of the Agency for Meteorology, Climatology and Geophysics (BMKG) recorded large aftershocks, while the Pijay-Net temporary local stations, deployed jointly by BMKG, German Research Center for Geosciences (GFZ-Potsdam), Institut Teknologi Bandung (ITB), and Syiah Kuala University (UNSYIAH), also recorded aftershocks including those having small magnitudes. Moreover, global networks, such as the United States Geological Survey (USGS) seismic stations, also detected some aftershocks. However, the aftershocks from these three catalogs show an ambiguous fault plane orientation, as depicted in Figure 1.
Two catalogs show a northeast-southwest aftershocks trend (Muzli et al., 2018;USGS), while one catalog depicts a northwestsoutheast aftershocks trend (Supendi et al., 2017). By using the fault instability method, Sahara et al. (2019) conclude that the Pidie Jaya fault plane orientation has a northeast-southwest aftershock trend which is more unstable in Sumatra regional stress.
FIGURE 1 | Mainshock and aftershocks (BMKG, Pijay, and USGS network) of the Pidie Jaya 2016 earthquake. Focal mechanism (mainshock) is from BMKG. The relocated BMKG aftershocks (December 7-19, 2016) from Supendi et al. (2017) are plotted with a yellow pentagon. Relocated Pijay aftershocks (December 14, 2016-January 15, 2017) from Muzli et al. (2018) are plotted with blue circle. USGS aftershocks (December 6-17, 2016) are plotted with a gray square. The relocated Pijay network aftershocks are creating planar fault that dips 63°toward southeast (Muzli et al., 2018); white stars are relocated past events from Hurukawa et al. (2014). Samalanga-Sipopok Fault is taken from Barber et al. (2005). Panteraja Fault is taken from Muzli et al. (2018). Seulimum and Aceh Faults are taken from the National Center for Earthquake Studies of Indonesia (PUSGEN). The figure was modified from Kusumawati et al. (2019b) and Sahara et al. (2019 Along the east side of the 2016 Mw 6.5 Pidie Jaya earthquake lies a nearly 180°strike fault, the Samalanga-Sipopok Fault (Barber et al., 2005). The historical 1967 mb 6.1 earthquake took place between the Samalanga-Sipopok Fault and the Pidie Jaya earthquake (Hurukawa et al., 2014;Muzli et al., 2018). Hurukawa et al. (2014) relocate past events in the Great Sumatran Fault (GSF), including the 1967 mb 6.1 event, using Modified Joint Hypocenter Determination. After the relocation, the 1967 mb 6.1 event moved closer to the Samalanga-Sipopok Fault, as seen in Figure 1. This suggests that the fault is active; thus, it should not be overlooked.
Based on the BMKG report, the suspected causative fault for the 2016 Pidie Jaya earthquake is the Samalanga-Sipopok Fault (Djatmiko, 2016). This fault has been well documented in prior studies (Keats et al., 1981;Cameron et al., 1983;Genrich et al., 2000;Barber et al., 2005). Coming later, however, Muzli et al. (2018) did not associate the event with this fault due to the fact that neither of its focal mechanism strikes aligns with the strike of the Samalanga-Sipopok Fault. Instead, Muzli et al. associated it with an unidentified sinistral fault, which they suggest as being either the newly inferred fault in the west, the Panteraja Fault, or the same fault that was responsible for the 1967 mb 6.1 event. Nevertheless, some obscurity has arisen: the relocated 1967 mb 6.1 event is in the proximity of the Samalanga-Sipopok Fault, suggesting a link to that fault. In order to elucidate this obscurity, we proposed a fault instability analysis of the Pidie Jaya earthquake's possible causative faults to analyze their relation to the 1967 mb 6.1 event.
As previously mentioned, the ΔCFS could show the likely region of subsequent events by the stress increase area. However, previous studies have shown that some aftershocks were located in the stress shadow areas (Kusumawati et al., 2019a;Kusumawati et al., 2019b). Despite that, stress increases along the order as small as 1 bar were reported, which could trigger subsequent events (Harris, 2000). The triggered events are assumed to be optimally oriented planes (King et al., 1994). Optimally oriented planes mean that the fault planes are unstable in the current regional stress field. The fault instability method could quantify fault plane stability. Thus, it might be insightful to perform a fault instability analysis prior to ΔCFS. Such an analysis might give insight into why ΔCFS "succeeded" in one area or "failed" in another area.
The second objective of this study is to calculate the fault instability of fault planes and then compare these to ΔCFS of the major earthquake as resolved on those planes. We applied the second objective to the M6.5 Pidie Jaya earthquake and its aftershocks.
To achieve these study objectives, we first describe fault instability (Vavryčuk, 2011) and ΔCFS (King et al., 1994) methods. Fault instability depends on the regional stress field, especially its orientation. Thus, preceding to the main analysis, we discuss the Sumatra regional stress and its stress perturbation possibilities. Then, we analyze the fault instability on faults related to the first objective and on aftershock fault planes related to the second objective. The results and discussion are included in the last section.

Fault Instability and Coulomb Failure Stress (ΔCFS)
Brittle failure of a rock under triaxial stress is governed by the Coulomb failure criterion (Sibson, 1990): 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): CFS |τ| + μ σ ″ + P − C.
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).
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: 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 Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 8 | Article 559434 4 Figure 2). σ 1 , σ 2 , and σ 3 (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: where σ and τ are 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- 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 Parameter
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).

Parameters of Panteraja and Samalanga-Sipopok Faults
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.

The M6.5 Pidie Jaya Earthquake Slip Model and the Focal Mechanisms of the Aftershocks
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)  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).

Sumatra Regional Stress Analysis
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.

Fault Instability of Panteraja and Samalanga-Sipopok Faults
We calculated the fault instability of the Panteraja and Samalanga-Sipopok Faults using approximated strikes of 43.6°F IGURE 4 | Strike and dip estimation of Panteraja Fault based on the 2016 Pidie Jaya aftershock distribution (Muzli et al., 2018 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:  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 and ΔCFS in M6.5 Pidie Jaya Aftershocks
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.

DISCUSSION
The seismic history of the northern Sumatra region was analyzed by Hurukawa et al. (2014), who gathered and relocated historical earthquake data along the Sumatran Fault up to 2012. In the northern part of Sumatra, events from 1935 to 2012 were detected in the Hurukawa Catalog. The 1967 mb 6.1 earthquake occurred near the 2016 Pidie Jaya earthquake. The 1967 earthquake moved closer to the Samalanga-Sipopok Fault after relocation. We suggest the 1967 event to be related to Samalanga-Sipopok Fault activity. Other than the 1967 event, there is no historical M ≥ 6.0 earthquake recorded by the catalog around the Samalanga-Sipopok or Panteraja Faults (Supplementary Figure S1).
Based on GPS observations, the northern part of GSF has slip rates ranging between 16 and 20 mm/year (Ito et al., 2012). By using these slip rates as estimated slip rates for the Samalanga-Sipopok and Panteraja Faults, we were able to estimate the accumulated slip on these faults. The mb 6. It is worth noting that, given the limited seismic catalog we have for northern Sumatra, 1.30-1.62 m is the possible lower limit of slip accumulation in the Panteraja Fault. Nonetheless, from these calculations, we could see that Panteraja Fault has a larger accumulated slip than Samalanga-Sipopok Fault; hence, it is building up more stress than the Samalanga-Sipopok Fault, making it more susceptible for an earthquake occurrence.
A fault with a higher instability value is more susceptible to undergo slip. From Figure 7, we see that the Samalanga-Sipopok Fault has higher instability values than the Panteraja Fault. However, the Panteraja Fault has a unique fault instability feature. The instability curve peaked at a dip range between 40°and 70°, reaching the maximum values. This dip range might be revealing the most optimum dip for the Panteraja Fault to have an earthquake; the reactivated structure for the Pidie Jaya 2016 earthquake has a dip value of 61.2°± 1.6°, which coincides with the optimum dip range. Based on the seismic history and the fault instability analysis, we suggest that Panteraja was the causative fault plane of the 2016 Pidie Jaya earthquake.
In attempting to resolve the ΔCFS on the aftershock focal mechanism nodal planes, one should infer the fault plane carefully because if one nodal plane has a positive stress change, the other may have a negative stress change (Toda et al., 2011). Thus, Toda et al. (2011) suggest that the fault plane should be inferred, based on independent information such as seismic alignment. While for fault instability calculation, local redistribution of Coulomb stress might activate a cluster of events with slightly low instability values, as suggested by Vavryčuk et al. (2013) in their study using swarm events. Indeed, analysis of the Pidie Jaya aftershock fault instability and ΔCFS comparison should be carried out cautiously. Therefore, we used the aftershock alignments to further confirm the results.
Recalling the instability result in Figure 9B, most of the NE-SW nodal planes have higher instability, except for events 9 and 12 which have slightly smaller instability at the NE-SW orientation. Then, we tried to plot the aftershocks focal mechanisms with the NE-SW nodal planes highlighted, in Figure 10. We also colored the focal mechanisms in this figure based on the ΔCFS result for the NE-SW nodal plane. Particularly for events 9 and 12, we plotted the NW-SE nodal plane and their ΔCFS results. Faults with higher instability are interpreted as faults which are more susceptible to failure. In the case of the Pidie Jaya aftershock focal mechanisms, we suggest that they failed in the NE-SW nodal plane. This is because the main fault of the mainshock ruptured in the NE-SW direction, as shown by the sharp alignment of local events and GPS observation. Gunawan et al. (2020) found that, by using the NE-SW fault plane orientation that extends to the offshore as the mainshock, the synthetic and observed displacement in nearby GPS stations fits very well. However, for events 9 and 12, FIGURE 8 | Instability of Panteraja Fault at dips of 40°-70°, using different friction coefficients and SH magnitudes. The Sv, Sh magnitudes, SH orientation, and strike of Panteraja Fault are kept constant: Sv 240 bar/km, Sh 180 bar/km, N39°E, and 43.6°, respectively. The stress magnitudes refer to Hennings et al. (2012), while the stress orientation refers to Mount and Suppe (1992 the failure plane might be flipped to NW-SE. The first week aftershocks recorded by the BMKG regional network show sparse NW-SE event orientation. Approximately 78.4% of these regional aftershocks are concentrated in depths of 0 -15 km. Events 9 and 12 are located at the northern tip of the local aftershocks; they have shallow depths and occurred eight days after the mainshock. We suggest that there may have been a possible conjugate failure at a shallow depth in the NW-SE direction (aside from the NE-SW main fault failure), in which events 9 and 12 might have failed. This might confirm the identified cracks which are once thought to be associated with the possible NW-SE mainshock fault plane (right-lateral faulting) (National Center for Earthquake Studies of Indonesia (PUSGEN), 2017) but later suspected as the secondary effect of the NE-SW mainshock fault plane (Gunawan et al., 2020).
The ΔCFS and fault instability comparison for each aftershock (marked as aftershock event numbers) is plotted in Figure 11. From the comparison, we found that fault instability and ΔCFS do not have similar trends. However, the comparison showed that some of the events FIGURE 11 | Comparison of instability and ΔCFS of Pidie Jaya aftershocks nodal planes. All events, except events 9 and 12, are plotted using their NE-SW instability and ΔCFS results. Both instability and ΔCFS calculation used a friction coefficient of 0.32.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 8 | Article 559434 11 with negative ΔCFS have high instability values (above 0.7). This might give insight regarding the occurrence of events in the negative ΔCFS. The Coulomb stress is not sufficient to trigger other events in the negative ΔCFS area, yet some events occurred there, as can be seen in Figure 10. These events might exist as they are optimally oriented beforehand, which is shown from their high instability value.

CONCLUSION
We have conducted a fault instability analysis of the Pidie Jaya earthquake and described its relation to ΔCFS distribution. By using the fault instability method for the Pidie Jaya causative fault analysis, we concluded that neither the fault responsible for the 1967 mb 6.1 event nor the Samalanga-Sipopok Fault is the causative fault; however, we strengthened the Panteraja Fault as the causative fault, due to its dip which coincides with the optimum dip range.
Aftershock events located in the stress shadows were found to have a quite high fault instability, which indicates that these are the reactivation of critically stressed fractures ( Figure 11). Therefore, we showed that fault instability of the preexisting fractures plays a role in inhibiting or promoting the occurrence of aftershock events.
The GSF is a well-documented active Fault in Sumatra. Besides the GSF, northern Sumatra is also exposed to off-GSF faults, i.e., the Panteraja, Samalanga-Sipopok, and Lhoksemauwe Faults (Barber et al., 2005). Several historical earthquakes, as well as recent earthquakes, were located around these faults, indicating that these off-GSF faults are also active. Moreover, our fault instability calculation has shown that the off-GSF (Panteraja and Samalanga-Sipopok) faults are unstable in the Sumatra regional stress setting. This is shown by the high instability value, thereby making it more susceptible for slip movement, stress build-up, and finally fault failure (earthquakes), which could occur in its return period (Hardebeck and Okada, 2018).

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
DK and DS contributed to modeling, result analysis, and manuscript writing. MM contributed to data preparation. SW, AN, MM, II, NP, and ZZ reviewed the manuscript. All authors have read and approved the final manuscript.