Research ArticleGEOPHYSICS

Repeating earthquakes record fault weakening and healing in areas of megathrust postseismic slip

See allHide authors and affiliations

Science Advances  05 Aug 2020:
Vol. 6, no. 32, eaaz9317
DOI: 10.1126/sciadv.aaz9317


Repeating earthquakes (REs) rupture the same fault patches at different times allowing temporal variations in the mechanical behavior of specific areas of the fault to be interrogated over the earthquake cycle. We study REs that reveal fault weakening after a large megathrust earthquake in Costa Rica, followed by fault recovery. We find shorter RE recurrence intervals and larger slip areas immediately following the mainshock that both gradually return to pre-earthquake values. RE seismic moments remain nearly constant throughout the earthquake cycle. This implies a balance between fault weakening (reducing slip) and transient embrittlement (increasing rupture area by converting regions from aseismic to seismic slip), induced by the increased loading rate following the mainshock. This interpretation is consistent with positive, negative, and constant moment versus RE recurrence interval trends reported in other studies following large earthquakes and with experimental work showing slip amplitudes and stress drop decrease with loading rate.


Subduction zones produce large and catastrophic earthquakes and tsunami. The recurrence interval and rupture characteristics of these megathrust events are strongly influenced by fault strength (resistance to shear stress before slip) and frictional properties of the plate interface (1).

Monitoring the interseismic fault strength and its variability (e.g., fault weakening and fault restrengthening or frictional healing) during the earthquake cycle is critical for earthquake forecasting and hazard assessment. Although an increasing amount of laboratory experiments [e.g., (2)] and numerical modeling (3) have been conducted to understand fault strength behavior during the earthquake cycle, observations in nature are limited and not in full agreement with experimental or theoretical results, hindering our understanding of fault physics (1, 3, 4).

The increasing quantity and resolution of seismological and geodetic data at subduction zones provide a heightened opportunity for monitoring and understanding earthquake processes along megathrust faults. Multiple studies have been conducted in this environment to better understand how these events propagate and affect the adjacent off-fault areas and to characterize the mechanical properties (e.g., roughness and frictional behavior) of the asperities where these events nucleate. Plate interface motion generates a broad range of events that span different spatial and temporal scales, including small- and intermediate-magnitude earthquakes (M ≤ 6.5), repeating earthquakes (REs), megathrust and tsunami earthquakes (M ≥ 7.0), low-frequency earthquakes, tectonic tremor, and slow slip events (SSEs).

REs, earthquakes with nearly identical waveforms, magnitudes, and locations have been previously recognized and studied in the laboratory and observed in multiple natural fault systems in California (58), Japan (914), Mexico (15), Greece (16), Taiwan (17), and Costa Rica (18). They are thought to represent the repeated rupture of the same fault asperity at a relatively constant recurrence interval (tr) due to continuous loading by surrounding aseismic slip (19, 20). These events are ideal to investigate temporal variations in source parameters and changes in the mechanical properties of the fault zone through the megathrust earthquake cycle. Large earthquakes substantially perturb the occurrence of REs on the same fault. Observations show that the frequency of existing families increases (reduction of recurrence interval, tr), and new sequences appear following large earthquakes (3, 6, 13). Variation in tr of postseismic repeaters has been widely documented with a range of relationships between Mo and tr reported. REs along the San Andreas fault system and off-Kamaishi, NE, Japan have revealed instances where Mo scales positively (Mo increases when tr increases), negatively (Mo decreases when tr increases or vice versa) or remains invariant with tr (7, 13, 21), and this variable scaling may be depth dependent (6). Intriguingly, observational and theoretical expectations of seismic moment (Mo) and tr scaling of REs differ, even during periods unaffected by stress perturbation from large earthquakes that occur on the same fault. Observations of REs in multiple locations find trMo1/6(5, 7) during periods of constant loading rate, but this is incompatible with the assumption that stress drop is independent of moment, which results in the relationship: trMo1/3 (see the Supplementary Materials) (8, 22, 23). Numerical models show that this conundrum can be solved if the rupture area of the REs involves an increasing fraction of aseismic slip as moment decreases (3, 23), but they can only match observed rates of recurrence by including enhanced dynamic weakening or patches of elevated normal stress into rate and state friction models (1, 2, 24).

