## Gravity gets into the earthquake game

Earthquakes generate large movements of mass, which slightly change the gravitational field. Unlike the elastic waves that propagate from the earthquake, the gravity perturbations travel at the speed of light. Vallée *et al.* have finally observed these gravity perturbations in seismometer records from the great Tohoku earthquake in Japan in 2011. The signal would have allowed an accurate magnitude estimation in minutes, rather than hours, for this catastrophic earthquake.

*Science*, this issue p. 1164

## Abstract

After an earthquake, the earliest deformation signals are not expected to be carried by the fastest (*P*) elastic waves but by the speed-of-light changes of the gravitational field. However, these perturbations are weak and, so far, their detection has not been accurate enough to fully understand their origins and to use them for a highly valuable rapid estimate of the earthquake magnitude. We show that gravity perturbations are particularly well observed with broadband seismometers at distances between 1000 and 2000 kilometers from the source of the 2011, moment magnitude 9.1, Tohoku earthquake. We can accurately model them by a new formalism, taking into account both the gravity changes and the gravity-induced motion. These prompt elastogravity signals open the window for minute time-scale magnitude determination for great earthquakes.

Earthquakes involve the displacement of large amounts of mass, which modifies the gravity field. This effect is not restricted to a permanent gravity change due to the final mass redistribution [e.g., (*1*–*3*)] but is also induced by the transient density perturbations carried by seismic waves. During the wave propagation, an observer feels attracted by the compressed parts of the medium and repelled by its dilated parts, with a global net effect depending on the earthquake mechanism. The gravity perturbations are transmitted at the speed of light (3.10^{5} km/s), far faster than the first-arriving (*P*) elastic waves that travel at 6 to 10 km/s in the crust and upper mantle. Additionally, at distances close to a large earthquake, it is difficult to estimate the event magnitude from the information provided by elastic waves, even when the area is densely instrumented by seismometers. In the case of the 2011, moment magnitude (*M*_{w}) 9.1, Tohoku earthquake, the near-real-time magnitude provided by the authoritative Japan Meteorological Agency (*4*) was 7.9 and was corrected only 3 hours later to 8.8 (*5*). This underestimation is due to the fact that real-time local magnitudes are generally derived from instrumental peak amplitudes, which are poorly correlated with moment magnitude when the earthquake is large. Detection of the gravity perturbations would provide a much faster method for estimating the size of fault ruptures.

The theoretical relations between the elastic and gravitational fields are well known [e.g., (*6*)], and analytical computations predicted the expected gravity change Δg* ^{P}* before the arrival of the

*P*waves in full-space (

*7*) and half-space (

*8*) models. The amplitude of Δg

*increases with increasing elastic deformation of the medium, and this growth is faster when the earthquake seismic moment rises quickly. Therefore the large-magnitude and short-duration earthquakes offer the best observation potential. As shown by (*

^{P}*7*) and (

*9*), the optimal distances for detecting Δg

*are not the closest ones to the earthquake. As long as the earthquake is in its accelerating phase, with seismic moment growing faster than quadratically with time, the gravity acceleration expected just before the arrival of the hypocentral*

^{P}*P*wave grows with the distance from the earthquake. For a very large earthquake of magnitude 9, for which the accelerating phase lasts on the order of 100 s, the gravity signal is expected to increase as a function of distance at least up to 800 km from the earthquake. This effect stems from the fact that an earthquake source is not instantaneous and from the increasing duration of the pre-

*P*observation window with distance. It is, however, not the only reason that close distances can be unfavorable to observe the early gravity signals with seismometers or gravimeters coupled to the ground. A previously overlooked phenomenon arises from ground accelerations themselves induced by the gravity changes, which tend to reduce the observability of the signal early in the earthquake (“early times”).

