In this study, we use the stochastic finitefault method with a dynamic corner frequency concept as implemented by Motazedian and Atkinson (2005) in the EXSIM software. This is used to model the given earthquake and then estimate the ground motion at KUMA. The fault plane is modeled as a combination of smaller subfaults which are all treated as stochastic point sources. The acceleration spectra from each point source (i,j) are expressed as follows in the frequency domain:
$${A}_{ij}\left(f\right)=\frac{C{M}_{{0}_{ij}}{H}_{ij}{\left(2\pi f\right)}^{2}}{1+{\left(\frac{f}{{f}_{{c}_{ij}}}\right)}^{2}}G\left(R\right){e}^{\frac{\pi f{R}_{ij}}{Q\beta }}Z\left(f\right){e}^{\left(\pi \kappa f\right)}$$
1
$$C=\frac{{\mathfrak{R}}^{\theta \gamma }\bullet FS\bullet PRTITN}{4\pi \rho {\beta }^{3}}$$
2
where M0 is the seismic moment in dyne∙cm, Rij is the distance from the observation point to the ijth subfault, β is crustal shearwave velocity in km/s, Q is frequencydependent quality factor term, Z(f) is the soil amplification function, kappa (κ) is a factor which models the linear decay in higher frequencies of Fourier Amplitude Spectrum of Swave portion of the acceleration records in semilogarithmic space, FS is the free surface amplification factor whose value is generally assumed to be 2, PRTITN is a factor applied to reflect the effect of shearwave energy partitioning into two horizontal components and its value is generally assumed to be 1/√2, ρ is the crustal density in g/cm3, and Hij is a frequencydependent scaling factor for highfrequencies. The radiation pattern constant Rθγ is mostly taken as 0.55 for shear waves. In this method, an average horizontal component is obtained. The corner frequency of a subfault, fcij in Eq. (1) is defined as:
$${f}_{{c}_{ij}}={N}_{R}{\left(t\right)}^{\frac{1}{3}}\text{*}4.9\text{*}{10}^{6}\beta {\left(\frac{\varDelta \sigma }{{M}_{0\_ave}}\right)}^{\frac{1}{3}}$$
3
where, NR (t) is the cumulative number of ruptured subfaults at time t, ∆σ is the stress drop and M0_ave is the average seismic moment on the fault.
The deterministic acceleration spectra defined in Eq. (1) is combined with random phases and transferred into the time domain for each point source on the fault plane. The contribution from each subfault is summed up in the time domain to yield the following total acceleration:
$$a\left(t\right)=\sum _{i=1}^{nl}\sum _{j=1}^{nw}{a}_{ij}\left(t\varDelta {t}_{ij}{T}_{ij}\right)$$
4
where, aij is the acceleration time series corresponding to the ijth subfault and a(t) is the acceleration of the entire fault. The terms nl and nw are, respectively, the number of subfaults considered along the length and width of the rectangular fault plane. Tij is the ratio of the subfault radius to the rupture velocity, ∆tij is the time delay of each subfault to the observation point.
Simulation Input Parameters For Analyses In Esg6 – Blind Trial Steps 2 And 3
Step 2
Accurate simulated ground motions simulations rely on welldefined simulation parameters which effectively represent the regional characteristics. In this study, the source parameters of the Step 2 earthquake (Mw = 5.5) are obtained from ESG6 – Blind Trial Step 2 documents. The regional parameters for anelastic attenuation, geometrical spreading, and duration are obtained from Zhang et al. (2016). The simulation parameters of the Mw = 5.5 ESG6 – Blind Trial Step 2 event are shown in Table 1.
Table 1
Input parameters for ground motion simulations of the Step 2 Mw = 5.5 event at SEVO and KUMA stations
Parameter

Value

Moment Magnitude

5.5

Hypocenter Location

32.9638°N 131.0868°E

Hypocenter Depth

5 km

Fault Orientation

Strike:209\(^\circ\), Dip: 60\(^\circ\)

Slip distribution

Random slips

Fault Dimensions

8 x 4 km (from empirical equations of Wells and Coppersmith, 1994)

Attenuation Model

