This effect is observed in the velocity time histories of many strong-motion earthquakes, such as for the 2015 Nepal earthquake in Figure 3A (Parajuli and Kiyono, 2015), and the near-source records of Mian Zu Qing Ping (MZQ, at a distance of about 1.7 km from the fault plane), during the 2008 Wenchuan earthquake, China (Figure 3B), which exhibits a two-sided velocity pulse (Lu et al., 2010). The intensive impulsive velocity effect is induced from the rupture directivity process and occurs in the direction vertical to the fault plane.

The fling effect is induced from the permanent tectonic offset of a rupturing fault. For example, Figure 4 shows the typical step-like displacements records of Si Fang Ba Jiao (SFB, at a distance of 1.2 km from the fault plane) during the 2008 Wenchuan earthquake and the TCU052 (at 1.7 km to the fault plane) and TCU072 (about 14 km to the fault plane) during the 1999 Chi-Chi, Taiwan, earthquake (Kramer, 1996). These records can be explained by elastic rebound theory. When the shear stress reaches the shear strength of the rock along the fault, the rock fails, and the accumulated strain energy is released. As this rupture progresses, it causes ground motions and induces the permanent displacement of the surface (Lu et al., 2008).

The third significant feature is the hanging-wall effect, whereby ground motions at sites located on the hanging wall of a dip-slip fault are larger than at sites located on the footwall at the same distance. Hanging-wall effects have been observed in the records of the 1994 Northridge earthquake and the 1999 Chi-Chi, Taiwan. The main reason for this effect is the wave propagation distance and the multi-reflection and refraction of the propagating waves between the surface and the fault planes (Liu et al., 2006).

First, based on the representation theorem, the displacements given by a kinematic fault model can be expressed in the frequency domain as: (1)Uk(Y,ω)=∫STik(X,Y;ω)Di(X;ω)dS where U_k is the kth component of displacement in Cartesian coordinates at an observation point Y, X is the source point on the fault plane, ω is the circular frequency, S is the fault plane, T_ik is the traction Green’s function, and D_i is the ith component of the fault slip. The slip rupture model and parameters are shown in Figure 5.

To simulate theoretical strong ground motions for near-faults, Hisada and Bielak (2003) introduced an efficient method for the fault integration of the representation theorem. This method evaluates the fault integration of the dynamic and static terms separately as: (2)Uk(Y,ω)=∫S{TikD(X,Y;ω)−TikS(X,Y)}Di(X;ω)dS+∫STikS(X,Y)Di(X;ω)dS where TikD and TikS are the dynamic and static traction Green’s functions of the layered half-space. On the right-hand side of Eq. 2, the first integral contains the dynamic terms, and the second integral includes the static terms (which means ω = 0).

The stochastic Green’s function method is very popular, flexible, applicable to structural designs, and can also generate dynamic ground displacements in the vicinity of the fault. However, it is difficult to simulate fault movements such as a “static” fling-up (permanent displacement) of the fault. The theoretical Green’s function method can generate both dynamic and static ground motions and has the advantage that the displacement near the fault almost directly reflects the source characteristics. However, its application is limited to idealized layered media and time-consuming. Thus, to incorporate the displacements induced by active faults, the hybrid method of modified stochastic (for dynamic calculations) and theoretical (for static calculations) Green’s functions is proposed for synthesizing the near-fault ground displacement. We shall show that this approach is much faster than some other simulation methods.

As the observation point is close to the fault plane, the waveforms exhibit sharp peaks of short duration in the region close to the observation point (Hisada and Bielak, 2003) Even if the distance to the fault is very small, the dynamic ground motion can be calculated by superposing small element waveforms. Thus, the stochastic Green’s function method (Irikura, 1983; Irikura and Miyake, 2006) is adopted to calculate the dynamic terms, as shown in Eq. 2. The basic principle of the original statistical Green’s function method is as follows: a large earthquake is composed of a series of small earthquakes and statistically calculated small earthquakes (namely the statistical Green’s function) are properly selected as the ground response caused by small-area sources. That is, statistical Green’s functions are overlaid in a specified manner to obtain the time history of a strong earthquake. Equations 3 and 4 describe the main procedure of this method: (3)U(t)=∑m=1NL∑n=1NWrrmn× [u(t−tmn)+∑k=1(ND−1)n′1n′ u (t−tmn−(k−1)·τ(ND−1)n′)] (4)tmn=rmn−roVS+ξmnVR where U(t) is the synthetic main-shock ground motion displacement, u(t) is the observed small ground motion, NL, NW, and ND are the ratios of fault length, fault width, and slip values between large and small events, and t_mn is the delay time of the point source (m, n) on the rupture surface. The parameter τ denotes the rise time of a small earthquake, and Vs, V_R denote the S-wave velocity near the earthquake source and rupture velocity, respectively. ξ_mn is the distance from (m, n) located on the fault plane to the starting point, as shown in Figure 6, and n' is an integer that weakens the artificial periodicity of n so that the tick interval represents the sampling rate. The other notation is defined in the schematic diagram of the superposition of small events in Figure 6.

