Research ArticleGEOPHYSICS

Characterizing large earthquakes before rupture is complete

See allHide authors and affiliations

Science Advances  29 May 2019:
Vol. 5, no. 5, eaav2032
DOI: 10.1126/sciadv.aav2032


Whether earthquakes of different sizes are distinguishable early in their rupture process is a subject of debate. Studies have shown that the frequency content of radiated seismic energy in the first seconds of earthquakes scales with magnitude, implying determinism. Other studies have shown that recordings of ground displacement from small to moderate-sized earthquakes are indistinguishable, implying a universal early rupture process. Regardless of how earthquakes start, events of different sizes must be distinguishable at some point. If that difference occurs before the rupture duration of the smaller event, this implies some level of determinism. We show through analysis of a database of source time functions and near-source displacement records that, after an initiation phase, ruptures of M7 to M9 earthquakes organize into a slip pulse, the kinematic properties of which scale with magnitude. Hence, early in the rupture process—after about 10 s—large and very large earthquakes can be distinguished.


An outstanding question in the field of earthquake dynamics and hazards is whether there is any determinism in the rupture process of large events. That is, when, within a minute-long rupture process, is a very large earthquake (~M9) distinguishable from an only large one (~M8)? The existence of any determinism (or lack thereof) is important for understanding the physics that govern rupture. Furthermore, it will define the minimum possible time at which characterization of an event and its resulting hazards (strong shaking or tsunami) can be accurately inferred.

One end member view poses that the characteristics of rupture onset at nucleation or shortly thereafter (1 to 3 s) can be used, in a predictive sense, to estimate the final size (the magnitude) of an event. Observational evidence from seismic data of this strong determinism has been proffered for events in the M3 to M8 range (13). These results are surprising considering that rupture durations for intermediate- to large-magnitude events (~M6 to ~M8) can be tens of seconds to several minutes. Studies supporting the contrasting view that there is no determinism whatsoever in rupture onsets use the same seismic observations (45). This paradigm suggests that one must simply wait until rupture has progressed a substantial amount or even until it is complete (67) to ascertain an earthquake’s magnitude and dimensions and to characterize the associated hazards.

In the past, the datasets best suited to addressing the question of how rapidly we can know how big an earthquake will become have been near-field seismic recordings of moderate to large events. Today, however, insights on rupture dynamics can be gained from other observations as well. Geodetic measurements from high-rate Global Navigation Satellite Systems (HR-GNSS), which are now becoming commonplace, show that, for M7+ events, if stations are close to a rupture, then it is generally possible to determine the final magnitude of an event before rupture is complete (8). Similarly, systematic studies of teleseismic slip inversions (6, 9) find that moment release patterns are distinguishable at some point before the end of rupture. The overarching question is when. As a result of these studies, a third view on determinism has emerged: It has been hypothesized that there exists a weak determinism (6, 9, 10); in this model, the first 1 to 3 s of rupture are too similar between events of different final magnitude for any measurable differences to be observed, but at some point, following nucleation (10 to 20 s), rupture organizes into a steady-state slip pulse (11) whose characteristics (average rise time and pulse width) are diagnostic of an earthquake’s final magnitude (9). However, evidence of that model has been scant; here, we will provide two very strong lines of evidence in support of this weakly deterministic model of large earthquake ruptures.

First, we analyze the rupture properties of large earthquakes using source time functions (STFs) (Fig. 1)—these depict the moment release of an earthquake as a function of time—from several independent databases. The primary STF data used below come from the analysis of moderate to large earthquakes using a synthetic Green’s function deconvolution approach (12). This database contains STFs for more than 3000 earthquakes of M6 and larger since the early 1990s and, as such, provides a rich dataset with which we can search for systematic properties of earthquake rupture growth. We confirm the properties revealed in these data using two finite fault databases (13, 14), which contain hundreds of rupture histories and associated STFs for large earthquakes (M7+) since about 1990, many of them for events common to the larger dataset.

