Examination of shallow and deep S-wave velocity structures from microtremor array measurements and receiver function analysis at strong-motion stations in Kathmandu basin, Nepal

The Himalayan collision zone, where the Indian Plate subducts beneath the Eurasian Plate at a low angle, has caused many devastating earthquakes. The Kathmandu basin, situated in this region, is surrounded by mountains on all sides and is filled with distinct soft lake sediments with a highly undulating bedrock topography. The basin has been experiencing rapid urbanization, and the growing population in its major cities has increased the vulnerability to seismic risk during future earthquakes. Several strong-motion stations have recently been deployed in the Kathmandu basin. It is expected that the data captured by this strong-motion station array will further enhance our understanding of site amplification in sedimentary basins. Clear P-to-S converted waves have been observed in the strong-motion records. In this study, we investigate the medium boundary that generated these converted waves. First, we estimate the shallow velocity structures, which correspond to the topographic slopes or surface geology, beneath the strong-motion stations. We then apply a receiver function analysis to the strong-motion records. The receiver function indicates that the interface between the soft sediment and seismic bedrock serves as a boundary that generates converted waves. The obtained results can be used for tuning three-dimensional velocity structures.


Introduction
The Kathmandu basin is an oval intermontane basin located in the Lesser Himalayas.It is approximately 30 km long in the east-west direction and approximately 25 km wide in the north-south direction and surrounded by Shivapuri Mountain to the north and Chandragiri range to the south.The average elevation is approximately 1340 m and the Bagmati River flows through the basin (Dhital 2015).Paleo-Kathmandu Lake disappeared from the basin around twelve thousand years ago due to the draining of the lake water (Sakai et al. 2016).The maximum sediment thickness at the central part of the basin is 650 m, as determined from gravity measurements.Exposed bedrock with local gravity highs can be seen at different sites within the basin (Moribayashi and Maruo 1980).The basin floor is overlain by unconsolidated sediments of Plio-Pleistocene origin.There are fluvio-deltaic facies in the north, lacustrine facies in the center, and alluvial fan facies in the south (Fujii and Sakai 2002).Kalimati clays, meaning black clay of lacustrine origin with rich organic matter are distributed with a thickness of 200 m in the central part (Katel et al. 1996).
The Himalayan collision zone, where the Indian Plate subducts beneath the Eurasian Plate at a low angle, has caused devastating earthquakes and damage to the basin (e.g., in 1833, 1934, and 2015) (Bilham et al. 2001;Chaulagain et al. 2018;Zilio et al. 2021).The population within the Kathmandu basin, consisting of Kathmandu, Lalitpur, and Bhaktapur districts, is 3.03 million in 2021 (National Statistics Office 2021).Rapid urbanization and population growth in the major cities within the basin have also increased the vulnerability to seismic risk during future earthquakes (Mesta et al. 2022).The most recent M w 7.8 earthquake on 25 April 2015 initiated around 80 km west of the Kathmandu basin in the Barpak region of Gorkha District.The large slip area was located in the north of the Kathmandu basin, and the fault rupture propagated toward the east (Kobayashi et al. 2016), causing severe damage in central Nepal.Although a few seismic stations acquired strong-motion records, long-period ground motions with a period of 2-5 s were observed in the basin (Bhattarai et al. 2015;Takai et al. 2016).The building damage around the strong-motion stations was less than anticipated for an earthquake of this magnitude.Severe building damage was concentrated in a few areas, such as Bhaktapur and Balaju in the basin (Bijukchhen et al. 2017a).Several soil liquefaction events were reported, indicating that the Kalimati Formation is susceptible to liquefaction (Gautam et al. 2017).The sedimentary layers enhance the amplification of the strong ground motion.An evaluation of the site amplification characteristics based on the strong-motion records is essential for assessing seismic hazards in the basin.
Meanwhile, the shallow S-wave velocity structure beneath the stations was not sufficiently clear, we conducted small-sized microtremor array measurements to investigate the S-wave velocity structure up to approximately 30 m beneath the stations.Besides, P-to-S converted (Ps) waves clearly appear in the accumulated strong-motion records for the Kathmandu basin.Ps waves are commonly used to estimate the interface depth of velocity structures.In this study, we apply receiver function analysis to the records acquired by the strongmotion station, including the stations installed in the period 2016-2018, to investigate the generation of converted waves.Finally, we compared the Ps time obtained through the receiver function analysis with the corresponding time derived from the shallow to deep velocity structure.

