Research ArticleGEOPHYSICS

Systematic deficiency of aftershocks in areas of high coseismic slip for large subduction zone earthquakes

See allHide authors and affiliations

Science Advances  14 Feb 2018:
Vol. 4, no. 2, eaao3225
DOI: 10.1126/sciadv.aao3225


Fault slip during plate boundary earthquakes releases a portion of the shear stress accumulated due to frictional resistance to relative plate motions. Investigation of 101 large [moment magnitude (Mw) ≥ 7] subduction zone plate boundary mainshocks with consistently determined coseismic slip distributions establishes that 15 to 55% of all master event–relocated aftershocks with Mw ≥ 5.2 are located within the slip regions of the mainshock ruptures and few are located in peak slip regions, allowing for uncertainty in the slip models. For the preferred models, cumulative deficiency of aftershocks within the central three-quarters of the scaled slip regions ranges from 15 to 45%, increasing with the total number of observed aftershocks. The spatial gradients of the mainshock coseismic slip concentrate residual shear stress near the slip zone margins and increase stress outside the slip zone, driving both interplate and intraplate aftershock occurrence near the periphery of the mainshock slip. The shear stress reduction in large-slip regions during the mainshock is generally sufficient to preclude further significant rupture during the aftershock sequence, consistent with large-slip areas relocking and not rupturing again for a substantial time.


Plate tectonics involves relatively steady long-term motions of lithospheric rock masses bounded by contact surfaces, or faults, across which the plates are offsetting (Fig. 1). Friction resists slip on the plate boundaries, resulting in accumulation of shear strain in the adjacent elastically deformed rock as plate motions proceed. When shear stress on the fault builds up sufficiently to overcome friction, the fault ruptures. The static stress drop lowers the fault shear stress to below a level sufficient to drive further slip, and the fault locks up again, awaiting stress buildup from further plate motions to repeat the process. Seismic waves radiated during rapid fault sliding can be used to determine the space-time distribution of earthquake slip, providing an estimate of the static stress drop. However, absolute stress levels are not resolved by seismic waves, and it has long been unclear how much shear stress remains within the rupture zone for large earthquakes and the associated seismic potential (1).

Fig. 1 Schematic illustration of the 1 April 2014 Mw 8.1 North Chile (Iquique) megathrust earthquake sequence.

The mainshock involved thrust faulting on the plate boundary between the underthrusting Nazca and the overriding South American plates. Locations of Mw ≥ 5.0 interplate and intraplate foreshocks within the preceding 60 days (blue) and aftershocks in the following 14 days (red) are projected on the plate boundary (37). Contour lines indicate an actual mainshock coseismic slip model (38). Bathymetry is from ETOPO1 (39). Vertical exaggeration is ×3.5. Megathrust depth regimes domains A, B, and C from the study of Lay et al. (23) are indicated.

Many studies of individual events have indicated that regions of large mainshock slip tend to have little seismicity after the mainshock, suggesting that the shear stresses are significantly reduced to well below the failure level (28). In certain cases, faults in the surrounding medium reverse from prior compressional activity to subsequent extensional activity, which requires total shear stress reduction on the megathrust. For example, the 2011 Mw (seismic moment magnitude) 9 Tohoku (Japan) earthquake substantially relaxed and rotated the local and regional deviatoric stress field (9, 10). Substantial stress drop within the main slip zone is consistent with the observation that larger aftershocks locate farther away than smaller aftershocks relative to the mainshock slip (11). Evidence of high friction on faults in the laboratory and the consistency of fault orientations with high frictional states are often cited to argue for high residual shear stress on the fault surface after an earthquake, which would indicate that the coseismic stress drop is a small fraction of the absolute shear stress (12, 13). In addition, it is known that earthquakes can be triggered by small changes in stresses and pore pressure in environments where the in situ stress is near critical (14). However, the situation can be different on a fault plane that has experienced a large and rapid slip. The extensive slip motion can produce a complex damage zone with strong heterogeneity in slip and stress. If the average residual stress remains high, then numerous aftershocks can occur at local high-stress spots. Coseismic stress changes can also modulate aftershock activity in surrounding near-critical stress state regions, either increasing or decreasing activity (15).

