Determination of the source parameters of the 2011 Tohoku-Oki earthquake from three-component pre-P gravity signals recorded by dense arrays in Japan

Dynamic earthquake rupture causes mass redistribution around the fault, and the emitted propagating seismic waves are accompanied by bulk density perturbations. Both processes cause transient gravity changes prior to the arrival of P-waves. Such pre-P gravity signals have been detected in previous studies of several large earthquakes. However, the detections were limited to the vertical component of the signal owing to the high noise level in the horizontal records. In this study, we analyzed dense tiltmeter array data in Japan to search for the horizontal components of the signal from the 2011 Mw 9.1 Tohoku-Oki earthquake. Based on the synthetic waveforms computed for a realistic Earth model, we stacked the horizontal records and identified a signal that evidently exceeded the noise level. We further performed a waveform inversion analysis to estimate the source parameters. The horizontal tiltmeter data, combined with the vertical component of the broadband seismometer array data, yielded a constraint on the dip angle and magnitude of the earthquake in the ranges of 11.5°–15.3° and 8.75°–8.92°, respectively. Our results indicate that the analysis of the three components of the pre-P gravity signal avoids the intrinsic trade-off problem between the dip angle and seismic moment in determining the source mechanism of shallow earthquakes. Pre-P gravity signals open a new observation window for earthquake source studies.


Introduction
Earthquake faulting causes mass redistribution in the vicinity of the fault, which results in a post-earthquake static gravity change. Gravity changes have been observed in large earthquakes and compared with model predictions (Imanishi et al. 2004;Matsuo and Heki 2011;Wang et al. 2012a, b). During dynamic rupture, bulk density perturbations are temporarily caused by propagating compressional waves. Both the mass redistribution due to faulting and the propagating compressional waves produce a transient gravity change, which commences before the P-wave arrival (Harms et al. 2015). Moreover, gravity change δg induces acceleration motion ü of the observation site; therefore, the observed acceleration s is the difference between δg and ü (Heaton 2017;Vallée et al. 2017;Kimura 2018;Kame and Kimura 2019;Kame 2021). The "pre-P gravity signal" refers to observable acceleration change s equivalent to ü − δg. "Pre-P gravity signal" was named "prompt elastogravity signal" by Vallée et al. (2017).
To date, the detection of the pre-P gravity signal has been reported in several large earthquakes, most notably in the 2011 Mw 9.1 Tohoku-Oki earthquake in Japan (Montagner et al. 2016;Vallée et al. 2017;Kimura et al. 2019a;Vallée and Juhel 2019). However, the detections were limited to the vertical component, even though the pre-P gravity signal is a three-component vector. The horizontal records are generally noisier than the vertical records (Peterson 1993); hence, signal detection remains challenging in the horizontal components. In fact, Kimura et al. (2019a) analyzed the tiltmeter records of the High Sensitivity Seismograph Network in Japan (Hinet), which works as horizontal accelerometer array, but failed to detect a signal. Their failure may also be attributed to the theoretical model used for waveform stacking analysis. They selected the data and reversed the polarity based on the synthetic amplitudes expected from the pioneer model of Harms et al. (2015), which only simulated gravity change δg in an infinite homogeneous non-selfgravitating medium.
Recent progress in the theoretical modeling of pre-P gravity signals has enabled the synthesis of signal waveforms for realistic Earth structure models. Zhang et al. (2020) developed a numerically stable and low-computation-cost method to solve the fully coupled elastogravity equations, i.e., the equation of motion of the continuum and Poisson's equation of the gravity change, for a spherically symmetric Earth. The method was based on Wang et al. (2017) and was developed to avoid numerical artifacts originating from the large amplitude of direct P-waves . Currently, Zhang et al. (2020) have proposed the most sophisticated methods for synthesizing pre-P gravity signals (Kame 2021). Zhang et al. (2020) also determined source parameters using their signal simulation method. They synthesized pre-P waveforms for a point source located at the hypocenter of the 2011 Tohoku-Oki earthquake, determining the source parameters that best fit the vertical broadband seismometer data observed in and around Japan. The estimated parameters were consistent with seismologically estimated values, such as the centroid moment tensor (CMT) solution, and the strike and rake angles were robustly determined (Fig. 7 of Zhang et al. 2020). However, relatively large uncertainties remained in the moment magnitude and dip angle, similar to the wellknown trade-off problem between the two parameters determined from long-period seismic waves (Kanamori and Given 1981).
In this study, we analyzed the Hi-net tiltmeter array data to detect the horizontal components of the pre-P gravity signal from the 2011 Tohoku-Oki event. We stacked the observed waveforms based on the synthetic waveforms using Zhang et al. (2020). The tiltmeter array data and vertical component of the F-net (Full Range Seismograph Network of Japan) broadband seismometer array data were inverted for point-source parameters. The sensitivities of the horizontal and vertical components to the dip angle and seismic moment were investigated. As a result, the parameters were well constrained to fit the observed three-component data.

