Open AccessArticle

Stochastic Green’s Function Method Considering Non-Uniform Rise Time Distribution to Simulate 3D Broadband Ground Motion

Longfei Ji

¹,

Xu Xie

^1,*

and

Xiaoyu Pan

College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, China

Department of Roads and Bridges, Zhejiang Institute of Communications, Hangzhou 311112, China

Author to whom correspondence should be addressed.

Appl. Sci. 2024, 14(21), 9796; https://doi.org/10.3390/app14219796

Submission received: 5 August 2024 / Revised: 19 October 2024 / Accepted: 24 October 2024 / Published: 26 October 2024

Download

Browse Figures

Versions Notes

Abstract

The stochastic Green’s function method has been widely used in the field of ground motion simulation in recent years. It is generally assumed that the rise time of each subfault is the same in this method. Since the rise time significantly influences the amplitude of simulation results in the intermediate frequency band, to improve the accuracy of stochastic Green’s function method for near-fault broadband ground motion simulation, referring to the numerical simulation results of Day, the rise time is assumed to be non-uniformly distributed on the fault, and an improved approximate expression of rise time on a rectangular fault considering that the rupture starting point may be at any position and the aspect ratio may be arbitrary is proposed. Additionally, the contributions of P, SV and SH wave are considered, respectively, and an improved stochastic Green’s function method is proposed for 3D broadband ground motion simulation. Taking the 1994 Northridge earthquake in America and 2013 Lushan earthquake in China as examples, under different subfault division numbers, the synthesized source spectra are compared with the omega-squared theoretical source spectra of the large earthquake, and the simulated ground motions at observation points are compared with observed records to verify the effectiveness of the improved method. The results show that when the Northridge earthquake fault and Lushan earthquake fault are divided into 9 × 10 subfaults and 11 × 7 subfaults, respectively, the simulation results obtained using the improved method are close to the observed records in the broadband frequency range. Therefore, the improved method can effectively simulate the 3D ground motion in near-fault regions.

Keywords:

3D ground motion simulation; stochastic Green’s function method; rise time; 1994 Northridge earthquake; 2013 Lushan earthquake

1. Introduction

1.1. Background

In recent years, the stochastic ground motion simulation method considering the effect of the rupture process of faults, the propagation path of seismic waves and the local site amplification has developed rapidly in the field of seismology. Several methods have been proposed, such as the stochastic Green’s function method [1,2] based on the scaling law between large and small earthquakes, the stochastic finite-fault method [3,4,5] that directly generates ground motions of subfaults and superimposes them to synthesize the ground motion of the large earthquake, and the stochastic summation of Green’s function method [6]. The basic idea of the stochastic method is to divide the fault into several subfaults. The stochastic point source method [7] is used to generate the ground motion time histories of small earthquakes in subfaults. After that, the ground motions of small earthquakes are superimposed according to the arrival time of the seismic waves. Then the ground motion time history of the large earthquake at the observation point caused by the rupture of the entire fault can be obtained.

In the stochastic Green’s function method, the rise time is an important parameter that affects the conversion function between small and large earthquakes. The commonly used formula for calculating the rise time τ at present is τ = α_rW/V_R, where α_r is an empirical coefficient between 0.25 and 0.6, generally taken as 0.5, W and V_R represent the width of the fault and the propagation velocity of rupture on the fault plane, respectively. Day [8] applied the 3D finite difference method to numerically simulate the earthquake slip function and analyzed the distribution of rise time on the rectangular fault with different aspect ratios (fault length/fault width) when the rupture starting point is located at the center of the fault. The results show that the rise time at a point on the fault plane is related to its position on the fault plane. Near the rupture starting point, the rise time is longer, while near the edge of the fault, the rise time is shorter. The commonly used formula mentioned above is approximately valid along the center line of a long, narrow fault; however, its applicability is limited if the aspect ratio of the fault changes and the fault does not belong to a long, narrow fault. In addition, the rise time at different points on the fault plane is different, and the actual earthquake slip cannot be reflected if the ground motion in each subfault is calculated using the same rise time. From the perspective of energy released during earthquakes, when the slip of the fault is constant, the longer the rise time, the slower the slip velocity of the fault, and the slower the energy is released during the earthquake. Therefore, it is necessary to establish a reasonable formula for calculating the rise time to reflect the actual earthquake slip and improve the accuracy of simulating and predicting the ground motion using the stochastic Green’s function method.

In this study, on the basis of improving the expression for the rise time, the contributions of the P wave, SV wave, and SH wave radiated from the source are considered, respectively [9,10], and an improved stochastic Green’s function method for 3D broadband ground motion simulation is proposed. Taking the 1994 Northridge earthquake in the United States and the 2013 Lushan earthquake in China as examples, under different subfault division numbers, the source spectra synthesized using the proposed approximate expression for the rise time are compared with the ω⁻² theoretical source spectra of the large earthquake, then the ground motions at the observation points are simulated using the improved method, and the simulation results are compared with the observed records, verifying the effectiveness of the improved stochastic Green’s function method proposed in this study.

1.2. Literary Review

In recent years, Otarola and Ruiz [11], Ruiz et al. [12], Wang et al. [13], Ojeda et al. [14], Li et al. [15], Qiang et al. [16], Luo [17], and Ji et al. [18] have conducted a series of studies on 3D ground motion simulations using stochastic methods. However, the simulation results cannot fully reflect the characteristics of near-fault ground motion, and the simulation accuracy needs to be improved.

Day [8] and Miyatake et al. [19] provided an approximate expression considering the non-uniform distribution of the rise time on the fault plane; however, this expression was derived under the condition that the fault belongs to the long, narrow fault and the rupture starting point is located at the center of the fault. On this basis, Hikita [20] proposed an approximate calculation formula for the rise time considering that the rupture starting point may be at any position on the fault; nevertheless, the application of this formula is still limited by the condition of the long, narrow fault.

In order to improve the accuracy of near-fault ground motion simulation using the stochastic Green’s function method, in this study, an improved approximate expression for the rise time on the rectangular fault considering that the rupture starting point may be at any position on the fault and the aspect ratio may be arbitrary is proposed referring to the numerical simulation results of Day [8].

2. Stochastic Green’s Function Method for 3D Ground Motion Simulation and Its Improvement

The principle of the stochastic Green’s function method is to divide a fault with non-uniform slip distribution into several subfaults, each of which slip sequentially in the rupture direction of the fault and radiate seismic energy to the surrounding area. The ground motion of the large earthquake at the observation point is the result of the superposition of the ground motions of the small earthquakes radiated from the subfaults [1,2], as shown in Figure 1. The following are the steps for simulating 3D ground motion in this study. The Fourier amplitude spectra of the P wave, SV wave, and SH wave of the small earthquake in each subfault are generated by the stochastic point source method. Then the phase spectra of the small earthquakes in subfaults are generated according to the similarity between the phase difference spectrum and the ground motion time history envelope curve. The acceleration time histories of three components

u_{ij}^{P} (t)

u_{ij}^{SV} (t)

and

u_{i j}^{SH} (t)

of the small earthquakes in the subfaults are obtained by inverse Fourier transform, and the P wave, SV wave, and SH wave of the large earthquake are constructed by superposing the small earthquakes in accordance with the type of wave. After that, the time histories of the three components of the large earthquake are projected into the NS, EW, and UD directions. Then the local site amplification effect of the horizontal components of the large earthquake is included, and the simulation results of the ground motion at the surface can be obtained.

2.1. Fourier Amplitude Spectra of Small Earthquakes in Subfaults

Based on the ω⁻² model, the acceleration Fourier spectrum

A_{i j}^{K} (f)

excluding the local site amplification effect of the small earthquake at the observation point caused by the ij subfault can be expressed as [7,9,10]:

A_{i j}^{K} (f) = \frac{F_{s}^{K} R_{θ ϕ i j}^{K}}{4 π ρ v^{3}} \frac{{(2 π f)}^{2} M_{0 i j}}{1 + {(\frac{f}{f_{c i j}^{K}})}^{2}} {[1 + {(\frac{f}{f_{\max i j}})}^{2}]}^{- \frac{1}{2}} {(\frac{ρ v}{ρ_{z} v_{z}})}^{\frac{1}{2}} N_{i j}^{K} (f) \frac{1}{r_{i j}} \exp [\frac{- π f r_{i j}}{Q_{K} (f) v}]

