There are a number of approaches in literature for simulating the ground motion time histories from a near seismic source, which are categorized into three main groups: 1) Physics-based approaches, 2) Stochastic methods, and 3) Hybrid approaches. A summary review of them with their main features can be found in the study of Douglas and Aochi (2008). In the present study, the stochastic finite-fault approach with dynamic corner frequency developed by Motazedian and Atkinson (2005) is employed. This approach integrates the characteristics of the earthquake source, path, directivity and site effects through a statistical approach. This method has an appropriate computation efficiency and widely used in simulating the ground motion waveforms in areas where past earthquake data is scarce. Figure 2 shows a general analytical framework of the stochastic finite-fault method. As depicted, in the first step, the rupture area is divided into sub-faults which are considered as point sources. In the second step, the ground motion waveform of sub-faults is derived. Finally, the ground motion from the rupture of the whole of fault is obtained by superimposing the waveform of sub-faults with an appropriate time delay. Eq. (1) is a mathematical representation of above explanations.

$$a\left(t\right)= \sum _{i=1}^{nl}\sum _{j=1}^{nw}{a}_{ij}(t+{\varDelta }_{ij})$$

1

In the above equation, \(nl\) and \(nw\) represent the number of sub-faults along the strike and dip of the fault, respectively. The acceleration time history from each element is denoted by \({a}_{ij}\) and \({\varDelta }_{ij}\) is the delay time which is a function of the distance between the source and observation point.

The waveform of sub-faults \({a}_{ij}\) are derived from the approach proposed by Motazedian and Atkinson (2005). The middle box of Fig. 2 shows a schematic analytical flow of this approach. As depicted, a Gaussian white noise signal is first generated. This is multiplied by a time window function. Generally, the time window function proposed by Saragoni-Hart is used in the simulation due to its appropriate representation of the average envelope of the squared acceleration time series. It is not, however, the ideal shape for earthquakes with two or more consecutive ruptures (Zhou and Chang, 2019). Then, the normalized Fourier spectrum of the time window signal is multiplied by the source amplitude in frequency domain. The source spectrum amplitude comes from the integration of source, path and site response. Finally, the simulated ground motion waveform is obtained by converting the signal to the time domain.

The mathematical form of Fourier amplitude of sub-fault is presented in Eq. (2). A detailed description of parameters in the following equation is presented in Table 1.

Table 1 A summary description of parameters for simulating the ground motion waveform in equation (1)

$${A}_{ij}\left(f\right)=C{{H}_{ij}M}_{0ij}\frac{{\left(2\pi f\right)}^{2}}{1+{\left(\frac{f}{{f}_{0ij}}\right)}^{2}}{e}^{-\pi {{f}_{k}}_{0}}\frac{{e}^{-\frac{\pi f{R}_{ij}}{Q\left(f\right)\beta }}}{{R}_{ij}}G\left({R}_{ij}\right)$$

2

Proper calibration of influencing factors such as source, rupture area, path and site parameters is an important step in stochastic finite-fault approach. Due to the absent of instrumentally records on the NTF, parameters are taken from literature. As illustrated in Fig. 3, the NTF seismic scenario has an approximate length of 110 Km with two main segments (Berberian et al, 1983). Here, the western segment, which is closer to Karaj, is considered as worst-case seismic scenario. Based on the empirical relation of Wells and Coppersmith (1994), the western segment with an approximate length of 57 Km has the potential of producing an earthquake with magnitude of 7.1 (Mw). This moment magnitude corresponds to a rupture plane with length of 46 and 26 Km along the fault's strike and dip, respectively. The fault plane is schematically depicted in Fig. 3. According to Samaie et al (2012) the strike and dip of the NTF fault are 305o and 35o, respectively. Because there are no instrumental strong earthquake data on the NTF, providing an accurate estimation of stress drop is associated with uncertainties. In this study, the value of stress drop is fixed to 50 bar, as proposed by Zafarani et al (2009). The quality factor in the present study is also taken from the study of Zafarani et al (2009). They proposed the form of \(Q\left(f\right)=87{f}^{1.47}\) for earthquakes in the north of Iran. This value is derived from the Manjil-Roudbar earthquake (7.4 Mw, 1990). Similarly, Motazedian (2006) developed a trilinear geometric spreading function based on data from the Manjil-Roudbar earthquake. Here, the same model is used. The proportion of pulsing area is assumed to be 50%, implying that during sub-fault rupture, no more than 50% of all sub-faults can be active, while the remainder are passive. Matazedina and Atkinson (2005) proposed this value based on the concept of self-healing provided by Heaton (1990). As there is no information on the slip distribution, a random normally slip distribution is employed for analysis.

In the present study, the analysis is carried out using the EXSIM code written in FORTRAN (Motazedian and Atkinson, 2005). The analysis is performed in a 0.5x0.5 Km grid cell. Figure 4 shows the results based on the above parameters on the bed-rock within Karaj's 12 municipal districts. As shown, the PGA values range from 0.1 to 0.42g, with the highest values in the north of Karaj nearby the NTF fault. Figure 5 presents the simulated time history waveform in 5 locations of Karaj (yellow triangle in Fig. 4). As shown, time history waveform of points in the north of Karaj has greater amplitude; while, time history waveform of points in the south of Karaj (which are far from NTF fault) has a weaker amplitude. This is in accordance with our expectations. Figure 5 also shows the pseudo response acceleration for the same points. Similarly, points located in north of Karaj has higher pseudo response acceleration.

To obtained the soil surface acceleration the VS30 map of region is required. The VS30 map is extracted from a microzonation study done by International Institute of Earthquake Engineering and Seismology (IIEES, 2013). Figure 6 shows the VS30 map of Karaj. As depicted, the northern part of the city, near the foothill of Alborz Mountains, has higher VS30; while, the southern part of the city has lower VS30. Figure 7 shows distribution of PGA values on the soil. As depicted, by employing the local site condition, the acceleration in the southern and central parts of the city is increased (i.e., district number 1, 2, 4 and 7). The existence of soft soil in these regions is the primary cause of such amplification.

To compare the results with GMPEs, a deterministic analysis is performed using the open-source OpenQuake program (Pagani et al, 2014). The epistemic uncertainty of GMPEs is considered in analysis using logic tree. To do so, five NGA-west2 relations with equal weights are applied (Bozorgnia et al, 2014). Figure 8 shows distribution of PGA on the soil surface for the same seismic scenario. As depicted, there is appropriate consistency between the values of stochastic finite-fault model and GMPEs. Figure 9 shows a grid-by-grid comparison of results from GMPEs and simulation. As depicted, the PGA values from the stochastic finite-fault are slightly higher in short distances; while, in far distances the finite-fault approach shows lower acceleration in relation to GMPEs.