Parametric Study of Local Site Response for Bedrock Ground Motion to Earthquake in Phuentsholing, Bhutan

Abstract

Earthquakes, when it comes to natural calamities, are characteristically devastating and pose serious threats to buildings in urban areas. Out of multiple seismic regions in the Himalayas, Bhutan Himalaya is one that reigns prominent. Bhutan has seen several moderate-sized earthquakes in the past century and various recent works show that a major earthquake like the 2015 Nepal earthquake is impending. The southwestern city of Bhutan, Phuentsholing is one of the most populated regions in the country and the present study aims to explore the area using geophysical methods (Multispectral Analysis of Surface Waves (MASW)) for understanding possibilities pertaining to infrastructural development. The work involved a geophysical study on eight different sites in the study region which fall under the local area plan of Phuentsholing City. The geophysical study helps to discern shear wave velocity which indicates the soil profile of a region along with possible seismic hazard during an earthquake event, essential for understanding the withstanding power of the infrastructure foundation. The acquired shear wave velocity by MASW indicates visco-elastic soil profile down to a depth of 22.2 m, and it ranged from 350 to 600 m/s. A site response analysis to understand the correlation of bedrock rigidness to the corresponding depth was conducted using EERA (Equivalent-linear Earthquake Site Response Analysis) software. The amplification factors are presented for each site and maximum amplification factors are highlighted. These results have led to a clear indication of how the bedrock characteristics influence the surface ground motion parameters for the corresponding structure period. The results infer that the future constructional activity in the city should not be limited to two- to five-story buildings as per present practice. Apart from it, a parametric study was initiated to uncover whatever effects rigid bedrock has upon hazard parameters for various depths of soil profile up to 30 m, 40 m, 60 m, 80 m, 100 m, 120 m, 140 m, 160 m, 180 m and 200 m from the ground surface. The overriding purpose of doing said parametric study is centered upon helping the stack holders who can use the data for future development. Such a study is the first of its kind for the Bhutan region, which suffers from the unavailability of national seismic code, and this is a preliminary step towards achieving it.

1. Introduction

Although the Himalayan region is considered as highly seismic, the Bhutan region has not faced an extreme earthquake (M > 6.5) for the past six decades [1]. On the contrary, the surrounding regions around Bhutan have seen major earthquakes, like the northeastern part of India (M 6.7 in Sikkim, 2011 and M 6.7 in Manipur, 2016) and the disastrous Nepal earthquake (Gorkha) in 2015 with M 7.8. Although the Nepal region saw a major earthquake in 2015, such a disastrous event in Bhutan has not occurred in the last century owing to variations in the segments of the Himalayas. These segments of Himalaya have divergent topography, geology and earthquake patterns [2]. There is a widespread belief that major earthquakes would not strike in Bhutan; however, it was reported that a major earthquake greater than a magnitude of 8 occurred in May 1713 which left in its wake widespread infrastructural damage and a staggering loss of human life. Similarly, a major earthquake happened in 1100 and it is highly likely that other major seismic events may have occurred in the past [3]. The geophysical and geological aspects of the Bhutan Himalayas are quite different primarily due to the out-of-sequence fault of thrusts and its significance on distortion along the prominent Shillong Plateau [1].

One of the most popular methods to understand the possible hazards of seismic events on infrastructure development is by using geophysical methods. These methods have been used for several years by researchers and practitioners related to soils and foundations applications [4,5] and highway site investigation [6,7,8,9]. Geophysical methods not only provide means for probing the soil properties, sediments and rock outcrops but also to calculate the dynamic properties of soils like soil compression and shear wave velocities [5]. Similarly, the local site effects can be related to the conditions of surface and sub-surface geology and topography of a particular area along with the influence of an earthquake and its secondary damage to the engineered infrastructure [10].

