Introduction

Earthquake hazard assessment is governed by three parameters such as source generation characteristics, wave propagation and ground response. To achieve this, analyses ought to be riveted on the analyses of earthquake history, earthquake source characteristics in the region and the attenuation relationships (^{Maiti et al., 2017}). Seismic hazard can be analyzed deterministically by constraining to a typical earthquake scenario Viz., the ever maximum earthquake hazard that has taken place in the terrain or probabilistically by taking into account the uncertainties involved in earthquake size, location, and time of occurrence, etc. over a considerable period of time to have data control and adequacy(Ministry of Earth Sciences, 2016).

A robust seismic design hinges on an understanding of critical parameters such as Peak Ground Acceleration (PGA), the Peak Ground Velocity (PGV), the Peak Ground Displacement (PGD) and the bracketed duration (^{Omidvar & Kivi, 2016}; ^{Takhirov et al., 2017}). Data of these parameters can be gathered from three types of strong-motion record such as seismogram (ground displacement), velocigram (ground velocity) and accelerograms (ground acceleration). The bracketed duration is the time duration of the ground shaking, defined as the elapsed time between the first and last acceleration excursions greater than some reference value (usually taken as 0.05 g, where g is the gravity acceleration). Recent studies have shown that the dominant periods of the Fourier Spectrum of the ground motion (displacement, velocity, and acceleration), are also very important parameters to be considered in the seismic design because ground shaking is more destructive on structures having a natural period around some of these dominant periods(^{Biasio, 2016}; ^{Agency, 2016}). Therefore, it is important to consider also the values of these dominant periods as important input parameters for the seismic design. It is essential to obtain relationships between the above-mentioned parameters (PGA, PGV, PGD, and duration) and the earthquake magnitude, to predict the possible values that these important parameters can take for future strong motions (^{Bommer, 2017}). Large-scale engineering projects, such as dams, tunnels, powerhouses and oil storage structures/ refineries etc., require estimates of the earthquake ground motions to which the construction may be subjected in a daily basis. Even though it is held that underground structures are significantly less vulnerable to earthquake shaking than structures on the surface, in view of the importance of underground structures and also it is reported that underground structures have suffered considerable damages due to earthquakes, a seismic design of structures underground structures is vital (^{Chian, 2016}) and therefore a detailed seismic hazard analysis is required and attempted in this work.

Relevant data sets have been gathered from the past earthquakes that have taken place in the country in the last 500 years. This also includes historical data prior to instrumentation and deductions made thereof. Time series methods use historical seismic data as the basis for estimating future seismicity and routinely updated on a regular schedule to incorporate additional data. Once the precursory signals are identified and studied then, modelling can identify the underlying factors that might influence the variable seismic signals and changes in the earthquake propagation zone. If the source parameters are understood, projections of the influencing variables can be made and used in the simulation; Probabilistic forecasting can be attempted.

Literature Review

The area of study falls in Northwest Himalayas, i.e. the highest seismic zone of the country (Zone - V; BIS) and to cite but a few of the important seismological studies carried out by various authors in the country and in various seismic zones, the following are listed:

The seismic hazard assessment studies for India have been carried out by ^{Baranwal et al. (2005}), ^{Chopra et al. (2010}, 2013), ^{Iyengar & Ghosh (2004}), ^{Jaiswal & Sinha (2007}), ^{Khattri et al. (1984}), ^{Kolathayar et al. (2012}), Kolathayar & Sitharam (2012), ^{Mhaske & Choudhury (2010}),^{Nath &Thingbaijam (2011}, 2012), NDMA (2011), ^{Parvez et al. (2003}), ^{Rao et al. (2002}, 2011), ^{Raghukanth (2011}), ^{Ravikumar & Bhatia, (1999}), Sitharam & Anbazhaga (2007), ^{Thaker et al. (2012}),^{Walling & Mohanty (2009}), ^{Yadav et al. (2008}). Each of these studies refers to specific regions and adopts a different methodology based on different set of assumptions to assess seismic hazard.