The recent global occurrence of major earthquakes and advanced slip-imaging methods (1619) now provide many well-determined mainshock fault slip distributions, allowing a systematic examination of observed faulting mechanism and location of aftershocks with respect to the mainshock slip (Fig. 1). We consider 101 large (Mw ≥ 7) subduction zone plate boundary earthquakes between 1990 and 2016 (fig. S1) with self-consistently determined seismic moments (M0) and coseismic slip models (table S1). We consider all large (Mw ≥ 5.2) aftershocks with focal mechanisms in the Global Centroid Moment Tensor (GCMT) catalog, examining the spatial relationship between the mainshock slip and the following seismicity both on and off the megathrust. Our goal is to use the aftershock seismicity with known faulting geometry for a large number of events to systematically evaluate the extent to which the aftershocks overlap mainshock large-slip zones and to thereby evaluate the relative degree of shear stress reduction in the mainshock slip zone for large subduction zone events.


Establishing the relationship between coseismic slip and aftershock distribution is challenging because there are uncertainties in both the slip models and the aftershock locations. In some cases, slip distributions are relatively well constrained spatially because of combined analysis of seismic, geodetic, and tsunami observations, but in many cases, only remote seismological data are available to determine the source slip distribution, and this leads to substantial uncertainty in the slip models, primarily associated with limited resolution of rupture expansion speed (20). We use source models with a suite of solutions for varying kinematic rupture expansion speed to address the uncertainty in the source slip areas. Independently, routine earthquake catalog locations have uncertainties of tens of kilometers in absolute locations of events, with shared overall regional biases and interevent scatter caused by the complex Earth structure and varying recording station distributions. Given that our slip models are tied to specific mainshock hypocenter locations, we reduce the relative location uncertainties by performing master event relocation (21) relative to the mainshock hypocenter for all aftershocks in our spatiotemporal windows. Three examples of the relocation results are shown in Fig. 2. In general, the location uncertainties are small compared to the slip model uncertainties, but throughout the analysis, we use only relocated positions of the aftershocks.

Fig. 2 Three examples of master event relocation of aftershocks.

The maps show the 14-day aftershock sequences at the initial National Earthquake Information Center (NEIC) epicenters (blue dots) and the relocated positions (red dots), connected by black lines. The mainshock NEIC hypocenters are shown by open circles. Loc, located; Reloc, relocated.

We measure the minimum in-plane distance (Δri) of the ith relocated aftershock from a contour of a specified fraction of the peak mainshock slip, χ, in a given source model (Fig. 3). Negative values of Δri indicate positions of aftershocks that locate within the slip contour, and positive values indicate aftershocks outside the slip contour. The enclosed area of the slip contour is used to compute the radius, Rj, of a circle with an equivalent area for the jth mainshock, which is used to normalize the measured distances. This scales all events for effective source dimensions, allowing us to combine data from all events to display the overall distribution of the aftershocks internal and external to the slip region. The calculated scaled event distances for different rupture area dimensions are a proxy for the radial stress changes from the finite sources.

Fig. 3 Schematic of the measurement of aftershock locations relative to mainshock slip zones.

Hypocenters for i magnitude ≥ 5.2 aftershocks are projected to each mainshock fault plane, and in-plane minimum distances (Δri) are measured from the slip contour for a given fraction of the peak slip (χ) (colored zones are shown for values of χ = 0.80, 0.50, and 0.15). Negative values are for aftershocks enclosed by the contour, and positive values are external to it. When combining measurements for different events into a composite plot, these distances are normalized by the source dimension Rj, taken as the radius of a circle with an area equivalent to that of the selected slip contour for the jth event.

