Sensitivity analysis of input ground motion on surface motion parameters in high seismic regions: a case of Bhutan Himalaya

Historical earthquakes demonstrate that strong motion characteristics and local soil condition, when coupled, significantly influence seismic site response. Interestingly, most of the Himalayan earthquakes depicted anomalous behavior per the site conditions historically. Being one of the most active seismic regions on earth, the eastern fringe of the Himalaya has observed many devastating earthquakes together with non-uniform damage scenarios. To quantify such anomalies, we evaluate surface motion parameters for a soft soil deposit located in the city of Phuentsholing in western Bhutan. Using one-dimensional site response analysis, the sensitivity of ground motion variation is estimated. This study accounts for the earthquakes of moment magnitudes 6.6 to 7.5 with a wide variation in peak ground acceleration (PGA). To dissect the characteristics of six inputted ground motions on eight local ground conditions, a sensitivity analysis is performed statistically. The statistical correlation of the response datasets and the linear regression model of the bedrock outcrop and the surface motion spectral acceleration along the stratified depth are examined to quantify the variation in surface motion parameters. The results highlight that the strong motions with PGA greater than 0.34 g demonstrate greater sensitivity, leading to some anomalies in response parameters, especially amplification. Similar results were obtained for the low PGA range (<0.1 g).

1 Introduction

Bhutan is located in the eastern fringe of the Hindu Kush Himalaya. Historical earthquakes that occurred in the Hindu Kush Himalayan region have resulted in enormous losses and damage (Gautam et al., 2016). Akin to the historical earthquakes, the impending earthquakes are certain to strike the region and result in detrimental consequences. The eastern fringe of the Himalaya, i.e., Bhutan, and neighboring areas were strongly shaken by significant earthquakes in the past; however, most of the earthquakes that occurred up to the 18th century are not well documented. The most recent events such as on 5 April 2021 (M_w=5.0) in Samtse (South Bhutan) and the September 2009 Mongar earthquake (M_w=6.7) in eastern Bhutan resulted in widespread damage to Bhutan and neighboring regions. These earthquakes caused major damage in the eastern parts of Bhutan and considerably affected the other parts of the country (Chettri et al., 2021a, b). All the past earthquakes highlighted anomalous damage patterns to structures and infrastructures in various parts of the country, especially in the plain areas. This evidence indicates the likely local site effect in Bhutan. So far, few studies on local seismic response have been in Bhutan, using a single strong motion record, but the reported studies mainly focus on the role of bedrock depth in ground response parameters (Tempa et al., 2020; Tempa et al., 2021a, b). The ground motion response analysis may not adequately address the accuracy in predicting the response when the information is limited regarding site characteristics and their variations within the same soil column (Stevens et al., 2020). In the case of data-scarce regions such as Bhutan, the variation in terms of material characteristics can be possibly accounted for using sensitivity analysis. For this reason, this study quantifies the characteristics and effects of several strong ground motions. Seismic ground response analysis fall in the grade III approach of microzonation studies (e.g., ISSMGE, 1999; Licata et al., 2019). It is a widely used method by researchers for various applications in order to capture local ground effects or site conditions that can affect the estimate of ground motion characteristics (Chavez-Garcia et al., 1990; Lopez-Caballero et al., 2012; Gautam and Chamlagain, 2016; Sil and Haloi, 2018). The outcomes of such studies aim to provide local seismic (Gautam, 2017) hazard parameters which can be adopted for the design of structures and infrastructures (Douglas, 2006). Ground response parameters typically characterize the complex nature of strong motion accelerograms using a simple expansion of predictive relationships. Two prominent approaches, deterministic and probabilistic, are widely used for seismic hazard studies. Tempa et al. (2021a) recommended the use of the deterministic approach that can estimate the parameters under various earthquake occurrence scenarios. Notably, selecting a single ground motion to consider amplitude for seismic site response analysis may not be a reliable approach to estimate site amplification. The selection of a wide amplitude range and the assessment of likely fluctuation scenarios for Bhutan is not done yet. Hence, ground motion parameters that are related to the amplitude are investigated to examine and predict the variability, often regarded as sensitivity, concerning mean values and associated scatter.

