Stress accumulation and earthquake activity on the Great Sumatran Fault, Indonesia

The seismically active Sumatra subduction zone has generated some of the largest earthquakes in the instrumental record, and both historical accounts and paleogeodetic coral studies suggest these were large enough to transfer stress to the surrounding region, including the Great Sumatran Fault (GSF). Therefore, evaluating the stress transfer from these large subduction earthquakes could delineate segments of elevated stress along the GSF where large earthquakes may potentially occur. In this study, we investigated eight megathrust earthquakes from 1797 to 2010 and resolved the accumulated Coulomb stress changes onto 18 segments along the GSF. Additionally, we also estimated the rate of tectonic stress on the GSF segments which experienced large earthquakes. We considered two cases, with: (1) no forearc sliver movement, and (2) the forearc sliver movement suggested by recent studies. Based on the historical stress changes of large earthquakes and the increase in tectonic stress rate, we analyzed the time evolution of stress changes on the GSF. The Coulomb stress changes on the GSF due to megathrust earthquakes between 1797 and 1907 increased the Coulomb stress mainly on the southern part of GSF, which was followed by four major GSF events during 1890–1943. The estimation of tectonic stress rates using case (1) produces a low rate of stress accumulation and long recurrence intervals, which would imply that megathrust earthquakes play an important role in promoting the occurrence of GSF earthquakes. When implementing the arc-parallel sliver movement of case (2), the tectonic stress rates are much higher than case (1), with an observed slip rate of 15–16 mm/yr at the GSF consistent with a recurrence interval for full-segment rupture of 100–200 years. The case (2) result suggests that the occurrence of GSF earthquakes is dominantly controlled by the rapid arc-parallel forearc sliver motion. Furthermore, the analysis of the evolution of stress changes with time shows that some segments such as Tripa (North and South), Angkola, Musi and Manna, which have experienced full-segment rupture and are therefore likely locked, appear to have returned to stress levels similar to those prior to previous historical events, suggesting elevated earthquake hazard along these GSF segments.


Introduction
The Sumatran subduction zone is the world's most classic example of oblique, partitioned subduction (Fitch 1972;Jarrad 1986), in the sense that it involves a consistent regime of oblique convergence, with trench-parallel component partitioned between the Sumatra megathrust and the Great Sumatran Fault (GSF), extending over 1900 km. The GSF is a mature fault, with total offset of at least 20 km and about 15 mm/yr dextral slip (Sieh and Natawidjaja 2000), which cuts along the length of the island of Sumatra. The GSF has a history of generating large earthquakes, with over 19 earthquakes of magnitude 6.5 or greater having occurred since 1892 (Hurukawa et al. 2014), the largest being a multi-segment rupture of magnitude 7.6. However, no earthquake of magnitude greater than 7 has occurred since 1950, prior to Indonesia's explosive growth in population during the latter half of the twentieth century. Because the GSF is in some places proximate to dense population centers, GSF earthquakes pose a significant threat to lives and infrastructure, and it is therefore important to understand the fault's seismogenic potential.
The large GSF earthquakes in the early part of the twentieth century followed a suite of great Sumatra megathrust earthquakes that occurred mainly in the nineteenth century. While any connection between Sumatra megathrust earthquakes and GSF seismicity has not been firmly established, connections between megathrust and forearc activity have been established in other oblique subduction zones like SW Japan (Hori and Oike 1996;Mitogawa and Nishimura 2020;Shikakura et al. 2014), where the stress perturbations associated with the megathrust earthquake cycle can create temporal "windows" during which the stress regime is more favorable to the occurrence of forearc earthquakes. Such a mechanism was proposed for earthquakes along the northern part of the GSF following the Great 2004 Sumatra-Andaman megathrust earthquakes (McClosky et al. 2005). Despite the subsequent occurrence of several great earthquakes along the Sumatran megathrust during 2005-2010, this sequence has not yet been followed by any GSF earthquakes of magnitude greater than 6.6.
In this study, we considered the Coulomb stress accumulation on the GSF prior to the 2004 Sumatra-Andaman earthquake. We implement heterogeneous coseismic slip models of each large megathrust earthquake to calculate the coseismic stress perturbation, as well as a model for tectonic stress due to the coupling of the obliquely convergent plate motion at the Sumatra megathrust. We consider two cases: (1) no movement of the forearc sliver, so as to facilitate comparison with other forearcs in obliquely convergent margins and (2) with movement of the forearc sliver that is consistent with a recent analysis of Global Positioning System (GPS) measurements of surface velocities in the forearc . Using the historical stress changes of large earthquakes and the rate of increase in tectonic stress, we estimate the time evolution of stress changes on the GSF, with a view to better understanding its potential to produce large, damaging earthquakes.