Examples of spatially and radially relocated aftershock distributions relative to slip models are shown in Fig. 4 for three of the earthquakes (all of the aftershock distributions and reference slip models are shown in fig. S1). Aftershocks are classified as either interplate (shallow dipping thrust events near the plate boundary) or intraplate (for example, events with any faulting geometry off the plate boundary) using the focal mechanisms provided by the GCMT catalog (Fig. 4). Events having focal mechanism solutions with pressure, tension, and null axes all within 30° of the mainshock values are identified as interplate events. This focal mechanism information is valuable for distinguishing events located on the megathrust from those located above or below the megathrust that may appear to overlap the mainshock slip in map view. We have analyzed other seismicity catalogs that lack focal mechanism information, which intrinsically provide more limited results, but we focus here on the results for the GCMT catalog. It is clear in Fig. 4 that the number and radial distribution of aftershocks are highly variable between events. The total number of aftershocks increases with rupture area (22), so the number within an annulus should scale as the source radius. We therefore normalize the aftershock counts in each annulus by the mainshock’s radius Rj when combining the aftershock spatial patterns from all events into composite distributions.

Fig. 4 Examples for three mainshocks and the associated 14-day aftershocks plotted with the calculated distances from the selected slip contour.

(A to C) Maps of three mainshocks showing all GCMT regional Mw ≥ 5.2 aftershocks, plotted at the relocated epicenters and showing the GCMT focal mechanism. Each mainshock epicentral location is from the NEIC hypocenter indicated by the black focal mechanism. The region(s) within the χ = 0.15 contour of the slip model is (are) plotted by the light pink polygon. The areas in which aftershocks are considered are indicated by the magenta circles. Interplate and intraplate aftershocks are distinguished by the color of the compressional quadrants, red and gray, respectively. The lower right corner of each map shows a lower-hemisphere stereographic plot of the distribution of the compressional (P), tensional (T), and null (B) principal stress axes of the aftershocks with respect to the P, T, and B axes of the mainshock (solid diamonds). Events having focal mechanism solutions with P, T, and B axes all within 30° of the mainshock values are identified as interplate events. (D), (E), and (F) show corresponding histograms of numbers of aftershocks at varying minimum distances from the χ = 0.15 contour. Interplate events are plotted along the positive y axis (red bars), and intraplate events are plotted along the negative y axis (gray bars). Negative values of Δr indicate events within the slip contour. The dashed vertical line marks the position of the reference contour.

Composite radial distributions of aftershocks for all mainshocks are shown in Fig. 5, with the x axis scaled by the χ = 0.15 contour dimension. Most relocated aftershocks in the GCMT catalog are positioned near or beyond the outer edge of the mainshock rupture zones defined by the contour for χ = 0.15 for our preferred rupture models (Fig. 5A). These external aftershocks show a progressive decrease in number of events with increasing distance, with activity decaying to low levels within two source dimensions.

Fig. 5 GCMT seismicity locations relative to mainshock slip zones.

Composite distributions of relocated aftershocks in the GCMT catalog with Mw ≥ 5.2 at their minimum distance from the χ = 0.15 contour for (A) all mainshocks, (B) domain A events, and (C) domain B events, with normalization by Rj of both the radial positions, Δri, and the number counts, Nij. The horizontal dashes within the histograms separate contributions from each event. (D), (E), and (F) show the corresponding subset distributions of interplate seismicity (positive bars) and the intraplate seismicity (negative bars). The blue curves superimposed on the interplate seismicity distributions are for a circular rupture area with seismicity decaying proportional to (Δr/R)−0.5 with respect to the reference contour (Δr = 0; both internally and externally) with a base level equal to the average of the observations in the two intervals adjacent to Δr = 0. This matches the basic distribution of the composite sequence well, for Δr/R from −1 to 1, capturing the reduction of seismicity both within the reducing area inside the contour and in the expanding area outside the contour.