Fig. 1 The calculation of STF moment acceleration, using the 2010 Maule, Chile M8.8 earthquake as an example.

Moment acceleration is calculated as the average slope of the STF over various time spans: 2 s (green), 5 s (cyan), 10 s (orange), and 20 s (red). The full STF for this event is shown in the inset.


To assess the early rupture properties of earthquakes, we calculate the “moment acceleration,” i.e., the rate of growth of the STF over the first 1 to 20 s following the earthquake origin time (Fig. 1). To focus on the broad early characteristics of the STF rather than the variable nature of the function at a specific time, we average moment acceleration from the beginning of the STF to the time in question. We find that individual events of a given magnitude demonstrate substantial variability (Fig. 2), a reflection of natural differences between earthquake ruptures. However, their median behavior over a range of magnitude bins shows systematic differences that scale with magnitude. On average, earthquakes with a larger final magnitude grow faster early on in the source process.

Fig. 2 Seismic observations of weak determinism from a database of STFs (12).

The moment acceleration of earthquake STFs (gray circles) at 2 s (A), 5 s (B), 10 s (C), and 20 s (D). Stars represent median values for each 0.2 magnitude-unit bin. Light gray dashed lines represent the theoretical limit to the data, given by M0/T2, where T is the time window in question. Dark gray dashed lines in the bottom right of each figure show the power “n” for a variety of scaling relations given by M¨0M0n, where n varies from 1, to 1/2, to 1/3, as shown in (A).

In the very early stages of rupture (at 2 s after origin time (OT); Fig. 2A and figs. S1 to S3), there is only a weak correlation between moment acceleration and magnitude. In other words, very early in the rupture process, minimal or no scaling is observed and earthquakes of different sizes cannot be distinguished from one another. However, as time goes on (at 5, 10, and 20 s after OT; Fig. 2, B to D), a clear log-linear relationship develops between moment acceleration and magnitude. The implication of this is that, at some short interval after OT (5 to 10 s), well before these ruptures are complete, large events can be distinguished from even larger earthquakes.

We believe that this change from the early stages of rupture (1 to 2 s) to the time period at which events are distinguishable (5 to 10 s) represents the time it takes for ruptures to organize into a slip pulse (11) with properties that are different, even this early on, depending on the final size of that event (9). The larger the event is, the broader the slip pulse, the larger the area of the fault rupturing at any one time, and the greater the rate at which slip occurs (i.e., the greater the moment acceleration).

When the slip pulse model was first proposed (11), it predicted that displacement records close to large earthquakes should show a rapid growth, following the initial P wave, to the peak ground displacement (PGD). Furthermore, if the slip pulse properties are diagnostic of final magnitude, then the PGD should always occur in a time much shorter than the source duration and its growth rate and the final amplitude should be indicative of the final magnitude (9). However, broadband ground displacement records from HR-GNSS close to the source of large events have, until recently, not been systematically collected or made available.

We study displacement records from an HR-GNSS database of M6.0 to M9.1 earthquakes (15). We pay special attention to recordings in that database in the near-source region of large events (Fig. 3 and fig. S4). HR-GNSS is not sensitive enough to record the earliest onsets of small-amplitude P waves (typical noise levels are in the 1- to 5-cm range), so one cannot tell when the record begins. In addition, the point at which the displacement grows above the HR-GNSS noise will depend on an event’s size and on the radiation pattern and distance to the source. To ameliorate this, we use P-wave arrival times at more sensitive seismic sensors and interpolate the P-wave arrival times to the locations of the HR-GNSS sites. Thus, a further constraint in the event selection is that a good-quality seismic network be available as well. We find that 12 earthquakes in the M7.1 to M9.1 range in Chile, Costa Rica, Ecuador, and Japan have enough usable data for analysis. Because we will study the early development of the displacement waveforms, it is important to understand any errors associated with the P-wave arrival time interpolation procedure. A bootstrap analysis shows that the median error in predicting the onset times at the HR-GNSS sites is 0.5 to 1.5 s depending on the earthquake (figs. S5 to S6).

