Simulation of Broadband Strong Motion Based on the Empirical Green’s Spatial Derivative Method

Abstract

This study sought to simulate strong broadband seismic motions beyond the corner frequencies used, for the same events, in previous studies. To correct discrepancies among the corner frequencies of events, a scaling law based on the ω−2 model was assumed and the spectral amplitude decay beyond the corner frequency was compensated. The observations were also corrected for location, focal mechanism, and time of occurrence. After estimating the empirical Green’s tensor spatial derivative (EGTD) from 11 aftershock events, using 0.2–10 Hz band-pass-filtered waveforms, the strong motion records for the mainshock and aftershocks were simulated. In the simulation of each event, the EGTD elements were multiplied by the moment tensor elements followed by summation and then corrected in their spectral amplitude taking the corner frequency of each event into account. At the closest epicentral distance, and for most events, an acceptable agreement was found between calculated and observed waveforms. The results were also compared with the outcomes of simulations using the empirical Green’s function method. As EGTD elements are determined by local underground structures, they could prove useful for future structural studies.

1. Introduction

Prediction of the strong ground movements produced by large earthquakes has been demonstrated using the empirical Green’s function (EGF) method, proposed by Hartzell [1] and extended by Irikura [2]. This practical technique is most effective when the focal mechanism of a small event is similar to that of a targeted event. When the focal mechanisms differ more significantly, the empirical Green’s tensor spatial derivative (EGTD) method, proposed by Plicka and Zahradnik [3], is more appropriate. This can predict the ground motion for events with diverse focal mechanisms. Using single-station inversion of waveform data of several small events whose focal mechanisms and source time functions are well determined, the EGTD elements can be estimated, with the expectation of stable and accurate ground motion prediction. However, discussion in the literature of EGTD applications has been limited [4, 5, 6, 7, 8, 9, 10, 11].

In recent studies [9, 10], the authors sought to simulate strong broadband motions, beyond the corner frequencies, for the same events as in previous studies [6, 8]. This chapter extends these recent studies and applies the results to other stations. The scaling law based on the ω⁻² model [12] was used to correct differences between the corner frequencies of events, assuming spectral amplitude decay at higher frequencies. The 0.2–10 Hz band-pass-filtered waveforms of 11 aftershock events were used to estimate the EGTD and then to simulate the strong motion records for main and aftershocks. The EGTD elements were multiplied by the moment tensor elements followed by summation. The spectral amplitude was adjusted by taking the corner frequency of each event into account. Agreement between the simulated waveforms, at the closest epicentral distance, and observation were acceptable.

2. Locations of the epicenters and stations

Figure 1 shows epicenter locations, determined by the Japan Meteorological Agency (JMA), of the mainshock (M_J5.4) and 25 aftershocks (M_J3.1–4.7) of the 2001 Hyogo-ken Hokubu earthquake at the target station, HYG004, one of the K-NET stations operated by the National Research Institute for Earth Science and Disaster Prevention (NIED). The fault zone is roughly 4 km east-west and 6 km north-south. HYG004 is on a nearby rock site (6–10 km from the fault zone). Acceleration data for the mainshock and 15 of the aftershocks (M_J3.5–4.7) came from the K-NET. Two other stations, HYG001 and HYG007, are shown in the inset of Figure 1; they were also tested to verify the applicability of the EGTD method.

Figure 1.

Map showing the epicenter locations of the mainshock and aftershocks of the 2001 Hyogo-ken Hokubu earthquake. The mainshock and 15 aftershocks analyzed are indicated by a solid circle and solid triangles, respectively. The 10 aftershocks not used in the present study are represented by open diamonds. The target station HYG004 is shown by an inverted solid triangle. The inset at the top right shows the surrounding region with the target area (solid rectangle) and the other stations, HYG001 and HYG007. (After Ohori and Hisada [8] with slight modification).

3. Source model

The EGTD method is most accurate if the focal mechanism and source time functions are known as accurately as possible. The source model described below, and illustrated in Figure 2, has been reevaluated from previous work [6]. Strike, dip, and rake of a double-couple point source, and source depth, were estimated using the grid search technique [13]. The observed acceleration records at HYG004 for the mainshock and 15 aftershocks were integrated into velocity waveform data with a band-pass filter of 0.2–1.0 Hz. A total of 5 s of data, which included the P-wave arrival and S-wave main portions, was inverted. The theoretical Green’s function for the layered underground structure [14] was calculated assuming a smoothed ramp function with a rise time of 0.32 s. The searching ranges of the strike, dip, and rake angles were set within 20° of the solutions determined by the F-net of the NIED. The source depth was estimated as 8–12 km. Seismic moments, released by five sequential slips with intervals of 0.16 s, and a source time function were then determined by the least-squares method with nonnegative constraints [15]. Events 3, 17, 19, and 26 were excluded from the EGTD inversion because of a relatively large discrepancy in waveform matching between the observation data and synthesis.