We subdivide the mainshocks based on the depth extent of their ruptures because this proves to be an influential factor for the relative distribution of interplate and intraplate aftershocks. We adopt the nomenclature from Lay et al. (23) of designating domain A for events that rupture (at least) the shallowest (≤15 km deep) portion of the megathrust (we include great events that rupture the entire megathrust extending to the trench in this group) and domain B for events that rupture from 15 to 35 km deep or deeper (we include events that rupture domain C from 35 to 55 km in this group). These domains are indicated in Fig. 1. The overall radial distributions are consistent for mainshocks rupturing domain A (Fig. 5B) and domain B (Fig. 5C). These composite patterns are robust with respect to the choice of varying aftershock time windows, as demonstrated in fig. S2.

The composite aftershock patterns have a simple spindle-shape distribution, with event numbers decaying asymmetrically away from the selected rupture contour, moderately externally and more rapidly internally. This corresponds to combined effects of areal extent at each distance from the contour and the decay of aftershock-triggering stress away from the contour. Before fully exploring the relative importance of these effects, we first quantify the overall pattern by considering a circular rupture with radius R with stress varying with radial distance from the perimeter Δr. For uniform thickness annuli, the area increases linearly with Δr/R external to the perimeter and decreases linearly with Δr/R interior to the perimeter. A spatial decay of seismicity with (Δr/R)−0.5 multiplied by the linear variations provides a reasonable match (blue curve) to the overall behavior for interplate aftershocks (Fig. 5D) for Δr/R from −1 to 1. For values of Δr/R > 1, the seismicity distribution falls off faster than the uniform circular distribution with a radial decay of (Δr/R)−0.5, either because the long-range spatial decay is stronger than apparent in the close-in data or because the finite width of the megathrust constrains the area where interplate aftershocks can occur to a more rectangular, one-dimensional (1D) shape. We do not attach any specific physical interpretation to the seismicity decay rate; it is just used here to account for the basic spindle shape of the observations. The shape of the observed patterns is similar for domain A and domain B subsets and for interplate and intraplate subsets, as well as for varying choices of aftershock time windows (fig. S2).

The depth extent of rupture during the mainshock has a strong effect on the relative numbers of interplate and intraplate aftershocks. Mainshocks that rupture the shallowest (≤15 km deep) near-trench part of the megathrust have a dominant percentage of intraplate aftershocks (Fig. 5E), whereas the mainshocks that rupture deeper primarily produce interplate aftershocks (Fig. 5F). The domain A interplate aftershock observations are depleted relative to the simple reference model as well (Fig. 5E). It has been proposed that this distinct behavior is likely due to rupture proximity to the free surface stress conditions (24), but for the current study, the key issue is the extent of any deficiency of interplate activity within the mainshock slip zone along with the overall majority of events falling outside the slip zone. The yet deeper (35 to 55 km deep) domain C portion of the megathrust is poorly sampled by our mainshocks but is characterized by having low aftershock activity that almost exclusively involves interplate events (fig. S3).

We seek to determine whether there is a deficient number of aftershocks in the innermost portion of the mainshock slip zones relative to the available area. To evaluate the deficiency within the rupture area, we calculate a reference distribution curve based on the actual geometry of each slip zone (fig. S1) event by event (see Materials and Methods for details about the calculation of the reference curve). Each reference distribution is normalized to the observed number of events in the interval from −1< Δr/R < 0.25 to focus on the area within the slip zone and to avoid the influence of uncertain external seismicity decay rate with distance. We then compare individual event aftershock distributions with the corresponding reference event distributions for each observed rupture topology.

