The hidden simplicity of subduction megathrust earthquakes

See allHide authors and affiliations

Science  22 Sep 2017:
Vol. 357, Issue 6357, pp. 1277-1281
DOI: 10.1126/science.aan5643

Universal scaling for big quakes

The amount of energy released as a large fault ruptures provides some clues about the overall size of an earthquake. Meier et al. looked at this energy release for more than 100 large earthquakes. Although the overall size of an earthquake cannot be predicted from the rate of energy release, a minimum size can be estimated. Estimating this minimum size could add valuable seconds to early earthquake warning algorithms, helping to avoid damage and save people from injury or death.

Science, this issue p. 1277


The largest observed earthquakes occur on subduction interfaces and frequently cause widespread damage and loss of life. Understanding the rupture behavior of megathrust events is crucial for earthquake rupture physics, as well as for earthquake early-warning systems. However, the large variability in behavior between individual events seemingly defies a description with a simple unifying model. Here we use three source time function (STF) data sets for subduction zone earthquakes, with moment magnitude Mw ≥ 7, and show that such large ruptures share a typical universal behavior. The median STF is scalable between events with different sizes, grows linearly, and is nearly triangular. The deviations from the median behavior are multiplicative and Gaussian—that is, they are proportionally larger for larger events. Our observations suggest that earthquake magnitudes cannot be predicted from the characteristics of rupture onsets.

Advances in earthquake research are often sparked by observations of “peculiarities” of individual earthquakes [e.g., (13)], implying that there is a “normal” or average earthquake behavior from which such events deviate. This expected normal behavior, however, is not well defined from an observational perspective; more often, our expectations are driven by conceptual rupture models—for example, that of self-similar cracks (4) or pulses (5)—that may or may not be supported by observational evidence. For large earthquakes in particular, whose rupture processes are obviously complex, it is questionable whether such simple models are meaningful.

Having a realistic and data-driven model for rupture evolution is important beyond earthquake source physics because the question of rupture predictability hinges on it. Although there is virtually unanimous consensus that the occurrence of earthquakes currently cannot be predicted (6), numerous studies have suggested that the evolution of earthquake ruptures themselves may follow a predictable pattern. These studies suggest that the characteristics of the rupture onset can be used to predict the final size of a rupture. In particular, it has been suggested that ruptures that start impulsively may be more likely to grow into large events (7) and, conversely, that larger ruptures start more slowly, with a longer (and seismically observable) nucleation phase (8). However, observational reports on correlations between rupture-onset characteristics and final magnitudes (9, 10) have been controversial (11, 12), and there is strong observational evidence that small and large shallow continental crust earthquakes start in indistinguishable ways (12).

Earthquake early warning (EEW) algorithms (1315) could provide substantially longer warning times if future rupture evolution could indeed be predicted. Longer warning times would allow alert recipients to perform more or better real-time damage mitigation procedures, such as shutting down gas lines and slowing trains (16). In the absence of a rigorous statistical description for the typical rupture evolution, however, it remains unclear whether this scenario is realistic.

Here we attempt to accurately describe the typical temporal rupture behavior of large subduction earthquakes, using three extensive source time function (STF) catalogs. The primary catalog considered here was compiled by Ye et al. (17), who used a classical kinematic teleseismic source inversion scheme (18, 19) to infer finite source models for 116 of the largest subduction zone earthquakes of the past three decades (Fig. 1A). Our results are further confirmed by analyzing two independent STF catalogs (20, 21).

Fig. 1 STF data set and magnitude binning.

(A) The 116 STFs from the kinematic source models of Ye et al. (17) colored by moment magnitude, Mw. (B) A Shapiro-Wilk (SW) normality test for absolute (black) and log-transformed (red) amplitudes shows that at most points in time, the STF amplitudes are consistent with a log-normal distribution. The P value for the amplitude values on a linear scale (black line) are below 0.1% for almost the entire time range, indicating that the amplitude values are grossly inconsistent with a normal distribution. (C) Magnitude binning, where in each bin we choose the 20 nearest neighbors in terms of Mw—i.e., the 20 STFs whose Mw is closest to the target magnitude (red line). For the larger-magnitude bins, data scarcity leads to saturation, and events from a large size range are included in the same bin. Median magnitude in each bin is given by the blue line.

