Determination of path attenuation and site characteristics of the North-west Himalayan region and adjoining regions within the Indian Territory using Generalized inversion method

North-west Himalayas and its adjoining regions have been experiencing deadly earthqaukes from time to time and are home for a large portion of population of Indian subcontinent. Knowledge of regional path attenuation and site parameters are prerequisite while attempting seismic hazard studies towards minimizing damages during future earthqaukes for a region. Present work focuses on the determination of path attenuation and site characteristics of earthqaukes recording stations, located in the north-west Himalayas and its adjoining regions, within India. It is done using two-step generalized inversion technique. In the first step of inversion, non-parametric attenuation curves are developed by constraining attenuation to be a smooth decaying function with hypocentral distance. Q s = (105 ± 11)f (0.94 ± 0.08) as S wave quality factor is obtained indicating that the region is seismically active having high degree of heterogeneities in the crustal medium. In the second step of generalized inversion, site amplification curve, at each recording station, is computed as the ratio of site spectral amplitude of horizontal and vertical components. In addition, based on Horizontal to vertical spectral ratio (HVSR) method, predominant frequency of each recording station is calculated. Values of predominant frequencies based on HVSR and generalized inversion are found matching for each of the recording station. Based on obtained predominant frequency, site class of 101 recording stations, which at present are absent, are determined in this work. Determined path attenuation as well as site parameters can be collectively used for developing regional ground motion models and subsequently for seismic hazard studies for the selected region.


Introduction
The Himalayan arc, extending approximately 2500km between Kashmir and Arunachal Pradesh within India, is one amongst the most seismically active regions across the globe. High level of seismic activity of this region can be understood based on induced damages witnessed primarily during 4 great earthquakes (EQs) including 1897 Shillong EQ, 1905Kangra EQ, 1934Bihar-Nepal EQ and 1950 Assam EQ, which occurred in the last 120 years. Based on seismotectonic characterization, the entire Himalayan belt was subdivided into three distinct segments namely; the North-western, the Central and the Eastern Himalayas [Bungum et al., 2017]. The Himalayan region and its foothills, within India come under seismic zone IV and V as per IS 1893: [2016], indicating regions of high to very high seismicity. Therefore, the necessity of precise seismic hazard assessment of this region is of great importance.
Intensity of ground shaking during an EQ, at a particular site is a collective effect of source, path and site parameters.
Source parameters include magnitude, fault mechanism, stress drop and rupture process. On the other hand, path parameters account for the reduction in the energy of seismic waves (called attenuation of seismic waves) during its propagation in the medium. Attenuation of seismic waves is caused due to scattering of elastic waves in the heterogeneous medium and anelasticity of the earth medium, and is measured in terms of a dimensionless parameter known as quality factor (Lay and Wallace, 1995). Another important parameter which influences the nature of ground motion during an EQ is the site parameter. It accounts for the modification in amplitude, frequency content and duration of the incoming seismic wave by subsurface medium. The determination of aforementioned EQ parameters at regional level is important for the development of region specific synthetic ground motion models, which can further be used for region/site specific seismic hazard assessment [Baro et al., 2018]. Above EQ parameters are estimated in this work from regional ground motion records using an inversion approach based on generalized inversion method [e.g. Andrews, 1986;Castro et al., 1990;Oth et al., 2009].
In order to understand the ongoing seismicity of various regions within India, the Government of India had installed a number of EQ recording stations in different parts of the country. All these recording stations and corresponding ground motion records, since 2004, are handled by PESMOS (Program for Excellence in Strong Motion Studies). At present, PESMOS manages EQ records from 300 recording stations which are distributed in the northern and northeastern parts of India as well as in the Andaman and Nicobar Islands [Kumar et al., 2012]. Thus, PESMOS is considered as the most significant resource of ground motion records in India. Along with ground motion data, PESMOS provides information regarding magnitude and localization of EQs. However, it must be highlighted here that PESMOS database is lacking in terms of accurate information about subsurface for majority of recording stations [Kumar et al., 2012;Harinarayan and Kumar, 2018]. Site class (SC) given by PESMOS is based on physical description of surface materials and local geology following Seismotectonic Atlas of India (GIS 2000) and Geological Maps of Indian [GSI, 1998] and hence not based on actual field investigation (Kumar et al. 2012). In the absence of accurate information about local soil, utilizing EQ records from PESMOS database for region specific seismic studies is a major challenge. Subsurface exploration studies on some of the recording stations in north-west Himalayas by Pandey et al. [2016a; had highlighted the above limitation in SC given by PESMOS. There are recording stations which are classified to be located independent S-wave quality factor (Q s ) for Delhi region using the coda normalization method based on ground motion records from 9 recording stations. In another study, Negi et al. [2015] estimated Q s for the Garhwal Himalayas (western part of Uttarakhand) using extended coda normalization method based on 40 ground motion records from 8 recording stations. Using a similar method as done by Negi et al. [2015], Singh et al. [2012] estimated Q s for the Kumaun Himalayas (eastern part of Uttarakhand) based on ground motions from 23 EQ events from 9 recording stations. In another study, Mukhopadhyay et al. (2014) estimated Q s for the Uttarakhand region using Multiple Lapse Time Window Analysis (MLTWA) method considering ground motions from 30 EQ events from 5 recording stations. Results of the aforementioned studies will be discussed in conjunction with those from the present study.
Similar to path attenuation studies, very few studies on the determination of site characteristics from EQ records exist for this region. Nath et al. [2002] computed site terms using the aftershocks of the 1999 Chamoli EQ, obtained from 5 recording stations located in the Uttarakhand region. Similarly, Sharma et al. [2014] estimated site parameters for the Garhwal region of Uttarakhand using EQ records in context with generalized inversion and HVSR. In another work,  reported a comparative study on site characteristics that were computed using EQ records from Tarai region of Uttarakhand using multiple analytical approaches. Similarly, Harinarayan and Kumar [2018] computed site parameters for recording stations in the north-west Himalayas in terms of predominant frequency (f peak ) alone using Horizontal to vertical spectral ratio (HVSR) based on the response spectra of the entire accelerogram.

Study area
Present study area includes the states of Himachal Pradesh, Uttarakhand, Punjab, Haryana and New Delhi.
It covers an area between 28º N to 34º N latitude and 75.8º E to 80.5º E longitude. According to 2011 Census, this region has a population of 96 million. From seismicity point of view, major seismotectonic features of the study area consist of three north-dipping thrust systems such as the Main Central thrust (MCT), the Main Boundary thrust (MBT) and the Himalayan frontal thrust (HFT) [Valdiya, 1981]. The MCT and the MBT run parallel to each other within north-western Himalayan region and were produced during the Cenozonic shortening [Malik and Nakata, 2003]. The HFT is the youngest active thrust separating the Himalayan region and the Indo-Gangetic alluvial plain [Kumar et al., 2009]. The HFT, the MBT and the MCT so far have generated numerous major EQs in this region [Philip et al., 2014]. Other tectonic features include the Jhelum Balakot fault, the Drang thrust, the Lesser Himalayan Crystalline Nappes, the Jammu thrust, the Vaikrita thrust, the Karakoram fault, the Jwala Mukhi thrust, and the Ramgarh thrust. As a result of the presence of above mentioned seismic sources, the region has been experiencing repeated EQs. Two most damage inducing EQs in the region in the last 120 years include 1905 Kangra-Himachal Pradesh EQ (Ms=7.8) [Ambraseys and Douglas, 2004] and 2005 Muzzafarbad-Kashmir EQ (Mw=7.6) [Avouac et al., 2006]. Both of these EQs had caused large number of casualties as well as damage to property. 1905 Kangra EQ killed more than 20,000 people and caused 15cm uplift in Dehradun region located 250km from the epicentre [Ambraseys and Bilham, 2000]. The 2005 Muzzafarbad-Kashmir EQ claimed more than 80,000 lives and caused extensive damages to buildings in Jammu and Kashmir regions [Kamp et al., 2008]. The 2005 Muzzafarbad-Kashmir EQ also triggered several landslides along the Jhelum and Neelum valleys in Kashmir region [Kamp et al., 2008].

Nelliparanbil Hareeshkumar Harinarayan and Abhishek Kumar
4 Further, ground motion recordings are done in trigger mode during each EQ with a sampling rate of 200 per second.
For the present analysis, ground motion records of EQs happened between 2004 and 2017 in the earlier discussed region, as available on PESMOS are used. For estimating site characteristics, 341 ground motion records from 86 EQs, with magnitudes ranging from Mw=2.3 to 5.8 and having focal depths ranging from 2 to 80km are used. Further, these ground motion records were recorded at 101 recording stations, located in the hypocentral distance range of 9 to 355km. Coordinates of each of the recording stations, used in this work are listed in Table 1. Further, details of EQs used for estimating site characteristics are summarized in Table 2.
For estimating path attenuation however, only those EQs which are recorded in at least two recording stations within which at least one recording station is located within hypocentral distance equal to or less than a reference distance (section 4.2 for further details) are considered in the analysis. Thus, out of 341 ground motion records which are used for estimating site parameters as mentioned earlier, only 207 ground motion records satisfy the above mentioned reference distance criteria are used in determining path attenuation. Thus, the database for estimating path attenuation consists of 207 ground motion records from 32 EQs having magnitude in the range of 3.1 to 5.5(Mw), with focal depths in the range of 3 to 55km and hypocentral distance in the range of 9 to 200km, recorded at 69 recording stations. Table 3 summarizes the details of the dataset used for estimating path attenuation.
In addition, Figure 1 shows source-to-recording station path of the dataset used in the present study for estimating path attenuation.

Data processing
Before a ground motion record is used for path attenuation or site characteristics determination, it is processed.
All the EQ records are corrected by removing the baseline by a 5% cosine taper and a band-pass filtering using a Butterworth filter, between frequency range of 0.25Hz and 15Hz. Later, S-wave portion of the accelerogram is separated. The beginning of S-wave arrival is manually picked based on visual inspection. Time window duration of S-wave portion of the accelerogram is determined starting from 0.5s before the beginning of the S-wave and ending when 90% of the total seismic energy of the EQ record is reached [Bindi et al., 2009;Ameri et al., 2011]. Typical lengths of the S-wave time windows to be used for further analysis vary from 4 to 15s. For some of the records, where the S-wave window length obtained is longer than 15s, it is fixed to 15s in order to minimize coda wave energy in the analysing time window (Oth et al. 2008). Later, from the extracted windows, the Fourier amplitude spectra is calculated for each EQ record which is smoothened by applying the Konno and Ohmachi [1999] algorithm, with the smoothing parameter (b) as 20. Further, signal to pre-event noise (all having equal window length) ratio (SNR) for all the ground motion records is computed (similar to Wang et al. 2019). Figure

Path attenuation
In the first step of inversion, path attenuation curves (path attenuation in logarithmic scale versus hypocentral distance) are developed using a non-parametric approach where attenuation is constrained to decay smoothly with hypocentral distance. Detailed discussion can be found in following sub-sections.

Methodology
Following Castro et al. [1990], observed spectral amplitude [U ij (f,R ij )] of the horizontal portion of the accelerogram (obtained as the root mean square average of the east-west and north-south components), of EQ j at recording station "i", and at a particular frequency "f" can be modelled linearly as: Here, M j (f) is a scalar, which is governed by the magnitude of the EQ (one value for each EQ), A(f,R ij ) is the attenuation function which is independent of the magnitude of the EQ. In Eq. (1), A(f,R ij ) incorporates both geometric spreading and anelastic attenuation variation with hypocentral distance. Moreover, A(f,R ij ) in Eq. (1) is not limited to a particular functional form, instead, is assumed to decay smoothly with hypocentral distance (R ij ) and take the value A(f,R 0 ) = 1 at a reference distance (R 0 ) [Castro et al., 1990;1996;2003]. Representation given by Eq.
(1) has no factor representing site effect separately, and thus site effect is contained in both and. Any rapid undulations in A(f,R ij ) are assumed to be due to the absorbed site effects (Oth et al. 2008). Two weighing factors, W 1 and W 2 are incorporated in the Eq. (1) following Castro et al. [1990]. W 1 is used to smoothen the path attenuation curve by supressing the undulations and thereby transferring any absorbed site effects in A(f,R ij ) to M j (f). Further, W 2 is used to impose A(f,R 0 ) = 1 constraint, as mentioned earlier. The value of W 1 is chosen such that the site effects in the term A(f,R ij ) are supressed but the change in the attenuation characteristics with hypocentral distance can be clearly observed [Oth et al. 2008].
Eq. (1) can be expressed as a linear system of the form Ax = B, where B is the data vector having the term lnU ij (f,R ij ), x is the vector having model parameters [lnM j (f) and lnA(f,R ij )], and A is the system matrix that relates x and B. The matrix formulation of Eq. (1) after incorporating the weighting factors W 1 and W 2 (Castro et al. 1990) takes the form as follows: ( 2) The hypocentral distances of the data set are discretized into number of bins of equal widths and the value of A(f,R ij ) is computed for each bin. The width of the bins is selected such that there are sufficient numbers of data points

Nelliparanbil Hareeshkumar Harinarayan and Abhishek Kumar
in every bin. Figure 3 shows the number of EQ records in each bin of the present data set. It can be observed from Figure 3 that there are very less EQ records available beyond hypocentral distance of 115km. For this reason, EQ records with hypocentral distance up to 115km are only considered for the determination of path attenuation curve. This hypocentral distance range from 15 to 115km is divided into 11 bins (10 km wide). Further, constraint A(f,R 0 ) = 1 is applied at R 0 = 15 km [Bindi et al. 2004], irrespective of the frequency. Further, path attenuation curves are computed separately for each frequency, solving Eq. (2) in a least square sense, using singular value decomposition method [Menke, 1989]. For each bin, path attenuation curves are computed at 18 frequencies ranging from 0.5Hz to 15Hz (see Table 4).

Spectral attenuation with distance
Variation of lnA(f,R ij ) with hypocentral distance, obtained from the above analysis, for selected frequencies is depicted in Figure 4. Based on Figure 4, a general trend in which path attenuation curves exhibit decay with hypocentral distance up an hinge distance of 105km can be observed. The change in the attenuation rate in correspondence of the hinge distance can be clearly detected especially at lower frequencies (<5. However, detailed study in this direction is required and is beyond the scope of the present work. After observing the attenuation curves at different frequencies as given in Figure 4, it can be concluded that for the present study region, higher frequency content (>5.5Hz) of seismic wave, decays more rapidly between the source and the site in comparison to lower frequency content. This observation is consistent with the findings by Castro et al. [2003] for Guadeloupe (France) and Oth et al. [2011] for Japan.

7
Path and Site Characteristics In order to assess the stability of the inverted path attenuation curves, bootstrap analysis [similar to Wang et al., 2018] is performed at each of the selected frequency points. 17 number of ground motion recordings (which accounts for approximately 8% of the total number of the ground motion records considered for analysis) are removed randomly from the data base, and the remaining ground motion records are considered as a new data set used for inversion. The above discussed procedure is repeated 100 times and inversion (discussed in section 4.1) is carried out for each of the 100 bootstrap samples. Figure 5 shows the inverted path attenuation curves obtained from 100 bootstrap analysis. It can be observed from Figure 5 that the deviation from the path attenuation curve obtained considering the entire data set is significantly small, indicating that the path attenuation curves obtained from inversion are stable.

Nelliparanbil Hareeshkumar Harinarayan and Abhishek Kumar
8 Figure 4. S wave spectral attenuation versus hypocentral distance. Note that ln A(f,R 0 ) at reference distance is zero.

Quality factor (Q S ) estimation
In order to estimate, attenuation curves within the hypocentral distance in the range 15km to 105km, (where a monotonic decrease of ln A(f,R ij ) with hypocentral distance is observed in Figure 4 is considered. Each of the attenuation curves shown in Figure 4 is modelled in terms of geometric spreading [G(f,R ij )] and quality factor (Q s ) [as per Castro et al. 1996] as: ( 3) Where, f is the frequency under consideration and is the mean shear wave velocity in the crustal medium which is taken as 3.5km/s for the study area after Mukhopadhyay and Kayal, [2003]. It has to be noted here that a For each of the 17 selected frequencies (see Table 4), Eq. (3) is linearized by taking logarithm on both sides and corrected for the effect of G(R ij ) as given in Eq. (4). (4) Rearranging Eq. (4) gives: Ascribed to Castro et al. [2003], Eq. (5) is written in the form; a(R) = -m R Where a(R) and m are given as; Where, m in Eq. (6) represents the slope between a(R) and R which for each of the selected 17 frequencies, is calculated based on a linear least-square fit. Further, Q values are estimated for each of the selected frequencies by substituting the value of m in Eq. (8). Table 4 list the values of m and Q respectively. In order to develop the frequency dependent relationship of the form; Q = Q 0 f n , above estimated values of Q are fitted as a function of frequency using a power law. In this expression, Q 0 is the Q at 1Hz frequency and n is the frequency dependent coefficient, which is approximately equal to 1 and varies with heterogeneity of the medium [Aki, 1980]. Variation of Q versus f as illustrated in Figure 6 gives the following correlation for the northwest Himalayas as; Q = (105 ± 11) f (0.94 ±0.08) It has to be emphasised here that the Q estimated in the present study is based on the assumption of G(f,R ij ) = 1/R ij , and Eq. 9 is only valid under this assumption.
Values of Q 0 and (in the expression Q = Q 0 f n )are attributed to the level of tectonic activity and degree of medium heterogeneity respectively, present in the region. Aki, [1980] concluded higher values of for tectonically active regions in comparison to that of stable regions. Similarly, low value of Q 0 (<200) is an indication of larger degree of heterogeneities in the medium [Joshi, 2006]. The values of n (= 0.94) and Q 0 (= 105), obtained in this study indicate that the present study region is tectonically active, characterized by higher degree of heterogeneities. Previous study by Kumar et al. [2005] also reveals a higher degree of heterogeneity in the crustal medium of north-west Himalayas based on studying the lapse time dependence of the Coda Q (Q C ) in the frequency range 1.5Hz to 18 Hz.

Comparison with regional and global attenuation characteristics
As previously discussed, numerous studies exist where path attenuation of different parts of the present study area were attempted in the past. Comparison of present results with those obtained by the previous researchers for the north-west Himalayas and Delhi region is attempted as shown in Figure 7. It can be seen from Figure   Himalaya: Q S = 151f 0.84 [Negi et al. 2015]; Delhi and NCR region: Q S = 98f 1.07 [Sharma et al. 2015]. India: Q S = 71f 1.32 [Sharma at al., 2007]; Umbria-Marche Italy: Q S = 40f 1.33 ; South Central Alaska: Q S = 96f 1.06 [Dutta et al. 2004]; South Eastern Region of Tibet: Q S = 151.2f 1.06 ].

Site effects
As highlighted earlier, site characteristics of the recording stations are determined in the second step of inversion, which is discussed in this section. These are estimated based on both inversion and HVSR methods in this work and are compared. Detailed discussion on the GINV and HVSR is given in the subsequent sections.
Since then, various forms of this technique have been developed and used for estimating the seismic site characteristics by various researchers [Castro et al., 1990;Boatwright et al., 1991;Oth et al., 2008 etc.]. Methodology used in the present study for estimating site characteristics is as follows; As per Andrews, [1986] The path attenuation term can be removed from the spectral content of the record following Andrews, (1986) as; ( 11) The value of A(f) ij in Eq. (10) can be estimated using Eq. (3) and above calculated Q (refer to Eq. 9). Further, Eq.
(11) can be linearized [as per Andrews, 1986] by taking natural logarithms on both sides giving; ln U A (f) ij = ln S(f) j + ln G(f) i For a particular recording station, let: ln S(f) j = s j (f), ln G(f) = g(f) and ln U A (f) j = d j (f), Eq. (12) in accordance with Joshi et al. (2010) and following the notations of Menke, [1989] in the matrix form can be written as; (13) Eq. (13) represents a purely under-determinate system since there are (n + 1) × m unknowns for "m × n" data (here m is the number of sample frequency and is the number of EQs recorded at a particular recording station).

Nelliparanbil Hareeshkumar Harinarayan and Abhishek
In this work, Eq. (13) is solved using Moore-Penrose matrix inversion procedure (minimum norm inversion) given by Penrose, [1955] to determine g(f) for each of the recording stations.

HVSR
HVSR method is an extension of Nakamura, (1989) technique, which is widely used to assess the subsoil characteristics using recorded ambient noises. Nakamura, (1989) technique is based on the assumption that the soil amplification effects are retained only in the horizontal component whereas the source and the path effects are maintained both in vertical as well as horizontal components of ground motion. Hence, the ratio of horizontal and vertical components gives an estimate of site amplification. Lermo and Chavez-Garccia, [1993] extended Nakamura, (1989) technique to S-wave part of the accelerograms and studied the theoretical basis of the technique by numerical modelling of SV waves. Later, HVSR method was applied to EQ recordings worldwide [Luzi et al., 2011;Yaghmaei-Sabegh and Tsang, 2011;D'Alessandro et al., 2012;Kumar, 2017a, b, 2018 etc.] to obtain the site characteristics.
Comparative studies between HVSR and other methods of evaluating site parameters reported by Field and Jacob, [1995], Parolai et al. [2004], Shoji and Kamiyama [2002], Harinarayan and Kumar, [2017b] etc. show that, HVSR can provide good and reliable estimate of predominant frequency. However, above literatures also point out discrepancies in amplification levels obtained from HVSR and other methods. In order to compare the site amplification functions obtained from HVSR and GINV methods, same S-wave portion of the accelerogram is used in HVSR as well as in GINV in this paper. HVSR curve for each recording station is determined [similar to Harinarayan and Kumar, 2018] using the following steps; 1) Calculate the response spectra considering 5% damping, for all the three components (north-south, east-west and vertical) of ground motion records.
2) Obtain the geometric mean of the two horizontal response spectra components (H) using Eq. (14) given below; 3) Calculate the ratio of H to V (H/V).
Where, H EW and H NS are the response spectrum of the horizontal east-west and north-south components respectively and V is the response spectrum corresponding to vertical component of ground motion. Then, the HVSR at each of the recording station can be estimated as; (15) Here, N i is the number of events recorded at recording station "i" and (HVSR) i indicates the average HVSR value for a particular recording station "i". The f peak is the value of frequency corresponding to a maximum value of HVSR i (denoted by A peak ) at the recording station "i" which is also known as the predominant frequency of the recording station.

Site parameters
Site spectral amplitude (SSA) curves are developed using GINV for the horizontal (GINV H) and the vertical (GINV V) components separately. Figure 9 (a-f) shows typical amplification curves obtained for GINV H (indicated by dashed lines) and GINV V (indicated by firm lines) for 6 typical recording stations based on present analysis. In general, obtained amplification values from GINV H are higher than GINV V for all frequencies (see Figure 9a Identification of the value of f peak from GINV H/V curve is carried out in accordance with , who classified GINV H/V curves into three categories. In the first category, GINV H/V curves have only one peak and the frequency at the peak is classified as f peak . In the second category, GINV H/V curves have more than one peak and the smaller peak is classified as f peak . In the third category, GINV H/V curves appears flat, indicating peak is unidentifiable. Site classification is not attempted for recording stations having flat GINV H/V curves. Based on the above discussion f peak for all the recording stations are identified in this study. For each recording station, values of f peak obtained from GINV H/V and HVSR methods exhibit 1:1 (see Figure 11) matching (typical observations can be made from Figure 10a-i). However, there is difference in terms of A peak values obtained based on GINV H/V and HVSR. A peak values obtained using HVSR are found to be higher in comparison to those obtained using GINV H/V.
Similar observations were also reported for other regions [Sharma et al., 2014;Field and Jacob, 1995]. Values of f peak obtained using GINV H/V and HVSR are tabulated in Table 1 respectively. Similarly, the values of A peak obtained using GINV H/V and HVSR are tabulated in Table 1 respectively.

Nelliparanbil Hareeshkumar Harinarayan and Abhishek Kumar
Where, Vz is the shear wave velocity, H is the thickness which is considered as 30m in this study following NEHRP site classification scheme. The ranges of f peak for various SC are tabulated in Table 5. Based on the values of f peak obtained using H/V earlier for each of the 101 recording stations, SCs are proposed as listed in Table 1. Based on SC determination from this work, it can be summarized that out of total 101 recording stations, 2, 10, 29, 33, and 27 numbers of recording stations belong to NEHRP SC A, B, C, D and E respectively.

Empirical correlation between f peak and V S30
SCs for the recording stations, as discussed in section 5.3, are done assuming the depth of overburden to be 30m following NEHRP site classification scheme. In actual site however, this depth of overburden, whose response is observed in Figure 10 Table 6. These values of V S30 and corresponding f peak are presented in Figure 12. Another set of values of V S30 , assuming 30m as soil column depth, as obtained from Eq. (16), for same 27 recording stations are also presented in Figure 13 with corresponding f peak which is used earlier to determine SC (section 5.3). Comparison of two set of data in Figure 11 clearly indicates that V S30 values in absence of actual subsoil depth information, as obtained from can be used for SC determination for the present study area.
Another correlation based on V S30 as per Pandey et al., [2016a,b] and f peak from present work, as available for above mentioned 27 recording stations is attempted as shown in Figure 13. In absence of any functional form of correlation between the two in literature, following empirical correlation is obtained based on highest R 2 value as;
(a-f) Site spectral amplitude curves obtained using GINV for horizontal component and vertical component. log f peak = (1.5 ± 0.18)(logV S30 ) -(3.3 ± 0.5) It has to be highlighted here that Eq. (17) is applicable for sites having f peak in the range 1.8 to 6Hz and V S30 in the range 198m/s to 565m/s.

Conclusion
In the light of on-going seismicity and catastrophic damages witnessed during earlier EQs, of path attenuation as well as site characterization of PESMOS handled EQ recording stations, located in the north-west Himalayas and adjoining area of is attempted. Dataset for the analysis consists of ground motion records from EQs happened between 2004 and 2017. Entire analysis is done based on two-step inversion. While determining path attenuation in first step, a kink in path attenuation curve at 105km hypocentral distance is observed. Referring to similar observations from other regions, presence of Moho discontinuity is proposed in the region. However, such finding must be validated based on detailed study which is beyond the scope of present work's objective. Further, based on attenuation curve obtained till 105km hypocentral distance and over wide range of frequencies (0.05-15Hz), Q s = (105 ± 11)f (0.94 ± 0.08) is obtained for the present study area, clearly indicating that the region is heterogeneous and seismically active. It has to be highlighted here that dataset used in present analysis covers a much larger range of ground motion records in the analysis in comparison to earlier studies attempted for some part of present study area.
In absence of proper geological information of PESMOS recording stations as highlighted by numerous recent studies, predominant frequency of each of the 101 recording stations are calculated based on GINV as well as HVSR and are found matching for all the stations. Based on predominant frequency values for each recording station, SCs of 101 recording stations located within present study area, are proposed. Based on proposed site class, out of total 101 recording stations, 2, 10, 29, 33, and 27 numbers of recording stations are found belonging to NEHRP SC A, B, C, D and E respectively. It has to be highlighted here that unless the SC of recording station is known, ground motion records from a recording station cannot be used confidently for development of regional ground motion prediction equation or for region specific ground response analysis.
Both path attenuation and site characteristics, determined in the present work, are the key factors in developing regional ground motion models and are important inputs for seismic hazards studies. Without proper knowledge of above two parameters on regional level, suitability of developed ground motion prediction equation as well as seismic scenario to capture correct seismicity will be debatable.