Variations in Forearc Stress and Changes in Principle Stress Orientations Caused by the 2004–2005 Megathrust Earthquakes in Sumatra, Indonesia

Coseismic changes in principal stress orientation in the northern Sumatra subduction zone due to two giant megathrust earthquakes there in 2004 and 2005 are estimated to investigate the in-situ stress. The two megathrust earthquakes, the 2004 Sumatra-Andaman and the 2005 Nias-Simeulue events, are both among the 11 largest earthquakes ever recorded. Previous studies have shown that these giant earthquakes perturbed the stress field in the Sumatra subduction zone enough to alter the principal stress directions there, and here we investigate whether these changes can be used to better understand spatial variations in stress along the subduction zone. We used 330 previously published focal mechanisms to estimate pre- and post-mainshock principal stress orientations in 3 outer forearc segments and assessed whether orientation differences were resolved and what they imply about the pre- and post-mainshock stress fields. Our results agree with previous studies in establishing that coseismic changes in stress orientation in the forearc are resolvable, and consistent with a low level of stress in the outer Sumatran forearc before the earthquake, with almost all the shear stress on the megathrust relieved in the 2004 and 2005 earthquakes. In this study, we reveal that both the stress orientations and coseismic changes in them exhibit along-strike variations, with a decrease in both the pre-mainshock stress and stress drop found in the rupture area of 2005 relative to that of the 2004 earthquake. The forearc segment between the 2004 and 2005 rupture areas, which coincides with a well-known megathrust rupture barrier beneath the island of Simeulue is observed to have a characteristic signature, with lower shear stress relative to the pre-mainshock stress field and higher shear stress relative to the post-mainshock stress field in the adjacent segments.


