Research ArticleGEOLOGY

Do slow slip events trigger large and great megathrust earthquakes?

See allHide authors and affiliations

Science Advances  31 Oct 2018:
Vol. 4, no. 10, eaat8472
DOI: 10.1126/sciadv.aat8472


Slow slip events have been suggested to trigger subduction earthquakes. However, examples to date have been poorly recorded, occurring offshore, where data are sparse. Better understanding of slow slip events and their influence on subsequent earthquakes is critical for hazard forecasts. We analyze a well-recorded event beginning 6 months before the 2012 Mw (moment magnitude) 7.6 earthquake in Costa Rica. The event migrates to the eventual megathrust rupture. Peak slip rate reached a maximum of 5 mm/day, 43 days before the earthquake, remaining high until the earthquake. However, changes in Mohr-Coulomb failure stress at the hypocenter were small (0.1 bar). Our data contradict models of earthquake nucleation that involve power law acceleration of slip and foreshocks. Slow slip events may prove useful for short-term earthquake forecasts.


Since their discovery more than two decades ago, it has been suggested that slow slip events (SSEs) may trigger subduction zone earthquakes, perhaps by stress loading of adjacent sections of the fault (14). Alternatively, SSEs reduce the probability of large earthquakes by relieving strain and reducing the magnitude of coseismic slip (5, 6). It is also possible that both are true: SSEs limit rupture area, reducing the long-term risk from earthquakes, but elevate the short-term probability of a seismic event through perturbations of the stress state. Unfortunately, it has been difficult to test these hypotheses as, in subduction zones, the critical region occurs offshore, where geodetic networks have limited sensitivity. Offshore geodetic techniques exist, but their deployment has been limited because of high cost and lower precision (7).

In Japan, an SSE may have triggered the giant earthquake of 2011 (810). While Global Positioning System (GPS) did not record the offshore SSE, eight offshore pressure sensors recorded deformation of the seafloor associated with an SSE (10). However, the signal was close to the noise level of the technique. Seismic activity associated with SSEs was observed to propagate toward the epicenter (9). Similar observations were made, leading up to the moment magnitude (Mw) 8.1 Iquique, Chile earthquake in 2014 (11, 12), with indicators preceding the earthquake by 8 months (13). An SSE was also observed several days before the Mw 6.9 Valparaiso, Chile earthquake in 2017 (14). Table S1 summarizes the current database for earthquakes associated with SSEs. In most cases, slow slip is postulated on the basis of migrating foreshocks or a few geodetic stations. Hence, it has been difficult to investigate the physical mechanism linking slow slip to earthquakes. In many subduction zones, SSE repeat times are short (one to several years) compared with earthquake recurrence intervals (30 to 500 years or longer), making it possible that their correlations are coincidental. In particular, it has been difficult to show that SSEs migrate in the vicinity of the earthquake nucleation point based on geodetic measurements.

In late February 2012, an SSE began in northern Costa Rica, 6 months before the 5 September 2012 Mw 7.6 earthquake. Both the SSE and earthquake were well recorded due in part to a peninsular region that allows instrumentation immediately above the seismogenic zone (15, 16). Preliminary analysis of nine GPS stations noted that the change in Mohr-Coulomb failure stress (MCS) associated with the SSE at the nucleation site of the earthquake was small (6). Since that time, data from 11 additional stations have become available (figs. S1 to S6), allowing considerable refinement of our knowledge about this episode, including the relationship with foreshocks (17), migration pattern, and a more accurate estimate of MCS.