Six representative examples of doing this for the χ = 0.15 contour are shown in Fig. 6 (A to F). It is evident that, in the first two to three intervals (−0.75 < Δr/R < 0.25) within the rupture at the larger-scaled distances from the contour (negative values), there are systematic deficiencies of aftershocks relative to the available area. This comparison between the modeled and observed seismicity requires sufficient aftershocks to establish both interior and overall behavior. For the 25 events with at least seven total aftershocks, deficiencies for the central region of the slip zone are found in every case, within the elevated seismicity at the rupture edge (fig. S4). Stacked observed and predicted distributions for varying thresholds of aftershock numbers demonstrate that the systematic deficiency of aftershocks is apparent for the full data set (fig. S5). Figure 7 shows that the collective deficiency, as measured by the percent summed difference between the reference curve and the observed aftershock numbers, is apparent for all threshold values on number of aftershocks, becoming larger and better defined as the number of observed aftershocks increases; thus, this is not an artifact of selection. For the region (−1.0 < Δr/R < 0.25) for our preferred rupture models, we determine the cumulative deficiency of aftershocks relative to available area to range from 15 to 45%, increasing with number of aftershocks in each sequence. We emphasize that the general behavior corresponds to the observed distributions directly apparent in map view in Fig. 4 and fig. S1; thus, the processing does not produce any artificial pattern.

Fig. 6 Examples of aftershock deficiency within the rupture zone.

Normalized distributions of relocated aftershocks, as shown in Fig. 5, for specific events for χ = 0.15. The green curves with circles are reference distributions for models calculated using the specific rupture area for each event with randomly distributed seismicity with respect to the reference contour (Δr = 0; both internally and externally) normalized so that the sum of the reference values (green) over the first five intervals from −1.0 to 0.25 is equal to the sum of the observations (red) over the same intervals. The dashed curves are SDs in the reference contours based on bootstrap analysis. The gray arrows show the deficiency in aftershocks from the predicted model values for a random distribution of the events within the slip contour. Corresponding plots for 25 events with at least seven total aftershocks are shown in fig. S4.

Fig. 7 Aftershock deficiency measurements.

Estimation of cumulative deficiency of aftershock activity within varying distance intervals from the individual earthquake rupture contours as a function of minimum number of interplate aftershocks considered. While deficiency is observed for all thresholds, it is more pronounced as the event population increases. The results presented in Fig. 6 are for a minimum number of aftershocks of seven.

Increasing the value of χ allows further examination of aftershock distribution with respect to higher-slip regions of each model. Figure S6 shows the results for χ = 0.5, with much reduced percentages of all events located within the slip contour. As the area contracts with increasing χ, the distribution of events within the slip zone will tend to approach a circular random distribution, but the number of events inside the contour quickly becomes too low to perform meaningful individual event comparisons for the specific rupture topologies.

The radial dimensions of the slip models tend to scale linearly with the assumed rupture velocity used in the inversions, producing substantial uncertainty in the rupture dimensions. To account for uncertainties in the rupture area that are affected by the uncertainties in rupture velocity, we radially scale the reference slip contour by ±20% from the original reference rupture velocity cases (fig. S7). Systematical increase of the slip dimensions increases the number of enclosed aftershocks, whereas shrinking the source dimensions decreases the number. Keeping in mind the fact that most aftershocks are distributed outside the mainshock slip zone, we focus on the internal region, computing the average ratio of the number of aftershocks within the contour versus the total number of aftershocks for that event for a range of χ values (fig. S8).

For χ = 0.5, the contours for all events that encompass slips larger than half of the peak slip, only 4 to 15% of the total GCMT aftershock populations plot within the contour. The numbers reduce even further to less than 10% when only interplate events on the megathrust are considered for the same contours. While the area of encompassed high slip decreases as (1 − χ) decreases very similarly (fig. S8), the rapid drop-off of number of enclosed events indicates that, across the population of events, very few aftershocks rerupture areas of large slip in the mainshocks.