The velocity of sound travelling through the soil depends on material composition and compaction. The seismic waves transmitted from a source at ground undergoes refraction at boundaries between different layers and return to the surface [5]. Evaluation of ground response is the most common problems encountered by geotechnical earthquake engineers [11,12,13]. Ground response analyses are crucial to substantiate accurate estimation of ground surface motions, which in turns aids in building design response spectra for proper evaluation of dynamic stresses and strains [14]. As mentioned above, the ground response usually depends on geological evolution and geotechnical characteristics of soil profiles. The various stiffness and damping characteristics of soil influence either elastic wave energy reflected and/or transmitted, with bedrock working all along as a fixed boundary [5,15,16,17,18]. In situ, geological and geotechnical parameters affect the ground shaking intensity and thus influence the extent of damage resultant of an earthquake [19]. It may also change the characteristics of surface seismic response, consequently impairing present structures constructed over such grounds.

The present study is a ground response analysis for bedrock characteristics using geophysical data and laboratory experiments on soil collected from the sites. Using the Equivalent-linear Earthquake Site Response Analysis (EERA) computer program as a tool for analysis, a layered visco-rigid soil profile is considered to analyze the effects of bedrock rigidness with changing depth. During geophysical site exploration using Multispectral Analysis of Surface Waves (MASW), it was concluded that the shear wave velocity measurements indicate soil profile down to a depth of 22.2 m. Next, a parametric study was initiated to study the effects of rigid bedrock on hazard parameters by considering the depth at 30 m, 40 m, 60 m, 80 m, 100 m, 120 m, 140 m, 160 m, 180 m and 200 m. For each analysis, it is assumed that there is a horizontal bedrock and a parallel soil layer due to its simplicity and practicality. In all the cases, the soil has been assumed as ground type B as per EURO Code EC8, 2003 as in the first case (up to 22.2 m), since due to topographical difficulties, borehole data could not be acquired. In the parametric study, the rigid bedrock with reference to EURO Code EC8, 2003 is accounted from 22.2 m to downward depth until the peak ground acceleration of varying depth of rigid bedrock has become approximately constant. The amplification factors are presented for each site and maximum amplification factors are highlighted. This parametric study has been carried out so that the same findings can be utilized for the other locations of the Chukha Dzongkhag where future development can be taken up by the Royal Government of Bhutan.

2. Study Area

Bhutan lies in the eastern end of the Himalayas. To the east, west and south it is bound by India while the Tibetan plateau stands to its North [20]. The nation is divided into eight subdivisions (known as “Dzongkhags” in Bhutanese) with an elevation ranging from 100 m to 7500 m [21]. The present study focusses on the Phuentsholing region situated in Chukha subdivision, which falls along the southwestern borders that Bhutan shares with India (Figure 1). The geographical coordinates, elevation and the profile directions of the sites collected for conducting the analysis are presented in

Table 1. Locations of the Sites (after NGI, 2009).

Site No	Site Name	Coordinates	Elevation	Remarks
1	Dhamdhara I	N 26 52 15.1 E 89 23 8.1	371 m	Profile NS Direction, 48 m
2	Dhamdhara II	N 26 52 29.8 E 89 23 01.4	309 m	Profile N 30 E Direction, 48 m
3,4	Toorsa riverbed (I & II) (Toorsa, Kabyetar)	N 26 52 33.1 E 89 22 45.0	192 m	Two Profile NS Direction, 96 m each
5	CST* hostel (Toribari)	N 26 50 57.8 E 89 23 42.1	439 m	Profile NE –SW Direction, 48 m
6	CST* football ground (Kharbandi)	N 26 51 01.3 E 89 23 35.8	423 m	Profile NS Direction, 96 m
7	Phajoding school (Pasakha)	N 26 51 4.7 E 89 23 50.1	406 m	Profile NS Direction, 48 m
8	Kharbandi monastry/Goenpa	N 26 50 57.0 E 89 23 22.7	378 m	Profile NS Direction, 48 m

Note: * College of Science and Technology (CST).

after Norwegian Geotechnical Institute (NGI), 2009. The Phuentsholing city is one of the main business hubs of Bhutan as Phuentsholing-Thimphu highway (or Asian Highway), is a primary road network of Bhutan which holds the key to bustling trade with the neighboring countries by road [22]. This has led to rapid urbanization and industrialization, playing a pivotal role in infrastructure development in this region.