The diverse behavior between Mo and tr is likely due to the increased loading rate following the mainshock, which introduces some competing factors: (i) It may enlarge rupture area due to a phenomenon termed transient embrittlement (6), an aseismic-to-seismic transition in frictional properties within fault patches, described by rate and state friction laws, and (ii) it may reduce the amount of slip and healing time (or tr) due to fast loading from afterslip. To resolve these effects, it is essential to separate out the variations in slip and rupture area that are combined in the moment measurements. Although previous work has attempted this, isolating the source process from wave propagation effects at relatively high frequencies (small-magnitude earthquakes) has proven difficult (8). The robustness and reliability of results depend on making observations with good azimuthal instrumental coverage and proximity to the source as earthquake magnitude reduces. In this study, we use extremely well-recorded REs on the Costa Rica megathrust to investigate these competing factors.


The advantageous location of the Nicoya Peninsula, in northern Costa Rica, extending seaward over the seismogenic zone is an ideal locality to study megathrust earthquake source processes (Fig. 1). A regional seismic and continuous Global Positioning System (GPS) network has recorded megathrust activity with outstanding detail over the past 15 years. The relatively rapid subduction (~85 mm/year) of the Cocos Plate beneath the Caribbean Plate along the Nicoya Peninsula has generated M7+ earthquakes every 50 years, with the last of these events, a moment magnitude (Mw) 7.6 earthquake, occurring on 5 September 2012 (25). Tectonic tremor and SSEs have been recorded every 2 years or less since 2003, occurring both along the up-dip and down-dip edges of the seismically active plate interface (26).

Fig. 1 Tectonic setting and instrumentation.

Map view showing the spatial distribution of the Nicoya Seismic Cycle Observatory and Volcanological and Seismological Observatory of Costa Rica (OVSICORI-UNA) seismic stations used in this study (orange triangles). The color contours illustrate the coseismic slip area of the 5 September 2012, Mw 7.6, Nicoya Peninsula earthquake, slip increments of 0.5 m (47). The yellow regions up-dip and down-dip of the plate interface highlight the accumulated slow slip (in mm) between 2007 and 2017 (26). The black dashed line represents the East Pacific Rise (EPR)–Cocos Nazca Spreading (CNS) center boundary. Stars show the epicentral location of the RE clusters analyzed in this study color coded by cluster identity (id). Adjacent focal mechanisms show the faulting geometry of the target event in each cluster. Table S1 shows the corresponding hypocenter location and origin time for the five families. The focal mechanism shows the fault geometry and epicentral location of the 2012 Nicoya Peninsula earthquake (25). INDI, Punta Indio broadband seismic station.

Using REs, we monitor the temporal evolution of the frequency content of the earthquake source spectra and estimate the stress drop at the same locations through the earthquake cycle. We observe a decrease in stress drop and infer a reduction in yield strength on RE fault patches after the 5 September 2012, Mw 7.6 Nicoya Peninsula Costa Rica megathrust earthquake, followed by recovery over a time scale consistent with the postseismic slip (27, 28).

Identifying REs

A previous study (18) using waveform matching to increase the number of events in the foreshock/aftershock catalog of the 5 September 2012, Mw 7.6 Nicoya Peninsula earthquake found 53 RE clusters between July and December 2012 located in areas of modeled coseismic and postseismic slip (27, 28). We select REs from this catalog supplemented by additional aftershocks (29) and use a similar waveform matching technique to greatly extend the time period of REs from 2007 to 2018 (see Materials and Methods for selection criteria and processing). We identify five RE families, located on the plate interface (18, 29) and consisting of (i) large or target events (magnitudes between M2.6 and M3.3) that repeat with a nearly constant magnitude between two to six times during the 12 years analyzed and (ii) several smaller events (at least one order of magnitude smaller than the largest event) with nearly identical waveforms to the large event that cluster around it in time. Figures 1 and 2 show faulting geometries and the spatiotemporal distribution of these families labeled C60, C55, C32, C23, and C16, maintaining the former nomenclature (18).

Fig. 2 Temporal distribution of the repeating families analyzed in this study.