The low level of aftershock activity overlapping the mainshock rupture zone is even more pronounced if we consider the cumulative moment ratio rather than event counts because larger-moment aftershocks tend to locate well outside the mainshock slip zone (fig. S9). The mean values of the cumulative moment of interplate aftershocks within the χ = 0.5 contour divided by the corresponding mainshock moments are very small (0.015, or 1.5%) and maintain a small ratio even when we consider a larger slip area, χ = 0.15 (0.035, or 3.5%; fig. S9A). The interplate activity is also found to be low in the coseismic slip area if we consider the average ratio of the cumulative moment of interplate aftershocks within the contour to the cumulative moment of all interplate aftershocks for each event (fig. S9B). Essentially, the aftershocks that do overlap the mainshock slip zones involve negligible additional slip for all but a handful of events. Note that these numbers are for the small fraction of events that actually have any aftershocks within the slip contours; most events do not, making the lack of aftershock moment even more pronounced.

There are relatively few intraplate aftershocks directly above or below the large-slip zone of the mainshock. Rather, most activity is distributed in an annular volume with off-plane events mainly located just beyond the margins of the slip zone. This observation provides additional evidence for a significant mainshock stress drop that modifies the relative magnitudes of the principal stresses in the volume surrounding the fault. The preferential activation of intraplate activity for ruptures of domain A demonstrates this volumetric effect; rupture to the seafloor at the trench has the potential to achieve a near-complete shear stress drop on the shallow megathrust due to the free boundary condition (24). This allows the occurrence of outer rise activity with diverse mechanisms produced by ambient stresses that are modulated by the interplate strain accumulation and release cycle. For domain B ruptures, the up-dip locked portion of the megathrust inhibits activation of outer rise intraplate activity.


The observation here that the lack of aftershock activity within the mainshock rupture zone is a robust characteristic for many large events (rather than just for a handful of instances, as has been reported in individual event studies) has important implications. The notion that large events release much of the stress accumulated since prior large ruptures implies that it will take time for sufficient stress to rebuild before having another rupture, making the region of the fault with the largest slip relatively safe after large earthquakes. The lack of aftershock activity within the mainshock high-slip zone (Fig. 6) suggests that the mainshock released sufficient accumulated shear stress in the relatively high-slip area of the rupture to bring the stress down to well below the critical stress. Pragmatically, this observation suggests that an observation of intense aftershock activity located within the high-slip rupture area would imply that a larger earthquake may still be possible. Only a very small subset of our mainshock ruptures overlap areas that had prior earthquakes that are of comparable size or larger within the seismological record. Thus, we do not strive to relate the degree of aftershock deficiency to the irregular seismic cycle experienced in each area because the sampling is too limited, and the pattern we interpret as favoring substantial stress drop in the large-slip areas emerges directly from the overall behavior.

However, competing observations of earthquake triggering and catalog statistics have called the intuitive notion of significant stress reduction based on observation of aftershocks not lying in large-slip zones into question, prompting a major debate about the validity of seismic gaps (25) (previously seismogenic fault regions where no large earthquake has been observed for a long time) as an indicator of future hazard (26). The observation that the spatial distribution of subsequent activity is influenced by the location of a large slip in the mainshock (fig. S1) reinforces the idea of stick-slip strain accumulation and release with a significant stress drop. The facts that the dynamic and static stress changes from the ensuing aftershock sequence fail to drive rerupture of most of the large-slip region of the mainshock and that only very low total moments are involved for the small fraction of events that do overlap indicate that the shear stress reduction is large and pervasive over the rupture surface. This possibility is supported by the 2011 Mw 9.0 Tohoku earthquake, as discussed above, where deviatoric stress resolved on the fault appears to reverse sign, although the coseismic stress change is itself low (4, 810, 27, 28).