To identify the typical temporal behavior of earthquakes of different size ranges, we grouped the STFs by magnitude with a nearest-neighbor approach. For seven target moment magnitudes Mw = 7.0, 7.25, 7.5, ..., 8.5, we found the STFs from the 20 events with most similar magnitudes (Fig. 1C) and, for each point in time, computed the median STF amplitude in each bin (Fig. 2A). The STF amplitudes are log-normally distributed according to a Shapiro-Wilk normality test (Fig. 1B). We chose the median rather than the arithmetic mean because using arithmetic means on log-normally distributed data may lead to a biased interpretation [compare synthetic test (19)].

Fig. 2 Median STFs in seven magnitude bins and for the 12 STFs with Mw ≥ 8.

(A) Median of STFs on physical scale and (B) median of STFs normalized by twice the centroid duration. The normalization is done for each individual STF separately, rather than on the median STF. The approximate range of when the median normalized STFs reach peak amplitudes is indicated in gray (35 to 55% of full-rupture duration). (C and D) Same as (B) but for two independent STF data sets (20, 21); see (19) and fig. S3 for details. All three data sets have the same near-triangular and scalable median shape. The median STFs on the physical scale are accurately approximated with the functional form yfit = μtexp[–1/2(λt)2] (dashed lines). Diamonds in (A) show centroid times for fitted functions. A larger version of (A) is shown in fig. S11.

Despite the large variability between individual STFs, the median behavior in the seven magnitude bins is fairly regular, with near-triangular and almost symmetrical shapes (Fig. 2A). The STFs grow linearly until the peak moment rate is reached and then start decaying with a similar rate. The average growth rate is ~1.25 × 1018 Nm/s2 and does not vary systematically with magnitude. Only the smallest-magnitude bin has a slightly lower initial slope, although this could be an artifact from the low source inversion resolution during the first 10 s (19). The medians for larger magnitudes are less regular than those of smaller magnitudes because the high-magnitude bins combine a much wider range of rupture sizes (compare Fig. 1C). We provide a comparison of each individual STF with the median STF as supplementary material (19).

To ascertain that the features we observed are not biased by the properties of the kinematic inversion procedure, we compare the results of this study with two alternative STF data sets (19) that were generated with a different slip parameterization (20) and a synthetic Green’s function deconvolution approach (21). This analysis demonstrates that the STFs are generally well constrained and change only in terms of second-order aspects across different methods. Furthermore, we performed synthetic tests regarding the slip parameterization, using a self-similar pulse (5) as input model (19). We found that the temporal rupture evolution, and in particular the observed linear rupture growth, is not an artifact of the slip parameterization assumed in the source inversion procedure by (17).

We approximated the median STF with the functional form yfit = μtexp[–1/2(λt)ϕ], with ϕ = 2 (Fig. 2A). This form corresponds to the prediction for the STF shape from discrete avalanche models that have been applied to earthquake ruptures (2225) and has been used to describe ground motion envelopes (26). The form essentially describes a triangle with a rounded hat and a short tail. Setting ϕ = 1 yields the classical Brune’s model (27), which gives rise to the often-reported Brune ω−2 displacement amplitude spectrum, but which has a poor fit with the median STF (fig. S11). The form we used with ϕ = 2 also has a ω−2 fall-off rate, but the spectral shape is more complicated (fig. S11C) and does not have a simple analytical expression.

Owing to the sparsity of high-magnitude events, none of the median curves accurately represent the events with Mw ≥ 8 (compare Fig. 1C). The median STF of the 12 events with Mw ≥ 8 breaks away from the average moment rate growth trajectory at ~12 s and does not appear to follow the regular shape of the smaller events. All of the Mw ≥ 8 events have extended periods in which they grow at a rate greater than four times the median rate of ~1.25 × 1018 Nm/s2 (fig. S7). This raises the question of whether these largest events differ from the smaller ones in terms of their physical rupture processes.

To facilitate comparison of the STFs across the large range of event sizes, we normalize the STFs with respect to their rupture durations. Rupture durations, however, are often poorly constrained because, at a late rupture stage, the direct waveform phases used to infer the STF overlap with indirect phases that the source model does not capture (17, 18). Deciding which phases should be modeled and which phases are secondary is difficult and involves subjective judgment. Consequently, STFs often include badly constrained low-amplitude tails (19). Using the time when the STF reaches zero may therefore lead to a systematic overestimate of rupture durations (28). Centroid times, by contrast, defined as Embedded Image, where Embedded Image is the moment rate, quantify the first temporal statistical moment of an STF and provide an objective estimate of rupture half-durations. Centroid times are much less affected by poorly constrained STF tails and subjective modeling choices (28) and are reported on a regular basis for large earthquakes (29). For the STF normalization, we therefore use centroid times.