Figure 2.

Map showing focal mechanisms and source time functions of 16 events (the mainshock and 15 aftershocks) from the work of Ohori and Hisada (2006). Single-station inversion for each event was performed using strong motion records recorded at the K-NET station HYG004. (After Ohori and Hisada [8]).

4. Estimation of the EGTD

The EGTD estimation method has been fully explained by Ohori and Hisada [6, 8]. It is applicable to simulation of strong motion in a frequency range below the corner frequency. The process is summarized below. Methods for compensating the spectral amplitude decay, beyond the corner frequencies, and simulating the broadband ground motion follow.

4.1. Basic equations

Ground motion displacement u_i(x_o, t) (i = x,y,z) excited by a double-couple point source is theoretically expressed as the convolution of moment tensor elements M_pq(x_s,τ) (p,q = x,y,z) and Green’s tensor spatial derivative elements G_ip,q (x_o, t | x_s,τ):

$$u_i\left(x_o,t\right)=M_{\mathit{pq}}{\left(x_s,\tau \right)}^{\ast }{G}_{\mathit{ip},\mathit{q}}\left(x_o,t|x_s,\mathrm{\tau }\right)$$E1

Hereafter, u_i(x_o, t), M_pq(x_s,τ) and G_ip,q (x_o, t | x_s,τ) are abbreviated as u_i, M_pq, and G_ip,q. Explicit expressions of M_pq for a double-couple point source are in the literature (e.g., see works by Aki and Richards [16]).

Considering symmetrical conditions (M_pq = M_qp) and no volume change [M_xx = −(M_yy + M_zz)] of the moment tensor elements, Eq. (1) can be rewritten as

$$u_i=\sum _{j=1}^5{M}_j\ast G_{\mathit{ij}}$$E2

where M_j (j = 1,2,…,5) is defined by M₁ = M_xy, M₂ = M_yy, M₃ = M_yz, M₄ = M_xz, and M₅ = M_zz, and G_ij (j = 1,2,…,5) is defined by G_i1 = G_ix,y + G_iy,x, G_i2 = G_iy,y - G_x,x, G_i3 = G_iy,z + G_iz,y, G_i4 = G_ix,z + G_iz,x, and G_i5 = G_iz,z - G_ix,x.

In the moment tensor inversion undertaken for a particular event using data of all possible components at all possible stations simultaneously, u_i and G_ij are given and M_j are the unknowns to be solved in a least-squares sense. Conversely, in the EGTD inversion computed for each component at each station using data from several events simultaneously, u_i and M_j are given and G_ij are the unknowns. While the moment tensor elements are determined from source parameters, the Green’s tensor spatial derivative elements express the underground structure of the area surrounding the source and the station.

4.2. Correction of the focal mechanisms

To compensate for differences in the locations of the main and aftershocks, the focal mechanisms are horizontally and vertically rotated, as described in the literature [4, 5], such that each event can be treated as a point source at the same location. Through the horizontal rotation, shown in Figure 3(a), the station azimuths of the mainshock and aftershocks are set to 90°, as measured from north [6]. Thus, the number of Green’s tensor spatial derivative elements is reduced to three (G_i1 = G_i4 = 0) for the radial component (i = y) and the vertical component (i = z) and two (G_i2 = G_i3 = G_i5 = 0) for the transverse component (i = x). Through the subsequent vertical rotation, as shown in Figure 3(b), the discrepancies in the take-off angles between the mainshock and aftershocks are corrected. The moment tensors are evaluated after the horizontal and vertical rotations.

Figure 3.

