Toward comparable relative locations between the mainshock slip and aftershocks via empirical approaches

The relative locations between mainshocks and their aftershocks have long been studied to characterize the mainshock–aftershock relationships, yet these comparisons may be subjected to biases inherited from various sources, such as the analysis method, data, and model parameters. Here, we perform both a relocation analysis of interplate events to obtain accurate relative centroid locations and a slip inversion analysis of the mainshock slip relative to the relocated events, with some of the relocated events used as empirical Green’s functions, to retrieve the spatiotemporal slip features of the mainshock relative to all of the relocated events. We perform these analyses on the large (6.0≤MW≤7.3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$6.0 \le M_{W} \le 7.3$$\end{document}) interplate earthquakes that occurred near four giant (MW≥8.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{W} \ge 8.5$$\end{document}) megathrust earthquakes: the 2007 MW8.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{W} \;8.5$$\end{document} Bengkulu (Indonesia), 2005 MW8.6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{W} \;8.6$$\end{document} Nias (Indonesia), 2010 MW8.8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{W} \;8.8$$\end{document} Maule (Chile), and 2011 MW9.1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{W} \;9.1$$\end{document} Tohoku-Oki (Japan) earthquakes. Most of the spatiotemporal slip features of the mainshocks are consistently recovered using different empirical Green’s functions. We qualitatively and quantitatively demonstrate that the large interplate aftershocks within 5 years of the four analyzed mainshocks are largely located on the periphery or outside of the large-slip regions of these four giant megathrust earthquakes.


Introduction
The relationship between the mainshock rupture and its subsequent aftershock sequence has been widely discussed. Many studies have associated aftershocks with the coseismic Coulomb stress changes caused by the mainshock rupture (e.g., Nakamura et al. 2016;Reasenberg and Simpson 1992;Seeber and Armbruster 2000;Stein 1999) or afterslip in the surrounding region (e.g., Hsu et al. 2006;Lange et al. 2014). Other mechanisms have also been proposed, such as the effect of fluid pressure pulses generated by the coseismic release of trapped CO 2 (Miller et al. 2004), the maximum change in shear stress, von Mises yield criterion, and the sum of the absolute value of the stress change tensor, which perform best in explaining neural-network-predicted aftershock locations (DeVries et al. 2018). If there are fewer aftershocks in the large mainshock slip region and an increase in aftershocks in the peripheral regions, then this could be indicative of a larger stress release in the large coseismic slip region and stress buildup along the periphery of this large coseismic slip region (Wetzler et al. 2018).
Numerous comparisons between the mainshock slip region and aftershock locations have been performed to characterize the mainshock-aftershock relationship. For example, Das and Henry (2003) and Mendoza and Hartzell (1988) qualitatively compared and summarized the slip locations for a number of mainshocks and their aftershock distributions, with both studies highlighting that the common feature among most of the mainshock events was the lack of aftershocks in the large-slip region. Quantitative studies have also been performed. Agurto et al. (2012) found that the M W > 4.0 interplate aftershocks of the 2010 M W 8.8 Maule Earthquake occurred mainly in regions where the coseismic slip was 20-70% of the peak slip, with these aftershock regions surrounding the patches that experienced the largest coseismic slip. Wetzler et al. (2018) conducted a slip inversion and aftershock relocation analysis for 101 M W ≥ 7.0 subduction-zone events, and found that the aftershocks occurred mainly along the periphery of the large mainshock slips. The occurrence of small numbers of aftershocks in the large-slip regions has also been reproduced by Yabe and Ide (2018) via the introduction of frictional heterogeneity in numerical simulations.
Accurate relative locations between the mainshock slip region and associated aftershock distribution are required to better understand this mainshock-aftershock relationship. Both the mainshock slip and aftershock locations can be determined more accurately with better station coverage, particularly for inland events (e.g., Ross et al. 2017;Woessner et al. 2006). However, the comparison can be challenging for subduction zone megathrust earthquakes due to the one-sided configuration of seismic networks on land, along with the limited data coverage and availability over the wide rupture area (e.g., Madariaga et al. 2010). For example, Rietbrock et al. (2012) pointed out the difficulty in discussing the relative locations between the 2010 M W 8.8 Maule Earthquake and its aftershocks due to the large discrepancy between the slip models determined by different studies. Therefore, care needs to be taken when attempting such a comparison since the different data types (e.g., seismic or geodetic), distribution of recording stations (e.g., nearfield or far-field), and applied parameters (e.g., velocity model) used to derive the slip models may contain their own biases.
Here we propose a new empirical approach to address these potential sources of bias, and analyze four giant megathrust earthquakes: the 2007 M W 8.5 Bengkulu (Indonesia), 2005 M W 8.6 Nias (Indonesia), 2010 M W 8.8 Maule (Chile), and 2011 M W 9.1 Tohoku-Oki (Japan) earthquakes. We use seismic data from the same teleseismic (30°-95°) network, and apply the same velocity structure for the travel time calculations to determine both the mainshock slip and aftershock locations under the same assumption that the Green's function from a given station to any location on the shallow (≤ 70 km) subduction zone is the same, except for the travel time difference. We first perform a centroid relocation analysis of the large ( 6.0 ≤ M W ≤ 7.3 ) interplate earthquakes around the rupture regions of the mainshocks using the teleseismic network correlation coefficient (NCC) method (Chang and Ide 2019), which is an extension of the NCC method (Ohta and Ide 2008), to relocate events with similar focal mechanisms that are distributed over a few hundred kilometers in subduction zones, assuming waveform similarity among the events. We then utilize the event waveforms of some of the larger relocated interplate events as empirical Green's functions (EGFs), again assuming the same Green's function between the mainshocks and the surrounding interplate events at a given station, for the slip inversion analysis. The reliable portions of the solutions are determined by applying different EGFs individually and comparing the results. We obtain the relocation and slip inversion results under a common coordinate, where we use the same earthquake catalog for both the initial locations in the relocation analysis and the hypocenter locations for the mainshocks in the slip inversion analysis. We lastly perform a quantitative comparison between the relocated aftershocks and slip inversion results, similar to the approaches of Agurto et al. (2012) and Woessner et al. (2006), and show a lack of large interplate aftershocks in the large-slip regions for most of the analyzed mainshocks.