Strong-motion network
Four permanent strong-motion stations with three components, installed in collaboration with Hokkaido University, Japan, and Tribhuvan University, Nepal, located along the east-to-west profile (Takai et al. 2016), have been collecting strong-motion data since September 2011.The KTP station is located at an exposed rock site in the basin.Strong-motion records from the 2015 Gorkha Nepal earthquake were acquired at these stations.For three months after this earthquake, four temporary stations, located along the north-to-south profile and in heavily damaged areas, were added for observing aftershock activity (Ichiyanagi et al. 2016;Shigefuji et al. 2022).Strong-motion records of aftershock sequences of the 2015 Gorkha earthquake were obtained at these eight stations.The one-dimensional velocity structures beneath these stations were discussed by strong-motion records of moderate-sized earthquakes (Bijukchhen et al. 2017b).However, the shallow parts of these structures have not been well examined.In addition, ten strongmotion seismometers were installed from November 2016 to May 2018 at the rock sites and low-density observation areas on the sediments (Takai et al. 2021) under the Science and Technology Research Partnership for Sustainable Development (SATREPS) project.These stations are managed by the Department of Mines and Geology (DMG), Government of Nepal. Figure 1 and Table 1 show information on these stations.Figure 1c shows the location of the stations on engineering environmental geological map from the DMG (Shrestha et al. 1998).

Shallow S-wave velocity structures
Microtremor array measurements Takai et al. (2015) and Shigefuji et al. (2018) carried out a multi-channel analysis of surface waves (MASW) close to the permanent strong-motion stations to estimate the S-wave velocity.In both studies, 24-channel vertical geophones with a natural frequency of 4.5 Hz (Geospace GS-11D) were laid out in a linear array, with a spacing of 1 or 2 m.This configuration provided satisfactory results to a depth of approximately 10 m.To extend the velocity structure to a depth of 30 m, small-sized microtremor array observations were carried out close to the 12 strong-motion stations on 17-20 September 2022, 5 November 2022, 3 May 2023 (temporary stations, and the KPN and SNK stations managed by the DMG were excluded).The microtremor array technique is the most frequently used method to estimate the sedimentary S-wave velocity structure (e.g., Karagoz et al. 2015;Yokoi et al. 2018;Takai et al. 2019;Chimoto et al. 2023).
In the present study, instruments for the microtremor array included three-channel geophones and a 24-bit A/D data logger (McSEIS-AT, OYO Corporation).The sensors were arranged at the vertices and center of an equilateral triangle.The side length of the microtremor array ranged from 1 to 46 m depending on the available space at a given station.The observation information, including for the MASW, is listed in Table 2.The sampling frequency was 250 Hz.Continuous records were acquired for more than 25 min for each array with a side length of 40 m, more than 10 min for each array with a side length of 4-30 m, and more than 4 min for each array with a side length of 1-3 m.The phase velocity of the Rayleigh waves was calculated by applying the Spatial Autocorrelation (SPAC) method (Okada 2003) to the  microtremor records.Each record was divided into segments of 8.92-65.536s, band-pass filtered at 0.1-100 Hz, and smoothed using a Parzen window with a bandwidth of 0.2-2.0Hz, depending on the array size.
Figure 2 shows the estimated fundamental Rayleigh wave phase velocity dispersion curves.The phase velocities calculated using SPAC method were combined with those obtained from the MASW at the corresponding site.The phase velocities in the overlapping frequency range almost corresponded and were continuous.Incidentally, the phase velocities of SNK and KPN are exclusively derived only from the MASW.The sites were categorized into three groups in Fig. 2 based on their locations: first one located on the bedrock, second and third located on the sediments either above or below an average elevation of 1340 m in the Kathmandu basin floor (Table 1).The average elevation corresponds to the upper boundary with the younger deposits in the lower level of the basin, which are the Thimi and Patan terraces (Sakai et al. 2006;Gautam et al. 2009).The sites below the average elevation were mainly located in the lacustrine Kalimati Formation (Shrestha et al. 1998).The phase velocities for the rock sites were steep in the frequency range of 5-80 Hz.The values for the SNG and JHR stations, located outside the basin, were more than 400 m/s at 5 Hz.The value for the KTP station, located on an isolated bedrock, was 700 m/s at 20 Hz.The phase velocity at the sites above the average elevation was obtained in the frequency range of approximately 1-30 Hz.It was found to be 200-400 m/s at 5 Hz.At the sites below the average elevation, low phase velocities of 70-200 m/s were found over a wide frequency range of approximately 2-30 Hz.The THM station had the lowest phase velocity (approximately 80 m/s at 10 Hz).