In this paper, a sensitivity analysis of site response for specific soil conditions in Phuentsholing, Bhutan, is explored by a statistical correlation function of the ground motion parameters for different earthquake shaking intensities. The study area is one of the major urban and commercial hubs in Bhutan Himalaya, and seismic site effects on existing structures may have detrimental consequences due to inherent vulnerabilities of structures and infrastructures, as well as due to likely phenomena such as amplification in loose soil deposits. To quantify the seismic site effects in terms of amplification of amplitude parameters, a range of time histories is selected, and site response parameters are estimated.

Figure 1 Geology and seismicity and the study area: (a) geological map of Bhutan reproduced from McQuarrie et al. (2013) and seismicity, (b) location of Phuentsholing and geology of the area, and (c) study area showing surveyed site using MASW (modified from © Google Earth Pro 2021).

Figure 2 Typical borehole stratigraphy in Toorsa and Dhamdhara (Zone I) and Rinchending (Zone II).

2 Seismicity and geology of the study area

The Himalaya is one of the most seismically active regions on earth, which observes both large- and moderate-sized events frequently (Drukpa et al., 2006). Bhutan is located in the eastern Himalaya formed due to the subduction of the Indian Plate beneath the Eurasian Plate and spans from the low-lying Brahmaputra Plain to the high Tibetan Plateau. Most of the land area of Bhutan is underlain by the Main Himalayan Thrust (MHT), which runs along the entire length of the Himalayan arc. A historical earthquake catalog (see Fig. 1a) indicates that Bhutan has experienced several earthquakes of moment magnitude greater than 5 since the early 1900s, among them, the 1915 Trashigang (M_w=6.6), 1954 Trashiyangtse (M_w=6.4), and the 2009 Mongar (M_w=6.1) earthquakes are the most notable ones. The 2011 Sikkim–Nepal earthquake (M_w=6.9) also caused noticeable damage to building stocks in Bhutan (Chettri et al., 2021a). The earthquakes in the vicinity of the study area (Phuentsholing) include the 1981 Dagana (M_w=5.1) earthquake and the 2003 Haa earthquake (M_w=5.5). The most recent event occurred in Samtse in 2021 (M_w=5.1), affecting Phuentsholing and the neighboring areas with an intensity level of IV in the Modified Mercalli Intensity (MMI) scale (Gautam et al., 2022). The continuity of seismic activities in Bhutan is attributed to the presence of major shear zones such as the MHT, Main Boundary Thrust (MBT), Main Central Thrust (MCT), and the South Tibetan Detachment System (STDS) (Long and McQuarrie, 2010) as shown in Fig. 1a. The study area is within the Phuentsholing Formation of the Buxa group of the Lesser Himalaya, mainly characterized by highly weathered dark grey to black slate and phyllite and thin interbedding of limestone with a substantial amount of cream-colored dolomite and fine-medium quartzite, additionally consisting of fine- to medium-grained conglomeratic quartzite interbedded with phyllite and dolomite towards the Rinchending area of Zone II. Hence, the lithological characteristic of the area indicates weak and highly unstable geology in the region. The presence of thrust faults in the proximity of the study area along the entire belt of the Lesser Himalayan units and the quaternary sediments in the south depict the area to be seismically active with the majority of the historical earthquake events concentrated within these geological units. In particular, this study focuses on the city of Phuentsholing in the Chukha district in Bhutan (Fig. 1c). The city is one of the major commercial hubs for trade with India. The study area is undergoing rapid infrastructure development activities and urban expansion for residential, commercial, and industrial purposes. Phuentsholing city covers an area of 15.6 km² and is located at 89.39^∘ N and 26.86^∘ E. The city has a population of 27 658, mostly distributed towards the peripheral international border area with a total of 2263 residential and commercial buildings per the 2020 statistics (http://www.pcc.bt/index.php/, last access: 15 June 2021). The seismic site characterization includes eight locations in the regions of Dhamdhara, Toorsa, and Rinchending in Phuentsholing, Bhutan. In this study, the sites are grouped into two main zones based on the geographical location and immediate availability of survey locations. These two zones also refer to the local area plan (LAP) of Phuentsholing. The zones are Zone I: Dhamdhara I, Dhamdhara II, Toorsa I, and Toorsa II; and Zone II: College of Science and Technology (CST) football ground, CST hostel, Phajoding, and the Monastery area. Among the eight LAPs, Dhamdhara and Toorsa (Zone I) are in the same region in the western part of the city, and Rinchending (Zone II) is in the east.

3 Materials and method

3.1 Geotechnical site characterization

The geotechnical reports collected by Phuentsholing municipality have 29 stratigraphic logs. From these records, the depth of the water table (GWT) was demarcated first. Drilling log data showed the highest depth of the water table in the Dhamdhara area at 12.5 to 16 m, whereas the groundwater table in the Rinchending area is at 5 m, followed by the Toorsa area at 0.5 and 3 m, which is located near the riverbed. The depth of the water table is one of the essential input parameters used for one-dimensional ground response analysis. Three drill holes are presented to illustrate the typical underground stratigraphy (Fig. 2). Table 1 presents a summary of soil properties from laboratory testing of in situ samples collected from the drill holes. The number of samples in each zone represents the total number of samples collected from all drill logs at various stratigraphic depths. All laboratory tests have been verified according to the Indian standard codes. Testing included physical identification, Atterberg limits, grain size distribution, and direct shear testing. Field tests such as standard penetration resistance (SPT) and core cutter test were performed to determine resistance to penetration (SPT-N) and field density, respectively.

Table 1 Average soil parameters in the study area.

As shown in the stratigraphic logs, the upper stratum comprises predominantly mixed coarse-grained soils characterized by a considerable fraction of sand. The soil classification of the Phuentsholing area carried out by sieve analysis highlighted that most soils consist of 22.74 % gravel, 74.89 % sand, and 2.37 % silt and clay. The sieve analysis results for the respective zones are shown in Fig. 3. The soils in Toorsa are non-plastic, as coarse-grained soils dominate the particle distribution, while the soils in Rinchending and Dhamdhara have low plasticity with a plasticity index (PI) of 6.5 and 10, respectively. The bulk density is 1.8 g cm⁻³ in Toorsa, 1.64 g cm⁻³ in Dhamdhara, and 1.33 g cm⁻³ in Rinchending. The shear strength parameter, cohesion (c), ranges between 0 and 0.18 kg cm⁻², while the angle of internal friction (ϕ) in the study area is up to 35^∘.

Figure 3 Representative grain size distribution curve for the study area.

Shear wave velocity profiles from eight locations in the study area based on the multispectral surface wave analysis (MASW) and empirical correlation developed by Tempa et al. (2021a) are used for input parameters. According to the shear wave velocity profile, engineered bedrock (V_s>800 m s⁻¹) lies at a depth of 150 to 400 m as shown in Fig. 4. According to the parametric analysis carried out by Tempa et al. (2020), the site condition in the study area is classified as ground type B per Eurocode EC-08 and the National Earthquake Hazards Reduction Program (NEHRP) with the majority of shear velocity (V_s,30) values falling between 380 and 470 m s⁻¹, except for Phajoding, which has shear wave velocity of 584.76 m s⁻¹ (Table 2).

Table 2 Site classification as per Eurocode EC-08.

Figure 4 Shear wave velocity profile of study locations in Phuentsholing, Bhutan.

Dynamic properties of soils are influenced by shear modulus and damping and are defined by the respective degradation models, regarded as the backbone curves. Figure 5 represents the dynamic soil model for sand used in this study. Degradation models are well established by many investigators for different types of soils (see, e.g., Seed and Idriss, 1970; Vucetic and Dobry, 1991; Darendeli, 2001; Dobry and Vucetic, 1982; Seed et al., 1986). A damped linear elastic model of the soil system is used for the analysis. Due to soil nonlinearity for which the shear modulus is strain dependent, ProShake performs an iterative process on the linear model until both the moduli and damping ratios are compatible with the average strains and convergence is achieved at the last iteration (Shafiee et al., 2011; Puri et al., 2018). The nonlinear and hysteretic stress–strain behavior of soils under cyclic loading is approximated as a function of G_sec and G_max. The predetermined estimation of G_sec or G and G_max is attributed to unit weight or bulk density, ρ, and shear wave velocity, $V_{\mathrm{s}}(G_{\max }=\mathit{\rho }{V}_{\mathrm{s}}^{\mathrm{2}})$. Similarly, damping ratios are predicted as a function of G_sec or G values. This estimation is achieved using an iterative procedure in the ProShake 2.0 program (EduPro Civil Systems Inc., 2017).

Figure 5 Average modulus reduction ratio and damping ratio adopted for sand (Seed and Idriss, 1970).

Figure 6 Strong motions and corresponding Fourier amplitude plots of the input ground motions.

3.2 Selection of input motion

The definition of the input motion that is considered for site response analysis of an area requires both subsurface characterization and careful selection of acceleration time histories. In Bhutan, records of acceleration time histories are very rare, if not absent. In the absence of a national seismic code, Bhutan is assumed to fall under Indian seismic zones IV and V, with an expected maximum peak ground acceleration (PGA) of 0.24 and 0.36 g for design purposes. For these two zones, the PGA for earthquakes with a return period of 475 years is expected to be half of the maximum considered earthquake (MCE), i.e., 0.12 and 0.18 g. Notably, the Global Seismic Hazard Map (GSHAP) depicts the PGA range between 0.2 and 0.28 g with an increasing trend towards the east of the country. Considering the variations in expected PGA, we selected six acceleration time histories as input motions with PGA ranging from 0.067 to 0.422 g, considering the lowest and the highest range of possible earthquake scenarios (Table 3). The acceleration time histories used for the one-dimensional ground response analysis are shown in Fig. 6 in ascending PGA order using the ProShake 2.0 computer program. In the ProShake 2.0 program, input motion and soil profile are denoted as “M” and “P”, respectively, and are annotated in the subsequent sections (Table 3). The amplitude and frequency content of the bedrock level motion are particularly the most important parameters (Kirtas et al., 2015; Kramer, 1996). To understand the strong ground motion characteristics, we plotted the Fourier amplitude versus period in the frequency domain, representing the Fourier amplitude spectra of the input motions, as shown in Fig. 6. The effect of local soils is indicated at a much higher frequency range in all the investigated sites.

Table 3 Selected strong motion records for ground response analysis.

Figure 7 The typical profiles of normalized peak ground acceleration (PGA) at (a) Toorsa II in Zone I and (b) CST football ground in Zone II.

3.3 One-dimensional ground response analysis

One-dimensional-equivalent linear analysis is performed at eight sites in Phuentsholing, Bhutan, to estimate local site effects using the ProShake 2.0 program. In this study, six strong motion records are used to represent low, medium, and high acceleration categories. The ProShake 2.0 program provides the flexibility to input ground motions and soil profiles and is useful for estimating the outcrop responses to input ground shaking. The improved shear wave velocity profiles down to the engineered bedrock depth (150 and 400 m) from eight sites are used. The deep shear wave profiles used in this study incorporate the effects of depth and soil type of visco-elastic soil layers above the predicted engineering bedrock. The one-dimensional ground response analysis accounts for wave propagation from the bedrock outcrop through the visco-elastically stratified soil deposit and provides an estimate of the surface motion parameters. The complex response method is solved by the equation of motion in the frequency domain. Nonlinear soil response is estimated by an iterative quasi-linear procedure in which successive linear analyses are performed while updating the shear modulus and damping ratio based on the shear strain level obtained from the preceding iteration. Iterations continue until the strain-compatible modulus and damping converge.

Figure 8 Typical spectral acceleration of bedrock and ground surface motion at Toorsa II in Zone I corresponding to the respective input motions.

Figure 9 Typical spectral acceleration of bedrock and ground surface motion at CST football ground in Zone II corresponding to the respective input motions.

Figure 10 Examples of amplification factors for various earthquakes at (a) soil profile P1 at Toorsa II in Zone I and (b) soil profile P4 at Dhamdhara I in Zone I.

Figure 11 Examples of amplification factors for various earthquakes at (a) soil profile P1 at CST football ground in Zone II and (b) soil profile P3 at Phajoding in Zone II.

4 Results and discussion

4.1 Seismic site effects

Figure 7 shows normalized PGAs on the surface at two typical locations of the investigated zones. The chart shows PGAs of 1.2 to 1.5 g for low-PGA earthquakes and 0.7 to ∼1.1 g for medium- and high-PGA earthquakes. Response parameters can be defined and characterized based on the amplitude parameters of the ground motion and the severity of the ground motion excitation in nearby structures. This, in turn, is a function of the amplitude or intensity, the frequency content, and the duration of the ground motion (Bradley, 2011). Natural periods or frequency domain parameters are related to the seismic behavior of structures and indirectly reflect the ground motion characteristics (Zafarani et al., 2020). Hence, to commensurate this relationship, the response spectra of bedrock and surface motion are presented in Figs. 8 and 9, respectively. The results of various input ground motions indicate the higher spectral acceleration of the soil profile in the period range between 0.3 and 3.0 s, with the peak spectral acceleration range of 0.14 to 1.62 g. Thus, the structures with similar fundamental vibration periods are likely to be exposed to greater peak spectral acceleration.

Figure 12 Box and whisker plot for ground motion parameters of soil profile at P1 Toorsa II in Zone I.

Figure 13 Box and whisker plot for ground motion parameters of the soil profile at P1 CST football ground in Zone II.

Figure 14 Linear regression model for bedrock and surface spectral accelerations for Toorsa II (Zone I).

Figure 15 Linear regression model for bedrock and surface spectral accelerations for CST hostel (Zone II).

Figures 10 and 11 show the results of typical amplification factors at two locations in the study area. The amplification factors range from 0.7 to 2.7, 0.6 to 2.6, 0.75 to 2.5, and 0.7 to 3.2 for Toorsa II, Dhamdhara I, CST football ground, and Phajoding, respectively, for 0.01 to 0.1 s natural period. In the period range from 0.1 to 1.0 s, the amplification factors are in the range from 1.1 to 3.6, 0.7 to 4.2, 1.0 to 3.7, and 1.2 to 5.2 for Toorsa II, Dhamdhara I, CST football ground, and Phajoding, respectively. In the natural period range, the amplification factors are 5.0, 6.2, and 5.8 for Toorsa II, Dhamdhara I, and CST football ground, respectively. However, in the Phajoding the amplification factor is ∼1.7 due to a much stiffer soil deposit ($V_{\mathrm{s},\mathrm{30}}=\mathrm{584.76}$ m s⁻¹) and shallow engineering bedrock at 150 m.

4.2 Correlation analysis

The main objective of this study is to demonstrate the sensitivity of input motion amplitudes to predict the variability in seismic site effects due to local ground conditions. We examined the potential trends, patterns, and relationships between datasets for the numerical results. Using statistical analysis, variation in amplitude parameters is projected by box plots (Figs. 12 and 13). Statistical correlations are fitted between peak ground acceleration (PGA), peak ground velocity (PGV), peak ground displacement (PGD), and spectral acceleration (S_a) to determine the correlation between the effects of strong ground motion and the local soil conditions. As anticipated, the 1992 Petrolia earthquake with 0.422 g PGA (M_w=6.6) led to the greatest response. However, the 1994 Northridge earthquake with a PGA of 0.329 g (Mw=6.7) shows greater variability in spectral acceleration compared to other earthquakes. This is because the spectral acceleration corresponds to the interaction between the ground and the shaking intensity of an earthquake. Therefore, from the perspective of seismic site effects the box plot of the spectral acceleration (period or frequency domain) is highly scattered with the outliers, confirming uncertainty in the ground response characteristics in both regions. The El Centro and Petrolia earthquakes, with the highest PGAs, also appear to be closely associated with spectral acceleration.

Primarily, propagating energy waves (outcrop motion) act on each stratified soil layer that amplifies or de-amplifies the ground motion response parameters at each layer. The sensitivity of the input motion parameters is critically monitored, and enhanced correlations are developed. To outline this, a linear regression model for bedrock outcrop motion and the predicted motion parameters as a function of bedding depth is developed. Regression analysis is performed for one particular soil profile from two zones (Toorsa II and CST hostel) to substantiate the sensitivity analysis (Figs. 14 and 15).

The 95 % confidence interval (CI) shows a linear relationship for the Loma Prieta 2, Taft Kern County, and Northridge earthquakes, indicating a closer impact on surface motion that corresponds to the outcrop motion. In this case, the predominant frequency content of the input motion is between 1 and 10 Hz. In contrast, the Loma Prieta 1, El Centro, and Petrolia earthquakes, with a predominant frequency between 0.3 and 1.2 Hz, exhibit typical nonlinearity throughout the spectral range, indicating possible damping of the spectral responses of the soil deposits.

Figure 16 Sensitivity of input ground motion in Zone I. (a) Peak ground acceleration, (b) response spectrum intensity, (c) Arias intensity, and (d) mean frequency. Soil profiles: P1: Toorsa II; P2: Toorsa 1; P3: Dhamdhara II; and P4: Dhamdhara I.

Figure 17 Sensitivity of input ground motion in Zone II. (a) Peak ground acceleration, (b) response spectrum intensity, (c) Arias intensity, and (d) mean frequency. Soil profiles: P1: CST football ground; P2: CST hostel; P3: Phajoding; and P4: Monastery area.

Since all analysis sites are in type B sites, the trend of ground motion variation to surface is very similar, so the average values may be crucial for the better implementation of the scenario-based seismic risk in the study area. Ground response parameters such as the PGA and response spectrum intensity including the Arias intensity show linear variation for aggregated values while increasing intensity of earthquake shaking corresponding to a given soil profile. The mean, median, and standard deviation of the output parameters are computed. The response spectrum intensity is computed based on Housner's approach (Housner, 1959) as an integral from 0.1 to 2.5 s of the pseudo-velocity spectrum that provides an indication of the average velocity for most civil engineering structures. The plot of sensitivity of various input motions on amplitude parameters to different local soils for the two zones is shown in Figs. 16 and 17.

The standard deviation is lower for a set of predominant natural periods for a soil profile compared to the response spectrum dataset, and the deviation from the mean value indicates a stronger soil response to the single degree of freedom systems, as shown in Tables 4 and 5. Soil nonlinearity often shows a significant scatter in spectral acceleration at higher and lower periods, and therefore the practical reliability of the result is that it prompts more analysis with many input motions to predict the mean (or median) response with some level of confidence (Kramer et al., 2012). The sensitivity of input motion is shown in Figs. 14 and 15 from two investigated locations. The results of the correlation analysis and the sensitivity plots indicate that the input motion M4 (Northridge) has a significant influence on most of the response parameters. The additional ground response parameters are provided in Tables S1 and S2 in the Supplement.

Table 4 Descriptive statistics for averaged ground response parameters in Zone I for all four soil profiles and six input ground motions.

Table 5 Descriptive statistics for averaged ground motion parameters in Zone II for all four soil profiles and six input ground motions.

The PGAs of M4 (Northridge) are mapped to show the spatial variability in two zones as shown in Fig. 18. The PGA in Zone I is distributed between 0.37 and 0.42 g. The variability of PGA in Zone II is higher compared to Zone I as the PGA range for Zone II is 0.33 to 0.47 g. The resulting interplay of strong ground motion parameters with local soil conditions primarily highlights the importance of input motion characterization.

Figure 18 PGA distribution map of input motion M4 (Northridge) earthquake at (a) Toorsa and Dhamdhara in Zone I and (b) Rinchending in Zone II.

5 Conclusions

Using one-dimensional site response analysis, we performed sensitivity of various input motions. Ground motion parameter sensitivity for soft soil deposits is assessed considering a typical eastern Himalayan setting. Aiming to quantify the variation in input motion characteristics, we assessed several ground motion parameters. The conclusions of the study can be depicted as follows.

• The trend in the variation in ground motion parameters such as PGA, PGD, PGV, and SA projects an increasing order with ground motion intensity as expected. However, the ground motions with input PGAs greater than 0.34 g and less than 0.1 g are more sensitive than the others. The conclusion is that sensitivity is more prominent in low and high PGA ranges than the moderate shaking scenario (0.1–0.34 g).

• For loose soil sites characterized as type B ground, peak spectral acceleration is prominent between 0.3 and 3 s, and this implies that the structures with their fundamental vibration period between 0.3 and 3 s will observe greater peak spectral acceleration. The consideration of earthquake-resistant designs for the structures with a fundamental vibration period requires additional attention due to the severity in peak spectral acceleration occurrence.

• In general, the peak amplification factor is obtained up to 6.2 for the study area. The lower amplification factor coincides with the occurrence of bedrock early. Meanwhile, the soil columns with greater depth of loose soil deposits have reflected greater amplification. The spatial variation in amplification factor is quite significant even in a small area. Thus, more rigor is necessary for site response analysis and microzonation studies in soft soil deposits to incorporate the spatial variation in soil columns. If soil stiffness is increased, the amplification factor can be checked; thus, soil improvement may be required to assure foundation performance in loose soil deposits.

This study uses various strong motions to depict the variability in ground motion characteristics. Although this is one of the first studies in the area, the results are still preliminary, and detailed investigation using sophisticated soil characteristics and approaches could be effective in obtaining more reliable results.