GPS-based slip models of one M_w 7.2 and twenty moderate earthquakes along the Sumatran plate boundary

Abstract

Earthquake-induced deformation along the Sumatran plate boundary has been monitored by the Sumatran GPS Array (SuGAr) since 2002. This continuous GPS network recorded the coseismic deformation of 10 earthquakes with moment magnitude (M_w) larger than 7 and 20 with M_w in the range of 5.9–7 from 2002 to 2013. Among all these recorded events, one large M_w 7.2 event and most of the moderate ones (5.9 ≤ M_w < 7) have yet to be modeled with available GPS data. This is partially due to the limited number (≤ 4) of stations that recorded each event. In this paper, we explore the possibility of using the limited observations to derive sensible slip models for these “forgotten” Sumatran events. We model each event as a single rectangular patch of uniform slip and constrain most of the patch parameters using external information based on slab geometry and global teleseismic catalogs. For each event, we use a grid-search approach to find the preferred location of slip patches, which we present along with contours of error-weighted variance explained to indicate the uncertainties. We compare the center locations of our final slip patches with the centroid locations from the global Centroid Moment Tensor (gCMT) catalog and the epicenter locations from four other global catalogs. Our results show that the gCMT centroid locations for the 21 Sumatran earthquakes are systematically biased toward the southwest relative to the centers of our slip patches, while the epicenter locations from the four other catalogs are all consistently shifted toward the northeast. Although the available data have no resolving power for other source parameters, we find that simple forward modeling based on sparse but reliable near-field GPS data generally provides less biased and more accurate locations than global teleseismic catalogs along the Sumatran plate boundary. The catalog of slip models we present will have particular utility in the event of other significant earthquakes being generated by the same or proximal areas of the Sunda megathrust.

Introduction

The Sumatran plate boundary has experienced a surge of seismic activity in the years since the 26 December 2004 M_w 9.2 Sumatra–Andaman earthquake and is currently one of the most seismically active convergent plate boundaries in the world (Feng et al. 2015). The deformation caused by seismic events has been monitored continuously by the Sumatran GPS Array (SuGAr), which was first established in 2002. Not only did the SuGAr record large events (M_w ≥ 7), but it also captured many moderate (5.9 ≤ M_w < 7) events that occurred close enough to the SuGAr stations. Although most of the large earthquakes and two moderate events (the 10 April 2005 and 16 August 2009 M_w 6.7 events (Wiseman et al. 2011; Wang et al. 2018)) have been modeled using the available SuGAr data, one M_w 7.2 event and the majority of the moderate events remained unmodeled geodetically. This paper therefore aims to model these “forgotten” events with the available SuGAr data and presents a catalog of coseismic slip models for the M_w 7.2 event and 20 recorded moderate ones (5.9 ≤ M_w < 7) that include 17 thrust and 3 strike-slip events (Fig. 1).

Fig. 1

Map of all the events modeled in this study. Rectangular colored patches or short lines show the final slip models for one M_w 7.2 event and 20 moderate (5.9 ≤ M_w < 7) events obtained from this study. Focal mechanisms are placed at the gCMT centroid locations (Dziewonski et al. 1981; Ekström et al. 2012). Yellow stars indicate the Advanced National Seismic System (ANSS) epicenters (NCEDC 2016). Solid colored lines outline published 1-m slip contours for the 2004 M_w 9.2 Sumatra–Andaman earthquake (Chlieh et al. 2007), the 2005 M_w 8.6 Nias–Simeulue earthquake (Konca et al. 2007), the 2007 M_w 8.4 Bengkulu earthquake and its two M_w 7.9 and 7.0 aftershocks (Tsang et al. 2016), the 2008 M_w 7.2 North Pagai earthquake (Salman et al. 2017), the 2010 M_w 7.8 Mentawai earthquake (Hill et al. 2012), the 2002 M_w 7.3 and 2008 M_w 7.4 Simeulue sibling earthquakes (Morgan et al. 2017), and the 2010 M_w 7.8 Banyaks Island earthquake (Morgan et al. 2015). Red circles are the SuGAr stations (Feng et al. 2015). Black dashed boxes delineate the inferred subducted fracture zones underneath Simeulue (Franke et al. 2008) and the Batu Islands (Pesicek et al. 2010) based on seismicity. Brown lines are slab contours at 20 km, 40 km, and 60 km intervals from Slab1.0 (Hayes et al. 2012)