We divided the time of each individual STF by 2Tc and multiplied the STF amplitudes by a constant such that the normalized STF integrates to 1. We then plotted the median of these normalized STFs in the seven magnitude bins (Fig. 2B). The median normalized STFs show a high degree of similarity. If we take the average of the seven median STFs as a reference, the root mean square deviation from this reference lies between 0.06 and 0.11 for the seven bins on the normalized amplitude scale—i.e., the average deviation from the reference is less than 5% of its peak value. The typical STFs are fairly symmetric in shape, and peak moment rates are reached between 35 and 55% of the rupture duration. This relatively broad range may in part be a consequence of the large uncertainties in rupture duration estimates.

The high similarity of the median normalized STF suggests an underlying general pattern that is typically shared by large earthquakes, despite the large variability of individual earthquakes. The general pattern is near-triangular with roughly linear moment rate growth and decay. The Mw ≥ 8 events no longer appear anomalous when normalized. Their longer rupture durations compensate for their higher amplitudes, and they follow the same shape as the smaller events, with linear moment rate evolution. That all median STFs have the same normalized shape suggests that the larger events can be represented as scaled-up versions of the smaller events. However, this scalability is not equivalent to the classical concept of earthquake similarity in the sense of Aki 1967 (4), which requires that the ratio of slip over length be constant. STFs do not constrain this ratio because they are proportional to the product of slip and rupture area. Classical similarity predicts a quadratic STF onset, which is at odds with the linear growth found here. In earthquake seismology, the similarity conditions of (4) are often referred to as “self-similar,” even though this designation has a different meaning in other fields of physics. Here we use “self-similar” to mean similar in the sense of (4).

To quantify the deviations from the typical STF shape, we fitted the functional form yfit = μtexp[–1/2(λt)2] to each individual STF, by optimizing μ and λ in a least-squares sense (Fig. 3). The μ values, which quantify the initial slopes of the STFs, vary strongly between events but do not show a systematic trend with magnitude (Fig. 3A). The lower values for the smallest events are likely affected by resolution limits of the inversion method. This suggests that rupture-onset characteristics are not diagnostic of final event sizes. The same observation is made with the two independent data sets used in the robustness analysis (fig. S4). The characteristic time scale of the ruptures is quantified by 1/λ and shows a power-law dependency with seismic moment (Fig. 3B).

Fig. 3 Fitted parameters and source duration estimates for the 116 individual STFs (black dots) and for the median STFs (diamonds).

Vertical bars for median values give 5th and 95th percentiles from 1000 bootstrap resamplings. (A) Initial slope μ versus seismic moment M0; solid line gives median μ value of fits to individual STFs. (B) The inverse characteristic time scale, λ, versus M0, with corresponding least-squares fit (solid line). (C) Centroid time versus seismic moment with least-squares fit to the durations of the median STFs (solid line). Dashed line gives best fit with imposed slope of 1/3 for reference.

Source durations have been widely reported to scale with seismic moment as T ~ M01/3 for a large range of earthquake sizes, and this scaling is well established, with abundant data for earthquakes with M ≲ 7 (4, 3034). However, the observed linear and magnitude-independent moment rate growth that we observe implies a T ~ M01/2 scaling, because the growth slope is proportional to ~M0/T2. Our observations, therefore, suggest a scaling break somewhere between the smaller earthquakes for which the T ~ M01/3 scaling is well constrained and the M > 7 events analyzed in this study.

Testing the hypothesis of such a scaling break is difficult, because seismological duration estimates usually do not cover a sufficiently wide magnitude range (19, 29) (fig. S9), or because they have been made in the frequency domain and are not necessarily consistent with the time-domain estimates discussed here (19, 33, 35) (fig. S10). Given the large scatter of the duration observations, the covered Mw range of our data set is not large enough to confirm or rule out either a M01/3 or M01/2 scaling (fig. S8). Both trends are consistent with the data, with coefficients of determination (R2) of 0.76 and 0.70, respectively. A study of observed source spectra from 942 thrust earthquakes with Mw ≥ 5.5 recently reported the same scaling break, from M01/3 to M01/2 at Mw ~ 7.0, but it had to use two different methods for small and large earthquakes (36). A scaling break to T ~ M0 was reported for crustal events with narrow seismogenic width (<30 km) (37), based on a database of finite-fault source models (38). In the present study of subduction earthquakes with much larger seismogenic width, the scaling break is a consequence of the typical STF shape.