The Tohoku earthquake in Japan (11 March 2011, *M*_{w} = 9.1) was a suitable event to search for such prompt gravity-induced signals. The earthquake shares a similar magnitude with the 26 December 2004 Sumatra earthquake but benefits from a shorter source duration and a better seismic station coverage. We retrieved all the regional broadband seismic data available at the Incorporated Research Institutions for Seismology (IRIS) data center (*10*) at distances up to 3000 km from the earthquake, as well as the broadband data from the F-net Japan network (*11*). Vertical waveforms are cut at the *P*-wave arrival time, deconvolved from the instrument acceleration response, and bandpass filtered between 0.002 and 0.03 Hz in order to get rid of most of the oceanic noise. We use two-pole high-pass and six-pole low-pass causal Butterworth filters; the simple data processing procedure is provided in data file S1 (*12*). We hereafter denote the observed signals as . We further select waveforms based on a signal-to-noise criterion, requiring that remains in the ±0.8 nm/s^{2} range in the 30-min-long window preceding the earthquake. Most of the nine regional stations thus selected (Fig. 1) are stations from the IRIS and GEOSCOPE global networks (*10*, *13*) known for their high quality (*14*). This data set is complemented by two stations from the F-net network [Fukue, Japan (FUK) and Shari, Japan (SHR)], selected to improve azimuthal and distance coverage without adding redundancy. The range of distances considered here, from 400 to 3000 km, extends the analysis made by Montagner *et al*. (*9*) based on gravimetric and seismic data located about 500 km from the epicenter. Their results were promising because they show that a signal is very likely to be present (from a statistical point of view), even at these short distances. We show at the selected stations, including the 30-min-long presignal window used to evaluate data quality (Fig. 1B). In the time frame between the earthquake’s origin time and the *P*-wave arrival, most stations show a consistent downward acceleration trend, particularly pronounced at stations located 1000 to 2000 km west of the earthquake [Mudanjiang, China (MDJ); FUK; Incheon, Korea (INCN); Zhalaiteqi Badaerhuzhen, China (NE93); and Baijiatuan, China (BJT)], where it reaches 1.6 nm/s^{2}. Even if this recorded acceleration is smaller by a factor of more than 10^{5} than the following elastic *P*-wave train (fig. S1), it remains above the seismic noise owing to the large earthquake size.

The modeling of such signals first requires clarifying the relation between and the physical fields. After correction from its response in acceleration, a seismometer is essentially sensitive to the difference between the gravity perturbation and the ground acceleration [e.g., (*6*)]. Combining the upward seismometry convention for with the downward convention for the pre-*P* gravity perturbation and ground acceleration , this leads to . We neglect additional contributions from the free-air and Bouguer gravity changes on the basis of their smaller amplitudes compared with the two other terms (*12*). originates from the space-time evolution of the displacement generated by the earthquake (*6*–*8*) and can therefore be modeled in a realistic Earth model by an elastic wave propagation simulation. We use here the AXITRA code (*15*), based on a discrete wave-number summation (*16*), as further detailed in the supplementary materials (*12*). In Fig. 2, is shown for two stations at different distances, Inuyama (INU) in Japan and MDJ in northeast China. In both cases, the perturbation is initially positive, because at early times the contributions to come from the volume elements located below the earthquake, which are compressed by the *P* waves. At MDJ station, the sign of the perturbation changes due to the increasing effect of the volume elements located closer to the station, which are dilated by the *P* waves (fig. S2). The same effect explains the minor inflexion observed at INU station when approaching the *P*-wave arrival.

The discrepancy between and , in particular at INU station, indicates that the ground acceleration cannot be neglected. Such pre-*P* ground acceleration exists because the gravity perturbation Δg* ^{P}*, occurring simultaneously with the earthquake rupture, itself acts as a secondary source of elastic deformation in the whole Earth (

*17*). We first calculated Δg

*not only at the station but at all locations where this secondary source can create waves arriving at the station before the hypocentral*

^{P}*P*-wave arrival (in an homogeneous medium, this would be a ball centered on the station with a radius equal to the distance between station and hypocenter). We then applied to each of these secondary source locations a body-force equal to ρΔg

*(where ρ is the density) and computed their radiated elastic waves with a seismic wave simulation method. Their overall wavefield provides (*

^{P}*12*) (figs. S3 and S4). This new approach accounts for both gravity changes and gravity-induced motion and offers a concrete method able to reproduce . It also explains why can be referred to as the elastogravity signal preceding the

*P*-wave arrival.

We found that tends to compensate at early times, as shown in Fig. 2 for both INU and MDJ stations. This effect, predicted by Heaton (*17*), can also be understood from the insightful infinite medium configuration in which there is a theoretical full cancellation of Δg* ^{P}* by (

*12*), which our modeling approach reproduces very well (fig. S4). Earth’s surface, and to a lesser extent the heterogeneities inside Earth, break the symmetries leading to this full-space cancellation (

*12*), and realistic simulations show that and diverge before the hypocentral