^{Ghosh et al. (2012}) aimed to study the current zoning maps of India. The Indian seismic code was prepared based on earthquake information available up to 1993 and was not updated since 2002. Accordingly, the study findings revealed that the Design Based Events (DBE) level seismic hazard is underestimated in many regions. The author states this finding is particularly important for the prevalent low-rise construction typically found in India. In future, the study has insisted on further collaborative research which is needed to update the seismic hazard map of India and to develop earthquake codes. The study presented by ^{Erdik et al. (2012}), homogeneously computed the seismic threats in the regions combining both local as well as the international expertise of data and provided the inputs for assessing and mitigating the risk of earthquakes in the areas of physical damage and socio-economic impacts. Similarly, the study by ^{Joshi et al. (2013})and Daniel et al. (2011) generated a seismic zonation map on the basis of deterministic modelling of the rupture plane of Uttarakhand Himalaya. The prepared map indicates that most of the area of Uttarakhand Himalaya falls under 10% or 20% exceedance in the probability of the PGA of 100 gals and 200 gals, correspondingly. This interprets that the vulnerability to seismic threats is highly prone in the regions of Uttarakhand Himalaya and its surroundings. Further, this study explains the need for seismic data in hazardous active regions and this study can also be used for the purpose of seismic zonation. Moreover, in the study done by ^{Sitharam & Kolathayar (2013}) presented the analysis of probabilistic seismic hazards in India by considering regional seismic zones and attenuation relations by taking into account of varies tectonic regions. The region under study was divided into minor grids with size 0.1. Additionally, for the periods of 0.1 and 1, Peak Horizontal Acceleration (PHA) as well as spectral acceleration was determined and also presented the contour maps displaying the spatial variation for the same. Through the presented study it can be interpreted that in the peninsular shield the vulnerability for seismic hazard is low, while it is high in major parts of the north as well as Northeast India. When the proposed study was compared with the recent study made by Iyengar (2010) in which the PHA values presented are greater in most parts except the regions with low seismic activity. Further, they suggested that strong motion records will aid in redefining of the seismogenic sources, which would be required for developing attenuation relations based on specific regions.

The research work by ^{Rout et al. (2015}) presented an improved technique for magnitude conversion for converting various magnitude scales to moment scales. In this study, the region of study and its surrounding areas was divided into seismogenic zones of 22 parts on the basis of seismicity, tectonics and geology including the source mechanism applicable to those regions. For the assessment of seismic hazard, attenuation equations which are region specific are utilized. Furthermore, the PSHA (Probabilistic Seismic Hazard Analysis) is conducted by adopting standard procedures, and from experimental results, it was found that for the probabilities of 2% and 10%, the values of PGA varied from 0.11 to 0.65g and 0.06 to 0.36g, respectively by taking into account of b-value. Moreover, greater PGA values were observed for the southeast region around Kaurik Fault System (KFS) as well as in the western side of Nepal. likewise, the study conducted by ^{Selva et al. (2016}) they developed a procedure for the quantification of uncertainties in Probabilistic Tsunami Hazard Analysis (PTHA), and additionally, special focus was given to uncertainties which are associated to numerical modelling in Seismic PTHA (SPTHA), and also on the process of subduction. Further to illustrate the procedure, they presented a case study by considering the sources in the Ionian Sea (central-eastern Mediterranean Sea), and by taking account of the coasts of Southern Italy as the target zone. From the experimental results, it can be interpreted that the quantification is efficient as well as complete even when considering a large number of sources and model formulations. Furthermore, the proposed model also allows the analysis of sensitivity and recombination which illustrates the applicability of the operational assessments.

The study conducted by ^{Ismail-Zadeh et al. (2017}) presented a model for block-and-fault dynamics (BAFD), which are utilized in simulating earthquakes due to the dynamics of the lithosphere and studies the effects of fault network properties as well as the regional movements of seismic patterns. Consequently, the performance analysis of the proposed model is done based on the generation of features of observed seismicity like earthquake clusters, the occurrence of events, fault slip rates, the relationship of frequency-magnitude and its mechanisms. Furthermore, they investigated new methods for the analysis of seismic hazard based on the historical, instrumentally recorded and BAFD-modeled earthquakes. Finally, the perceptions in the modelling of occurrences in earthquakes and the improvements in hazard assessment are discussed. The research done by ^{Galupino et al. (2017}) provided a review on the analysis of Probabilistic Seismic Hazard, and qualitative classification was done on the basis of the priority of low, medium and high Geographic Information System (GIS). Additionally, the authors guarantee that the developed structures resist a certain level of shaking of the ground, and also maintaining the required level of maintenance by the usage of provided maps. Moreover, this information is used as the basis for provisions related to seismicity in budget allocations in reducing risk, building codes and in risk models used in municipalities. In future works, suggested for respond hazard curve for the analysis of reliability for the periods like 0.2, 0.5, 3 or 4 secs.