The source parameters of these events are available from global teleseismic catalogs, but their solutions might be poorly constrained or biased. Our near-field GPS observations, though the number is limited, might provide some extra independent information on source parameters. We show in the rest of the paper that sparse near-field GPS data can be useful for determining and sometimes improving the location of moderate earthquakes. The improved locations will potentially be useful for providing a more complete and accurate slip history of the Sunda megathrust, particularly in the event of significant future earthquakes occurring in the vicinity.

Data and methods

Data

We obtained both the horizontal and vertical static coseismic offsets for the M_w 7.2 event and 20 moderate ones from Feng et al. (2015). Feng et al. (2015) derived these coseismic offsets, which are usually small in magnitude, by simultaneously estimating many other and larger signals that dominate the SuGAr daily position time series. These other signals include (1) long-term rates, (2) annual and semiannual signals, and (3) coseismic offsets and postseismic decays for the great and large earthquakes that are not modeled in this study.

Methods

Given the relatively small magnitude (5.9 ≤ M_w < 7) of the moderate earthquakes, each event was recorded only by a few (up to 4) SuGAr GPS stations. As the data from these stations were too sparse for an inversion, we conducted iterative grid-search forward modeling instead. We modeled each earthquake as a single rectangular patch of uniform slip using the Okada dislocation model (Okada 1985, 1992). To limit the number of unknown parameters, we imposed constraints on the length, width, strike, and dip of the slip patch while exploring the location (longitude, latitude, and sometimes depth) of the slip patch. We additionally varied the rake to test the sensitivity of the location results.

The dimensions of the slip patch can be constrained using empirical scaling relations that provide characteristic rupture lengths and widths for given earthquake magnitudes (Wells and Coppersmith 1994; Strasser et al. 2010; Blaser et al. 2010). We chose to use the relations from Blaser et al. (2010) as they included all types of events with a special focus on the subduction zone environment, so as to provide a consistent method to constrain the rupture dimensions. Once the rupture dimensions were constrained, we calculated the required slip for a given magnitude by assuming a rigidity of 30 × 10⁹ Pa.

We constrained the strike, dip, and depth of the slip patch based on the type of event, i.e., whether it was a megathrust event, a thrust event not on the megathrust (non-megathrust thrust event), or a strike-slip event. It is sometimes difficult to tell whether a thrust event was on the megathrust or not. Combining evidence from GPS horizontal displacements, focal plane solutions, focal depths, and whether the event occurred within narrow seismic bands that have previously been suggested to occur on the megathrust, Feng et al. (2015) suggested that 13 out of the 17 moderate thrust events were likely to be megathrust events. We thus modeled these 13 events as megathrust events. Feng et al. (2015) also suggested that two other events (the 15 April 2009 M_w 6.3 and 4 January 2008 M_w 6.0 events) could be either on the megathrust or within the overriding plate. As the data for these two events were too sparse to distinguish between the two scenarios, we still modeled them as megathrust events but noting that they could be shifted shallower. In contrast to the seaward horizontal displacements generated by the abovementioned 15 events, the remaining two thrust events (the 10 April 2005 M_w 6.7 and 16 August 2009 M_w 6.7 events) generated large landward horizontal motions and are therefore thought to be thrust events not on the megathrust.

For megathrust events, we constrained the strike, dip, and burial depth (top of the slip patch) to the interpolated geometry obtained from the Slab1.0 subduction interface model (Hayes et al. 2012). For other events, the fault geometry is generally unknown. As such, we varied burial depths every 5 km between 0 and 25 km to find the best-fit depth for these events when constraining the strike and dip to the global Centroid Moment Tensor (gCMT) solution (Dziewonski et al. 1981; Ekström et al. 2012). The gCMT catalog is one of the most commonly used products for moment tensor solutions.

For all events, we varied rake ± 90° in steps of 0.1° from the initial rake of the event type to find the best-fit rake. The initial rake is 90° for thrust events, 0° for left-lateral strike-slip events, and 180° for right-lateral strike-slip events.