\(Q\left(f\right)=85.5{f}^{0.68} \text{S}\text{E}\text{V}\text{O}\)
\(Q\left(f\right)=120{f}^{0.64} \text{K}\text{U}\text{M}\text{A}\)

Duration Model

0.0 (R: 0–10 km)
+ 0.16 (R: 10–70 km)
−0.03 (R: 70–130 km)
0.04 (R > 130 km)

Geometric Spreading

\({R}^{1} R\le 100 km\)
\({R}^{0.5} R>100 km\)

Crustal Shear Wave Velocity

3600 m/s SEVO
3700 m/s KUMA

Crustal Density

2800 kg/m3

Stress Drop

10 bar

Pulsing Area Percentage

50%

Amplification Factors

Generic rock amplification factors at SEVO Thereotical 1D transfer functions at KUMA

In the simulations, at SEVO station, generic rock conditions (Boore and Joyner, 1997) are assumed and generic rock site amplification factors are used along with a zerodistance kappa value of κ0 = 0.035. After validation at the reference station SEVO, the ground motion at KUMA station is amplified in the frequency domain by generic site amplification factors using the site class NEHRPD (Boore and Joyner, 1997). This site class has been assigned according to the Vs30 value computed from the 1D soil model described at KUMA. The soil model at KUMA (Fig. 2) is based on the velocity model obtained from Asten et al. (this issue) as part of ESG6 – Blind Trial Step 1 result. At KUMA station, a κ0 value of 0.051 is adapted from the regional study of Zhang et al. (2016).
Step 3
Step 3 involves simulations of the 2016 Kumamoto foreshock (Mw = 6.5) and Kumamoto mainshock (Mw = 7.0).. In this step, similar to Step 2, at SEVO station, generic rock conditions are assumed and generic site amplification factors (Boore and Joyner, 1996) are used along with a κ0 = 0.035. For both the foreshock and the mainshock of the 2016 Kumamoto earthquake, after a coarse validation at the reference station SEVO (despite the problematic data at SEVO for the Mw = 7.0 event), the hard rock ground motion at KUMA is computed with the validated simulation parameters. Then, it is propagated through the 1D soil model for KUMA (Fig. 2) to obtain the surface acceleration at KUMA with an equivalent linear approach (Kramer, 1996). The proposed soil model given in Table 2 at KUMA is formed using the shearwave velocity model obtained from Asten et al. (this issue) as part of ESG6 – Blind Trial Step 1 result. The model shown in Fig. 2 is combined with information provided in Step 3 (Kumamoto Earthquake Ground Structure Survey). The simulation parameters of the Mw = 6.5 Kumamoto foreshock and Mw = 7.0 Kumamoto mainshock are summarized in Table 3.
Table 2 Proposed Soil Model at KUMA station
Layer Number

Material Name

Thickness
(m)

Unit Weight
(kN/m3)

Vs
(m/sec)

G Max
(MPa)

Soil Model ID

PI (%)

1

Bank (Gravel)

2.00

17.7

161.0

46.7

Rollins et al. (2020)

0.00

2

Sandy Silt

12.00

18.6

181.0

62.3

Vucetic and Dobry (1991)

47.90

3

Silt with Sand
(top layer incudes clay with sand)

12.30

19.6

170.0

57.8

Vucetic and Dobry (1991)

120.50

4

Gravel followed by Tuff breccia

75.00

20.6

470.0

464.0

Rollins et al. (2020)

0.00

5

Rock

251.00

23.0

980.0

2252.5

Rock (Idriss)

0.00

6

Rock

400.00

23.0

1210.0

3433.8

Rock (Idriss)

0.00

7

Rock

498.00

24.0

2500.0

15295.8

Rock (Idriss)

0.00

8

Rock

0.00

24.0

2700.0

17841.0

Rock (Idriss)

0.00

Table 3 Input parameters for ground motion simulations of Step 3 Kumamoto Mw=6.5 foreshock and Mw=7.0 mainshock at SEVO and KUMA stations
Parameter

Value (Foreshock)

Value (Mainshock)

Moment Magnitude

6.5

7.0

Hypocenter Location

