Three-dimensional velocity structure models in and around the Kathmandu Valley, Central Nepal

We analyzed observations of microtremors and seismic ground motion to determine the depths of sedimentary interfaces and basement surface in the Kathmandu Valley. The results of these analyses and a reflection survey were combined to produce a data set of seismic depths for sedimentary interfaces and basement surface in the valley. Moreover, we introduced a data set of gravimetric depths, which had been produced for the basement surface from gravity observations. Because both sets of data differ in accuracy and number, we adopt the Sequential Gaussian Co-Simu-lation to generate the depth distributions of the sedimentary interfaces and basement surface. This method includes the cokriging approach, wherein the seismic and gravimetric depths are used as primary and secondary data, respectively. A three-dimensional velocity structure model of the Kathmandu Valley was constructed from the obtained depth distributions, which showed that the sedimentation in the valley is mainly related to the Paleo-Kathmandu Lake. Subsequently, we incorporated this model into a previous regional-scale model, and revised the latter to solve the issues of the extension of the underthrusting lithosphere and the S-wave velocity in its adjacent area.


Introduction
Nepal is an earthquake-prone country and its capital is located in the Kathmandu Valley.It is therefore important to accurately assess the seismic hazards there, and ground motion simulations based on accurate velocity structure models in and around the valley are essential for such assessments.
The collision of the Indian plate with the Eurasian plate is causing the Himalayan orogeny, which has formed many valleys and basins, as shown in Fig. 1a.The largest among them is the Kathmandu Valley.This valley is filled with sediments from the Bagmati River (thick blue lines in all panels of Fig. 1) and its tributaries (medium blue lines in Fig. 1b, c).The Paleo-Kathmandu Lake was formed 1 million years ago by the damming of the Bagmati River, and dried 12,000 years ago (Sakai et al. 2016).This geological history has resulted in the velocity structure of the Kathmandu Valley being three-dimensionally complicated.
The collision of the plates also causes the Indian lithosphere to underthrust beneath Nepal along the Main Himalayan Thrust, which reaches the ground surface at the Main Frontal Thrust (MFT; red lines with triangles in Fig. 1a).This underthrusting generates earthquakes and three-dimensionally complicates the velocity structure around the Kathmandu Valley.Therefore, in this study, we construct the three-dimensional velocity structure models in and around the Kathmandu Valley using diverse data sets and information.