(A to E) Magnitude distribution with respect to the time of occurrence of the 5 September 2012, Mw 7.6 Nicoya Peninsula earthquake (red dashed lines). Repeating earthquake clusters (RECs) for each family are formed by a target event (stars) and several smaller-magnitude events or EGFs (circles). (F) Spatial distribution of the REC with respect to the short-term afterslip, indicated by the white contour lines (27), and the long-term afterslip of the Nicoya Peninsula earthquake (28), indicated by the colored map. The amount of slip for both the short- and long-term afterslip is indicated in meters. Note the spatial decorrelation between the areas of maximum long-term afterslip and the location of the repeating families.

All RE families locate along the up-dip margin of the rupture zone of the 2012 Nicoya Peninsula earthquake, at the frictional transition between unstable sliding (seismic) down-dip, to stable sliding (hosting slow slip), up-dip (Fig. 1). This region also experienced maximum early postseismic slip (Fig. 2F) (27, 28) and positive Coulomb static stress changes (29) following the Nicoya mainshock. Three of the five earthquake families (C32, C23, and C16) had no activity for at least 5 years before the 2012 mainshock, followed by target events that repeated two to three times in the subsequent 3 years (Fig. 2). One repeating family occurred once about 4 years before the mainshock and twice in the following 4 years. Assuming that the interseismic recurrence interval for families C32, C23, and C16 is longer than 5 years (the time span of pre-mainshock data), the recurrence time, tr, of the target events in these families decreases following the Nicoya mainshock, consistent with several other studies that report decreases in tr following large earthquakes (19, 20). C60 continues throughout the entire period and exhibits no systematic variation in recurrence interval. The sequences that persist longest after the Nicoya earthquake (C32, C55, and C60) occur in regions with long-term afterslip (slip ≥1.4 m), and the sequences that stop soonest are located in regions of little long-term afterslip (slip ≤0.4 m). We measure moment for the events in the sequences (see Materials and Methods and Fig. 2) and observe no systematic or substantial variation in Mo with tr.

Temporal variation in stress drop, rupture duration, slip, and area

We model the shape of the earthquake source spectra to estimate rupture duration and area and combine these measurements with the Mo to estimate slip and stress drop (∆σ) of the target events (stars in Fig. 2, A to E) in each recurring earthquake cluster (REC) (see Materials and Methods). In theory, calculating stress drop from the spectral shape should be relatively simple, but the intrinsic and complex nature of the earthquake source, combined with the difficulty of isolating it from attenuation effects, makes it a property that is hard to quantify robustly (8, 30, 31). Earthquakes that occur before the 2012 mainshock or after the postseismic slip have slowed (fault healed) and contain more power at higher frequencies than those from the time period immediately after the 2012 megathrust earthquake (fig. S1). A similar variation is seen in colocated smaller earthquakes. This variation could be revealing a change in source duration (and stress drop) or alternatively could simply reflect increased attenuation surrounding the interface following the large earthquake, as has been observed elsewhere (32). To separate the effects of attenuation and source, we follow an empirical Green’s function (EGF) approach to remove the path effects from the earthquake spectra.

Given the waveform similarity across the seismic network, indicative of identical focal mechanism and hypocenter location, as well as the magnitude distribution with respect to the corresponding target event (at least one order of magnitude smaller than the largest event), we consider all small-magnitude repeating events in a REC as EGFs (circles in Fig. 2, A to E). Source spectra for target events are obtained using the EGF approach [e.g., (8, 31)], taking advantage of the numerous EGFs in a cluster that are available for each target event. An example of the waveform deconvolution procedure is shown in fig. S2 and described in Materials and Methods.

We only use EGFs that occur within a few months before or after the target events (and both target and EGFs must all occur either before or after the megathrust event) to isolate temporal changes in the physical properties of the recurrent fault patches and the evolution of the strength of the fault zone, without introducing possible time variations in propagation properties of the medium. Although the seismic moments of the target events remain invariant throughout the time period of this study, we observe a notable increase in high-frequency content with recurrence interval in their spectral ratios. Figure 3 compares the earthquake source spectral ratios and source time functions (STFs) between the megathrust-affected target event (blue) in 2012 and the healed target event (green) in 2015 for three of the repeating families (C16, C60, and C32). Healed repeating events contain more power than the megathrust-affected events at frequencies between 6 and 18 Hz, where the spectral ratio is dominated by the shape of the larger earthquake and least affected by the smaller one.