The Nicoya Peninsula lies along the Middle America trench where the Cocos plate is subducting beneath the Caribbean plate at a rate of ~9 cm/year (18). The peninsula extends toward the trench, with the plate surface of ~15 km beneath the coastline. SSEs have been identified both updip and downdip of the peninsula, with recurrence times of about 22 months (19, 20). On 5 September 2012, a Mw 7.6 earthquake took place within a region that had been previously identified as a locked patch (Fig. 1) (16, 21). Using newly available GPS data, a large geodetic signal consistent with an SSE (southwest motion) is visible mid-February, persisting until the day of the Mw 7.6 earthquake (Fig. 2). The signal appears first on stations located southeast of the peninsula in February (CIQU, JACO, PUNT, and RIDC). The initial phase lasts for 4 months, with interseismic-like GPS velocities reappearing in July, 1.5 months before the Mw 7.6 earthquake. A second pulse of motion is observed in early August and continued until the earthquake. Coastal stations (GRZA, EPZA, and SAJU) show a final increased south-westward movement 2 weeks before the mainshock.

Fig. 1 Map of the study area in Nicoya Peninsula, Costa Rica.

Light pink areas are regions with interseismic SSEs (6). Red arrow represents Cocos-Caribbean convergence direction (39). Dashed line marks the transition between oceanic crust from Cocos-Nazca spreading (CNS) center and East Pacific Rise (EPR). Blue contours mark the slab depth (18). Yellow triangles mark the GPS stations. Mainshock focal mechanism is indicated by red beach ball (15). Red star marks the epicenter of the 2012 El Salvador earthquake, 450 km to the northwest of the Nicoya Peninsula.

Fig. 2 Average slip rate and GPS time series.

Average SSE slip rate for different time periods (left), color coded to match individual site displacement time series (right). Black lines showing modeled fit due to fault slip, and scatter showing the horizontal displacements. Red triangle (left) marks the epicenter of the 2012 earthquake (15).


We invert the GPS displacements for the time-dependent slip on the megathrust (22, 23). We include the vertical GPS displacements in the datasets, weighting them ~3 times less than the horizontal position estimated. Incorporating the vertical GPS displacements is critical for estimating the downdip limit of slip, and tests removing the vertical GPS time series were unstable. Model results indicate that the initial transient started southeast of the peninsula, under Herradura (Fig. 2 and movie S1). Slip rates were highest during this initial transient, reaching 5 mm/day (fig. S6). This transient is similar to other deep Nicoya SSEs, with slip magnitudes peaking at ~6 cm (Mw 6.5) (5, 20). The transient migrates to the northwest where it slowly decays beneath the locked zone. About 3 weeks before the earthquake, a second shallow slip pulse appears in the vicinity of the locked zone, which is the rupture area of the 5 September earthquake. Kinematic modeling indicates that the earthquake nucleated near or immediately downdip of the region of shallow SSEs (15). The migration of slip toward the seismogenic zone is constrained by the displacement of coastal stations, particularly stations SAJU and GRZA. Peak slip rate occurs on this shallow section, immediately before the 2012 El Salvador earthquake, which occurred 10 days before and 450 km to the northwest of the Nicoya earthquake. Slip rates decay following this event despite a brief increase in seismicity associated with the El Salvador event (17). In the month preceding the Nicoya earthquake, there is a good correlation between foreshock productivity (defined as any earthquakes in the 30 days before the mainshock) and slip rate, with foreshocks clustering updip of the region of the largest aseismic slip (Fig. 3). Earthquake productivity in the Nicoya region during the interseismic period, before the events mentioned here, is 25 (±5) earthquakes per day.

Fig. 3 Slip rate and foreshock activity before the Nicoya earthquake.

Gray bars are daily earthquake counts (17). Blue line is the modeled slip rate. Purple vertical line is the 2012 El Salvador earthquake. Red vertical line is the 2012 Nicoya earthquake.