INTRODUCTION
The island of Sumatra, Indonesia, is adjacent to an active plate boundary, where oceanic lithosphere of the Indian and Australian Plates subducts beneath the continental lithosphere of the Sunda Plate ( Figure 1). The Sumatra subduction zone experiences 40-60 mm/yr of oblique convergence, the accommodation of which is partitioned between mostly thrust slip along the Sumatra Trench megathrust, and dextral slip along the Great Sumatra Fault that parallels the island's southwestern coast about 250 km landward of the Sumatra Trench axis (Hamilton, 1979;McCaffrey, 1992;McCaffrey et al., 2000;McCaffrey, 2009;Malod et al., 1995;Moeremans and Singh, 2015).
In Figure 1 the relative motions of the Indian and Australian Plates with respect to the Sunda Plate are indicated for the GSRMv2.1 plate model (Kreemer et al., 2014). Because the oceanic lithosphere subducting off northern Sumatra is part of the wide, diffuse boundary between the Indian and Australian Plates that is actively deforming (Wiens et al., 1985;Gordon et al., 1990;Kreemer et al., 2014), GPS-derived plate models do not necessarily provide useful estimates of the obliquity of plate convergence. We believe that GSRMv2.1 is the most appropriate estimate of relative plate motion to use in our study of forearc stress, because it accounts for internal deformation of both the diffuse Indian-Australian plate boundary as well as that of the Sunda Plate that includes Sumatra. As can be seen in Figure 1, the obliquity of convergence of both the Indian and Australian with respect to the Sunda Plate exhibited by GSRMv2.1 is pronounced, at 45-52°. This is consistent with the 20-25 mm/yr dextral slip inferred on the Sumatra Fault by GPS and geologic studies Sieh et al., 1999, respectively).
The oblique plate convergence and the accompanying strain accumulation along the Sumatra forearc and the diffuse Indian-Australian plate boundary have resulted in the occurrence of major earthquakes off northern Sumatra during 2004-2012, including the 2004 Mw ∼9.3 Sumatra-Andaman (Stein and Okal, 2007) and 2005 Mw 8.6 Nias-Simeulue (Konca et al., 2007) megathrust earthquakes, as well as the 2012 Indian Ocean earthquake doublet (Mw 8.6 and Mw 8.2, Satriano et al., 2012). Three of these earthquakes are among the largest 11 earthquakes ever recorded (U.S.Geological Survey, 2017), and the stress transfer caused by such large earthquakes can change the stress field in the Earth's lithosphere as they relieve and redistribute stress over a wide area (Pollitz et al., 2006b;Gunawan et al., 2014).
Stress changes caused by large earthquakes can also promote or inhibit seismic activity in the surrounding area (Freed, 2005), and indeed changes in seismic activity following the 2004 Sumatra-Andaman and 2005 Nias-Simeulue earthquakes have been linked to such changes in the stress field (Wiseman and Bürgmann, 2011;Sevilgen et al., 2012). Some studies of earthquake-induced stress changes have also described how they can influence the earthquake cycle of major faults (Michael, 1987;Lin et al., 2011;Hasegawa et al., 2012;Hardebeck and Okada, 2018). (McCloskey et al., 2005) showed that the stress changes due the 2004 Sumatra-Andaman earthquake likely brought forward the next rupture of the adjacent Nias section of the Sumatra megathrust, which subsequently did rupture in the 2005 earthquake, and of the northern part of the Great Sumatran Fault (GSF), which has yet to rupture in a major earthquake. The suggestion of most of these studies is that many subduction zone megathrusts are subject to relatively low background differential stress that can experience near-complete relaxation due to the stress drop associated with large megathrust earthquakes (Hardebeck, 2012).
Most studies reporting stress changes after large earthquakes infer them from rotations of principal stress axis orientations estimated before and after the earthquakes (Michael et al., 1990;Hardebeck and Hauksson, 2001;Hardebeck, 2012;Yoshida et al., 2015). The orientation of principal stress axes can be estimated by inverting a set of focal mechanisms or moment tensors (Michael, 1987) from a large number of earthquakes. Calculating the stress orientations before and after a significant earthquake can reveal how the earthquake perturbed the pre-seismic stress field.
In the present study, we first apply stress inversion to yield the orientation of principal in-situ deviatoric stress axes before and after the 2004 and 2005 megathrust earthquakes. We subdivide the megathrust segmentation of Hardebeck (2012) on the northern part of Sumatra Island according to the seismicity pattern to infer the stress perturbation along strike. We use the iterative joint inversion approach of Vavryčuk (2014) to estimate principal stress orientations and apply it to different segments of the Sumatra subduction zone's outer forearc, to understand spatial variations in the pre-and post-mainshock stress fields. The bootstrap resampling technique is applied to assess the uncertainty of the solution, and results for coseismic changes in stress orientation are compared to previous studies of forearc stress and megathrust earthquake activity (Briggs et al., 2006;Meltzner et al., 2012).