The city has a population of around 27,000 people residing in 504 masonry, and 584 Reinforced Cement (RC) buildings (Standard and Quality Control Authority, 2009). Due to the ever-growing population, the shortage of housing is a major problem faced by the Phuentsholing City Corporation (PCC). To fulfill the housing problem, an additional urban land of 170 km² is required, which is more than four times the land currently occupied by human settlements (Norwegian Geotechnical Institute, 2009). To address such housing shortages, the PCC has come up with several Local Area Plans (LAP) since 2007 to extend its boundary in the east from Kharbandi and Toribari to Pasakha and in the west from Dhamdhara to Torsa, and in the north to Kabraytar (Geological Survey of Bhutan, 2010). The objective of such extension of the city involves the construction of buildings and facilities (public, commercial or private) and the PCC intends to develop all kinds of modern amenities in the LAP.

A general plot of earthquake data for the period 1937–2000 procured from the National Earthquake Information Center (NEIC) catalog covering the Bhutan region is depicted in Figure 2. As mentioned earlier, though Bhutan has not faced any major earthquake in recent times, the neighbouring regions of two countries India and Nepal have faced few considerable magnitude earthquakes. Apart from it, considering Indian Standard IS 1893:2002, Bhutan lies in between seismic zones IV and V. Again, Phuentsholing is a region which the Royal Government of Bhutan has identified for future development. The safety of such urbanization becomes one of the priorities to be addressed and implemented. One of the major challenges for this safety is understanding the effect of ground conditions on which urban development would begin [22]. Therefore, the present study on local site response analysis to an earthquake in Phuentsholing has become of utmost importance for a possible finding which can be utilized in a new practice for future infrastructure development in the region.

The geology and the geotechnical properties of a study region are significant as they enable to understand the spatial variation in ground shaking, amplitude, frequency and duration during the earthquake event. Similarly, lithology can also be a significant factor [23], therefore understanding the local geological sections (Figure 3) is significant for classifying area into different hazard zones. The study area lies in Phuentsholing formations of the Baxa Group, with formations established upon grounds of varying lithology and/or color [24]. The Phuentsholing formation is marked by substantially fractured and weathered phyllites, slates, and schists with an excessive amount of clay minerals [21]. It is comprised of slate- from dark gray to black, phyllite having thin interbeds of limestone as well as the cream colored dolomite with thin beds of dark-brown, fine-medium grained quartzite [25].

3. Methodology

The methodology adopted for the present study is depicted in Figure 4.

3.1. Geophysical Site Exploration

The study area is explored using a non-invasive MASW method to determine the shear wave velocity (V_s) profile of the region using wave energy. In this method, the dispersive nature of surface waves is manipulated to calculate the shear wave velocity in various layers, derived from the logic that shear wave velocity changes between soil layers due to change in medium properties. For the present study, multiple geophones and receivers were used to detect the shear wave velocity. The sledgehammer (source) was hit on a metal plate to produce the seismic waves. The source is offset at a predetermined distance from the near geophone and the Rayleigh wave velocity structure is obtained by utilizing the theoretically proven relation between Rayleigh and shear wave velocity as V_R = 0.9 V_S.

Table 2. Geophone array setup.

Site No.	No. of Shots	No. of Geophones Used	Spacing between Geophones (m)	Remarks
1, 2, 5, 6, 7 and 8	25	24	2	Restrictions in space
3 and 4	25	24	4	Flat surface, more spacing will give information of soil profile up to deeper depth

shows the geophone array setup in the present study. At all the study locations, the number of geophones used is 24, with a spacing between the geophones as 2 m except locations at sites 3 and 4, where spacing is kept at 4 m. The justification for the spacing of geophones is mentioned in the remark column of Table 2.

A typical time-distance relationship is shown in Figure 5 for ready reference, where the source was at 0 m. To determine shear waves, all the raw data acquired or waveforms generated after the shots for a site have to be opened at once using SeisImager software. Next, the shot points and geophone interval should be edited for all the files opened. The collected seismograms were processed to obtain the surface wave dispersion curves and then common midpoint (CMP) cross-relation points were generated where yellow dots are the receiver locations, i.e., the position of the geophones and the blue dots are the locations at which the hammer was hit. A typical step is shown in Figure 6. The next step includes the generation of a dispersion curve which displayed the graph of phase velocity to frequency. From the dispersion curve, 1-D shear wave velocity is generated as shown in Figure 7.