For each rake and each depth (when applicable), we conducted a grid search over the longitude and latitude of the top-left corner of the slip patch. In the first coarse-step (0.05°) search, we varied the longitude and latitude in a 2 × 2° area centered at the gCMT location. In the second fine-step (0.01°) search, we searched within a smaller area of 1° × 1° centered at the best-fit location from the previous coarse step.

As a result of the grid search, we iteratively ran millions of forward models for each rake and depth (when applicable) of each event. We quantified the goodness-of-fit by calculating the percentage of the error-weighted data variance that can be explained by each forward model, which we name the error-weighted variance explained (ve) and its equation is given as follows:

$${\text{ve}} = \left( {1 - \frac{{\sum\nolimits_{i} {\left( {\frac{{d_{i} - m_{i} }}{{e_{i} }}} \right)^{2} } }}{{\sum\nolimits_{i} {\frac{{d_{i}^{2} }}{{e_{i}^{2} }}} }}} \right) \times 100\% ,$$

where d represents the observed coseismic offsets (both horizontal and vertical components) recorded by the GPS stations, e the errors in these observations, and m the coseismic offsets predicted by forward models. The reason why we used ve instead of the traditional root-mean-square misfit is because ve does not depend on the magnitude of offsets as misfit does, thus allowing the direct comparison between events for consistent uncertainty assessment.

For megathrust events, the contours of error-weighted variance explained for each set of grid search often have two or three local maxima that can be chosen as the preferred model. For most events, the global maximum was closer than other maxima to the locations given by global teleseismic catalogs, in which case the global maximum was chosen as the preferred model (Fig. 2a). However, for two events (the 5 July 2005 M_w 6.6 Nias and 9 May 2010 M_w 7.2 Simeulue events), the global maximum was not located on the lobe closest to the teleseismic locations in which case the local maxima closest to the teleseismic locations was then chosen as the preferred model (e.g., Fig. 2b). For other events, the contours usually show only one maximum, which is simply the preferred model.

Fig. 2

Grid-search results for (a) the 26 February 2005 M_w 6.7 Simeulue event and (b) the 5 July 2005 M_w 6.6 Nias event with contours showing the error-weighted variance explained (ve) of numerous forward models. Black boxes outline the surface projection of the preferred rupture plane. Red circles are the SuGAr stations that had been installed before the event, while white circles are those installed after the event or decommissioned before the event. Green and red vectors represent the observed vertical and horizontal displacements, while black vectors represent the displacements predicted from our preferred model. Yellow stars, orange triangles, blue diamonds, and red inverted triangles represent the epicenters from the ANSS, ISC, ISC-EHB, and ISC-GEM teleseismic catalogs, respectively. Focal mechanisms are placed at the gCMT centroid locations (Dziewonski et al. 1981; Ekström et al. 2012). Both events have three local maxima in their contours. The global maximum for the Simeulue event is the closest to the teleseismic locations; thus, the global maximum is chosen as the preferred model. However, the global maximum near LHWA for the Nias event is only the second closest to the teleseismic locations; thus, the local maximum closest to the teleseismic locations is chosen as the preferred model. Brown lines are slab contours at 20 km, 40 km, and 60 km intervals from Slab1.0 (Hayes et al. 2012)

We obtained one preferred model for each rake and each depth (when applicable). The final model is the preferred model when we use the best-fit rake and best-fit depth (when applicable).

Geodetic results and comparisons with other datasets

Among the 21 events that we model in this paper, 11 were recorded by only one GPS station and the other 10 were recorded by no more than four stations that are often located on the same side (either trenchward or landward) of the event. Given the data limitations, we ask the question of whether one station or a sparse and one-sided group of stations are sufficient for determining the location of an earthquake. To answer this question, we first conduct additional tests for megathrust events to quantify the change in location due to the change in certain model parameters; we then compare the center locations of our final slip patches (Table 1) with geologic field observations, global teleseismic catalogs, and local seismic catalogs.

Table 1 Summary of the preferred model for all events modeled in this study

Sensitivity tests