Fig. 3 Geodetic observations of weak determinism.

(A) Sample waveform from station CNBA (67 km from the hypocenter) during the M8.3 Illapel earthquake. Shown are the three components of motion and the total displacement (red) and the time evolution of PGD (PGDt) waveform (dashed black). (B) PGDt waveform for 12 events colored by distance to the hypocenter. The thick black line represents the cumulative moment from the STF [database of Hayes (2017) (13) for completeness purposes] for the event normalized for presentation. Gray areas are the noise level of HR-GNSS within which the data cannot be reliably interpreted. STFs for two events, Ibaraki and Iwate, are not available because the earthquakes occurred in quick succession after the Tohoku-oki mainshock.

We analyze the time evolution of PGD (PGDt; Fig. 3). PGDt is the maximum value of displacement, at any given instant, observed since the inferred P-wave onset. It is obtained from the total displacement, which is the epoch-by-epoch Cartesian norm of the three components of motion (Fig. 3). We compare the growth of PGDt to the temporal evolution of the moment release as measured by the STFs. For these events, we find that the near-source records finish growing to their final value (or very close to it) in a time much shorter than the event duration. This behavior is a strong function of distance to the source. While the records will be affected by complexities in the STFs, such as multiple moment rate peaks, at close distances, and at long periods, HR-GNSS records will be mostly unaffected by elastic wave propagation and are dominated by the ramp-like growth to PGD and to the static offset. Thus, they more accurately reflect the behavior of slip at adjacent portions of the fault (10). At greater distances, this relationship breaks down and HR-GNSS records reflect mostly S waves and surface waves.

We further study whether events of different magnitudes exhibit systematically different early PGDt growth patterns (Fig. 4). We stack waveforms in four magnitude bins (M7.1 to M7.4, M7.6 to M7.8, M8.1 to M8.3, and M8.9 to M9.1) and two hypocentral distance bins (30 to 150 km and 150 to 250 km). Before this, we apply a geometric spreading correction to a reference distance in the middle of each of the two intervals (90 and 200 km). For the close-in-distance interval, we expect a significant contribution from near-field terms and apply a 1/r2 correction. For the further-away-distance interval, we apply a 1/r correction. Because events have different numbers of waveforms for each magnitude and distance bin (fig. S7), we first calculate the mean PGDt waveform for each earthquake and then the mean of all average event waveforms. To add a measure of uncertainty, we bootstrap the analysis by removing 20% of the waveforms for each earthquake and repeat the analysis 100 times. The results for both distance ranges indicate that, very early on, the near-source displacements show different growth rates as a function of final magnitude. By 10 s, the M8.8 to M9.1 magnitude bin begins to display faster growth than smaller magnitude events. Earthquakes in the M8.1 to M8.3 and M7.6 to M7.8 magnitude bins are indistinguishable from one another in the first 10 s but, soon after, diverge and can be distinguished. These results echo and are consistent with the moment acceleration patterns seen in the STFs (Fig. 2) and, critically, offer a potential path forward for observing these differences in real time. We note that, by studying the PGDt growth following P-wave onset, we have eliminated an extra delay that will slow down an algorithm that leverages these observations. This delay is from the finite travel time of elastic waves from the earthquake to the HR-GNSS stations. For the events studied here (Fig. 2) at the closest sites, this travel-time delay can be anywhere between 5 and 25 s depending on the network and the location of the hypocenter.

Fig. 4 Average behavior of the PGDt waveforms for four magnitude and two hypocentral distance (Rhyp) bins.

For each earthquake, we average the PGDt waveforms to create a single-event waveform. We then averaged all the event waveforms within a given bin. Shaded regions represent the one σ uncertainties from bootstrap analysis. The gray area is the noise level of HR-GNSS within which the data cannot be reliably interpreted. Before averaging, for each waveform, we apply a geometric spreading correction to a reference distance in the middle of the interval (90 and 200 km). Dashed lines are centroid times for Mw 7.5 and Mw 8.0 earthquakes from a scaling law (9).