The shear wave velocity (V_s), arrived at by following the steps as mentioned above for all the present study locations, is considered for site classification in the present study. The international codes—Eurocode-8 (2003), National Earthquake Hazards Reduction Program (NEHRP) (BSSC, 2003) and International Building Code (IBC, 2009), stipulate site classification established on the average V_s in the top 30 m soil column (V_S,30). However, it is mandatory to consider near-surface Vs of soil when coming into contact with rock within a depth of about 30 m. Otherwise, the obtained value of V_s,30 will be greater for the velocity of a rock mass [26]. For the present study, a mean value of Vs down to 30 m (V_S,30), as described by NEHRP, is considered grounds for site classification referring to

Table 3. Mean shear wave velocity and soil parameters of the study area.

Site No	Site Name	VS,30 (m/s)	Ground Type	N Value	qult (kN/m2)
1	Dhamdhara I	386.43	B	38	1017.6
2	Dhamdhara II	435.92	B	50	1358.9
3	Toorsa riverbed I	439.54	B	51	1386.1
4	Toorsa riverbed II	464.30	B	57	1581.0
5	CST Hotel	426.76	B	48	1291.3
6	CST Football Ground	426.61	B	47	1290.3
7	Phajoding School	584.76	B	96	2750.1
8	Monastry/Goenpa	446.20	B	53	1437.1

, EURO Code EC8. The mean value of shear wave velocity for the top 30 m, V_S,30 is computed by using Equation (1) as follows:1

$$V_{S,30}=\frac{30}{{{\displaystyle \sum }}_{i=1,N}()\frac{h_i}{V_{S,i}}} [\mathrm{m}/\mathrm{s}],$$

(1)

where h_i and V_S,i represent the thickness and shear wave velocity of the i-th formation or layer, in a total of N, existing in the top 30 meters. Table 3 shows V_S,30 of different sites on which the local seismic response study was conducted. To calculate the N value and ultimate bearing capacity, q_ult, Equations (2) and (3) below are used, proposed respectively by [4,27], with values listed in Table 2.

$$\mathrm{N}={\left(\frac{\mathrm{Vs}}{76.55}\right)}^{2.24719}$$

(2)

$${\mathrm{q}}_{\mathrm{ult}}=\frac{{\mathrm{Vs}}^{2.4}}{1590}[\mathrm{kN}/{\mathrm{m}}^2]$$

(3)

3.2. Ground Response Analysis

The ground response analysis aims to ascertain the effects of soil conditions of the particular region on the intensification of seismic waves, thus estimating the ground response spectra for the sake of design purposes. Taking into account the range of approaches utilized in the modeling and analyzing of site effects, the 1-D ground response analysis method is especially flexible and uncomplicated in its application to estimate site amplification under earthquake shaking [19,28]. One of the most popular techniques to perform 1-D ground response is the use of the Equivalent-linear Earthquake Response Analysis (EERA) program developed by [29]. Input data required to run the program are full time history of an earthquake, soil profile and dynamic soil properties. Damping ratio and shear modulus degradation curves influence the dynamic soil properties, and are a representation of equivalent damping ratio regarding secant shear modulus and strain [29].

The shear modulus of the soils is a property bearing immense effect on the evaluation of a suitable response spectra and study of shear stresses. In equivalent linear approximation of nonlinear response, the shear modulus and damping ratio change with strain and this takes place according to equivalent linear soil properties. For this reason, a known time history bedrock motion is considered for fast Fourier transformation (FFT) of Fourier series. It is required to determine the transfer function with respect to each layer using presently available properties of a soil profile. Non-linearity of the soil conduct is notable, which is why it is important for a linear methodology to adapt in order to produce sensible evaluations of ground response concerning practical issues of intrigue. Propagation of ground motion is determined by initially estimating the shear modulus (G) and damping ratio for available ground motion time history. The equivalent linear soil properties can be approximated from nonlinear hysteresis of the stress-strain behavior of cyclically loaded soils [11]. The equivalent linear shear modulus, G, is in general considered as secant shear modulus, the equivalent damping ratio, ξ, as damping ratio which gives rise to the same energy loss in a single cycle in an actual hysteresis loop of irreversible soil behavior (Figure 8).