Schematic diagram explaining how the focal mechanisms are rotated to reduce the number of unknowns in the empirical Green’s tensor spatial derivative (EGTD) and to compensate for the different locations of the mainshock and aftershocks. In the horizontal rotation (a), the station azimuths for the mainshock and aftershocks are set to 90°, measured from north, so that the number of EGTD elements is reduced to three for the radial and vertical components and two for the transverse component. In the vertical rotation (b), following horizontal rotation, the discrepancies in the take-off angles between the mainshock and aftershocks are corrected. In the top right panel, take-off angles for the mainshock and aftershocks are denoted ξ_m and ξ_a, respectively. The difference between take-off angles, ξ_m-ξ_a, is slightly exaggerated in the bottom panel for clarity. (After Ohori and Hisada [8]).

4.3. Correction applied to the waveform data

Time shifts were estimated from the 0.2–1.0 Hz band-pass-filtered velocity waveforms and then applied to the observation data of aftershocks to fit their S-wave arrival time with that of the mainshock. In addition, the observation data was deconvolved to remove discrepancy in the source time function. The observed waveforms, used in the estimation of the EGTDs, were corrected such that the source time function has a constant seismic moment (1.0 × 10¹⁵ Nm, nearly equal to M_w4.0) and a single-isosceles slip velocity function with a rise time of 0.32 s.

Assuming the source spectrum obeys the ω⁻² model [12], the corner frequency of the mainshock was about 1.0 Hz, while those of 11 aftershocks were between 1.2 and 3.0 Hz, as shown in Figure 4. To simulate the broadband ground motion beyond the corner frequency, the effect of the differences among corner frequencies must be removed. Again, assuming a ω⁻² scaling law, and compensating for spectral amplitude decay beyond the corner frequency, each event could then be treated as having the same corner frequency as that of the mainshock.

Figure 4.

Relationship between the corner frequency and the seismic moment. The numeral next to solid circle represents the sequential event number in the present study. According to the scaling law, the seismic moment is inversely proportional to the cube of the corner frequency when the stress drop is constant.

4.4. EGTD estimation

Using the focal mechanisms of aftershocks, rotated and time corrected, as described above, simultaneous linear equations for each component were solved for each of the sampling data sets. No smoothing or minimization for unknown parameters was included in the EGTD estimation. Figure 5 illustrates the transverse component elements of an EGTD for HYG004. Each element is scaled, the same as the mainshock, for an event with a seismic moment of 1.0 × 10¹⁵ Nm and a corner frequency of 1.0 Hz. As the Green’s tensor spatial derivative elements are determined by the local underground structure, the EGTD elements could be useful for future structural studies.

Figure 5.

An example of the EGTD transverse component elements with normalized peak amplitude.

5. Simulation of the strong ground motion using the EGTD

Figure 6 compares the observed radial, transverse, and vertical component velocity waveforms for HYG004, normalized for time, seismic moment, and corner frequency, with a 0.2–10 Hz band-pass filtering and corresponding synthesis calculated from the EGTD. The top trace for the mainshock (Event 1) is not included in the EGTD inversion. Considering the complexity of the high-frequency components, the broadband synthesis, using EGTD, acceptably reproduced the observed waveforms.

Figure 6.

Comparison of 0.2–10 Hz band-pass-filtered observed velocity waveforms for HYG004 (red lines) and corresponding syntheses calculated using the EGTD (black lines). For each trace, the source time function, seismic moment, and corner frequency of each event are taken into account. Absolute peak amplitudes are given at the end of each trace: Upper numerals are observed amplitudes and lower are synthesized. Note that the top trace for the mainshock (Event 1) is not included in the EGTD inversion. (a) Radial-comp; (b) Transverse-comp; (c) Vertical-comp.

Three frequency ranges of the band-pass filter are compared in Figure 7: 0.2–1.0, 1.0–10, and 0.2–10 Hz. From this figure, it is evident that simulation results for HYG004 from EGTD are most accurate in a frequency range of 0.2–1.0 Hz, except that the radial components of Event 14 are somewhat overestimated. The results for a frequency range of 1.0–10 Hz are acceptable and very similar to those for a frequency range of 0.2–10 Hz. For most of the events, the maximum amplitude ratio between the synthesized velocity waveforms and observed data is between 0.5 and 1.5. On the whole, the simulation of broadband ground motion, using the EGTD method, successfully reproduced the observed waveforms.

Figure 7.

Comparison of the maximum amplitude ratio between the synthesis and observed data for HYG004. The frequency ranges of the band-pass filters are 0.2–1.0 Hz (upper), 1.0–10 Hz (middle), and 0.2–10 Hz (lower).

