Near-surface seismic refraction tomography and MASW for site characterization in Phuentsholing, Bhutan Himalaya

It is essential to understand the soil characteristics of the subsurface layers for any engineering construction. In difficult terrains like hilly areas, conventional methods of investigation are expensive and difficult to conduct. It calls for nondestructive testing methods to get reliable estimates of subsurface properties. In the present study, seismic refraction tomography (SRT) technique and multichannel analysis of surface waves (MASW) methods were carried out along five selected profiles in Phuentsholing region of Bhutan Himalaya. The profile length ranges from 37 to 81.5 m, and depth of imaging down to 10 m. While the SRT data imaged the P-wave velocity (Vp) structures, the MASW imaged the shear wave velocity (Vs) structures. The P-wave images provide a fair knowledge of geological layers, while the MASW images provide S-wave velocity structures (Vs). These results are useful to estimate soil parameters, like the density, Poisson’s ratio, Young’s modulus, shear modulus, N-value and the ultimate bearing capacity. The seismic images reveal the presence of sand, sandy clay, gravels and shale layers below the selected sites. Bhutan Himalayas being seismically vulnerable, the obtained results in terms of shear wave velocity were accustomed to categorize the sites as per NEHRP site classes, and a ground response analysis was performed to determine the reliable amplification factors. From the study, it is suggested that the engineering construction is feasible at all the sites except in one site, where an indication of saturated soil is observed which is vulnerable for liquefaction, and ground needs to be improved before construction at that site.


Introduction
Before commencement of any engineering construction, it is important to check soil profiles and soil characteristics in order to ensure that the site is acceptable for making such construction. In Bhutan Himalaya, the terrain poses difficulties in conducting conventional soil investigations using shell and auger, open-pit, etc. These methods are noneconomic and take time to perform. The conventional methods, viz. cone penetration test (CPT), provide required parameters of the soil but are costly and invasive, particularly in hilly terrains [16]. Seismic refraction tomography technique (SRT) as well as the multichannel analysis of surface waves (MASW), on the other hand, is time efficient and cost effective. Seismic refraction tomography (SRT) provides velocity structure and the required soil parameters without disturbing the natural ground condition [2,5,18,25]. Over the past years, seismic refraction methods have been used for civil engineering purposes in mountainous terrain [10,14,28]. They have been used to locate faults and in determining the soil parameters [4,16]. The SRT is a geophysical method which uses inversion technique to obtain 2D velocity depth profile [17,22,23] that images the cross-sectional picture along the profile by using the soils' response to the energy from the external source [6,7,11].
The multichannel analysis of surface waves (MASW) is a method used to estimate the shear wave velocity at a shallow depth and also help to determine the elastic parameters of the soil [14,30]. The SRT and MASW data are used for estimating seismic velocities and elastic parameters which are needed to estimate the bearing capacity of the soil mass and for the design of structures which is going to be constructed over it [15]. The damages due to seismic activity are controlled by local geology and soil parameters in addition to construction characteristics of structures. Seismic ground response at a particular location is strongly affected by in situ conditions of local geology and soil characteristics [19][20][21]. The soil condition has an important effect on the ground vibration. This gives the attention toward the importance of site classification in seismically vulnerable areas. Seismic site enactment due to shear wave velocity will indicate the mechanics of site amplification. The geophysical methods are seen as the most reliable for seismic site characterization [2,26]. This paper focuses on site characterization at five locations at Phuentsholing in Bhutan Himalayas using the SRT and MASW to assess the feasibility of engineering constructions and to classify the sites as per NEHRP site classification.

Study area
Bhutan is a small kingdom located in the eastern Himalayas, nestling between India and China with a total area coverage of 38,394 km 2 (Fig. 1). Our study areas are located in Phuentsholing Gewog under Chukha Dzongkhag covering an area of 15.6 km 2 with an average elevation of 293 m which serves as the administrative seat of the Dzongkhag [24]. Phuentsholing, about 172 km away from Thimphu (Capital of Bhutan), is one of the commercial centers of Bhutan and is a developing border town with India. Due to the importance of trade and commerce, the town has to expand with upcoming infrastructures [8]. In this study, five sites were selected in Phuentsholing area: Kabreytar 1, Kabreytar 2, Pipaldhara 1, Pipaldhara 2 and Phuentsholing town as shown in Fig. 1. Elevation of the sites varies from 210 m (Phuentsholing town) to 392 m (Kabreytar 1), and it is mentioned in Table 1 and shown in Fig. 2.
The sites were chosen based on three factors, such as easy accessibility, terrain and vegetation cover. All the sites were accessible and near to the road available. The sites Pipaldhara 1, Pipaldhara 2 and Phuentsholing town were on plane grounds, whereas Kabreytar 1 and Kabreytar 2 were on sloping grounds with the angle of elevation with horizontal road surface of 37.38° and 21.32°, respectively. Most of the sites have vegetation cover except Phuentsholing town.
The Phuentsholing area forms a part of the Himalayan foothill belt with a youthful topography steep topography and high precipitation. The Phuentsholing formation belongs to Baxa Group of Rocks consisting of quartzite, greenish gray, variegated and carbonaceous phyllite. The formation falls over the Main Boundary Thrust. The rocks in the area are crushed and weathered. The Himalaya being youngest mountain belt is still evolving due to northward movement of the Indian Plate toward the Eurasian plate. Hence, Himalaya is seismically active with large number of earthquakes. The region has weak intact rock strength having uniaxial compressive strength (UCS) value ranging from 0.60 to 1.50 MPa. There are also infrequent interbands of quartzite and dolomite that are strong rock with (UCS) values ranging from 34 to 55 MPa [29].

