Multiple Observations of the Prompt Elastogravity Signals Heralding Direct Seismic Waves
Abstract
The recent first observations of the prompt elastogravity signals (PEGS) induced by the 2011 Mw = 9.1 Tohoku megathrust earthquake generated interest in how these signals might best be observed, especially for lower-magnitude events. Simulations of these signals preceding the direct P wave, for different depths and focal mechanisms, first reveal that shallow strike-slip earthquakes offer a better detection potential than subduction megathrust earthquakes. Consistently, clear PEGS are observed at several broadband seismometers during the 2012 Mw = 8.6 Wharton Basin earthquake. Due to their short source durations, large deep earthquakes are then shown to have an even larger detection potential, confirmed by the successful seismological observations for the 2018 Mw = 8.2 Fiji and 1994 Mw = 8.2 Bolivia earthquakes. Detection is even improved when an earthquake is recorded by a number of good-quality stations, allowing for stacking techniques. Thanks to the deployment of the USArray network across Alaska, the recent 2018 Mw = 7.9 off-Alaska earthquake (strike slip) is thus observed with an excellent signal-to-noise ratio. Array stacking is also shown to reveal the PEGS induced by the large 2010 Mw = 8.8 Maule megathrust earthquake, for which individual observations are impeded by the long-lasting radiation generated by a distant large earthquake. As a whole, we show new observations and successful modeling of the PEGS for five earthquakes in the 7.9–8.8 magnitude range. These findings demonstrate that, even without considering promising future instruments, the PEGS detection is not restricted to exceptional events, confirming their potential for magnitude and focal mechanism determination within the few minutes following a large earthquake.
Key Points
- We show new seismic observations of prompt elastogravity signals (PEGS) for five Mw7.9-8.8 earthquakes
- Magnitude, focal mechanism, depth, and source time function control the signal shape and amplitude
- Single-station or array-derived PEGS are a new promising data type for large earthquake monitoring
1 Introduction
When an earthquake occurs, the Earth's density is perturbed in a volume whose size increases with time as the elastic waves propagate. The Newtonian Poisson's equation then predicts that the gravitational field is immediately perturbed everywhere and not only inside this time-evolving volume perturbed by the elastic waves. Even if we know from Einstein theory that the gravitational perturbation is in fact not felt immediately (but at the speed of light), this distinction can easily be neglected at the scale of the Earth. An apparatus extremely sensitive to the gravity changes would therefore be able to detect any distant earthquake as soon as it started (Harms et al., 2015).
These prompt elastogravity signals (PEGS) are readily present in the self-gravitating equations governing the earthquake-induced motion (e.g., Dahlen & Tromp, 1998, p. 84), but their observations have only been provided recently by analyzing local and regional records of the 2011 Mw = 9.1 Tohoku earthquake (Montagner et al., 2016; Vallée et al., 2017). Vallée et al. (2017) show that acceleration signals in the 0.002- to 0.03-Hz frequency range, recorded by excellent broadband seismometers located 1,000 to 2,000 km away from the earthquake, offer the most striking evidence. At these stations, where the P wave arrives 2 to 3 min after the earthquake origin time, an early signal is clearly measurable some tens of seconds before this first seismic wave arrival. As the PEGS can be accurately modeled and are shown to be strongly magnitude dependent (Juhel, Montagner, et al., 2018; Vallée et al., 2017), their detections have the potential of providing a highly valuable rapid estimate of the earthquake size.
This motivates further characterization of these signals, both by modeling and additional observations. In a first section, we numerically explore the influence of the earthquake focal mechanism and depth on the expected signal amplitude, using realistic scenarios of earthquakes with Mw = 8.5. This reveals that at a given magnitude, shallow strike-slip and deep earthquakes have a better potential to be recorded than megathrust earthquakes (thus compensating for the larger magnitudes of megathrust earthquakes). Additionally, these synthetic cases indicate in each configuration where to expect the largest signal.
Helped by these simulations, we analyze in the next sections the broadband records of large earthquakes of the last 25 years and first find three additional events with direct observational evidence: the 11 April 2012 Mw = 8.6 Wharton Basin earthquake (shallow strike slip), the 19 August 2018 Mw = 8.2 Fiji earthquake (depth ≃ 560 km), and the 9 June 1994 Mw = 8.2 Bolivia earthquake (depth ≃ 640 km).
As the tiny PEGS are close to the ambient and/or instrumental noise level, we expect the signal-to-noise ratio (SNR) to be enhanced by stacking techniques, allowing for detection of smaller or nonoptimally recorded earthquakes. The recent 23 January 2018 Mw = 7.9 − 8 off-Alaska earthquake benefits from the excellent coverage of the USArray deployment in Alaska, and we will show in this case how several stacking approaches result in an excellent extraction of the PEGS. Another successful application is shown for the 27 February 2010 Mw = 8.8 Maule megathrust earthquake, where individual observations suffer from the long-lasting low-frequency ground motion generated by a Mw = 7 earthquake in Japan, 10 hr before.
In order to illustrate all the conditions to be fulfilled for a precise observation of these signals, we will discuss in the last section some interesting undetected large events. The first case is the 24 May 2013 Mw = 8.3 Okhotsk deep earthquake. The reason here is a noise level even higher than for the Maule earthquake, with a Mw = 7.4 earthquake occurring 12 hr before in the Fiji region. The second case is the 2004 Mw = 9.1–9.2 Sumatra megathrust earthquake. Even if its magnitude is similar to the Tohoku earthquake, its source duration is much longer, and even more importantly, the moment function grows much slower in the first 60–80 s. PEGS with detectable amplitude levels therefore arrive too late with respect to the direct P wave.
In all these examples, we will show that waveform modeling always supports the observations, confirming that the approach taking into account both the direct gravity effect and the elastic gravity-induced effect (Vallée et al., 2017; see also Juhel, Montagner, et al., 2018) is an accurate formalism. The global scope of this study is therefore twofold: showing new observations of PEGS not restricted to an exceptional event and demonstrating the potential of retrieving earthquake source parameters from these early signals.
2 Expected Influence of the Focal Mechanism and Depth of the Earthquake
Three types of fault rupture historically led to large-magnitude events (Mw > 8), hence to a large redistribution of internal masses: shallow megathrust earthquakes, shallow strike-slip earthquakes, and deep earthquakes. In order to explore the influence of the depth and earthquake focal mechanism (hereafter parametrized by its strike ϕ, dip Δ, and rake λ) on the elastogravity perturbation, we select four fault geometries that illustrate these large-magnitude events. We consider two 20-km-deep reverse dip-slip events, accounting for realistic values of the dip angle in megathrust environments ((ϕ = 180°, Δ = 10°, λ = 90°) and (ϕ = 180°, Δ = 20°, λ = 90°)), a 20-km-deep vertical strike-slip event (ϕ = 0°, Δ = 90°, λ = 0°), and a 650-km-deep horizontal dip-slip event (ϕ = 0°, Δ = 0°, λ = 90°). These events are all scaled to the same magnitude (Mw = 8.5) and excited by the same source time function (STF; isosceles triangle with half-duration of 45 s), in the PREM Earth model (Dziewonski & Anderson, 1981).
Vallée et al. (2017) showed that the PEGS at frequencies higher than ≃ 0.002 Hz can be computed as the difference of two terms, a direct gravity perturbation Δg and a gravity-induced elastic acceleration 
. Such formalism enabled us to accurately model the perturbation induced by the Tohoku earthquake at single station locations, using the discrete wavenumber AXITRA code (Cotton & Coutant, 1997). The normal-mode formalism enables a faster computation of Δg in a straightforward mode summation, such that Δg surface maps of thousands of grid points can be computed at low computational costs (Juhel, Montagner, et al., 2018). Δg can thus be efficiently computed for the 12,000 surface points of our target grid (one point every 40 km for all locations closer to 20° from the epicenter location). The normal-mode computation of the induced acceleration 
however still requires a two-step approach, hence a larger computation time. Indeed, for each surface grid point, it first requires to compute the body force term ρΔg in the surrounding volume (discretized by hundreds of thousands of points for epicentral distances close to 20°) and then to sum all the associated secondary elastic waves (Figure S1). In order to make the numerical approach more tractable, the elastic Green's functions relating the secondary sources to any point at the Earth surface are precalculated, as further detailed in Supporting Information S1.
We show the vertical PEGS amplitudes at P wave arrival time for all four mechanisms in Figure 1. These amplitudes are extracted after bandpassing the acceleration waveforms between 0.002 and 0.03 Hz, consistently with the signal processing used in Vallée et al. (2017) and in the following sections. The perturbations are decomposed into their direct and induced components in Figures S2 to S5 of the supporting information. For a given magnitude and source duration, the shallow strike-slip and deep dip-slip earthquakes induce larger perturbations than the subduction megathrust earthquakes. The best observation potential is reached with the strike-slip event, where the PEGS amplitude exceeds 1.25 nm/s2 along the maximum tensional (T) and compressional (P) axes at regional distances (∼1,000 km away from the epicenter). Such distances offer optimal observations since the rupture is well developed at P wave arrival time and because of the distance dependence of the gravity perturbation (Harms et al., 2015). Amplitudes above 1.1 nm/s2 are observed around the P and T axes of the deep earthquake, that is, about 500 km away from the epicenter in the slip direction. We further expect that such values are lower bounds for real deep earthquakes, as they usually have shorter source durations than shallow earthquakes (Frohlich, 2006; Houston & Williams, 1991; Vallée, 2013), which increase the PEGS. In other words, amplitudes of 1 to 2 nm/s2 will be likely reached for deep earthquakes in the 8–8.5 magnitude range.
The PEGS induced by the subduction megathrust earthquakes do not exceed 0.6 and 1.0 nm/s2 for the 10° and 20° dipping scenarios, respectively. These lower amplitudes are first caused by the fact that such earthquakes generate elastic waves with amplitudes proportional to 
at low frequency (Kanamori & Given, 1981; Tsai et al., 2011). Additionally, at regional distances, such earthquakes cannot be observed at the Earth's surface close to their P or T axes.
These synthetic cases provide guidance for the search of elastogravity signals generated by earthquakes of the last 30 years: Even if the maximum observed magnitude is 8.6 for strike-slip earthquakes (11 April 2012 Wharton Basin earthquake) and 8.3 for deep earthquakes (24 May 2013 Okhotsk earthquake), these types of events offer an observation potential similar to the magnitude ≃ 9 megathrust earthquakes.
3 New Individual Observations and Models of the PEGS
3.1 The 11 April 2012 Mw = 8.6 Wharton Basin Strike-Slip Earthquake
3.1.1 Data
The 11 April 2012 Mw = 8.6 Wharton Basin earthquake is the largest strike-slip event ever recorded by broadband seismometers. We retrieved all the broadband data publicly available at distances closer than 2,000 km from the earthquake epicenter and processed them as in Vallée et al. (2017): The 1-hr-long pre-P vertical data are converted into acceleration using the instrumental response and bandpass filtered to remove both unreliable very low frequencies and microseismic noise. We use a two-pole causal Butterworth high-pass filter at 0.002 Hz (except at station IPM where low frequencies require to be filtered with a four-pole filter) and a six-pole causal Butterworth low-pass filter at 0.03 Hz. The five stations shown in Figure 2 (belonging to GE, MY, PS, and IU networks) are the ones where the obtained signal remains below ±0.8 nm/s2 in the 600-s-long time window preceding the earthquake origin time.
In Figure 2a, we show the PEGS expected amplitudes at P wave arrival time, using the methodology of section 2 and the Global Centroid Moment Tensor (GCMT) parameters of the earthquake (Ekström et al. (2012); Table 1). The GCMT STF, parametrized by an isosceles triangle with half duration of 45 s, is shown in Figure 3a. The two stations PSI and IPM are expected to have a clear positive signal (larger than 1 nm/s2 at IPM) due to their distance to the earthquake and their locations close to the T axis. On the other hand, the PEGS should be weaker and negative at CISI and MALK stations, located at further distances and closer to the nodal planes, and almost 0 at the nodal CHTO station. Data shown in Figure 2b (red lines) well agree with these expectations, both in polarity and relative amplitude, with clear PEGS at PSI and IPM stations.
| Earthquake | Source type | Origin time | Moment (N.m) | Latitude | Longitude | Depth (km) | Strike | Dip | Rake |
|---|---|---|---|---|---|---|---|---|---|
| 11 April 2012 Mw = 8.6 | GCMT | N/A | 9.14 × 1021 | 2.35 | 92.82 | 45.6 | 289 | 85 | 166 |
| Wharton Basin | SCARDEC | 08:38:36 | 1.36 × 1022 | 2.33 | 93.06 | 24.0 | 289 | 83 | 175 |
| 19 August 2018 Mw = 8.2 | GCMT | N/A | 2.52 × 1021 | −17.86 | −177.85 | 555.0 | 210 | 21 | −74 |
| Fiji | SCARDEC | 00:19:38 | 2.72 × 1021 | −18.16 | −178.11 | 572.0 | 223 | 21 | −66 |
| 9 June 1994 Mw = 8.2 | GCMT | N/A | 2.63 × 1021 | −13.82 | −67.25 | 647.1 | 302 | 10 | −60 |
| Bolivia | SCARDEC | 00:33:16 | 2.85 × 1021 | −13.84 | −67.55 | 631.0 | 310 | 10 | -53 |
| 23 January 2018 Mw = 7.9 | GCMT | N/A | 9.60 × 1020 | 56.22 | −149.12 | 33.6 | 257 | 64 | 4 |
| Gulf of Alaska | SCARDEC | 09:31:40 | 1.07 × 1021 | 56.00 | −149.17 | 21.0 | 262 | 72 | 23 |
| 27 February 2010 Mw = 8.8 | GCMT | N/A | 1.86 × 1022 | −35.98 | −73.15 | 23.2 | 19 | 18 | 116 |
| Maule | SCARDEC | 06:34:11 | 1.91 × 1022 | −36.12 | −72.90 | 24.0 | 24 | 18 | 117 |
- Note. For each earthquake, we provide the GCMT-type and SCARDEC-type source parameters. The GCMT-type configuration uses the full GCMT solution (Ekström et al., 2012), with the goal to predict the PEGS amplitude distribution at the P wave arrival time (Figures 2a, 4a, 5a, 6a, and 8a). Note that this simulation does not require any hypocentral time and that the GCMT double-couple solution is only shown here for information. The SCARDEC-type configuration uses the SCARDEC focal mechanism (Vallée & Douet, 2016) associated with the STFs shown in Figure 3. Epicentral coordinates and origin time are not determined in the SCARDEC method, and this information here comes from USGS. For shallow earthquakes, depths are from SCARDEC (as the use of surface reflected body wave makes this method very sensitive to depth), and for the two deep Bolivia and Fiji events, this information is taken from USGS. The SCARDEC-type configuration is used to model all the PEGS waveforms shown in this study (Figures 2b, 4b, 5b, 6b, 7, 8b, and 9). PEGS = prompt elastogravity signals; STFs = source time functions; GCMT = Global Centroid Moment Tensor; USGS = U.S. Geological Survey.
3.1.2 Modeling Results
Expected amplitudes shown in the map of Figure 2a are insightful for determining the most promising stations. The GCMT source description is however too simplistic to model the elastogravity waveforms when earthquakes have a temporal development not well represented by the prescribed isosceles triangle. This is the case for the 2012 Wharton Basin earthquake whose STF from the SCARDEC method (Vallée & Douet, 2016) accelerates faster than the GCMT triangular function (Figure 3a).
In order to simulate the PEGS waveforms, the 2012 Wharton Basin earthquake is thus modeled by a point source (located at the hypocentral coordinates) whose focal mechanism and STF are determined by the SCARDEC method (Table 1 and Figure 3a). We then follow the procedure described in Vallée et al. (2017): The direct elastic displacement wavefield u is computed in the Mantle PREM model (Dziewonski & Anderson, 1981) using the AXITRA method (Cotton & Coutant, 1997) and Earth flattening corrections (Muller, 1977). This provides access to the direct gravity perturbation Δg through integration over the time-evolving volume affected by u (Dahlen & Tromp, 1998). As shown by Vallée et al. (2017), the PEGS are the difference between Δg and the gravity-induced elastic ground acceleration 
, the latter term being calculated with the AXITRA method by considering that Δg also acts as a body force in the medium surrounding the observation point. The resulting synthetic vertical PEGS are finally processed with the same filtering as the observed waveforms and shown with the black line in Figure 2b.
The positive perturbations are well reproduced at the two large SNR stations (IPM and PSI), with a higher accuracy at IPM station. We expect the point source modeling to be less precise at PSI station, due to its closer distance to the earthquake (650 km); as the 2012 Wharton Basin earthquake propagated for more than 200 km west of the epicenter (Duputel et al., 2012; Meng et al., 2012; Satriano et al., 2012), predicted PEGS have an overestimated amplitude for close stations located east of the earthquake. CISI and MALK stations have a noisier signal than the other three stations, but the shape and amplitude of the negative trend preceding the P wave arrival are also well reproduced. Finally, CHTO station is the best station in terms of intrinsic quality, but, as also predicted in Figure 2a, no signal can be detected there due to the nodal location of the station with respect to the focal mechanism.
3.2 The 19 August 2018 Mw = 8.2 Fiji Deep Earthquake
3.2.1 Data
The 19 August 2018 Mw = 8.2 Fiji earthquake is the most recent of the three very large deep earthquakes (Mw≥ 8 and depth ≥ 500 km) recorded by modern broadband seismometers (the other two being the 9 June 1994 Bolivia and 24 May 2013 Okhotsk earthquakes). We retrieved all the broadband data publicly available at distances closer than 2,000 km from earthquake epicenter and processed them as in section 3. The four stations shown in Figure 4 (belonging to IU, II, and G networks) are the ones where the vertical acceleration signal remains below ±0.9 nm/s2 in the 600-s-long time window preceding the earthquake origin time.
Like the other Mw ≥ 8 deep earthquakes, the Fiji event has a dip-slip mechanism on a shallow-dipping plane, and we expect relatively strong PEGS at the locations shown in the deep earthquake scenario (see section 2). Using the same approach with the earthquake GCMT parameters (Table 1 and triangular STF shown in Figure 3b) provides the GCMT-expected amplitudes at the P wave arrival time, shown in the map of Figure 4a. MSVF station appears to be optimally located to record a negative signal. At the other three stations, we expect weaker PEGS, and likely difficult to detect at RAR station due to its distance to the earthquake. In terms of sign, we expect a positive polarity at NOUC and a negative polarity at AFI. Data shown in Figure 4b (red lines) well agree with these expectations, both in polarity and relative amplitude. Absolute amplitudes are however larger than predicted (reaching in particular almost 2 nm/s2 at MSVF station), indicating that moment rate grows faster than the GCMT STF.
3.2.2 Modeling Results
The GCMT-expected amplitudes shown in Figure 4a use a generic isosceles triangle whose duration is only magnitude dependent. As most deep earthquakes, the Fiji earthquake has however a shorter and more impulsive source process, as indicated by the SCARDEC STF shown in Figure 3b. To reproduce the PEGS at the four station locations, we thus follow the methodology detailed in section 4, with focal mechanism and STF determined by the SCARDEC method (Table 1 and Figure 3b). The resulting synthetic vertical PEGS (black lines in Figure 4b) first show that the obvious signal at MSVF station is well reproduced in shape and amplitude. At AFI and NOUC stations, the SNR is less favorable, but both the rapid acceleration decrease at AFI and the long-lasting acceleration increase at NOUC are well modeled. The modeled signal at RAR is well below the noise recorded at the station and is therefore consistent with the absence of detection at this station.
3.3 The 9 June 1994 Mw = 8.2 Deep Bolivia Earthquake
3.3.1 Data
The 9 June 1994 Mw = 8.2 Bolivia earthquake was recorded by only two permanent broadband seismometers whose distance to the earthquake and intrinsic quality offer the possibility to track the PEGS: the LPAZ station (GT network) located close to La Paz and the NNA station (II network) in Central Peru. We add to this data set one station in southern Bolivia (YUNZ) from the temporary BANJO experiment (Broadband ANdean JOint network, XE code). When high-pass filtered at 0.002 Hz with a four-pole causal Butterworth filter (to remove unreliable low frequencies), this station has a pre-event noise similar to the LPAZ station. The three considered stations are indicated in the map of Figure 5a, and their processed waveforms (see section 3) are shown in Figure 5b. In Figure 5a, we also add the expected PEGS amplitude at the P arrival time, based on the GCMT parameters of Table 1. This simulation uses the GCMT boxcar STF shown in Figure 3c and predicts that LPAZ station is located in an almost optimal location for PEGS detection. On the contrary, signal is expected to be small at NNA and YUNZ stations, making any detection unlikely there. The observed waveforms confirm these expectations, with a clear positive acceleration signal at LPAZ and no visible PEGS at YUNZ and NNA, even if NNA has a low pre-event noise.
3.3.2 Modeling Results
The source process of the 1994 Bolivia earthquake has been carefully analyzed by several studies (e.g., Ihmlé, 1998; Kikuchi & Kanamori, 1994; Zhan et al., 2014) that retrieved moment rate functions very similar to the SCARDEC STF shown in Figure 3c. Additionally, the source extent has been found very compact (about 40 km × 40 km) and the point source approach is therefore expected to be precise for all the stations considered here. We use the SCARDEC focal mechanism and STF (Table 1 and Figure 3c), and follow the procedure described in section 4 to model the waveforms. Simulations (black lines in Figure 5b) well reproduce the clear positive signal at LPAZ and confirm that the expected signal is below the noise of the stations at NNA and YUNZ.
4 Additional Observations by Waveform Stacking
4.1 Stacking Strategies
(1)
are the observed acceleration, its predicted waveform, and the P travel time at sensor i, respectively. The reference time for ai and si is chosen as the origin time of the event, such that So is defined for negative times and So(0) is the weighted sum of the individual PEGS at their respective P arrival time. Equation 1 assumes that the theoretical signals si in P reference time only differ by a proportionality constant (i.e., 
, where 
is a waveform independent of the sensor location). 
is the variance of the zero mean observed waveform ai in the pre-event time window (thus associated with a “noise” estimation):
(2)
(3)
(4)
will be computed together with a simple stack Ss, defined as
(5)
(by replacing So by Ss in equations 3 and 4), for comparison with the detection gain provided by the optimal stack strategy.
4.2 The 23 January 2018 Mw = 7.9 Gulf of Alaska Strike-Slip Earthquake
4.2.1 Data and Individual Station Modeling
The 23 January 2018 Mw = 7.9 Gulf of Alaska earthquake is a large shallow strike-slip earthquake (Figure 6a) recorded by a very dense regional network of broadband sensors. The recent deployment of the stations of the USArray project across Alaska and Northwest Canada, together with the preexisting broadband stations, result in more than 300 instrumented locations in the 0- to 1,800-km epicentral distance range. We process all data as in section 3 and keep the sensors where σi (equation 2) computed in the 800 s preceding origin time is below 5 nm/s2. This loose selection criterion is intended to eliminate the sensors that may not work properly, without removing any waveforms that can significantly contribute to the stack (taking into account that the best sensors have σi≃ 0.3 nm/s2, a sensor with σi= 5 nm/s2 is down-weighted by a factor ≃ 300 according to equation 1). After this selection, we now consider the data of 262 sensors (Figure 6), among which 167 belong to the USArray Transportable Array (network code TA), 64 to the Alaska Regional Network (network code AK), and the remaining ones to the AT, AV, US, II, IU, and CN networks.
coefficient: The 13 stations with the highest values are filled with yellow, the 165 stations (168 sensors) with the lowest values are filled with blue, and the remaining 79 stations (81 sensors) are filled with orange. (b) Observed (red) and modeled (black) pre-P vertical accelerations of the 13 sensors with best observation potential (filled with yellow in (a)). Waveforms are filtered in the 0.002- to 0.03-Hz frequency range and represented in a time window starting 5 min before the earthquake origin time and terminating at the P wave arrival time at each station (1 nm/s2 scale is shown to the right). Modeling is based on a double-couple point source simulation using focal mechanism and STF from SCARDEC (Table 1 and Figure 3d). Network, name, azimuth θ, and epicentral distance (in kilometers, following Earth's surface) of the stations are shown to the left of each signal. PEGS = prompt elastogravity signals; GCMT = Global Centroid Moment Tensor; STF = source time function.
Due to the magnitude of the earthquake, we do not expect the individual PEGS amplitudes to be directly detectable at any of the stations. The GCMT-expected amplitudes (Figure 6a, using the GCMT parameters of Table 1 and Figure 3d) confirm that the maximum values remain below 0.5 nm/s2. All the individual predicted waveforms si have also been modeled using the SCARDEC STF (Figure 3d) and focal mechanism (Table 1), following the methodology of section 4. After ordering the sensors by decreasing values of the coefficient 
, we show in Figure 6b the observations and simulations for the 13 first sensors. No conclusion can be drawn from these individual observations at the best potential stations, but the good prediction of the waveforms polarity at the P wave arrival encourages the use of stacking approaches.
4.2.2 Stacking Results
We first apply the optimal stack strategy (equation 1) to the 262 sensors. As shown in Figure S6, the PEGS now clearly appear, with 
reaching a SNR value of 7.4 at t = 0, corresponding to the P arrival time at each station. Under the hypotheses mentioned in section 13, equation 1 predicts that SNR value is optimized when the maximum number of sensors are used. However, these hypotheses tend to be less and less respected when expected PEGS |si| are low and pre-event noises σi are high. Observation of signals with large σi indeed reveals that high pre-event noise is frequently associated with linear trends or large oscillations at long period, which break the hypotheses that the noise is stationary and gaussian. On the other hand, stations with low expected |si| are along the nodal planes with azimuths ≃ 75° and ≃ 345°. The amplitudes 
are difficult to predict accurately there, as they are sensitive to the details of the source mechanism and of the Earth structure.
We therefore order the sensors by decreasing values of the coefficient 
and explore the optimal number N of sensors to be considered in the stack. We find that the SNR increases with the number of the sensors up to N≃90–100, where it reaches a value of 10, before slowly decreasing to the 7.4 value obtained when considering all the sensors. We thus show in Figure 7a the optimal stack with the first 94 sensors, whose locations are shown in Figure 6a by the yellow and orange triangles. Besides the high SNR of 10, Figure 7a also shows that the expected stack signal (by replacing ai by si in equation 1) is very close to the observed one. PEGS are also clearly detected when using the normalized simple stack 
applied to the same 94 sensors: Figure 7b shows that the SNR value is only slightly lower compared to the optimal stack (9 instead of 10). In acceleration units, the average value of the PEGS absolute amplitudes at the P wave arrival time Ss(0) peaks to 0.37 nm/s2, close to the predicted value of 0.39 nm/s2. Unambiguous detection at such low acceleration values cannot be reached with any single broadband seismometer today and thus illustrates the power of the stacking strategies when applied to a dense and high quality seismic network.
(yellow and orange triangles in Figure 6). The resulting observed stack 
(red line) reaches a SNR value of 10 at t = 0; 
is compared with the modeled stack (black line) where PEGS synthetics are summed instead of data. (b) Simple stack strategy, using equation 5 with the same 94 sensors as in (a) (red line). Acceleration scale (to the right) and SNR scale (to the left) refer to Ss and 
, respectively. At t = 0, 
peaks to 9 and Ss peaks to 0.37 nm/s2, the latter value being very similar to the one predicted by the modeled simple stack (black line). PEGS = prompt elastogravity signals; SNR = signal-to-noise ratio.
4.3 The 27 February 2010 Mw = 8.8 Maule Megathrust Earthquake
4.3.1 Data and Individual Station Modeling
The 27 February 2010 Mw = 8.8 Maule megathrust earthquake (Figure 8) is the third largest earthquake ever recorded by modern broadband seismometers (after the 11 March 2011 Tohoku earthquake and the 26 December 2004 Sumatra earthquake). We retrieved all the broadband data publicly available at distances closer than 2,100 km from earthquake epicenter and processed them as in section 3. For the same reasons as for the 23 January 2018 Alaska earthquake, we then keep the sensors where σi (equation 2) computed in the 600-s preceding origin time is below 5 nm/s2. After this selection, we now consider the data of 31 sensors (Figure 8), among which 13 belong to the IPOC network (Integrated Plate Boundary Observatory Chile, network code CX), 10 to the temporary network of the “Lithospheric Structure above the variably dipping Nazca Slab” project (network code XH), and the remaining ones to the IU, GT, II, and C networks.
coefficient: the 8 stations with the highest values are filled with yellow, and the remaining 21 stations (23 sensors) are filled with orange. (b) Observed (red) and modeled (black) pre-P vertical accelerations of the eight sensors with best observation potential (filled with yellow in (a)). Waveforms are filtered in the 0.002- to 0.03-Hz frequency range, and represented in a time window starting 5 min before the earthquake origin time and terminating at the P wave arrival time at each station (2-nm/s2 scale is shown to the right). Modeling is based on a double-couple point source simulation using focal mechanism and STF from SCARDEC (Table 1 and Figure 3e). Network, name, azimuth θ, and epicentral distance (in kilometers, following Earth's surface) of the stations are shown to the left of each signal. PEGS = prompt elastogravity signals; GCMT = Global Centroid Moment Tensor; STF = source time function.
We show in Figure 8a the expected PEGS amplitudes when using the GCMT source parameters and STF (Table 1 and Figure 3e). As predicted for this mechanism type (Figure 1), PEGS are maximal in the direction orthogonal to the trench and at epicentral distances of the order of 1,000 km. However, none of the sensors located in the vicinity of this area (TRQA and stations from the temporary XH network) are optimal in terms of pre-event noise: at TRQA and at the four quietest XH stations (ESP10, ES07, ESP03, and ESP12), σi ranges between 0.5 and 0.7 nm/s2. Such values at the high quality TRQA primary sensor are related to the occurrence, 10 hr before the Maule earthquake, of a Mw = 7 earthquake in the Japan region. This large-magnitude earthquake excited the Earth for many hours through multiple Earth-cycling surface waves (R7 Rayleigh waves are expected to arrive in the Chile area at the time of the Maule earthquake). This unfavorable configuration implies that PEGS are unlikely to be above ≃ 3σi at any of the stations, and their individual detections are thus expected to be difficult.
These expectations are confirmed when using the SCARDEC STF (Figure 3e) and focal mechanism (Table 1), and modeling the individual waveforms with the methodology of section 4. After ordering the sensors by decreasing values of the coefficient 
, we show in Figure 8b the observations and simulations for the eight first sensors. As in the Alaska earthquake case, polarities observed at the P arrival time are generally consistent, but pre-event noise prevents from directly using these individual observations.
4.3.2 Stacking Results
We show in Figure 9a the optimal stack for all the 31 considered sensors. At t = 0, 
peaks at a SNR value of 7.6, meaning that the PEGS of the Maule earthquake are unambiguously detected. Figure 9a further shows that the expected stack signal (black line) is close to the observed one, with a SNR value peaking at 8.2. Contrary to the Alaska earthquake case, 
does not peak to a higher SNR value when waveforms with the lowest values of the 
coefficient are removed from the stack. This can be understood by the fact that most of the latter waveforms were recorded by high-quality sensors located in northern Chile (LVC and the stations belonging to the CX network). Such waveforms with low σi are much less prone to have a non-Gaussian and/or nonstationary pre-event noise, and due to their along-trench position, the PEGS amplitudes 
are small but do not depend critically on their exact position (as it was the case for the nodal stations of the Alaska earthquake).
(red line) reaches a SNR value of 7.6 at t = 0; 
is compared with the modeled stack (black line) where PEGS synthetics are summed instead of data. (b) Simple stack strategy (red line), using equation 5 with the best eight sensors in terms of coefficient 
(yellow triangles in Figure 8). Acceleration scale (to the right) and SNR scale (to the left) refer to Ss and 
, respectively. At t = 0, 
peaks to 6.7 and Ss peaks to 1.4 nm/s2, the latter value being similar to the one predicted by the modeled simple stack (black line). PEGS = prompt elastogravity signals; SNR = signal-to-noise ratio.
PEGS are also clearly detected with the normalized simple stack 
applied to the eight best potential sensors shown in Figure 8b. As shown in Figure 9b, the SNR value is lower but remains at a very significant value of 6.7 at t = 0. In acceleration units, the average value of the PEGS absolute amplitudes at the P wave arrival time Ss(0) peaks to 1.4 nm/s2, close to the predicted value of 1.3 nm/s2. Such acceleration values would have been very likely individually detected, in particular at TRQA station, without the occurrence of the preceding Mw = 7 Japan earthquake. Stacking strategies are here shown to reveal PEGS for a large earthquake recorded in a nonoptimal configuration, both in terms of station distribution and pre-event noise.
5 Discussion and Conclusions
5.1 Undetected Large Earthquakes
Based on simulations and observations shown in the previous sections, we might expect two additional earthquakes to have detectable PEGS: the 26 December 2004 Mw = 9.1–9.2 Sumatra megathrust earthquake and the 24 May 2013 Mw = 8.3 Okhotsk deep earthquake.
The Okhotsk earthquake is the largest deep earthquake recorded by modern broadband seismometers but was preceded 12 hr before by a Mw = 7.4 earthquake in the Fiji Islands. This configuration leads to even higher pre-event noise than for the Maule earthquake: for the two best located sensors, PET in Kamchatka and MA2 in Eastern Russia, σi is around 0.9 nm/s2 (while it is close to 0.2 nm/s2 in a quiet period). This means that even strong PEGS amplitudes of the order of 2 nm/s2 cannot be reliably detected at these sensors. Additionally, as the available station distribution is sparse in this area, we do not expect stacking approaches to be efficient for this earthquake.
Even though it shares similar magnitude and mechanism type with the 11 March 2011 Tohoku earthquake, the 2004 Sumatra earthquake, however, has a very different rupture time history. Its total duration is 3 to 4 times longer (≃ 500 s compared to ≃ 150 s; Guilbert et al., 2005; Krüger & Ohrnberger, 2005; Ni et al., 2005) and even more importantly, its moment rate remained low for the first 70 s of the rupture (Ammon et al., 2005, 2006). At that time, moment rate function indicates that the accumulated seismic moment was about 3.1021 N·m (Mw = 8.3), 10 times smaller than the Tohoku earthquake (Mw = 8.9). For epicentral distances in the range 1,000–2,000 km, which were shown to be optimal in the case of the Tohoku earthquake (Vallée et al., 2017), measurable signals would therefore arrive after the arrival of the elastic waves and cannot be detected. At further distances, the moment rate history is less critical, but signals become close to the noise level even at the quietest stations.
5.2 Summary of the Observations
We show in this study that detection of PEGS with broadband seismometers is not restricted to exceptional (Mw > 9) earthquakes. PEGS of much smaller earthquakes can be detected with individual instruments when their source parameters (mechanism, depth, and moment rate evolution) as well as the recording configuration (well located high-quality sensors and low seismic noise level) are favorable. Our simulations show that for a given magnitude, the optimal earthquake type is a shallow strike-slip earthquake with a fast moment rate evolution. In contrast, thrust earthquakes on shallow-dipping interfaces and slow moment rate evolution are the most difficult to detect. PEGS are consistently clearly observed for the 11 April 2012 Mw = 8.6 Wharton Basin earthquake (strike slip) while they are undetected for the 26 December 2004 Mw = 9.1–9.2 Sumatra megathrust earthquake. If considering only their focal mechanism, deep earthquakes generate slightly weaker PEGS amplitudes than shallow strike-slip earthquakes. However this lower detection potential is usually compensated by the faster moment rate evolution of deep earthquakes, explaining why the PEGS of both the 19 August 2018 Mw = 8.2 Fiji and the 9 June 1994 Mw = 8.2 Bolivia earthquakes are successfully observed by regional broadband seismic sensors. For all these detected earthquakes, observed PEGS waveforms are in good agreement with modeled PEGS provided the appropriate source parameters (location, focal mechanism, and moment rate function) are used. This simultaneously supports the observations, confirms the accuracy of the modeling approaches described in Vallée et al. (2017) and Juhel, Montagner, et al. (2018), and shows the PEGS sensitivity to earthquake characteristics.
Accurate modeling is also useful for PEGS detection when the searched signals are close to the instrumental or seismic noise at a number of sensors. Waveforms can then be stacked in an optimal way, taking into account both the intrinsic sensitivity of each sensor and the expected amplitudes at each location. When such an approach is applied to a dense and high quality seismic network, the PEGS detection threshold is largely improved: we show that thanks to the recent efforts made to instrument the whole Alaska region, the 23 January 2018 Mw = 7.9 Gulf of Alaska earthquake is observed with a high SNR, meaning that even smaller magnitude events could have been detected in this area. Stacking strategies are finally shown to be efficient to extract the PEGS of large earthquakes recorded in less optimal recording configurations: waveform stacking reveals the PEGS of the 27 February 2010 Mw = 8.8 Maule earthquake, in spite of high pre-event noise (due to the recent occurrence of another large earthquake) and much sparser station distribution than in the Alaska earthquake case. As a whole, this study shows PEGS observations for five earthquakes in the magnitude range 7.9–8.8. These observations are summarized in Figure 10, in which original regional observations made for the 2011 Tohoku earthquake (Vallée et al., 2017) are also mentioned.
5.3 Near-Future Potential of PEGS for Earthquake Monitoring
These foregoing observations show that PEGS can be detected on a regular basis with broadband seismometers, even with network configurations not specifically designed to record them. In a rapid response perspective, PEGS also have the advantage to be recorded at regional distances far from the largest ground motions, therefore avoiding potential difficulties related to sensors saturation, instrument damages, or data transmission gaps. Since we also have documented their intrinsic sensitivity to key source parameters (focal mechanism and moment rate function), the PEGS potential for earthquake monitoring applications is concretely demonstrated. We discuss here how these signals can be used in this respect in the near future, even without considering the dramatic progress (Juhel, Ampuero, et al., 2018; McManus et al., 2018) that could be provided by next generation instruments currently explored in the gravitational-wave physics community (Ando et al., 2010; Harms et al., 2013; McManus et al., 2017; Paik et al., 2016).
A straightforward implementation for regional seismic networks is to extract the time window preceding each P wave arrival (manually or automatically picked). Vertical acceleration at the P arrival time can then be measured in the 0.002- to 0.03-Hz frequency range, and a SNR can be evaluated in real-time by using the amplitude levels of earlier time windows (that have been previously computed and stored). Obtaining several sensors with significant SNR is a direct and early evidence of the occurrence of a very large earthquake (Mw > 8). Together with real-time local GNSS measurements (Allen & Ziv, 2011; Blewitt et al., 2009), PEGS therefore offer a standard way to detect the presence of a large growing earthquake in the minutes following its initiation. At this very early stage, a simple observation can also be very informative for tsunami hazard assessment: the existence of even one signal with significant positive PEGS proves that the event is not a subduction megathrust earthquake (for which PEGS are negative everywhere).
More ambitious goals can be pursued once the earthquake has been located, even approximately. In well-instrumented areas, this does not slow down the use of PEGS, because earthquake location can be obtained just after local stations were triggered, and thus before the P wave arrival at the stations with best PEGS potential (located at regional distances, further than ≃ 400 km). In a similar way as done today with seismic waves or static displacements, a source model can then be obtained by optimizing the agreement between data and synthetics. In a near-real time perspective (less than a few minutes after earthquake origin time) such source characterization can theoretically only be provided by modeling local seismic and geodetic data (Allen & Ziv, 2011; Blewitt et al., 2009; Delouis, 2014; Li et al., 2013) or regional PEGS. But modeling the local data of a large earthquake requires to invert for many source parameters because the extended nature of the seismic process cannot be neglected; on the contrary, regional PEGS modeling benefits from the approximate validity of the point source representation. A robust magnitude estimation can therefore be obtained by a standard source parameterization, simply involving the focal mechanism and a simplified moment-rate evolution.
The technical implementation of the associated inverse problem is out of the scope of the present study, but several important insights about the concrete PEGS potential for seismic early warning can be inferred from scenarios where earthquake location and focal mechanism are known. In Figure 11, we show the example of shallow strike-slip earthquakes, for magnitudes representative of large events expected on major transform boundaries (such as the San Andreas or North Anatolian faults). For each of the two considered magnitude scenarios (Mw= 7.8 and Mw= 8), three source durations are considered: The central value comes from the empirical magnitude-duration relation determined by Houston (2001), and the two extreme ones use the standard error values of the lognormal duration distribution observed by Courboulex et al. (2016). For each of these six cases, Figure 11 shows both the P wave arrival time and the corresponding PEGS amplitudes, therefore providing the information whether and when a PEGS-derived magnitude can be obtained. Based on the observations for the Gulf of Alaska earthquake and assuming a seismic network with similar density and quality, the detection threshold is considered to be reached if PEGS amplitudes are in some places above ≃ 0.4 nm/s2 (dashed yellow isocontour). In terms of early warning, our simulations show that key areas are then located along the P and T axes, at epicentral distances between 450 and 700 km from the earthquake, corresponding to P arrival times between 60 and 90 s. At a time of 90 s after the earthquake origin, all the broadband signals located at distances closer than 700 km can therefore be optimally stacked, in order to reveal the PEGS. As shown by Figure 11a), a positive detection at this very early time is a direct evidence of an earthquake with magnitude larger than 7.8. Lower magnitudes can indeed be ruled out, because only the shortest duration Mw = 7.8 scenario generates areas with PEGS amplitudes above ≃ 0.4 nm/s2. This directly evidences how PEGS information is complementary with classical approaches relying on local seismic waves, for which larger magnitudes make early determinations increasingly challenging. It has to be mentioned, however, that a negative PEGS detection does not prove that the earthquake has a magnitude smaller than 8, as illustrated by the longest duration Mw = 8 scenario, where the 0.4 nm/s2 level is not reached. At this postorigin time of 90 s, Figures 11a and 11e finally show that the PEGS amplitude pattern can be similar between a short Mw = 7.8 earthquake and a classical Mw = 8 earthquake. This ambiguity however disappears with time and distance and an exact magnitude determination, rather than only a lower bound, is expected with the progressive integration of additional sensors (in the present case, up to distances of about 1,100 km, which requires to wait 135 s after earthquake origin time).
, where σT = 0.35 is the empirical standard error of the lognormal duration distribution observed by Courboulex et al. (2016). PEGS are simulated using the normal-mode approach using a point source at 20-km depth with an isosceles triangular moment rate function. PEGS amplitudes at P wave arrival time are both color coded and shown by the two isovalue contours at 0.2 nm/s2 (dash-dotted yellow line) and 0.4 nm/s2 (dashed yellow line). P wave arrival times are themselves shown by isovalue contours (dotted black lines). For each quadrant, the white cross indicates the location of the maximal PEGS amplitude. PEGS = prompt elastogravity signals.
Multiple PEGS observations provided in this study, together with these last concrete scenarios for strike-slip earthquakes, demonstrate the interest of integrating these signals in earthquake monitoring protocols. In all the areas threatened by Mw > 7.5 strike-slip earthquakes or Mw > 8.5 megathrust earthquakes, PEGS detection is a direct evidence that such earthquakes have just happened. And even if only rare events will be detected, their early and reliable quantification is a key information for rapid response and tsunami warning. PEGS-based monitoring approaches could today be implemented in a lot of earthquake-prone areas but with a detection potential that may suffer from the sparsity of the seismic networks at regional distances. For example, the monitoring of the South America megathrust earthquakes would benefit from a denser instrumentation of the sub-Andean domain, and future events along the San Andreas fault would be better recorded with additional stations in Nevada and Arizona. We highlight that the existence of these signals is an additional motivation to instrument locations at 500–1,500 km from the potential sources, as it is today the case in Alaska.
Acronyms
-
- PEGS
-
- prompt elastogravity signals
Acknowledgments
This study was inspired by original theoretical and observational works on PEGS, in which Jean Paul Ampuero, Jean-Paul Montagner, Matteo Barsuglia, Jan Harms, and Pascal Bernard were particularly involved. We are grateful to Olivier Coutant for upgrading the AXITRA code, allowing for waveform simulation in a multisource and multireceiver configuration. Robert Busby, Kasey Aderhold, and Jean Paul Ampuero incited us to look at the data of the 2018 Gulf of Alaska earthquake. Jean Paul Ampuero, Robert Busby, Jean-Paul Montagner, and Susan Beck provided useful insights on the manuscript or on the data used in this study. We thank both reviewers for their comments and in particular Barbara Romanowicz for her advice to further illustrate the PEGS potential for early warning. Such an observational study deeply relies on engineering and technical efforts made to install and maintain high-quality broadband seismic sensors: The following seismic networks specifically contributed to the PEGS observations shown in this study (alphabetic order): Alaska Regional Network (AK, https://doi.org/10.7914/SN/AK), National Tsunami Warning Center Alaska Seismic Network (AT, https://doi.org/10.7914/SN/AT), Chilean National Seismic Network (C), Plate Boundary Observatory Network Northern Chile (CX, https://doi.org/10.14470/PK615318), GEOSCOPE (G, https://doi.org/10.18715/GEOSCOPE.G), GEOFON (GE, https://doi.org/10.14470/TR560404), Global Telemetered Seismograph Network (GT, https://doi.org/10.7914/SN/GT), IRIS/IDA Seismic Network (II, https://doi.org/10.7914/SN/II), Global Seismograph Network (IU, https://doi.org/10.7914/SN/IU), Malaysian National Seismic Network (MY), Pacific21 network (PS), USArray Transportable Array (TA, https://doi.org/10.7914/SN/TA), United States National Seismic Network (US, https://doi.org/10.7914/SN/US), Broadband Study of the Altiplano and Central Andes network (XE, https://doi.org/10.7914/SN/XE_1994), and Argentinean Slab network (XH, https://doi.org/10.7914/SN/XH_2008). We also thank the latter seismic networks for the public distribution of their data that were retrieved through the IRIS-DMC (http://ds.iris.edu/ds/nodes/dmc/), GEOFON (https://geofon.gfz-potsdam.de/waveform/), and IPGP(http://datacenter.ipgp.fr/) data centers. The SAC (http://ds.iris.edu/ds/nodes/dmc/software/downloads/sac/) free software was used for data processing. Most numerical computations were performed on the S-CAPAD platform, IPGP, France. We acknowledge the financial support from the UnivEarthS Labex program at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02), from the Agence Nationale de la Recherche (grant ANR-14-CE03-0014-01), and from the GEOSCOPE program (funded by CNRS-INSU and IPGP).