Estimated shallow S-wave velocity structures
The S-wave velocity structures were estimated by fitting the obtained dispersion curve and theoretical phase velocity dispersion curve of the fundamental mode of Rayleigh waves using the genetic algorithm inversion technique (Yamanaka and Ishida 1996).The P-wave velocity and density were calculated from the S-wave velocity relationship (Ludwig et al. 1970).The search ranges for the S-wave velocity and thickness are listed in Table 3.Since the range of observed phase velocities varies from station to station, the number of layers was set for each station.The search ranges for each station are shown in Supplementary Table S1.The number of generations, population size, crossover probability, and mutation probability of the inversions were set to 100, 20, 70%, and 1%, respectively.We performed 100 calculations with different random numbers and selected the structural model with the minimum mean squared error between the observed and theoretical values.Figure 3 compares the observed and theoretical phase velocities for the fundamental Rayleigh wave at the SNG station.
Figure 4 shows the estimated S-wave velocity structures, which were classified into three groups (see Fig. 2).The S-wave structures for each station are shown in Supplementary Table S2.Several S-wave velocity profiles do not reach the engineering bedrock with S-wave velocities of 300-700 m/s.The S-wave velocity in the uppermost layer is approximately 200 m/s for the stations situated on bedrock and sediments above the average elevation, while it ranges from 80-200 m/s in sites below   the average elevation.The KTP station, located on the isolated bedrock, had an S-wave velocity of 900 m/s at a depth of 10 m, in contrast to the other rock sites the JHR and SNG stations where S-wave velocities range from 300-400 m/s.Although these stations are located at a higher elevation and the geological map shows that these are the bedrock sites, the obtained S-wave structure shows that soft layers have been deposited on the surface.These soft layers might consist of weathered rock deposits or additional soil used in constructing school buildings.We must continue to examine the details.The S-wave velocities for sites located above the average elevation were 200-350 m/s at a depth of 10 m, while for those below the average elevation were 100-200 m/s at a depth of 10-30 m.These observed shallow S-wave velocities appear to be correlated with the elevation, and it was confirmed that no distinct velocity contrast in the shallow sedimentary layers.