As shown in Figure 8a, from Equation (4), G_s at the ends of symmetric strain-controlled cycles is:

$${\mathrm{G}}_{\mathrm{s}}=\frac{{\mathsf{\tau }}_{\mathrm{c}}}{{\mathsf{\gamma }}_{\mathrm{c}}}$$

(4)

where τ_c and γ_c are the shear stress and strain amplitudes, respectively. In site response analysis, the material behavior is usually represented as depicted in Figure 8b. G_s-γ curves do not appear as arbitrary shapes, being derived from a stress-strain curve (Figure 8a). Figure 8b illustrates G_sec and damping to strain curves, which shows the uncertainty in degradation properties for a specimen site. Often average G_s and damping curves should be acquired from a few laboratory tests or from generic results, and uncertainties must be adopted to reflect unknown conditions. The degradation curves have a negative correlation and this means a higher than average G_s curve is likely associated with a lower than average damping curve [14].

4. Data Used

4.1. Input Motion and Bedrock

The Equivalent-linear Earthquake Response Analyses (EERA) program allows a ground response analysis by examining the attributes of a soil or rock layer exposed to an acceleration time history. For convenience, the Los Angeles University station’s record of the acceleration time history of the 1994 Northridge earthquake (M_w = 6.7) is the one taken up for the current ground response study (Figure 9).

Findings derived from the record indicate peak ground acceleration (PGA) of ground (input) motion as 0.235 g, approximately. The recommended value of PGA from the Global Seismic Hazard map for Phuentsholing city (GSHAP) was 0.20–0.28 g. Therefore, for the present study, the response spectra are considered as per illustrated in Figure 10a,b. Currently, Bhutan does not have any seismic zonation map and therefore, the seismic zonation map of India (IS 1893: 2002) is being followed for making any building earthquake resistant. As mentioned above, the MASW analysis in the present study indicates shear wave profile down to 22.2 m, PGA is constant at this level and depth of the bedrock is accounted for at 23 m. According to Euro Code EC8 (2003), a rock or rock-like formation describes ground type A with V_s,30 > 800 m/s. So, for simplicity, in the present case, the average shear wave velocity of the bedrock is considered as 1000 m/s with 0.5% damping.

4.2. Soil Properties

In order to execute a thorough seismic response analysis specific to the study site, sub soil information is one of the most essential data. The geological parameters of a layer up to a depth of the first 30 m are significant in amplification for earthquake shaking [30,31].

4.2.1. Static Properties

Based on sieve analysis and a shear strength parameter test performed in the laboratory on the soils collected from the sites, it is clear that in Toribari, Kharbandi and Pasakha areas, the gravel and sand contents are 25–50% and 25–55%, respectively and c = 8–22 kPa and ϕ = 19°–28°. Similarly, for the Dhamdhara and Torsa areas, the gravels (> 2 mm size) are more than 50% and medium-fine coarse sand to coarse gravels with less fine contents are < 5% with cohesion intercept c = 10–24 kPa and ϕ = 18°–26°.

4.2.2. Dynamic Properties

Soil characteristics that influence the propagation of wave and other low-strain phenomena include density, stiffness, damping and Poisson’s ratio. Out of the mentioned criteria, stiffness and damping stand to be predominantly important; the other two, in contrast, have little influence and belong within relatively narrow ranges [11]. With respect to the equivalent linear 1-D analysis in EERA, the dynamic soil properties denote damping ratio to the shear modulus degradation curves. Based upon extensive experimental studies, curves for cohesionless soil (sand) and cohesive soil (clay) were put forward by [32,33]. One for gravel was submitted by [34]. For the current study, shear modulus degradation and damping curves used are in accordance with [32] for sand since soil parameters of the specific site concerned satisfy ground type B as shown in Figure 11a,b, respectively. Three curves of sand are represented—lower bound, upper bound and average values, showing that lower bound indicates a less stiff with less damping ratio in comparison to the upper bound at the same shear strain. For both cases, in the present study, average value curves are used. The average values are considered by keeping in mind that the soil available in the sites are coarse-grained soils (Table 2) with plasticity index PI = 0 [35].