To explore possible triggering by the SSE, we used Coulomb failure stress (CFS) analysis (24). CFS represents the relative contributions of normal stress change and shear stress change, resolved onto the fault in the direction of fault slip. Increased CFS implies that a fault is closer to failure. We find that, while positive in the time-dependent case, CFS due to the SSE is less than 0.1 bar at the nucleation point (Fig. 4), well below the threshold typically found in static earthquake–triggering cases [~1 bar (25)] but above the level for modulating SSEs (26). While CFS is affected by the slip gradient, which is dependent on regularization, it is unlikely to have been increased by an order of magnitude. We find that reduction of regularization constraints requires more slips in the coseismic region, leading to further decrease in the Coulomb stress near the earthquake nucleation point. To explore the effects of this regularization, we also modeled the SSE using a static approach (Fig. 5). We used the cumulative offset estimated through fitting a cubic spline through the GPS time series. The time period spans the same period analyzed in the time-dependent modeling. This approach gives qualitatively similar results to the time-dependent inversion mentioned earlier in the region of the peninsula and requires slip within the cosesimic region regardless of choice of regularization (figs. S7 and S8). Static modeling requires slips to the south of the peninsula, in a region of poor geodetic resolution. We attribute this slip to differences in network geometry for the static inversion, where discontinuous time series cannot be included. Some signal at station JACO is identified as a bench mark motion in the Network Inversion Filter (NIF) modeling because of its early onset, and its displacement might be overestimated in the static modeling. However, slip is required in this region even when JACO is excluded.

Fig. 4 Cumulative slip and Coulomb stress.

Left: Cumulative slip for the SSE, with chartreuse line marking the 5-mm contour. Right: CFS on the megathrust fault associated with the cumulative slip. White and cyan contours mark the 1, 2, and 3 m of coseismic slip, and red focal mechanism marks the 2012 earthquake epicenter (15). Gray circles are foreshocks in the 30 days before the earthquake, with denser concentrations appearing black (17).

Fig. 5 Static inversion modeling of GPS offsets associated with SSE. Slip exists within the co-seismic region.

Left: Static inversion of GPS time series with preferred regularization weighting. Right: Coulomb stress associated with the static inversion modeling results.

Perhaps a more dynamic process is responsible for triggering. The presence of slow slip, both updip and downdip of the seismogenic zone during most Nicoya interseismic SSEs (6, 20, 27), suggests that slow slip is able to cross the frictional barrier between seismic slip and aseismic slip. The magnitude of shallow slip is above our uncertainty levels, although resolution decreases offshore (fig. S9). At 60-km depth, where the SSE initiated, the presence of fluid (28) is thought to promote aseismic slip, as evidenced by the large signal inland. Perhaps the SSE allows fluid to migrate updip toward the seismogenic zone at 15- to 25-km depth. While conditions within the seismogenic zone are unfavorable for aseismic slip, fluids could weaken this region through dilatancy hardening (2) or other mechanism. Further updip, conditions once again favor slow slip and do so persistently (6, 20). If the seismogenic zone is close to failure, as was the case for 2012 earthquake when the Nicoya segment was more than 20% past its average 50-year characteristic recurrence time, then additional fluids driven by the SSE were sufficient to trigger the earthquake.

Migrating seismicity before large earthquakes has been reported (814). It has been hypothesized that this migration is indicative of aseismic slip behavior. Our results provide strong evidence that foreshocks can be temporally associated with the slip rate of SSEs. Precursor microseismicity rates may therefore be a reasonable proxy for aseismic slip behavior, albeit at lower slip rates than previously reported (29). Our data also allow us to rule out at least one model for earthquake triggering by SSE, whereby rupture is initiated through power law acceleration of slip (3032). In this case, both slip rates and foreshock rates were higher in the weeks before the rupture.

While there seems to be a strong case for temporal correlation between foreshocks (seismic behavior) and slow slip (aseismic behavior), spatiotemporal patterns of seismicity do not necessarily track SSE behavior. Notably, foreshocks cluster near the megathrust rupture but predominantly outside the SSE region (Fig. 4).