Modeling in Kathmandu Valley (data)
In the Kathmandu Valley, Yokoi et al. (2021) performed broadband array observations and analyses of microtremors, and Takai et al. (2021) recorded ground motions from earthquakes with strong motion seismometers.The locations of microtremor observations and strong motion seismometers are plotted in Fig. 1b, c with red and magenta circles, respectively.The details of these observations will be reported in separate papers.The gravity observations were also performed as described by Pradhan et al. (2018).The violet dots in Fig. 1b, c indicate these observation sites including those of Moribayashi and Maruo (1980).Kawasaki et al. (2019) performed a seismic reflection survey for the first time in the Kathmandu Valley along a line extending adjacent to the runway of the Tribhuvan International Airport (TIA).The airport and survey line are shown in Fig. 1b, c with dark green areas and light green lines.The details of the gravity and seismic reflection surveys will also be reported in other separate papers.We first evaluated the Rayleigh-wave phase velocities and horizontal-to-vertical spectral ratios of microtremors, which were provided by Yokoi et al. (2021), for S-wave velocity (V S ) structures using the joint inversion method of Suzuki and Yamanaka (2010) and the assumption of Arai and Tokimatsu (2000).The assumption is that microtremors are the summation of surface waves from point sources distributed at all directions and distances.Microtremor data do not have resolution for deep parts of the V S structure.Therefore, we set the V S of the two layers in the basement as 2.2 and 2.6 km/s referring to crustal structures in the Eurasian plate, although we determined both the V S and the thicknesses of the five sedimentary layers above them.Additionally, we assumed fixed structures of the crust and upper mantle based on CRUST 1.0 (Laske et al. 2013).The results at the TIA site are shown in the center of Fig. 2, and all results are shown in Additional file 1: Fig. S1.
This site is located close to the reflection survey line as shown in Fig. 1c, so that the elevation profile obtained from the survey by Kawasaki et al. (2019) is shown in the background of Fig. 2. Their interpretation of the basement surface is indicated by the green line.The black line indicates the basement surface along the reflection survey line obtained from the gravity survey performed by Pradhan et al. (2018).The shapes of the basement surface are similar to each other and they agree with the interface between the sediments (red) and basement (green) in the result of the microtremor analysis.
We then analyzed the observed ground motions from earthquakes for V S structures at seismometer stations.
The observed north-south and east-west components were converted to radial and tangential components using the earthquake epicenter location determined by the United States Geological Survey or Horiuchi et al. (2021).A Fourier transform was performed on the coda parts of this radial component and the observed vertical component to obtain the R/V spectral ratio.For example, seven earthquakes with determined epicenters were observed at the SGL station (Fig. 1b, c), and their R/V spectral ratios were calculated as indicated by the gray lines in Additional file 1: Fig. S2c.The average of these R/V spectral ratios was considered as the observed R/V spectral ratio (black line).This station is located close to the TIA site, as shown in Fig. 1c.Therefore, starting with the velocity structure at TIA (Additional file 1: Fig. S2a), the velocity structure at SGL was improved through trial and error (Additional file 1: Fig. S2b) to match the peak of the theoretical R/V spectral ratio of the Rayleigh wave to that of the observed ratio (magenta and black lines in Additional file 1: Fig. S2c).In this process, because of limited performance of R/V spectral ratio data, trial and error was performed by changing only the thicknesses of the sedimentary layers.The layer thicknesses at SGL were determined to be 80% of those at TIA (two-headed lines in Additional file 1: Figs.S2a, b).All the results of this analysis are shown in Additional file 1: Fig. S3.
The depths of the layer interfaces were determined from the obtained V S structures shown in Additional file 1: Figs.S1 and S3.The interface at the top of the 2.2 km/s layer is the basement surface, and the interfaces shallower than this are called sedimentary interfaces.The  2018) is also shown with the black line.The embedded column represents the S-wave velocity structure obtained from our microtremor survey at the TIA site basement surface was also obtained from the results of the seismic reflection survey (Fig. 2) and borehole survey results of JICA ( 2018).The depths of the basement surface from the microtremors, ground motions, reflection survey, and borehole surveys are highly accurate but limited in number, whereas those from gravity observations by Pradhan et al. (2018) are numerous, as shown in Fig. 1b, c.The latter are called gravimetric depths and are plotted along with the former, collectively called seismic depths, in Fig. 3.This figure shows that the gravimetric depths have fair accuracy, although there are some variations.
According to the similarity and accuracy shown in Figs. 2 and 3, we use both the seismic and gravimetric depths for the modeling of the three-dimensional velocity structure in the Kathmandu Valley.It is noted that the microtremors and ground motions also provided us with the seismic depths of the sedimentary interfaces.In addition, to further utilize the advantage of the large amount of gravity data, a grid was defined over the area including the Kathmandu Valley, and Bouguer anomalies at its nodes are obtained by interpolation and extrapolation of the observations using the minimum curvature algorithm (Smith and Wessel 1990).The gravimetric depths calculated from them (Pradhan et al. 2018) were used in this study.
We assumed seven layers for the sediments and basement based on the V S structures in Additional file 1: Figs.S1 and S3.We adopted the accurate V S structures obtained from the joint inversions of the microtremor data (Additional file 1: Fig. S1), but excluded Chobar, because its shallow V S structure is very different from those of other stations.V S of an assumed layer was determined by averaging the values of the corresponding layers in the V S structures, as shown in Additional file 1: Table S1.The P-wave velocity (V P ), density (ρ), and S-wave and P-wave quality factors (Q S and Q P , respectively) for each layer were obtained using the empirical relationships by Kitsunezaki et al. (1990), Ludwig et al. (1970), and Kawabe and Kamae (2008).