5. Results and Discussion

5.1. Geophysical Site Exploration

For the present study, the V_s profiles of the sites are interpreted based on the 2-D V_s profile images. Three to four shear wave profiles are reported for each site. Each shear wave profile is interpreted, and the average interpreted V_s profile is used for 1-D response analysis (Figure 12).

In the current study, average shear wave velocity over the upper 30 m (V_S,30) is arrived at as per Equation (1), fitted to the definition of the NEHRP. The shear wave velocity for sites 1, 2, 3, 4, 5, 6, 8 indicates that the soil in all the sites broadly exhibited similar properties with the majority of shear velocity ranging between 380–470 m/s (Figure 12 and Table 2), except site no. 7 which has a shear velocity of 584.76 m/s. As we know, the properties of each layer of soil beneath the ground surface may vary, due to which a major influence on the dissemination of ground motions from the bedrock to the ground surface may become evident. In the case of Phajoding School, the type of soil may be of very dense sand or gravel than the other locations where the soil profile may consist of dense or medium-dense sand. However, as a whole, considering all study locations, the value of Vs ranges from 360–800 m/s, which satisfies the ground type B as per Euro Code EC08.

Again from Table 3.1, Euro Code EC08, in the parameters column, it is mentioned that to define the ground type, along with Vs, the standard penetration test (SPT) N-value will help to reinforce the claim of soil type. So the N-value is determined based on Equation (3), proposed by [4]. It is seen that apart from site Dhamdhara I, the rest of the N-values are either >50 or nearer to 50 (Table 2) which almost satisfies the above claim of ground type B listed in Euro Code EC08.

The bearing capacity of a soil determines whether the soil is competent for any building construction or not. Therefore, in the present study, the ultimate bearing capacity (q_ult) of soil is calculated for future reference using Equation (3), proposed by Imai and Yoshimura (1976) and listed in Table 2. From Equation (3), as q_ult is directly proportional to V_s, so among all the locations, the bearing capacity of Phajoding school is the highest (2750.1 kN/m²), which is almost more than double the average ultimate bearing capacity of the rest (1337.47 kN/m²).

5.2. Ground Response Analysis

As discussed above, the 1-D ground response analysis was carried out using EERA, referring to the historic 1994 Northridge earthquake as the input ground motion. In executing this particular analysis, shear modulus and damping curves were generated in accordance with [32]. A sequence of input parameters and output results for Dhamdhara I is graphically listed in Figure 13. The acceleration time history of the 1994 Northridge earthquake is illustrated in Figure 9. As an input parameter, the PGA 0.235 g extracted as recorded at station Los Angeles University on the day of the event. Figure 13a depicts variations in shear modulus (G_max) with relation to depth, thereby establishing that shear modulus increases with depth, which indicates that the soil becomes stiffer and can bear more load for construction. The value of G_max is extracted from Figure 11 as suggested by [32]. The value of G_max when it comes to Dhamdhara I ranges between 117–304 MPa, from ground surface to bedrock level at the depth of 22 m.

Figure 13b illustrates the variation of shear wave velocity (V_s) in relation to depth, which is used as one of the input parameters in EERA. Here too, it is easy to observe that the value of V_s ranges from 233–360 m/s starting from the ground surface. One should note that the present Vs is showing the variation with respect to changes of soil profile due to an increase in depth, and it is different from average V_s in the top 30 m soil column (V_S,30) calculated using Equation (1) as defined by NEHRP and which is used for ground type determination. Yet another input parameter applicable to ground response analysis is the unit weight of the soil, adopted to be 21.08 kN/m³ (Figure 13c).