Receiver function analysis
Figure 5 shows the radial and vertical component velocity waveforms recorded at the THM station on 22 July 2015 for the earthquake (M 4.3, depth: 7 km) that occurred near the Kathmandu basin.In the radial component, Ps wave is clearly observed between the direct P and S waves with approximately 1 s difference between the arrival of Ps and P waves (Ps-P time).The presence of distinct velocity contrast between the seismic bedrock and sediment (Koketsu et al. 2024), but not within the shallow sedimentary layers, suggests that a clearly noticeable Ps wave was generated at the interface of the sediment and seismic bedrock.We used a receiver function analysis (Langston 1979) to further investigate the Ps wave and its multiple reflections utilizing strong-motion waveforms from 45 earthquakes (Fig. 1b) with a magnitude of more than 3.5 from 2015 to 2019.Earthquake information was obtained from the United States Geological Survey (USGS) National Earthquake Information Center (NEIC) catalog, Horiuchi et al. (2021), Adhikari et al. (2015), and the National Earthquake Monitoring and Research Center (NEMRC), DMG catalog (Supplementary Table S3).These earthquakes were mostly aftershocks of the 2015 Gorkha earthquake.Due to the power-related issue at the KPN station (managed by the DMG), records were unavailable, consequently, data from 17 stations were utilized for the analysis.Although the operation periods differ among stations, more than five earthquakes were recorded at each station.
For the analysis, time windows of 5 s (starting from P-wave onset) or shorter than 5 s (from P-wave to S-wave onset) from radial and vertical components were selected.Both windows were zero-padded to 20.48 s with a cosine taper of 0.5 s on each end.The receiver function was calculated from a deconvolution of the radial component by the vertical component.The water level method (Clayton and Wiggins 1976) was applied to fill the spectral hole with a parameter c = 0.01, following Langston (1979).The obtained spectrum was then bandpass filtered at 0.2-5.0Hz.The average incident angle to the sediment from the seismic bedrock with an S-wave velocity of 2.2 km/s, taken from Koketsu et al. (2024), was approximately 45 degrees.The average ray parameter was 0.19 s/km.The phase-weighted stack (Schimmel and Paulssen 1997) was applied to the receiver functions for each station to reduce incoherent noise and obtain clear signals.
Figure 6 shows the receiver function for the THM station for the back azimuth.Clear phases appear at around 1.0 and 1.6 s.The approximately 1-s signal corresponds to the Ps-P time.The subsequent signal time corresponds to the difference in the arrival time of the P-wave and the wave that was reflected as a P-wave and converted to an S-wave at the interface (PpPs-P time).The theoretical Ps-P and PpPs-P times at the sediment-bedrock interface are 1.2 and 1.7 s, respectively, based on the velocity structure reported by Koketsu et al. (2024).The observed Ps-wave converted at the sediment-seismic bedrock interface was confirmed.The arrival time variation for each Ps phase is smaller regardless of the back azimuth or ray parameter values.The variation for the PpPs phase is larger.The stacked receiver function is also shown in Fig. 6.The amplitude of the Ps phase is clear, whereas that of the PpPs phase is somewhat unclear.
Figure 7 shows the receiver functions for all stations.The Ps-P times range from 0.1 to 1.0 s at the sedimentary sites, indicating that the basin has a complex bedrock topography.Positive Ps phases appear between approximately 0.6 and 1 s for the stations located in the center of the basin.The maximum Ps-P time was 1.0 s at the THM station.The Ps-P time decreases toward the edges of the basin.The Ps phases appear at less than 0.1 s at rock sites.Notably, although the KTP station (located on exposed bedrock) and the TVU station (located on the sediment in the center of the basin) are only 1.6 km apart, their Ps-P times differ by approximately 0.7 s, indicating a steep bedrock topography around these stations.The PpPs phase of the stacked receiver function for the THM and SGL stations (in the center of the basin) appears 0.6 s after the Ps phase.The Ps and PpPs phases approach each other and eventually overlap as the sedimentary layer at the stations becomes thinner.The arrival time difference between the Ps and PpPs phases for the RNB and TVU stations are 0.2-0.3s.The arrival times of the PpPs phase at other stations varies greatly and do not clearly appear in the stacked receiver functions.