Modeling in Kathmandu Valley (methods)
The shapes of the basement surface (basement/sediment interface) and the four sedimentary interfaces separating the five upper layers were then determined in the form of a depth distribution.This determination is a variant of the common problem of mapping a regionalized variable from limited sample of data.We considered the depth distribution of an interface as a regionalized variable z(x) at unsampled location x , and assumed that seismic depths z i were obtained at sampled locations x i (i = 1, 2, . . ., n) .To recover z(x) from z i , we use the 'kriging' method developed by Matheron (1963), which was named after pioneering work of Daniel Krige in the early 1950s.z(x) was statistically distributed in region A, which was the Kathmandu Valley, with a probability density function.We adopted the normal score transformation of Deutsch and Journel (1997) such that z(x) was converted to y(x) , the probability density of which is nor- mally distributed with zero mean and unit variance.We note that the expectation value of a function g(x) includ- ing y(x) is defined as where f (x) is the normal distribution with zero mean and unit variance.
It is assumed that z i is equal to the true variable at x i , which is z(x i ) or y(x i ) in the normal score space, and that y(x) is represented by a linear combination of y(x i ) as Because the integral in (1) is a weighted sum of g(x) in A, the variance (i.e., the weighted sum of squared residuals) for the observation equation ( 2) can be written as Fig. 3 Plot of gravimetric depths against seismic depths for the basement surface.The blue, green, and red dots represent gravimetric depths at the microtremor sites and seismometer stations, those along the reflection profile, and those at the boring sites, respectively.They are located at the horizontal positions corresponding to the seismic depths from the microtremor, ground motion, or boring data and those from the result of the reflection survey The covariance and variogram of y(x) are defined as They are related to each other by where γ (h) is a half of the variogram and called a semi- variogram.S is then rewritten using C in (3) as (3) S is minimum for the solutions of the normal equations ∂S/∂w i = 0 , i = 1, 2, . . ., n.
This formulation is extended to 'cokriging' where data with different attributes were included (e.g., Doyen 1988).We included the gravimetric depths (2) then yielded the observation equation: and w ′ i ′ is the linear combination coefficient for y ′ x ′ i ′ .The weighted sum of squared residuals for (4) is

Using
(2), C ′ (h) = E y(x)y ′ (x + h) , and which are related to each other by the variance S in ( 5) is rewritten as This yields the partial derivatives and the normal equations The number of covariances included in this system of linear equations is as many as n + n ′ 2 + n + n ′ , which requires much computation.However, because the volume support of the gravimetric depth is larger than that of the seismic depth, it can be assumed that the statistical dependence of seismic depth on gravimetric depth is limited to the co-located depths.This approximation is called the Markov Model 2, and results in linear relations and where ρ 12 is the correlation coefficient inferred from co-located data y(x) and y ′ (x) , and ρ R (h) is any permis- sible covariance (Journel 1999).If C ′′ is calculated from the gravimetric depths, all the covariances in (8) can be obtained using ( 9) and ( 10).Therefore, we can solve (8) and calculate the cokriging estimate and variance using the solution of (8), the observation equation ( 4), and the normal score back-transformation.It is additionally noted that because the probability density function for y(x) and y ′ (x) is the normal distribution with zero mean and unit variance.
To obtain the depths at the nodes of a grid as the depth distribution of an interface, we adopted the procedure of the Sequential Gaussian Co-Simulation (Remy et al.

2009) as follows:
1.Each node x of the grid on an interface is defined.