To test how the change in rake affects the location of megathrust events, we alternatively fix the rake at the gCMT solution and compare the results with those from the best-fit rake models. The best-fit rakes differ from their corresponding gCMT rakes in a wide range (Fig. 3a), which suggests that rake is not well-constrained. Even though we use the best-fit rake in our final models, the location difference is not significant when switching from the best-fit rake to the gCMT rake. Constraining the rake to the gCMT solution causes an average absolute location shift of 3.4 ± 3.3 km along the east–west direction and 4.0 ± 3.7 km along the north–south direction (Fig. 3c). These values are relatively small compared to the dimensions of a moderate earthquake (11 km × 8 km for M_w 6.0; 42 km × 23 km for M_w 7.0), implying that the location is not sensitive to rake. In addition, we observe no systematic bias toward any direction as shown by the near-zero average location shift, though the location shift tends to occur along the northwest–southeast direction parallel to the trench (Fig. 3b).

Fig. 3

Sensitivity tests for all modeled megathrust events. Changes in a rake, b location, and c absolute location, when rake is constrained to the initial solution as opposed to the best-fit value. Changes in d rake, e location, and f absolute location, when strike and dip are constrained to the gCMT solutions as opposed to the interpolated Slab1.0 interface

To test how the change in strike and dip affects the location of megathrust events, we alternatively constrain the strike and dip to the gCMT solution and compare the results with those from the Slab1.0 models. The average absolute east–west and north–south shifts are 1.0 ± 0.9 km and 2.1 ± 3.0 km, respectively (Fig. 3f). The average east–west and north–south shifts are 0.0 ± 1.4 km and 0.2 ± 3.8 km, respectively (Fig. 3e). These small values suggest no significant systematic shift in location when the strike and dip are constrained to either the gCMT solutions or the interpolated Slab1.0 interface. Therefore, small changes in the strike and dip do not seem to significantly change the center locations of our final slip patches.

Comparison with geologic field observations

The vast majority of the modeled events occurred under the ocean; thus, it is difficult to find direct evidence for the exact location of these events. The only exception is the 6 March 2007 M_w 6.4 earthquake and another M_w 6.3 event 2 h later, for which twin-surface ruptures were documented along the Sumatran fault (Daryono et al. 2012). They found unequivocal evidence of tectonic fault rupture at 10 localities (tan squares in Fig. 4) along the Sianok and Sumani segments that straddle a prominent releasing stepover, providing the most accurate locations that we can compare with our results.

Fig. 4

Grid-search results for the 6 March 2007 M_w 6.4 and 6.3 Sumatran fault earthquake doublet when the burial depth is fixed at surface. Tan squares show 10 localities where twin-surface ruptures associated with this event were documented (Daryono et al. 2012). The final patch is located over the pull-apart basin between the Sianok (dark brown) and Sumani (dark blue) segments to the northwest and southeast, respectively. See Fig. 2 caption for description of other features

Because the two events occurred so close in time, we cannot separate their individual effects using the coseismic offsets derived from daily GPS positions by Feng et al. (2015). Thus, we model the total displacements caused by the doublet using one single event of moment magnitude M_w 6.56 that combines the energy released from both events, and constrain the strike and dip of the plane to the gCMT solution of the first larger shock. Searching from 0 to 25 km every 5 km, we find the best-fit burial depth at 0 km (Fig. 4), which is consistent with the fact that the twin ruptures reached the surface. The horizontal location of our final model is located over the pull-apart basin between the Sumani and Sianok segments, consistent with the general locations of the mapped surface ruptures. The contours of error-weighted variance explained show two local maxima with one concentrated along the Sumatran fault and the other located approximately between the two GPS stations (Fig. 4).

Comparison with global teleseismic catalogs

We also compare the centers of our final slip patches with locations from five widely used global teleseismic bulletins or catalogs. Those include the global Centroid Moment Tensor (gCMT) catalog (Dziewonski et al. 1981; Ekström et al. 2012), the Advanced National Seismic System (ANSS) composite catalog (NCEDC, 2016), the Bulletin of the International Seismological Centre (ISC) (International Seismological Center, 2016), the ISC-Engdahl–Van der Hilst–Buland (ISC-EHB) Bulletin (Engdahl et al. 1998; Weston et al. 2018), and the ISC-Global Instrumental Earthquake (ISC-GEM) catalog (Storchak et al. 2013, 2015). The three ISC products suit different research needs: the ISC Bulletin is a comprehensive global summary of natural and anthropogenic events; the ISC-EHB Bulletin is a relocated catalog of selected events from the ISC Bulletin with the focus on using teleseismic depth phases (pP, pwP, and sP) to improve depth estimation; and the ISC-GEM is an extensive list of moderate to large global earthquakes selected from the ISC Bulletin with the focus on homogeneous estimates for location and magnitude (Weston et al. 2018).