We visualize the essence of aftershock spatial behavior with respect to a mainshock coseismic slip contour by considering simple 2D scenarios (Fig. 8). Characterization of the radial aftershock pattern is based on the number of earthquakes enclosed by expanding annuli of constant width from the mainshock epicenter/rupture perimeter. In cases where aftershocks are randomly distributed over the fault with no influence of the mainshock stress changes, we expect a quasi-linear variation of seismicity with respect to distance from the mainshock epicenter/rupture perimeter, concave upward inside the irregular rupture area (Fig. 8A). If the variable mainshock slip reduces stress in internal large-slip areas but increases it on the margins and external to the slip zone, then we expect a deficiency of aftershock activity within the core of the coseismic slip area, surrounded by a concentration near the perimeter that decays with increasing distance due to stress decay, out to a distance of about two times the source radius (Fig. 8B). If the stress drop is near total over the entire slip area, then the internal deficiency should be close to complete out to the rupture perimeter (Fig. 8C). The difference between Fig. 8B and 8C is certainly hard to resolve given the uncertainty in slip areas and aftershock locations. If both are very well determined, as for the Mw 5.2 Borrego Springs earthquake (5) or for the great 2011 Mw 9.1 Tohoku earthquake (6), then a stronger case can be made for total stress drop, but this requires very extensive, nearby station distributions. In our situation, when considering individual events and composite behavior for many events together, some of which could be in any of the three categories, it is clear that the overall behavior is best represented by the case in Fig. 8B, with the central region of the slip zones likely having a significant stress drop to well below the critical stress, and in some cases, total stress drop may have essentially occurred, while the margins experience stress increases that produce a perimeter (or torus, if intraplate events are considered) of activity that may or may not rerupture low slip regions of the mainshock.

Fig. 8 Schematic of possible aftershock distributions relative to a mainshock slip zone.

Seismicity is shown by the black dots, and the coseismic slip area is marked by the magenta contour and shaded regions. (A) Randomly distributed aftershocks unaffected by the specific mainshock slip distribution. (B) Aftershock distribution relative to a substantial stress drop mainshock that still has residual stress within the margins of the slip zone. (C) Aftershock distribution relative to a total stress drop mainshock that has no shear stress remaining in the slip zone. (D to F) The corresponding schematic radial distributions of aftershocks relative to the margins of the slip zones for the three cases. The black dashed line is the radius of the slip zone, as used in all previous measurements (that is, Fig. 5), and the green dashed line is the predicted level of possible activity for a random distribution of aftershocks within the slip contour.

Absolute stress levels remain very difficult to resolve, but the systematic behavior displayed by the 101 large megathrust sequences considered here does support the basic notion that large shear stress reduction occurs in the mainshock rupture zones, and it will require time for reaccumulation of sufficient stress to rerupture those interior regions. Thus, immediate rerupture of the precise large-slip zone is unlikely (this is convincingly supported by observations), although regional clustering may occur due to stress increments on regions outside the main slip zone. This supports a simple fundamental behavior of stick-slip sliding for very large events; however, the heterogeneity of slip and stress change is undeniable, as is the existence of triggering interactions that advance or delay the time to rerupture, giving rise to observed overall complexity of large earthquake sequences.


The mainshock data set comprises 101 large (Mw ≥ 7) subduction zone plate boundary earthquakes between 1990 and 2016 with self-consistently determined seismic moments (M0) and coseismic slip models (16, 20, 2935). The slip models and seismicity sequences are shown in fig. S1, a bundle of figures for all mainshocks. We excluded a handful of mainshocks that have very close spatial or temporal proximity to preceding large-magnitude earthquakes. We considered all large (Mw ≥ 5.2) aftershocks from the GCMT catalog, for which we used initial hypocentral and phase information from the U.S. Geological Survey Preliminary Determination of Epicenters catalog.

The uncertainty of the absolute hypocenter locations was influenced by the use of a 1D velocity model, network azimuthal coverage, and number of detected phases. We were mainly interested in the position of our events relative to our slip models, which were similarly pinned to uncertain mainshock hypocenters. Thus, we adopted a master event–relative relocation procedure by minimizing the differences between the differential travel times (predicted − observed) for each aftershock and the corresponding station differential times for the mainshock (21). Each event was relocated by a grid search over a range of possible epicenters, separated by spacing dx and dy, accounting for ray parameter changes with distance from each station. Then, we found a minimum value forEmbedded Imagewhere tO − tP denotes the residual time for the mainshock (tm) and relocated event (te). The minimum value over the grid defines the new epicenter location.