Modeling in Kathmandu Valley (results)
The depth distributions of the five interfaces, which comprised the basement surface (basement/sediment interface) and the four sedimentary interfaces, were determined using the data and methods described in the previous two sections.In Step 1 of the above procedure, a grid with 200 m intervals was applied to the rectangular region including the Kathmandu Valley.This grid and its nodes ( x k , k = 1, 2, . . ., K ) were commonly used for all the five interfaces.In Step 2, we took the normal score transforms of the gravimetric depths for the basement surface ( ) and those of the seismic depths for an interface ( y(x i ), i = 1, • • • , n ), and constructed the initial data set by combining y ′ x ′ i′ and y(x i ) .The correlation coefficient ρ 12 of the normal score transforms of the co-located seismic and gravimetric depths were calculated in advance.
In 3.1 of Step 3, we proceeded to a random node at x k in the grid.We first calculated the variogram γ ′′ (h) of y ′ x ′ i′ since we used the approximation of Markov with a range a = 8000 m fit them.C ′′ (h) s in (8) were cal- culated using ( 6), ( 11) and γ ′′ (h) = γ s (h) in ( 12) with the range.C ′ (h) s in (8) were calculated using ρ 12 calculated in Step 2 and (9).We assumed ρ R (h) in ( 10) to be 1 − γ s (h) , where γ s (h) is the spherical model in ( 12) with a = 8000 m.C(h) s in (8) are calculated using this assumption and (10).Thus, the system of linear equations in (8) was set up and solved for w 1 , . . ., w n and w ′ 1 , . . ., w ′ n ′ .From these solutions and the observation equation ( 4), we obtained the cokriging estimate (estimate by the cokriging) for y at x k : Using the solutions together with (7), we obtained the cokriging variance (variance by the cokriging) S(x k ) for y at x k .
In 3.2 of Step 3, as reported in Pyrcz and Deutsch (2014), we drew a random residual R(x k ) that followed a normal distribution with zero mean and variance of S(x k ) using a pseudo-random number generator for a given "seed".We added the cokriging estimate and the residual to get the simulated value y s (x k ) = y * (x k ) + R(x k ) .At the end of the for loop in Step 3, we add y s (x k ) to the (12) initial data set to ensure that the covariance with this value and all future predictions is correct.
We revisited the beginning of the for loop with a next random node at x k ′ in the grid.The processes of 3.1 and 3.2 were repeated for this node and the data set after the previous iteration to obtain the simulated value y s (x k ′ ) and added this to the data set.We then visited all K nodes in a random order to generate the simulated values.In Step 4, all the simulated values in the normal score units were finally back-transformed to the distribution of depths of each interface, which was collectively called a realization.Steps 3 and 4 were repeated 50 times with different random number seeds, and the average of the 50 realizations was used as a model for the depth distribution of each interface.Figure 4 shows the models for the four sedimentary interfaces and the basement surface.To show the level of confidence for them, we drew the distribution of the standard deviation for the depth of the basement surface in Additional file 1: Fig. S5.
Comparing Figs.1b and 4, we found deep interfaces, which correspond to thick sediments, in the center of the valley and along the Bagmati River and its tributaries.In the southern half of the valley, the thickest sediments (deepest interfaces) are not distributed along the Bagmati River but extend south from where the river runs from east to west.The southern extension is located in an area where lacustrine sediments and delta deposits are distributed (Sakai et al. 2016).This suggests that sedimentation in the valley is predominantly related to the Paleo-Kathmandu Lake.
To observe the effect of the valley on seismic ground motion, we simulated ground velocities from an M 5.4 earthquake in the east at points A and B, which represent the edge and center of the valley (left panel of Fig. 5).We used a 3D finite difference method with the fourth-order difference scheme (Levander 1988), a staggered grid (Virieux 1986), and viscoelastic formulation of Robertsson et al. (1994).A bandpass filter from 0.1 to 0.5 Hz was applied to the simulation results.The ground velocities at point A had simple waveforms consisting mainly of body waves, whereas those at point B had complex waveforms with long duration and large amplitudes, including surface waves that developed while propagating through the valley (right panel of Fig. 5).Koketsu et al. ( 2016) constructed a regional-scale velocity structure model around the Kathmandu Valley to simulate the ground motions caused by the 2015 Gorkha earthquake.The target region is indicated by the large rectangle in Fig. 1a.This model was based on the geometry of the Indian and Eurasian plates by Sapkota et al. (2013) and the crustal structures in CRUST 1.0 developed by Laske et al. (2013).However, CRUST 1.0 provides no basement on the Nepalese side along the MFT; therefore a layer of V S 2.1 km/s was introduced as on the Indian side.The shape of the Kathmandu Valley was defined according to and Igarashi (1984), and the valley was filled with uniform sediments of V S 0.6 km/s.The V S of the Indian lithosphere was assumed to be 3.2 km/s.Apart from these, the V S, V P , and ρ of each layer were obtained from CRUST 1.0.Q S was obtained from V S using the empirical relationship of Brocher (2008), and Q P was assumed to be twice Q S .