Method and data
We adopt the teleseismic NCC method to relocate the large ( 6.0 ≤ M W ≤ 7.3 ) interplate earthquakes around giant megathrust events, and obtain accurate relative centroid locations. The original NCC method is outlined by Ide (2008, 2011) and has since been extended by Chang and Ide (2019) to relocate teleseismic events (30°-95° distances). The relative locations between event pairs are first determined by maximizing the sum of the correlation coefficients between the waveforms from all of the available stations and components throughout the entire network under the appropriate arrival time shifts. This is followed by an inversion that yields the final locations, with the initial catalog locations as prior information and the relative locations being weighted using the summation of the correlation coefficients Ide 2008, 2011). We account for the waveform differences due to the systematic increase in rupture duration with increasing event magnitude during the cross-correlation to determine the relative location between two events by convolving the waveforms of one event with an assumed triangle source-time function (Meier et al. 2017) that possesses the rupture duration of the other event (Chang and Ide 2019), as follows: where t d is the rupture duration, R is the radius of an assumed circular source with a 3 MPa stress drop (Kanamori and Anderson 1975;Ye et al. 2016), and V R is the rupture velocity (2.5 km/s; Chounet et al. 2018;Ye et al. 2016). (1) Here we relocate the large, shallow (6.0 ≤ M W ≤ 7.3, depth ≤ 70km) interplate earthquakes that occurred near four giant earthquakes ( M W ≥ 8.5 ; subplots in Fig. 1) between 2001 and 2016, taken from the Global Centroid Moment Tensor (GCMT) catalog (Dziewoński et al. 1981;Ekström et al. 2012). We define interplate earthquakes as those where the Kagan angle (Kagan 1991(Kagan , 2007) between a given earthquake and its neighboring giant megathrust earthquake is ≤ 30 • , which is determined using the focal mechanisms provided by the GCMT catalog (Dziewoński et al. 1981;Ekström et al. 2012). A complete list of the M W ≥ 6.0 interplate events in each of the four regions, along with a detailed explanation if an event is not relocated, is provided in Additional file 1: Tables S1-S4. We use the locations from the International Seismological Centre (ISC) reviewed bulletin (International Seismological Centre 2019) as the initial locations for the relocation analysis; we also employ this hypocenter catalog for the mainshock hypocenters in the subsequent slip inversion analysis. The analysis performed here largely follows that of Chang and Ide (2019), with the main differences being a longer analysis period in this study (2 years longer) and a different initial earthquake catalog, as Chang and Ide (2019) used the ISC-EHB catalog (Weston et al. 2018) for comparison with a catalog containing better depth constraints.
We use three-component broadband teleseismic (station-event distance of 30°-95°) waveform data from the Global Seismographic Network (GSN) and International Federation of Digital Seismograph Network (FDSN) to perform the relocation analysis, assuming P and S waveforms on vertical and horizontal components, respectively. The raw waveforms are down-sampled to 10 Hz and detrended without removing the instrument response or filtering. We then perform the convolution, as well as signal-to-noise ratio screening (Additional file 1: Text S1). We use 44-s waveforms that begin 4 s prior to the catalog-predicted arrival time for the crosscorrelations, as in Chang and Ide (2019). Event pairs with < 30 station components are not included in the relative location analysis, with this threshold being slightly larger than the 20-component threshold of Chang and Ide (2019). The source-to-station travel time is calculated using the IASP91 (Kennett and Engdahl 1991) velocity structure.

Results and discussion
The relocated centroids of the interplate events generally exhibit mild dipping trends along each subduction direction ( Fig. 1). We retrieve the standard error (SE) of the location results and examine the consistency of the relative locations determined via bootstrapping (e.g., Efron and Tibshirani 1994). Here the same number of relative location pairs as those used in the actual analysis are randomly selected from all of the available pairs, with replacement, and the inversion is performed to determine the final locations, with this process run 5000 times (Chang and Ide 2019). The SE is then calculated as follows: where K is the number of realizations (5000), x k is the location determined by the k th realization, and x is the mean of all of the realizations. The mean SEs are 1.0, 0.7, and 1.2 km in the E-W, N-S, and depth directions, Centre 2019) and the relocated events in this study (circles) color-coded by depth. Gray squares denote the events that are not relocated (see Additional file 1: Tables S1-S4 for further details). The beach ball plots represent the focal mechanisms of the selected EGF (pink) and mainshock (red) events, with the locations of the mainshocks and focal mechanisms from the EGF events and the mainshocks taken from the GCMT catalog (Dziewoński et al. 1981;Ekström et al. 2012). The blue rectangles denote the assigned fault planes for the slip inversion analysis respectively, for all of the analyzed events. This indicates that the relative locations determined via the cross-correlation analysis are robust and consistent with each other. We then select some of the larger relocated events in each study region to be used as EGFs (Fig. 1, Additional file 1: Tables S1-S4) and perform a slip inversion for the four great megathrust earthquakes using these EGFs (see the next section). This process yields the slip relative to the event used as the EGF, and also the mainshock slip relative to all of the interplate earthquakes in the surrounding region using the accurately determined relative locations.