(1)

In which f is the frequency; the superscript K represents the P wave, SV wave, or SH wave; and the subscripts i, j indicate the serial number of the subfault along the length and the width.

F_{s}^{K}

is the amplification factor of the free surface. As seismic waves near the ground surface propagate in a nearly vertical direction,

F_{s}^{K}

is taken as 2.0 for all types of waves.

R_{θ ϕ i j}^{K}

is the radiation pattern considering the frequency and distance dependence [18]. M_0ij and r_ij represent the seismic moment and the source distance of the small earthquake, respectively.

f_{c i j}^{K}

and f_maxij are the corner frequency [21,22] and high-frequency cutoff frequency [23] of the small earthquake, respectively. ρ, v represent the medium density and the P-wave or S-wave velocity near the source, respectively. ρ_z, v_z represent the medium density and P-wave or S-wave velocity at the seismic bedrock, respectively.

N_{i j}^{K} (f)

is the correction term including the influence of the near-field and intermediate-field terms of the ij subfault [9,10,24]. Q_K(f) is the quality factor, which represents the anelastic attenuation characteristics of the propagation path.

In this study, in order to improve the accuracy of low-frequency ground motion simulation using the stochastic Green’s function method, according to the study of Pitarka et al. [25] and Kotha et al. [26] on the frequency and distance dependence of radiation pattern,

R_{θ ϕ i j}^{K}

in Equation (1) is calculated using the formula that considers the effects of both frequency and distance [18]:

R_{θ ϕ i j}^{K} (f, r_{i j}) = \{\begin{matrix} R_{θ ϕ f i j}^{K} (f) \\ \frac{[\log (r_{2}) - \log (r_{i j})] R_{θ ϕ f i j}^{K} (f) + [\log (r_{i j}) - \log (r_{1})] R_{θ ϕ m}^{K}}{\log (r_{2}) - \log (r_{1})} \\ R_{θ ϕ m}^{K} \end{matrix} \begin{matrix} , \\ , \\ , \end{matrix} \begin{matrix} r_{i j} \leq r_{1} \\ r_{1} < r_{i j} \leq r_{2} \\ r_{i j} > r_{2} \end{matrix}

(2)

where

R_{θ ϕ f i j}^{K} (f) = \{\begin{matrix} R_{θ ϕ 0 i j}^{K} \\ \frac{[\log (f_{2}) - \log (f)] R_{θ ϕ 0 i j}^{K} + [\log (f) - \log (f_{1})] R_{θ ϕ m}^{K}}{\log (f_{2}) - \log (f_{1})} \\ R_{θ ϕ m}^{K} \end{matrix} \begin{matrix} , \\ , \\ , \end{matrix} \begin{matrix} f \leq f_{1} \\ f_{1} < f \leq f_{2} \\ f > f_{2} \end{matrix}

(3)

In Equations (2) and (3), r_ij represents the source distance of the ij subfault, r₁ and r₂ represent the distance range where the radiation coefficient transitions linearly with the logarithm of the source distance, with r₁ = 40 km and r₂ = 100 km, respectively; f₁ and f₂ represent the frequency range in which the radiation coefficient transitions linearly with the logarithm of the frequency. In this study, it is assumed that f₁ = 0.5 Hz and f₂ = 5 Hz.

R_{θ ϕ m}^{K}

is the average radiation coefficient at high frequencies.

R_{θ ϕ m}^{P}

and

R_{θ ϕ m}^{S}

are taken as 0.52 and 0.63, respectively [27].

R_{θ ϕ 0 i j}^{K}

is the theoretical radiation coefficient at low frequencies, taking the relative position relationship between the source and the observation point and the slip direction into account. The expressions are [28]:

\begin{array}{l} R_{θ ϕ 0 i j}^{P} & = \cos λ \sin δ \sin^{2} i_{ξ i j} \sin 2 (ϕ_{i j} - ϕ_{s}) - \cos λ \cos δ \sin 2 i_{ξ i j} \cos (ϕ_{i j} - ϕ_{s}) \\ + \sin λ \sin 2 δ \times [\cos^{2} i_{ξ i j} - \sin^{2} i_{ξ i j} \sin^{2} (ϕ_{i j} - ϕ_{s})] + \sin λ \cos 2 δ \sin 2 i_{ξ i j} \sin (ϕ_{i j} - ϕ_{s}) \end{array}

(4)

\begin{array}{l} R_{θ ϕ 0 i j}^{SV} & = \sin λ \cos 2 δ \cos 2 i_{ξ i j} \sin (ϕ_{i j} - ϕ_{s}) - \cos λ \cos δ \cos 2 i_{ξ i j} \cos (ϕ_{i j} - ϕ_{s}) \\ + \frac{1}{2} \cos λ \sin δ \sin 2 i_{ξ i j} \sin 2 (ϕ_{i j} - ϕ_{s}) - \frac{1}{2} \sin λ \sin 2 δ \sin 2 i_{ξ i j} [1 + \sin^{2} (ϕ_{i j} - ϕ_{s})] \end{array}

(5)

\begin{array}{l} R_{θ ϕ 0 i j}^{SH} & = \cos λ \cos δ \cos i_{ξ i j} \sin (ϕ_{i j} - ϕ_{s}) + \cos λ \sin δ \sin i_{ξ i j} \cos 2 (ϕ_{i j} - ϕ_{s}) \\ + \sin λ \cos 2 δ \cos i_{ξ i j} \cos (ϕ_{i j} - ϕ_{s}) - \frac{1}{2} \sin λ \sin 2 δ \sin i_{ξ i j} \sin 2 (ϕ_{i j} - ϕ_{s}) \end{array}

(6)

In Equations (4)–(6), ϕ_s and δ represent the strike and the dip angle of the fault plane, respectively, λ is the rake angle, ϕ_ij is the azimuth angle from the center of the ij subfault to the observation point, and i_ξij is the takeoff angle, which is defined as the included angle between the direction of the seismic ray leaving the source and the vertical downward direction. Assuming that the propagation medium is homogeneous, the seismic wave travels along a straight line, and the takeoff angle refers to the angle between the line connecting the center of the ij subfault and the observation point and the vertical downward direction. In order to avoid underestimating the amplitude of the simulated ground motion, the lower limit of the absolute value of

R_{θ ϕ 0 i j}^{K}

is set to 0.1 [9,10,27].

2.2. Phase Spectra of Small Earthquakes in Subfaults

The random phase spectrum is generally used in the stochastic Green’s function method. Due to the fact that the ground motion time history corresponding to a random phase spectrum is a random wave, the directivity effect of the near-fault ground motion cannot be reflected in the synthesized results, and the simulation accuracy is lower in the low-frequency range. Therefore, based on the similarity between the phase difference spectrum and the ground motion time history envelope curve [29,30], the ground motion of the small earthquake in the first subfault satisfying the envelope curve is generated as follows:

(1): Using the time history envelope function which is based on scattering theory and proposed by Satoh [31,32], set the phase difference probability density distribution consistent with the shape of the envelope curve;
(2): Generate the phase difference randomly according to this phase difference probability density distribution;
(3): Derive the phase spectrum of the small earthquake in the first subfault according to the phase difference;
(4): After combining the phase spectrum with the amplitude spectrum, the ground motion of the small earthquake in the first subfault $F_{11}^{K} (f)$ in the frequency domain can be obtained.

Nozu et al. [33] found that seismic waves with the same propagation path have similar phase characteristics by analyzing the observed seismic records. Therefore, the idea of the modified empirical Green’s function method is used for reference, and Equation (7) is applied in this study to obtain the ground motions of the small earthquakes in other subfaults in the frequency domain:

F_{i j}^{K} (f) = A_{i j}^{K} (f) \frac{F_{11}^{K} (f)}{|F_{11}^{K} (f)|}

(7)

where,

F_{11}^{K} (f)

and

F_{i j}^{K} (f)

represent the ground motion (complex form) of the small earthquake caused by the rupture of the first subfault and the ij subfault, respectively;

|F_{11}^{K} (f)|

is the amplitude of

F_{11}^{K} (f)

and is smoothed using a Parzen window, with a bandwidth of 0.05 Hz.

Thus, the phase spectra of small earthquakes among different subfaults within a large fault can be correlated and have a certain degree of discreteness, and the directivity effect of the near-fault ground motion can be reflected. After performing the inverse Fourier transform on the complex Fourier spectra of small earthquakes in subfaults, the time histories of small earthquakes in subfaults can be obtained.

2.3. Small Earthquakes in Subfaults Superimposed to Form the Large Earthquake

After the P wave, SV wave, and SH wave of the small earthquakes in subfaults are independently generated, the time delay of their propagation in the crust and the rupture propagation on the fault plane is considered, and the three components of the waves are superimposed to generate the P wave, SV wave, and SH wave of the large earthquake, respectively. In this study, a ground motion synthesis method based on the scaling law is used [34,35], and the ground motion of the large earthquake U_K(t) at the bedrock is expressed as:

U_{K} (t) = \sum_{i = 1}^{N L} \sum_{j = 1}^{N W} f (t) * u_{i j}^{K} (t - t_{i j}^{K})

(8)

where

f (t) = δ (t) + \frac{a}{n' [1 - \exp (- a)]} \times \sum_{k = 1}^{(N D - 1) n'} \{\exp [- \frac{a (k - 1)}{(N D - 1) n'}] \cdot δ [t - \frac{(k - 1) τ_{i j}}{(N D - 1) n'}]\}

(9)

where δ(t) is Dirac function; a is taken as 1.0; τ_ij is the rise time of the ij subfault; NL, NW, and ND represent the division number of the length, width, and slip of the fault, respectively; n′ is the further subdivision number of the slip; and

t_{i j}^{K}

is the time delay when the ground motion of the ij subfault reaches the observation point. In order to avoid the occurrence of dominant frequencies due to the close interval between the rupture start times of subfaults, a random value is set for the rupture start time of each subfault [36]. The time delay is calculated according to Equation (10):

t_{i j}^{K} = \frac{r_{i j}}{v} + \frac{ξ_{i j}}{V_{R}} + \frac{(λ - 0.5) Δ W}{V_{R}}

(10)

where, ξ_ij is the distance from the rupture starting point on the fault plane to the center of the ij subfault; V_R is the velocity of rupture propagation on the fault plane; λ is a random number between 0 and 1; ΔW is the width of the subfault.

The rise time τ is generally calculated according to the commonly used Equation (11) [8]:

τ = 0.5 \frac{W}{V_{R}}

(11)

However, the rise time at a point on the fault plane is related to its position. Near the rupture starting point, the rise time is longer, while near the edge of the fault, the rise time is shorter. Equation (11) is only approximately true for the long, narrow faults [8]. According to the numerical simulation results of Day [8], for faults with various aspect ratios, there is a law that the longest rise time occurs at the rupture starting point, which equals approximately the width of the fault divided by the propagation velocity of the rupture. As the calculation point gradually moves away from the rupture starting point and approaches the edge of the fault, the rise time gradually shortens, and this attenuation law approximates a parabola.

Therefore, in order to improve the accuracy of the stochastic Green’s function method for ground motion simulation, based on the above laws, an improved approximate expression for the rise time on the rectangular fault considering that the rupture starting point may be at any position and the aspect ratio may be arbitrary is proposed in this study:

τ_{i j} = \frac{W}{V_{R}} [1 - {(\frac{x_{i j} - x_{R}}{C_{1}})}^{2}] [1 - {(\frac{y_{i j} - y_{R}}{C_{2}})}^{2}]

(12)

where, (x_R, y_R) and (x_ij, y_ij) stand for the local coordinates of the rupture starting point and the center of the ij subfault on the fault plane, respectively. The local coordinate system of the fault plane is oxy presented in Figure 1, o is the upper left corner of the fault plane, x and y are the direction of fault length and width, respectively. C₁ is the distance from the rupture starting point to the left or right edge of the fault along the direction of fault length when the center of the ij subfault is located on the left or right of the rupture starting point, respectively; C₂ is the distance from the rupture starting point to the top or bottom edge of the fault along the direction of fault width when the center of the ij subfault is located on the top or bottom of the rupture starting point, respectively:

C_{1} = \{\begin{matrix} L - x_{R}, & x_{i j} \geq x_{R} \\ x_{R}, & x_{i j} < x_{R} \end{matrix}

(13)

C_{2} = \{\begin{matrix} W - y_{R}, & y_{i j} \geq y_{R} \\ y_{R}, & y_{i j} < y_{R} \end{matrix}

(14)

In Equations (12)–(14), L and W represent the length and width of the fault, respectively. Equations (12)–(14) represent a functional form describing that the rise time would be longer when a subfault is closer to the hypocenter.

According to the above method, the time histories of the P wave, SV wave, and SH wave of the large earthquake are obtained. Assuming that seismic waves propagate in a straight line, the time histories of the P wave, SV wave, and SH wave of the large earthquake are projected and converted to the NS, EW, and UD directions according to Equation (15), respectively:

\{\begin{matrix} U_{NS} (t) = U_{P} (t) \sin i_{ξ} \cos ϕ + U_{SV} (t) \cos i_{ξ} \cos ϕ + U_{SH} (t) \sin ϕ \\ U_{EW} (t) = U_{P} (t) \sin i_{ξ} \sin ϕ + U_{SV} (t) \cos i_{ξ} \sin ϕ + U_{SH} (t) \cos ϕ \\ U_{UD} (t) = U_{P} (t) \cos i_{ξ} + U_{SV} (t) \sin i_{ξ} \end{matrix}

(15)

where, U_P(t), U_SV(t), and U_SH(t) represent the time histories of the P wave, SV wave, and SH wave of the large earthquake, respectively; U_NS(t), U_EW(t), and U_UD(t) represent the time histories of the large earthquake in NS, EW, and UD directions, respectively; ϕ is the azimuth angle from the center of the fault to the observation point; and i_ξ is the average takeoff angle. Since it is assumed that the seismic wave travels along a straight line, the average takeoff angle refers to the angle between the line connecting the center of the fault and the observation point and the vertical downward direction.

The ground motion obtained from the above simulation is the ground motion at the bedrock below the observation point. On this basis, it is necessary to consider the local site amplification effect to obtain the ground motion at the surface. The horizontal-to-vertical component spectral ratio (HVSR) method [37,38,39] is applied in this study to consider the local site amplification effect. The theory of the HVSR suggests that vertical ground motion consists of amplification from part of the source effects that need to be eliminated to achieve a suitable site amplification [38,40]. Therefore, the ratio of the Fourier spectrum of the horizontal component to that of the vertical component of the seismic records at the local site is taken as the local site amplification of the horizontal ground motion. To reduce the impact of the randomness of a single ground motion on the calculated results, it is necessary to take the average of the calculated results from three or more seismic records. The frequency range for site amplification is set to 0.2 Hz to 10 Hz [41,42,43]. The horizontal components of the simulated ground motion at the bedrock are multiplied by the Fourier spectral ratio, and the site amplification of the vertical component of the ground motion is ignored. Then the simulation results of the ground motion at the surface can be obtained.

The program code of the improved stochastic Green’s function method for 3D ground motion simulation is developed based on FORTRAN 90.

3. Validation of Effectiveness

3.1. Ground Motion Simulation of the 1994 Northridge Earthquake

The 1994 Northridge earthquake in the United States is taken as an example for analysis. The moment magnitude of the earthquake is Mw 6.7. The epicenter was located at 34.211° N, 118.546° W, and the focal depth was 17.5 km. Referring to the inversion source model of Wald et al. [44] and Wang [45], the fault parameters used in this study are listed in Table 1. The quality factor Q_K(f) is expressed as follows [46]:

Q_{K} (f) = 150 f^{0.5}

(16)

The fault location and station location of the 1994 Northridge earthquake are shown in Figure 2. In the figure, OXYZ is the overall coordinate system, O is the projection of the upper left corner of the fault plane on the surface, X, Y, and Z are the north, east, and vertical downward direction, respectively, oxy is the local coordinate system of the fault plane, o is the upper left corner of the fault plane, x and y are the direction of fault length and width, respectively, Q is the rupture starting point on the fault plane, and d_t is the top depth of the fault plane.

To test the effectiveness of the improved method proposed in this study, firstly, the source spectra of the 1994 Northridge earthquake are synthesized using the following two methods with subfault division numbers of 9 × 10, 18 × 20, and 36 × 40, respectively:

(1): Assuming that the rise time is a constant on the fault plane, the stochastic Green’s function method where the commonly used Equation (11) is used to calculate the rise time;
(2): Considering the non-uniform distribution of rise time on the fault plane, the stochastic Green’s function method where the proposed expression (12) is applied to calculate the rise time.

The source spectrum here refers to the Fourier amplitude spectrum of the moment-rate function (without considering the influence of high-frequency cutoff frequency) [36].

According to the ω⁻² model, the theoretical source spectra of the small and the large earthquake are represented by Equations (17) and (18), respectively [21]:

|{\dot{M}}_{0 i j} (f)| = \frac{M_{0 i j}}{1 + {({f / f}_{c i j}^{S})}^{2}}

(17)

|{\dot{M}}_{0} (f)| = \frac{M_{0}}{1 + {({f / f}_{c}^{S})}^{2}}

(18)

In which M₀ and

f_{c}^{S}

are the seismic moment and the corner frequency of the large earthquake, respectively.

The source spectrum of the large earthquake synthesized by superimposing the small earthquakes is expressed as:

{\dot{M}}_{0} (f) = \sum_{i = 1}^{N L} \sum_{j = 1}^{N W} G (f) {\dot{M}}_{0 i j} (f) \exp \{2 i π f [\frac{ξ_{i j}}{V_{R}} + \frac{(λ - 0.5) Δ W}{V_{R}}]\}

(19)

where G(f) is the conversion function of the source spectrum from the small earthquake to the subfault earthquake, expressed as:

G (f) = 1 + \frac{a [1 - \exp (- a - 2 i π f τ_{i j})]}{n' [1 - \exp (- a)] \{1 - \exp [\frac{- a - 2 i π f τ_{i j}}{(N D - 1) n'}]\}}

(20)

The source spectra of the large earthquake synthesized using the above two methods are compared with the theoretical ω⁻² source spectrum, as shown in Figure 3.

The results indicate that when assuming a uniform distribution of rise time on the fault plane, the sag phenomenon of the source spectrum in the intermediate frequency band becomes more and more significant as the subfault division number gradually increases, which is similar to the research results of Hisada [36]. However, when using the improved expression for the rise time, as the subfault division number increases, the amplitude of the source spectrum changes very little, and the sag phenomenon is relatively insignificant.

In order to further verify the effectiveness of the improved expression for the rise time proposed in this study, six stations are randomly selected around the 1994 Northridge earthquake fault. The rupture distances of these six stations range from 16 km to 47 km, with two stations located on the hanging wall and four stations located on the foot wall of the fault. In the cases of subfault division numbers of 9 × 10 and 36 × 40, the ground motions at these stations are simulated using the two methods mentioned above. The location, rupture distance R_rup, time-averaged shear-wave velocity of the upper 30 m V_s30, and the National Earthquake Hazards Reduction Program (NEHRP) site classification of these stations are listed in Table 2. According to the NEHRP code [47], the site types of these stations are class C sites, which are soft rock or very dense hard soil sites. The basic information and observed seismic records of these stations are from the strong ground motion database of the Pacific Earthquake Engineering Research Center.

In order to minimize the impact of randomness on the simulation results of the stochastic Green’s function method and verify the effectiveness of the improvement proposed in this study, 30 acceleration time histories are randomly generated for each simulation, and the 5% damped acceleration response spectra and the displacement Fourier amplitude spectra of the 30 acceleration time histories are calculated. The average value of the 30 acceleration response spectra and the average value of the 30 displacement amplitude spectra are taken as the simulated acceleration response spectrum and the simulated displacement amplitude spectrum, respectively, and the acceleration time history corresponding to the response spectrum closest to the average response spectrum is taken as the simulated acceleration time history.

As the acceleration time history of the ground motion is mainly concerned with the engineering seismic design, in this study, the acceleration time histories, acceleration response spectra, and displacement amplitude spectra of the synthesized 3D ground motions are compared with those of the observed records of the 1994 Northridge earthquake. Due to the uncertainty of fault parameters, propagation characteristics, and local site conditions, the simulated acceleration response spectra of the horizontal components are also compared with the prediction results obtained by the ground motion prediction equation (GMPE) ASK14 proposed by Abrahamson et al. [48], as shown in Equation (21):

\begin{array}{l} \ln S a (g) = f_{1} (M_{w}, R_{r u p}) + F_{R V} f_{7} (M_{w}) + F_{N} f_{8} (M_{w}) + F_{A S} f_{11} (C R_{j b}) + f_{5} ({\hat{S a}}_{1180}, V_{s 30}) + \\ F_{H W} f_{4} (R_{j b}, R_{r u p}, R_{x}, R_{y 0}, W, d i p, Z_{T O R}, M_{w}) + f_{6} (Z_{T O R}) + f_{10} (Z_{1.0}, V_{s 30}) + R e g i o n a l (V_{s 30}, R_{r u p}) \end{array}

(21)

The specific meanings represented by the parameters in the equation can be found in the literature [48].

The comparison results are shown in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9.

In Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, the comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the red line represents the GMPE, the yellow line represents the simulation result with the constant rise time when the subfault division number is 9 × 10, the green line represents the simulation result with the modified rise time expression when the subfault division number is 9 × 10, the brown line represents the simulation result with the constant rise time when the subfault division number is 36 × 40, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 36 × 40. Comparing the simulation results with the observed records at the stations, it can be seen that if the rise time is constant, the amplitude of the simulated response spectra and amplitude spectra will decrease in the intermediate frequency range. The amplitude decrease becomes more significant as the subfault division number increases. After using the modified rise time expression, the sag phenomenon that appeared in the intermediate frequency range is alleviated, and the simulation results are more consistent with the observed records. As the subfault division number increases, the decrease in amplitude is not as significant as the decrease when the rise time is constant. The simulation results of the horizontal components are compared with the results of the GMPE, and the results show that the response spectra simulated by the improved stochastic Green’s function method are generally close to the results predicted by the GMPE, further verifying the rationality of the simulation results.

From a sensitivity perspective, for the L4B station, the sensitivity of the simulation results to the subfault division number and the rise time expression is roughly similar in the period range of approximately 0.05 s to 0.5 s, while the simulation results show stronger sensitivity to the rise time expression than the subfault division number in the period range shorter than 0.05 s and longer than 0.5 s. For the L09 station, the simulation results show stronger sensitivity to the rise time expression than the subfault division number in almost the whole period range. For ATB, CHL, and LV1 stations, the simulation results show stronger sensitivity to the subfault division number than the rise time expression in the short period range, while the simulation results show stronger sensitivity to the rise time expression in the long period range, and the boundary is approximately 0.5 s, 2 s, and 0.8 s, respectively. For the HOW station, the sensitivity of simulation results to the subfault division number is stronger than the rise time expression in almost the whole period range. The change in sensitivity may be related to the azimuth and site effects of the station.

In order to measure the simulation accuracy of the improved method, the simulation misfit is defined as the average logarithm of the ratio of the simulated response spectrum to the response spectrum of the observed record over all stations in this study. The simulation misfit of the ground motion of the Northridge earthquake in the NS, EW, and UD directions is shown in Figure 10.

The results show that when the subfault division number is 9 × 10, the simulation accuracy is improved in the period range longer than 0.4 s after using the improved method. In the period range of 0.01 s to 6 s, the misfit values of the NS, EW, and UD components using the improved method are within ±0.35, ±0.3, and ±0.25, respectively. When the subfault division number is 36 × 40, the simulation accuracy is improved in the period range longer than 0.2 s after using the improved method. In the period range of 0.01 s to 6 s, the misfit values of the NS, EW, and UD components using the improved method are within ±0.5, ±0.35, and ±0.45, respectively. Overall, the simulation accuracy is higher when the subfault division number is 9 × 10 than when the subfault division number is 36 × 40.

3.2. Ground Motion Simulation of the 2013 Lushan Earthquake

In order to further verify the effectiveness of the improvement proposed in this study, the ground motion of the 2013 Lushan Ms 7.0 earthquake in China is simulated below. The epicenter of the Lushan earthquake was located at 30.3° N, 103.0° E, with a focal depth of 10.2 km. According to the inversion results of Wang et al. [49], the fault parameters used in this study are listed in Table 3. The quality factor Q_K(f) is expressed as follows [50]:

Q_{K} (f) = 274.6 f^{0.423}

(22)

The fault location and station location of the 2013 Lushan earthquake are shown in Figure 11.

Firstly, using the two methods described in the previous section, the source spectra of the 2013 Lushan earthquake are synthesized with subfault division numbers of 11 × 7, 22 × 14, and 44 × 28, respectively, and they are compared with the theoretical ω⁻² source spectrum, as shown in Figure 12.

The results indicate that when assuming a uniform rise time distribution on the fault plane, as the subfault division number increases, the sag phenomenon of the source spectrum in the intermediate frequency band becomes more and more significant, and the frequency range where the sag phenomenon occurs is also expanding. After improving the rise time expression, as the subfault division number increases, the amplitude change in the source spectrum is not significant.

Next, six stations are randomly selected around the 2013 Lushan earthquake fault, and the ground motions at these stations are simulated with subfault division numbers of 11 × 7 and 44 × 28, respectively. According to the results of Huang [51], Table 4 lists the basic information of these stations. The rupture distances of these six stations range from 4 km to 44 km, with four stations located on the hanging wall and two stations located on the foot wall of the fault. According to the NEHRP code [47], the site types of these stations are all class C sites.

In this study, the acceleration time histories, acceleration response spectra, and displacement amplitude spectra of the simulated 3D ground motions are compared with those of the observed records at these stations in the 2013 Lushan earthquake. The comparison results are shown in Figure 13, Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18.

In Figure 13, Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18, the black line represents the observed record, the yellow line represents the simulation result with the constant rise time when the subfault division number is 11 × 7, the green line represents the simulation result with the modified rise time expression when the subfault division number is 11 × 7, the brown line represents the simulation result with the constant rise time when the subfault division number is 44 × 28, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 44 × 28. The results indicate that if the rise time is constant, the sag phenomenon of the simulated response spectra and amplitude spectra in the intermediate frequency band will not be obvious when the subfault division number is small, whereas the amplitude in the intermediate frequency band will show a significant decrease when the subfault division number is large. After using the modified rise time expression, the sag phenomenon that appeared in the intermediate frequency range is alleviated, and the simulation results are closer to the observed records. As the subfault division number increases, the decrease in amplitude in the intermediate frequency range is not as significant as the decrease when the rise time is constant.

From a sensitivity perspective, for 51BXZ and 51BXM stations, the sensitivity of simulation results to the subfault division number is stronger than the rise time expression in almost the whole period range. For the 51YAM station, the simulation results show stronger sensitivity to the rise time expression than the subfault division number in the period range shorter than 0.4 s, while the simulation results show stronger sensitivity to the subfault division number in the period range longer than 0.4 s. For the 51LSF station, the simulation results show stronger sensitivity to the rise time expression than the subfault division number in almost the whole period range. For the 51YAL station, the sensitivity of the simulation results to the subfault division number and the rise time expression is roughly similar in the whole period range. For the 51HYT station, the sensitivity of the simulation results to the subfault division number and the rise time expression is roughly similar in the period range shorter than 0.3 s, while the simulation results show stronger sensitivity to the subfault division number in the period range longer than 0.3 s. The change in sensitivity may be related to the azimuth and site effects of the station.

In order to measure the simulation accuracy of the improved method, the simulation misfit of the ground motion of the Lushan earthquake in the NS, EW, and UD directions is calculated, as shown in Figure 19.

The results show that when the subfault division number is 11 × 7, the simulation accuracy of the NS and EW components is improved in the period range longer than 1.5 s after using the improved method, and the simulation accuracy of ground motion in the UD direction is improved in the whole period range. In the period range of 0.01 s to 6 s, the misfit values of the NS and EW components using the improved method are within ±0.2, and the simulation misfit of the UD component using the improved method is within ±0.3. When the subfault division number is 44 × 28, the simulation accuracy of the three components is significantly improved in most period ranges after using the improved method. In the period range of 0.01 s to 6 s, the misfit values of the NS, EW, and UD components using the improved method are within ±0.35, ±0.4, and ±0.6, respectively. Overall, the simulation accuracy is higher when the subfault division number is 11 × 7 than when the subfault division number is 44 × 28.

Overall, after using the improved rise time expression proposed in this study, the accuracy of the simulated response spectra and amplitude spectra has been improved in the broadband frequency range. When the subfault division number of the Northridge earthquake fault is 9 × 10, and the subfault division number of the Lushan earthquake fault is 11 × 7, the simulation results are in good agreement with the observed records. Therefore, the improved stochastic Green’s function method proposed in this study can better reproduce the ground motion in the near-fault region. However, due to the uncertainty of source parameters, propagation paths, and local site conditions, as well as the approximation of the simulation method, there are still some differences between the simulation results and the observed records [52,53].

In addition, it seems that the simulation results show stronger sensitivity to the rise time expression at stations with forward directivity effects, such as L09 station in the Northridge earthquake and 51LSF station in the Lushan earthquake, while the simulation results show stronger sensitivity to the subfault division number at backward stations, such as HOW station in the Northridge earthquake and 51BXZ and 51BXM stations in the Lushan earthquake. It needs to be verified with more earthquake events in the future.

4. Conclusions

In order to improve the accuracy of the stochastic Green’s function method in simulating near-fault broadband ground motion, based on the numerical simulation results of Day on the earthquake slip function, an improved approximate expression for the rise time on a rectangular fault considering that the rupture starting point may be at any position and the aspect ratio may be arbitrary is proposed in this study. On this basis, considering the contributions of P wave, SV wave, and SH wave to the 3D ground motion, respectively, an improved stochastic Green’s function method for simulating 3D broadband ground motion is proposed. Taking the 1994 Northridge earthquake in the United States and the 2013 Lushan earthquake in China as examples, under different subfault division numbers, the source spectra synthesized using the improved approximate expression for the rise time are compared with the theoretical ω⁻² source spectra of the large earthquake, then the ground motions at the observation points are simulated using the improved method, and the simulation results are compared with the observed records. After that, the effectiveness of the improved stochastic Green’s function method proposed in this study is verified. The following conclusions can be drawn:

(1): If it is assumed that the rise time is uniformly distributed on the fault plane, as the subfault division number gradually increases, the sag phenomenon of the source spectrum of the large earthquake synthesized by superimposing small earthquakes using the stochastic Green’s function method in the intermediate frequency band will become more and more significant. After improving the expression for the rise time, as the subfault division number increases, the amplitude change in the synthesized source spectrum of the large earthquake is not significant.
(2): If the rise time is assumed to be constant, the amplitude of the response spectrum and amplitude spectrum of the simulated ground motion at the observation point will decrease in the intermediate frequency range. As the subfault division number increases, the amplitude decrease becomes more pronounced. After using the improved rise time expression, the sag phenomenon that appeared in the intermediate frequency range is alleviated, and the simulation results are closer to the observed records. As the subfault division number increases, the decrease in amplitude in the intermediate frequency range is not as significant as the decrease when the rise time is constant.
(3): After using the improved rise time expression proposed in this study, when the subfault division number of the Northridge earthquake fault is 9 × 10, in the period range of 0.01 s to 6 s, the simulation misfit values of the NS, EW, and UD components are within ±0.35, ±0.3, and ±0.25, respectively; when the subfault division number of the Lushan earthquake fault is 11 × 7, in the period range of 0.01 s to 6 s, the simulation misfit values of the NS and EW components are within ±0.2, and the simulation misfit of the UD component is within ±0.3. The waveform of the time histories, the response spectra, and the amplitude spectra in the broadband frequency range simulated by the improved stochastic Green’s function method are in good agreement with those of the observed records, indicating that the improved method can simulate the near-fault broadband ground motion more accurately.

Author Contributions

Conceptualization, L.J. and X.X.; formal analysis, L.J. and X.P.; methodology, L.J.; resources, X.X.; software, X.X.; supervision, X.X.; validation, L.J. and X.P.; writing—original draft, L.J.; writing—review and editing, X.X. and X.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Natural Science Foundation of China grant number 52178174.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.13234766 (accessed on 18 October 2024).

Acknowledgments

Data for this study are provided by the Institute of Engineering Mechanics, China Earthquake Administration. The seismic records of the 1994 Northridge earthquake were downloaded from the Pacific Earthquake Engineering Research Center (https://peer.berkeley.edu/peer-strong-ground-motion-databases, last accessed on 1 December 2020).

Conflicts of Interest

The authors declare no competing interests.

References

Kamae, K.; Irikura, K.; Fukuchi, Y. Prediction of strong ground motion based on scaling law of earthquake by stochastic synthesis method. J. Struct. Constr. Eng. 1991, 56, 1–9, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Kamae, K.; Irikura, K.; Pitarka, A. A technique for simulating strong ground motion using hybrid Green’s function. Bull. Seismol. Soc. Am. 1998, 88, 357–367. [Google Scholar] [CrossRef]
Beresnev, I.A.; Atkinson, G.M. Modeling finite-fault radiation from the ωⁿ spectrum. Bull. Seismol. Soc. Am. 1997, 87, 67–84. [Google Scholar] [CrossRef]
Beresnev, I.A.; Atkinson, G.M. FINSIM—A FORTRAN program for simulating stochastic acceleration time histories from finite faults. Seismol. Res. Lett. 1998, 69, 27–32. [Google Scholar] [CrossRef]
Motazedian, D.; Atkinson, G.M. Stochastic finite-fault modeling based on a dynamic corner frequency. Bull. Seismol. Soc. Am. 2005, 95, 995–1010. [Google Scholar] [CrossRef]
Sedaghati, F. Simulation of Strong Ground Motions Using the Stochastic Summation of Small to Moderate Earthquakes as Green’s Functions; The University of Memphis: Memphis, TN, USA, 2018. [Google Scholar]
Boore, D.M. Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra. Bull. Seismol. Soc. Am. 1983, 73, 1865–1894. [Google Scholar] [CrossRef]
Day, S.M. Three-dimensional finite difference simulation of fault dynamics: Rectangular faults with fixed rupture velocity. Bull. Seismol. Soc. Am. 1982, 72, 705–727. [Google Scholar] [CrossRef]
Satoh, T. Generation method of broadband statistical green’s functions of horizontal and vertical ground motions. J. Struct. Constr. Eng. 2008, 73, 1087–1094, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Satoh, T. Improvement of generation method of three components of statistical Green’s functions considering near-field and intermediate-field terms. J. Struct. Constr. Eng. 2009, 74, 629–638, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Otarola, C.; Ruiz, S. Stochastic generation of accelerograms for subduction earthquakes. Bull. Seismol. Soc. Am. 2016, 106, 2511–2520. [Google Scholar] [CrossRef]
Ruiz, S.; Ojeda, J.; Pastén, C.; Otarola, C.; Silva, R. Stochastic strong-motion simulation in borehole and on surface for the 2011 Mw 9.0 Tohoku-Oki megathrust earthquake considering P, SV, and SH amplification transfer functions. Bull. Seismol. Soc. Am. 2018, 108, 2333–2346. [Google Scholar] [CrossRef]
Wang, H.W.; Ren, Y.F.; Wen, R.Z. A modified stochastic finite-fault method for simulating ground motions in three dimensions: A case study of Ludian earthquake. Earthq. Eng. Eng. Dyn. 2021, 41, 181–191, (In Chinese with English Abstract). [Google Scholar]
Ojeda, J.; Akinci, A.; Tinti, E.; Arriola, S.; Ruiz, S. Hybrid broadband strong-motion simulation to investigate the near-source characteristics of the M6.5, 30 October 2016 Norcia, Italy earthquake. Soil Dyn. Earthq. Eng. 2021, 149, 106866. [Google Scholar] [CrossRef]
Li, C.G.; Wang, H.W.; Wen, R.Z.; Qiang, S.Y.; Ren, Y.F. Three-component ground motion simulations based on the stochastic finite-fault method for the 2021 Maduo Ms 7.4 earthquake, Qinghai Province. Seismol. Geol. 2021, 43, 1085–1100, (In Chinese with English Abstract). [Google Scholar] [CrossRef]
Qiang, S.Y.; Wang, H.W.; Wen, R.Z.; Li, C.G.; Ren, Y.F. Three-dimensional ground motion simulations by the stochastic finite-fault method for the Yangbi, Yunnan Ms 6.4 earthquake on May 21, 2021. Chin. J. Geophys. 2021, 64, 4538–4547, (In Chinese with English Abstract). [Google Scholar] [CrossRef]
Luo, Z.F. Research on Stochastic Simulation Method of High Frequency Ground Motion Based on Improved Stochastic Finite Fault Method; Tianjin Chengjian University: Tianjin, China, 2022; (In Chinese with English Abstract). [Google Scholar]
Ji, L.F.; Xie, X.; Pan, X.Y. Improvement of stochastic Green’s function method for 3D broadband ground-motion simulation. Seismol. Res. Lett. 2023, 94, 331–349. [Google Scholar] [CrossRef]
Miyatake, T.; Yoshimi, M.; Terasaka, M. An approximate expression of slip rate for near-fault strong ground motion simulation, II. Application to strike slip fault and revision of the expression. Zisin (J. Seismol. Soc. Jpn. 2nd Ser.) 2003, 56, 125–139, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Hikita, T. Evaluation of earthquake ground motion by stochastic Green’s function method using rectangular crack model. J. Jpn. Assoc. Earthq. Eng. 2012, 12, 62–79, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Brune, J.N. Tectonic stress and the spectra of seismic shear waves from earthquakes. J. Geophys. Res. 1970, 75, 4997–5009. [Google Scholar] [CrossRef]
Brune, J.N. Correction to “Tectonic stress and the spectra of seismic shear waves from earthquakes”. J. Geophys. Res. 1971, 76, 5002, Erratum in J. Geophys. Res. 1970, 75, 4997–5009. [Google Scholar] [CrossRef]
Faccioli, E. A study of strong motions from Italy and Yugoslavia in terms of gross source properties. In Earthquake Source Mechanics; Geophysical Monograph Series, 37, Maurice Ewing; AGU: Washington, DC, USA, 1986; pp. 297–309. [Google Scholar] [CrossRef]
Nozu, A. A simple scheme to introduce near-field and intermediate-field terms in stochastic Green’s functions. In Proceedings of the 12th Japan Earthquake Engineering Symposium, Tokyo, Japan, 3–5 November 2006; JEES: Tokyo, Japan, 2006; pp. 190–193, (In Japanese with English Abstract). [Google Scholar]
Pitarka, A.; Somerville, P.; Fukushima, Y.; Uetake, T.; Irikura, K. Simulation of near-fault strong-ground motion using hybrid Green’s functions. Bull. Seismol. Soc. Am. 2000, 90, 566–586. [Google Scholar] [CrossRef]
Kotha, S.R.; Cotton, F.; Bindi, D. Empirical models of shear-wave radiation pattern derived from large datasets of ground-shaking observations. Sci. Rep. 2019, 9, 981. [Google Scholar] [CrossRef] [PubMed]
Boore, D.M.; Boatwright, J. Average body-wave radiation coefficients. Bull. Seismol. Soc. Am. 1984, 74, 1615–1621. [Google Scholar] [CrossRef]
Aki, K.; Richards, P.G. Quantitative Seismology, 2nd ed.; University Science Books: Sausalito, CA, USA, 2002. [Google Scholar]
Ohsaki, Y. On the significance of phase content in earthquake ground motions. Earthq. Eng. Struct. Dyn. 1979, 7, 427–439. [Google Scholar] [CrossRef]
Yamane, T.; Nagahashi, S. A study on a generation of simulated earthquake ground motion considering phase difference characteristics: Part 4 Generation method of simulated earthquake ground motion with the Fourier amplitude prescribed by the ω-square model and the Fourier phase selected based on the standard deviation of the phase difference. J. Struct. Constr. Eng. 2005, 70, 45–53, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Satoh, T. Study on envelope model of ground motions based on inversion of group delay time and scattering theory. J. Struct. Constr. Eng. 2004, 69, 71–78, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Satoh, T. Envelope model of acceleration ground motions for earthquakes caused by subducting Philippine Sea plate based on scattering theory. J. Struct. Constr. Eng. 2006, 71, 75–82, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Nozu, A.; Nagao, T.; Yamada, M. Simulation of strong ground motions using empirical site amplification and phase characteristics. Doboku Gakkai Ronbunshuu A 2009, 65, 808–813, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Irikura, K.; Kagawa, T.; Sekiguchi, H. Revision of the empirical Green’s function method. Programme Abstr. Seismol. Soc. Jpn. 1997, 2, B25. [Google Scholar]
Nozu, A. Necessity for a generalized expression of slip velocity correction function for use in empirical Green’s function method. Zisin (J. Seismol. Soc. Jpn. 2nd Ser.) 2002, 55, 233–238. (In Japanese) [Google Scholar] [CrossRef]
Hisada, Y. Causes and remedies of Fourier amplitude sags at intermediate frequencies for simulating broadband strong ground motions using stochastic Green’s function methods. J. Jpn. Assoc. Earthq. Eng. 2020, 20, 46–68, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Atkinson, G.M.; Boore, D.M. Earthquake ground-motion prediction equations for eastern North America. Bull. Seismol. Soc. Am. 2006, 96, 2181–2205. [Google Scholar] [CrossRef]
Nakamura, Y. A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface. Q. Rep. Railw. Tech. Res. Inst. 1989, 30, 25–33. [Google Scholar]
Yamazaki, F.; Ansary, M.A. Horizontal-to-vertical spectrum ratio of earthquake ground motion for site characterization. Earthq. Eng. Struct. Dyn. 1997, 26, 671–689. [Google Scholar] [CrossRef]
Nakamura, Y. What is the Nakamura method? Seismol. Res. Lett. 2019, 90, 1437–1443. [Google Scholar] [CrossRef]
Hata, Y.; Ichii, K.; Nozu, A.; Furukawa, A.; Tokida, K.I. Strong motion estimation at the embankment of the tomei expressway damaged by the 2009 suruga bay earthquake based on empirical site amplification and phase effects. Doboku Gakkai Ronbunshuu A 2010, 66, 673–690, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Hata, Y.; Nozu, A.; Kamai, T. Seismic waveform estimation at the Higashitakezawa landslide induced by the 2004 Mid Niigata Prefecture Earthquake based on empirical site amplification and phase effects. Landslides 2012, 49, 51–60, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Nagao, T.; Yamada, M.; Nozu, A. A study on estimation method of amplification factor at sites with thin sedimentary layer. J. Jpn. Soc. Civ. Eng. Ser. A1 (Struct. Eng. Earthq. Eng. (SE/EE)) 2010, 66, 1–11, (In Japanese with English Abstract). [Google Scholar] [CrossRef]
Wald, D.J.; Heaton, T.H.; Hudnut, K.W. The slip history of the 1994 Northridge, California, earthquake determined from strong-motion, teleseismic, GPS, and leveling data. Bull. Seismol. Soc. Am. 1996, 86, S49–S70. [Google Scholar] [CrossRef]
Wang, H.Y. Finite Fault Source Model for Predicting Near-Field Strong Ground Motion; Institute of Engineering Mechanics, China Earthquake Administration: Harbin, China, 2004; (In Chinese with English Abstract).
Hartzell, S.; Leeds, A.; Frankel, A.; Michael, J. Site response for urban Los Angeles using aftershocks of the Northridge earthquake. Bull. Seismol. Soc. Am. 1996, 86, S168–S192. [Google Scholar] [CrossRef]
Building Seismic Safety Council. NEHRP Recommended Seismic Provisions for New Buildings and Other Structures; National Institute of Building Sciences: Washington, DC, USA, 2015. [Google Scholar]
Abrahamson, N.A.; Silva, W.J.; Kamai, R. Summary of the ASK14 ground motion relation for active crustal regions. Earthq. Spectra 2014, 30, 1025–1055. [Google Scholar] [CrossRef]
Wang, W.M.; Hao, J.L.; Yao, Z.X. Preliminary result for rupture process of Apr. 20, 2013, Lushan Earthquake, Sichuan, China. Chin. J. Geophys. 2013, 56, 1412–1417, (In Chinese with English Abstract). [Google Scholar] [CrossRef]
Hua, W.; Chen, Z.L.; Zheng, S.H. A study on segmentation characteristics of aftershock source parameters of Wenchuan M8.0 earthquake in 2008. Chin. J. Geophys. 2009, 52, 365–371, (In Chinese with English Abstract). [Google Scholar] [CrossRef]
Huang, X.T. The Near-Field Strong Ground Motion Characteristics of Lushan Earthquake; Institute of Engineering Mechanics, China Earthquake Administration: Harbin, China, 2014; (In Chinese with English Abstract).
Beresnev, I.A.; Atkinson, G.M. Stochastic finite-fault modeling of ground motions from the 1994 Northridge, California, earthquake. I. Validation on rock sites. Bull. Seismol. Soc. Am. 1998, 88, 1392–1401. [Google Scholar] [CrossRef]
Wang, H.W.; Ren, Y.F.; Wen, R.Z. Source parameters, path attenuation and site effects from strong-motion recordings of the Wenchuan aftershocks (2008–2013) using a non-parametric generalized inversion technique. Geophys. J. Int. 2018, 212, 872–890. [Google Scholar] [CrossRef]

Figure 1. Stochastic Green’s function method for 3D ground motion simulation. (OXYZ is the overall coordinate system, O is the projection of the upper left corner of the fault plane on the surface, X, Y, and Z are the north, east, and vertical downward direction, respectively, oxy is the local coordinate system of the fault plane, o is the upper left corner of the fault plane, x and y are the direction of fault length and width, respectively, Q is the rupture starting point on the fault plane, and d_t is the top depth of the fault plane, red arrows indicate the vibration direction of waves).

Figure 2. Fault location and station location of 1994 Northridge earthquake. (OXYZ is the overall coordinate system, O is the projection of the upper left corner of the fault plane on the surface, X, Y, and Z are the north, east, and vertical downward direction, respectively, oxy is the local coordinate system of the fault plane, o is the upper left corner of the fault plane, x and y are the direction of fault length and width, respectively, Q is the rupture starting point on the fault plane, and d_t is the top depth of the fault plane. The black star marks the hypocenter, and the black triangles mark the locations of the stations.).

Figure 3. The source spectra of the 1994 Northridge earthquake. (The red solid line represents the source spectra of the large earthquake synthesized using the constant rise time, the blue solid line represents the source spectra of the large earthquake synthesized using the modified rise time expression, and the black dashed line represents the theoretical ω⁻² source spectrum.).

Figure 4. Simulation results of the L4B station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the red line represents the GMPE, the yellow line represents the simulation result with the constant rise time when the subfault division number is 9 × 10, the green line represents the simulation result with the modified rise time expression when the subfault division number is 9 × 10, the brown line represents the simulation result with the constant rise time when the subfault division number is 36 × 40, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 36 × 40).

Figure 5. Simulation results of the L09 station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the red line represents the GMPE, the yellow line represents the simulation result with the constant rise time when the subfault division number is 9 × 10, the green line represents the simulation result with the modified rise time expression when the subfault division number is 9 × 10, the brown line represents the simulation result with the constant rise time when the subfault division number is 36 × 40, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 36 × 40).

Figure 6. Simulation results of the ATB station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the red line represents the GMPE, the yellow line represents the simulation result with the constant rise time when the subfault division number is 9 × 10, the green line represents the simulation result with the modified rise time expression when the subfault division number is 9 × 10, the brown line represents the simulation result with the constant rise time when the subfault division number is 36 × 40, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 36 × 40).