*P*arrival time (Fig. 2). As a consequence, the difference between these two quantities offers a much larger observation potential than predicted by (

*17*) using an oversimplification of Earth’s response. For the INU station, we show that is stronger than , which leads to the negative signal recorded at this station. A ground-attached gravimeter would record the same, implying that it is paradoxically less sensitive to the gravitational field perturbations than to its induced effects. Recording exclusively the positive (and stronger) gravity perturbation would require an instrument uncoupled to the ground. The global feedback effect is less drastic when the station is farther away from the earthquake, as illustrated in Fig. 2 with the MDJ station. In this case, for about 60 s we observed an almost full cancellation of the signal, but the last 100 s before the

*P*-wave arrival were increasingly dominated by . This explains why the observed signal is very pronounced at MDJ station and at the neighboring INCN, FUK, NE93, and BJT stations.

We gathered the observed elastogravity signals and present our modeling results for the 11 regional stations (Fig. 3). The data and synthetics are systematically in very good agreement. The oscillations present at the shortest considered period (33 s), also in pre-earthquake time windows, are related to Earth’s natural noise. The modeling helps us understand several important features of the waveforms. The observability of the signals reaches a maximum at about 1000 to 1500 km from the earthquake, as illustrated by the high-quality MDJ station. These distances are favorable because the hypocentral *P* wave arrives 130 to 190 s after the origin time, when the Tohoku earthquake already had released most of its seismic moment (*12*). Closer to the earthquake, the pre-*P* time window is short relative to the earthquake duration, and the canceling effects between and (Fig. 2) reduce even more the window where the signal can be observed. At further distances [e.g., Ulaanbaatar, Mongolia (ULN) and Xi’an, China (XAN) stations], where the latter two effects marginally affect the signal, the weaker amplitude is simply due to the distance dependence of (*7*). Finally, independently of the distance, we observed and modeled the strong azimuthal effect due to the earthquake’s focal mechanism. Because the Tohoku earthquake is a thrust event occurring on the north-south-striking, shallow-dipping subduction interface, is very small for the stations to the north [Magadan, Russia (MA2) and SHR].

These observations of the elastogravity signals preceding the *P* waves, and their successful quantitative modeling, strongly motivate their use for an early magnitude estimate of megathrust earthquakes. Based on the proportionality in infinite space between and the second temporal integral of the moment time function *m* of the earthquake (*7*), we expect a strong dependency of on the earthquake magnitude. We show in Fig. 3 the synthetic signals for a realistic Tohoku-type earthquake, down-scaled to *M*_{w} = 8.5. Keeping the assumption of a triangular moment rate function (*12*), scaling relationships (*18*, *19*) predict such an earthquake to have growing half as fast as and with half the duration of the Tohoku earthquake. As expected from the respective moment time evolutions, the signal amplitudes of the simulated *M*_{w} = 8.5 earthquake are about half the ones of the Tohoku earthquake at early times [INU and Matsushiro, Japan (MAJO) stations] and become increasingly smaller at late times, approaching the moment ratio value of 1/8. Even in this simulation, where the *M*_{w} = 8.5 earthquake lasts 70 s (which is short for such a magnitude) (*19*, *20*), never exceeds 0.5 nm/s^{2}. This simple test therefore shows that detection of pre-*P* acceleration amplitudes reaching 1 nm/s^{2} is a direct evidence of an earthquake with a seismic moment at least twice as large (*M*_{w} > 8.7), hundreds to thousands of kilometers away.

This estimate can be refined when a large earthquake has been detected and its epicenter has been located with local data (which can be done in the tens of seconds after the origin time). In this case, based on the theoretical or empirical (with classical triggering techniques) *P*-wave arrival time at regional stations, it is straightforward to extract the pre–*P*-wave arrival time window. Compared with the usual post–*P*-wave time window recording the complex regional elastic wavefield, the former window provides both an earlier and a simpler way to evaluate how large the earthquake was. In this respect, Fig. 3 can be directly used to get a reliable lower bound of a megathrust earthquake magnitude. As the Tohoku earthquake has a short duration compared with its magnitude (*12*, *19*, *20*), it is unlikely that a smaller magnitude earthquake generates larger values at a given distance. Observing values of about 1300 km from the earthquake (as in the case of MDJ, FUK, or INCN stations), just before the *P*-wave arrival time, is therefore evidence of the occurrence of a *M*_{w} > 9 earthquake. If such an approach were followed for the Tohoku earthquake, using these stations where the *P* arrival times are less than 180 s, a lower bound of its huge magnitude would have been reliably detected 3 min after the origin time. Using additional elastogravity signals recorded at further distances (like ULN or XAN stations) delays the time at which a first magnitude can be provided but offers the potential to provide an exact magnitude determination. Such data can indeed better detect that the earthquake has stopped growing (*7*), a necessary condition to move from a lower bound estimation to an exact determination.