Method and data
We perform a slip inversion analysis to retrieve the rupture process from recorded waveforms via the Green's function. There are generally two approaches to acquire the Green's function: calculate the theoretical waveform using generalized ray theory (Langston and Helmberger 1975) or other methods, or treat the waveforms of smaller events (EGF events) that are located near the target earthquake and possess similar focal mechanisms as EGFs (Hartzell 1978). One may approximate an EGF event that is smaller than the target mainshock by about 1.5 or more in magnitude and possesses a similar radiation pattern as a point source that possesses path propagation, attenuation, and site effects similar to those of the mainshock (e.g., Baltay et al. 2014;Dreger et al. 2007).
There are also multiple ways to implement this EGFbased analysis. One common approach involves spectral deconvolution, where the relative moment rate spectrum is extracted via division of the mainshock spectrum by the EGF event spectrum in the frequency domain. The relative moment rate spectrum is converted to the time domain via an inverse Fourier transform, yielding the relative moment rate function, which is then inverted to determine the rupture process (e.g., Dreger et al. 2007;Kim et al. 2016;Mori and Frankel 1990;Mori and Hartzell 1990). Ide (2001) and Ide et al. (2011) simplified this procedure by first convolving the mainshock waveforms with triangle functions that possess a duration equal to the assumed source-time function of the EGF event to account for its rupture duration, and then expressing the Green's function simply by the appropriately shifted waveforms of the EGF event.
We implement the approach of Ide (2001) to keep the processing similar to that employed for the relocation analysis, using the same assumption for the triangle source-time function duration t d (Eq. 1) that is used for the convolution in the relocation analysis. We invert for the spatiotemporal slip rate distribution �u(ξ , τ ) on an assumed fault plane, which is expressed as follows (Ide and Takeo 1997;Okuda and Ide 2018): where a lmn are the expansion coefficients, which are determined via the inversion, and φ 1 l (ξ 1 ) , φ 2 m (ξ 2 ) , and ψ n (τ ) are the assumed basis functions on the 2-D fault in the strike, dip, and time directions, respectively. Here we assume the basis functions to be triangle functions with a 40-km length and 20-s duration. We utilize bootstrapping to regularize the inversion solutions, as in Okuda and Ide (2018), and retrieve the most robust features of the slip inversion results instead of applying commonly implemented smoothing constraints. We randomly select the same number of waveform data from the available dataset, which include those from different stations and components, with replacement, and perform a slip inversion, with this process run 300 times. We then take the mean of the model parameters and examine various features, such as the waveform recovery and slip distributions given by this averaged solution, and estimate the model error for each feature based on the standard deviation of each feature produced in each realization (Hartzell et al. 2007;Hayes 2011;Okuda and Ide 2018). We implement the non-negative constraint by applying the non-negative least squares algorithm (Lawson and Hanson 1974), as well as the causality constraint, which only allows initial slip expansion in space and time from the hypocenter after an assumed rupture velocity of 3 km/s. We use the ISC catalog locations (International Seismological Centre 2019) for the mainshock hypocenters, and set the fault planes as rectangles that intersect the hypocenters and expand in the positive and negative strike and dip directions, generally following the mainshock strike and dip as listed in the GCMT catalog (Dziewoński et al. 1981;Ekström et al. 2012). As the slip models we obtain are relative in time, location, and amplitude to the centroids of relocated EGF events whose waveforms are used as the Green's functions, we believe that the usage of un-relocated mainshock hypocenter from ISC catalog (International Seismological Centre 2019) would not affect the relative location between the slip model and the relocated neighboring seismicity noticeably.
We select several EGF events for inclusion in the analysis for each mainshock (Fig. 1) based on the following two criteria. (1) It is a relocated interplate event with a magnitude difference dM W between corresponding mainshocks that is within 1.5 ≤ dM W ≤ 2.2 , such that it can be considered a point source relative to the mainshock while retaining enough high-quality waveforms.
(2) It is located within 150 km of the mainshock hypocenter. We use one EGF at a time when we perform the inversion, such that the Green's functions for the entire fault plane of the mainshock are the same except for appropriate time shift. We can understand the variability in the results by comparing the slip inversion results for the mainshock using different EGF events, as well as extract the stable portion of the results. We use the same broadband teleseismic (30°-95°) body waveforms from the GSN and FDSN that were used in the relocation analysis to perform the inversions. We taper, detrend, and remove the instrument responses from the original data, transferring them to velocity waveforms. We include the P and SH waveforms in the analysis, using the vertical component for the P phase and the horizontal transverse component for the SH phase, which is obtained by first rotating the horizontal components of the seismograms to the great circle paths and then manually adjusting the rotation by examining their particle motion. We further down-sample the waveforms to 1 Hz, apply a 0.005-0.030 Hz bandpass filter for the Tohoku-Oki Earthquake and a 0.008-0.030 Hz bandpass filter for the other three earthquakes, and include the high-quality waveforms after checking their signal-to-noise ratios and comparing their radiation pattern to those of the mainshock and EGF event for each station (Additional file 1: Text S2). Furthermore, we exclude the data that are later than the theoretical PP-or SS-phase arrivals in each waveform from the analysis. We use the IASP91 velocity structure (Kennett and Engdahl 1991) for the travel time calculation, as in the relocation analysis, to minimize the bias that arises from different velocity models.