Fig. 3 Temporal variability of RE source.

Earthquake source spectra and time function for three of the RE families: (A) C16, (B) C32, and (C) C60 obtained through the deconvolution procedure (fig. S2). Each panel shows the comparison between the megathrust-affected target event (blue) in 2012 and the healed target event (green) in 2015. The red and black dashed lines indicate the best-fit Boatwright model for each event, and the circles correspond to the corner frequency for each event. Note that for the three families, the healed (green) event shows more power at frequencies between 6 and 18 Hz, and its corner frequency is more than twice the corner frequency of the 2012 event. Avg., average.

Stress drop as a function of time for each target event within the five families studied here is shown in Fig. 4A. Our results indicate that target events experience a reduction in ∆σ after the Nicoya earthquake, followed by a systematic increase through 2016 when all repeaters cease. Consistent with many previous studies that have inferred a direct relationship between stress drop and fault material properties (33, 34), we interpret the temporal changes in stress drop as representing changes in frictional strength. The increase in frictional strength with tr for the five repeating families (Fig. 4B) is consistent with laboratory experiments that show that frictional healing increases logarithmically with tr (1, 2, 24).

Fig. 4 Stress drop variability as a function of time.

(A) Stress drop variability as a function of time with respect to the occurrence of the 5 September 2012, Mw 7.6 Nicoya Peninsula earthquake (black dashed line). Circles represent the target event for each of the RECs within the five families, color-coded by family id (see Materials and Methods for calculation of error bars). The gray bars mark periods of slow slip events from 2008 to 2017 (26), and the orange and green regions indicate the duration of the short- and long-term Nicoya earthquake afterslip (27, 28). Note that no repeaters were identified after 2016. (B) Stress drop for each family as a function of recurrence time (tr) from previous event (circles) and time delay of first event after the Mw 7.6 Nicoya Peninsula earthquake (stars). Note that stress drop increases with log (tr), consistent with laboratory experiments (2, 19, 20).

Like those of the target earthquakes, the spectra of the EGFs that occurred immediately after the Nicoya megathrust earthquake in 2012 contain less high-frequency power than those that happened 3 years after this event (fig. S3), consistent with a decrease in stress drop, an increase in attenuation, or both. In the deconvolution process, we use only EGFs from the same time period as the target events to ensure reliable correction for attenuation. If we use the same EGFs (from any chosen time period), or a stack from all time periods, to correct all the target events, then high-frequency power of the events most affected by the mainshock would be even lower, leading to a larger difference in stress drop between the two sets of events. Our results are consistent with previous studies showing similar stress drop variations following large transient events (8, 35, 36).


Frictional weakening and restrengthening

The reduced tr, exhibited by the target events in four of the five earthquake families directly after the 5 September 2012, Mw 7.6 Nicoya Peninsula megathrust earthquake, is consistent with previous observations (13, 18, 19) and modeling (3, 7), demonstrating that repeating sequences rupture more frequently after large crustal events. The reduced tr is coincident with a systematic decrease in ∆σ for all the RE families as predicted by rate and state friction laws (7) and observed in both laboratory experiments (13, 7) and RE sequences along the San Andreas fault system in California (1, 6, 24). We attribute this to an increased loading rate from afterslip hastening failure, reducing the healing time and the effective area of contact, Aeff (area of the fully locked region), within each rupture patch that consists of small, weak asperities. This interpretation is supported by our observation that RECs that persist longest after the Nicoya earthquake occur in regions with maximum long-term afterslip, whereas those sequences that stop soonest locate in regions of little long-term afterslip.

A weakened fault allows repeaters, driven by surrounding aseismic creep, to occur more frequently, because less accrued stress is required before seismic failure. Laboratory measurements and numerical simulations of earthquakes have demonstrated that at low interseismic rates, frictional strength and therefore ∆σ of faults show minor variations, whereas strong frictional weakening is observed at faster loading rates such as those during afterslip episodes (2, 37). As loading rate decreases, tr and healing time increase and so does fault strength and ∆σ. Our results indicate that ∆σ attains pre-mainshock values after ~1.5 years following the 2012 Nicoya Peninsula earthquake (Fig. 4). Variations in earthquake stress drop, which correlate with P-wave velocity reductions, along the Gofar transform fault (34) were interpreted as reflecting differences in the degree of fault zone damage and the ability of the fault to accumulate strain energy. A similar mechanism may be operating along the Nicoya megathrust in concert with frictional changes induced by afterslip.