Oblique subduction and forearc slivers: Sumatra compared with SW Japan
Almost all the world's subduction zones involve some degree of oblique convergence. These often exhibit convergence partitioning, with arc-perpendicular and some arc-parallel convergence taken up on the megathrust, while the remaining component of arc-parallel convergence is taken up by strike-slip faults in the upper plate, often near the volcanic arc (Fitch 1972;Jarrad 1986;McCaffrey 1992). The crustal block between this strike-slip fault system and the trench is called a forearc sliver and may migrate along strike in the direction of arc-parallel convergence. Almost all forearc slivers are buttressed at their leading edge, however, which hinders along-strike movement and leads to a dominantly trenchparallel compressive stress field in the forearc sliver (Beck et al. 1993;Wang 1996).
The Nankai Trough, where the Philippine Sea Plate converges with SW Japan at a rate of 59-65 mm/yr and obliquely of 20°, has one of the world's best-studied forearcs. Although similar in some ways to Sumatra, it is also a useful benchmark that is typical of oblique subduction zone forearcs worldwide. The arc-parallel component of convergence at the Nankai Trough is partitioned between the trough and a dextral fault system that cuts through SW Japan over a distance of over 500 km, from Kyushu to the Kii Peninsula, known as the Median Tectonic Line (MTL). The forearc sliver between the trough and the MTL migrates westward in the direction of arc-parallel convergence at a rate of 5-10 mm/ yr (Tabei et al. 2003). The subduction plate boundary rotates counter-clockwise offshore Kyushu, creating a geometric buttress (Beck et al. 1993) that is consistent with the arcparallel compression in SW Japan inferred from focal mechanism data (Wang 2000). The MTL is the longest onshore fault in Japan, consisting of 12 segments separated by gaps, bends, etc., and paleoseismic evidence indicates a recurrence interval of 1000-3000 years, with relatively high slip (5-8 m) events suggestive of multi-segment rupture (Tsutsumi and Okada 1996). Given that almost all forearcs are buttressed (Beck 1991;Wang 1996) and there are very few examples of large, strike-slip earthquakes on forearc-sliver-bounding faults worldwide, we regard this low rate of arc-parallel movement and long earthquake recurrence time as typical of forearc slivers in obliquely convergent subduction zones.
While the Sumatra subduction zone is similar to SW Japan in being an ocean-continent subduction zone with a similar rate and obliquity of convergence (about 41-62 mm/yr and 20°, respectively), its forearc sliver and the GSF that bound it differ in ways that make it unique among forearc slivers worldwide. At 1900 km length, it is the longest active strikeslip fault in the world, and its 18 + segments and total offset of 20 km (or possibly more, see Sieh and Natawidjaja 2000) fall in the middle of the range of major strike-slip faults (Stirling et al. 1996). Unlike almost all other forearc slivers, the Sumatran forearc is unbuttressed with its movement along the arc-parallel direction of convergence accommodated by spreading in the Andaman Sea (Curray 2005;Kimura 1986). This is supported by the forearc stress orientations estimated by Rafie et al. (2021), who found that the principle horizontal compressive stress in the Sumatra forearc is oriented in the direction of plate convergence, rather than being arc-parallel as it would be for a buttressed forearc. It was long thought that the slip rate along the GSF decreases from 25 mm/yr near the equator to only 6 mm/yr near the Sunda Strait (McCaffrey 1991), which seemed to be consistent with both geodetically and geologically measured slip rates (Genrich et al. 2000;Sieh and Natawidjaja 2000, respectively) and implied arc-parallel extension of the forearc sliver. However, more recent work on both geologically and geodetically determined slip rates 1 3 now suggest that these data can be fit with a uniform slip rate of about 15 mm/yr Natawidjaja et al. 2017). In particular, Bradley et al. (2017) show that the geodetic data are consistent with along-arc movement of the forearc sliver by 15 mm/yr, 1.5-3 times faster than the SW Japan sliver.
The greatest contrast between the GSF and the MTL, and other forearc slivers worldwide, is in its level of earthquake activity. While the MTL has not experienced a large, destructive earthquakes since at least 1600 AD (Tatsumi and Okada 1996) and possibly longer, the GSF has experienced at least 20 earthquakes of magnitude ≥ 6.5 (Table 1 and Hurukawa et al. 2014) since 1892, many of which caused widespread damage. When the greater fault length but shorter historical record is taken into account, this amounts to more than an order of magnitude greater number of large earthquakes per year and km of fault length on the GSF vs. the MTL. In Table 1, we have used the average slip-and rupture area-magnitude relationships of Wells and Coppersmith (1994) to estimate slip and rupture area for these large, historical GSF earthquakes. For those earthquakes that involve full-segment rupture (fractional area > 60%), the average slip varies between 1.8 and 3.3 m. Considering the 15 mm/yr average slip rate obtained by Bradley et al. (2017), a slip-predicable earthquake recurrence model would suggest recurrence times of 120-220 years, again an order of magnitude more frequent than paleoseismic estimates of recurrence times for the MTL.