Figure 7. Simulation results of the HOW station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the red line represents the GMPE, the yellow line represents the simulation result with the constant rise time when the subfault division number is 9 × 10, the green line represents the simulation result with the modified rise time expression when the subfault division number is 9 × 10, the brown line represents the simulation result with the constant rise time when the subfault division number is 36 × 40, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 36 × 40).

Figure 8. Simulation results of the CHL station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the red line represents the GMPE, the yellow line represents the simulation result with the constant rise time when the subfault division number is 9 × 10, the green line represents the simulation result with the modified rise time expression when the subfault division number is 9 × 10, the brown line represents the simulation result with the constant rise time when the subfault division number is 36 × 40, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 36 × 40).

Figure 9. Simulation results of the LV1 station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the red line represents the GMPE, the yellow line represents the simulation result with the constant rise time when the subfault division number is 9 × 10, the green line represents the simulation result with the modified rise time expression when the subfault division number is 9 × 10, the brown line represents the simulation result with the constant rise time when the subfault division number is 36 × 40, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 36 × 40).

Figure 10. The simulation misfit of the ground motion of the Northridge earthquake. (The yellow line represents the simulation misfit with the constant rise time when the subfault division number is 9 × 10, the green line represents the simulation misfit with the modified rise time expression when the subfault division number is 9 × 10, the brown line represents the simulation misfit with the constant rise time when the subfault division number is 36 × 40, and the blue line represents the simulation misfit with the modified rise time expression when the subfault division number is 36 × 40).