Competing process of fault weakening and transient embrittlement

The target events off Nicoya repeat with nearly constant seismic moment, unaffected by changes in tr, during the 12 years analyzed (Fig. 2). Because Mo scales with rupture area and average slip, RE seismic moments that do not vary with tr imply a balance between competing processes that reduce slip (fault weakening) and increase the rupture area (transient embrittlement that increases rupture area by converting frictional properties of fault regions from aseismic to seismic), both induced by increased loading rate from the mainshock. We invoke the balance between these competing processes to explain the Nicoya RE moment invariance. Similar results were reported for REs following the 2004 Parkfield earthquake (30, 35). The trade-off between these two processes can explain all observed moment and recurrence time trends. A positive correlation between Mo and tr (with both decreasing) implies that slip reduction from the increased loading rate dominates, and little, if any, transient embrittlement occurs. A negative Mo and tr relationship (Mo increases as tr decreases) implies that transient embrittlement dominates over slip reduction. Figure 5 summarizes our interpretation with a schematic representation of the space-time evolution of the RE fault patches before and after the 5 September 2012, Mw 7.6 Nicoya Peninsula earthquake.

Fig. 5 Space-time evolution of the RE fault patches before and after the 5 September 2012, Mw 7.6 Nicoya Peninsula earthquake.

The figure shows the time-dependent behavior of frictional heterogeneities within the fault area of REs during the (A) late interseismic period (before the occurrence of the mainshock), (B) some days following the occurrence of the large event, and (C) ~1 year after the mainshock. The white regions of the fault patch represent areas with velocity-strengthening frictional properties, continuously creeping and elastically loading the edges of the unstable sliding asperities; the black areas with velocity-weakening (seismic) frictional properties. Following the occurrence of the mainshock in 2012, the slip area increases (from R1 to R2), but the effective area of contact, the black regions on the plate interface elastically coupled, reduces, decreasing the amount of slip in the subsequent rupture. As loading rate increases because of afterslip following the main event, asperities with conditionally stable sliding properties experience transient embrittlement, a transition from aseismic (light gray) to seismic (black) slip mode. After ~1 year following the Nicoya earthquake, afterslip amplitude reduces and some of the unstable areas (black) within the conditionally stable (dark gray) patches switch slip mode and behave aseismically again. Slip area reduces (from R2 to R3), but the effective area of contact increases.

During the late interseismic period (before the occurrence of the Nicoya mainshock), sections of the fault plane governed by velocity-strengthening (aseismic creep) frictional properties (regions in white) are steadily loading the asperities that otherwise behave seismically (shown in black). When loading exceeds the frictional strength of a single asperity or two, small earthquakes in a REC occur. Stress interactions between these small earthquakes, or similar failure thresholds that are exceeded for several neighboring asperities, eventually cause the target fault plane with radius R1 to fail in target events with recurrence time tr. Following the 2012 Nicoya Peninsula earthquake, increased loading rate from afterslip (i) reduced the effective area Aeff of the asperities; (ii) decreased their slip in subsequent rupture (4, 37, 38): (iii) converted regions with conditionally stable, aseismic frictional behavior (light gray) to seismic (black), resulting in an increase in the total rupture area from R1 to R2 (transient embrittlement); and (iv) reduced tr. As the increased loading rate from afterslip decayed, the fault plane gradually returned to its pre-mainshock configuration, asperity Aeff, slip, and tr increased while total rupture area decreased.


Our results represent the first well-documented evidence that these competing factors may be the mechanisms responsible for generating fault weakening along a subduction zone megathrust. After ~1 year following the Nicoya megathrust earthquake, afterslip amplitude reduced, and some of the frictionally unstable areas (black) within the conditionally stable (dark gray) patches switched slip mode and behaved aseismically again. Slip area reduced (from R2 to R3), but the effective area of contact increased, which we interpret as healing or frictional restrengthening. The time-dependent restrengthening of the unstable fault patches observed in Fig. 4 and schematically represented in Fig. 5 suggests that the plate interface healed within ~1.5 years after the 5 September 2012, Nicoya Peninsula earthquake, a time scale that is consistent with the period of largest afterslip (28).