Seismicity Data
The change in seismic activity before and after the 2004-2005 Sumatra megathrust earthquakes is illustrated by the cumulative seismicity curve in Figure 2. We divided this activity into two parts: the pre-mainshock seismicity which has a constant rate (from 1907 to 2004), and the post-mainshock activity characterized by rapid increases in seismic activity following the 2004 earthquake. We remove the 2005 Nias and 2012 Indian Ocean earthquakes from our analysis of the postmainshock activity, which includes the aftershock sequences associated with each earthquake. We regard the pre-mainshock seismicity as reflecting failure of faults in response to the stress field prevailing in the latter part of the megathrust earthquake cycle, while the post-mainshock seismicity reflects fault failure in response to stress in the early part of the cycle, following release of stress on the megathrust by a great earthquake.
The interpretation that the pre-2004 seismicity reflects a premainshock stress field that has gradually accumulated over the interseismic period seems straightforward, since no earthquake large enough to perturb the regional stress field has occurred in northern Sumatra for at least 100 years prior to the 2004 Sumatra-Andaman event (no event with Mw ≥ 6 was ever recorded in northern Sumatra prior to 2004, see Storchak et al., 2015). The interpretation of the post-2004 seismicity as reflecting the stress field in the early part of the earthquake cycle is more open to question, since the release of stress during the 2004 Sumatra-Andaman and 2005 Nias-Simeulue earthquakes occurred in a manner that was heterogeneous in time and space. The 2005 Nias-Simuelue earthquake occurred only 3 months after the 2004 Sumatra-Andaman, so the aftershocks that occurred prior to the 2005 event may reflect a different stress field from that which prevails after the 2005 event, and a stress inversion that includes them may be biased (see Figure 2). However, we believe this influence is mitigated by these potentially biased aftershocks having occurred well north of the rupture area of the 2005 event, so any bias due to this "transitional" stress field between the times of the 2004 and 2005 events should be minor. A similar consideration applies to the 2012 Indian Ocean earthquakes.
An additional complication is the significant rotation of the stress orientations observed during the months following the 2004 event (Hardebeck, 2012). This was due to a combination of the nearly complete relaxation of stress following the 2004 event-if stress is close to zero even a small stress recovery can change the principal stress directions-and post-seismic motion (Pollitz et al., 2006a;Gunawan et al., 2014). Moreover, significant post-seismic relaxation continued at least until 2010 and is likely ongoing (Gunawan et al., 2014). While it seems certain that these effects influence our estimate of the post-mainshock stress field using aftershocks during 2005-2017, we believe this estimate is nevertheless a useful measure of the average stress at the beginning of the interseismic part of the earthquake cycle, since the 13 years over which our estimate is averaged is much smaller than the previous interseismic intervals -550 (Rubin et al., 2017) and 143 (Newcomb and McCann, 1987)  Moment Tensor (GCMT) Project (Dziewonski et al., 1981;Ekström et al., 2012) as well as NIED CMT solutions (Kubo et al., 2002) and solutions from the United States Geological Survey as published in the Bulletin of the International Seismological Centre (Lentas et al., 2019). These data were divided into two-time windows:  Figure 2 is 5.0. Following Triyoso et al. (2020), the magnitude completeness (M c ) in Sumatra subduction is also 5.0. Therefore, we are sure that the catalogue completes at this level for the whole period of time. The potential bias in the figure refers to the aftershocks of the 2004 earthquake that could affect the stress field on the 2005 earthquake rupture zone (or in this case segment C).
FIGURE 3 | Focal mechanism distributions in northern Sumatra for the pre-mainshock (A) and post-mainshock seismicity (B). The yellow polygons show the segments used to group focal mechanisms for estimation of stress parameters in outer forarc. Background bathymetry/topography is from GEBCO (Kapoor, 1981 The spatial distributions of focal mechanisms used for premainshock and post-mainshock stress inversions are shown in Figure 3. Because in this study we attempt to resolve arc-normal variations in stress, we have segmented the data as indicated in Figure 3. The choice of segments used to group the data was guided by the clustering of aftershocks, and the rupture areas of the 2004 and 2005 mainshocks as described by finite fault rupture models (we considered Ammon et al., 2005;Banerjee et al., 2007;Konca et al., 2007;Tanioka et al., 2006). Spatial segmentation of the data is more consistent with the assumption in stress inversion that the stress field is uniform in the region over which focal mechanisms are inverted, but this comes at the expense of increased uncertainty due to the reduction in the number of events per segment. For this reason, we believe it is not practical to subdivide the data further in either time or space.

Stress Inversion
The methods for estimating principal stress parameters from focal mechanism data assume that: 1) the tectonic stress is uniform in the region over which focal mechanisms are inverted, 2) earthquakes occur on pre-existing faults with varying orientations, 3) the fault slip vector points in the direction of shear stress on the fault (Wallace, 1951;Bott, 1959), and 4) the earthquakes do not interact with each other and do not disturb the background tectonic stress (Vavryčuk, 2014). As discussed above, by grouping focal mechanisms into segments along the subduction zone forearc, we have reduced the region over which focal mechanisms are inverted and should thereby better satisfy the condition of uniform stress. Additionally, by dividing the data into pre-mainshock and post-mainshock intervals, and removing the mainshocks themselves, we have greatly reduced the potential for earthquake interaction and excluded earthquakes that could significantly disturb the overall stress field.
A simple approach for stress inversion is the method proposed by (Michael, 1984) which employs the expressions for shear traction on a fault in the following form: Where the sum over indices i, j and k is implied, δ ik is the Kronecker delta, N is the unit vector in the direction of shear stress resolved along the fault plane, τ is the deviatoric stress tensor, n is the unit vector normal to the fault plane and τ is the magnitude of shear stress on the fault. Eq. 1 can be expressed in matrix form as: Where t is a vector of the five independent components of the deviatoric stress tensor τ (i.e., τ 33 −(τ 11 + τ 22 )), A is a 3 × 5 matrix calculated from the fault normal n, and s is the shear traction on the fault plane. It is further assumed that this shear traction has close to the same value for all studied earthquakes. In the inversion the magnitude of s is normalized to have the value 1, since the method cannot determine absolute stress values. This equation can be solved using generalized linear inversion in the L2-norm as follows (Lay and Wallace 1995, their section 6.4): The linear system of Eq. 3 is inverted for all earthquakes simultaneously, resulting in estimated values for the five independent components of τ subject to the constraint that |s| 1. The four remaining stress parameters can be solved for three angles expressing the mutually orthogonal directions of the principle stresses, along with the "shape ratio" R: where σ 1 , σ 2 and σ 3 are the maximum, middle and minimum principal stresses, respectively. However, one basic disadvantage of applying this method is the lack of knowledge of which of the fault mechanism nodal planes is the actual fault plane. An incorrect nodal plane can degrade the accuracy of stress estimates (Vavryčuk, 2014). Therefore, (Vavryčuk, 2014), introduced a modified stress inversion method that iteratively inverts for both stress and fault orientation, with the latter determined using the fault instability criterion of (Lund and Slunga, 1999). (Vavryčuk, 2014) shows that this approach results in improved stress estimates, particularly for the shape ratio R.
Because the choice of fault plane made during the course of iterative inversion is based on the frictional instability criterion of (Lund and Slunga, 1999), the resulting stress orientation is dependent on a value for fault friction coefficient chosen a priori. In order to assess the sensitivity to and make an optimal choice for the fault friction coefficient, the inversion was run repeatedly with friction coefficient values ranging from 0.0 to 1.0, in steps of 0.05, and the optimal value was chosen as that which produced the highest overall fault instability (see Vavryčuk, 2014).
Finally, we calculate the uncertainty using the bootstrap method with 250 resampled datasets, as this approach is widely used to assess uncertainty in stress inversion (see e.g. Michael, 1984, andMichael, 2006). The best value is the mean of data distribution, and we express uncertainty as half the width of the 95% quantiles of the distributions associated with the bootstrap resampling.