From the review of the literature, seismological investigations recommended that concentration of stress and earthquakes are probable for occurrence when boundaries of three plates meet, like in NE-Himalayas and Kashmir. Also, it is the fact that major earthquakes are happening in these regions. Moreover, areas of lesser risks exist in the Himalayas in between Assam and Jammu as only two plates are interacting in the normal direction. Thus, attention should be focused on micro-seismic zonation.

Methodology

Seismotectonics of the Northwest Himalayas

The Himalayan region has undergone significant development and to ensure safe and secure progress in such a seismically vulnerable region, there is a need for hazard assessment. For seismic hazard assessment, it is important to assess the quality, consistency, and homogeneity of the seismicity data collected from different sources. The study pertains to a proposed road tunnel between Srinagar and Sonamarg that will pass through an assemblage of varied landforms such as gentle undulating plains, moderate to steep slopes of different mountain ranges, snow avalanche tracts, landslide zones, deeply incised gorges etc. with Sind River constituting the principal drainage. The area beyond Kangan receives heavy snow precipitation during the winter months that keeps the motor road presently blocked at several points for a considerable part of the year. Hung valleys south of Sonamarg and Zoji La pass are the main snow accumulation zones that severely impede the traffic movement every year. In order to considerably increase the period of smooth flow of traffic, a road tunnel designated as Z-Morh tunnel has been proposed in lieu of existing road which incidentally will avoid the heavy snow accumulation zones of Hung gorge section and Zozila. The tunnel will pierce through rock formations such as jointed slate, quartzite, and pyroclastics traversed by dolomite, shear zones and faults. There are several drainage crossings observed across the alignment, and some of the streams are perennial. Variation in lithologies, the presence of moderately to highly jointed rock mass, contacts of rock types with different permeability characteristics is likely to result in water ingress into the tunnel of varied quantum. Springs reported in the area confirm intersection of pheratic gradient.

Seismic performance of the underground tunnel as discussed by ^{Chen & Gui (2011}) states that seismic waves in underground tunnels are influenced by the stiffness of the tunnel and the geometrical layout .The tunnel lining and the soil surrounding structure also plays a significant part in determining the magnitude of the seismic wave. It is observed that deeper in the tunnel the lesser the earthquake impact on the tunnel lining. The soil pressure is redistributed after the earthquake which minimizes the bending moment of the tunnel lining. The seismic behaviour of the underground structures can also be seen in the study of tunnels by ^{Ertugrul (2010}). With large inertia of the ground, the seismic response is dominated by the surrounding soil of the underground structure. In underground structures, the deformation of the structure controls the seismic effect rather than the forces or stresses. So while designing the underground structures considering seismic hazard free-field deformation of the ground and its interaction with the structure are considered. This study mainly focuses on the seismic performance of the proposed tunnel which is assessed by investigating seismic vulnerability of underground tunnels. Seismic Atlas SEISAT (2000) has brought out a number of epicentres of Earthquakes of varying magnitudes close to the Main Boundary Thrust (MBT), Jhelum fault and Shinkiari fault which are in the vicinity to the project. As a part of the study, we have devised methods of using seismic databases to model and analyses. The data sets collected from ISS (International Seismological summary), BCIS (Bureau Central International Seismologic), EHB (Engdahl, van der Hilst, and Bulandi. e., Bulletin of the Seismological Society of America), CGS (California Geological Survey) and MOS (Marine Observation Satellite) are used as the database for conducting experiments in Matlab. Also, the data on earthquake was acquired from sources like the India Meteorological Department (IMD) from the year 1500 to 2012 was considered for the analysis. This type of analysisis unique because there is no need for any expert knowledge of system operation. But the software shows the propagation and distribution of seismic waves like time, magnitude and the energy associated with travel could be observed.