The possibility of detecting, 3 min after the origin time, that the Tohoku earthquake had a magnitude larger than 9 has to be compared with our current ability to quantify large earthquakes’ magnitudes. The determination of the moment magnitude in the minutes after an earthquake is possible with local data [e.g., (*21*, *22*)], but for large-magnitude events, this is complicated by finite-source effects. Currently, moment magnitudes are more efficiently determined at distances far from the source (*23*–*25*), with a fundamental limitation imposed by the time needed for elastic wave propagation. Even the fastest available methods (*24*) are unlikely to provide a reliable magnitude estimate within the first 20 min after the earthquake.

Synthetic signals of the *M*_{w} = 8.5 earthquake show that maximum amplitudes are lower than 0.5 nm/s^{2} everywhere, making individual detection difficult, even with excellent broadband seismometers located in quiet sites. We therefore emphasize the strong benefit of installing and maintaining high-quality sensors at regional distances from potential large earthquakes, such that stacking or coherence techniques can be applied to detect early gravity signals from earthquakes in the 8 to 9 magnitude range. At lower magnitudes, the potential detection of such signals depends on our ability to separate the gravity signal from the background seismic noise. This can be done in principle by measuring the gradient of the gravity perturbation between two or more seismically isolated test masses. Relevant technologies are being developed in the context of low-frequency gravitational-wave detectors, with concepts such as torsion bars antennas (*26*, *27*), superconducting gravity gradiometers (*28*, *29*), and atom interferometers (*30*, *31*). In the first two concepts, the test masses are linked to the ground by a common frame; the displacements driven by the seismic noise and affecting the gravity measurement can be made very similar for the two masses, and they are hence rejected by the differential measurement. In an atom interferometer, the phase of a laser beam is sensed by its interaction with two or more atomic clouds, giving an intrinsic partial immunity to the background seismic noise. The gravity gradient is, however, much weaker than the gravity itself, and making its measurement feasible should motivate further research to overcome additional challenges besides the suppression of seismic noise.

## Supplementary Materials

This is an article distributed under the terms of the Science Journals Default License.

## References and Notes

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵Information on materials and methods is available in the supplementary materials.
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
**Acknowledgments:**We are grateful to the developers of AXITRA, and more specifically to O. Coutant for making his code publicly available, both in the moment tensor and in the single-force configurations. We thank engineers involved in the installation, maintenance, and data distribution of broadband seismic stations. High-quality signals come from the Global Seismograph Network (https://doi.org/10.7914/SN/IU), New China Digital Seismograph Network (https://doi.org/10.7914/SN/IC), National Research Institute for Earth Science and Disaster Resilience /F-net, GEOSCOPE (https://doi.org/10.18715/GEOSCOPE.G), and Northeast China Extended Seismic Array (https://doi.org/10.7914/SN/YP_2009) networks, all publicly available at the IRIS data management center (http://ds.iris.edu/ds/nodes/dmc/), Institut de Physique du Globe de Paris (IPGP) data center (http://datacenter.ipgp.fr/), or the F-net data center (www.fnet.bosai.go.jp/top.php). Global Centroid-Moment-Tensor Project (www.globalcmt.org/) source parameters of the Tohoku earthquake were very valuable for this study. The careful reviews of three anonymous reviewers were very valuable to clarify observational and theoretical aspects of the initial manuscript. We acknowledge T. Heaton for stimulating discussions that pushed us to carefully study the role of the gravity-induced acceleration. J. Harms provided valuable comments on the relationships between gravity and elastic perturbations. Most numerical computations were performed on the S-CAPAD platform, IPGP, France. We acknowledge the financial support from the GEOSCOPE program (funded by CNRS Institute for Earth Sciences and Astronomy and IPGP), the UnivEarthS Labex program at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02), and the Agence Nationale de la Recherche through grant ANR-14-CE03-0014-01. Seismic Analysis Code (http://ds.iris.edu/ds/nodes/dmc/software/downloads/sac/) and Generic Mapping Tools (http://gmt.soest.hawaii.edu/projects/gmt) free software packages were used for data processing and figure preparation.