The gCMT location is a centroid that refers to the center of the seismic moment release in space (Ekström et al. 2012), while the ANSS, ISC, ISC-EHB, and ISC-GEM locations are epicenters where the rupture nucleates. These two types of location do not necessarily coincide with each other or with the center of a slip patch, especially for large earthquakes. But as the rupture areas of moderate earthquakes are relatively small (~ 100 km² for a M_w 6.0; ~ 1000 km² for a M_w 7.0), we assume they can be treated as point sources in teleseismic distances, and thus we expect patch center, centroid, and epicenter to be located relatively close to each other.

We summarize the average location shifts for the five global catalogs relative to our patch centers in Fig. 5. The gCMT centroids show an overall systematic shift toward the southwest with respect to our patch centers (Fig. 5a), while the ANSS, ISC, ISC-EHB, and ISC-GEM epicenters show a systematic shift toward the northeast relative to our patch centers (Fig. 5b–e). This general pattern is not true for every event, but it is valid for the bulk of the events studied. The average north–south shifts for all the five global catalogs are on the similar order of 10 km, while the average east–west shift for ISC-EHB is the smallest (2.8 ± 11.6 km) with less than half of the values for the other four catalogs.

Fig. 5

Location shifts for a gCMT centroids, b ANSS epicenters, c ISC epicenters, d ISC-EHB epicenters, and e ISC-GEM epicenters relative to the centers of our final models. $\bar{E}$ and $\bar{N}$ show average shifts and their standard deviations in the east–west and north–south directions, respectively. The gCMT centroids are biased toward the southwest with respect to our patch centers, while the ANSS, ISC, ISC-EHB, and ISC-GEM epicenters are biased toward the northeast with respect to our patch centers

Although the five global catalogs use different seismic phases and periods, they all rely heavily on the Global Seismographic Network (GSN) that provides a worldwide monitoring of global seismicity with over 150 modern seismic stations. But the station distribution of the GSN is not uniform. The Sumatran subduction zone happens to reside in a region where the GSN station coverage is less dense than other seismically active subduction zones such as Japan (Ammon et al. 2010). In addition, the Sumatran subduction zone faces the Indian Ocean to the west, where on land instrumentation is not possible, so most stations lie to the east resulting in an unbalanced station geometry (Engdahl et al. 2007). All the five catalogs use a 1-D Earth velocity structure that does not account for lateral variations in seismic velocities. This unmodeled 3-D bias is extremely large, even with good station coverage, in subduction zones where high-velocity subducting slabs exist (Bondár et al. 2004). Therefore, the Sumatran subduction zone suffers from both the unfavorable station coverage and unmodeled 3-D velocity structure, which we suggest are the causes for the observed average regional bias in location reported in global teleseismic catalogs relative to our GPS-based locations. We also suggest that the smaller east–west bias in the ISC-EHB epicenter determination is possibly due to the inclusion of depth phases and the stricter selection of seismic stations in the EHB algorithm (Weston et al. 2018).

Comparison with local seismic catalogs

While someone could argue that the regional bias is due to the sparsity and one-sided geometry of our GPS data, we show that our GPS-based locations are overall less biased than the five global teleseismic catalogs by comparing the centers of our final slip patches with locations from three local seismic catalogs. The three local catalogs were all derived from a dense temporary local seismic network (Tilmann et al. 2010; Lange et al. 2010; Collings et al. 2012); thus, they should provide more accurate epicenters than global teleseismic catalogs (Bondár et al. 2004). Among the three temporary networks, two comprising both land stations and ocean bottom seismometers were deployed in the Simeulue segment between October 2005 and March 2006 (Tilmann et al. 2010), and the Nias-Batu segments between April 2008 and February 2009 (Lange et al. 2010); the other comprising only land stations was deployed in the Mentawai segment between December 2007 and October 2008 (Collings et al. 2012).