Figure 14 presents the output characteristics of Dhamdhara I as given by EERA. Figure 14a illustrates the variation of maximum acceleration with reference to depth, showing how ground acceleration at the bedrock level is made stronger by properties of sediment available in between bedrock and surface level. Ground acceleration at bedrock rises, growing from 0.081 g to 0.148 g at surface level. Figure 14b shows the amplification spectrum between the bedrock and ground surface. For Dhamdhara I, the maximum amplification factor was found to be 1.38 for the frequency 3.8 Hz. Figure 15 represents a time history of spectral response spectrum of surface level and maximum spectral acceleration (g) is 0.469 g with 5% critical damping.

Table 4. Seismic response output characteristics from Equivalent-linear Earthquake Site Response Analysis (EERA) for all study locations.

Site No.	Sites	PGA at Surface Level (g)	Maximum Strain (%)	Maximum Stress (kPa)	Amplication Factor	Frequency at Amplication (Hz)	Maximum Spectral Acceleration (g)
1	Dhamdhara I	0.158	0.0014	1.553	1.38	3.8	0.469
2	Dhamdhara II	0.181	0.0016	1.866	1.50	10.0	0.596
3	Toorsa riverbed I	0.172	0.0012	1.734	1.38	4.0	0.534
4	Toorsa riverbed II	0.171	0.0010	1.735	1.35	4.8	0.525
5	CST hostel	0.159	0.0010	1.620	2.46	2.8	0.562
6	CST football ground	0.126	0.0006	1.200	1.26	1.8	0.366
7	Phajoding school	0.212	0.0013	2.211	1.74	6.8	0.666
8	Monastry/Goenpa	0.128	0.0007	1.220	1.06	1.8	0.388

supplies a number of seismic response characteristic figures from EERA for all the study locations. It is seen that PGA for all the locations ranges from 0.126–0.212 g and the acceleration at a bedrock level for all the locations ranges from 0.081–0.110 g. High ground acceleration at the surface level is detected in Phajoding school site while the opposite is found in Dhamdhara I site. It is because the thickness and characteristics of soil profile formation are essential criteria above the bedrock. They influence the frequency band of ground motion and might be significantly amplified by in situ soil conditions. From the analysis, it is seen that the time required to reach PGA ranges from 3.94–3.97 s with Phajoding school taking the shortest time while, Dhamdhara I, CST hostel and CST football ground taking the longest time.

Table 4 also shows the maximum level of strain in the upper layer during seismic events which varies from 0.06–0.16% with CST football ground having the least strain, and Dhamdhara I the highest. Apart from it, the maximum pressure propagated by seismic events will be generated in the range of 1.200 kPa (CST football ground) to 2.211 kPa (Phajoding School). Outcomes indicate seismic pressure on the soil is tied more closely to the depth of the soil than its elasticity. The thicker the soil layer with different properties, the more the pressure and strain developed. Fourier spectra demonstrate that the basic frequency due to seismic events will fall within the range of 5.52 Hz to 10.08 Hz. Taking all the study locations into consideration, the highest amplification factor is determined to be 2.46 in Goenpa with frequency at an amplification of 2.8 Hz, whereas the lowest 1.06 is found at CST hostel with a frequency of 1.8 Hz. The frequency at an amplification of 10.0 Hz is the highest in Dhamdhara II.

In the present study, a response spectrum analysis is carried out with a critical damping ratio of 5%. Study outcomes point to a bigger spectral acceleration on the surface, 0.67 g with a period of 0.09 s (Phajoding school) and the smaller 0.37 g with a period of 0.07 s (CST football ground). From the analysis, it is seen that all the locations have different response spectral acceleration profiles for the different characteristic of a soil profile, and it is from medium-dense sand to very dense sand or gravel (type B, Euro Code EC089).

5.3. Parametric Study

