Mapping long-period soil resonances in the Kathmandu basin using microtremors

This study reports the geostatistical analysis of a set of 40 single-station horizontal-to-vertical spectral ratio (HVSR) passive seismic survey data collected in the Kathmandu basin (Nepal). The Kathmandu basin is characterized by a heterogeneous sedimentary cover and by a complex geo-structural setting, inducing a high spatial variability of the bedrock depth. Due to the complex geological setting, the interpretation and analysis of soil resonance periods derived from the HVSR surveys is challenging, both from the perspective of bedrock depth estimation as well as of seismic-site effects characterization. To exploit the available information, the HVSR data are analyzed by means of a geostatistical approach. First, the spatial continuity structure of HVSR data is investigated and interpreted taking into consideration the geological setting and available stratigraphic and seismic information. Then, the exploitation of potential auxiliary variables, based on surface morphology and distance from outcropping bedrock, is evaluated. Finally, the mapping of HVSR resonance periods, together with the evaluation of interpolation uncertainty, is obtained by means of kriging with external drift interpolation. This work contributes to the characterization of local seismic response of the Kathmandu basin. The resulting map of soil resonance periods is compatible with the results of preceding studies and it is characterized by a high spatial variability, even in areas with a deep bedrock and long resonance periods.


Introduction
Kathmandu, the capital of Nepal (Fig. 1), has been exposed to the effects of various earthquakes. The last one, the M7.8 Gorkha earthquake (25 April 2015), struck central Nepal causing severe damages in the city, including the monumental heritage (Kathmandu is listed in the UNESCO World heritage). The urban area is sprawl in an intermountain sedimentary basin, once partially covered by the waters of an ancient lake (Igarashi et al. 1988;Sakai et al. 2001). The sedimentary cover, characterized by variable thickness and stratigraphic heterogeneity, and the complex geo-structural setting influence and diversify spatially the seismic response Takai et al. 2016). From this perspective, the Kathmandu sedimentary basin has some analogies with the geological structure and seismic response of Mexico City, including the occurrence of seismic amplifications at low soil resonance frequencies. In fact, also Mexico city is built on a sedimentary basin once occupied by a lake and it is characterized by a complex seismic behavior, as testified by the various researches on seismic local effects, fueled by the experience of the Michoacán M8.1 September 19, 1985 earthquake (e.g., Singh et al. 1988). However, relevant aspects differentiate the two basins. One difference is related to the temporal evolution of the lacustrine environment. In Mexico City, some lakes of limited extent and with strong anthropic influence are still present; differently, the Kathmandu valley dried up in the Late Pleistocene (Igarashi et al. 1988). Another relevant difference is related to the stratigraphy: Mexico City sedimentary basin has a maximum depth of approximately 80 m, with relatively low Vs (150 m/s) and it is mainly constituted of clays with high water content Sáncez-Sesma et al. 1988). Differently, Kathmandu sedimentary basin thickness exceeds 500 m and it is characterized by a high lithological heterogeneity (Moribayashi and Maruo 1980;Khadka 1993) and by higher Vs of the sedimentary cover, with an average Vs of the first 30 m and 200 m, respectively, of 200 m/s and 250 m/s (Molnar et al. 2017;Sandron et al. 2019). Consequently, the seismic amplifications at low resonance frequencies (or long periods) are mainly due to the low Vs of the sediments in the case of Mexico City and to the high thickness of the sedimentary cover in the case of Kathmandu.
The Kathmandu area is peculiar and challenging from the seismological perspective (e.g., Wesnousky et al. 2017), in relation to the geo-structural setting and to the heterogeneity of the sedimentary cover. Kathmandu was hit by numerous historical earthquakes (Ambraseys and Douglas 2004;Sandron et al. 2019) and many of these events caused relevant damages, even if the epicenters were located at distances longer than 200 km. The most important recent seismic event in the area was the 1934 Nepal-Bihar M8.2 earthquake (Sapkota et al. 2013), that, with an epicenter distance of approximately 200 km, caused damages in the capital, with a calculated intensity MSK (Medvedev Sponheuer Karnik Intensity) of VIII (Ambraseys and Douglas 2004). The Shillong M8.1 earthquake of 1897 is still more emblematic; it caused an intensity of V with an epicentral distance of 600 km. Notwithstanding, the epicenter of Gorkha 2015 M7.8 earthquake was closer than in the Nepal-Bihar 1934 M8.2 earthquake, historical heritage and the vernacular dwellings sustained comparable damages. In addition, a relevant spatial variability of the macroseismic intensity has been observed; some sites of the valley reached VIII, but most of the assessments were in the range of VI-VII . A similar seismic behavior has been observed in Mexico City, in the case of the Michoacán M8.1 earthquake; for example, even if the epicenter of this event was very far from the City, approximately 350 km, the lacustrine sediments amplified the seismic waves (long wavelengths) and phenomena of soil liquefaction caused foundations failure of many buildings. Some areas of the Kathmandu basin seem to amplify seismic shaking with periods longer than 3 s, as testified by the strong motion data recorded during the 2015-2016 sequence (Dixit et al. 2015;Bhattarai et al. 2015;Takai et al. 2016;Rupakhety et al. 2017). In the cited studies, a sensor network based on 6 accelerometers and at least 3 GPS stations has been considered. These studies highlight the appearance of ground plastic deformation as response to the main shock, mainly related to the low-frequency component of the accelerometer spectrum. Later Takai et al. (2016), analyzing the data from three accelerometers located in the center of the basin and one located near them but on the bedrock, identified some long-period (T 3-5 s) and high-amplitude peaks on the spectrum and considered these an inherent seismic characteristic of the basin. This was confirmed also in Rupakhety et al. (2017) study, which emphasizes the presence of long-period peaks in the instrumental records; in addition, the polarization analysis of T 4 s pulse of accelerometers on sediments against the one on bedrock seems to confirm the presence of a characteristic amplification in the band of 3-5 s. Moreover, Dhakal et al. (2016) show that the attenuation relations for the Himalayan region fail in the over mentioned band of periods, while elsewhere correctly predict the ground motion periods. In conclusion, the similar response to low-frequency seismic waves of Mexico City and Katmandu is directly related to the combination of sedimentary cover thickness and seismic velocity.
Given the complex seismic behavior of the Katmandu area, it is worth to analyze the spatial distribution of possible long-period soil resonances, both for the study of seismicsite amplifications as well as for gathering local stratigraphic information. In this context, passive seismic techniques based on seismic noise recordings are particularly well suited for applications in urban environments (D'Amico et al. 2008;Avila et al. 2016;Trevisani and Boaga 2018), being this methodology inherently non-destructive and easily applicable with lightweight equipment. Among these methodologies, the single-station seismic noise approach based on Horizontal to Vertical Spectral Ratio (HVSR, Nakamura 1989) is particularly popular. As reported in subsect. "The HVSR methodology", the HVSR analysis permits to detect soil resonance frequencies and an indication of the seismic impedance contrast, furnishing relevant information on seismic response and on local stratigraphy. In this context, the studies of Pandey (2000) and Paudyal et al. (2012Paudyal et al. ( , 2013 are noteworthy for the Kathmandu area. The former study is based on single-station seismic noise recordings, with 60 measurements, acquired in 1985-86; the signals have been analyzed considering the power spectrum and the ratio with a reference bedrock spectrum. The latter, conducted by a joint team of Nepali and Japanese researchers, is based on an extensive survey of 172 single-station measurements in the northern/central part of the city. In both studies, the duration of recording was shorter than 5 min and the signals have been analyzed with window lengths of approximately 20 s; consequently, the reliable frequency response band lower limit is roughly of 1 Hz or more (Sesame 2004). Unfortunately, with this kind of experimental setting, it is possible to detect only short periods of HVSR peaks, mainly related to shallow impedance contrasts internal to the sedimentary cover. For example, Paudyal et al. (2012Paudyal et al. ( , 2013 recorded only one site with a period longer than 2.0 s (0.5 Hz) (site number 116), and all the other recordings detected higher frequency peaks (Sandron et al. 2019). In fact, as reported in Paudyal et al. (2013), the detected HVSR peaks are referred to seismic reflectors located in the sedimentary sequence, hundreds of meters above the effective bedrock surface. This is also highlighted by the study of Gilder et al. (2020), where a preliminary analysis of an extended stratigraphic database confirms that the "bedrock" surface detected by Paudyal et al. (2012) is much shallower than the expected depth of lithological bedrock. Paudyal et al. (2012Paudyal et al. ( , 2013 interpret the detected seismic reflectors as a seismic bedrock and not as the geological bedrock. However, as reported above, the recordings equipment and surveys setting adopted by Paudyal et al. (2012) and by Pandey (2000) prevent to detect deep seismic reflectors. Accordingly, the derived seismic surface of Paudyal et al. 2012, even if significative for local seismic response, is not representative of the seismic bedrock nor of the lithological bedrock (Sandron et al. 2019). Unfortunately, the deep stratigraphy and the effective bedrock top morphology is still not well known. Moreover, most of the deep boreholes (e.g., Gilder et al. 2020) have been conducted for ground water purposes by means of a non-core drilling, and consequently the stratigraphic records have a high level of uncertainty and potential interpretative ambiguities. Even the seismic analyses conducted by means of passive methods (e.g., Molnar et al. 2017;Poovarodom et al. 2017;Sandron et al. 2019) for deriving Vs depth profiles have a limited areal coverage and often a shallow investigation depth.
To characterize long resonance periods, this study considers a set of 40 HVSR analyses specifically designed for the detection of deep seismic reflectors. 31 HVSR data (Sandron et al. 2019) have been collected during two surveys (2015 and 2018) promoted by the National Society for Earthquake Technology NSET (Nepal), in collaboration with National Institute of Oceanography and Experimental Geophysics (OGS, Italy). To integrate the number of recordings capable to detect long-period peaks, other 8 data are derived from the work of Molnar et al. (2017). Then, a measurement located on rock (referred as KTP, Fig. 2, see Takai et al. 2016) has been included as reference for outcropping bedrock. The combined dataset of 40 measurements makes feasible the characterization of the western part of the basin, covering the downtown of Kathmandu and important areas of the city (Lalitpur). The studied spatial domain covers also rural areas affected by a rapid urban sprawl. This area is characterized by a very heterogeneous geological structure, conditioned also by the dendritic morphology of the southern edge of the basin, with the main drainage valley of the river Bagmati.
The aim of this work is to analyze the spatial distribution of soil resonance periods referred to deep reflectors, highlighting potential relationships with the geological and geomorphological setting. To achieve this main goal, a geostatistical approach is adopted, by means of three intermediate objectives. A first objective is the analysis of the spatial continuity of resonance periods, also evaluating its compatibility with available geo-structural knowledge. A second objective is to test the possible use of secondary data to improve the spatial mapping of resonance periods. Then, the final objective is to build a map of resonance periods related to deep seismic reflectors, together with a measure of prediction uncertainty. The produced map would be useful for improving the knowledge concerning the seismic response and the geological structure of the Kathmandu basin.

Geology of Kathmandu basin
The geo-structural setting of Kathmandu basin is complex, both regarding the bedrock, strongly deformed during the orogeny, as well as concerning the thick sequence of sediments filling the Kathmandu basin (e.g., Suresh et al. 2018;Gilder et al. 2020). The Kathmandu valley is an intermountain basin, developed within a Nappe structure located in the Lesser Himalaya belt, once occupied by a lake until the Late Pleistocene. The bedrock is composed mainly of metamorphic rocks (Precambrian-Palaeozoic; Paudyal et al. 2012;Gilder et al. 2020): meta-sedimentary, in the southern and central part, gneiss, schists and granite in the north. Lesser Himalaya (Fig. 1) is a geological region included between Main Central Thrust (MCT) in the north and Main Boundary Thrust (MBT), two important tectonic features' surficial expression of the Main Himalayan Thrust (MHT). Fluviolacustrine sediments and alluvial fan deposits have filled the basin (Sakai et al. 2006(Sakai et al. , 2015, in a period between Late Pliocene till the Quaternary, with more than 600 m of sediments (Moribayashi and Maruo 1980) in the center of the Basin. Near Ratna Park, in the city center, (Fig. 2, borehole BR), close to Bhrikutimandap Exhibition Hall, a borehole encounters the basement rock at 550 m of depth (Khadka 1993). As reported in Khadka (1993) and from the analysis of the stratigraphic database presented in Gilder et al. (2020), the transition from the sedimentary cover to the bedrock is often not sharp, and tens of meters of weathered and fractured rock can be found before the fresh bedrock. The sedimentary architecture is very complex and heterogenous (Dill et al. 2001) with the presence of fluvial paleo-channels and, on the margins of the basin, heterogeneous slope deposits (Fig. 2). Moreover, active tectonics profoundly influenced sedimentary and geomorphological processes, including the drainage system (e.g., Moribayashi and Maruo 1980;Asahi 2003;Sakai et al. 2006Sakai et al. , 2015. This is for example evident in western sector, where the active Chandragiri and Chobar faults ( Fig. 2) conditioned the sedimentary sequence (Saijo et al. 1995) and the bedrock top morphology.
The more relevant sedimentary formations ( Fig. 2) are, in ascending order, Tarebhir (older than 2.8 Myr), Lukundol (2.8 to 1 Myr) and Kalimati (2.8 to 30 Kyr) (see Paudyal et al. 2012, Table 1 modified after Sakai et al. 2001). The Tarebhir Formation (Sakai et al. 2001 is mainly composed of fluvial gravel and sand deposits. Then in the south, Lukundol and Kobgoan Formations consist of muddy sediments. Finally, the Kalimati Formation (from 2.8 to 30 Kyr) is mainly composed of black organic clay, silts, sand and gravel. It is the most important formation, locally referred as "black mud". The Kalimati Formation is overlaid by the Chapagoan Formation (up to Pleistocene, Sherstha and Shah 2014).These formations are overlaid by the Patan formation (from 14-19 to 10 Kyr, last glaciation, Paudayal 2011), consisting of fluvial deposits (sand, silt and gravel). In the north, the Gokarna Formation (Middle Pleistocene, 1 to 30 Kyr) is also present and is composed of gravel and sands. Gokarna Formation is then overlaid by the Thimi Formation, mainly composed of silt and clays and in the north-west sector by Tokha Formation, mainly sands (Sherstha and Shah 2014). The spatial distribution of the geological Formations is complex, with heteropic indentations (for more details see Dill et al. 2001;Sakai et al. 2001Sakai et al. , 2006Sakai et al. , 2008Sakai et al. and 2015. The sedimentary architecture is complex even in relation to the seven climatic oscillations, glacial and interglacial fluctuations, occurred during its deposition from middle to   (Kuwahara et al. 2004;Igarashi et al. 1988;Sakai et al. 2001).
These oscillations generated an alternation of draining and filling of water of the Paleo Kathmandu Lake. For example, the Patan formation was affected by a period of cold and drier climate (glacial) with a drainage of the lake, as confirmed by the studies of Mukunda and Sakai (2008) and Mampuku et al. (2008), based on radiocarbon dating.

The HVSR methodology
The seismic environmental noise, that for a long time in seismic reflection/refraction survey methodologies was considered a disturbance, is currently exploited in environmental seismology and extensively adopted in seismic hazard-related studies. The seismic noise sources (e.g., Longuet-Higgins 1950;Bonnefoy-Claudet et al. 2006) are considered to be natural (oceanic waves, meteorological perturbations) for very low frequencies (< 0.5 Hz) and anthropic (urban and industrial related activities, e.g., Poli et al. 2020) for higher frequencies (> 1 Hz). The seismic waves have amplitudes in the range of 10 -6 -10 -5 m. Nogoshi and Igarashi (1970a, b) introduced the use of seismic noise in geophysics; it was initially devised to achieve the resonance frequency (f 0 ) of a stratigraphic layer in presence of an impedance contrast (e.g., a soft layer overlying the bedrock). Then, Nakamura (1989) implemented the Horizontal to Vertical Spectral Ratio (HVSR) method, based on the recording of seismic noise by means of three components highly sensitive geophones (tromograph). The key idea is to plot the curve of the ratio between the spectra of horizontal and vertical velocimetric components versus the frequency (or the period) and individuate one or more peaks in the curve. The peak (e.g., Fig. 3) in the HVSR curve (Field and Jacob 1993;Lachet and Bard 1994) identifies the subsoil resonance frequency (RF) or, considering the reciprocal, the resonance period (RP). This approach assumes that the vertical component of the ground motion, differently from the horizontal ones, is unaffected by amplification phenomena related to the local stratigraphy. Consequently, the HVSR method uses the vertical velocimetric component as reference signal, avoiding the use of a second reference seismic station located on a site without local effects. The goal of the method is to detect the vibration modes and the dynamic amplification factors of one or more geologic layers through which the seismic waves propagate. The method proved to be an efficient tool for the detection of subsoil impedance contrasts (e.g., Lermo and Chavez-Garcia 1994;Konno and Omachi 1998;Bodin et al. 2001;Bonnefoy-Claudet et al. 2008;Boaga et al. 2013;Trevisani and Boaga 2018). In particular, the method is widely adopted in the context of seismic microzonation (e.g., Lermo et al. 1988;Chatelain et al. 2007;Gosar 2007;Bonnefoy-Claudet et al. 2009;D'Amico et al. 2008;Pilz et al. 2010;Foti et al. 2011;Teza et al. 2015) and for studying the thickness of sedimentary cover over bedrock (e.g., Ibs Von Sehts and Wohienberg 1999; Trevisani et al. 2017).
Briefly, considering a simple 1D layered stratigraphy (the theory can be extended to multiple layers), with a soft layer overlying a seismic bedrock, the subsoil RF (or RP −1 ) is controlled by the thickness (H) and the shear wave velocity (Vs) of the soft layer according to: Even if the equation is expressed in terms of shear wave velocity, it should be noted that the origin of HVSR peaks is still debated (Lermo and Chavez-Garcia 1994;Bonnefoy-Claudet et al. 2006;Castellaro and Mulargia, 2009), because it cannot be related purely to constructive interference of vertically bouncing of horizontal shear waves.
(1) RF = Vs∕4H Many authors suppose a contribution from surface wave ellipticity (Malischewsky and Scherbaum, 2004;Chatelain et al. 2007;Tuan et al. 2011;Boaga et al. 2013). Currently, the technique of single-station passive recordings and the HVSR analysis have been standardized (Sesame 2004). It should be noted that in the present study, the analysis is performed in terms of RP, both for spatial statistical motivations (see Trevisani et al. 2017) as well as because this notation is widely adopted in the context of seismic microzonation and structural engineering.

Passive seismic surveys in Kathmandu
The passive seismic surveys considered in this study have been designed to detect possible site amplifications with long resonance periods (e.g., Fig. 3), likely related to the impedance contrast of the bedrock-sedimentary cover transition. Given the expected high variability of bedrock top surface (Dill et al. 2001) and the heterogeneity of the sedimentary cover, some spatial variability of RPs is expected. As already mentioned, a portion of HVSR data have been collected during two surveys, one shortly after the Gorkha earthquake (November 2015) and another in November 2018, in the framework of a collaboration between NSET and OGS. The surveys have been designed to integrate the work of Paudyal et al., 2012 and, with long-duration records capable to detect long RPs. In fact, the results of the first campaign (Sandron et al. 2019) revealed that the RPs reported in the cited studies, considering also the stratigraphy of some boreholes (BH1, BH2, and BR, see Fig. 2) were referring to shallower reflectors and not to the deep reflectors, likely related to bedrock/sedimentary cover transition.
Some specifically selected data (8) from the study of Molnar et al. (2017) have been included, considering that the acquisition parameters were designed to detect lowfrequency resonance peaks. Then, a measurement located on rock (KTP, see Takai et al. 2016) has been included as reference of outcropping bedrock (Fig. 4).
Globally, 40 values of RPs (Table s1 in Supplementary Material) have been considered (Fig. 4), covering the western portion of the basin.

Geostatistical analysis of resonance periods
RP spatial statistical analysis has been performed following a standard geostatistical approach (e.g., Goovaerts 1997). Geostatistics is characterized by a well-known and welldeveloped theoretical framework and offers a rich suit of tools for the analysis and interpolation of spatial of data. The geostatistical approach is currently a standard tool for GIS specialists and quantitative earth scientists (e.g., Burrough and McDonnel 1998). A distinctive characteristic of geostatistics is the possibility to interpolate available spatial data objectively, i.e. taking into account their spatial continuity structure, and to evaluate prediction uncertainty (Isaaks and Srivastava 1989;Gandin 1965;Herzfeld 1996). It is worth to highlight that the geostatistical methodology is not a black-box approach; differently, a relevant part of the methodology relies on the Explorative Spatial Data Analysis (ESDA), including the study of spatial continuity. The ESDA is a fundamental step capable of shedding light on many aspects of studied phenomena, stimulating the users to find connections between computed spatial statistics and expert knowledge. From this perspective, the ESDA should be conducted not only when adopting geostatistical algorithms for interpolation but also when using alternative approaches, for example, based on Machine Learning (e.g., Kanevsky and Maignan 2004). In brief, a geostatistical analysis focused on interpolation is based on 3 main steps: (1) ESDA; (2) Inference of spatial continuity function (generally the variogram); (3) Interpolation and accuracy assessment via cross-validation. For the ESDA of RPs, statistical summaries and spatial statistical tools, such as classed post maps, (Isaaks and Srivastava 1989), variogram cloud (Ploner 1999), and directional variograms and variogram maps (Isaaks and Srivastava 1989), have been computed. The analysis has been carried out using the Gstat library (Pebesma 2004), implemented for the R statistical programming environment (R Core Team 2009), and some exploratory tools of the "Geostatistical analyst" package of the GIS software Arcmap 10.7 (Esri). To derive a preliminary mapping of the data, a Natural Neighbor (NN) interpolation approach (Sibson 1980;Sambridge et al. 1995;Watson 1999) has been adopted. NN is a convenient deterministic interpolation approach because of does not requires user-based decisions, being based on the Voronoi tessellation, and is capable to take into account the spatial clustering of the data. The main objectives of this step of the analysis are: the detection of local and global outliers; the detection of potential non-stationarity in spatial variability (e.g., presence of a trend); the evaluation of the necessity to transform the data due to highly skewed distributions (e.g., log and Box-Cox transformations); and the detection of anisotropies in spatial continuity. In particular, the variogram can be estimated correctly from the data only under, at least, intrinsic stationarity (stationarity of the mean and of the variogram). Consequently, during ESDA, it is crucial to detect departures from stationarity related to the presence of trends and non-stationarity of local variance.
The second step of the analysis is focused on the characterization and inference of the spatial continuity structure of RP. More frequently, the variogram (Eq. 2) is adopted as spatial continuity index, which is estimated using the experimental variogram estimator (Eq. 3).
The reported equations should be interpreted according to the geostatistical perspective, i.e. the available data are considered the outcome of a spatial random function. Consequently, the variogram ( ) (Eq. 2) is the semi-variance of the difference between two random variables Z(u) and Z(u+h) ; where is a vector defining a specific location in space and is a vector defining a given spatial separation. Under second order or, at least, intrinsic stationarity (Goovaerts 1997), the variogram can be estimated by the available data z . The experimental variogram is the half of the mean squared differences between the N( ) data pairs z − z( + ) spatially separated by the vector . Considering various discrete values of , it is possible to estimate the corresponding values of the theoretical function γ( ).
The experimental variogram estimated with Eq. 3 is unsuitable for interpolation for two main reasons: (1) a continuous function ( ) is required for determining the value of variogram for all possible values of ; and (2) the negative of the variogram should be a positive definite function (Goovaerts 1997). Accordingly, it is necessary to fit the experimental variogram values computed with Eq. 3 with a theoretical function, selected among a family of authorized models, by means of a least square or maximum likelihood approach (Goovaerts 1997;Pebesma 2004). In this study, the variogram models are fitted from experimental variograms by means of a weighted least squares approach, as implemented in Gstat software (Pebesma 2004). The weight given to a specific experimental variogram value is proportional to the related number of sample pairs and is inversely proportional to the squared distance.
The third step is focused on interpolation using Krigingbased algorithms (Goovaerts 1997). Kriging interpolators are members (Journel 1989) of a wide family of interpolators that are referred as Best, Linear, Unbiased, Estimators (BLUE) and are based on the normal equations of linear regression. The peculiarity of these linear interpolators is that the weights given to neighboring data values are defined with the objective of minimizing estimation variance, under the constraint of an unbiased error. The minimization of the variance under unbiased conditions is accomplished using the variogram or covariance function; from this perspective, the weighting scheme is objective, because it is defined by the spatial statistical structure of the data (Herzfeld 1996). The kriging algorithms have also two additional specific benefits: (1) estimation of prediction uncertainty; (2) possibility to interpolate on non-punctual spatial supports (block kriging), permitting, at the cost of reduced spatial resolution, a reduction of estimation variance (Goovaerts 1997;Isaaks and Srivastava 1989;Burrough and Mc Donnel 1998).
In this case study, given the spatial statistical characteristics of the data set, two kriging algorithms have been tested (Goovaerts 1997); one is the well-known ordinary kriging (OK) and the other is kriging with external drift (KD). In both cases, a block interpolation has been used, with block dimensions of 250 m, chosen as the balance between reducing estimation variance and resolution of the maps produced. The quality of the interpolated maps and the consistency of the inferred variogram model were evaluated by means of leave-one-out cross-validation (CV) diagnostics (Isaaks and Srivastava 1989), including the root mean square error (RMSE), the root mean square standardized error (RMSSE), and the statistical and spatial distribution of the errors. The CV consists of performing kriging interpolation in correspondence of sample data, removing temporary the sample value when performing the interpolation at specific point.
While the univariate interpolation approach OK is well known, the KD algorithm is less popular. KD is an interpolation approach that can be adopted when the primary variable of interest (e.g., RPs) is related to one or more secondary variables (e.g., the depth of the bedrock). This relation can be taken into account with specific algorithms (Goovaerts 1997;Hengl et al. 2007;Trevisani et al. 2017), such as co-kriging and kriging with external drift (KD), or its generalization, i.e. regression kriging (RK). KD and RK, which are mathematically equivalent (Hengl et al. 2007), with the exception of some differences at the operative level, are specifically useful when the secondary information is known exhaustively in the spatial domain of study. An example is represented by near-surface air temperature data correlated to elevation, the last known exhaustively by means of a digital terrain model (e.g., Hudson and Wackernagel 1994).
The KD/RK algorithms rely on the stochastic conceptualization, inherent to geostatistical theory, of the spatial properties of interest, according to which the random variable Z(u) can be decomposed in a residual stochastic component R(u) (second-order stationary) and a deterministic trend component m(u) (Eq. 4) R(u) is a stationary random function with zero mean and stationary variance and variogram; m(u) is a deterministic component which generally represents the regional spatial patterns. A first benefit of this representation is the possibility to handle the non-stationarity of Z(u) filtering out the trend component m(u) : this approach is at the base of OK and of Universal Kriging algorithms (Goovaerts 1997). A second benefit is that the trend component m(u) can be linearly related to secondary variables, known exhaustively in the study domain. In KD, m(u) is derived via a linear regression between the primary and a secondary variable (Goovaerts 1997); RK represents a generalization of KD, and the regression could be multiple (i.e., with more than one secondary variable) and could be potentially non-linear (Hengl et al. 2007). Both for KD and RK the regression of the trend component should be computed via Generalized Least Squares (GLS) using an interactive process, which takes into consideration the spatial correlation of data values, using the variogram of residuals R(u) . However, in many circumstances, given the difficulties in the implementation and possible issues of numerical instabilities, the regression can be computed according to a simplified approach based on Ordinary Least Squares (OLS). In general, the simplified approach based on OLS is acceptable in the presence of a large amount of data not excessively spatially clustered (Goovaerts 1997;Hengl et al. 2007). Both RK and KD, in spite of the operative differences, permit an estimate (i.e., obtain maps) of the primary variable of interest Z(u) , the trend m(u) and of the residuals R(u).

Results and discussion
The 40 RP data available have a clustered and elongated sampling geometry and are characterized by a high shortrange spatial variability (Fig. 5). The high short-range spatial variability is detectable both with the analysis of classed post maps as well as with the analysis of spatial continuity, using the variogram cloud (Fig. 6b). Even the preliminary interpolation with NN highlights abrupt variations in RPs (Fig. 7). The high spatial variability is compatible with the complexity of the geo-structural and geomorphological setting of the Katmandu valley and with the findings of preceding studies (Paudyal et al. 2012(Paudyal et al. , 2013. For example (Fig. 5), it is not rare to detect long RPs near bedrock outcrops, indicating that the bedrock top is characterized by abrupt changes in morphology (e.g., Gilder et al. 2020). Moreover, the high spatial variability can be partially related to the lithological heterogeneity of the sedimentary cover and hence to variations in the Vs (Sandron et al. 2019).
The variogram analysis (Fig. 6b-d) has been conducted considering a base lag of 1000 m, corresponding approximately to the median nearest neighbor spacing of the data. The few data available and the anisotropic spatial domain of analysis, elongated in the North-South direction, prevent a reliable modeling of potential anisotropies in spatial continuity. Notwithstanding, a higher spatial variability is detectable approximately in East-West direction, more evident when considering distances larger than 5000 m; unfortunately, in this direction, it is not possible to calculate the variogram for distances longer than 6000 m, differently from the North-South direction. The higher variability detected in this direction can be related to the presence of a trend, represented by a general tendency of increasing RP values moving from the mountains toward the center of the sedimentary basin, and, hence, from West to East (Fig. 5). However, it cannot be excluded, that a contribution to the anisotropy in spatial continuity is related to the residual component R(u) (Eq. 4), and hence to local patterns in the spatial distribution of RPs, such as the elongated pattern in NW-SE direction (Fig. 7).
To integrate the information available and to model quantitatively the trend, some secondary variables, potentially related to site-specific geological conditions and indirectly to RP, have been considered. Given the characteristics of the trend, a first secondary variable considered is the distance from the outcropping bedrock (Fig. 8a), due to a likely direct proportionality between bedrock depth and distance from outcropping bedrock (e.g., Trevisani et al. 2017). The simplified outlines of the outcropping bedrock have been digitized considering the geological maps presented in Sakai et al. (2015) and Shersta and Shah (2014). Only a weak linear relationship between RPs and distance from bedrock is detectable (Fig. 8b).
Then, potential correlations between RP and basic geomorphometric variables have been explored; the geomorphometric variables have been derived from the ALOS DSM (Jaxa 2018), which is a quasi-global DEM free of use and capable to describe surface morphology at a relatively high detail (e.g., Florinsky et al. 2019). Basic morphometric variables (e.g., Florinsky 2017) that could be potentially related to site-specific geological conditions have been considered ( Fig. 9): elevation, slope, total curvature and residual topography (e.g., Trevisani and Rocca 2015). The selected geomorphometric features have been calculated from a smoothed DTM, derived via a multi-pass moving window approach, to filter out fine-scale details and potential noise. The analysis of the scatter plots (Fig. 10) reveals a lack of evident correlations between secondary geomorphometric variables and RPs. The lack of simple correlations (e.g., linear or quadratic) between RPs and geomorphometric variables can be in part related to the fact that these variables are inherently related to shallow geological conditions. In this regard, the RPs considered in the analysis have been derived from HVSR recordings and analyses specifically designed for the detection of long RPs; accordingly, high-frequency resonance peaks mainly related to shallow reflectors have been filtered out. In addition, the relationships between geomorphometric variables and RPs could be characterized by complex non-linear patterns, not detectable with the few data available.
Given the spatial statistical characteristics of the data, two kriging interpolation approaches have been tested: (1) OK, using an isotropic model variogram inferred from the directional experimental variogram, computed in the direction of maximum continuity (North-South); (2) KD, exploiting the correlation between RPs and the distance from the outcropping bedrock and inferring an isotropic variogram from the residuals (Eq. 4). With both approaches, given the low spatial density of data, the interpolations have been performed on a relatively large spatial support, with blocks of 250 m × 250 m, subjectively selected to obtain a balance between interpolation accuracy and mapping resolution. The results of the two interpolation approaches, both in regard to spatial distribution of RPs as well as in relation Fig. 7 Natural neighbor preliminary interpolation of RPs, hatched areas highlight the outcropping bedrock. It is evident the elongated pattern in NW-SE direction, parallel to the alignment of isolated bedrock outcrops and to the Chobar fault (Fig. 2)   Fig. 8 a Computed distance from outcropping bedrock (hatched areas); b Regression of RPs versus distance from bedrock and related diagnostics to cross-validation diagnostics, have small differences (Table 1). For both, the RMSE is relatively high, the error is symmetrical and is not spatially correlated. KD resulted in a more realistic representation of RP values near the outcrops and with some improvements in cross-validation diagnostics; in particular, a slightly lower RMSE and a RMSSE nearer to 1, indicating a lower underestimation in modeling variance respect to OK. Given the better results of KD, the discussion will focus on the analysis and the results of this interpolation approach. The details and the results of OK interpolation are reported in the supplementary material (s2) accompanying this paper.
KD is capable to exploit the relationship between the trend and a secondary variable, here the distance from the outcropping bedrock. Given the few data available and the high spatial clustering of the data, an iterative GLS approach has been adopted for modeling the trend and the residual variogram. Four iterations were sufficient to stabilize the experimental residual variogram, which has been fitted with a spherical model variogram (Fig. 11). It is worth to highlight, that, even if the few data do not permit to model the anisotropy, from the residual experimental variogram, it is possible to detect a slight anisotropy, with the maximum continuity approximately in NW-SE direction. This direction corresponds to relevant structural lineaments crossing the Kathmandu basin (e.g., Figs. 2, 7 and 12). For the interpolation, a global search approach has been adopted, i.e., all the data are used for interpolating each grid node.
In Fig. 12 b, the areas characterized by interpolations variance higher than 0.5 s 2 (corresponding to a standard deviation of 0.7 s) have been masked, being characterized by interpolation artefacts and too high uncertainty. The selection of the estimation variance threshold for blanking is subjective and can be a choice related to the objectives of the study or user defined cost-functions (Goovaerts 1997).
Interpreting RPs as a proxy of bedrock depth, both the map of RP (Fig. 12b) as well as the low correlation between RP and the distance from the bedrock (Fig. 8) are compatible with the high variability of bedrock depth. The general tendency toward longer RPs moving toward the center of the basin is evident, generating an elongated pattern in North-South direction. The area of longer resonance periods, in the center of the basin, is compatible with the borehole BR (Fig. 2), reporting the weathered bedrock top at 550 m. Moreover, it is noticeable the southwestern portion of the area, where an abrupt variation in RPs is characterized by an elongated pattern in Northwest-Southeast direction, compatible with the directionality of structural lineaments of the area (Sakai et al. 2006(Sakai et al. , 2015. The elongated pattern Fig. 9 Basic geomorphometric variables considered in the study, hatched areas outline bedrock outcrops extends the alignment of rocky outcrops and it is likely connected to the Chobar fault (Fig. 2).
The map of RPs has been compared with the results of the work of Moribayashi and Maruo (1980), in which the depth to the bedrock is estimated by means of gravimetric method. Moribayashi and Maruo derived bedrock top isobaths (Fig. 5 in the cited paper) by means of gravimetric inversion of 112 gravimetric measurements, considering a density difference between bedrock and sediments of 0.8 g/cm 3 . Their results were compatible with the evidence of isolated bedrock outcrops and some deep boreholes. A limitation of their model relies on the use a homogenous value of density difference, that in some places seem to overestimate bedrock depth (Gilder et al. 2020). However, the spatial patterns represented are generally consistent with bedrock top morphology (Gilder et al. 2020). Consequently, it is worth to compare the results of two complementary geophysical methods both furnishing potential information on bedrock top morphology. The spatial distribution of RPs depicted in our map is compatible with the contour lines reported in the cited study (Fig. 13). Respect to the map of RPs of Paudyal et al. (2012 and, the RP mapping performed in the current analysis is characterized by much longer periods. In our study, where the estimated bedrock depth is more than 500 m, the RP block values are up to 3.7 s; in the cited studies, where the authors performed the interpolation performed on a punctual spatial support, are only up to 2.05 s. Accordingly it is evident that the maps of Paudyal et al. (2012 and, as confirmed by Sandron et al. (2019) and Gilder et al. (2020), outline shallower seismic reflectors. The analogies of spatial patterns derived from the gravimetric survey and the ones from HVSR mapping are promising. This suggests that a combined passive seismic and gravimetric approach could  be effective for reconstructing the bedrock top morphology of the Kathmandu basin. However, this can be achieved confidently only using high-quality stratigraphic and geomechanical data, including deep Vs profiles. In fact, part of the spatial variability of resonance periods and of gravimetric anomalies can be related to the heterogeneity of sedimentary cover and/or of the bedrock-sedimentary cover transition.

Conclusion
The geostatistical approach adopted, both with the spatialexplorative data analysis as well as with the mapping of resonance periods, furnishes useful information from the perspective of seismic-site effects related to long-period soil resonances. The analysis of spatial continuity is coherent with the available geological knowledge of the area. With this dataset, the use of geomorphometric secondary information seems not feasible to improve the interpolation results; only the distance from outcropping bedrock improves the geostatistical analysis. Despite the few data available, it is possible to map long-period resonance frequencies, even if for a limited extent and with a low resolution due to the high prediction uncertainty. The analysis of the produced map and the findings from the spatial data exploration reveal multiple interesting aspects.
• The produced map outlines evident sharp variations in RPs, with elongated spatial patterns aligned to structural lineaments, as in correspondence of the Chobar fault. • The high variability and the presence of anisotropic structures are also evidenced by the spatial continuity analysis. This suggest that even RPs related to deep reflectors, despite the larger volume of subsoil interested by the seismic elastic perturbation, can have an unexpected high spatial variability. • The RP map is compatible with the results of a preceding study based on a gravimetric survey and could represent a potential proxy of bedrock morphology.
These findings are promising but the available data limit the possibility to calibrate long resonance periods as predictor of bedrock depth. Given the spatial variability of RPs and the complex geo-structural setting of Kathmandu basin, it is first necessary a spatially dense set of passive seismic surveys, designed with long-duration recordings, to obtain a detailed and accurate mapping of long-period RPs. The spatial variability of RPs can be partially related also to spatial variations of Vs in the sedimentary cover; accordingly, in future studies, the analysis of Vs profiles and stratigraphies could improve the interpretation of the derived map. The combined use of HVSR and gravimetric methodologies seems a promising approach for improving the knowledge on the deep structure of the Kathmandu basin.
Funding Open access funding provided by Istituto Nazionale di Oceanografia e di Geofisica Sperimentale within the CRUI-CARE Agreement.
Code availability Not applicable.
Availability of data and material (data transparency) See data table in supplementary material.

Conflict of interest Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.