As the conventional stochastic Green’s function method is mainly adopted to generate far-fault ground motions, it only considers the far-field terms of the S-wave. For regions very close to the source, the ground motion features are rather complicated and the attenuation relationship is very different from that of the far-fault ground motions. Thus, the complete waveforms, including the near-, intermediate-, and far-field items of P-waves and S-waves, should be taken into consideration. The complete waveforms have been calculated using a finite-difference method (Dreger et al., 2007), whereas Onishi and Horike (2004) improved the stochastic Green’s function by introducing theoretical radiation coefficients to the P-, SV-, and SH-waves and using a ray tracing technique in the layered half-spaces. In this paper, we introduce a simple model in which the complete waveforms are calculated based on the method proposed by Nozu (2006), which used the ratio of the Fourier transform of the total wave to introduce near-field and intermediate-field terms of P- and S-waves with respect to that of far-field S-wave. This process is described by Eq. 5 [the derivation and more details can be found in Nozu (2006)]. Furthermore, the calculation for far-field P-wave is the same, except that the P-wave velocity is used instead of the S-wave velocity in the radiation pattern.

The stochastic Green’s function method does not consider the static displacement, because the statistically calculated small earthquake does not contain permanent displacements. To obtain the near-fault time history, the static terms (second integral on the right-hand side of Eq. 2) describe the attenuation of the slip function due to the static traction of the Green’s function, and the theoretical Green’s function is calculated by the wavenumber integration method (Hisada, 1994, 1995; Hisada and Bielak, 2003). In this study, we only use the second term of Eq. 2, and calculate the static displacement as: (6)Uk(Y;ω)=∫STikS(X,Y)Di(X;ω)dS where Uk(Y; ω) is the static displacement (in this situation, ω = 0), and D_i is the same as in Eq. 1. Clearly, it is very easy to calculate Eq. 6; the details can be seen in Figure 5.

First, a simple and pure strike-slip model of surface faulting (see Figure 7) was simulated by taking into account both the static and dynamic terms of the near-field ground motions. In Figure 7A, the fault size is 10.0 km in length and 5.0 km in width, with strike angles of N0°E, a dip angle of 90°, and rake angle of 0°, and the seismic moments were set to M_o = 3.825 × 10²⁵ dyne cm. Surface observation points (numbered 1–12) were calculated along a line perpendicular to the fault plane. Note that points 1 and 2 are only 100 m and 500 m away from the fault trace, respectively. Similarly, points 3–12 are 1.5~10.5 km away from the fault trace at intervals of 1.0 km. The characteristic fault model is shown in Figure 7B, which was divided into 10 × 5 sub-faults, the asperity has a slip of 1.0 m, stress drop of 141.0 bar and slip of 0.6285 m, stress drop of 28.8 bar for the background. Another important parameter is the slip velocity function; we assume it has the form of an isosceles triangle with 1.0-s duration here. The homogeneous half-space with physical properties is also described in Figure 7A. The static terms are calculated by the proposed representation theorem, Eq. 6, and the dynamic terms are simulated by a modified statistical Green’s function.

The symmetry of the model indicates that this is a pure strike-slip fault. The fault-normal components (namely along the EW direction) are very small, and the vertical components are negligible compared with the fault-parallel components. Thus, only the calculated velocity and displacement time histories along the fault-normal direction are shown in Figures 8A–F. As mentioned above, the near-fault ground motion is rather complicated.

For dynamic terms, in Figure 8A, the velocity time histories for observation points 1~4 (at distances of 0.1~3.5 km from the fault traction) exhibit random dynamic peak values, whereas at sites 4~12, the peak dynamic terms exhibit a degree of attenuation. Somerville (1998) obtained an empirical relationship for the variation of the PGV (peak ground velocity) as a function of the moment magnitude and the distance to the causative fault, assuming that PGV varies as 1/R for distances larger than 3.0 km. The simulated velocity time histories exhibit this tendency. The displacement time histories of dynamic terms in Figure 8D also illustrate that the peak displacements of dynamic terms not only exhibit pure time delay but also have a random distribution located very near the causative fault.

For the static terms, Figures 8B,E show that the peak values of velocity and displacement attenuate with distance from the surface fault traction. Compared with the dynamic terms, this decrease is more pronounced. In this situation, almost all the static terms are larger than the dynamic displacements in Figure 8E, which confirms that the permanent displacements should not be neglected when designing near-fault bridges or pipelines.