Field operation for SRT and MASW
The equipment for SRT and MASW comprises source, detector and a seismograph. A sledgehammer of 8 kg was used to produce the seismic waves. Twelve geophones with a frequency of 4.5 Hz were used. Seismograph of DAQ LINK 3 was used to record the signals received from the geophones. An alignment was chosen along which the geophones were inserted firmly into the ground at regular intervals. The seismic cables were used to connect the geophones and the seismograph. The sledgehammer was hit on a metal plate to produce the seismic waves (P-wave). These waves travel into the ground through refraction and reflection. Only the critically refracted waves are recorded by the seismograph. The sledgehammer was hit three times on the metal plates at regular intervals of the geophones and at the mid of two geophones. Taking a stack of three, the number of stacks collected for the sites Kabreytar 1, Kabreytar 2, Pipaldhara 1, Pipaldhara 2 and Phuentsholing town was 68, 27, 28, 33, and 30, respectively. Due to limited open plain area, the geophone interval for Kabreytar 1 was kept as 2 m and for the rest of the sites, it was 2.5 m. The lengths of the profiles obtained including the offset distances were 37 m, 44.5 m, 53.5 m 81.5 m, and 43.5 m for Kabreytar 1, Kabreytar 2, Pipaldhara 1, Pipaldhara 2, and Phuentsholing town, respectively.

Analysis of primary waves
A set of refracted P-arrival data are illustrated in Fig. 3. The SRT data were analyzed using the software SeisImager/2D. The waveform are displayed as a graph between time and distance. The first arrival time of P-wave was picked to generate the time-distance curve. The time-distance curve provides an initial velocity model (Fig. 4). The initial model is used for inversion; iteratively, the model is reduced with minimum root mean square (RMS) error for the observed and calculated travel times. Finally, the depth velocity model is presented in 2D (Fig. 5).

Analysis of secondary waves (MASW analysis)
For the analysis of shear waves, all the waveforms generated for a site have to be opened at once using Pick win, SeisImager/SW. Similar to the primary waves, the shot points and geophone interval should be edited for all the files opened. Next, the common midpoint (CMP) crossrelation gathers have to be generated (Fig. 6). CMP shows the receiver locations the shot point locations. The 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. The next step includes the generation of the dispersion curve using the WaveEq module. The dispersion curve is shown as a graph between phase velocity and the frequency. From the dispersion curve generated, 1D shear wave velocity profile (Fig. 7) was generated. The 1D shear wave velocity profile is shown as a plot between depth and velocity. The final step is to perform inversion analysis to get the 2D shear wave velocity, which is obtained using Geoplot (Fig. 8).

Results and discussion
Since no prior soil tests were done in the past and no information on the soil types was available for the sites, the velocity data were interpreted referring to Appendix B; Clause 3.2 of the IS: 1892-1979 [13]. Soil parameters such as density (ρ), Poisson's ratio (σ), Young's modulus (E), Shear modulus (μ), N-value and ultimate bearing capacity (q ult ) were calculated using equations as shown in Table 2. From the velocity of primary waves (V p ), the types of soil found in these sites are sand, sandy clay, gravels and shales which are presented in Table 3.

Ground response analysis
The results of the geophysical method indicate that all the five sites belong to site class 'D' (180 ≤ V s ≤ 360) as per NEHRP soil classification. In general, sand, silt and stiff clays with SPT value of 15 ≤ N ≤ 50 and shear strength of 50 kPa ≤ S ≤ 100 kPa belong to Site Class D.
The amplification factors proposed by previous researchers [22,27] can be directly used to assess the ground surface spectral acceleration if the input motion at bedrock is known. However, equivalent ground response analysis based on the subsurface profile at the site will give a more reliable estimate of site amplification value for the surface level and the variation of acceleration with varying depth. DEEPSOIL v7 [9] was used to perform equivalent ground response analysis using the shear wave velocity profiles and soil properties estimated from the geophysical methods. Figure 9 presents the typical subsurface soil profile The ground response analysis was carried out for the site considering a scenario earthquake with input motion from the Chi-Chi earthquake. The acceleration time history for the input motion of Chi-Chi earthquake is shown in Fig. 10.
The analysis in DEEPSOIL considers each and every layer of the soil profile, and within a layer, the properties like V s , 30 were varied for the site. It was seen that the acceleration value gets amplified as the waves travel to the top of the surface. The plot showing the variation of peak  Imai and Yoshimura [12] ground acceleration (PGA) with varying depth is shown in Fig. 11a for the site Kabreytar 1. The amplification factor was estimated as 1.36 by taking the ratio of surface PGA to the PGA at the bottom layer. The ground response analysis considering the varying soil profile gives a more accurate and reliable estimate of the amplification factor. The response spectrum showing spectral acceleration corresponding to various natural periods (of the structure) at the site Kabreytar 1 for an input motion from a scenario earthquake of Chi-Chi earthquake is shown in Fig. 11b. Similar analyses were carried out for other sites as well, and the results are presented from Figs. 12, 13, 14 and 15.
The surface level PGA values and estimated amplification factors for the considered sites in Bhutan Himalaya for input motion of Chi-Chi earthquake at the bottom layer are presented in Table 5. It was observed that amplification      9. Based on the weighted average shear wave velocity for top 30 m, the identified five sites were classified as site class D as per NEHRP site classification 10. Site-specific ground response analyses were performed for all the five sites based on shear wave velocity profiles and soil parameters at different depths. The amplification factors estimated for the sites Kabreytar 1, Kabreytar 2, Pipaldhara 1, Pipaldhara 2 and Phuentsholing town were 1.36, 1.14, 1.08, 1.11 and 1.07, respectively.