At face value, these results may seem to be in contrast with recent analyses of STFs from Meier et al. (6). In that study, after normalizing the STFs by centroid time and peak moment rate, it was concluded that large (M7+) earthquakes exhibit universal early growth. The critical difference between the two studies is the time scale of analysis—Meier et al. imply that all earthquakes start the same, while the work presented here shows that, after such “simple” nucleation, earthquakes of eventually different final magnitudes rapidly develop distinguishable characteristics.

Recent work on displacements from inertial recordings found no scaling in the earliest P-wave onsets (5). As we noted, these signals are high-pass–filtered by necessity, so it is possible to argue that there are longer period features that could foretell that a large event is likely. However, because of the noise levels typical of HR-GNSS waveforms (approximately 1 to 2 cm), it is very difficult to see these earliest-onset P waves; thus, whether such a signal exists remains unresolved. Nonetheless, because our analysis of STFs and HR-GNSS data is consistent, we favor the notion that, in the earliest stages of these earthquakes (1 to 3 s), events cannot be reliably distinguished from one another (5). However, at some subsequent time, distinct features emerge with different moment accelerations and rates of growth of near-field displacement.

A physical model is necessary to explain the relationship between final magnitude and moment acceleration evident in Fig. 2 and in the growth of the PGDt waveforms (Fig. 4). We propose that there is an initial chaotic growth phase that, at first, has the same characteristics for all events, irrespective of their eventual final magnitude. This initial phase likely goes on for longer time periods for larger earthquakes (1), and subsequently, the ruptures organize into slip pulses with measurable properties diagnostic of final magnitude. Here, we note that two slip pulse models are admissible: a self-similar pulse model where the pulse width grows linearly with time during rupture at a rate that is magnitude dependent and a steady-state pulse model where the pulse width remains unchanged as the rupture evolves. Recent work has shown that average pulse widths scale with magnitude (9) but could not distinguish between these two models. However, scaling of the moment accelerations (Fig. 2) perhaps favors the self-similar pulse model. Since pulse width and slip for these large events scale with the cubed root of final moment (M0) (9, 16), one expects that in a steady-state pulse, where the width remains unchanged, the slipping area grows in the same way (proportional to the power of 1/3). Thus, for a particular time window since rupture onset, moment, and thus moment rate and moment acceleration, is proportional to slip and rupture area. This then predicts scaling proportional to M02/3. However, if area grows at a faster rate because the pulse width is changing (as in the self-similar model), then scaling closer to M01 is possible. The median moment accelerations (Fig. 2, blue stars) are consistent with this transition from a chaotic growth phase to a self-similar pulse. The figure suggests that, early on (e.g., at 2 s; Fig. 2A), moment acceleration scales proportional to moment raised to the power of 1/2 or 1/3 At later times, moment acceleration scales faster, closer to linearly, with moment. The transition between these two regimes is magnitude dependent. For example, at 5 s, events with M < 7.5 have transitioned to linear growth, while larger events retain the slower scaling. By 15 s, the same is true for events with M < 8, and last, by 20 s, all events show linear scaling. This is also consistent with what has been observed recently in near-field P waves from strong motion records for moderate events (5). It is in this second stage when the pulse has developed (+10 s) that easily distinguishing features in the near-field HR-GNSS data become evident. Here also, two magnitude-dependent growth patterns are evident (Fig. 3): The time it takes for the moment acceleration or PGDt features to stabilize is longer for larger magnitude events (fig. S8). Note that, because of finite-source effects and depending on the source to station distance, a particular HR-GNSS site will be sensitive to the slip and rise time of some finite portion of the fault closest to it (17). As a result, the rise time observed at a particular PGDt record is an upper bound of the slip rise time (10). Nonetheless, both the teleseismic and HR-GNSS results imply that the time it takes for the pulse width to stabilize is also magnitude dependent; the larger an event is, the longer it takes for its rupture pulse to develop, grow, and stabilize. This time scale is much shorter than the time scale of the rupture itself, thus revealing deterministic features.