The use of foreshocks and SSEs for earthquake forecasting remains challenging, as most subduction zones lack the necessary monitoring. In particular, identifying foreshock sequences that culminate in a large earthquake in real time has not been possible. When the precursor SSE in Nicoya initiated, it was indistinguishable from interseismic SSEs, which often begin with a high slip under the Gulf of Nicoya. Both the timing of the SSE [22-month recurrence (19, 20)] and earthquake [50-year recurrence (16)] were consistent with historical records. This suggests that near-term hazard forecasts should incorporate information about the timing of SSEs, particularly as a fault enters the later stage of the earthquake cycle (4). We note that SSEs are not required for a nucleation of a megathrust earthquake. Better measurements of the offshore region of subduction zones will be required to separate precursor activity from normal interseismic behavior.


GPS processing

GPS data were processed using GIPSY-OASIS 6.4 software, with orbits and clock estimated provided by the NASA Jet Propulsion Laboratory. Daily solutions were calculated using the precise point positioning method (33). Phase ambiguity resolution was performed (34), and ocean loading was removed using FES2004 (35). Tropospheric delay was estimated using the VMF1 (Vienna mapping function 1) mapping function (36), and ionospheric delay was estimated using the Ionex model (37). Stations were processed in the IGb08 reference frame.

Seasonal signal removal

Seasonal signals were removed from the time series via a least-squares fit of a function that included annual, semiannual, a linear term, and H function (offset). Periods in which SSE was occurring were masked, and a Heaviside function was used to remove the offset. Parameters were estimated using the Levenberg-Marquardt algorithm, as implemented in the Python package lmfit ( The seasonal terms were then subtracted from the time series, and the estimated trend was included as an a priori trend during the subsequent time-dependent inversion for fault slip.

Fault slip inversion using the NIF

We used a modified version of the NIF (38). We generated Green’s function using a high-resolution model of the subducting interface of the Cocos plate (39) discretized into triangular elements (40), with approximate spatial dimensions of 20 km2. The modeled times series were a function of the slip rate on the plate interface, network error, random walk error, and common-mode error. Regularization was imposed on both the temporal and spatial smoothness and was enforced via maximum likelihood (29). The slip rake direction is constrained to be parallel to the relative motion between the Cocos and Caribbean plates (29° east of North) (17). Modeled displacement estimates were compared to observed position time series in fig. S1.


Supplementary material for this article is available at

Fig. S1. Station time series and model fits for CABA, CIQU, ELVI and EPZA.

Fig. S2. Station time series and model fits for GRZA, HATI, HORI, and HUA2.

Fig. S3. Station time series and model fit IND1, JACO, LAFE, and LEPA.

Fig. S4. Station time series and model fits for LIBE, LMNL, NICY, and PUJE.

Fig. S5. Station time series and model fits for PUNT, RIDC, SAJU, and VERA.

Fig. S6. Comparison of moment rate and slip rate.

Fig. S7. L curve for choice of regularization parameter.

Fig. S8. Comparing slip distributions and CFS from different regularization enforcement parameters.

Fig. S9. Uncertainty estimates for cumulative slip.

Table S1. Comparison of aseismic precursors to large subduction earthquakes.

Movie S1. Movie of the 2012 SSE slip rate history.

References (4245)

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 UNAVCO for continued station maintenance and the University of Costa Rica (UCR) for providing access to their database. We thank the Instituto Geográfico Nacional de Costa Rica for access to additional GPS stations. Figure 1 was created using Generic Mapping Tools (41). Funding: This work was supported by grant NSF-EAR 1345100 to T.H.D. Research at the Jet Propulsion Laboratory, California Institute of Technology, was supported by NASA’s Earth Surface and Interior focus area. Author contributions: N.V. analyzed the data and wrote the first draft of the manuscript. T.H.D. suggested the study and wrote the original proposal that funded the work. R.M. processed the GPS data. M.P. collected the GPS data and maintained the network. Z.L. generated the analysis tools. All authors contributed to writing and editing the text. 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 used in this study are available through the UNAVCO data archive without restriction and from UCR upon request. Software developed for this study may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article