The fits of the functional form yfit = μtexp[–1/2(λt)2] to the individual STFs allow us to analyze quantitatively if the Mw ≥ 8 events deviate from the general pattern of the smaller events. The residual ratios, Δyr (t′) = yobs (t′)/yfit (t′), as a function of normalized time t′ = t/(2Tc), are very similar in different-magnitude bins, and their standard deviation does not vary systematically as a function of magnitude (Fig. 4A). The empirical residual distributions from different-magnitude bins are nearly identical and conform with a normal distribution (Fig. 4D) with a standard deviation of σ = 0.38 (Fig. 4D). Here, we have excluded three outlier events that feature multiple episodes of strong moment rate growth and therefore are not well described with the simple fitted model (fig. S12).

Fig. 4 Residual ratios between fitted and observed STF.

(A) Residual ratios, ∆yr (t′) = yobs (t′)/yfit (t′), for the 116 individual STFs in the seven magnitude bins. (B) Shapiro-Wilk normality test for the amplitude distributions of the 116 residual ratios at each time step. Except at very early and very late times, the residual ratios conform with a normal distribution. (C) Example STF illustrating how residual ratios are calculated. (D) Concatenated residual ratios between 10 and 90% of t′ in each magnitude bin, and normal distribution with standard deviation of 0.38 for reference (red dashed line). Near the start and end of the rupture, when either yobs or yfit are close to zero, the residual ratios may be unstable.

The residual ratios represent higher-frequency oscillations around the smooth model, yfit. These moment accelerations and decelerations create seismic radiation throughout the rupture process and contain a substantial fraction of the total signal energy. The amplitude spectra of the smooth models saturate at higher frequencies in that all magnitudes have similar amplitudes at frequencies above ~0.03 Hz (fig. S13B). This stands in contrast to the amplitude spectra of the observed STFs, yobs, in which the energy content in this frequency range increases with magnitude (fig. S13A), and demonstrates that the observed source spectra at higher frequencies are strongly affected by the oscillations around the smooth models—i.e., by the irregularities of the rupture process. This result is noteworthy, because commonly used simple source models assume that high-frequency radiation is only generated during rupture initiation and termination [e.g., (27, 39)]. The residual ratios are characterized by a distinct 1/ω decay in terms of their Fourier amplitude spectra (fig. S13C).

The similarity of the residual ratios for all earthquake sizes suggests that all STFs, including the Mw ≥ 8 events, are compatible with the same universal base model for moment rate evolution, with normally distributed interevent variability superimposed multiplicatively, yobs(t′) = yfit(t′) × [1 + ε(t′)], where ε ~ N(0,0.382). That the deviations from the base model are multiplicative, rather than additive, implies that the deviations are proportional to the absolute STF amplitudes. Consequently, the Mw ≥ 8 events feature the largest absolute deviations from the base model. This explains why, for these events, a simple base model seems least appropriate if individual STFs are compared at face value, and the largest events may appear exceptional. The observed multiplicative nature of the residuals, furthermore, is likely of physical importance. The multiplicative excursions from the simple model may require a dynamical model in which fluctuations of the rupture process are exacerbated proportionally to the current rupture size. One explanation could be that larger ruptures, with larger energy-release rates, can rupture further into zones with unfavorable rupture conditions (e.g., areas of low stress), giving rise to larger deviations.

The typical STF shape that we observed challenges standard models of earthquake rupture. In particular, self-similar pulse and crack models predict that moment rate onset evolves as ~tη, with η = 2 [e.g., (5, 40, 41)], in contrast to the linear growth that we observed. The implications for rupture models can be understood through simple scaling relations. For a rupture with area growing as Embedded Image and average slip growing as Embedded Image, seismic moment Embedded Image scales as Embedded Image and moment rate grows with an exponent of Embedded Image. The exponent that we observed is Embedded Image such that Embedded Image. This constraint is incompatible with standard self-similar rupture models, which have Embedded Image and Embedded Image. One admissible end-member model is a circular pulse growing with constant speed and constant slip (i.e., Embedded Image and Embedded Image). Such behavior emerges in steady-state pulse models (42, 43) and in models of rupture over disconnected slip patches with areas of systematic slip deficits (2225, 39). If slip is not constant (i.e., Embedded Image), the rupturing area must grow more slowly than predicted by these models (Embedded Image). An admissible model with Embedded Image is that of a rupture expanding in one predominant direction [e.g., (36, 44)] with linearly increasing slip (i.e., Embedded Image and Embedded Image). Combined with independent scaling relations (e.g., between moment and slip), the condition Embedded Image derived here from STFs is a valuable constraint on rupture models.