*Seismic hazard assessment*

There are two basic methods for assessing the seismic ground motion hazard in a particular region or at a specific site, namely deterministic methods and probabilistic methods. A full description of these methods is given in ^{Reiter (1992}). A Probabilistic Seismic Hazard Assessment (PSHA) combines seismic source zoning, earthquake recurrence and the ground motion attenuation to produce “hazard curves” in terms of the level of ground motion and an associated annual frequency of being exceeded. The elements used to carry out the PSHA in this study are as follows:

A definition of the seismotectonic source zones that define the geographical variation of earthquake activity.

A model of earthquake recurrence with respect to earthquake magnitude. There are generally more small (low-magnitude) earthquakes than large (high-magnitude) earthquakes.

A ground-motion prediction equation (GMPE), or attenuation relationship, is required, which defines what ground motion should be expected at location A due to an earthquake of known magnitude at location B.

The basic methodology adopted is based on that originally proposed by ^{Cornell (1968}), which includes integration over the aleatory variability of the ground-motion prediction equations. The probabilistic seismic hazard calculations were carried out using the Matlab software.

Ground motion database

In this study, horizontal components of the spectral acceleration values have been simulated for moment magnitude (Mw) ranging from 5 to 10 in increments of 0.8 units, at 30 values of hypo-central distances ranging from 1 to 200 km. To capture the finiteness of the source, the ground motions are also simulated for eight azimuths ranging from 0° to 270° in increments of 45°. Thus, a total number of 200 distance samples are considered for each magnitude. In all, there are 2000 pairs of magnitudes and distances. Since the stress drop, focal depth, dip, radiation coefficient, and pulsing percentage area are random variables, also the uncertainty aroused from these parameters are also included. The uncertainty associated with the local site amplification is included in the ground motion relations. The acceleration time histories are first simulated for class A site using the amplification function reported in (Working Committee of Experts (WCE), 2011). The simulated time histories combined with the developed shear wave velocity profiles are used as input to Matlab to obtain the surface level response spectrum. Thus, a database of 1000 PGA and Sa samples corresponding to 1600 simulated earthquake events are generated assuming 2% damping. This synthetic database is developed separately for both Kashmir and Himalayan regions with different site conditions.

Estimation of ground motion

The relationships that are predictable for the parameters like peak acceleration which decreases with increasing distance are generally named as attenuation relationships. For the utilization in engineering applications, PGA and its response spectra are required for the considered location. For the analysis of seismic hazard for a particular location and to construct a map for seismic zoning, the determination of attenuation and other scaling characteristics of several parameters with geological condition, distance and earthquake size are needed. The choice of suitable attenuation relationship influences the determination of ground motion for a particular location; hence the selection plays a significant aspect. Its accuracy depends on the data considered, function chosen and techniques used to determine it. Further, the relationships are classified based the zone and region considered.

Estimation of displacement Response Spectrum

A displacement response spectrum is the plot of maximum displacement of a single degree of freedom system (SDOF) to a particular ground motion as a function of the natural frequency and damping ratio of the SDOF. Derivation of the displacement response spectrum forms the basis for deriving other spectra. The time response of a linear SDOF system to a given ground acceleration is governed by the differential equation:

Where m is mass of SDOF, u is the displacement relative to the ground, üg is the ground acceleration; c,k are damping and stiffness of SDOF respectively. Further, the above equation can be written as follows:

Where ωn = 2Π/ Tn is the natural frequency, and Tn is the natural period of vibration. To get the response of SDOF system, the above equation is integrated with respect to time for a given damping, for a given seismic input once the natural period of vibration is known. If the response for different SDOF systems is calculated, maintaining the damping ratio and varying the natural period of vibration, and represent the maximum value of each response against the natural period of vibration, we get the Displacement Response Spectrum of the given seismic data for a particular damping ratio.

*Pseudo Velocity and Pseudo Acceleration Response Spectra*

The pseudo-velocity and pseudo-acceleration are given by:

The respective plots are called pseudo-velocity and pseudo-acceleration response spectra.

Results