In this part, a parametric study was performed for determining the effects of bedrock on hazard parameters. For this, the depth of bedrock is considered only after 30 m, 40 m, 60 m, 80 m, 100 m, 120 m, 140 m, 160 m, 180 m and 200 m of a soil layer. For simplicity and practicality of analysis, the assumptions considered are—horizontal bedrock, parallel soil layer, ground type B as per EURO Code EC8, 2003, and for the bedrock, the unit weight is considered as 25 kN/m³. The result of the site-specific response analysis for all the investigated sites shows that the PGA of the ground surface motion decreases as the depth of bedrock increases (Figure 16a,b). Figure 16a illustrates the fluctuations in PGA with respect to depth of bedrock for the Monastry, CST hostel and football ground and Phajoding school. It can be noted that the value of PGA ranges from 0.182–0.201 g at the depth of 30 m for Monastry, CST hostel and football ground. An exception observation was noted for Phajoding School showing a PGA value of 0.263 g at the same depth of bedrock as the former three. The PGA at a specific depth of bedrock can be predicted more appropriately through the exponential equation provided in the same figure for the above four locations. Since three sites—Monastry, CST hostel and football ground, are showing identical variations of PGA with a depth of bedrock, only one exponential equation is proposed for these three sites taking the average of PGAs with respect to depths. An exception observation of PGA was noted for Phajoding school, so a separate exponential equation is proposed. Similarly, Figure 16b depicts fluctuations in PGA with respect to the depth of bedrock in the sites Dhamdhara I and II, and Toorsa I and II. From the figure it can be noted that the PGA value ranges from 0.191–0.221 g at the depth of 30 m for all the above sites of Dhamdhara and Toorsa. Here too, since all the four sites are identical variations of PGA with respect to depth of bedrock, only one exponential equation is proposed taking the average of PGAs with depths to predict the PGA at a specific depth of bedrock.

Figure 17a–d shows the change of amplification factor with the depth of bedrock for 5% critical damping for CST hostel, CST football ground, Phajoding school and Monastery, respectively. Similarly, Figure 18a–d shows the change of amplification factor with the depth of bedrock for 5% critical damping for Dhamdhara I and II, and Toorsa I and II, respectively. The site amplification was computed as the ratio of the response spectrum of the surface motion divided by that of the bedrock. From given figures, the response analysis due to local ground conditions with varying depth of rigid bedrock shows a significant amplification on surface ground response parameters. At bedrock depth of 30 m to 60 m, amplification factors are observed to be as high as 2.81 between 0.35 s to 0.72 s which generally holds good for two- to eight-story buildings based on Indian standard IS 1893:2002. As of now, the competent body of the Royal Government of Bhutan allows construction of up to five-story buildings in Phuentsholing city. So, from the present findings, it can be proposed that for future building constructional activity in the city, one should not be limited to two- to five-story buildings as per present practice.

6. Conclusions

In the present study, eight different locations were selected in Phuentsholing to explore the area using geophysical methods (MASW) for understanding the possibilities pertaining to future infrastructural development. A site response analysis to understand the correlation of bedrock rigidness to the corresponding depth was conducted using EERA software. The amplification factors are determined to signify the maximum amplification factors. The concluding remarks are listed as below:

• The present study is a parametric ground response study for bedrock characteristics using the geophysical data and laboratory experiments on soil collected from the sites.
• The N-value and ultimate bearing capacity values are calculated on the basis of the shear wave velocity of a soil profile.
• Average shear wave velocity indicates soil in each site broadly exhibits similar properties with a majority of V_s,30 ranging between 360–800 m/s which according to the Euro Code EC08, 2003 falls in soil type B, which means the soil is very dense sand or gravel.
• In the present study, no amplification due to the site-specific response analysis concerning the specific site in the study shows that PGA near ground surface ranges from 0.191 to 0.221 g for all the study locations at 30 m depth of bedrock.
• At bedrock depth of 30 m with amplification factor 2.81 between 0.35 to 0.38 s generally holds good for two- to eight-story buildings based on Indian standard IS 1893:2002.
• Topographical features are considered. So, during the construction of building in the study area, the designer should consider other cataclysmic events which are regular features in Bhutan due to its orography, like landslides, slope instability near the slope edges, etc. This section is not mandatory, but may be added if there are patents resulting from the work reported in this manuscript.