Further insight into the deviation from self-similar rupture can be gained by combining our observations with peak ground motion displacement observations from short-distance recordings of shallow crustal earthquakes (12). Peak displacements are a proxy for moment rate and were observed to be consistent with self-similar two-dimensional growth (η = 2) for ≤1 s after rupture onset, after which growth was substantially slower. The shallow crustal earthquakes also had the same magnitude-independent onset behavior in that small and large earthquakes grow with equal rates at first until the smaller ones branch out as they approach their peak moment rates. The STFs analyzed here, which become reliable only after several seconds (19), have a well-resolved linear trend (η = 1). This suggests that the ruptures at some point transition from initially self-similar to systematically slower growth. This transition must occur before the linear trend is well constrained, at source dimensions much smaller than the seismogenic width, which is ~200 km for subduction earthquakes (45). Therefore, whereas small earthquakes with durations <1 s (i.e., M ≲ 5.0) may be accurately described as self-similar ruptures, the nature of larger ruptures seems to change after they have reached a critical size that is not controlled by seismogenic width. The change could occur, e.g., because plausible weakening mechanisms, such as fault melting or thermal pressurization, need a minimum energy level to be triggered.

The multiplicative nature of the STF residuals has direct consequences for long-period (T > 1 s) strong earthquake ground motion. Large events with high moment rates have proportionally high moment rate fluctuations around the smooth model. Moment rates are directly proportional to long-period ground motion, which explains why ground motion variability at such periods grows strongly with magnitude, when measured on a linear scale. In ground motion prediction equations, this increasing variability is commonly captured with a log-normal residual description [e.g., (46, 47)]. The multiplicative nature of STF residuals implies that large-magnitude events can feature extreme values of long-period ground motion, not as outliers, but as random realizations of a multiplicative process.

The shape that we observed for typical moment rate evolution precludes strong rupture predictability as smaller and larger earthquakes are statistically indistinguishable until they reach their peak moment rate. This result has implications for EEW. Rupture-onset observations cannot diagnose the final rupture size. However, the fairly symmetrical nature and the gradual decay of STFs imply a weak rupture predictability: At any point during moment rate growth, the final seismic moment is statistically at least twice as large as the moment already generated. This makes a large difference for large ruptures. For instance, an earthquake that reaches Mw = 7.2 during its growing phase will at least double its size to Mw = 7.4. This strongly increases the area experiencing dangerous ground motion. Estimating moment rates in real time with EEW algorithms could therefore substantially lengthen warning times.

Individual source time functions of large shallow subduction zone earthquakes are often complicated and have a diverse range of forms, reflecting complex physical rupture processes. Despite this large diversity, however, there is a well-defined temporal pattern that such earthquakes typically follow. This pattern attests to the universality of the physical mechanisms that dominate rupture process in such earthquakes. The simple base model and the multiplicative deviations from it constitute a comprehensive characterization of first- and second-order universal features of the macroscopic behavior of large earthquakes. They provide powerful observational constraints for our attempts to understand the physics of earthquake rupture.


Supplementary Text

Figs. S1 to S13

References (4850)

References and Notes

  1. Materials and methods are available as supplementary materials.
  2. Acknowledgments: The STF data used in this study were generated by Ye et al. (2016), and can be accessed from The finite fault models and STFs from Hayes (2017) can be obtained through the U.S. Geological Survey National Earthquake Information Center Combined Catalog ( The STFs from Vallée et al. (2011) can be downloaded from This study has been partially funded by the Swiss National Science Foundation and the Gordon and Betty Moore Foundation. J.P.A. acknowledges funding from NSF CAREER award EAR-1151926. We thank L. Ye, H. Kanamori, K. Dahmen, N. Lapusta, H. Owhadi, L. Rivera, D. Sornette, and J.-P. Avouac for discussions, and M. Vallée and G. Hayes for providing their STF data sets.
View Abstract

Navigate This Article