In order to predict the impact of earthquakes as well as to perform waveform analysis, we have used the Matlab software. The data set is utilized for a period of 500 years from the year 1501 to 2000. The data sets have the values of Mw, Ms, and Mb in which the Mw represents the moment of magnitude, Mb represents the body wave magnitude, and Ms represents the surface wave motion. The mathematical modelling is done systematically with the help of Matlab without user intervention.

Estimation of seismic hazard parameters

Magnitudes of completeness (Mc), Gutenberg- Richter (G-R) recurrence parameter (b and a) values have been estimated for each seismogenic zone from the homogenized earthquake catalogue. The G-R relation (Gutenberg & Richter, 1944) is the relation between logarithms of the annual rate of exceedance versus earthquake magnitude and expressed as:

where λ_{m} is the mean annual rate of exceedance of magnitude m, a is the general level earthquake activity in a given area during the study period and depends on the seismicity of the region, b the slope describes the relative likelihood of large to small earthquakes(^{Rout et al., 2015}).

Estimation of PGA, PGV, PGD

In order to find the correlation coefficient, PGA, PGV, and PGD values are considered. The PGA value is obtained by calculating the b-value. The b-value is calculated with the help of dataset available in the U.S. Geological Survey (USGS) website from 1501 to 2012. The empirical formula for the b-value is:

The PGA value is estimated with the help of the following formula:

The value b1, b2, b3 is found with b-value formula, which is given in equation (6). In the derived equation (7), M is the magnitude; R is considered as a distance observed from a point to the earthquake focus.

The PGV value is estimated with the following formula:

In which Ts is the time series which is calculated using SD_{s} and SD_{1}

Whereas SD_{S} is the short period and the SD_{1} is the 1-second period in the considered time series T_{S}:

The site modified response spectrum (SM_{S}) can be obtained as,

Where, S_{S} and S_{1} are spectral response accelerations. F_{V} and F_{A} represent the function of site class when velocity and acceleration are considered respectively.

When PGA is known, PGD is calculated from:

Where α is the amplification factor that takes a value of 0.036, 0.054 and 0.065 for rock, stiff soil and soft soil sites respectively (^{Bommer et al., 2017}).

The following figures are plotted on the basis of the obtained Table 1 and the response spectrum of acceleration, velocity along with the input databases considered are plotted. In Table 1, only few values are considered due to page restrictions. Also, comparison with previous results is difficult, since different return periods are used.

Seismic input data acquisition consists of gathering and recording of continuous seismic signals from seismic stations. The accelerogram is plotted as shown in the above graph, with the use of seismic input data acquired, by considering time period (Sec) in x-axis and acceleration (g) in the y-axis. The obtained accelerogram shows that the input seismic data is high for a particular time period of 15secs and it slowly dampens.

The collected seismic input data is used to plot displacement in metres (in y-axis) with respect to the time period in secs (in x-axis). Displacement is a vector quantity that refers to "how far out of place an object is"; it is the object's overall change in position. The response of SDOF (single degree of freedom) system with 2% damping and natural period of 5s to strong ground motion is depicted. From the resultant graph, it is seen that the Peak Ground Displacement (PGD) is in the time interval of 25 to 35 secs.

Displacement response spectrum is the plot of maximum displacement of a single degree of freedom system (SDOF) to a particular ground motion as a function of the natural frequency and damping ratio of the SDOF. The displacement response spectrum is obtained by plotting the maximum displacement of the system along y-axis and time period along the x-axis. The obtained graph shows that the response spectrum rises in the time interval of 9 to 10 secs.

A response spectrum is merely a plot of the peak or steady-state response (displacement, velocity or acceleration) of a series of oscillators of varying natural frequency, which are forced into motion by the same base vibration or shock. In the above graph, displacement (m) Verses Time period (s) is plotted. Through resultant figure, it can be interpreted that the response increases exponentially.

The above graph shows the Pseudo Velocity Response Spectrum of Earthquake for 2% damping. Pseudo Velocity (m\s) is considered along with the y-axis and Time period (s) is considered along with the x-axis. The resultant graph shows that the response spectrum of velocity decreases in the time interval 2 to 5 secs and approximately remains constant over the time interval 7 to 10 secs.