In the results given by the theoretical Green’s function method, the attenuation is different: very near the fault traction, the static terms are dominant, whereas further away from the fault traction, the dynamic terms prevail. As different analysis methods have been applied here, these differences are acceptable because the hybrid method considers the complete waveforms of the dynamic terms; according to Nozu (2006), the near-fault ground motion is heavily dependent on the frequency and velocity. This illustrates that the near-fault ground motion is rather complicated.

Next, we calculated the near-fault strong motions using a dip-slip surface faulting model. The fault model is shown in Figure 9, with strike angles of N0°E, a dip angle of 45°, and rake angle of 90°(reverse fault), and the seismic moment was M_o = 3.825 × 10²⁵ dyne cm. Here, 24 observation points were positioned on the free surface along the line perpendicular to the fault traction; points 1~12 were located on the footwall, and points 13~24 were on the hanging-wall side. Points 12 and 13 were closest to the fault trace, 100 m on either side, with points 11 and 14 of 500 m away from the fault trace. Points 1~10 and 15~24 were arranged at intervals of 1.0 km. The characteristic fault model with slip distribution and the stress drop was schematically the same as example 1 shown in Figure 7B, and slip function was also assumed as triangle with duration of 1.0 s. The material properties of the layered half-space are also shown in Figure 9.

As this reverse-fault model considers a pure dip-slip fault, the simulated results along the fault-normal (EW) and up–down directions are analyzed. Figures 10A–F show the results along the fault-normal direction, and Figures 11A–F showed the simulation time histories along the up–down direction, including the dynamic terms, static terms, and total velocity and displacement at the 24 observation points.

For the EW direction, the velocity time histories in Figure 10A indicated that there are no major differences between the sites located on the hanging wall or the footwall side; the dynamic terms are all larger than the static velocity as compared with Figures 10A,B, whereas the static displacements (Figure 10E) are all larger than the dynamic displacements in Figure 10D. Thus, the permanent displacements should be considered when conducting a seismic analysis of near-fault spatially expanded structures.

Along the up–down direction in Figure 11, the results illustrate the hang-wall effect, i.e., the observation points on the hang-wall indicate larger motions than the observation sites on the footwall, especially the static terms of both velocity in Figure 11A and displacement in Figure 11D.

As the distance from the fault trace increases, the static displacements attenuate rapidly, and the dynamic components become dominant, the same attenuate tendency as in strike-slip fault model.

The proposed hybrid method was used to simulate the time histories of near-fault ground motions on the 1999 Chi-Chi earthquake in Taiwan, including the permanent displacements.

The 1999 Chi-Chi earthquake (Mw 7.6; September 20, 1999, 14:47:15.9 UTC; hypocenter located at 23.853°N, 120.816°E at a depth of 7.5 km) inflicted severe regional-scale damage on Taiwan. The surface fracture trace ran along the Chelunpu fault, with strike angles of N5°E and a dip angle of 30°. According to Wu’s (Wu et al., 2001) fault model A and Hisada’s simulation (Hisada, 0000), the fracture had a length and width of 82 and 42 km, respectively. The characteristic fault model with two asperities is plotted in Figure 12B. The fault plane was divided into 21 × 11 sub-faults, and the seismic moments were set to M_oa1 = 8.2315 × 10²⁵ for asperity 1, M_oa2 = 6.3874 × 10²⁵ for asperity 2, and M_ob = 5.2829 × 10²² for the background of each sub-fault. In addition, the 16-layer soil conditions are listed in Table 1.

We used a rupture velocity of 3.0 km/s, rake angle of 60°, and calculated the ground motions at four observation sites on the hanging wall which shown in Figure 12A, TCU052 (1.7 km away from the surface fault traction), TCU072 (14.2 km away from the surface fault traction), TCU089 (16.4 km away from the surface fault traction), and TCU078 (16.2 km away from the surface fault traction). For the slip velocity, we considered five time windows, based on Nakamura and Miyatake (1997). Note that band-pass filters (0.0~0.01 Hz and 10.0~15.0 Hz) and baseline correction were applied in processing the results.

For reasons of limited space, we only present the results from two stations, TCU052 and TCU072, including the velocity in Figure 13 and displacement in Figure 14 time histories and comparisons with the observed recordings.

For the velocity time histories in Figure 13, the simulated velocities along three directions are smaller than the observed waves at TCT052 and generally fit well with the observed velocities at TCU072. As the near-fault ground motions are highly dependent on the frequency in the proposed method, it is clear that high-frequency regions affect the simulated results.

In the displacement time histories shown in Figure 14, the vibrations of the simulated results are not so obvious in the waveforms, which suggest that the dynamic parts are smaller than the static terms. The simulated displacement along the NS direction for both TCU052 and TCU072 are notably different to the observed recordings. The displacement along the UD direction at TCU052 in Figure 14A is a little smaller than the observed movement. The other simulated displacements are in good agreement with the observations, which showed that the proposed hybrid method could effectively simulate the near-fault strong ground motions.