Identification of REs

We use EQcorrscan (39), a template-matching detection algorithm in Python, to extend the published (July to December 2012) catalog of RE clusters generated for the aftershock sequence of the 5 September 2012, Mw 7.6, Nicoya Peninsula earthquake (18) to cover a period of about 12 years. Starting with 53 well-located repeating clusters in this catalog, we select those that contain events with at least one order of magnitude difference in earthquake size, reducing the catalog to 26 clusters. For each repeating cluster, we select the largest magnitude event (M2.2 to M3.30, termed target event) as a template to scan through continuous seismic data from 2007 to March 2018 to find matches. We use the vertical component at all available stations from the Nicoya Seismic Cycle Observatory (Fig. 1) to compute waveform cross-correlations. We apply a 1- to 15-Hz Butterworth band-pass filter to the data and use a 5-s time window starting 0.5 s before the manually picked P-wave arrival. All events with a network-averaged cross-correlation coefficient (CC) ≥ 0.95 are considered repeaters, and those with CC ≥ 0.85 are considered near repeaters. In total, we identify five earthquake clusters whose target events repeat at least twice between 2008 and 2016. Figures 1 and 2 show the spatial and temporal distribution of the repeating families C60, C55, C32, C23, and C16. Each of these families is formed by several clusters composed by a target event (represented by stars in Fig. 2) whose magnitude is constant or nearly constant with time and a minimum of three smaller events (at least one order of magnitude smaller than the largest event) that obey our definition of repeating or near-repeating events.

Calculation of seismic moment

Local magnitudes previously determined for the aftershock sequence of the Nicoya Peninsula earthquake (18) are used for computing the Mo for events that occurred in 2012. We use the singular value decomposition (SVD) approach (40) to compute relative moments for all the repeating events that occur in other years. This method takes advantage of the highly coherent waveforms of REs and yields more precise and accurate descriptions of earthquake size than standard catalog techniques. For the SVD analysis, we use the vertical component of all the seismic stations available for a given year that recorded the events. We low-pass filter seismograms with a Butterworth filter with a corner of 10 Hz to enhance the signal-to-noise ratio. Signals are aligned using a cross-correlation on the P-wave. Because of the proximity of the stations to the earthquake source, we use a total window length of 4.5 s that begins 0.5 s before the P-wave arrival.

Corner frequency determination

We analyze temporal changes in the corner frequency and stress drop of all target events in these RE clusters. We compute corner frequencies using the EGF spectral ratio approach (8, 30), summarized here. We classify the largest event of each repeating sequence as the target or main event and the smaller magnitude earthquakes within that REC as EGFs. Each EGF-target waveform pair is windowed using an empirical magnitude-dependent duration, that is of the order of 10 times the expected pulse duration of the target earthquake, assuming a constant stress drop scaling (31).

