## 3 Physics-Based Stochastic Finite-Fault Simulation  
### Stochastic Finite-Fault Method  
In this study, the stochastic finite-fault simulation method proposed by [19] was adopted. The basic assumption of this method is that the near-fault ground motion of a target site [\(a\)(t)] can be simulated by the summation of \(e\) and \(m\) subfaults along the rupture fault length and width, respectively, with a proper time delay in the time domain. \[a(t)=\sum_{i=1}^{n}\sum_{j=1}^{m}a_{ij}(t+\Delta t) \tag{14}\]  
where \(a_{\rm ij}(t+\Delta t_{ij})\) is the acceleration time series generated by the \(ij\)th subfault to reach the target site with a relative delay time \(\Delta t_{ij}\). The acceleration Fourier spectrum for a subfault [\(A_{ij}(f)\)] at a hypocentral distance \(R_{\rm ij}\) is modeled as a point source with a Brune's source spectrum [50; 51]. \[A_{ij}(f)=H_{ij}E_{ij}(f)P\big{(}R_{ij},f\big{)}S(f) \tag{15}\]  
in which  
\[H_{ij}(f)=\sqrt{\frac{N\Sigma\{f^{2}/[1+(ff/f_{0})^{2}\}}{\Sigma\{f^{2}/[1+(ff /f_{0 ij})^{2}\}}} \tag{16}\]  
\[E_{ij}(f)=\frac{(2\pi f)^{2}C\Omega_{0 ij}}{[1+(f/f_{0 ij})^{2}]} \tag{17}\]  
\[P(R_{ij},f)=G(f)\exp[\frac{-\pi fR_{ij}}{Q(f)\beta}] \tag{18}\]  
\[S(f)=O(f)\exp(-\pi f\kappa_{0}) \tag{19}\]  
where \(H_{ij}(f)\) is a scaling factor to conserve the high-frequency spectral level of subfaults, \(E_{ij}(f)\) is the Brune's source spectrum, \(P(R_{\rm ij},f)\) is the propagation path attenuation, and \(S(f)\) is the site effect of the target site. \(N\) is the ratio of the total energy at high frequencies from the entire fault and from the \(ij\)th subfault. \(f_{0}\) is the source corner frequency of the earthquake, given by \(f_{0}\) = 4.9 \(\times\) 10\({}^{6}\beta(\Delta\sigma/M_{0})^{1/3}\), in which \(\Delta\sigma\) is the stress drop. \(f_{0 ij}\) is  
model and geometrical spreading on the simulation is ignorable. Thus, we used this \(Q\)(\(f\)) model and geometrical spreading without considering its uncertainty. As shown in Figure 7a, the crustal velocity model of the target site is constructed via the borehole data and the crustal velocity model [56]. The \(V_{\text{S30}}\) of the target site is 510 m/s, and the soil sediment thickness is 12 m. Thus, the site class of the large site is site class B following the GeoMatrix site classification proposed by [PERSON] et al. [57]. Then, the crustal amplification function \(O\)(\(f\)) is calculated by the quarter-wavelength approximation method proposed by [PERSON] and [PERSON] [53] (Figure 7b). Previous studies have proven that the \(O\)(\(f\)) is a first-order approximation of site amplification without considering the site response fluctuation induced by refraction, reflection, and diffraction of seismic wave energy in the layered site medium [55; 58]. Because there is no ground-motion record of the target sites, the \(\kappa_{0}\) is estimated by eight site class B stations in [PERSON] et al. [55] (for detailed site information, please refer to their Table 1). To obtain the uncertainty of \(\kappa_{0}\) value of site class B, the \(\kappa_{0}\) value of each record is calculated by  
\[\kappa_{0 pq}=\kappa_{q}-m_{p}R_{pq} \tag{20}\]  
where \(\kappa_{q}\) is the high-frequency attenuation parameter of \(q\)th record, \(m_{\text{p}}\) is the average distance-related parameter of the \(p\)th station, and \(R_{pq}\) is the epicentral distance. Then, the \(\kappa_{0}\) value of site class B is estimated by  
\[\kappa_{0}=\frac{\sum_{N_{r}}\kappa_{0}}{N_{r}} \tag{21}\]  
where \(N_{r}\) is the total number of records. Figure 8 shows the histogram distribution of the \(\kappa_{0}\) values of site class B (in 0.005 s bins). The result shows a positively skewed distribution. However, because only 88 accelerograms recorded at eight stations are available, in this study, we simply adopted the average \(\kappa_{0}\) value with SD, that is, 0.0288 \(\pm\) 0.00686 s. The weight of the average \(\kappa_{0}\), the \(+1\) SD, and the \(-1\) SD are 0.65, 0.14, and 0.21, respectively, depends on the data distribution.