The above graph shows the Pseudo Acceleration Response Spectrum of the seismic zone for 2% damping. The graph is obtained by plotting Pseudo Acceleration (m/s^{2}) versus Time period (s). The resultant figure shows that the response spectrum decreases exponentially and later reaches zero.

The above graph demonstrates the Acceleration Response Spectrum of Earthquake for 2% damping. The graph is plotted with Acceleration (m/s^{2}) along the y-axis and Time period (s) along the x-axis. The resultant graph shows that the response spectrum of acceleration decreases exponentially.

Discussion

In this, the software validation is done with the help of Matlab software, and the data validation is done with the PGA, PGV, and PGD value estimation. As it is observed that the seismic behaviour of the underground road tunnels is impacted by the ground motion, its parameters estimation and validation is done. ^{Trombetti et al. (2008}) discussed the characteristics of these parameters. Peak Ground Acceleration (PGA), Peak Ground Velocity (PGV) and Peak Ground Displacement (PGD) are both necessary and sufficient to exhaustively characterize the seismic input. The software validation is done with the help of the dataset considered. The collected data sets are analyzed using Matlab software and the metrics which include magnitude, acceleration, velocity, displacement, are predicted on an average basis and the predicted values are plotted in graphs (Figure 3, 4 and 5).

The geological and geotechnical site conditions greatly influence the strong ground motion at a site. The basement topography represented by the soil overburden thickness is connected to the site-specific hazard, especially if contrast exists in geophysical properties between the basement and the soil deposits. The sediment thickness implicates rebounding of the seismic waves leading to site amplification, and therefore, has direct implications. Likewise, undulating topography (elevations) would produce scattering, focusing, or defocusing of incident seismic waves. Observed data contain in database satisfy the following conditions that hypocentral distance is ranged from 5 to 100 km and AVS30 (Average shear-wave velocity in the upper 30m) of observed sites is ranged from 600m/s. Another hazard factor is the soil liquefaction, which is connected to loose soil and geo-hydrological conditions, and therefore, is a determinant geotechnical hazard, especially at the reclaimed sites previously of natural water bodies, and swampy tracts. In hilly terrains, on the other hand, seismic hazards are a foremost determinant induced hazard. The underlying primary hazard is generally represented by peak ground acceleration, a short period ground motion parameter signifying damage potential to the buildings enabling an overall quantitative basis for the design codes and construction practices. However, when the analysis is targeted at specific buildings, period-specific spectral acceleration would also be a better index to account for applicable resonance frequencies. The typical case examples in this paper employ attenuated frequency in the multi-criteria assessment to derive generic hazard conditions. Period-specific analyses based on peak ground accelerations are envisaged for building typological specific and ward-wise studies. The amalgamation of the different hazard factors is aimed at the better representation of the local specific seismic hazard variation in the study regions.

Conclusions

In this study we have collected and analysed seismic data from 1501 to 2012 (500 years) from India Meteorological Department(IMD) and also utilized the databases of ISS (International Seismological summary), BCIS (Bureau Central International Seismologic), EHB (Engdahl, van der Hilst, and Bulandi.e., Bulletin of the Seismological Society of America), CGS (California Geological Survey) and MOS (Marine Observation Satellite)e and have validated the experiments by the estimation of PGA, PGV and PGD values and applied the same for a proposed road tunnel in the Kashmir Himalayas (GPS coordinates of Latitude 32.20° to 37.20° N, Longitude 72.50° to 80.50° E). The experimental results and the graphs plotted established that the datasets used in the evaluation for Kashmir Himalayas is reliable and matched with that of data set of USGS (National Earthquake Hazards Reduction Program, 2018). The prediction of the earthquake and its seismic intensity area is made with the values of different probability. This probabilistic approach used the Matlab software which predicted the Velocity, time and energy associated with seismic wave propagation and also proved them with dataset model analysis carried out. The PGA, PGV and PGD studies were correlated based on these datasets and this may be used in preparing the Seismic hazard map as well as building code design as a source. Experimental modelling of scale models with PGA, PGV and PGD may provide further opportunities for refining studies of Seismic Hazard Microzonation and therefore updating Seismic Zonation of the Country.