One pronounced cluster of seismicity in the Simeulue and Nias-Batu catalogs is a continuous coast-parallel seismic band situated immediately seaward of Simeulue and Nias (Figs. 6 and 7). The band of seismicity became active mainly after the 2005 M_w 8.6 Nias-Simeulue earthquake (Pesicek et al. 2010) and formed the vast majority of the 2005 aftershock sequence (Tilmann et al. 2010). This band has a narrow width of 20–30 km, with most seismicity occurring on the megathrust (Tilmann et al. 2010; Lange et al. 2010). Although none of our events were recorded by the two local catalogs, several common events were found in the Simeulue catalog and global catalogs (Tilmann et al. 2010). For these common events, Tilmann et al. (2010) found a significant seaward bias in the gCMT locations and a lesser degree landward bias in the EHB locations relative to their locally determined locations.

Fig. 6

Map of seven modeled earthquakes in the Simeulue section. Rectangular colored patches represent our final slip models, and the models elaborated upon in the main text are highlighted with a gray border. Yellow stars, orange triangles, blue diamonds, and red inverted triangles represent the epicenters from the ANSS, ISC, ISC-EHB, and ISC-GEM teleseismic catalogs, respectively. Black dots indicate seismicity between October 2005 and March 2006 located by Tilmann et al. (2010). See Fig. 1 caption for description of other features

Fig. 7

Map of seven modeled earthquakes in the Nias section. Black dots indicate seismicity between October 2005 and March 2006 located by Tilmann et al. (2010). Gray dots indicate seismicity between April 2008 and February 2009 located by Lange et al. (2010). See the captions of Figs. 1 and 6 for description of other features

This bias is consistent with the regional bias we find in the gCMT and EHB locations relative to our locations (“Comparison with global teleseismic catalogs” section). This regional bias is best illustrated by the 10 thrust events (as shown in Fig. 16b–k of Feng et al. 2015) that occurred in the vicinity of the narrow seismic band. While their locations from different global catalogs are scattered, it is reasonable to assume that their actual locations fall within the seismic band. Compared to the location of the seismic band, the gCMT centroids of the 10 events tend to shift toward the southwest (seaward), while their epicenters tend to shift northeastward (landward) (Figs. 6 and 7). On the contrary, our slip patches are located mostly within this seismic band (Figs. 6 and 7) with only one exception (the 25 July 2012 M_w 6.4 Simeulue event, Fig. 6), indicating that our locations are generally more consistent and less biased than those reported in the global catalogs.

We find only one event (the 4 January 2008 M_w 6.0 event, Fig. 8) recorded in both our catalog and the Mentawai catalog (Collings et al. 2012). Similar to the pattern observed for the 10 thrust events in the narrow seismic band, the gCMT centroid and all the epicenters are biased toward the southwest and the northeast, respectively, while our slip patch is very close to the locally determined epicenter (Fig. 8). It is worth noting that our model with depth fixed at the slab interface significantly underestimated the coseismic offsets (Additional file 1: Section S1.9). As the locally determined depth (27 km) (Collings et al. 2012) is shallower than the slab depth (~ 40 km), it is possible that the event occurred shallower in the overriding plate. Even though the depth of our model might be misplaced and the magnitude is significantly underestimated, the horizontal location seems to be still relatively well-constrained with only one station. Examining all the thrust events recorded by only one station, we find that their final models are all aligned along the direction of their horizontal coseismic offsets (see Additional file 1). Thus, we suggest that reliable coseismic offsets, particularly reliable coseismic direction estimates, play a key role in good location determination for moderate events recorded by sparse GPS data.

Fig. 8

Map of six modeled earthquakes in the Mentawai section. Black dots indicate seismicity between December 2007 and October 2008 located by Collings et al. (2012). Gray star indicates the epicenter of the 4 January 2008 event from Collings et al. (2012). Backthrust fault traces are taken from (Singh et al. 2010). See the captions of Figs. 1 and 6 for description of other features

Discussion