Modeling around Kathmandu Valley
The underthrusting Indian lithosphere was modeled to a depth of 8.4 km in this model (Fig. 6a), but extended beyond that depth to the northern edge of the modeled region, as shown in Fig. 6b.This extension was performed following well-known reviews on Himalayan tectonics (e.g., Avouac 2003;Bilham 2019).The uniform sediments in the Kathmandu Valley (Fig. 6a) were replaced with the velocity structure model obtained in the previous sections, as shown in Fig. 6b.If we use the old model in the Kathmandu Valley, for example, the calculated phase velocities and R/V spectral ratios at SDB (green lines in Additional file 1: Fig. S6) do not match the observed ones (gray circles), but those calculated with our new model (violet lines) well fit to the observed ones.Kobayashi et al. (2016) found that the source region of the Gorkha earthquake was located in the V S 3.5 km/s layer of their velocity structure model.Because large historical earthquakes occurred even on the shallow part of the underthrusting lithosphere, the V S 3.5 km/s layer on the Nepalese side was raised to below the MFT, as shown in Fig. 6b.

Conclusions
We analyzed the microtremors observed by Yokoi et al. (2021) and the seismic ground motions observed by Takai et al. (2021) to obtain horizontally layered structures at their observation sites.Kawasaki et al. (2019) conducted a reflection survey to obtain a two-dimensionally layered structure along a survey line.From these results, we obtained seismic depths of four sedimentary interfaces and basement surface (basement/sediment interface).Pradhan et al. (2018) performed gravity observations to obtain the gravimetric depths of the basement surface.The accuracy of the seismic depths is higher than that of the gravimetric depths, while the number of the former is smaller than that of the latter.To generate the depth distributions of the interfaces from these depths with such different attributes, we adopted the Sequential Gaussian Co-Simulation (Remy et al. 2009) including the cokriging approach (e.g.Doyen 1988).We combined the results of this method and assumptions for deeper parts (Additional file 1: Table S1 suggests that the sedimentation in the valley is predominantly related to the Paleo-Kathmandu Lake.We then incorporate this model into a regional-scale model for ground motion simulation (Koketsu et al. 2016).We also revised the regional-scale model to solve the problems of the extension of the underthrusting lithosphere and V S in its adjacent area.

Fig. 1
Fig. 1 Topographic maps for Nepal and the Kathmandu Valley (a and b).c is a close-up map around the valley center and airport (dark green area).The large black box and red lines with triangles in a represent a wide area around the valley and MFT, respectively.The light green line labeled X and Y indicates where the cross sections in Fig. 6 are taken.The thick lines in blue indicate the Bagmati River in all maps, while thin lines in medium blue indicate other major rivers in a or tributaries of the river in b and c.The Kathmandu Valley and its internal sediments are surrounded by the basin boundary (brown) and sedimentary boundary (black) of Shrestha et al. (1998) in b.In b and c, the sites for microtremor observations, stations of strong motion seismometers, and locations of boring survey results are plotted with red, magenta, and orange circles, respectively.The violet dots are points of gravity observation, and a seismic reflection survey is conducted along the very thick lines in light green

Fig. 2
Fig.2Elevation profile obtained from the seismic reflection survey byKawasaki et al. (2019).Their interpretation for the basement surface is shown with the green line.The basement surface obtained from the gravity survey byPradhan et al. (2018) is also shown with the black line.The embedded column represents the S-wave velocity structure obtained from our microtremor survey at the TIA site Modeled depth distributions of the four sedimentary interfaces (a to d) and the basement surface (e).The brown line in each distribution indicates the basin boundary of Shrestha et al. (1998) Model 2. The calculated values are plotted as red crosses in Additional file 1: Fig. S4, and the black line in this figure indicates that a spherical model:

Fig. 5
Fig. 5 Right panel compares the simulated ground velocities from an M 5.4 earthquake in the east at points A and B, which represent the edge and center of the basin as shown in the topographic map on the left.The brown and black lines in the map indicate the basin and sedimentary boundaries of Shrestha et al. (1998), respectively

Fig. 6
Fig.6Panels a and b show the cross sections of the regional-scale velocity structure models byKoketsu et al. (2016) and this study, respectively, along line XY in Fig.1a.The underthrusting Indian plate is modeled to a depth of 8.4 km in a, but extends beyond that depth to the northern edge of the modeled region in b.The sediments in the Kathmandu Valley are assumed to be homogeneous with V S of 0.6 km/s in a, while the velocity structure model obtained in the previous sections is applied in b.The V S 3.5 km/s layer on the Nepal side is raised to below the MFT (dotted line) based onKobayashi et al. (2016)