Influence of Sumatra megathrust earthquakes on GSF earthquake activity
In low-stress forearcs like SW Japan and Sumatra, large megathrust events can cause stress perturbations sufficient to alter the orientations of principle stress axes (Rafie et al. 2020;Wang 2000) and thereby alter earthquake activity in the forearc. In SW Japan, changes in shear stress on the inland faults caused by great megathrust earthquakes drive temporal variations in seismicity in the SW Japan forearc, modulating the occurrence of large shallow crustal earthquakes in the forearc so that their frequency is greatest 50 years before and 10 years after the occurrence of megathrust earthquakes (Hori and Oike 1999;Mitogawa and Nishmura 2020;Shikakura et al. 2014). Similar perturbations to the forearc stress field due to the Great 2004 Sumatra-Andaman earthquake caused changes in seismic activity in the Andaman Sea (Sevilgen et al. 2012), and likely brought forward the next rupture of the northern part of the GSF (McClosky et al. 2005).
Sumatra forearc stress is likely to be influenced by an extensive history of large Sumatra megathrust earthquakes. Since the beginning of the twenty-first century, the Sumatra subduction zone has recorded some of world's largest instrumentally recorded megathrust earthquakes. Historical accounts describe a similar level of earthquake activity extending over the past 225 years (Newcomb and McCann 1987;see Fig. 1), and coral paleogeodetic and paleotsunami studies extend a similar record to the early Holocene (Philibosian et al. 2014;Rubin et al. 2017). A series of historical events occurred during the colonial period: 1797 Mw ~ 8.7 , 1833 Mw ~ 8.9 ), 1861Mw ~ 8.5-8.6 (Briggs et al. 2006Meltzner et al. 2015;Natawidjaja et al. 2006), and 1907Mw ~ 8.2-8.4 (Martin et al. 2019. No large Sumatra megathrust earthquakes were experienced through the rest of the twentieth Century, until a series of very large events at the beginning of the 21st: 2004 Mw 9.1 Sumatra-Andaman (Ammon et al. 2005;Rhie et al. 2007), 2005 Mw 8.6 Nias-Simeulue (Konca et al. 2007), 2007 Bengkulu (Konca et al. 2008), and 2010 Mw 7.8 Mentawai (Yue et al. 2014).
As described above, the GSF also has an extensive history of large earthquakes (Table 1 and Hurukawa et al. 2014;Natawidjaja and Triyoso 2007;Newcomb and McCann 1987). Soon after the initiation of a sequence of major megathrust earthquakes between 1797 and 1907, several large GSF earthquakes were observed (Fig. 2). The delays between the occurrences of major events along GSF and nearby (within 500 km) megathrust earthquakes are presented in Fig. 3. Following Natawidjaja et al. (2006), we divided the GSF into two areas: Table 1 Segment geometry and average slip and stress drop estimates for large (M ≥ 6.5) GSF earthquakes.
Magnitude estimates are taken from Hurukawa et al. (2014). Slip is determined using the average slipmagnitude relationship of Wells and Coppersmith (1994). Fractional area reflects the ratio of rupture area estimated from Wells and Coppersmith (1994) divided the area of the segment, assuming a downdip rupture width of 20 km. Stress drops are calculated using the formula for strike-slip faults of Kanamori & Anderson (1975), again assuming a downdip width of 20 km and a rigidity of 32 GPa. Creeping segments are identified from Ito et al. (2012), Natawidjaja (2018a) and Tong et al. (2018). Where a segment has not experienced and earthquake, knowledge of its seismogenic potential is indicated by "creeping," "unknown" or, in the case of Barumun, "stress shadowed," indicating it has been affected by the 1892 earthquake in the adjacent Angkola segment a This 1943 earthquake was multi-segment and fully ruptured both Sumani and Suliti segments b Full downdip width assumed to have ruptured when the fractional area > 50% northern (Seulimeum to Angkola) and southern (Sumani to Kumering). The histograms of the time intervals between the occurrence of each GSF earthquake in Table 1 and all preceding megathrust events that occurred with 500 km of the GSF event show that, for the northern GSF earthquakes, 11 out of 16 time intervals are 150 years or less. For the southern GSF earthquakes, all 25 time intervals are in the range 50-250 years. This clustering of time intervals between GSF and megathrust earthquakes suggests that megathrust events modulate earthquake activity on the GSF.

Data and fault models
The best way to study the interaction of the effects of megathrust strain accumulation and rupture on inland faults is by modeling changes in forearc stress, as has been considered by studies of the Nankai Trough (Hori and Oike 1999;Mitogawa and Nishmura 2020;Shikakura et al. 2014), Cascadia (Wang 2000), Chile (Ryder et al. 2012) and Caribbean (ten Brink and Lin 2004) forearcs.
To analyze the occurrence of large earthquakes along the GSF and investigate whether, and if so how, they might be influenced by major megathrust earthquakes, we used the catalogue of Hurukawa et al. (2014) to provide precise hypocenter locations of large earthquakes M ≥ 7.0 and their fault planes since 1892. The GSF was segmented into 19 segments by Sieh and Natawidjaja (2000), and these were used by Hurukawa et al. (2017) and other studies of the GSF. Although the National Center for Earthquake Studies of Indonesia (PUSGEN 2017) detailed the GSF segmentation into 41 sub-segments, we prefer the simpler 19-segment model, and we have excluded the Sunda segment (Fig. 2). Also, we The 19 segments of the Great Sumatran Fault (GSF) and the GSF earthquakes as described by Hurukawa et al. (2014), along with GSF trace data and subduction megathrust segments from Irsyam et al. (2020). Note that the bifurcations of the GSF in the north and south, where the Seulimeum/Aceh and Semanko West/East segments overlap, respectively, have been simplified to consider only the Seulimeum and Semangko West (referred here simply as Semangko) segments, since they are most favorably oriented with respect to the stress field. Also, the Tripa segment of Hurkawa et al. (2014) has been split into Tripa North and Tripa South segments, to better conform to the surface trace of the GSF as given in Irsyam et al. (2020). GPS vectors of Australian and Indian Plate motion relative to the Sunda Black are from Kreemer et al. (2014) 1 3 have ignored the short-thrust segment in the northernmost 30 km of the Tripa segment as specified by PUSGEN (2017), so that we can consider the entire GSF as a vertical strikeslip fault. Further refinement of fault geometry could be considered in future work, but for the present study we want to focus on a simple model for GSF earthquake activity.