Figure 11. Fault location and station location of 2013 Lushan earthquake. (OXYZ is the overall coordinate system, O is the projection of the upper left corner of the fault plane on the surface, X, Y, and Z are the north, east, and vertical downward direction, respectively, oxy is the local coordinate system of the fault plane, o is the upper left corner of the fault plane, x and y are the direction of fault length and width, respectively, Q is the rupture starting point on the fault plane, and d_t is the top depth of the fault plane. The black star marks the hypocenter, and the black triangles mark the locations of the stations.).

Figure 12. The source spectra of the 2013 Lushan earthquake. (The red solid line represents the source spectra of the large earthquake synthesized using the constant rise time, the blue solid line represents the source spectra of the large earthquake synthesized using the modified rise time expression, the black dashed line represents the theoretical ω⁻² source spectrum.).

Figure 13. Simulation results of the 51BXZ station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the yellow line represents the simulation result with the constant rise time when the subfault division number is 11 × 7, the green line represents the simulation result with the modified rise time expression when the subfault division number is 11 × 7, the brown line represents the simulation result with the constant rise time when the subfault division number is 44 × 28, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 44 × 28).

Figure 14. Simulation results of the 51YAM station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the yellow line represents the simulation result with the constant rise time when the subfault division number is 11 × 7, the green line represents the simulation result with the modified rise time expression when the subfault division number is 11 × 7, the brown line represents the simulation result with the constant rise time when the subfault division number is 44 × 28, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 44 × 28).