These findings are relevant for hazard warning systems, in particular for earthquake early warning (EEW), where every second gained matters. The modern EEW paradigm centers on forecasting the shaking intensity at a site of interest some distance away from an earthquake. Recent studies have shown that uncertainty of an event’s final magnitude leads to a tradeoff between the possible warning time and the certainty of the ground motion forecast (7). At present, a fast warning can only be made with large uncertainty in the associated shaking. While waiting for more information about the source significantly reduces the uncertainty, this comes at the cost of shrinking the warning time and enlarging the blind zone where no warning is possible. Modern EEW systems rely on a combination of point source and finite fault algorithms, neither of which is designed to exploit weak determinism. Point source algorithms focus on the first 1 to 3 s of a P wave as measured by inertial seismic sensors (1820); thus, it is unlikely that they will overcome the familiar problem of magnitude saturation for the largest events. Of course, for the smaller events in our dataset, by 10 to 15 s, the earthquakes are in the later portion of their source processes. As a result, from an early warning perspective, the results here would still leave room for future improvements in speeding up warning times. Finite fault algorithms (2123) can more effectively deal with large magnitude events but rely on the entire rupture (or a significant portion of it) being complete before they produce reliable magnitude estimates.

Our findings argue that a different paradigm, both in terms of the observational platform and the algorithm design, is possible. While we are not currently able to generate real-time STFs for an evolving rupture (Fig. 2), rapid imaging of an unfolding large earthquake can be made from broadband geodetic measurements as close as possible to a large earthquake (Fig. 3). For large continental ruptures, this analysis can be performed with currently existing HR-GNSS technology, but for subduction zone environments, where most of the rupture is offshore, new ocean bottom technologies will be necessary. Seafloor HR-GNSS exists for monitoring long-term deformation (24, 25), but real-time high-rate collection does not, at this time, seem feasible. Alternative measurements such as absolute pressure (26) or ocean-bottom strain (27) will need to be further developed. Alternatively, far-field measurements of near-instantaneous gravity perturbations (28) are another potential path forward for identifying the very largest earthquakes. Thus, while our results augur a positive future for rapid earthquake hazard characterization, they also illuminate the pressing need for new and more diverse measurement techniques.


STF analysis

The STF database of Vallée et al. (12) was generated via a deconvolution of synthetic Green’s functions from globally distributed body waves. For each event, an “average” STF was generated by combining apparent STFs from each station recording the earthquake.

The finite fault database of Hayes (1) used a nonlinear wavelet domain inversion approach (29). For each earthquake, teleseismic body wave (P and SH, band-pass–filtered between 1 and 200 s) and surface wave (Rayleigh and Love, band-pass–filtered between 200 and 500 s) data were inverted. Rupture velocity can be fixed or allowed to vary, a decision largely based on data quality and/or evidence from other studies of the event in question. Fault planes were divided into a series of subfaults along the strike and dip directions, and the inversion solved for the slip amplitude, slip direction, rise time, and rupture initiation time of each subfault, where subfault STFs were modeled with an asymmetric cosine function (2931).

Ye et al. (1) used a similar kinematic source inversion approach to generate finite fault models for large subduction zone earthquakes from 1990 to 2015. They inverted teleseismic P waves (band-pass–filtered between 1 and 200 s) onto a planar fault surface broken into a series of subfaults and parameterized each subfault slip rate function as a series of overlapping triangles. Thus, the inversion is linear, solving for the amplitude of each triangle and minimizing the misfit between model and observations in a least-squares sense. Rupture velocity was fixed (although a series of models were run for each event to determine the best-fitting rupture velocity).

While there can be some differences between the three datasets for any one event, by and large, the STFs share the same gross features (fig. S9). For each of these three datasets, we calculated the moment acceleration at 2, 5, 10, and 20 s (Fig. 1) simply as the slope of a line joining the event origin time to the moment rate at each time of interest. We present the moment acceleration values individually for each event. We also calculated the median moment acceleration in discrete magnitude bins (e.g., Fig. 2 and figs. S1 to S3) and compare the results between datasets to examine them for internal consistency.