Estimation of change in Coulomb failure stress due to great megathrust earthquakes
The Coulomb failure criterion describes the stress condition under which a pre-existing fault can rupture, and relates the shear stress on the fault to the normal stress (King et al. 1994). The static Coulomb Failure Stress change (ΔCFS) is often implemented to explain the spatial pattern of aftershock occurrence or subsequent triggered events caused by a large earthquake (e.g., King et al. 1994;Stein 1999;Toda et al. 1998Toda et al. , 2005. ΔCFS is defined as: where Δ is the shear stress change on the fault (positive in the slip direction), is the fault friction coefficient, Δ is the change in normal stress (positive in extension) and ΔP is the pore pressure change. Δ and Δ arise from the perturbation to the regional stress field caused by a large earthquake. The change in pore pressure is often related to mean stress change and assumed to be proportional to the normal stress change, so that Eq. (1) is usually replaced by: (1) ΔCFS = Δ + (Δ + ΔP)  Table 1 and all preceding megathrust events that occurred with 500 km of the GSF event. The origin of the horizontal axis, year 0, corresponds to the occurrences of the great Sumatra megathrust earthquakes in 1450, 1797, 1833, 1861, 1907, 2004, 2005, 2007 and 2010. GSF earthquakes have been split into northern GSF events, consisting of ruptures of the Seulimeum, Renun, Tripa and Angkola segments, and southern GSF events, consisting of ruptures of the Sumani, Sianok, Suliti, Siulak, Dikit, Ketaun, Musi, Manna and Kumering segments where ′ is called the apparent or effective friction coefficient which includes changes in pore pressure. Assuming the secular change in tectonic stress is such that the stress on a fault will accumulate to a level resulting in rupture, a positive ΔCFS will bring the fault closer to failure than it was prior to the large earthquake, while a negative ΔCFS will lengthen the time it takes for the accumulation of stress to lead to rupture (Harris 1998). We use the COULOMB 3.4 code (Lin and Stein 2004;Toda et al. 2005) to calculate ΔCFS using the elastic dislocation approach (Okada 1985(Okada , 1992 for eight historical Sumatra megathrust events. The Coulomb stress calculations are then estimated for the specific orientation of each of the 18 segments of the GSF considered here. The ΔCFSs caused by the megathrust earthquakes were calculated based on the slip distributions of historical megathrust earthquakes. The rupture model of each large earthquake was obtained from previous studies ( Table 2). The slip models of the 1797 and 1833 earthquakes were taken from Natawidjaja et al. (2006) which was constructed based on coral microatoll measurements. We take the slip distribution of the 2005 Nias-Simeulue as an analog for that of the 1861 earthquake, since both earthquakes have similar patterns of coseismic uplift measured from corals (Meltzner et al. 2015). We calculated the ΔCFS caused by these eight historical megathrust earthquakes on each GSF segment, which we henceforth denote as ΔCFS m . As has been presented in previous studies (e.g., Mildon et al. 2016 andMildon et al. 2019), the Coulomb stress calculation is most sensitive to the variation of strike. Therefore, for bending fault of Angkola and Barumun segments, we treated the segments by dividing their geometry into two-parts: from initial fault to the bending point and from bending point to the other end of fault.
(2) ΔCFS = Δ + � Δ 1 3 Cattin et al. (2009) studied how variations in seismicity in the Sumatra-Andaman-Sagaing fault system might be caused by the stress perturbations due to the 2004 and 2005 megathrust earthquakes, but were unable to constrain the effective friction of the GSF. The effective friction coefficient value we used for all segments was 0.1, because this is consistent with many studies of the SW Japan forearc (Iio 1997;Hori and Oike 1999;Shikakura et al. 2014;Mitogawa et al. 2020), as well as with studies of major crustal strikeslip faults (Carpenter et al. 2015;Provost et al. 2003;Townend and Zoback 2004). We calculated ΔCFS m at 10 km depth. This choice of depth is supported by previous studies which showed that the stress changes on the GSF are relatively constant from 10 to 30 km depth (Cattin et al. 2009;Qiu and Chan 2019). From the calculated ΔCFS m , we investigated the effect of stress transfer with the occurrence of major earthquakes along the GSF. The results for ΔCFS m calculated on the GSF segments due to each large Sumatran megathrust earthquake are shown in Fig. 4.
The ΔCFS m distribution on the GSF due to the 1797 earthquake (Fig. 4a) has negative values around the central part of the GSF, especially the Toru, Angkola and Barumun segments. Positive values for ΔCFS m of > 0.1 MPa are found southward from Sumpur to Dikit. As for the 1833 earthquake (Fig. 4b), similar positive and negative stress transfers southward and northward are found, respectively. The 1907 earthquake (Fig. 4c)

Estimation of tectonic stress rate
Apart from the ΔCFS m imparted by megathrust events, the tectonic stress rate also plays an important role in determining the Coulomb Failure Stress acting on a fault at any time. Here, we denote the tectonic stress contribution to CFS on the GSF as ΔCFS t . The higher the tectonic stress rate, the faster the fault can be brought back to its critical condition after rupture, and vice versa (Hori and Oike 1999). The tectonic stress acting on the GSF comes from two sources: (1) coupling of the oblique convergence at the megathrust, which transmits the stress acting on the megathrust through the forearc sliver and onto the GSF, and (2) arc-parallel motion of the forearc sliver itself, which exerts an additional traction on the GSF (Fitch 1972;Malod and Kemal 1996). Earlier studies (e.g., McCaffrey 2009) suggested there is a strong northward increase in the slip rate of the GSF, which would result in stretching of the forearc and imply a commensurate northward increase in tectonic stress rate. However, the more recent study of Bradley et al. (2017) uses GPS measurements of crustal velocities to show that the Sumatra forearc sliver undergoes little or no internal permanent deformation, but experiences arc-parallel movement that is consistent with an almost constant slip rate of 15 mm/yr along the GSF indicated by revised geologic estimates of the slip rate Natawidjaja et al. 2017).
We calculate ΔCFS t for two cases: (1) tectonic stress solely due to coupling of the oblique plate convergence at the megathrust, and (2) tectonic stress due to the megathrust coupling as well as sliver movement with respect to the Sunda plate as described by Bradley et al. (2017). The former, no-sliver-movement case, is useful for comparison with other forearcs, like SW Japan, that are buttressed and have less sliver movement. The latter, sliver movement case, is preferred because it better explains the available GPS data .
We derive coupling models for both cases, using the 3-block model of Bradley et al. (2017) for the transition between Indian and Australian plate motion with respect to the Sunda Plate, as illustrated in Fig. 5. For case (1) with no sliver motion, all of the oblique plate convergence between these blocks and the Sunda plate is applied as a "back-slip" dislocation on the megathrust (Savage 1983). For case (2), the arc-parallel motion of the forearc sliver is subtracted from the oblique plate convergence in order to determine the back-slip applied on the megathrust, and then, the forearc sliver movement is added to all the velocities calculated in the sliver. In both cases, we have allowed coupling in the four segments: Aceh-Nias, Batu, Padang and Enggano to vary, as indicated in Fig. 5. We used a grid search of coupling values over these four segments, with the preferred coupling values based on the minimum root mean square residual between calculated and observed GPS velocities (see Table S2 and Table S3). Following Bradley et al. (2017), we allowed  Hayes et al. (2018) are displayed, as is a comparison of modeled GPS velocities (green arrows) and observed velocities used for the kinematic modeling of Bradley et al. (2017, black arrows). (c) shows profiles of the megathrust coupling obtained by optimizing the fit of the GPS data for the two cases of no sliver movement and sliver movement in (a) and (b), respectively. In order to fit the data for each case, coupling values were varied in the Aceh-Nias, Batu, Padang and Enggano portions of the Sumatra megathrust the Aceh-Nias and Padang segments to be fully locked, but considered that the Batu and Enggano segments are likely to have low coupling , with coupling in fully locked and low-coupling segments allowed to vary in the ranges 0.0-1.0 and 0.0-0.4, respectively. Results for the preferred coupling values for each of cases (1) and (2) are shown in Fig. 5a, b, respectively. We note that in both cases, similar fits to the GPS data were obtained near the Padang and Batu segments, but only case (2), with sliver movement, results in a good fit for the data in the Aceh-Nias and Enggano segments.
In order to model the elastic deformation of the forearc sliver due to coupling of the plate convergence at the megathrust, an appropriate distribution of coupling is applied along the megathrust. We use the Triangular Dislocation Element (TDE) approach of Meade (2007) to calculate the elastic deformation due to back-slip applied on triangular elements used to tessellate the irregular shape of the SLAB2.0 model (Hayes 2018) for the Sumatra megathrust (see Fig. 5). We can estimate the stress tensor using Hooke's law, = tr( ) + 2 , where is the strain tensor caused by the back-slip dislocation, tr( ) is the trace of the strain tensor and is the identity matrix. The shear stress and normal stress can then be resolved on each segment of the GSF to calculate the associated ΔCFS t according to Eq. 2. For case 1, with no sliver movement, this corresponds to the total Coulomb stress change ΔCFS t 1 ( Table 3). As for case (2) with sliver movement, the whole calculation is similar to what we have used for no sliver movement; however, we use the forearc rigid block movement of Bradley et al. (2017) to calculate a slip deficit for each GSF segment (Δu s in Table 3). Then, we calculate the corresponding shear stress accumulation Δτ s based on the relation of Kanamori and Anderson (1975), and add this to the shear stress Δτ obtained from dislocation modeling for Sumatra megathrust coupling, to obtain the total ΔCFS t 2 for case (2) as: An important caveat to this approach is its failure to account for fault-normal sliver movement, as indicated in Table 3 by the angle θ between the forearc sliver and Sunda plate motion and at each segment. Except for the Tripa North segment, this angle is generally less than 20°, with mostly negative angles indicating extension in the south and positive angle indicating compression in the north, consistent the forearc sliver's pole location and sense of rotation . The fault-normal movement must be accommodated either by oblique motion on the GSF itself, or by faulting and/or folding in the forearc sliver. Since there is little or no information on the dip of the GSF segments and internal strain rates of the sliver are not detectable using currently available GPS data ), we do not account for this motion here.

Evaluation of ΔCFS m imparted by Sumatran megathrust earthquakes on the GSF
We estimated the accumulation of ΔCFS m imparted by the eight Mw > 8.0 megathrust events between 1797 and 2010 to the GSF (Fig. 4) and display the cumulative values of these on the different segments of the Sumatra fault in Fig. 6. The distribution of the cumulative ΔCFS m along the GSF exhibits high positive values in the north and south, but lower and even negative values in the central part of Sumatra. To analyze the effect of smaller (3) ΔCFS = Δ + Δ s + � Δ earthquakes (Mw ≤ 6.5) that occurred near the GSF segments, we calculated the ΔCFS of 1977 Mw 6.1 earthquake that ruptured on Angkola segment (see Figure S1 and Table S1). The ΔCFS was then imparted to the segments adjacent to the earthquake (i.e., Barumun and Sumpur segments). This event produces a maximum ΔCFS m of 5.7 kPa on a small area near the middle of the Barumun segment ( Figure S1), which might help trigger rupture if the segment is already stressed close to failure. However, the average ΔCFS m imparted by this event to the whole Barumun segment is only -0.7 kPa (-1.9 kPa for the Sumpur segment). This suggests that the influence of smaller events like this on the ΔCFS accumulation on whole segments is very small, and they are unlikely to be a major contributor to the segment-scale evolution of ΔCFS on the GSF.
The contribution to the cumulative ΔCFS m made by the major Sumatran megathrust events between 1797 and 1907 affected mainly the southern GSF, showing a cumulative ΔCFS m of 0.1 to 0.4 MPa on the Siulak-Manna segments, and a slight increase of 0.04 MPa on the northern GSF, particularly on the Renun and Toru segments. This high ΔCFS m in the southern part of GSF was followed by four major (magnitude ≥ 7) GSF events observed Table 3 Estimation of the rates of increase in tectonic stress on the GSF for the two cases considered: (1) no forearc sliver movement and (2) forearc sliver movement from Bradley et al. (2017) In each case, Δτ and Δσ are perturbations to shear and normal stress, respectively, on each GSF segment caused by the coupling of oblique convergence at the Sumatra megathrust. For case (1), these are combined according to Eq. (2) to yield ΔCFS t 1. For case (2), an additional term Δτ s , calculated from the slip deficit Δu s on each GSF segment associated with sliver movement according to Kanamori & Anderson (1975), is added to Δτ in order to calculate the change in Coulomb stress ΔCFSt 2 . θ is the angle between the forearc sliver and Sunda plate motion and at each segment from the Siulak to Manna segments between 1890 and 1943. During this same period, however, a similar number of major events occurred on the central and northern GSF, and also one on the Kumering segment far to the south, where ΔCFS m did not exceed 0.04 MPa and was in some places negative. Thus, while at first glance the correlation between the pre-1907 megathrust earthquakes and subsequent GSF earthquakes might seem high (e.g., see Fig. 2 Giacomo et al. 2018;Storchak et al. 2015) pattern following the 1833 Bengkulu earthquake. None of the twenty-first century Sumatra megathrust events have been followed by major GSF events, with the exception of the magnitude 6.6 earthquake in the Dikit segment, despite the expectation following the 2004 Sumatra-Andaman earthquake that a major GSF earthquake was imminent (McCloskey et al. 2005).
The values of ΔCFS m that we obtain for the northern GSF due to the 2004 and 2005 megathrust earthquakes are less than those of McCloskey et al. (2005) and Qiu and Chan (2019), who obtained maximum ΔCFS m of 0.9 MPa and 1.5 MPa, respectively. This is most likely due to a combination of two factors. First, our reported ΔCFS m of 0.4 MPa is the average for the Seulimeum-Aceh segment, which actually experiences a maximum ΔCFS of 0.7 MPa. Second, Qiu and Chan (2019) assumed a high effective friction values of � = 0.4 , while we used � = 0.1 ( McCloskey et al. 2005, do not say what value they used). As explained above, we prefer the low value ′ =0.1 along the GSF because it agrees with studies of the SW Japan forearc and of large, crustal strike-slip faults (Carpenter et al. 2015;Iio 1997;Hori and Oike 1999;Mitogawa and Nishimura 2020;Shikakura et al. 2014;Provost et al. 2003;Townend and Zoback 2004).

Time Evolution of ΔCFS on the GSF Without Sliver Movement
In order to model the time evolution of ΔCFS on any GSF segment, we need the results of ΔCFS m due to megathrust earthquakes, the tectonic stress rate ΔCFS t , and the stress drops of large intraplate earthquakes on nearby GSF segments. Table 1 shows the slip of each GSF earthquake calculated using the relationship between magnitude and Average Displacement (AD) of Wells and Coppersmith (1994) based on a linear regression of global magnitude and slip data. We also used the relation of Kanamori and Anderson (1975) to calculate stress drops of all large (M ≥ 6.5) GSF earthquakes, assuming a down-dip rupture width of 20 km. We note that several earthquakes are unlikely to have involved full segment rupture and may not have ruptured the full downdip width, so the stress drops in Table 1 may be underestimated for these events. We used the Wells and Coppersmith (1994) formula for rupture area vs. magnitude to estimate the fraction of the full segment area that ruptured. We consider 7 events that ruptured more than 50% of the segment area to have ruptured the full downdip width and find that these earthquakes have stress drops of 1.8-3.4 MPa, which we adopt as characteristic for the GSF. For the effects that GSF earthquakes have on adjacent GSF segments, we assume that only these seven earthquakes are significant and have neglected the contribution of smaller earthquakes. With these results as well as those of previous sections, we can model the time evolution of CFS changes on the GSF from 1750 to 2050 (see Fig. 7).
As shown in Table 3, the estimated increase in tectonic stress rates on the GSF with the case (1) of no sliver movement produces a tectonic stress rate of 0.5-2.0 kPa/yr. This low rate of tectonic stress increase is similar to that obtained by Shikakura et al. (2014) and Mitogawa and Nishimura (2020) for SW Japan's Nankai Trough forearc. If we consider 1.8-3.4 MPa as the threshold level for earthquakes to occur, this would imply a recurrence interval of 1000-7000 years, not dissimilar to the 1000-3000 suggest by paleoseismic studies of SW Japan's MTL (Tatsumi and Okada 1996). On the other hand, the characteristic slip of the 7 full-width earthquake ruptures from Table is 1.5-3.3 m, and the recently revised  slip rate for the GSF is about 15 mm/yr. Estimating a recurrence interval using a slip-predictable model (Shimazaki and Nakata 1980) would suggest a much shorter recurrence time of 100-220 years, about an order of magnitude shorter than suggested by the low tectonic stress rate of case (1), with no sliver movement.
A low tectonic stress rate implies a recurrence time of 1000 years or more for GSF earthquakes, much longer than the recurrence time of 130-450 years of major Sumatra megathrust earthquakes (Malik et al. 2019;Philibosian et al. 2014;Rubin et al. 2017). Given that the ΔCFS m such an earthquake imparts to the GSF can be as much as 1 MPa, and the typical stress drops of large GSF earthquakes can be 1.8-3.4 MPa, large Sumatra megathrust earthquakes should be able to bring forward fault failure times on the GSF by hundreds of years. For example, the Musi and Manna segments experience tectonic stress rates of 0.6 and 0.8 kPa/year, and the 1833 megathrust earthquake resulted in ΔCFS m of 0.4 and 0.2 MPa, respectively. If the GSF earthquakes in this segment followed a stresspredictable model, their rupture times would have been brought forward by 650 years for Fig. 7 Time evolution of Coulomb stress change on each GSF segment, for the case (a) without arc-parallel movement of the forearc sliver, and (b) with arc-parallel sliver movement. Coseismic changes in CFS due to the occurrence of megathrust (red) earthquakes, tectonic stress rate (green) and GSF earthquakes (magenta for same segment, blue for adjacent segments) are indicated (all from Table 3). The absolute value of the range of change in CFS from these sources, accumulated over the years 1750-2050 for each segment is indicated on the left of (a) and the right of (b). Note, however, that the poor availability of historical records in the early to mid-nineteenth century means that some segment ruptures may have been missed, possibly resulting in over-estimate of the accumulated ΔCFS Musi and 250 years for Manna. Since earthquakes did in fact occur only 60-70 years after the 1833 event, it would be logical to conclude that their occurrence was significantly influenced by the megathrust earthquake. Furthermore, while the low tectonic stress rate might lead one to expect it would take over 1000 years for stress to accumulate to the level required for earthquake recurrence, the 2007 Sumatra megathrust earthquake has moved these forward so that the Musi and Manna segments could rupture in a large earthquake even now.
On the other hand, despite the fact that the 2004 Sumatra-Andaman earthquake was widely regarded as having the potential to trigger a large GSF earthquake on the Seulimeum segment (McClosky et al. 2005), Fig. 7a suggests that, although the 2004/2005 earthquakes may have brought forward rupture there by 375 years, the ΔCFS m was not large enough to bring the CFS to the same level at which the segment ruptured in 1964. Case (1) nevertheless suggests that large megathrust earthquakes on the Sumatra fault can bring forward rupture of GSF earthquakes by hundreds of years, and in that sense, they could have a strong influence on the time-dependent seismic hazard in the Sumatra forearc.

Time evolution of ΔCFS on the GSF with forearc sliver movement
For case (2), including forearc sliver movement (Fig. 7b), the tectonic stress rate is 14-15 times greater than in case (1). This rate is much higher than in SW Japan, and this is almost certainly the reason that the GSF has a much higher level of earthquake activity than the MTL. In fact, for case (2) the tectonic stress rate is so high that the ΔCFSs due to both megathrust earthquakes and GSF earthquakes in adjacent segments are much less significant. The time evolution of the Coulomb stress on any segment of the GSF is essentially controlled by the high tectonic stress rate and occurrence of large earthquakes on that segment. If we consider that the range of stress drops for full-width rupture of past GSF earthquakes (Table 1) is 1.8-3.4 MPa, and that the average tectonic stress rate is 17 kPa/year, it takes 100-200 years for the CFS to recover its pre-earthquake level following a major event. A similar conclusion could be drawn from a consideration of the characteristics slip of 1.8-3.3 m and the ≈ 15-16 mm/yr slip rate determined from geodetic and geologic data Ito et al. 2012;Natawidjaja et al. 2017;Natawidjaja 2018a, b).
For this case, all segments on the northern part of the GSF have experienced large increases in CFS throughout the interval studied, and segments that have had large earthquakes in the past-the Seulimeum, Tripa North and South, Renun and Angkola segments-have already seen CFS recover to what it was prior to the respective earthquake. In the southern part of the GSF, the Sumani and Suliti segments appear to have not yet recovered from the previous large GSF earthquake. The Siulak segment on the other hand, despite experiencing decreases in CFS following the occurrence of a small GSF earthquake in 1995 and the 2007 megathrust earthquakes, has almost recovered from the large decrease in CFS it experienced during a large (Ms 7.3) earthquake in 1909. The Ketaun segment is similarly close to having recovered from a major earthquake it experienced in 1943, while the Kumering segment has yet to recover from the CFS decrease it experienced during a major earthquake in 1933.
Of considerable concern in Fig. 7b are the GSF segments that have never experienced a major earthquake since the high tectonic stress rate results in an increase in CFS on these segments of 5-6 MPa. These include the Aceh, Toru, Barumun, Sumpur, Sianok, and Semangko segments. First, it should be noted that the history of large GSF earthquakes is complete only since 1892 (Hurukawa et al. 2014), so any of these segments could have experienced a large earthquake at almost any time during the nineteenth century that could have substantially reduced its accumulated CFS over the 1750-2050 period. Second, some of these segments may be creeping. GPS (Ito et al. 2012) and Interferometric Synthetic Aperture radar (Tong et al. 2018) have confirmed that the Aceh segment is creeping, but geodetic studies on other segments have resolution that is too poor to establish whether they are creeping or not.
Whether we consider forearc sliver movement or not, the Musi and Manna segments appear to have accumulated enough stress that CFS has recovered to the level it was at prior to large earthquakes on these segments.