Figure 15. Simulation results of the 51LSF station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the yellow line represents the simulation result with the constant rise time when the subfault division number is 11 × 7, the green line represents the simulation result with the modified rise time expression when the subfault division number is 11 × 7, the brown line represents the simulation result with the constant rise time when the subfault division number is 44 × 28, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 44 × 28).

Figure 16. Simulation results of the 51BXM station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the yellow line represents the simulation result with the constant rise time when the subfault division number is 11 × 7, the green line represents the simulation result with the modified rise time expression when the subfault division number is 11 × 7, the brown line represents the simulation result with the constant rise time when the subfault division number is 44 × 28, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 44 × 28).

Figure 17. Simulation results of the 51YAL station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the yellow line represents the simulation result with the constant rise time when the subfault division number is 11 × 7, the green line represents the simulation result with the modified rise time expression when the subfault division number is 11 × 7, the brown line represents the simulation result with the constant rise time when the subfault division number is 44 × 28, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 44 × 28).

Figure 18. Simulation results of the 51HYT station. (The comparisons of acceleration time histories are shown in the top row, the comparisons of acceleration response spectra are shown in the middle row, and the comparisons of displacement amplitude spectra are shown in the bottom row. The simulation results in the NS, EW, and UD directions are shown from left to right. The black line represents the observed record, the yellow line represents the simulation result with the constant rise time when the subfault division number is 11 × 7, the green line represents the simulation result with the modified rise time expression when the subfault division number is 11 × 7, the brown line represents the simulation result with the constant rise time when the subfault division number is 44 × 28, and the blue line represents the simulation result with the modified rise time expression when the subfault division number is 44 × 28).