We calculate spectra and spectral ratios using a multitaper approach (41) that uses complex division, enabling the calculation of the mainshock STF. Deconvolution of a clear source pulse confirms the EGF assumption, because the approximation is good enough to work in both phase and amplitude. Time domain source functions or STFs are also useful for discriminating complex earthquakes that cannot be fully described by a simple source model (fig. S1). We fit each individual spectral ratio to calculate corner frequencies using the spectral source model (42, 43)Ṁ1(f)Ṁ2(f)=M01M02[1+(f/fc2)γn1+(f/fc1)γn]1γ(1)where f is the frequency, fc1 and fc2 are the corner frequencies of the target and the EGF, respectively, n is the frequency falloff (here, we assume n = 2), and γ is a constant controlling the shape of the corner. M01 and M02 correspond with the seismic moment of the target and the EGF, respectively. We try fitting the observed spectral ratio to the original Brune model (γ = 1) and a sharper corner or Boatwright model (γ = 2) and determined that, overall, the shaper corner provided a better fit to the majority of the spectral ratios analyzed. Before fitting, we log-sample the individual spectral ratios resulting from the deconvolution of each EGF-target earthquake pair for each station in the frequency domain to reduce the weighting toward the higher-frequency part of the spectra. We use a log10(∆f) sample rate of 0.05 and lastly select the frequency range where the signal, defined by the P-wave train, with length nsec, is at least three times the noise, defined by the pre–P-wave arrival with the same length of the signal. We average the resulting individual spectral ratios and STFs to obtain the mainshock spectrum and STF. Averaging across stations at multiple azimuths improves the stability and resolution of the corner frequency measurements and reduces the possibility of underestimating or overestimating the mainshock source duration (1/fc1) due to possible finiteness effects of the source (e.g., earthquake directivity). We fit the average source spectrum using the Levenberg-Marquardt algorithm in Non-Linear Least-Squares Minimization and Curve-Fitting for Python (LMFIT) (, a Python, open source library for nonlinear curve fitting and perform a gridsearch around the resulting value of fc1 to determine the range of fc1 and fc2 in which the fit between the data and the model given by the variance reduction (VR) is at least 85%VR=[1(datamodel)2(data)2]×100(2)

The gridsearch limits around the fc1 are defined as 2 × F, where F (F test) represents the confidence interval (uncertainty), and it is defined asF(Pfix,NP)=(χf2χ021)×NPPfix(3)where χ0 is the null model (which is the best fit we have found), χf is an alternate model where one of the parameters is fixed to an specific value, N is the number of data points, P is the number of parameters of the null model, and Pfix is the number of fixed parameters. To limit the effects of complexities in the source spectra, we use only measurements where the difference in amplitude of the high- and low-frequency levels in the fit is greater than 2. We do not use estimates of fc2 because for some cases, they are typically out of the available bandwidth.

Stress drop determination

To calculate stress drop (∆σ) from the previously determined seismic moment and corner frequency measurements, we assume a simple circular model (43)σ=716 Mo fc3k3 β3(4)where Mo is the seismic moment, β is the shear wave velocity (β = ~3.5 km/s), and k is a constant equal to 0.32 for P-waves (44). The estimated uncertainty for the mainshock corner frequency is used in Eq. 4 for estimating the stress drop SD, presented as error bars in Fig. 3. This widely used model is known to be an oversimplification (8), and the absolute values of stress drop estimated in this way are subject to large uncertainties, both random and systematic (30, 45). Application of this approach is most stable when applied to repeating sequences (8), and the relative values in our detailed analysis reflect real variation in the frequency content of the energy radiated by the earthquakes. Any temporal changes in the velocity, rigidity, or attenuation in the immediate source region caused by the occurrence of the megathrust earthquake may also contribute to the estimates of Mo, β, and ∆σ affecting the uncertainties in our final results. Such changes are indistinguishable from changes in earthquake source parameters given the spatial resolution and the source dimensions of the REs (M < 3.0) analyzed. Although current observations of postseismic decrease in velocity are only a few percent (46), they could be higher within a low-velocity fault zone in the immediate source region (34). This may be enough to affect the behavior (D¯,Aeff) and evolution of RE fault patches during the earthquake cycle, adding to the temporal variation in source parameters we observe.


Supplementary material for this article is available at

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 M. Protti and V. González from the OVSICORI-UNA (Observatorio Vulcanológico y Sismológico de Costa Rica Universidad Nacional) for installing and operating the Nicoya seismic network with critical support from A. Newman, D. Sampson, and Z. Peng. E.J.C. acknowledges support from all the scientists, engineers and administration team of OVSICORI-UNA, Consejo Nacional para Investigaciones Científicas y Tecnológicas de Costa Rica (CONICIT) and Ministerio de Ciencia, Tecnología y Telecomunicaciones de Costa Rica (MICITT). Discussion and comments from E. Brodsky, T. Lay, C. Cattania, and W. Frank greatly improved the manuscript. Funding: S.Y.S. and R.E.A. acknowledge support from NSF awards EAR-1551683 and EAR-1551758, respectively. E.J.C. acknowledges support from CONICIT award FI-0129-13. Author contributions: E.J.C. performed all the computations. E.J.C., S.Y.S., and R.E.A. analyzed the data, interpreted the results, and wrote the paper. S.Y.S. and R.E.A. contributed equally to the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All seismic data from the Nicoya Seismic Cycle Observatory, including all stations used in this work, are freely available at the Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC;, last accessed in April 2018). 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. Figures were generated using open source libraries written in Python, Generic Mapping Tools (GMT) and through the seismological community scientific library ObsPy.

Stay Connected to Science Advances

Navigate This Article