Conclusion
We have estimated and investigated the effect of the historical accumulation of Coulomb stress on the GSF from the late sixteenth to early twenty-first centuries, spanning the occurrence of eight great Sumatra megathrust earthquakes (1797 to 2010 earthquakes). This involved the estimation of the tectonic stress rate on the GSF, for which we have considered two cases: (1) no arc-parallel movement of the forearc sliver, to better compare with buttressed forearcs in other obliquely convergent margins, and (2) including the relatively fast (15-17 mm/yr) arc-parallel movement of the forearc sliver that is consistent with recent GPS studies of the Sumatra forearc .
For the case with no arc-parallel movement of the forearc sliver, the tectonic stress rate on the GSF (denoted ΔCFS 1 in Table 3) is determined solely by the coupling of obliquely convergent plate motion at the megathrust. In this case, the tectonic stress rate is very low, with all but the Tripa North segment experiencing a rate of between 0.3 and 2.0 kPa (the stress rate on the Tripa North segment is negative for this case, probably due to accommodation of thrust movement that was not considered here). This suggests a recurrence interval for full segment rupture of 1000 years or greater. Both the low tectonic stress rate and long earthquake recurrence interval are similar to those obtained for SW Japan's Nankai Trough (Hori and Oike 1996;Mitogawa and Nishimura 2020;Shikakura et al. 2014) and are likely typical of other forearcs whose arc-parallel motion is impeded by a buttressed leading edge (Beck et al. 1993). This case appears to have an order of magnitude less earthquake activity than is observed on the GSF, and the forearc deformation does not match GPS measurements of crustal velocity in the northern and southern parts of the Sumatra forearc sliver. In this case, coseismic changes in CFS can advance rupture by hundreds of years, so in that sense modulation of GSF earthquakes by megathrust events is important.
For the case including 15-17 mm/yr of arc-parallel sliver movement, which Bradley et al. (2017) show match the GPS data well throughout the length of the Sumatra forearc, the tectonic stress rate (ΔCFS 2 in Table 3) is 9 to 58 times higher than when there is no sliver movement. The recurrence interval for full-segment rupture is 100-200 years, roughly consistent with a slip-predicable model in which such ruptures experience 1.5-3.3 m slip (Table 1) at the observed slip rate of 15-16 mm/yr . In this case, changes in CFS imparted by megathrust earthquakes advance rupture by only a few years, and in this sense, they are not very important, though they may still modulate GSF event occurrence to some extent. In summary, seismic hazard along the great Sumatra fault is driven much more by the rapid arc-parallel forearc sliver motion than by megathrust earthquake activity. The Sumatra forearc may be unique in this respect, since instead of being buttressed at its leading edge, arc-parallel motion of the sliver appears to 1 3 be accommodated by seafloor spreading in the Andaman Sea to the north of Sumatra (Curray 2005).
A concern highlighted in this study is the potentially high stress accumulation on the GSF segments that have never experienced full segment rupture, like the Aceh, Toru, Barumun, Sumpur, Sianok, and Semangko segments. Although the Aceh segment is creeping (Tong et al 2018;Ito et al. 2012) and therefore unlikely to accumulate the high level of stress indicated in Fig. 7b, it is unknown whether the other segments are creeping or not. It is entirely possible that some of these segments may have ruptured in the early to mid-nineteenth century, when historical records are sparse, which would reduce the accumulated stress. A new compilation of historical data by Martin et al. (2018) suggests some of these segments are likely to have ruptured early in the nineteenth century, although the magnitudes are difficult to constrain. In any case, further study of creep on these segments would help constrain seismic hazard along the GSF. On the other hand, some segments that have experienced full-segment rupture and are therefore likely locked, such as the Tripa (North and South), Angkola, Musi and Manna segments, appear to be well advanced in their seismic cycles and could experience another large event soon.
Finally, we point out that our model of the GSF is simplified and likely does not account for much of the complexity on the actual GSF. Segment geometries have been simplified and are all assumed to be pure strike-slip, when in reality some segments likely involve some component of oblique motion. Although Bradley et al. (2017) find that the available GPS data are best fit with no internal deformation of the Sumatra forearc sliver, there may be some deformation that is simply not resolved by the relatively sparse station spacing, and such internal deformation would likely reduce the tectonic stress rates calculated in Table 3. Also, accurate calculation of the time evolution of ΔCFS along the GSF should account for viscoelastic, postseismic motion, which we have not considered here. However, Qui and Chan (2019) have shown that these can change the ΔCFS on the GSF by about 10%, so this is unlikely to alter our conclusions. In any case, we believe the model we have developed in this study for stress accumulation and release along the GSF will prove useful in guiding further seismic hazard studies in Sumatra.

Conflicts of interest
The authors declare that the research was conducted in the absence of any commercial and financial relationships that could be construed as a potential conflict of interest. The authors have no relevant financial or non-financial interests to disclose.

Human and animal rights
No animal studies are presented in this manuscript. No human studies are presented in this manuscript.
Informed consent There are no human subjects in this article, and informed consent is not applicable.