While we present the whole catalog of our final slip models for the M_w 7.2 and 20 moderate events as tables and figures in Table 1, and Additional files 1 and 2, in this section we discuss selected events in the Simeulue, Nias, and Mentawai sections, respectively, along the Sumatran subduction zone. This allows us to put our results in the broad context of regional seismicity.

Simeulue

The Simeulue section is known as a persistent rupture barrier against which several great earthquakes historically terminated from both the north and the south (Meltzner et al. 2012). It is also a rupture generator within which moderate to large earthquakes occurred (Morgan et al. 2017). Six of our moderate earthquakes occurred seaward of Simeulue (Fig. 6). Our results suggest five of them are typical megathrust events that occurred within the well-known narrow seismic band. The remaining event (the 25 July 2012 M_w 6.4 event) is located to be seaward of the seismic band. If our location is accurate, this event may present seismicity along the inferred subducted fracture zone (Franke et al. 2008). However, according to the contours of the error-weighted variance explained, it is possible to move the location landward so that this event may also fall into the seismic band (Additional file 1: Section S1.16).

The largest earthquake we model in this paper is the 9 May 2010 M_w 7.2 event north of Simeulue (Fig. 6). The magnitude of this event is large, but the available GPS data are too sparse for a rigorous study. Although our final location is much closer to the gCMT centroid at a shallower depth than the various epicenters, the location could be at a deeper depth closer to the epicenters (Additional file 1: Section S1.13). In any case, this event occurred near the southern downdip edge of the 2004 M_w 9.2 Sumatra–Andaman earthquake (Chlieh et al. 2007) and perhaps terminated against the inferred subducted fracture zone (Fig. 6).

Nias

In the Nias section, our results suggest four thrust events to be typical megathrust events that occurred within the seismic band seaward of Nias (Fig. 7). This seismic band is a continuation of the seismic band seaward of Simeulue.

The near-trench area in the Nias region is often characterized by extremely low seismicity (Tilmann et al. 2010); however, the 16 May 2006 M_w 6.8 event seems to occur in a locus of microseismicity (Lange et al. 2010) (Fig. 7). Based on our modeling results, we find that a left-lateral strike-slip earthquake on an N–S trending fault plane fits the data much better than a right-lateral strike-slip event on an E–W trending fault plane. Furthermore, the best-fit depth of the N–S trending plane is ~ 20 km, suggesting that this event most likely occurred within the subducting plate. Like some other strike-slip events within the Indian and Australian plates in the Wharton Basin (Robinson et al. 2001; Lay et al. 2016), this event is probably associated with the reactivation of an existing N–S trending fracture zone (Deplus et al. 1998).

Mentawai

In the Mentawai section, the 10 April 2005 and 16 August 2009 M_w 6.7 events are of particular interest. These events were initially suggested to occur on the shallowly seaward-dipping Mentawai backthrust faults (Wiseman et al. 2011); however, recent broadband waveform modeling refined the locations of these two sequences and suggested that they occurred on steeply (~ 60°) landward-dipping backstop faults (Wang et al. 2018). From our modeling, we find that the backthrust and backstop models fit the data equally well, so our GPS data are insufficient to distinguish the two fault planes. Given the recent evidence that favors backstop faults, we use our final backstop models for these two events.

Conclusion

Our model results show that even with a limited number of near-field GPS stations, we can constrain the horizontal locations of moderate earthquakes relatively well when coseismic offsets and directions are reliably estimated. Compared to our final locations, we observe a southwestward (seaward) bias in the gCMT centroids and a northeastward (landward) bias in the ANSS, ISC, ISC-EHB, and ISC-GEM epicenters for events along the Sumatran plate boundary. A joint inversion of teleseismic data with near-field GPS data could potentially reduce epicenter mislocation errors in teleseismic locations for regions like Sumatra that has a significant lack of seismic stations. Our catalog presented in Additional file 1 provides information that may be useful for compiling a more complete and accurate slip history of the Sunda megathrust, and for studying future large earthquakes if they occur in the vicinity of our events. Many of the events in our catalog occurred near the boundaries of either geologic or geometric structures or major earthquakes, so they may also be useful for providing insights into the behavior of these boundaries.