Average S-wave velocity in the upper 30 m
The average S-wave velocity in the upper 30 m (V S30 ) is correlated with the site amplification at a high frequency.It is typically used as a parameter that represents the site effect in ground motion prediction equations.We calculated V S30 from the estimated S-wave velocity structures.For S-wave velocity profiles that did not extend to 30 m, V S30 was directly calculated from the phase velocity, following Konno and Kataoka (2000), who proposed that the phase velocity of Rayleigh waves with a wavelength of 35 to 40 m corresponds to the V S30 .The V S30 for the rock sites was above 320 m/s.The maximum value was 730 m/s at the KTP station.The mean of V S30 for the sites above the average elevation was 276 m/s, and that for below the average elevation was 165 m/s.The minimum value was 135 m/s at the TVU station.According to the seismic site classification (BSSC 2004) of the National Earthquake Hazards Reduction Program (NEHRP), our sites located on bedrock are either on class C (very dense soil and soft rock) or on class D (stiff soil); Similarly, sites above the average elevation are on class D, and sites below the average elevation are on class D or on class E (soft soil).
V S30 is commonly correlated with geological data.Since some previous studies (e.g., Yoshida and Igarashi 1984;Shrestha et al. 1998;Sakai 2001) have proposed different geological maps for the Kathmandu basin, it is difficult to classify using these.On the other hand, V S30 has been shown to correlate with the topographic slope (e.g., Wald and Allen 2007).Figure 8 shows the relationship between the V S30 and the topographic slope.We used World Elevation Data (30 m mesh version) from the Japan Aerospace Exploration Agency (JAXA) (Tadono et al. 2016).The slope for the site below the average elevation in the center of the basin has low gradients and the slope becomes steeper toward the edges of the basin.The V S30 and topographic slopes show a positive relationship.As an exception, the KTP station (located on the isolated bedrock inside the basin) has a low gradient but a high V S30 .In future studies, we will conduct many shallow microtremor array measurements to clarify the relationship between V S30 and parameters such as the topographic slope and geomorphological classification, as well as the relationship Fig. 8 Relationship between V S30 obtained from MASW and microtremor measurements and topographic slope between V S30 and the site amplification factor of the shallow velocity structure.

Relationship between receiver functions and velocity structures
Previous studies have proposed deep S-wave velocity structure models for the Kathmandu basin based on strong-motion records.Bijukchhen et al. (2017b) constructed one-dimensional sedimentary S-wave velocity structure models beneath eight strong-motion stations using aftershock records of the 2015 Gorkha earthquake.Bijukchhen (2018) extended the one-dimensional S-wave velocity structure models to an initial three-dimensional sedimentary S-wave velocity structure model for strongmotion simulation based on available geological and borehole data.Mori et al. (2020) proposed a site correction term for ground motion prediction equations using the same aftershock records and found that this coefficient is correlated with the bedrock depth.Koketsu et al. (2024) constructed a three-dimensional velocity structure model based on observations of microtremors and strong-motions, reflection survey data, and gravity data.
The structure was tuned using the radial-to-vertical spectral ratios of the strong-motion records from permanent stations (Table 1).Similar to the receiver function analysis using the vertical and radial components of the strong-motion records, the horizontal-to-vertical spectral ratio in the frequency domain is widely used for estimating sedimentary structure.While the spectral ratio is useful for detecting average velocity values, it is not as effective for determining the depth of the layer boundary.In contrast, the Ps-P time provides valuable information about the depth of the boundary (Chong et al. 2017).In this study, we focus only on the Ps-P time of the receiver function because the velocities for each layer have been determined in a previous study.The bedrock depth of the S-wave velocity (2.2 km/s) reported by Koketsu et al. (2024) is shown in Fig. 9.The sedimentary layer consists of five layers with S-wave velocities of 0.13, 0.19, 0.29, 0.44, and 0.51 km/s, respectively.The velocity contrast between the sedimentary layer and the seismic bedrock is large.The central part of the basin is the deepest region and there are several exposed rocks, so the basin has an undulating topography and small basins.The distribution of the Ps-P time from the observed receiver functions is also shown in Fig. 9.The Ps-P times are around 1 s in the central part of the basin with a depth of over 450 m.The vertical sections for line a − d of the sedimentary velocity structure reported by Koketsu et al. (2024) and the stacked receiver functions for stations along each line are shown in Fig. 10.The changes in the Ps phase correspond to changes in the bedrock depth.Although the PpPs phase appears at stations in the central part of the basin, which has a smooth topography, such as in lines b and d, it does not clearly appear at the other stations.The rapid change in the surrounding bedrock topography beneath the stations seems to affect the arrival times of the PpPs phase and the amplitudes of the stacked PpPs phases cancel out.
Theoretical arrival times were calculated from the one-dimensional sedimentary velocity structure beneath the stations reported by Koketsu et al. (2024).The theoretical Ps-P time at the seismic bedrock interface (2.2 km/s) appears in the observed receiver functions.The incident angle was used for the average observed value.The variation of Ps-P time with the incident angle is negligible.The observed and theoretical Ps-P times are generally in agreement.The theoretical Ps-P times are overestimated at the temporary stations, which were not used for tuning in Koketsu et al. (2024).The seismic bedrock depth around these strong-motion stations may be shallower than this model.Continuous validation of velocity structure models using additional datasets is required for accurately assessing the strongmotions of future earthquakes in this region.In the future, we will investigate the S-wave velocity structure from shallow to deep by inverting the receiver function, including its amplitude.