Figure 19. The simulation misfit of the ground motion of the Lushan earthquake. (The yellow line represents the simulation misfit with the constant rise time when the subfault division number is 11 × 7, the green line represents the simulation misfit with the modified rise time expression when the subfault division number is 11 × 7, the brown line represents the simulation misfit with the constant rise time when the subfault division number is 44 × 28, and the blue line represents the simulation misfit with the modified rise time expression when the subfault division number is 44 × 28).

Table 1. Fault parameters of 1994 Northridge earthquake.

	Fault Parameter	Parameter Value
Macroscopic fault parameter	Strike, dip, rake	122°, 40°, 101°
	Fault length and width	18 km, 21 km
	Rupture area	378 km²
	Depth at the top of the fault	5.0 km
	Stress drop	5.0 MPa
	Seismic moment	1.3 × 10¹⁹ N·m
	Average slip	1.0185 m
	Medium density near the source	2.8 g/cm³
	S wave and P wave velocity near the source	3.7 km/s, 6.4 km/s
	Rupture propagation velocity	3.0 km/s
Microscopic fault parameter	Number of asperities	2
	Area of the first asperity	58.8 km²
	Slip of the first asperity	2.016 m
	Area of the second asperity	8.4 km²
	Slip of the second asperity	1.779 m
	Area of the background	310.8 km²
	Slip of the background	0.6618 m