Finally, Figure 8 shows the results for the other two stations. In this figure, because of data quality issues, only transverse components are shown. Generally, the waveform matching between synthesis and observatory data is inferior to that for HYG004 (see Figure 6(b)). For HYG002, the synthesis underestimated motion for Events 2, 5, 6, and 25. For HYG007, synthesis also underestimated motion for Event 6. On the other hand, for the mainshock (Event 1), the agreement between synthesized and observed waveforms was satisfactory.

Figure 8.

Comparison of 0.2–10 Hz band-pass-filtered observed velocity waveforms (red lines) and corresponding syntheses calculated from the EGTD (black lines) for transverse components of data from HYG001 (a) and HYG007 (b). Other conditions are the same as in Figure 6.

6. Comparison of the simulation results obtained using the EGTD and EGF methods

A previous study [8] examined the accuracy of simulation results using the EGTD method by comparing them with those obtained using the alternative empirical approach—the EGF method [1, 2]. In that study, velocity waveform data with a 0.2–1.0 Hz band-pass filtering was used. This study compared the two methods over a broader frequency range. The EGF method used was almost the same as in the previous study [8], except the differences in the corner frequencies between the mainshock and each aftershock were corrected using the scaling law of the ω⁻² model [12]. In the EGF approach, each aftershock was used in the simulation of the mainshock. To compensate for differences in radiation pattern coefficients, for each aftershock, the data were multiplied by the ratio of radiation pattern coefficients between the mainshock and the aftershock. The radiation pattern coefficients of the SV-waves were used for simulation of the radial components, while SH-waves were used for the transverse components. The source time function of the mainshock was then convolved.

Figure 9 compares the synthesized waveforms of the mainshock for HYG004 (Event 1) obtained using the EGF method, with the observation data and the synthesis obtained using the EGTD method. The EGF method is applied for each of the 11 aftershocks. The amplitude levels of the synthesized waveforms, obtained using the two methods, were mostly in agreement with the observed waveform over the full duration of the analysis. In the low-frequency range examined in the previous study [8], the EGTD method provided a better match with the data than the EGF method. In the broadband frequency range, in this present study, the EGTD method seems to give stable results but not the best match with the data. It is noteworthy that the EGF method, in the case of Events 9 and 16, acceptably reproduced the amplitude for the radial components but twice, and almost a third time, overestimated the amplitude for the transverse components. For simulation using the EGF methods in a low-frequency range, a water-level of 0.20 has best suppressed extraordinary overestimation [8]. A water-level of 0.20, as the minimum absolute value of radiation coefficient, was tested with the current data but could not suppress this overestimation. This may arise from the complexity of high-frequency source processes. The synthesized waveforms using the EGF method look different from event to event. It therefore seems that the EGTD method, by removing the dependency of the EGF method on aftershock selection, provides more stable results.

Figure 9.

Comparisons of the synthesized waveforms of the mainshock for HYG004 (Event 1) obtained using the EGF method and each aftershock with the observed data and the synthesis obtained using the EGTD method. The numeral after ‘EGF’ represents the sequential event number of the specific aftershock used in the simulation for the mainshock. (a) Radial-comp; (b) Transverse-comp.

7. Conclusions

The applicability of the EGTD method to simulating near-field strong motion seismic records has been demonstrated. Previous studies [6, 8] in EGTD estimation used the low-corner frequency of 1.0 Hz for the mainshock of the 2001 Hyogo-ken Hokubu earthquake (M_J 5.4) and targeted 0.2–1.0 Hz band-pass-filtered velocity waveforms. This study extended the upper limit of the target frequency range to 10 Hz, while the corner frequency of the events was in a range from 1.0 to 3.0 Hz. To correct corner frequency discrepancies, the scaling law based on the ω⁻² model [12] was used to compensate the spectral amplitude decay beyond the corner frequency. Agreement between the observed and calculated waveforms for the mainshock was satisfactory for a long duration, and there was a good match between amplitudes. Further data accumulation and investigation enhance the applicability of the EGTD method.

Acknowledgments

The strong motion data used in this study were recorded at the K-NET stations, provided by the National Research Institute for Earth Science and Disaster Prevention (NIED) on their website (http://www.kyoshin.bosai.go.jp/kyoshin/, last accessed February 2018). The Japan Meteorological Agency (JMA) unified hypocenter catalogue, and the F-net source parameters were also provided by the NIED on their website (http://www.fnet.bosai.go.jp/top.php, last accessed Feburary 2018). This study was partially supported by Grants in Aid for Scientific Research (C) (16 K01316).