HR-GNSS processing

We analyzed data from the HR-GNSS waveform database of Ruhl et al. (15). We selected large megathrust events for which there exists good regional seismic data as well. A total of 12 events were selected. The waveforms for each event were processed as follows. First, we obtained P-wave picks (arrival times) from regional networks. The strong motion data for these events are available for download from a database used to test EEW algorithms (32). We then fit a fourth-order surface to the geographically distributed P-wave picks and used that geometric surface to interpolate the expected P-wave arrival time at the HR-GNSS stations. We performed a bootstrap analysis where we randomly removed 20% of the seismic stations and performed the analysis and predicted P-wave arrivals at the removed sites. We repeated this procedure 100 times to estimate errors in the interpolation. The median pick time error is 0.5 to 1.5 s depending on the event. After defining the pick time at each HR-GNSS site, we cropped the waveform from this pick time to anywhere between 60 and 180 s after. Next, we use the three components of displacement (N, E, Z) to obtain the total displacement waveform (D) such that D = (N2 + E2 + Z2)0.5. From this, we obtained the time evolution of PGDt. PGDt is defined asPGDt(t)=max[D(τ)],tp<τ<twhere tp is the P-wave arrival time at a given station. The PGDt waveforms were color-coded by hypocentral distance and are presented in Fig. 2. We chose hypocentral distance instead of distance to the rupture because we are interested in the early evolution of rupture, which is more properly referenced to the point of rupture nucleation, the hypocenter. We analyzed the early PGDt growth patterns by stacking waveforms in four magnitude bins (M7.1 to M7.4, M7.6 to M7.8, M8.1 to M8.32, and M8.9 to M9.1) and two hypocentral distance bins (30 to 150 km and 150 to 250 km). Because events have different numbers of waveforms for each magnitude and distance bin, we calculated the mean PGDt waveform of each earthquake and then the mean of all average event waveforms. To add a measure of uncertainty, we bootstrapped the analysis by randomly removing 20% of the waveforms for each earthquake and repeated the analysis 100 times (these are shown as shaded regions in Figs. 2 and 3). An analysis of the errors in pick times from this process is also shown in figs. S5 and S6.


Supplementary material for this article is available at

Fig. S1. Seismic observations of weak determinism.

Fig. S2. Seismic observations of weak determinism.

Fig. S3. As in Fig. 1 and figs. S2 and S3.

Fig. S4. Same as Fig. 2 but color-coded by minimum distance to the rupture instead of by hypocentral distance.

Fig. S5. Bootstrap estimate of travel time interpolation errors for the 12 events with HR-GNSS data.

Fig. S6. P-wave travel times observed at seismic sites (blue), compared to a reference velocity model (PREM; dashed black line).

Fig. S7. Distribution of HR-GNSS stations used in this study.

Fig. S8. Time to maximum PGD for waveforms used in this study.

Fig. S9. Comparisons of STFs from the three STF databases used (1214), for the earthquakes used in HR-GPS analyses (Fig. 2).

Table S1. Events with HR-GNSS data analyzed in this work.

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 D. Shelly and two anonymous reviewers for comments that helped to improve this manuscript. We thank M. Vallée and L. Ye for providing STF datasets and D. Goldberg, T. Lay, M. Vallée, and W. Yeck for discussions. Funding: The authors acknowledge that they received no funding in support of this research. Author contributions: G.P.H. conducted the STF analyses. D.M. performed HR-GNSS analyses. D.M. and G.P.H. jointly interpreted results and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: STF data from Vallée et al. (2011) (12) can be downloaded from Finite fault models and STFs from Ye et al. (2016) (14) can be accessed from Finite fault models and STFs from Hayes (2017) (13) can be downloaded from the USGS Advanced National Seismic System Combined Catalog ( HR-GNSS data can be downloaded at 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.
View Abstract

Navigate This Article