The 2007 M W 8.5 Bengkulu Earthquake
We begin the discussion with the smallest of the analyzed mainshocks, the 2007 M W 8.5 Bengkulu Earthquake, which is the third M W ≥ 8.5 earthquake in the Sumatra subduction zone since 2004, following the 2004 M W 9.0 Sumatra-Andaman and 2005 M W 8.6 Nias earthquakes. Fujii and Satake (2008) and Lorito et al. (2008) analyzed tsunami data to demonstrate a smaller amount of slip near the hypocenter at ~ 4.5° S, with a large-slip region at ~ 3°-4° S, to the northeast of the hypocenter. Gusman et al. (2010) analyzed interferometric synthetic aperture radar (InSAR) and tsunami data, and Konca et al. (2008) analyzed teleseismic waveforms, global positioning system (GPS) measurements, InSAR data, and measurements of coral uplift, with both studies identifying large-slip patches at ~ 3.0°-4.2° S.
We first present the waveform recovery (Fig. 2a), station distribution (Fig. 2b), and snapshots of the slip (Fig. 2c) for the solution using EGF4, which is an aftershock that occurred 2 days after the mainshock. We observe excellent waveform recovery in Fig. 2a, with the variance reduction (VR), which is a measure of the fit between two waveforms, being 90% between the observed and recovered waveforms (Eq. 4): where u obs j,k (t) and u calc j,k (t) represent the observed and calculated velocity waveforms for station j and component k , respectively.
In addition to the total slip distribution, we use the spatiotemporal slip distribution information to assist in identifying the unreliable portion of the solution, such as the scattered slip at the very end of the recovered slip history, since there can be a tradeoff between space and time during the slip inversion (e.g., Wald and Heaton 1994). We observe some initial slip near the hypocenter, followed by downdip slip to the north at ~ 3.5° S, and then updip slip to the northwest (Fig. 2c). We compare the slip snapshots using three different EGFs to better understand the robust portion of the solution (Fig. 3), revealing similar spatiotemporal slip histories: initial slip around the hypocenter that then moves downdip to the north after a decrease in slip and terminates with updip slip to the northwest (pink arrows in Fig. 3). Our results are most similar to those of Konca et al. (2008) and Lorito et al. (2008), who identified a smaller amount of updip slip around the hypocenter, larger downdip slip to the north, and larger updip slip to the northwest of the hypocenter, although the northern limit of the large-slip region in our solution is slightly to the south compared with these two studies. The full slip histories and waveform recoveries using all of the selected EGFs other than EGF4 are provided in Additional file 2: Figures S1-S5.