Conclusions
We investigated the Ps waves observed by the strongmotion station array recently installed in the Kathmandu basin using receiver function analysis.We confirmed that the observed Ps-waves were converted at the sediment-seismic bedrock interface.The rapid variation of the Ps-P time distribution suggests an undulating basin topography.We confirmed that the relative Ps-P times of the observed receiver functions generally correspond to the theoretical values reported in a previous study.Additionally, small-sized microtremor array measurements were carried out to target the shallow S-wave velocity structure around the strong-motion stations.We estimated the phase velocities using the SPAC method and constructed velocity structure models up to 30 m in depth.The V S30 were 135-730 m/s at 14 sites, and 165 m/s at sites on the lacustrine Kalimati Formation.The relationship between the V S30 and topographic slopes is positive.
The trend for the isolated bedrock site in the basin was different from that for other sites.Future work should examine the site amplification characteristics of the seamless velocity structure from the seismic bedrock to the surface based on the strong-motion records.

Fig. 1
Fig. 1 Overview of study area in Nepal.a Location map of study area.Black box indicates the location of map shown in b. b Distribution of earthquakes used for receiver function analysis.Epicenters are shown as circles (indicating scale).Black box indicates location of map shown in c. c Location of strong-motion stations plotted on engineering environmental geological map of Kathmandu basin based on the DMG (Shrestha et al. 1998).Blue triangles indicate permanent stations installed in 2011 (Takai et al. 2016) and blue inverted triangles indicate temporary stations operated from May to August 2015(Shigefuji et al. 2022).Red triangles indicate permanent stations installed in 2018(Takai et al. 2021)

Fig. 2
Fig. 2 Observed phase velocity dispersion curves.a Sites on bedrock.b Sites above average elevation of 1340 m. c Sites below average elevation

Fig. 3 aFig. 4 Fig. 5
Fig. 3 a Comparison of observed and theoretical phase velocities at the SNG station.Circles indicate the observed phase velocities.Phase velocities for the fundamental mode of Rayleigh waves are shown by a solid black line for the best-fit model and a thin grey line for a 10% model.b S-wave velocity structure.The solid black line is the best-fit model and the solid grey thin lines are the 10% models

Fig. 6
Fig. 6 Receiver functions for the THM station.The receiver functions and ray parameters arranged by back azimuth are shown on the left in the top and bottom panels.The stacked receiver function is shown in the right panel.The arrival time difference between the Ps (PpPs) waves and the P-wave is shown as t Ps (t PpPs )

Fig. 7
Fig. 7 Receiver functions for all stations.Stacked and individual receiver functions are shown as a solid black line and thin grey lines, respectively

Fig. 9
Fig. 9 Distribution of top depth of seismic bedrock with the S-wave velocity of 2.2 km/s reported by Koketsu et al. (2024) and Ps-P time obtained from each receiver function at strong-motion stations (circles).The dotted lines correspond to the profile for the vertical cross-section of the sedimentary layers shown in Fig. 10

Fig. 10 S
AbbreviationsDMG Department of Mines and Geology MASW Multi-channel analysis of surface waves M w Moment magnitude SPAC Spatial Autocorrelation USGS United States Geological Survey UTC Coordinated Universal Time V S S-wave velocity V S30 Average S-wave velocity in the upper 30 m

Table 1
Location of the strong-motion stationsHU-TU and HU-TU-T are permanent and temporary networks of collaboration between Hokkaido University and Tribhuvan University, respectively.DMG-SATREPS is a permanent strong-motion network installed under the SATREPS project and operated by the DMG

Table 2
Sites of microtremor and MASW measurements, the side length of microtremor array, interval of MASW, and obtained V S30

Table 3
Search range for genetic algorithm