RESULTS
In this study, we applied the iterative joint inversion scheme of (Vavryčuk, 2014) to invert for principal stress orientations and shape ratio, using 47 focal mechanisms for the pre-mainshock period and 282 focal mechanisms for the post-mainshock period.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712144 Stress inversions were performed for the 3 outer forearc segments A, B and C for the pre-and post-mainshock periods-segments A directly above the rupture area of 2004 Sumatra-Andaman earthquake, B at the gap between rupture area of 2004 and 2005 megathrust earthquakes, and C in the rupture area of the 2005 Nias-Simeulue earthquake (see Figure 3). The accuracy of stress inversion depends on the number of focal mechanisms inverted, the noise level in the data (Vavryčuk, 2014), and the distribution of different focal mechanisms. The directions of maximum horizontal compressive stress (S Hmax ) for the pre-and post-mainshock intervals are plotted in Figures 4A,B, respectively, and the numerical results for principal stress directions and shape ratio for these are indicated in Table 1 and 2. For each stress tensor, we calculated the directions S Hmax and their 95% confidence intervals following (Lund and Townend, 2007). Plots of the focal mechanism P-T axes, principal stress directions, and shape ratios are displayed in Figures 5A,B for pre-mainshock and post-mainshock time periods in the outer forearc segments, respectively. Histograms of bootstrap resampling distributions for σ 1 plunge and S Hmax are also shown in Figure 6.
We considered testing for bias introduced into the postmainshock stress solutions for segment C due to aftershocks of the 2004 earthquake (indicated as "Potential bias" in Figure 2), by removing such events from the Segment C stress solution. However, we found that the focal mechanism dataset for Segment C, after removal of megathrust events as described below, included no events during this time period between the 2004 and 2005 events (see  Supplementary Figure S1). Segment B contains four events between the time of the 2004 and 2005 earthquakes, but their effect on the estimated stress orientations was negligible.
To ensure the quality of the stress inversion, the diversity of the fault orientation has to be assessed. As the events were located in the subduction zone, we checked whether the orientations of thrust events are similar to megathrust orientation or not. If the focal mechanism data are dominated by megathrust events, the results of stress inversion may reflect the megathrust geometry rather than the orientation of the forearc stress field (Vavryčuk, 2011(Vavryčuk, and, 2014McKenzie, 1969). To avoid this potential source of bias, we identified and removed megathrust events from our dataset.
To identify megathrust events, we first estimated the strike and dip of the Sumatra megathrust for each segment. The strike direction of the megathrust for each segment is estimated based on the direction of the Sumatra trench axis and the megathrust dip of each segment is calculated based on the slab geometry of SLAB2.0 (Hayes, 2018), which is divided into 0-20 km and 20-50 km depth ranges. The strikes of the megathrust segments A, B, and C are 324°, 314°, and 328°, respectively. Dips estimated for each segment vary from 5°to 8°for the 0-20 km and about 20°for the 20-50 km depth ranges, respectively. We then compared the strikes and dips for the most shallowly-dipping nodal planes of the thrust events in our dataset, and rejected events with strike and dip within ±10°and ±5°, respectively, of the strike and dip of the megathrust for the respective segment. Figure 7 depicts the strike and dip values of the most shallow-dipping nodal planes for the thrust events, and shows which were identified and TABLE 2 | Comparison of pre-mainshock and post-mainshock stress inversion results. Angles are in degrees and measured positive clockwise (S Hmax calculated from the stress orientation parameters in Table 1 according to Lund and Townend, 2007 The uncertainty for all values is estimated as half the 95% confidence intervals. FIGURE 5 | Stress Inversion results for outer forearc segments A, B and C. Displayed are results for the pre-and post-mainshock stress field (A) and (B), respectively, in terms of P-T axes for individual focal mechanisms, principal stress axes, and shape factors. The uncertainty of principal stresses and shape ratios is estimated as half the 95% confidence intervals of bootstrap resampling method.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712144 rejected as megathrust events. Vertical cross sections representing the different geometries of the slab including rejected and accepted fault plane solutions in our stress inversion as well as the obtained stress axes are shown Figure 8. The stress inversion results displayed in Figure 5 are in overall agreement with those of (Hardebeck's, 2012), who estimated stress directions in the rupture area of the 2004 Sumatra-Andaman earthquake before and after the event. The directions of σ 1 and σ 3 lie in a vertical plane roughly aligned with the direction of plate convergence, with σ 1 having shallow plunge directed opposite to the megathrust dip, and σ 3 plunging steeply in the opposite direction. (Hardebeck's, 2012). used an outer forearc segment A that is almost identical to the segment A used here, so results there can be compared in detail. A comparison of the top panel of our Figure 5 with the bottom panel of (Hardebeck's, 2012) Figure 5C shows that the estimated post-mainshock stress directions in segment A are almost identical. Although the pre-mainshock stress directions are very similar, the bootstrap resampled distributions of our estimated stress directions exhibit more scatter, and our σ 1 has a plunge shallower than that of (Hardebeck 2012).
Because there were so few pre-mainshock events, the bootstrap resampled stress inversion showed larger scatter for each premainshock result compared to the corresponding postmainshock result, both for σ 1 plunge and S Hmax (Figure 6). The combined segment A-B-C produced smaller uncertainties for pre-and post-mainshock stress orientations. Although the stress inversions for individual segments have greater uncertainty, results for segments A, B and C reveal a pattern of along-strike variation in stress orientations, which is most clearly evident in the variation of σ 1 plunge ( Figure 6A).

DISCUSSION
In order to better understand the stress field of the northern Sumatra subduction zone, and how it was influenced by the giant megathrust earthquakes that occurred there in 2004 and 2005, we inverted focal mechanisms for the pre-and post-mainshock principal stress orientations in several segments of the outer forearc. The shallow plunge of σ 1 in opposite direction to the megathrust dip, and the steep plunge of σ 3 , are consistent with numerical models of low-stress outer forearcs (Wang and He, 1999). A plunge of σ 1 of 20°for the combined segments A-B-C (Table 1) implies a dip angle of about 40°with the megathrust, which is quite similar to other subduction megathrusts exhibiting high seismic coupling, like SW and NE Japan, and South America (Hardebeck, 2015). On the other hand, we find that plunges for individual segments are 11°, 29°and 29°for segment A, B and C, correspond to angles of 31°, 49°and 49°between the megathrust and σ 1 . Following (Hardebeck, 2015) assessment of variations in stress orientations and megathrust coupling in subduction zones worldwide, these results imply that segment A is strongly coupled, while segments B and C have weak or intermediate coupling.
The post mainshock stress field of the separate segments in the individual outer forearc is less consistent. Segments A and C, the like combined segment A-B-C, show a steeper plunge of σ 1 and clockwise rotation of S Hmax (see Table 2), which is consistent with the release of shear stress on the megathrust.
The exception to the near-complete release of pre-mainshock stress in the Sumatra outer forearc is segment B. Segments B and C both have an 18°steeper pre-mainshock plunge of σ 1 than segment A, indicating a weaker megathrust. However, relative to its pre-mainshock stress parameters, the plunge of σ 1 in segment B has actually decreased by 7°, in contrast to the 22°and 7°i ncreases of σ 1 plunge in segments A and C, respectively. Also, FIGURE 6 | Histograms of (A) σ 1 plunge and (B) SH max for outer forearc segments A-B-C, A, B and C estimated in this study. Each histogram represents the variation in samples taken from bootstrap resampling method.
TABLE 3 | Stress rotation (Δθ), the pre-mainshock angle of σ 1 to the fault (θ), and inferred stress drop ratios (Δτ/τ) according to Hardebeck and Hauksson (2001). Δτ/τ could not be resolved for segment B. segment B has experienced only a 8°clockwise rotation of S Hmax , which is much smaller than in either segments A or C (38°and 12°, respectively). While individual these differences may not be significant at the 95% confidence level, we believe the consistency in the contrasting behavior of segment B relative to A and C cannot be ignored. Using the approach of (Hardebeck and Hauksson, 2001), we infer the stress drop ratio, i.e., the ratio of stress drop to initial megathrust shear stress, using the estimated changes in σ 1 plunge (Δθ) and the pre-mainshock angle of σ 1 to the fault (θ) (see Table 3). This stress drop ratio is computed using Eq. 4 in (Hardebeck and Hauksson, 2001), and we take into account the uncertainty of the σ 1 plunge in our estimates of Δθ and θ. The intersection of the 95% quantile of boostrap resampling distributions for Δθ and θ with the allowable areas of the stress drop plot in Figure 9 encompass a wide range of stress drop ratios for each segment. However, if we consider that the occurrence of large earthquakes in segments A and C should be associated with a positive stress drop, then the stress drop ratio for segment A is in the range 86-100%, while that for C is in the range 0-89%. If we allow for segment B to have experienced no or partial rupture, then its range of possible stress drop ratios is so wide as to be essentially meaningless. The principal stress orientations and their rotations due to the occurrence of the mainshocks exhibited a consistent variation along strike, in which the pre-mainshock plunge of σ 1 was shallower and the coseismic rotations of both the σ 1 plunge and S Hmax were larger in segment A, the southern part of the 2004 earthquake's rupture area, than in segment C in the rupture area of the 2005 earthquake. The shallower σ 1 suggests a higher level of stress (Hardebeck, 2015;Hardebeck and Loveless, 2017) on the megathrust and the larger coseismic rotations indicate a larger stress drop ratio of the 2004 earthquake relative to the 2005 event. This apparent along-strike variation in stress and stress drops correlates with the time intervals since the previous earthquakes in these segments: the predecessor to the 2004 earthquakes in segment A occurred in AD 1450 ± 3 (Meltzner et al., 2012), so stress had built up over about 550 years, while the predecessor to the 2005 earthquake in segment C occurred in 1861, only 144 years earlier. Hence, megathrust stress accumulated on segment A over a longer interseismic interval before reaching failure, but when it does fail the stress drop is larger than in segment C. However, we note that paleotsunami studies indicate a wide variation in recurrence intervals of large earthquakes on the Sumatra megathrust (Rubin et al., 2017), so such a model is likely an oversimplification.
Like segment C, segment B showed a steeper pre-mainshock plunge of σ 1 which suggests that the accumulated stress on the Sumatra megathrust prior to the 2004 and 2005 mainshocks was less on Segments B and C than it was on segment A. The decrease in plunge of -7°following the 2004-2005 mainshocks implies that the stress on the segment B megathrust has increased as a result of the 2004 and 2005 earthquakes.
This apparently anomalous behaviour of Segment B can be explained by its position at the "gap" between the 2004 Sumatra-Andaman and 2005 Nias-Simeulue rupture areas, as indicated by the juxtaposition of the outer forearc changes in stress orientation  (Banerjee et al., 2007) for comparison, since they use exactly the same methodology for both earthquakes, and they optimize the fit to coseismic displacement observations that should reflect changes in stress.
In order to better identify the edges of the rupture area, we normalized each model with respect to its maximum slip, since otherwise the regularization used in the finite fault inversions can spread the rupture area beyond any rupture barrier. That such a barrier exists in segment B is well established by coral microatoll measurements of coseismic uplift along the coast of Simeulue north of Nias, which clearly show a 70 km-wide "saddle" of reduced uplift (0.5-1.5 m) that separates the rupture areas of the 2004 and 2005 earthquakes (Briggs et al., 2006, whose data are reproduced in Figure 10A). Moreover, (Meltzner et al., 2012), more detailed analysis of the coral microatoll data revealed that this barrier is persistent: over the past 1,100 years, rupture of at least three predecessors to the 2004 event and two predecessors to the 2005 event that uplifted the north and south coasts, respectively, of Simeulue, did not extend beyond this barrier. (Briggs et al., 2006) hypothesise that this barrier is a section of megathrust where slip occurs during the interval between large earthquakes, either aseismically or in moderate earthquakes. Thus, when a large earthquake occurs on either side of the barrier, rupture would be unlikely to propagate through the de-stressed section. Our measurements of stress orientations in segment B, especially when compared to those in segments A and C which lie in the rupture areas of the 2004 and 2005 events, respectively, are consistent with this hypothesis for the rupture barrier beneath Simeulue. The steeper plunge of σ 1 in segments B and C suggests that at the time of the 2004 event the stress was lower there than in the adjacent segment A where σ 1 had shallower plunge. While the megathrust in segment C ruptured at this lower level of stress, segment B did not completely rupture. Instead, the occurrence of the 2004 and 2005 earthquakes in the adjacent segments appears to have increased the stress on segment B, as evidenced by its reduction in σ 1 plunge. If the (Briggs et al., 2006) hypothesis is correct, this additional stress will be relieved in future years either through aseismic slip or the occurrence of moderate earthquakes.
Such an interpretation is supported by the aftershock study of (Tilmann et al., 2010), who noticed that shallow aftershocks of the 2004 and 2005 earthquakes mark the transition between stable . The obtained principal stress solution is also displayed with P (red circles) and T (blue plus) axes of every focal mechanism on the right of cross section for each time domain. The red focal mechanisms are the data that accepted in the stress inversion while the gray focal mechanisms are the rejected data as they identified as megathrust event. The slab used in the cross section is obtained from SLAB2.0 (Hayes, 2018).
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712144 FIGURE 9 | Plot of stress rotations in the three segments analysed in this study, where θ is the pre-mainshock angle of σ 1 axis to the orientation of the megathrust interface and Δθ is the stress rotation between the pre-vs post-mainshock. The theoretical stress drop ratio following Hardebeck and Hauksson (2001) for each segment is plotted. Bars plotted along with the theoretical stress drop is the uncertainty estimated as half the 95% confidence intervals (see Table 3).
FIGURE 10 | Comparison of (A) coseismic uplift data from Coral microatolls on the coast of Simeulue (Briggs et al., 2006) measured before and after the 2004 and 2005 earthquakes with (B) profiles of normalized slip (maximum slip as a function of latitude divided by the maximum slip for the respective event, taken from the finite fault models of Banerjee et al., 2007), and (C) distributions for changes in the pre-vs post-mainshock changes in the plunge of σ 1 and azimuth SH max among the outer forearc segments A, B and C.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712144 11 sliding of the megathrust near the trench and unstable sliding in the rupture zones downdip. However, exactly at the Simeulue Saddle this narrow band of aftershocks undergoes a pronounced shift down-dip, suggesting that stable sliding occupies a greater downdip extent of the megathrust in segment B than elsewhere. (Tilmann et al., 2010). note that this aftershock gap, and the Simeulue Saddle, occurs where a fracture zone in the subducting Indian Ocean crust intersects the trench, and they speculate that this could be a zone of enhanced fluid release into the megathrust, reducing its effective normal stress and promoting stable, aseismic slip (Moore & Saffer 2001;Scholz 2019).

CONCLUSION
We have investigated along-strike variations in stess of the Sumatra forearc, along with its perturbation due to the occurrence of the 2004-2005 megathrust earthquakes. Prior to the occurrence of these earthquakes, the outer forearc had an S Hmax well aligned with the oblique convergence of the Indian and Australian Plates relative to the Sunda Plate. The plunge of maximum principle stress σ 1 was only 11°in segment A, the rupture area of the 2004 earthquake, indicating a strong megathrust. This plunge increased to 29°in segments B and C, indicating a pronounced decrease in megathrust strength south of the 2004 earthquake's rupture area.
After the two megathrust earthquakes occurred, S Hmax rotated counterclockwise, to a direction more normal to the plate boundary, and the plunge of σ 1 increased, at least in segments A and C, which cover the rupture areas of the 2004 and 2005 earthquakes, respectively. This behaviour is consistent with previous studies and is consistent with a low-stress forearc, in which almost all the stress on the megathrust was released in the earthquakes (Herdebeck, 2012). The stress rotations in segments A and C indicate a larger stress drop ratio of the 2004 earthquake (in a range of about 85-100%) relative to the 2005 (in a range of about 0-90%).
The most intriguing results of this study were obtained for outer forearc segment B, which straddles the boundary between the 2004 and 2005 earthquakes, known as the "Simeulue Saddle" because it marks a gap in megathrust rupture that appears to have persisted over many earthquake cycles (Briggs et al., 2006;Meltzner et al., 2012). The post-mainshock stress results for segment B show that the plunge of σ 1 has decreased, implying that stress on the segment B megathrust has actually increased over its pre-mainshock state. This is exactly what would be expected for a rupture barrier that experiences aseismic slip: its pre-mainshock accumulated stress is less than the adjacent segments because stress has been released in the preceding interseismic interval by aseismic creep. Rupture in each of the 2004 and 2005 earthquakes has "loaded" the stress in segment B, as indicated by the reduced plunge of σ 1 , even though the aseismic zone itself fails to rupture. Under this hypothesis, the increased stress on the megathrust in segment B will be gradually released by aseismic creep until its stress is less than that recovered in the adjacent segments during the next interseismic interval, so that it returns to the pre-mainshock state prior to the 2004 earthquake.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
MR, PC, and DS contributed to calculate all the stress estimation. PC, DS, SW, WT, and AN contributed to supervise and validate the results. All authors contributed to the writing and the preparation of the manuscript. All authors have read and approved the final manuscript.