The 2005 M W 8.6 Nias Earthquake
The 2005 M W 8.6 Nias Earthquake occurred about 3 months after the more widely known 2004 M W 9.0 Sumatra-Andaman Earthquake. Briggs et al. (2006) and Hsu et al. (2006) used coral records and continuous GPS records to infer the slip history, with both studies identifying two large-slip regions at ~ 2°-3° N and ~ 0.5°-1.5° N. Konca et al. (2007) also included teleseismic waveforms in their analysis, and identified large-slip regions at ~ 2.0°-2.5° N, ~ 0.5°-1.5° N, and close to 0.2° N.
We first show the waveform recovery (Fig. 4a), station distribution (Fig. 4b), and snapshot of slip (Fig. 4c) using EGF4. The waveform recovery is excellent (Fig. 4a), with 94% variance reduction for a broad station distribution (Fig. 4b). We observe a larger initial slip around the hypocenter at 2°-3° N, followed by a smaller slip at ~ 1.0°-1.8° N, with the slip terminating by a patch at ~ 0.5° N (Fig. 4c). We then compare (4) VR = t,j,k u obs j,k (t) −u calc j,k (t) the results using three different EGFs (Fig. 5), where we observe that both the large-slip region at ~ 2°-3° N to the northwest and smaller slip at ~ 1.0°-1.8° N to the southeast are recovered (pink arrows in Fig. 5), although there is some variability in the slip locations among the results. Our results generally agree with those from previous studies that identified two major patches of slip (e.g., Konca et al. 2007), although the northern slip patch showed greater slip than the southern patch in our study, which is the opposite to that determined by Konca et al. (2007). The full slip histories and waveform recoveries using all of the EGFs