32.7417°N 130.8087°E

32.7545°N 130.7630°E

Hypocenter Depth

11.39 km

12.45

Fault Orientation

Strike:212\(^\circ\), Dip: 89\(^\circ\)

Strike:205\(^\circ\), Dip: 72\(^\circ\)

Slip distribution

Asano and Iwata (2016)

Fault Dimensions

14 x 13 km

42 x 18 km

Attenuation Model

\(Q\left(f\right)=85.5{f}^{0.68} \text{S}\text{E}\text{V}\text{O}\)
\(Q\left(f\right)=120{f}^{0.64} \text{K}\text{U}\text{M}\text{A}\)

Duration Model

0.0 (R: 0–10 km)
+ 0.16 (R: 10–70 km)
−0.03 (R: 70–130 km)
0.04 (R > 130 km)

Geometric Spreading

\({R}^{1} R\le 100 km\)
\({R}^{0.5} R>100 km\)

Crustal Shear Wave Velocity

3600 m/s SEVO
3700 m/s KUMA

Crustal Density

2800 kg/m3

Stress Drop

25 bars

50 bars

Pulsing Area Percentage

50%

Amplification Factors

Generic rock amplification factors at SEVO Thereotical 1D transfer functions at KUMA

Ground Motion Simulation Results
For the Mw = 5.5 event in Step 2, initially, the horizontal ground motion component simulated at SEVO station is shown in comparison with the observed NS and EW components in Fig. 3 for both time and frequency domains. It is well known that the stochastic method does not accurately compute the arrival times as it contains only the Swaves. Thus, the time history is shifted to compare the Swave portions in the time domain. Figure 3 shows that the Swave portion is simulated closely in time domain. For FAS comparisons, the bottom left panel indicates the comparison of the Swave portions of the simulated and recorded time histories while the right panel shows the same comparison with the full waveforms. It can be observed that the Swave portion of the records is simulated effectively. For making not only visual but quantitative validations, the Goodness of fit (GOF) criteria of Olsen and Mayhew (2010) is used to evaluate the simulated motions at SEVO station can be found in Table 4. The definition of these GOF scores are as follows: 80–100 is excellent fit, 65–80 is very good fit, 45–65 is fair fit, and 35–45 is poor fit. In the computation of GOF scores, the main parameters including FAS, PGA, PGV, PGV/PGA, Housner Intensity, and Arias Intensity are employed. Table 4 indicates that the simulated motion of the Mw = 5.5 event at SEVO provides a GOF score of 79 which is close to being an excellent fit.
Event

GOF score

Table 4
Goodness of fit scores at SEVO station for the Mw = 5.5 target event, 2016 Kumamoto foreshock and the mainshock
Step 2: Target earthquake (Mw = 5.5)

79

Step 3: 2016 Kumamoto foreshock (Mw = 6.5)

45

Step 3: 2016 Kumamoto mainshock (Mw = 7.1)

56

Next, we simulate the horizontal component at KUMA station for the Mw = 5.5 event as shown in Fig. 4. The reliable frequency range for the simulated data with the stochastic approach is generally considered to be 0.2525 Hz.
As part of the ESG6 Blind Trial Step 3, initially, the horizontal ground motion components for the 2016 Kumamoto foreshock and mainshock are simulated at SEVO station. The results are shown in comparison with the corresponding observed NS and EW components in Figs. 5 and 6, for two events respectively, in both time and frequency domains. For the foreshock, the simulated motion underestimates the low frequencies but provide a closer fit for the higher frequency bands. The stochastic method cannot simulate the longperiod surface waves thus the discrepancy of the low frequency content can be attributed to this inherent limitation. Time domain comparisons are observed to simulate the Swave content closely. Similarly, the GOF criteria is presented in Table 4 for the foreshock and the mainshock. The GOF values indicate that the fits for the two events are both computed to be fair fits consistent with the visual inspections in time and frequency domains.
Finally, we simulate the horizontal components of the two events at KUMA station as shown in Fig. 7. It is observed that for the Kumamoto mainshock we observe longer durations with a broadband frequency content including higher frequencies when compared to the foreshock.