Data
The Tohoku-Oki earthquake (Mw 9.1) occurred off the coast of northern Honshu, Japan, on March 11, 2011. We retrieved two-component borehole tiltmeter records at 689 Hi-net stations and vertical broadband seismometer records at 70 F-net stations (Obara et al. 2005;Aoi et al. 2020). The station locations are shown in Fig. 1a-c. The tilt data provided by the National Research Institute for Earth Science and Disaster Resilience (NIED) are in radians. They were converted into the horizontal acceleration in m/s 2 through multiplication with the gravitational acceleration (− 9.8 m/s 2 ). The installation orientation of the tiltmeters was corrected at each station according to Shiomi et al. (2003) and Shiomi (2012). The linear trend was removed based on the slope 30 min before the event origin time. The horizontal noise spectrum is presented in Additional file 1: Fig. S1 and is a typical spectrum of an ambient seismic field with two peaks; the primary microseism occurs at approximately 0.07 Hz and secondary at approximately 0.3 Hz (e.g., Fig. 1 of Nishida 2017). The 0.002-Hz two-pole high-pass and 0.05-Hz sixpole low-pass causal Butterworth filters were applied to eliminate the lower tidal noise and higher microseisms. We followed the filter in Vallée et al. (2017), and the difference was only the low-pass corner frequency 0.05 Hz in the analysis of horizontal records. The frequency is shown in Additional file 1: Fig. S1. The horizontal signal processing is in the time domain and is free from acausal artifacts. The broadband seismometer records were processed fully following Vallée et al. (2017), and the data were terminated at visually picked P-wave arrival time t obs P and deconvolved from the instrumental acceleration response. Then, the 0.002-Hz two-pole high-pass and 0.03-Hz six-pole low-pass causal Butterworth filters were applied to eliminate tidal and microseismic noise, as well as artifacts originating from the frequency-domain deconvolution (Kimura et al. 2019b).

Waveform simulation
We used the QSSPPEGS code developed by Zhang et al. (2020), which was modified from the QSSP code (Wang et al. 2017), to synthesize the pre-P gravity signal waveforms. QSSPPEGS numerically solves the fully coupled elastogravity equations for a spherically symmetric Earth model: and calculates pre-P ground acceleration ü and gravity change δg(= ∇ψ) . Here, ρ = ρ(r) is the unperturbed density; u = u(r, θ , φ, t) is the particle displacement vector; σ = σ (r, θ , φ, t) is the incremental stress tensor in the Lagrangian description; ψ = ψ(r, θ, φ, t) is the incremental gravity potential in the Eulerian description; g = g(r) is the undisturbed static Earth gravity (downward positive); f = f (r, θ, φ, t) represents the seismic source; G is the gravitational constant; r, θ, φ are the radial distance, polar angle, and azimuthal angle with the origin at the center of the spherically symmetric Earth, respectively; and e r is the radial unit vector. The pre-P gravity signal measured by ground-based sensors was synthesized as s =ü − δg . We ignored the tilt terms −g ∂u z ∂x and − g ∂u z ∂y for the horizontal components s x and s y , and free-air term u z dg dz for the vertical component s z due to their insignificant contributions to the signal amplitudes Zhang et al. 2020). Here, the x, y, and z axes were considered as eastward, northward, and upward positives, respectively.
As detailed in "Detection of the horizontal components of the pre-P gravity signal" section, a point source was adopted for the waveform simulation of the 2011 Tohoku-Oki earthquake. It was located at the hypocenter (latitude = 38.19°N, longitude = 142.68°E, depth = 21 km) determined by Chu et al. (2011). Event origin time t eq , seismic moment M 0 , rupture duration T , and source mechanism (strike, dip, rake) were set to 05:46:23 UTC, 5.31 × 10 22 Nm, 140 s, and (203°, 10°, 88°), respectively (Global Centroid Moment Tensor: GCMT; Ekström et al. 2012). The moment rate function Ṁ (t) is described as a squared half-period sinusoidal function: A modified Earth structure model based on AK135 by Wei et al. (2012) was used, and the Green's functions were calculated to a frequency of 0.25 Hz. The synthesized signal waveforms were truncated at the P-wave arrival time t P .
is the smaller of the visually picked value t obs P and t theo P (2 s before the theoretical value calculated using the TauP toolkit; Crotwell et al. 1999). In "Waveform inversion" section, the dip angle and M 0 were chosen as the model parameters that could not be well constrained from vertical records (Zhang et al. 2020). The rupture duration T was assumed to obey the scaling law (Kanamori and Anderson 1975).

Point source correction
The pre-P gravity signals calculated for a point source tend to be larger than those calculated for a finite fault model with the same total seismic moment (Fig. S6 of Zhang et al. 2020). Additional file 1: Figure S2 shows that this overestimation is particularly prominent in the horizontal components, which affects the inversion results. To evaluate further acceptable signal amplitudes, in "Waveform inversion" section, we employed an empirical constant factor (estimated roughly to be 0.75 as in Additional file 1: Fig.  S2), which was multiplied by the synthesized waveforms in the horizontal components.

Waveform stacking
The amplitude of the pre-P gravity signal is comparable to or smaller than the background noise level To reduce the noise and enhance the signal, we used the waveform stacking approach (Kimura et al. 2019a;Vallée and Juhel 2019). Because the amplitudes of the pre-P gravity signals are in many cases expected to monotonically increase until the P-wave arrives ( Fig. 1 of Zhang et al. 2020), the waveforms are aligned and stacked with t P as: where s i (t) , a i (t) , and t i P are the synthetic acceleration, observed acceleration, and P-wave arrival time at the i th sensor, respectively; a(t) is the stacked trace; N is the total number of stacked sensors; and sgn( * ) is the sign function for the polarization reversal. The synthetic waveforms can also be stacked by replacing a i (t) with s i (t).

Inversion method
For the waveform inversion, we used the misfit function defined as where p is the variable parameter to be estimated, σ i is the standard deviation of a i (t) before t eq (here, we defined a time window t eq − 10min, t eq ), and s i (t; p) is s i (t) calculated using the parameter values of p . Because the amplitudes of pre-P gravity signals increase just before the P-wave arrival, we used the last quarter of the time window between t eq and t i p , i.e., t i 2 = t i P and t i 1 = t i P − 1 4 t i P − t eq . We defined the normalization constant as α = N i=1 t i 2 t i 1 dt such that R(p) reached unity if the noise level was time-invariant and the signal was completely fitted with the parameter values of p.

Detection of the horizontal components of the pre-P gravity signal
The horizontal components of the pre-P gravity signal from the 2011 Tohoku-Oki earthquake are expected to have an amplitude of approximately 1-2 nm/s 2 over the Japanese land area with a range of 0.002-0.05 Hz (Fig. 1a,  b). The tiltmeter data were stacked based on the synthetic polarities (Eq. 3). Here, from 1378 records (two components at 689 stations), 429 records were selected that met the following two criteria: (i) records with the ratio of the synthetic amplitude s i t i P to the noise level σ i greater than 0.7, and (ii) records that did not exceed 10 nm/s 2 30 min before t eq . The locations of the selected stations are shown in Fig. 1a As a result of stacking, a clear signal above the noise level was observed before the P-wave arrival (Fig. 2a). In the stacked trace, the range of origin times varying per station is indicated. For a station with t p , its origin time appears at −t p and for all stations origin times are in the range of [ −t max p , −t min p ], where t max p and t min p are the largest and smallest t p among the stations, respectively. The signal amplitude was 0.64 nm/s 2 , which was 9.5 times higher than the reduced noise level of 0.067 nm/ ]. The identified signal shows a characteristic duration of 1 min, which appears to reflect the window lengths [ −t p , 0 ] whose ends were masked by the large-amplitude P-wave arrivals, rather than the rupture duration, as theoretical syntheses show signal durations longer than 1 min away from Japan (Fig. 1 of Zhang et al. 2020).
In addition, the radiation pattern in the north-south (NS) component was confirmed; the NS-component records were stacked for the positive and negative synthetic value regions, respectively. In both cases, signals with signs identical to those of the synthesized waveforms were found to exceed the noise level (Additional file 1: Fig. S3).
However, the amplitude of the detected signal was less than half that of the stacked synthetic waveforms (1.65 nm/s 2 ), and their difference of approximately 1 nm/ s 2 exceeded the noise level. This result contrasts that of the vertical records (Fig. 2b). Among the 70 broadband seismometer records, 43 records with a ratio of the synthetic amplitude s i t i P (shown in Fig. 1c) to the noise level σ i greater than unity were selected and stacked. The stacked synthetic amplitude of 0.75 nm/s 2 was close to the signal amplitude of 0.68 nm/s 2 , and the difference was within the noise level range of 0.050 nm/s 2 .