The 2010 M W 8.8 Maule Earthquake
The 2010 M W 8.8 Maule Earthquake occurred in the Chilean subduction zone, a highly active seismic zone with many records of megathrust earthquakes (Ruiz and Madariaga 2018 Similar to many previous studies, especially those also using teleseismic data (e.g., Lay et al. 2010), we find that the slip to the north at ~ 34°-35° S dominates our solutions using different EGFs, with a smaller amount of slip near the hypocenter to the south. We show the latter part of the slip history in Fig. 6 to highlight these observations and compare the solutions using different EGFs, where the slip to the north of the hypocenter at 34°-35° S is consistently observed across the different solutions at a similar location and timing (pink arrows in Fig. 6). For all solutions, we additionally note that this large slip is occurring in the latter part of the rupture and occurring to the north of EGF1, which is consistent with Kiser and Ishii (2013). The full slip histories and waveform recoveries using all of the selected EGFs are provided in Additional file 4: Figures S9-S13.

The 2011 M W 9.1 Tohoku-Oki Earthquake
The 2011 M W 9.1 Tohoku-Oki Earthquake off the Pacific coast of Japan is well studied due to its magnitude. Lay (2018) provides a thorough comparison and Fig. 4 Slip inversion for the Nias Earthquake using EGF4. As in Fig. 2 discussion of the many studies of slip distributions for this earthquake, with a region of very large slip within 80 km of the trench being the common feature in most of the slip images (Sun et al. 2017). Among which, Ide et al. (2011) utilizes EGF-based slip inversion analysis as in this study, using only P wave data and with smoothing constraint, yielding large slip close to the trench that extends further downdip. We also observe a large-slip region very close to the trench across various solutions, with a comparison between the spatial and temporal slip using various EGFs highlighting that the shallow slip close to the trench is well recovered using different EGFs and at a similar timing (pink arrows in Fig. 7). Our result of large slip close to the trench is consistent with many previous studies compiled in Sun et al. (2017), while the downdip slip extends mainly to about 142° E, which is shallower than Ide et al. (2011). The full slip histories and waveform recoveries using all of the selected EGFs are provided in Additional file 5: Figures S14-S17.

Short discussion of the results
We now discuss the accumulated slip distribution across the slip history (Fig. 8) based on the comparison of the spatiotemporal slip distribution histories for the four mainshocks (Figs. 3, 5, 6, 7), and extract the stable portion of the solution, with the uncertainty information determined via a comparison of the different results for the same mainshock. The main slip features for each mainshock generally agree well among the solutions using different EGFs, including the C-shaped slip distribution during the Bengkulu Earthquake, the slip region at 2°-3° N during the Nias Earthquake, the large-slip region at 34°-35° S during the Maule Earthquake, and the strong slip close to the trench during the Tohoku-Oki Earthquake (Fig. 8). We also observe in more detail that the slip distributions using different EGFs do not exhibit a clear drift toward or away from the location of the EGF event, although we do find slip along the boundary in some of the solutions (slip outside of the dotted rectangles in Fig. 8a-d), which is considered the portion of the solution that is not well resolved by our analysis and needs to be excluded from the comparison between the slip distribution and aftershock locations. We then combine the slip distribution with the relocated aftershock distribution for each mainshock, and discuss their relative locations in a more comparable sense (see the following section).

Relative locations of mainshocks and aftershocks
We now discuss the relative locations of the mainshocks and their large interplate aftershocks based on our relocation and slip inversion analyses. We first plot the different slip-distribution solutions using different EGFs and all of the interplate earthquakes with magnitudes M W ≥ 6.0 in the GCMT catalog that occurred within the 5-year period following each mainshock (Fig. 8), including a few that are not relocated (gray circles in Fig. 8; Additional file 1: Tables S1-S4) due to either the waveforms possessing a low signal-to-noise ratio or the earthquakes being too large ( M W ≥ 7.8 ) to be assumed point sources compared with the other events during the relocation. For the Bengkulu Earthquake, a group of aftershocks occurred in a region that did not experience major slip during the mainshock for all of the different slip distributions (Fig. 8a). One M W 7.9 event occurred ~ 12 h after and to the northwest of the mainshock, and another M W 7.8 event occurred ~ 3 years after the mainshock. For the Nias Earthquake, most of the large aftershocks are distributed to the northwest and south of the main rupture region, with some occurring along the periphery of the large-slip region at 2°-3° N (Fig. 8b). We find that no large interplate aftershocks occurred close to the large-slip region to the north of the hypocenter of the Maule Earthquake at 34°-35° S (Fig. 8c). We find similar results of no large interplate aftershocks in the large-slip region of the Tohoku-Oki Earthquake, as most of the slip occurred close to the trench and weakened in the downdip direction. We present the aftershock locations relative to the mainshock slip in a quantitative sense, as in Woessner et al. (2006) and Wetzler et al. (2018), to better define the mainshock slip-aftershock location relationship. We first determine the ratio of the amount of mainshock slip at each aftershock location to the maximum mainshock slip, and then plot the normalized accumulated number of aftershocks located within each threshold of the mainshock slip (Fig. 9). We perform this analysis using the slip results determined using different EGFs for each mainshock that are located within the areas indicated by dotted rectangles in Fig. 8, noting that the slip in the regions outside of the dotted rectangles occurred along the boundaries, which are not well resolved. There are no aftershocks in the regions that experienced > 40% of the mainshock slip, with most of the aftershocks occurring in the regions that experienced < 20% of the mainshock slip for the Maule and Tohoku-Oki earthquakes (Fig. 9c, d). Similarly, we find no aftershocks in the regions that experienced > 50% of the mainshock slip, with most aftershocks occurring in the regions that experienced < 30% of the mainshock slip for the Bengkulu Earthquake (Fig. 9a). However, we find it difficult to similarly infer the mainshock slip-aftershock location relationship for the Nias Earthquake (Fig. 9b) since there are a few aftershocks along the periphery of the large-slip region, with the large mainshock slip region varying slightly among the results using different EGFs. Note that the exact percentage of slip may vary with the smoothness of the slip images (i.e., how concentrated or smeared out the mainshock slips are). Here we compare the aftershocks with the concentrated features of the mainshock slip since we use bootstrapping results, as opposed to smoothing constraints during the slip inversion analysis, to produce a robust image. Some events are not relocated due to either the waveforms possessing a low signal-to-noise ratio or the events being too large ( M W ≥ 7.8 ) to be relocated since they cannot be assumed point sources compared with the other events (Additional file 1: Tables S1-S4). We do not exclude these un-relocated events in this comparison, but since they did not occur in the regions with large mainshock slip, their inclusion does not change the observations (Fig. 8).
We further attempt to assess the uncertainty of the distribution in Fig. 9 from the individual slip images within the 300 times of bootstrap by first retrieving the same curves as in Fig. 9 for the 300 individual slip images using each EGF for each mainshock, then taking the 10th and 90th percentiles and also the median of the aftershock count at each threshold to maximum slip. The results for EGF1 of each mainshock is shown in Additional file 6: Figure S18. We note that the confidence intervals we have retrieved here represents that of the 300 individual slip inversions, which are each derived using a particular dataset from which a certain combination of error is inherited, hence tends to be more concentrated and therefore naturally more likely to be separated from the aftershocks, as can be observed from the black curves in Additional file 6: Figure S18 for averaged bootstrap results lying in most cases to the left of the median or even the 90th percentiles of the specific slip images. Despite the variation between each individual case, however, we found that the observation that there is no large interplate aftershock in large slipping region or at the periphery of it does not change.
We then attempt to determine the statistical significance of this distribution of fewer aftershocks in the large-slip region, following an approach similar to that of Agurto et al. (2012) and Wetzler et al. (2018), since the lack of aftershocks in the large-slip region could simply be the result of mainshock slip regions being fairly concentrated. We randomly scatter the same number of aftershocks as in Fig. 9 onto each slip image for each mainshock, calculate the same distribution as in Fig. 9, and run this process 5000 times. We then take the mean, and the 10th and 90th percentiles of all of those distributions from all of the EGFs, and compare them with the average of the actual distribution over different EGFs (Fig. 10). We found that the regions that accommodated the large mainshock slip during the Bengkulu, Maule, and Tohoku-Oki earthquakes are generally lacking in aftershocks compared with the mean case (i.e., regions where the blue curves lie above the red curves in Fig. 10a, c, and Fig. 8 Accumulated slip using all of the EGFs and interplate aftershocks in the 5 years following the a Bengkulu, b Nias, c Maule, and d Tohoku-Oki earthquakes. The total slip of each mainshock (contours represent 35% and 70% of the maximum slip) using all of the EGFs (same color as corresponding contours of slip, locations indicated by stars) are shown for each study region. Circles represent the interplate aftershocks that occurred in the 5 years following each mainshock (relocated in black and un-relocated in gray), and are sized by their 2R values (Eq. 1). Gray rectangles denote the assigned fault planes for the slip inversion analysis, and dotted dark gray rectangles indicate the regions where slip was considered in the quantitative analysis d). However, we cannot verify the statistical significance by comparison with the 10th and 90th percentiles. This could inherently result from only including M W ≥ 6.0 aftershocks in the analysis, which limits the sample size. While this magnitude threshold can be lowered by using more regional data (e.g., Agurto et al. 2012;Kato and Igarashi 2012), this approach may introduce a bias to the comparison, as demonstrated by Rietbrock et al. (2012).

Discussion and conclusion
Here, we retrieve the robust slip features by bootstrapping over results using different station and component combination, then extract the stable portion of the slip models by comparing results using several EGFs separately. This approach offers us the reliability information of the EGF-based slip inversion analysis, as we assumed that all EGFs behave similarly, and treated them in the same manner except for their timing, location, and magnitude differences. From Fig. 8, we can observe that there are some variations between the slip models, suggesting the existence of some EGF-specific properties that exceeds our assumption, which may be investigated further in future studies.
In this study, we conducted a careful study to empirically and relatively determine the mainshock slip and aftershock locations in order to minimize the bias between the mainshock and aftershock locations. One could argue that the relocation analysis did not appear to change the conclusion greatly, as the relocated locations in Fig. 1 do not deviate greatly from the original ISC catalog locations (International Seismological Centre 2019). This is because our entire analysis was based on the ISC-coordinate, where we used the ISC catalog locations as the initial locations for the relocation analysis, and then inverted for the mainshock slip using EGF events that were relocated from the ISC catalog, with mainshock hypocenters also taken from the ISC catalog. However, if we were to compare the aftershocks determined using one catalog and the slip inversion results using another, then a potential bias between the two Fig. 9 Quantitative comparison between the aftershock locations and normalized mainshock slip for the a Bengkulu, b Nias, c Maule, and d Tohoku-Oki earthquakes. The number of aftershocks, which is normalized by the total aftershock count, that occurs at the location that experienced mainshock slip within various thresholds (vertical axes), is plotted against different thresholds of slip during the mainshock, which are the amount of slip at the site of the aftershock normalized by the maximum mainshock slip (horizontal axes) catalogs is introduced, such as the shift between the ISC-EHB catalog (Weston et al. 2018) and the Japan Meteorological Agency hypocenter catalog (Japan Meteorological Agency 2019) demonstrated by Chang and Ide (2019), which would result in two incomparable locations that could strongly influence the conclusions derived from such a comparison. This is the very bias that prevented Rietbrock et al. (2012) from making a conclusive comparison between the mainshock slip location and aftershock distributions in their study. Therefore, our approach may have the potential to better handle such bias issues in future applications.
Here, we address the potential pitfall due to biases when comparing the relative locations of the mainshock slip and the aftershock distribution. We first relocate large interplate earthquakes under the assumption that they have similar Green's functions between a given seismic station. We then use some of these relocated events as EGF events, and perform a slip inversion of the mainshock to retrieve the mainshock slip relative to the EGF events. Similar spatiotemporal slip distributions are observed using different EGF events. We qualitatively and quantitatively demonstrate that the large interplate aftershocks occurred primarily either outside of or along the periphery of the region of large coseismic slip of each mainshock. These results are in agreement with many previous studies and were obtained using a common-coordinate approach for the mainshock slip and aftershock locations, which is an effective way to minimize the bias between the mainshock slip and its aftershock locations.
Additional file 1: Text S1. Detailed explanations of the data screening process for the relocation analysis. Text S2. The inversion analysis. Tables S1-S4. list every M W ≥ 6.0 interplate event in the study regions, excluding the mainshocks, for the Bengkulu, Nias, Maule, and Tohoku-Oki earthquakes, respectively.

Fig. 10
Statistical significance tests of the aftershock locations relative to the mainshock slip for the a Bengkulu, b Nias, c Maule, and d Tohoku-Oki earthquakes. The axes are the same as in Fig. 9. The figure shows the average of the observed distribution (red) for the slip inversion results using different EGFs, as well as the mean (blue), and 10th and 90th percentiles (black) of the combined distribution of the 5000 bootstrapping realizations of random scattering of the same amount of aftershocks as in the actual observation onto each slip image using each EGF