Table 2. Station information of the 1994 Northridge earthquake. (R_rup is the rupture distance, V_s30 is the time-averaged shear-wave velocity of the upper 30 m).

Station	Latitude	Longitude	R_rup (km)	V_s30 (m/s)	NEHRP Site Classification
L4B	34.650	−118.477	31.69	523.54	C
L09	34.608	−118.558	25.36	670.84	C
ATB	34.758	−118.361	46.91	572.57	C
HOW	34.204	−118.302	16.88	581.93	C
CHL	34.086	−118.481	20.45	740.05	C
LV1	34.594	−118.242	37.19	499.31	C

Table 3. Fault parameters of 2013 Lushan earthquake.

	Fault Parameter	Parameter Value
Macroscopic fault parameter	Strike, dip, rake	205°, 38.5°, 88.8°
	Fault length and width	66 km, 35 km
	Rupture area	2310 km²
	Depth at the top of the fault	0 km
	Stress drop	6.0 MPa
	Seismic moment	3.61 × 10¹⁹ N·m
	Average slip	0.43 m
	Medium density near the source	2.8 g/cm³
	S wave and P wave velocity near the source	3.6 km/s, 6.2 km/s
	Rupture propagation velocity	2.88 km/s
Microscopic fault parameter	Number of asperities	3
	Area of the first asperity	120 km²
	Slip of the first asperity	0.613 m
	Area of the second asperity	300 km²
	Slip of the second asperity	1.069 m
	Area of the third asperity	210 km²
	Slip of the third asperity	0.78 m
	Area of the background	1680 km²
	Slip of the background	0.26 m

Table 4. Station information of the 2013 Lushan earthquake. (R_rup is the rupture distance, V_s30 is the time-averaged shear-wave velocity of the upper 30 m).

Station	Latitude	Longitude	R_rup (km)	V_s30 (m/s)	NEHRP Site Classification
51BXZ	30.5	102.9	15.6	393.9	C
51YAM	30.1	103.1	13.3	600.4	C
51LSF	30.0	102.9	4.1	517.4	C
51BXM	30.4	102.7	22.5	398.0	C
51YAL	29.9	102.8	14.9	535.0	C
51HYT	29.9	103.4	43.8	437.2	C

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Ji, L.; Xie, X.; Pan, X. Stochastic Green’s Function Method Considering Non-Uniform Rise Time Distribution to Simulate 3D Broadband Ground Motion. Appl. Sci. 2024, 14, 9796. https://doi.org/10.3390/app14219796

AMA Style

Ji L, Xie X, Pan X. Stochastic Green’s Function Method Considering Non-Uniform Rise Time Distribution to Simulate 3D Broadband Ground Motion. Applied Sciences. 2024; 14(21):9796. https://doi.org/10.3390/app14219796

Chicago/Turabian Style

Ji, Longfei, Xu Xie, and Xiaoyu Pan. 2024. "Stochastic Green’s Function Method Considering Non-Uniform Rise Time Distribution to Simulate 3D Broadband Ground Motion" Applied Sciences 14, no. 21: 9796. https://doi.org/10.3390/app14219796

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu