Ground motion parameters for the 2015 Nepal Earthquake and its aftershocks

The 25th April 2015 Nepal earthquake was the first major event in the Himalayan orogeny to provide a relatively well-recorded dataset. This paper presents a comprehensive analysis of the mainshock and its five major aftershocks through 21 well-established ground motion parameters. The analysis is presented for near-field stations of the Kathmandu basin and far-field stations of the Indo-Ganga Plains, including the site response behavior with varying sediment depths. Moreover, Hilbert–Huang Transform is used to study site amplifications associated with sediment depth and soil class, especially at long periods. In addition, a new ground motion model (GMM) is derived for all 21 parameters using moment magnitude, source-to-site distance, site class and sediment depth as predictor variables. Further, to ensure that the developed GMM is not biased, an inter-event and intra-event residual analysis is performed. Furthermore, the GMM is compared with other well-established prediction models applicable to the study region. The efficacy of the current model is further addressed by the evaluation of ranking indices.


Introduction
The Himalayan tectonic zone of the Indian subcontinent has been constantly subjected to incessant seismic activities, wherein a total of eight major earthquakes of M w > 7 have been recorded in the past two centuries. These include the 1866 Nepal (M w 7.4), 1897 Shillong (M w 8.1), the 1905 Kangra (M w 7.8), the 1934 Bihar-Nepal (M w 8.1), the 1947 Arunachal Pradesh (M w 7.9), the 1950 Arunachal Pradesh (M w 8.6), the 2005 Muzaffarpur (M w 7.6) and the 2015 Nepal (M w 7.9) earthquakes. The most common and unfortunate aspect of all these events is that none of these earthquakes were recorded well until the M w 7.9 Nepal earthquake. Thus, a thorough study of this 2015 event gained extreme importance. The 25th April 2015 Nepal earthquake originated near the Main Frontal Thrust (MFT) of the Himalayas, which is at the subduction plate boundary between the Indian and Eurasian plates. The rupture originated in the Gorkha district (28.147 °N and 84.708 °E), which is nearly 80 km northwest of the capital city-Kathmandu and propagated eastwards for about 160 km. Ground failure was reported in many densely populated areas, where the main event triggered landslides and ground liquefaction in the areas surrounding Kathmandu. There was unprecedented damage to buildings and heritage structures, where the economic loss is nearly 10 billion USD, which is about half the Gross Domestic Product (GDP) of Nepal and the human loss includes over 9000 fatalities. Several moderate to strong aftershocks followed this main event. According to USGS, 227 earthquakes of M w > 4.1 were recorded within a 250 km radius of the main event epicenter, in a span of just three weeks surrounding the main event. Among these, there are 30 aftershocks of M w > 5 and five aftershocks of M w > 6 due to the M w 7.9 rupture. In the immediate aftermath, many studies were conducted on source rupture characteristics (Hayes et al. 2015;Galetzka et al. 2015;Grandin et al. 2015) and artificial simulation of ground motions Dhanya et al. 2017). These studies emphasized the significant influence of rupture directivity and basin effect on ground motions. There were many studies based on the ground motion dataset as well. Goda et al. (2015) attributed the main cause of structural damage and collapse of buildings to the short period peak, which is due to ground vibration induced by near-source rupture. Dhakal et al. (2016) confirmed that the basin sediments at the KATNP station amplified the long period components of ground motions and concluded that the prediction equations of a similar tectonic environment such as Japan underestimate the recorded peak and spectral acceleration (S a ) values. The ground motion study conducted by Takai et al. (2016) with the five time histories (KATNP, PTN, THM, KTP and TVU) in the Kathmandu valley has also emphasized the amplification due to valley response at sediment sites. These studies, therefore, emphasize the effect of the sediment basin response in the Kathmandu valley.
A catastrophic event in the Himalayan region would also have a disastrous effect on the Indo-Gangetic Plains (IGP). Ground motions of the Nepal earthquake and its aftershocks recorded by the Central Indo-Gangetic Plains Network (CIGN) and the Indian Institute of Technology, Roorkee (IITR), networks have presented a rare opportunity to study the same. Due to substantial alluvial deposits, significant amplifications at the IGP sites were noticed previously (Srinagesh et al. 2011;Singh et al. 2020). Therefore, it would be interesting to study and compare the site response of the Kathmandu valley and the IGP. Chadha et al. (2016) and Rajaure et al. (2017) have obtained site amplifications for some of the IGP and Kathmandu stations for the Nepal earthquake using a Fourier-based approach. However, an advanced non-stationary approach, including the evolutionary characteristics of the ground motion, would yield better results (Zhang 2006). For this purpose, Ensemble Empirical Mode Decomposition (EEMD) with the Hilbert-Huang transformation (HHT) technique has been employed to include site nonlinearity. Further, Sreejaya et al. (2021) derived a nonparametric artificial neural network ground motion model to estimate a total of 21 parameters for the Himalayan region, using a dataset comprising 138 earthquakes, which includes the 2015 Nepal earthquake as well. However, the characteristics of the Nepal earthquake are not explored in detail, and the prediction model does not constitute sediment depth (D s ), which is deemed one of the key parameters affecting the ground motion records of this earthquake. Moreover, higher-order parameters such as the Root-Mean-Square acceleration (a RMS ), Cumulative Absolute Velocity (CAV), Characteristic Intensity (I c ), Arias Intensity (I a ), etc., and spectral parameters extracted from the Fourier spectrum, Power Spectral Density (PSD) and Evolutionary Power Spectral Density (EPSD) functions would extract important information from the ground motions, which were not explored thoroughly. Therefore, the present study aims at analyzing the ground motions of the 2015 Nepal earthquake and its five aftershocks using a total of 21 Ground Motion Parameters (GMPs) representing various amplitude, frequency and duration characteristics of the ground motion.
Many local and regional ground motion models (GMMs) are available for the Himalayan region of India (Sharma et al. 2009;Nath et al. 2005Nath et al. , 2012Iyengar and Kanth 2006;Kanth and Iyengar 2007;Joshi et al. 2013, etc.). Among these models, very few GMMs are derived using strong motion datasets (Sharma et al. 2009;Bajaj and Anbazhagan 2019;Sreejaya et al. 2021) and many of these GMMs are not applicable to the current earthquake and its aftershocks. Moreover, among the few available GMMs that are derived using the Nepal earthquake data (Singh et al. 2017), most were derived for standard GMPs such as Peak Ground Acceleration (PGA), Peak Ground Velocity (PGV), Peak Ground Displacement (PGD) and S a . There are fewer global models available for higher order parameters, such as CAV and I a (Campbell and Bozorgnia 2010;Foulser-Piggott and Goda 2014). However, since these models are derived for a different tectonic environment, the applicability as well as the goodness of these models to the Indian scenario is debatable. Hence, there is a necessity to generate a new set of GMMs, especially for the higher-order parameters, to verify the goodness of the local and global models that were applicable to this major earthquake. Therefore, this paper also aims at studying the behavior of these GMPs with source-to-site distance (R), D s and site characteristics. Accordingly, a novel regressionbased GMM for the Nepal earthquake is desired, and the attenuation relations should be validated through residual analysis. In addition, the hierarchy of the current model against the entire local as well as the global models should be validated through a ranking system.
In obtaining the above objectives, the rest of the manuscript is organized into six sections. In Sect. 2, the strong motion database used in the current study is discussed, along with the seismo-tectonics and geology of the region. Detailed analysis of the ground motions through the first-, second-and third-order GMPs is given in Sect. 3. In this section, amplification studies are also carried out using an HHT-based method. In Sect. 4, a region-specific GMM for all the GMPs is developed and validated through residual analysis. In Sect. 5, the developed GMM is compared with previously established models that are applicable to the region. Further, the efficiency of the current model is emphasized through the evaluation of ranking indices in Sect. 6. Finally, significant conclusions are provided in Sect. 7.

Strong motion data
Despite being seismically active, a nationwide seismic monitoring network is not available for Nepal. However, the earthquake and its aftershocks were recorded by four different strong motion networks: The United States Geological Survey (USGS), Hokkaido University, Japan, in collaboration with Tribhuvan University, Nepal (HU-TU), Central Indo-Gangetic Plains Network (CIGN) and the Indian Institute of Technology, Roorkee (IITR). The first set of ground motion records was provided by USGS at the Kantipath (KATNP) station, Kathmandu (USGS 2015). Along with the main event, this station has also recorded accelerograms of seven aftershocks. Later, HU-TU provided four strong motion records of the 2015 mainshock at stations KTP, PTN, THM and TVU. Center for Engineering Strong Motion Data (CESMD) made these acceleration time histories accessible (Takai et al. 2016). The relatively new CIGN network comprises 26 strong motion stations, 19 of which are located in the IGP (Chadha et al. 2016), whereas the IITR seismic network consists of around 100 accelerographs in the IGP (Kumar et al. 2012). However, due to technical issues, the mainshock of the 2015 earthquake was recorded by just 15 CIGN stations and 13 IITR stations. The number of recording stations is even lower for the succeeding aftershocks. Based on the records available, ground motions of the M w 7.9 2015 Nepal earthquake and its five major aftershocks are considered in the present study. Table 1 gives the details of these aftershocks, which range between M w 5.3 and M w 7.2 and were dated within a month of the main event. Figure 1 shows the epicenters of all the earthquakes considered in the current study, along with the distribution of ground motion recording stations.

Seismo-tectonics and site geology of the region
The Himalayan front is divided into four major tectonic divisions from south to north, which were demarcated by the presence of physiographic transition zones in the form of prominent faults (Gansser 1964;Valdiya 1964). The four tectonic divisions include-Sub Himalaya, Lesser Himalaya, High Himalaya and Tethyan Himalaya. The southernmost division, i.e., the Sub Himalaya, is separated from the IGP by the Himalayan Frontal Thrust (HFT), which is sometimes also referred to as the MFT. Other prominent faults from south to north include Main Boundary Thrust (MBT), Main Central Thrust (MCT), and the South Tibetan Detachment System (STDS). Nepal is diversified with all four tectonic divisions, wherein the Kathmandu valley is located in the Lesser Himalayan belt (Sakai et al. 2002). It can be observed from the seismo-tectonic map given in Fig. 1 that the earthquake rupture originated in the MCT, and the effect of these tremors was felt even in the IGP, though the impact of these tremors was severe in the Kathmandu valley. It is to be noted here that each of the tectonic divisions is characterized by similar physiography and the thickness of sediment deposits usually decreases from north to south. Therefore, D s is considered as an important attribute of the ground motions considered in this study. Therefore, the D s for the recording stations in IGP are obtained from the spatial distribution of subsurface bedrock of the IGP provided by CSIR-NGRI as presented by Sreejaya et al. (2022). Thus, D s of each ground motion recording station is obtained from 1D velocity models given by Bijukchhen et al. (2017) and Dhakal et al. (2016). In addition, site characteristics are assumed to be represented by the average shear wave velocity in the top 30 m of the soil strata (V s 30). The V s 30 values for the IGP sites are calculated from the SPT-N studies of Anbazhagan and Bajaj (2020), and the V s 30 values for the sites in the Kathmandu basin are obtained from Bijuchhen et al. (2017) and USGS. Table 2 gives the V s 30-based site classification adopted in this study, which is recommended by IBC (2009). The details of all the 40 ground motion recording stations, along with their respective D s and site class, are given in Table 3. In addition, finite source rupture models proposed by Hayes et al. (2017) are considered to calculate rupture distance, which is the shortest distance between the station and the fault plane. Further, in the case of a few aftershocks for which rupture models are not available, hypocentral distance is considered as the appropriate distance measure.

Preliminary studies
In the present study, 155 ground motions recorded at 40 stations are used for the analysis. Since it is established that the current study includes analysis of ground motions in the IGP as well, the range of distance measure is considered from 10 to 1000 km. Figure 2 shows the distribution of strong motion records used in the present study with respect to M w and R. It can be observed that the dataset consists of near-field records with R as low as 12 km and far-field records with R as high as 896 km. Histograms showing the number of records with respect to M w , R and hypocentral depth (d) are given in Fig. 3.
The raw data recorded by CIGN and IITR stations with a sampling rate of 100 Hz are subjected to a baseline correction of 0.1-25 Hz. On the other hand, HU-TU data with a sampling rate of 200 Hz are corrected from 0.1 to 49.5 Hz. Overall, the highest PGA of the dataset is observed to be 0.266 g at the near-field station KTP. Interestingly, the observed PGA values are unexpectedly low for such an earthquake which was notorious for severe    KTP, which is situated on a rock site, has similar (if not higher) PGA values to that of other near field records that are observed at soft soil sites. Therefore, PGA alone is not a reliable indicator of the ground motion record. Further, S a plots of the respective ground  Fig. 5, to observe any underlying pattern with respect to the period. In the figure, predominant long period amplitudes are observed clearly for near-field stations-TVU, PTN and THM, which are located on soft soil sites. On the other hand, long period amplitudes are relatively much lower for the near field rock site KTP as well as all the other far field records of IGP. Moreover, among the far field IGP records, the AMT record with a C type site class has slightly higher long period S a values to that of the other three records that were of D type site class. Hence, D s and site class are considered as critical factors to which ground motion amplifications of the current dataset can be attributed. However, such long period amplifications at sediment sites in the IGP cannot be extracted solely from the S a plots. Therefore, a detailed analysis of these ground motions is conducted using various GMPs.

Ground motion parameters
Ground Motion Parameters (GMPs) are necessary for defining the unique characteristics of ground motions in a concise, quantifiable manner. Thus, for seismic design and linear structural analysis, parameters representing peak values such as PGA, PGV and PGD are more preferred by the engineers. But, peak values solely are not sufficient to characterize damage potential completely. Therefore, many other GMPs like CAV, Vertical-to-Horizontal (V/H) spectral ratio, a RMS , I a , etc., were developed and were correlated with structural damage. Bhargavi and Raghukanth (2019a, b) have presented 21 such parameters and rated the damage potential of a ground motion using principal component analysis. The same 21 parameters are selected to analyze the ground motions considered in the current study. Amplitude, frequency content and duration are the three significant characteristics of a ground motion (Kramer 1996). Therefore, all the GMPs are divided into three categories based on their dependence on these key ground motion characteristics: first-order, secondorder, and third-order parameters. Table 4 explains the basic definitions of all the GMPs derived from acceleration time histories.
Correlation between time and frequency, (

First order parameters
First-order parameters represent only one key characteristic: amplitude or frequency content or duration. Among these, amplitude-based GMPs include PGA, PGV, PGD and V/H ratio of PGA. However, seismic damage is complicated to interpret solely based on peak values such as PGA, PGV and PGD of the ground motion. A time history with medium amplitude content that was recorded over longer duration or near the natural frequency of engineering structures might be more damaging than that of a transient time history with large PGA. Therefore, studying frequency content and duration characteristics of the ground motion is equally important. Thus, frequency content-based parameters are estimated through analysis of the time history in the frequency domain using Fourier Amplitude Spectra (FAS) and Power Spectra. Among these, the central frequency (Ω) and shape factor (q) are the frequency GMPs that were derived from the PSD function using spectral moments (Vanmarcke 1976). Similarly, the predominant frequency (F p ) is characterized as the frequency corresponding to the maximum amplitude of the FAS, which is also an important frequency-based parameter. Further, the PGV/PGA ratio reflects the frequency content characteristics of ground motion. On the other hand, the duration of the ground motion indicates the time needed for the fault rupture to release accumulated strain energy. The most popular duration-based GMP is the significant duration (T sig ), which is defined as the time interval between 5 and 95% energy thresholds of the Husid plot (Husid 1969). The first-order parameters that were derived for two near-field stations-'KTP' and 'THM' in the Kathmandu basin and one far field station-'AMT' in the IGP are tabulated in Table 5. Extremely low values of the peak GMPs, including PGV and PGD, are noticed for all the IGP records compared to that of the Kathmandu basin records. The F p values of records at near-field sediment sites are of the order 0.2, whereas the F p obtained at the station KTP is of the order 3.4. Thus, these low F p values emphasize the presence of predominant long period amplitudes in the ground motions recorded at sediment sites of the Kathmandu basin. On the other hand, even though predominant frequencies as high as 3 and 7.5 are recorded at some sites in the IGP, no clear pattern is observed with the variation of these values with the D s . Interestingly, the Ω of the vertical component of ground motion records is greater (double in most cases) than that of the horizontal component for both the Kathmandu and IGP. However, the V/H (PGA) values are less than 1 for most cases. Further, the parameter 'q', which represents the bandwidth of the time history, is almost the same in both the components (horizontal and vertical), if not slightly higher for the horizontal component of the ground motions. Therefore, based on the definition of 'Ω' and 'q', it can be assumed that the horizontal component of these ground motions has registered almost the same amount of energy at lower frequencies to that of the vertical component at high frequencies. As expected, significant durations of near-field records of the Kathmandu basin are extremely low compared to far-field IGP records.

Second-order parameters
First-order parameters are not sufficient for complete characterization of the ground motion. Housner and Jennings (1964) proposed a RMS as an indicator of damage as it considers the acceleration time history over a significant duration of the Husid energy plot. I a was introduced to observe the cumulative effects of ground motion over the complete duration (Arias 1970). Benioff (1934) proposed response spectra to quantify the physical response of a structure due to ground shaking. Park and Ang (1985) introduced a new secondary parameter-I c , in terms of a RMS and significant duration, which is linearly correlated to the index of structural damage due to maximum deformation and absorbed hysteretic energy (Ang 1990). Acceleration Spectrum Intensity (ASI) and Velocity Spectrum Intensity (VSI) were developed by Von Thun et al. (1988) and Housner (1952) to characterize ground motion at higher frequencies. These spectrum intensities are referred to the area under the acceleration and pseudo-velocity response spectrum for 5% damping ratio at all natural periods. Electrical Power Research Institute (EPRI 1988) mathematically defined a ground motion intensity measure named CAV as an efficient indicator of structural damage. Reed and Kassawara (1990) proposed CAV for determining the exceedance of the Operating Basis Earthquake (OBE) after the occurrence of an earthquake at a nuclear power plant.
In the preliminary analysis, it was noted that the maximum PGA of 0.266 g was observed at the station KTP, despite KTP being on a rock site. It should be noted that all four stations of the HU-TU network including the KTP station are located at almost similar source-to-site distances. Therefore, this preliminary result would be considered an inconsistency from the common notion that soft soil sites or sediment sites usually register larger ground motion amplitudes. Upon the analysis of the second-order parameters, it was revealed that the sediment sites indeed recorded larger amplitudes. The highest a RMS value of the dataset 0.055 g is obtained for the station TVU, whereas the a RMS value for the highest recorded PGA component of KTP is 0.033 g. In fact, all the other second-order amplitude parameters such as I a , CAV and I c are higher for TVU compared to that of KTP. This is due to the reason that the significant duration of the former ground motion is greater than the latter. Results like these emphasize the importance of the second-order parameters over the first-order parameters such as PGA, PGV and PGD. Table 5 shows the second-order parameters, i.e., a RMS, I a , CAV, I c , ASI, etc., for stations KTP, THM and AMT.

Third order parameters
An acceleration time history is highly stochastic and non-stationary. Thus, the statistical properties vary in both temporal and spectral domains simultaneously. Therefore, parameters reflecting the non-stationarity of ground motions are defined as the third-order parameters. To obtain these parameters, Evolutionary Power Spectral Density (EPSD) is generated through EEMD with HHT. Huang et al. (1998) proposed that the EMD-HHT is an effective technique to understand nonlinear and non-stationary data, due to the adaptive basis functions. Then, Wu and Huang (2009) modified the EMD as an Ensemble EMD, which includes the shifting of an ensemble of white noise modulated signals to remove signal intermittency. After that, the procedure has been employed in several studies to obtain the third-order parameters from the EPSD of ground motions (Raghukanth and Sangeetha 2012;Bhargavi and Raghukanth 2019a, b). In this method, the acceleration time history is decomposed into simple empirical orthogonal oscillatory modes called Intrinsic Mode Functions (IMFs). An IMF is a time series where the number of extrema and zero crossings must either be equal or differ by one. Further, the mean value of the envelopes defined by the local maxima and minima of the IMF should be zero. Initially, white noise of finite amplitude is added to the acceleration time history. From the resulting signal, the consecutive maxima and minima are identified. These maxima and minima are connected as cubic splines to obtain the upper and lower envelopes. Thus, IMFs are obtained by removing the local median in every step through a shifting process. This process is continued until a monotonic function is obtained at the end. Once the IMFs are obtained, the Hilbert Transform of these IMFs is used to construct EPSD. An Analytic Function ( AF j ) is formed from the conjugate pair of the obtained IMF ( IMF j ) and its respective Hilbert Transform IMF ( hIMF j ) to calculate the corresponding amplitude ( C j ) and instantaneous frequency ( IF j ) of EPSD. The amplitude of energy in time is combined over different frequencies for the extracted IMFs to obtain the energy distribution in the time-frequency domain in the form of a surface plot of EPSD (G( , t)). (1) The sample EPSD plots obtained for the near-field stations in Kathmandu valley (KTP, PTN, THM and TVU) as well as the far-field IGP stations (VNS, AMT, BSR and TDR) are given in Fig. 6. The energy distribution that was varying with time and frequency simultaneously can be visualized here, wherein high concentrations are observed over 6-8 Hz  Figure A5 and Vertical Figure A6) for the near-field stations. The statistical properties of the EPSD function are further characterized using six GMPs: Total Energy (E acc ) of the EPSD, first and second moments of the frequency and time axes: Spectral Centroid (E w ), Spectral standard deviation (S w ) and Temporal Centroid (E t ), Temporal standard deviation (S t ), and the Correlation coefficient of the time and frequency (ρ). All the 21 GMPs thus calculated for the three stations-KTP, THM and AMT-are provided in Table 5.

Analysis using IMFs and instantaneous frequencies
The IMFs of the acceleration time history (EW) for four stations in the IGP: VNS, AMT, BSR and TDR and four stations in the Kathmandu basin: KTP, TVU, PTN and THM, along with their respective percentage variances, are shown in Fig. 7. The original signal can be expressed as the sum of these IMFs. The percentage of variance of each IMF is considered as the contribution of each IMF to the total variance of the data. In most cases, the original data can be reconstructed using the first 10 IMFs. For the VNS station, it can be observed that the first four IMFs are sufficient to recreate 93% of the ground motion, wherein the first IMF represents 52% of the original data. For the station AMT, the first five IMFs are required to regenerate 94% of the original time history, with the first and second IMF having a percentage variance of 25%. For station BSR, 95% of the ground motion can be obtained by adding up the first six IMFs. In this case, the second IMF is more dominant than the first IMF, with 24% of the total variance. For station TDR, IMF4 is predominant with 19% of the variance of ground motion, and the first six IMFs are required to reproduce 93% of the ground motion. Thus, it can be assumed that as the thickness of sediments in the IGP increases from south to north, the IMFs with higher modes become more predominant.
In the case of the KTP record, which was registered at a hard rock site in the Kathmandu basin, the first five IMFs are essential for rebuilding the original ground motion. Here, the maximum percentage variance is observed for the first IMF, representing 38% of the original data. On the other hand, for the ground motions registered at sediment sites (PTN, THM, TVU and KATNP), higher IMFs (third to sixth IMF) are more dominant. For instance, IMF3 constitutes 30% of the variance for station TVU. On the other hand, IMF4 represents 49% of the ground motion for station THM and 30% of the ground motion for station PTN. Moreover, IMF6 constitute 54% of the variance for station KATNP. Therefore, it is inferred that IMFs with higher modes are indeed predominant for ground motions recorded on sediment sites of the Kathmandu valley as well.
The estimated IFs from IMFs of acceleration time histories for both IGP and Kathmandu basin are shown in Fig. 8. Table A1 shows the mean and standard deviation of the IF, along with the percentage variances of the first 10 IMFs for all the ground motion records of the M w 7.9 earthquake. Due to the inherent principle involved in the extraction of IMFs, IF values would be highest for the first IMF and decrease for subsequent IMFs. Moreover, the IF corresponding to the dominant IMF, i.e., the IMF with the maximum percentage of variance, is considered the dominant frequency of the acceleration time history. Consequently, it is observed that the dominant frequencies of the IGP stations-VNS, AMT, BSR and TDR-vary between 3 and 16 Hz, whereas the dominant frequency range for the Kathmandu stations-PTN, THM, TVU and KATNP-varies between 1 and 5 Hz. Therefore, the dominant mean IFs are decreasing as the R increases and the thickness of the sediment increases. In order to observe this property, the dominant IMFs and their corresponding mean IFs are plotted against D s in Fig. 9. Even though there is no underlying pattern in the dominant IMFs of the data at intermediate R, it can be observed that the higher modes are prevalent for larger sediment depths at larger distances as well as at near-field distances. Further, a  1 3 thickness (from north to south of the basin). Consequently, significant amplifications due to sediment depth are expected at the site's fundamental frequency. However, the GMPs that are obtained so far do not indicate any significant amplifications in the IGP. Therefore, in order to obtain amplifications at each of these sites, an HHT-based procedure is adopted.

HHT-based site amplification
In this section, site amplification analysis is carried out to provide a greater understanding of the effect of sediment depth on ground motions over short and long time periods. Traditionally, the surface site response of a site is compared with that of the bedrock site response to estimate the respective site amplifications (Shahri et al. 2013; Tempa et al. 2020). However, when only surface data are available, a reference site of hard rock strata (near the site) is selected to extract relative amplifications at all the other sites within the vicinity of this reference site. Based on this principle, various methods such as the Horizontal-to-Vertical Spectral Ratio (HVSR), the Standard Spectral Ratio (SSR) technique and the Ratio of Source Spectrum (RSS) technique have been used previously to estimate the site amplification in the IGP and Kathmandu basins (Chadha et al. 2016;Singh et al. 2020;Rajaure et al. 2017). However, these studies consider Fourier-based amplification, which assumes linear behavior of the site. In Fourier-based estimation, frequency is constant for the entire length of the data window, whereas HHT-based frequency content analysis  Figure A3 and Vertical Figure A4) considers temporal non-stationarity as well. Therefore, the effect of site nonlinearity on the ground motion time histories can be characterized by estimating site amplification based on the marginal spectrum obtained from the HHT (Zhang 2006). It is defined as the ratio of marginal Hilbert amplitude spectra of sediment sites to hard rock sites, similar to the SSR technique involving Fourier-based analysis. The marginal Hilbert amplitude spectra (A( )) over a temporal duration T of the ground motion and the HHT-based amplification factor ( FH s ( ) ) are given in Eqs. 5 and 6, respectively.
Here, G( , t) is the EPSD; the subscripts EW and NS represent the two horizontal components, and subscripts s and r represent the target sediment site and reference hard rock site, respectively. The HHT-based amplifications thus obtained at sediment sites of IGP and Kathmandu basin for the M w 7.9 mainshock are shown in Fig. 10. Site amplification factors are estimated for VNS, AMT, BSR and TDR stations in the IGP, considering hard reference sites as ALB, ALM, ALM and RPG, respectively. Similarly, the KTP station is considered as a reference hard rock site for TVU, PTN, THM and KATNP in the Kathmandu basin. In order to compare these values with the Fourier-based amplifications, the corresponding plots are obtained for the IGP and Kathmandu basin sites using the procedure given by Chadha et al. (2016) and Rajaure et al. (2017), respectively. The HHT-based dominant frequency is similar to the fundamental frequency obtained using the Fourier-based method for most of the sites in IGP but slightly diminuend at some stations in the Kathmandu basin. However, a significant increase in the amplification values can be observed at these frequencies for the HHT-based method than the Fourier-based method. The amplifications obtained for the four far-field stations VNS, AMT, BSR and TDR of the IGP for the M w 7.9 mainshock and the two aftershocks of M w 6.7 and M w 5.3 are shown in Fig. 11. It can be observed that the HHT-based dominant frequency is the same for M w 7.9 mainshock and M w 5.3 aftershock for IGP stations. The amplitude for the fifth aftershock (M w 5.3) is much lesser than the M w 7.9 mainshock and its immediate aftershock (M w 6.7) for most of the stations in IGP. This is a fair observation since aftershocks are usually created along a ruptured fault and hence have lower amplitude than the mainshock. The increased site amplification at lower frequencies for the higher magnitude events can be considered as an indication of soil nonlinearity at these sites. Remarkably, the station VNS with much lesser sediment depth compared to the other sites shows weak or no site nonlinearity under the mainshock. However, the dataset used in the present study is considerably small and sparse to obtain a relation between sediment depth and site amplification. In fact, the amplification factors cannot be found for many sites in the IGP due to the lack of information regarding a sufficiently close hard rock site. Therefore, given a comprehensive dataset, the nonlinear behavior and the effect of sediment depth on the amplification of the site can be studied more comprehensively for insightful observations. Overall, the HHT-based method proved to be a better option for obtaining site amplification.

Ground motion model
In this section, a region-specific GMM is developed. A suitable GMM is most important for predicting the future seismic potential of the region. Due to the lack of a comprehensive dataset, developing a regional GMM for Himalayan orogeny is challenging. Despite this, Singh et al. (2017) have developed empirical equations for the Himalayan region for three amplitude parameters: PGA, PGV and PGD. However, developing novel prediction equations for conventional parameters as well as those essential for engineering applications, like a RMS , I a , CAV, etc., would be invaluable. Due to the small and sparse dataset, the nonlinear mixed effect model (Abrahamson and Youngs 1992) is used to develop GMM in the present study to incorporate both fixed and random effects. The algorithm for the estimation of model parameters and variances involves estimating the parameters first using the fixed effects regression method. The most general framework of the GMM model is where Y ij is a GMP, ij = f M i , R ij , is the attenuation equation, M is the earthquake source measure, R is the distance measure, is the vector of model parameters, i represents the inter-event residual for an event i, and ij represents intra-event residual for the jth recording of event i.

Ground motion model for GMPs
The predictor variables considered for the GMM include M w , R , flags for V s 30 of the site ( F V s 30 ) and sediment depth variation ( F sed ). In order to formulate attenuation relationships Fig. 11 HHT-based site amplifications for stations VNS, AMT, BSR and TDR against the reference sites ALB, ALM, ALM and RPG of the IGP, for the M w 7.9 mainshock and the M w 6.7 M w 5.3) aftershocks between a GMP and the predictor variables, the following simple functional form is selected for all the 21 GMPs: Here, f mag represents the magnitude term where the coefficient c 0 corresponds to the magnitude or absolute scale of the GMP (Eq. 8). The term ′ c 1 M ′ represents the effect on the ground motion concerning earthquake magnitude considered as source measure. The distance measure f dis represents the term � c 2 R + c 3 ln R + c 4 exp c 5 M � that corresponds to geometric and anelastic attenuation (Eq. 11). Since all the high amplitude records of the dataset are concentrated at distances less than 20 km from the source, near-field magnitude saturation due to the M w 7.9 mainshock is expected. In order to capture this property, the term ′ c 4 exp c 5 M ′ is incorporated in the model, wherein these nonlinear coefficients are obtained using nonlinear least square regression analysis. The f site represents the term ′ c 6 F ′ V s 30 incorporating the local site conditions where ′ F ′ V s 30 includes flags for the site classes as per V s 30 (Eq. 12). Further, it is already mentioned in the earlier sections that the sediment depth varies from 0.6 to 5 km in the IGP, whereas the sediment thickness ranges between 0.04 and 0.46 km in the Kathmandu basin. Based on sediment thickness variation, the sites are classified into three different classes. The sediment depth flags ( F sed ) are considered as 0, 1 and 2 for sediment depth ranges D s ≤ 1, 1 < D s ≤ 2.5 and D s > 2.5 respectively. The event-to-event (inter-event, i ) and within-the-event (intra-event, ij ) residuals are independent, identically distributed and jointly normal (Eq. 14). The coefficients c k (k = 0, 1, … , 7) in Eq. 8 and the corresponding inter-event and intra-event standard deviations are shown in Table 6. The estimates thus obtained are plotted against recorded data for all the 21 GMPs in Fig. 12. It can be observed that the functional form considered in the study is able to capture the trend in the behavior of GMPs with source and distance measures. The inter-event and intra-event residuals from the mixed effect regression analysis are plotted against the predictor variables in the model to assess the validity of the GMM developed in this study. The variation of the residuals is shown in Figure A7 (Electronic supplement). In these residual plots, a positive residual indicates the underestimation of the actual value by the predicted model, and a negative residual indicates the overestimation of the actual value by the predicted model. It is observed that neither inter-event nor intra-event residuals show a trend with respect to the predictor variables.

Comparison with established GMMs
In order to study the goodness of the model, the GMM developed here is compared with previously established GMMs that were applicable to the region considered in this study 1 3 (Douglas 2020). Subsequently, a comparison of these already established models with the present model is presented for the two events of magnitudes M w 7.2 and M w 7.9.

Peak Ground Acceleration (PGA)
Regionally  the mean horizontal component of PGA and PGV for the same database. Similarly, CB14 was also developed for the NGA-West2 database, which includes period-dependent magnitude saturation, faulting effects, magnitude-dependent scaling with hypocentral depth and fault dip, regionally dependent anelastic attenuation, hanging wall effects and so on. K06 includes two regression models for shallow and deep events with a V s 30 of 300 m/s. Z06 includes attenuation models considering tectonic source types and focal mechanisms. SSSA17 was proposed for IGP based on an approximate solution of a circular finite source model using the same database. The GMM is derived using an exponential integral function where M w and R are considered as the predictor variables. Figure 13 shows the comparison of different GMMs with the PGA values estimated using a flag value of 2 for both F V s 30 and F sed and for two moment magnitudes: M w 7.2 and M w 7.9, for the current model. It is observed from the figure that the global models are under-predicting the recorded data.
Remarkably, the regional prediction model proposed by SSSA17 is also able to capture the trend of the data. However, it can be observed from Table 6 that the inter-event and intraevent standard deviation of the present model is lesser than that of SSSA17. Moreover, the present model also incorporates site classes and sediment depth variations along with the M w and R in the functional form. Thus, the GMM derived in the present study provides better estimates than other models.

Peak Ground Velocity (PGV)
Similarly  Figure 14 shows the comparison of these GMMs, with the PGV estimated using a flag value of 2 for both F V s 30 and F sed and for two moment magnitudes: M w 7.2 and M w 7.9. The models BSSA14, CB14, and CY14 are under-predicting the records, whereas estimates of SSSA17 and K06 are closer to that of the present model.

Peak Ground Displacement (PGD)
The models given by Singh et al. (2017)-SSSA17 andBozorgnia (2008)-CB08 are selected to compare the GMM derived for PGD. The model CB08 was derived for crustal events and includes hypocentral depth, shallow linear and nonlinear site responses, 3-D basin response and hanging wall effects. Due to a lack of information, the present model does not consider fault mechanisms and basin responses. Figure 15 shows the comparison plot of the GMM using a flag value of 2 for both F V s 30 and F sed and for magnitudes M w 7.2 and M w 7.9. The model CB08 is over-predicting the estimated values  Campbell and Bozorgnia (2008) of PGD, whereas the model SSSA17 is giving sufficiently closer estimates to the recorded data.

Arias intensity (I a )
Very few GMMs are available for complex parameters like I a , a RMS and CAV. The GMMs derived by Sarma and Srbulov (1998) The model SS98 considered duration and number of pulses as a function of acceleration and the usual predictor variables such as magnitude and distance. Further, surface wave magnitude (M s ) was used as a source measure, and equations were derived for rock and soil site conditions. Model A09 also used surface wave magnitude to derive four different equations for different ground types and tectonic conditions for Iran. Therefore, correlation relations provided by Scordilis (2006) have been used to obtain appropriate estimates. The comparison of these GMMs with the present model for I a using a flag value of 2 for both F V s 30 and F sed is shown in Fig. 16. It is observed that A09 estimates are drastically deviating from the attenuation trend with respect to R. Moreover, the model SS98 is over-predicting the data for M w 7.2 while able to follow the trend for M w 7.9. FPG14 is under-predicting the I a values for far-field records and overestimating the same for near-field records. As expected, the regression equation derived in the present study provides better estimates compared to other models.

Cumulative Absolute Velocity (CAV)
In this section, models proposed by Campbell and Bozorgnia (2010) (2014) magnitude, source-to-site distance, site conditions, faulting style, hanging wall effect, and linear and nonlinear site response. However, both models do not follow the trend in the behavior of parameters and give under-predicting estimates for CAV values. However, the model FPG14 provides sufficiently close estimates of CAV for far-field data while overestimating the near-field data. Figure 17 shows that the present model is relatively good at predicting the estimates.

Ranking of GMMs
Since the model derived in the current study used a limited dataset of ground motion records, the GMM cannot be used directly for any broader applications. However, the same model can be used in a logic tree approach with other regional models, especially in hazard studies, as the current model gives the best regional estimates. In this regard, there are two types of uncertainty in GMMs: aleatory variability and epistemic uncertainty. Aleatory variability is irreducible as it is associated with intrinsic variability. Epistemic uncertainty is attributed to incomplete or limited information and can be reduced with more pertinent data. The epistemic uncertainty can be handled by obtaining estimates from multiple prediction models, in a logic tree approach. In order to support the above statement, a detailed qualitative analysis needs to be conducted. Therefore, a statistical analysis is performed to facilitate the ranking of the established GMM hierarchically for future applications such as hazard and risk assessment. Many recent studies have proposed various methods for numerically scoring and weighing a set of GMMs, in order to identify the best GMM for a region (Mak et al. 2017;Roselli et al. 2016). The most popular methods for this task are the likelihood (LH) method (Scherbaum et al. 2004), the log-likelihood (LLH) method (Scherbaum et al. 2009) and the Euclidean distance (EDR) method (Kale and Akkar 2013). The basic concept in the LH method involves normalization of residuals obtained between the recorded and estimated data, with the standard deviation of the GMM. The LLH technique determines  (2013) the average log-likelihood of the observed and estimated data, where both are considered as continuous log-normal probability density functions. On the other hand, the EDR method makes use of the Euclidean distance concept to obtain the raking of various GMMs. All the above three methods may sometimes give different rankings for a set of GMMs, as the LH method is sensitive to standard deviations. Therefore, the concept of normalized RI N is opted here, wherein all the above three methods are used to derive an unbiased ranking system for the GMMs (Kale 2019). Table 7 gives the ranks thus obtained for all the GMMs considered for comparison of PGA in this study. It can be observed that all the GMMs are placed in a different ranking order for each of the three methods. Similar observations are noticed in the cases of PGV, PGD, I a and CAV, which are given in Tables 8, 9, 10 and 11, respectively. Therefore, the ranking indices are normalized for each method, with respect to the best performing GMM. Further, the final RI N is obtained by taking a weighted combination (EDR-0.50, LH-0.25 and LLH-0.25) of the corresponding normalized rankings (EDR N , LH N and LLH N ). In addition, the weights associated with each GMM are estimated from the normalized ranking (RI N ) of the GMMs using Eq. 15 to address the epistemic uncertainty in the GMMs.
The final ranks thus obtained for the GMMs of PGA, PGV, PGD, I a and CAV are given in Tables 12 and 13, respectively. From the tables, it can be observed that the (15) w j = 2 − log 2 (RIN) ∑ k j=1 2 − log 2 (RINj) model derived in the present study is ranked number one for all the GMPs, the reason being the calibration of the model exclusively to the recorded data used in this study.
On that note, it can also be observed that the other locally derived model SSSA17 is ranked number 2 for PGA and PGD. All the other global models register much higher RI N values for PGA, verifying that the Nepal earthquake is a truly unique event of the Himalayan tectonic zone, which cannot be estimated effectively by many well-established global models. Interestingly, the RI N values are much less for PGV and PGD, especially for the K06 model that was derived for Japan, which signifies that some of the global models are effective in predicting these parameters. The same pattern can be obtained from the plots provided in Figs. 13, 14 and 15. The rankings obtained for I a and CAV also follow the same pattern observed from the plots of Figs. 16 and 17.

Conclusion
The 2015 M w 7.9 Nepal earthquake was not only the most catastrophic but also the first well-recorded event to have occurred in the Himalayan orogenic belt. This paper presents a detailed analysis of this event using 21 GMPs, representing amplitude, frequency content, duration and evolutionary characteristics of the ground motion. Key results of the preliminary analysis include demonstration of the effectiveness of a higher-order parameters such as a RMS , compared to that of a first-order parameter such as PGA. Moreover, the analysis of S a emphasizes the presence of long period spectral amplitudes in near field ground motions of the Kathmandu basin, only at sediment sites. The behavior of the parameters concerning the thickness of sediment is studied through IMFs using the HHT-EMD technique. Sites with larger sediment depths have recorded low dominant mean IFs, due to a shift in predominant IMFs from lower to higher mode with an increase in sediment depth. The effect of sediment depth on ground motions is further analyzed through the HHT-based site amplification study, wherein higher amplification values are observed compared to the Fourier-based method. Comparing the M w 7.9 mainshock with two of its aftershocks, the increased site amplification at lower frequencies for the higher magnitude events can be considered as an indication of soil nonlinearity at these sites. Additionally, in this study, a GMM is derived for the 21 estimated GMPs, using nonlinear mixed effect regression analysis to incorporate both fixed and random effects. Among the few models applicable to parameters such as PGA, PGV, PGD, I a and CAV, superior prediction is observed from the estimates of the current model. A similar class of estimates can be expected for all the second-and third-order GMPs of the current model, which is a novel endeavor for the Nepal Himalayan region. In addition, the evaluation of RI N also highlights the efficiency of the current model. Since all the earthquakes considered in this study are originated due to a similar source mechanism, the GMM given here is not only region specific but also source specific. Analysis using a wide range of data would further ensure the goodness of this model and its applicability. The GMM developed in this study may not be used directly in applications such as seismic hazard analysis due to the constricted dataset used in deriving the models. However, when used with a logic tree approach, it will contribute quite significantly. Author contributions Jahnabi Basu contributed to formal analysis, investigation and writing-original draft preparation, Bhargavi Podili contributed to writing-review and editing and supervision, S. T. G. Raghukanth contributed to conceptualisation, methodology and supervision, D Srinagesh contributed to data provision. All authors read and approved the final manuscript.

Funding
The authors did not receive support from any organization for the submitted work.
Availability of data and materials Strong-motion data of the Gorkha earthquake sequence from the Central Indo-Gangetic Network (CIGN) are provided by D. Srinagesh at srinagesh@ngri.res.in. CESMD (Center for Engineering Strong Motion Data) provided ground motion records for stations KTP, PTN, TVU, and THM (Takai et al. 2016), whereas USGS provided ground motion records for station KATNP (USGS 2015).
Code availability Custom codes in MATLAB are used in the current study.

Conflict of interest
The authors have no conflicts of interest to disclose that are relevant to the content of this article.

Ethical approval
No human participants and/or animals are involved in the research.

Consent to participate Not applicable.
Consent to publication Not applicable.