Because the activation zones of aftershocks scale with rupture dimensions (22), we defined the aftershock region using an adaptive spatial window that scales with the mainshock seismic moment. We used the Eshelby crack model (36) to capture aftershocks within a circle from the mainshock hypocenter with radius Embedded Image, where a uniform static stress drop, Δσ = 3 MPa, was assumed. We used 14-day time windows after each mainshock for aftershocks for the primary analysis. The GCMT catalog was formally complete down to Mw ~5.2 throughout the time period considered (24).

For the GCMT solutions, aftershocks were classified as interplate (shallow dipping thrust events on/near the megathrust) or intraplate (all mechanisms off the megathrust) based on their individual focal mechanisms. We designated this by comparing the angles between the P, T, and B principal stress axes for each earthquake with those of the mainshock. Earthquakes designated as interplate are required to have all three angles be less than 30° from the mainshock values; hence, they are shallow-dipping thrust events. Source depth was also considered, but due to fairly large uncertainty, this was a secondary criterion because there was scatter in depth estimates, but we excluded thrust events that were located too far away from the megathrust to be interplate ruptures. Figure S1 shows the GCMT aftershock mechanisms and relocated hypocenters relative to the mainshock slip distributions for all sequences for our preferred rupture models.

We determined the spatial distribution of the seismicity after each mainshock relative to the models for the coseismic slip distribution. For a given fraction of slip relative to the peak slip, χ, we defined a contour (or contours) on the slip model separating internal higher slip regions from external lower slip regions. As χ reduces, more of the mainshock slip region is enclosed. Given the discretization and damping of the slip inversions, values of χ lower than ~0.15 involved very poorly resolved areas of the model. We projected seismicity locations perpendicularly onto the megathrust fault plane and computed in-plane minimum distances of each event location from the slip contour (Δri), with negative values for events within the contour and positive values for events outside the contour (Fig. 3). To compare the spatial distributions for mainshocks with different magnitudes, the distances from the χ slip contour were normalized by Rj, the radius of a circle with an area equivalent to that enclosed by the slip contour for the jth event. Normalizing the distances by the source dimension produced a proxy for the static stress change outside the rupture zone. We then examined the measurements for all the mainshocks together, normalizing the individual event counts by Rj to balance the event population contributions for varying size events (Fig. 5).

To calculate the background reference seismicity distributions, we randomly scattered 10 events over the entire search area (magenta circle in fig. S1) for 500 realizations. This bootstrap procedure provided the reference seismicity curve (green line in Fig. 6) calculated from the mean distribution and the SD for each histogram bin (dashed lines in Fig. 6).


Supplementary material for this article is available at

table S1. Mainshock information.

fig. S1. Bundle of event sequences for 101 mainshocks.

fig. S2. Seismicity locations relative to mainshock slip zones for different time windows.

fig. S3. Seismicity locations relative to mainshock slip zones.

fig. S4. Bundle of 25 event-specific slip deficiency estimates.

fig. S5. Cumulative distributions for varying aftershock thresholds.

fig. S6. GCMT seismicity locations relative to mainshock slip zones.

fig. S7. Examples of rupture dimension modification with ±20% change in modeled rupture velocity.

fig. S8. Seismicity locations within the mainshock slip zone.

fig. S9. Cumulative seismic moment distributions.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank L. Ye for access to mainshock slip models, which are available from We thank two anonymous reviewers for helpful suggestions that improved the presentation. Funding: N.W. was supported by the Geological Survey of Israel at the Ministry of Energy and Water Resources. T.L. was supported by NSF grant EAR1245717. Author contributions: N.W. performed the computer analysis. T.L. conceived the analysis of aftershock patterns relative to mainshocks and contributed to the analysis design. E.E.B. and H.K. contributed to the analysis design. All authors contributed to the writing of the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. The GCMT earthquake catalog is obtained from
View Abstract

Stay Connected to Science Advances

Navigate This Article