Waveform inversion
We performed a waveform inversion analysis to resolve the inconsistency between the synthetic and observed records, as revealed in the previous section. We estimated the appropriate source parameters that fit the observations. Figure 3 shows the distribution of the misfit function for the dip angle and Mw in the range of 5°-25° and 8.4°-9.3°, respectively. The misfit tendencies differed between the horizontal and vertical components; the misfit function for the horizontal components was sensitive to the dip angle, whereas that for the vertical component exhibited an inverse proportionality between the dip angle and Mw. The sum of the misfit functions for the three components was minimal at a dip angle of 13° and Mw of 8.85 (Fig. 3c). Here, all 1378 tiltmeter records and 70 vertical broadband seismometer records were used to calculate the misfit functions. Note that the contribution of each record to the misfit functions depends on the noise level σ i (Eq. 4). Subsequently, estimation errors were calculated using the bootstrap method (Rubin 1981). The mean value and 90% confidence interval of the dip angle and Mw that minimized the misfit function were 13. 2° and [11.5°, 15.3°], and 8.85 and [8.75, 8.92], respectively. In this study, 689 stations were randomly selected from 689 Hi-net tiltmeter stations, and 70 stations were selected from the 70 F-net broadband seismometer stations, both allowing duplicates. The misfit function calculations were repeated 500 times. Figure 4 shows the resultant stacked traces for the improved parameters (dip angle, 13°; Mw, 8.85). The synthetic and observed waveforms were in quantitative agreement with the uncertainty of the background noise.

Discussion: Dip angle determination from pre-P gravity signal
As mentioned in "Introduction" section, long-period seismic waves have poor sensitivity to the seismic moment M 0 and dip angle δ , which are not well constrained for large shallow dip-slip events (Kanamori and Given 1981). This behavior is a result of all the excitation terms for the normal modes of the Earth becoming proportional to M 0 sin2δ to satisfy the traction-free boundary condition at the Earth's surface when the source depth becomes zero; therefore, M 0 and δ cannot be determined independently (Tsai et al. 2011).
Conversely, our analysis clearly showed that the horizontal and vertical components of the pre-P gravity signal have different dependencies on M 0 and δ ; therefore, the dip angle is successfully determined from the pre-P gravity signal. However, the excitation characteristics of the pre-P gravity signal remain difficult to understand. Because the excitation terms are common among elastic deformation and gravity change, it is questioned how the trade-off between M 0 and δ is resolved in the observation of the pre-P gravity signal. Further research is necessary to clarify this point.
Compared with the dip angle of the CMT solution, the estimated dip angle of 13.2° was a relatively large value, which may be related to the temporal properties of the pre-P gravity signal, which rapidly increases in amplitude with time, and the minute signal is readily masked by the P-wave arrival. Accordingly, the signal does not indicate the complete average of the moment release, but rather emphasizes the rupture in the region closer to the observation station. One key point is that the pre-P signal at each station is masked by P-wave arrival. If a large slip occurs in a region farther than the hypocenter from the observation station, the amplitude of the signal originating from the large slip will be overwhelmed by the initial P-wave that reaches the observation station from the hypocenter. Therefore, its contribution to the overall signal amplitude is limited. In the case of the 2011 Tohoku-Oki event, observation of the pre-P gravity signal on the Japanese land is thought to be heavily influenced by the rupture of the western, deeper, and large-dip regions of the fault, while that of the eastern, shallow, and small-dip regions has a limited effect.
The estimated moment magnitude was slightly smaller than the seismologically estimated values (9.0-9.1). Given that in most Japanese data the fault rupture was incomplete at the P-wave arrival time, this may mean that the total released moment of the earthquake was larger than expected from the beginning phase of the rupture to which the pre-P gravity signal has greater sensitivity.
The estimated dip angle and magnitude, different from the seismic ones, indicates that the pre-P signal provides a new observation window for the earthquake source, responding to spatiotemporal mass redistribution. This contrasts with conventional seismic observations, which measure the surface vibrations of the Earth.
Pre-P gravity signals have two distinct characteristics: almost simultaneous signal arrivals and different signal amplitudes from discrete slipping regions on a fault. The former comes from a virtually instantaneous propagation speed, and the latter comes from the Newton-gravity interaction with the inverse square law. These characteristics will help clarify the spatiotemporal characteristics of the source process, including fault finiteness, in future inversion analyses.

Conclusions
Using the synthetic waveforms calculated for a realistic Earth structure model as a reference, we detected the horizontal components of the pre-P gravity signal from the 2011 Tohoku-Oki earthquake in tiltmeter array data with a signal-to-noise ratio of 9.5. The discrepancy between the signal and synthetic amplitudes was resolved by changing the dip angle and Mw of the source, and the waveform inversion of the three-component records constrained the two source parameters in the ranges of 11.5°-15.3° and 8.75-8.92, respectively. The analysis of the observed pre-P gravity signal incorporated with the theoretical synthetics showed that the pre-P gravity signal opened a new observation window for the earthquake source, independent of the seismic observations.
Additional file 1: Fig. S1. Average noise spectrum of the Hi-net tiltmeter records. Fig. S2. Magnitude maps of the three-component pre-P gravity signal. Fig. S3. Stacked waveforms of the horizontal NS component data and synthetics.