Introduction

The Jammu and Kashmir region in northern India has witnessed several natural disasters like earthquakes, floods, landslides, and avalanches (Sana 2019; Ansari et al. 2021a, b). The continuing collision of the Indian and Eurasian plates has resulted in earthquakes of varying magnitudes in the Himalayan region including the deadliest Muzaffarabad earthquake (Mw 7.6) in 2005 and the most recent Mirpur earthquake (Mw 5.6) in 2019 (Ansari et al. 2022d; Gupta and Satyam 2022). The 2005 Muzaffarabad earthquake resulted in severe damage in Muzaffarabad, Uri, Kupwara, and Balakot causing the deaths of 72,763 people and injuring a further one lakh people (Sana and Nath 2016; Ansari et al. 2017; Sana 2019; Sana et al. 2019; Powali et al. 2020). The tremors of the 2019 Mirpur earthquake in Pakistan were felt in the Jammu Region (JR), Punjab region, Uttarakhand, and Delhi causing the deaths of 50 people and property damage in small towns of Jatlan, Manda, and Afzalpur near Mirpur city. The slope failure and ground shaking were observed in the Rajouri and Poonch districts due to the 2019 Mirpur earthquake. The post-effect of these earthquake events created the problem of liquefaction in different parts of Jammu and Kashmir (Ansari et al. 2022b). The consequently induced acceleration and internal strains may generate higher in situ shear stress along a potential failure line during earthquakes (Papathanassiou 2008; Zhan and Kanamori 2016). Furthermore, as a result of the ground shaking, the material may tend to densify, resulting in increased pore pressure in saturated loose soil and, as a result, a reduction in effective stresses. Liquefaction of soils occurs when a significant extra pore pressure builds up due to partially undrained deformation of loose soils, resulting in a lower undrained shear strength than drained shear strength. In general, each site’s liquefaction threat is assessed using a factor of safety (\({F}_{s}\)) based on the cyclic stress and resistance concept, and liquefaction potential index (\(\mathrm{LPI}\)). The sites with \({F}_{s}\) < 1 are regarded as liquefiable, while those with \({F}_{s}\) > 1 are considered as non-liquefiable. \({F}_{s}\) = 1 denotes the limiting equilibrium state.

Over the last few years, construction of mega projects including the world’s tallest Chenab bridge, India’s first cable-stayed Anji Khad bridge, Pir Panjal tunnel, and Udhampur Srinagar Baramulla Rail Link Project (USBRL) have progressed at a faster rate in the JR (Ansari et al. 2022a, b, c, d, e). The past historical records of the earthquake events and rapid infrastructural development push to assess the liquefaction hazard of the study area. In the present study, an attempt has been made to provide the liquefaction hazard maps in terms of \({F}_{s}\) against liquefaction, liquefaction potential index (\(\mathrm{LPI}\)), and probability of liquefaction (\({P}_{L}\)) at different depths for the JR based on geotechnical and geophysical approaches. The standard penetration test (SPT) and cone penetration test (CPT) are the most common method used to calculate cyclic resistance. The Seismic Dilatometer Marchetti Test (SDMT) is increasingly becoming more widely used in normal geotechnical investigations, allowing for the collection of a large number of data points. The SDMT calculates a horizontal stress index that has a strong link to soil liquefaction (Grasso and Maugeri 2008; Grasso et al. 2020).

The liquefaction hazard assessment is carried out using both SPT (Boulanger and Idriss 2014) and shear wave velocity (\({V}_{s}\))-based methods (Andrus and Stokoe 2000; Andrus et al. 2004). Extensive borehole data have been collected from geotechnical consultancies for the study area. The geophysical tests were conducted on natural land to examine the subsurface data for shear wave velocity at different depths based on the multichannel analysis of surface wave (MASW) technique. The ground rupture was observed during the field investigation in the Simbal village in the JR. Local people of Jatah (Samba district) observed sand blows near Degh Nala after the 2019 Mirpur earthquake. The liquefaction potential of the study area was assessed using the PGA at the surface level from seismic response analysis in DEEPSOIL (Hashash et al. 2008), geotechnical data from consultancies, and geophysical data from field tests at 243 site locations. Based on the possibility of liquefaction, the study area is classified into different zones using geotechnical, geophysical, and seismic response data.

The areas in the southern part of Jammu, Samba, and Kathua have been identified as liquefiable due to the existence of silty sand and fine sand, as well as a shallow water table. The sites in the northern districts like Ramban and Poonch, and the northeastern part of Doda and Kishtwar are not susceptible to liquefaction. An integrated liquefaction hazard map for the study area is also generated using the findings of both SPT- and \({V}_{s}\)-based techniques. It has been also reported that a few locations that are prone to liquefaction using a \({V}_{s}\)-based approach is not susceptible to liquefaction using SPT-based approach. The maps for liquefaction potential index (\(\mathrm{LPI}\)) and probability of liquefaction (\({P}_{L}\)) is also presented, which gives a better delineation of hazards in the study area. This data may be used as a tool for engineers to make risk-based design decisions for liquefaction mitigation in the study area.

Tectonics and local geology of the study area

The JR in the northwestern Himalayas is limited on the north by Kashmir valley and on the east by Ladakh, where the southwestern part has flat terrain, and the rest is made up of Jammu hills (Ansari et al. 2022a). The height of the Jammu hills ranges from 330 m in the Jammu city to 1638 m in the Kishtwar on the eastern side, which touches Ladakh on the India-China border (Paul et al. 2018; Pandita et al. 2019). The groundwater table depth varies between 2 and 22 m, with north to south flow, due to variations in subsurface condition and topographical changes from Reasi in the north to Kathua in the south (Bhanja et al. 2018; CGWB 2020; Qadir et al. 2020). The Chenab, Tawi, and Ravi are the main rivers flowing in the JR, which further enters the territory of Pakistan. The Chandra and Bhaga rivers, which meet near Keylong in Himachal Pradesh to form the Chandra-Bhaga or Chenab River, are two important tributaries of the Chenab River. After passing through the Kashmir Himalayas, the river emerges into the plains near Akhnoor in Jammu. The southern portions of Samba and Kathua, near Pathankot in Punjab fall under the basin of the Ravi River (Razzia and Dar 2021). The Tawi River passes through the Doda, Udhampur, and Jammu and then merges with the Chenab River near Sialkot in Pakistan.

The Jammu, Samba, Kathua, and southwestern part of Udhampur districts lie at the foothills of Siwalik, built over deep sedimentary basins made up of unconsolidated quaternary sediments of the Jammu formation, older and younger alluvium (Ansari et al. 2021a). Subsurface strata from Jammu to Reasi show Dolomite formation and Murre formation (Tiwari and Latha 2020). The northern part of Jammu city comprises of Nagrota formation while the southern part consists of younger alluvium near the banks of the Tawi River. On the western side of Udhampur, clast-supported deposits can be found at the top of the fan, sandy deposits can be found in the central area, and clay deposits within gravel bands can be found at the bottom of the fan (Mahajan et al. 2012). Because of its extensively weathered and fractured character, dolomite rock formation regions are the most predominant among the primary rock types found in the study area. The thickness of dolomite ranges from 150 to 200 m (Srinivas et al. 2014). The geographical location of the study area along with lithological formations is presented in Fig. 1.

Fig. 1
figure 1

Lithology map of the Jammu Region with benchmark locations and sites where liquefaction was observed during far-field earthquakes

Full size image

The JR has been impacted due to the near-field as well as far-field earthquakes in the Himalayan region extending up to Hindukush in Afghanistan (Lister et al. 2008; Asim et al. 2017). The major earthquake events in the northwestern Himalayas are listed in Table 1. The following Fig. 2 presents the total number of earthquake events that occurred in and around the JR for different magnitude classes in the last 500 years. The Jhelum Fault (JF), Attock Fault (AF), Reasi Fault (RF), Balakot-Bagh Fault (BBF), Deosai Fault (DF), Jwalamukhi Thrust Jwalamukhi Thrust Hanna Fault (HF), Batal Fault (BF), and Mawer Fault (MF) are few of the active faults that surround this area (Bilham and Wallace 2005). The Main Central Thrust (MCT) distinguishes the crystalline rocks of the higher Himalayas from the formations of the lower Himalayas (Gupta and Gahalaut 2014). Along the JR’s northern boundary, the Main Boundary Thrust (MBT) and Panjal Thrust (PT) run parallel. NNW-SSE and NW–SE trends are shown by the Jwalamukhi Thrust (JT) and Balapur Thrust (BT) respectively (Malik and Mohanty 2007; Alam et al. 2015; Sana and Nath 2017). The NE-dipping Kishtwar Fault (KF) and the NS-trending Jhelum Fault (JF) are the two major local strike-slip faults in the JR. The BBF, which is NE dipping in Pakistan, is the primary cause of the 2005 Muzaffarabad earthquake (Avouac et al. 2006; Pathier et al. 2006). Active Reasi Fault (RF) and Udhampur Fault (UF) pass through the core center section of the JR in addition to the MBT. As a result of the imbrication of the lower Himalayas into a deeper structural level, the Kishtwar Window (KW) developed inside the crystalline upper Himalayas (Singh 2010; Pandey et al. 2017; Pandita et al. 2021). In 2013, an earthquake of a magnitude of Mw 5.7 struck Kishtwar, with the epicenter 258 km from Jammu. The major tectonic features along with earthquake epicenters having magnitude Mw ≥ 4.1 are shown in the following Fig. 3.

Table 1 Damage and causalities caused due to historical earthquakes around the Jammu Region in the northwestern Himalayas
Full size table
Fig. 2
figure 2

Total number of earthquake events that occurred in the northwestern Himalayas in the last 500 years (1521 AD to 2021 AD)

Full size image
Fig. 3
figure 3

Seismotectonic map of the Jammu Region and surroundings presenting major tectonic features and epicenter of earthquake events with magnitude, Mw ≥ 4.1

Full size image

Post-earthquake effects reported during field investigation

The epicenters of the 2005 Muzaffarabad earthquake and the 2019 Mirpur earthquake are located around 232 and 122 km from Jammu (Fig. 4a). Jayangondapeumal et al. (2008) reported the well-developed liquefaction features at Simbal village near Jammu airport. Open cracks and sand blow on the ground owing to intense motion observed as liquefaction evidence. Through longitudinal fissures, the sand vented water to the ground surface (Rai and Murty 2006; Malik et al. 2007). In May to June 2021, a field survey in the Jatah village, which is next to Degh Nala in the JR’s Samba district. During the 2019 earthquake near Mirpur, a few villagers in Jatah village on the Indo-Gangetic Plain felt strong triggering and noticed the swaying of electric poles. Due to the shaking effect of the far-field earthquake, this event caused a shift in the water table. Local residents observed sand blows at three locations as evidence of liquefaction caused by the far-field earthquake, which had an epicenter in Mirpur, 137 km away. As shown in the following Fig. 4c, the largest sand blows appeared on the east side of Degh Nala. In Simbal, a 20-m long ground rupture was identified, which is created due to strong motions during the 2019 Mirpur earthquake (Fig. 4b). The observed liquefaction characteristics for these two sites which are located at a far distance from the epicenter of both the 2005 Muzaffarabad earthquake and 2019 Mirpur is due to dynamic stress transmission inside the subsurface strata. If the soil has previously liquefied as a result of an earthquake, it may do so again in future earthquake events. Prior seismic strain history has a major impact on the resistance of soils to liquefaction phenomena.

Fig. 4
figure 4

Liquefaction features observed in the Jammu Region (Jammu and Kashmir) during the field survey after the 2019 Mirpur earthquake. a Location map, b ground rupture at Simbal in Jammu, and c sand blows near Degh Nala in Jatah village in Samba

Full size image

Geotechnical, geophysical, and ground response investigation

For a detailed assessment of liquefaction susceptibility of any area in an earthquake-prone zone, subsurface soil data is very important to proceed with the analysis. For this purpose, subsurface geotechnical data is collected from geotechnical consultancies working on different projects related to infrastructure and town planning in Jammu and Kashmir. The standard penetration test (SPT) gives an idea about the N value, density, and plasticity index at different depths for a borehole at a particular site of interest to reflect the penetration resistance (Satyam 2006; Rao and Satyam 2007; Thaker and Rao 2014; Naik and Patra 2018). The soils in the southern sections of Jammu, Samba, and Kathua are generally clayey to clayey silt in the top few meters, followed by sands, and silty sands. The shear wave velocity (\({V}_{s}\)) is associated with the stiffness of soil deposits showing the overall behavior of soil under seismic conditions (Trivedi et al. 2009). For this purpose, the multichannel analysis of surface waves (MASW) method (Park et al. 1999, 2001; Rix et al. 2001; Satyam and Rao 2008; Naik et al. 2021) was used to get a high resolution of subsurface strata at shallow depth. Geophysical approaches have been used in the past to check the susceptibility of sites for liquefaction in tectonically active areas (Park et al. 2002; Yilmaz et al. 2006; Karastathis et al. 2010; Rathod 2011; Rao and Rathod 2014; Satyam and Rao 2014). The geophysical field testing was performed at 243 locations to cover the whole study area using Tromino (Micromed s.p.a. 2009) provided by the Geotechnical Division of National Institute of Technology Srinagar, Jammu and Kashmir, having a frequency range of 0.1–1024 Hz. The testing locations are marked in the map presented in Fig. 1. The recorded data is analyzed to produce Rayleigh wave dispersion curves, which essentially show how the phase velocity of seismic waves in the shallow subsurface changes with frequency. The average shear wave velocity (\({V}_{s30}\)) is calculated for the top 30 m which ranges between 185 m/s on the southern side to 935 m/s on the northern side of the JR.

Depth-wise variation of shear wave velocity (\({V}_{s}\)) and SPT (N) values are presented in Fig. 5. Although it is preferred to establish the shear wave velocity directly from field experiments, doing so at all site locations is generally not cost-effective. Various empirical relationships for determining shear wave velocity using SPT (N) values are known for the Indo-Gangetic area (Rao and Rathod 2014; Naik and Patra 2018). The geotechnical parameters obtained from consultancies and \({V}_{s}\) data obtained from geophysical field testing were integrated for the study region to generate a correlation for all types of soils using regression analysis, as follows:

Fig. 5
figure 5

Depth-wise variation of a SPT (N) value and b shear wave velocity (\({V}_{s}\)) for benchmark locations in the study area

Full size image
$${V}_{s}=211.44\; {N}^{0.2862}$$
(1)

The shear wave and other dynamic soil profiles of the JR can be characterized using the proposed empirical correlation between shear wave velocity and SPT blow counts. The primary soil type in the southwestern part of the study region is sandy silt or silty sand. The performance of the suggested correlation is evaluated by comparing it to the shear wave velocity (\({V}_{s})\) recorded during geophysical testing. The proposed Eq. (1) better fits the majority of the data for shear wave velocity estimated for examined soils. The existence of mechanically unstable soil formations was shown by the delineated clay and peat within the typical foundation depth.

For any site of interest, seismic waves transmitted through a soil layer maybe get amplified or deamplified as compared to bedrock motion. This effect can be studied by considering local site effects to understand the seismic response of soil deposits. Local site effect and extensive site characterization are necessary for places where shallow geology is diverse and strong topographic fluctuations prevail (Bonilla et al. 1997; Nardone and Maresca 2011; Singh et al. 2017; Ferraro et al. 2018; Putti and Satyam 2018; Naik et al. 2020).

In this study, input bedrock motion recorded at nine different seismic stations located in the central, northwestern, and northeastern Himalayas were considered and one-dimensional seismic response analysis for the equivalent linear case has been done using DEEPSOIL (Hashash et al. 2008; Thaker et al. 2012; Jishnu et al. 2013) at all sites where geophysical field testing was carried out within the study area. The bedrock PGA of 0.24 g is suggested by BIS-1893 (2016) for Jammu. In the present study, we used bedrock PGA 0.41 g for safe design purposes considering the maximum credible earthquake scenario (Ansari et al. 2022c). The equivalent linear ground response analysis used in this study was done with the premise that the site is stacked horizontally for all frequency components. Although the comparable linear technique is quick and efficient, it is always an approximate estimate of soil non-linear behavior. The damping versus shear strain curves and modulus reduction curves are used to represent the soil parameters of each soil layer. Seed and Sun (1989) and Seed and Idriss (1970) for clays and sands have employed shear modulus reduction curves (G/Gmax) based on the shear strain to define the characteristics of the soil using discrete points. For the response analysis, Idriss (1990) modulus reduction curves were used.

Nine seismic stations in the Himalayas recorded the motions (IGM1-IGM9) for the 1991 Uttarkashi earthquake, 1995 Chamba earthquake, 1999 Chamoli earthquake, 2010 South Hindukush, and the most recent 2013 Jammu and Kashmir earthquake. All these past earthquake events triggered the Jammu region JR. Among all nine input motions, the 1991 Uttarkashi and 1999 Chamoli earthquakes are very strong. In a way to consider the effect of both the strong and weak nature of motions, the mean value of PGA is preferred. As shown in Fig. 6, the maximum PGA of the mean spectra of these input motions is 0.52 g.

Fig. 6
figure 6

5% damped response spectra of selected input bedrock motions recorded at Himalayan seismic stations for past earthquake events

Full size image

Based on average shear wave velocity at 30 m depth (\({V}_{s30})\), local geology and topography, complete JR is divided into five different classes, class A, class B, class C, class D, and class E. Few sites in the Kishtwar area, in the northwestern part of JR, have hard rock and \({V}_{s30}\) more than 825 m/s; hence, they are classified as class A. Sites in Kishtwar, northern Reasi, Ramban, and Poonch that exhibit having granite rock at shallow depth are classified as class B, with \({V}_{s30}\) ranging from 734 to 820 m/s. All locations in Doda, northwestern side of Jammu, Rajouri, and Poonch show dense gravely sands, and are kept in class C for which \({V}_{s30}\) ranges between 603 to 723 m/s. Most of the sites in southern areas like Jammu, Samba, Kathua, and the central part of Rajouri fall in class D and class E due to sandy deposits. \({V}_{s30}\) fluctuates between 320 and 587 m/s for class D and 185 to 316 m/s for class E, respectively. Due to younger sedimentary formation, areas like Kunjwani, Rehari Colony in Jammu, and Katli and Salor in Kathua have an average shear wave velocity of 208 m/s. Table 2 shows the results of the statistical analysis used to determine the disparity of shear wave velocity for JR.

Table 2 Average shear wave velocity (\({V}_{s30}\))-based statistical parameters for site classification in the Jammu Region
Full size table

Figure 7 shows the comparison of spectral acceleration, Fourier amplitude ratio, and amplification factor for all these five site classes. At 0.3-s time period, spectral acceleration for class A is 1.56 g, which is highest among all sites taken under consideration during geophysical testing. Fourier amplitude ratio (FAR) is determined by taking the ratio of Fourier amplitude spectra of surface motion and Fourier amplitude spectra of input bedrock motion. Class E has the lowest FAR of 2.75 at 8.5 Hz for all locations. The spectral acceleration computed at the surface for any motion at bedrock during the seismic event determines the amplification ratio (AR).

Fig. 7
figure 7

a Spectral acceleration, b Fourier amplitude ratio, and c amplification factor for different site classes in the Jammu Region, based on the shear wave velocity measured during geophysical testing

Full size image

As one move towards site classes with lower shear wave velocity, the peaks of the amplification factor appear to shift towards a higher period. In general, the amplification factor increases with increasing period length, peaks at the lowest period, and then declines with increasing period length. In the southern section of Jammu, the amplification factor for rock sites (class A) is larger than that for alluvium sites (class E). The FAR plot demonstrates that alluvium sites have a maximum peak at lower frequencies, whereas rocky sites have a higher peak at higher frequencies. Over places with a higher shear wave velocity for a shorter period of time, the amplification is greater. For places with lesser shear wave velocity over a long period of time, the amplification is maximum.

For complete JR, the average trend for these parameters is highlighted using a black-dotted line which satisfies the selected boundary condition of input bedrock motions. The graphical values for these black-dotted lines act as a preliminary database to assess the dynamic behaviors of foundation soil in JR. In general, analysis of soil deposits needed for localized construction projects can employ the general trend plot at the initial stages. But site-specific projects need to follow the specific curves provided for various classes as per characterization and microzonation.

The spectral acceleration peaks are maximum for the rocky sites in the northern region as compared to alluvium sites in the southern region covered by Samba, Kathua, and south-eastern parts of Jammu. It has been noted that surface peak ground acceleration (PGAsurface) ranges between 0.13 to 0.55 g. The rocky sites in the eastern part of Kishtwar, Poonch, and Ramban on the northern side present PGA > 0.3 g. Due to young sediments and alluvium deposits near the Tawi River and Ravi River, maximum sites in Jammu, Samba, Kathua, and the central part of Udhampur exhibit very low surface PGA, as shown in Fig. 8. A maximum surface PGA of 0.56 g was reported at Siksha Niketan School at Jeevan Nagar and Govt. Degree College at Bishnah in Jammu. The surface PGA illustrated in Fig. 8 is used as an input parameter to calculate Cyclic Stress Ratio (CSR) while doing a liquefaction hazard assessment.

Fig. 8
figure 8

Spatial distribution of surface peak ground acceleration (PGAsurface) values for Jammu Region using 1D equivalent linear seismic response analysis

Full size image

Methodology

The earthquake-induced liquefaction leads to create catastrophic hazard potential for life killing and property damages. To overcome such natural disasters and associated failure risks, a liquefaction hazard assessment for an area located in a seismically active region is needed. The most commonly used field testing-based approaches for getting undisturbed soil samples for liquefaction hazard study are standard penetration test (SPT), cone penetration test (CPT), and seismic tests to measure the shear wave velocity of soil deposits. The first method evolved based on the concept of applied stress and resistance offered by subsurface soil layers for SPT tests (Seed and Idriss 1971). Robertson and Campanella (1985) were the first to suggest a CPT-related liquefaction model. Idriss and Boulanger (2008) established the liquefaction analysis methodology, and the expression for rd was developed from site response investigations by Boulanger and Idriss (2014). Using a maximum likelihood algorithm and an updated case history database, Boulanger and Idriss (2014) proposed probabilistic CPT-based liquefaction triggering the procedure for cohesionless soils. Several studies have significantly improved and altered the SPT-based approach, as the number of field liquefaction case histories using SPT data has grown (Seed and Peacock 1971; Seed and Idriss 1981; Iwasaki et al. 1982; Idriss 1999; Youd et al. 2001; Boulanger and Idriss 2004, 2012; Cetin et al. 2004; Idriss and Boulanger 2006, 2010). The sites having gravelly soils, cobbles even boulders, and boring becomes difficult for SPT and CPT methods. In such cases, the shear wave velocity (\({V}_{s}\))-based method is adopted for site characterization. The \({V}_{s}\) is a fundamental characteristic of soil material that is directly linked with shear modulus and is used to predict dynamic soil behavior (Kramer 2000). The liquefaction triggering potential may be calculated using existing correlations and in situ \({V}_{s}\) measurements. The shear wave velocity–based approach was used by several researchers to analyze the possibility of soil liquefaction (Lin et al. 2004; Beroya et al. 2009; Karastathis et al. 2010; Castellaro et al. 2015). Qualitative information on subsurface conditions was acquired through shear wave velocity–based frequency measurements to generate a hybrid hazard map for zonation of liquefaction sites (Huang and Tseng 2002; Beroya et al. 2009; Satyam and Priyadarsini 2021).

Based on \({V}_{s}\), a relationship to assess the cyclic resistance of subsurface strata was created by considering liquefaction case histories from 20 earthquakes and 193 liquefaction and non-liquefaction case histories throughout the world (Andrus and Stokoe 1997). For Andrus and Stokoe (2000) method, cyclic loading due to earthquake, overburden stress corrected shear wave velocity, and soil resistance against liquefaction are the main parts to calculate the \({F}_{s}\) for a site susceptible to liquefaction during strong seismic waves. Andrus et al. (2004) adjusted this relationship by introducing a correction factor for cemented and old soil deposits. In this study, the most recent technique suggested for SPT (Boulanger and Idriss 2014) and shear wave velocity (Andrus and Stokoe 2000; Andrus et al. 2004) is used to assess the liquefaction hazard in the study area.

The factor of safety (\({F}_{s}\)) is calculated by taking the ratio of cyclic stress ratio (CSR) and cyclic resistance ratio (CRR). The CSR depends on the surface PGA and can be calculated as follows:

$$\mathrm{CSR}=0.65 . \frac{{\sigma }_{v}}{{\sigma }_{v}^{^{\prime}}} . \frac{{a}_{\mathrm{max}}}{g} . {r}_{d}$$
(2)

In the above Eq. (2), \({\sigma }_{v}\) is the vertical total stress, \({a}_{\mathrm{max}}\) is the maximum PGAsurface, and \({r}_{d}\) is a shear stress reduction factor, defining the dynamic response of soil deposits. The cyclic resistance ratio (CRR) is associated with the SPT (N) value, which is corrected for overburden stress, delivered energy, rod diameter, energy ratio, and type of sampler (Boulanger and Idriss 2014). Duration of ground shaking is also one of the major factors to be considered in terms of magnitude scaling. The empirical equation to calculate CRR is given as follows:

$${\mathrm{CRR}}_{M,{\sigma }_{v}^{^{\prime}} }= {\mathrm{CRR}}_{M=7.5,{\sigma }_{v}^{^{\prime}}=1\; \mathrm{atm} }. MSF . {K}_{\sigma }.{K}_{\alpha }$$
(3)

Here, in this Eq. (3), \(\mathrm{MSF}\) is a magnitude scaling factor, \({K}_{\sigma }\) is overburden correction factor, and, \({K}_{\alpha }\) is the static shear stress correction factor. The correlation for \({\mathrm{CRR}}_{M=7.5,{\sigma }_{v}^{^{\prime}}=1 \mathrm{atm}}\) is as follows,

$${\mathrm{CRR}}_{M=7.5,{\sigma }_{v}^{^{\prime}}=1 \mathrm{atm} }=\mathrm{exp} \left\{\frac{{\left({N}_{1}\right)}_{60cs}}{14.1}+ {\left[\frac{{\left({N}_{1}\right)}_{60cs}}{126}\right]}^{2}- {\left[\frac{{\left({N}_{1}\right)}_{60cs}}{23.6}\right]}^{3}+ {\left[\frac{{\left({N}_{1}\right)}_{60cs}}{25.4}\right]}^{4}-2.8\right\}$$
(4)

where \({\left({N}_{1}\right)}_{60cs}\) is the equivalent clean sand value for cohesionless soils. The \(\mathrm{MSF}\) can be estimated using the following Eq. (5).

$$\mathrm{MSF}=6.9 . \mathrm{exp} \left\{\frac{-M}{4}-0.058\right\} \le 1.8$$
(5)

According to Andrus and Stokoe (2000) method, for shear wave velocity–based evaluation of liquefaction hazard, CRR can be calculated using Eq. (6), as presented below.

$$\mathrm{CRR}=\left\{a {\left[\frac{{K}_{c}{V}_{s1}}{100}\right]}^{2}+ b \left(\frac{1}{{V}_{s1}^{*}-{K}_{c}{V}_{s1} }- \frac{1}{{V}_{s1}^{*} }\right)\right\} . \mathrm{MSF}$$
(6)

In the above Eq. (6), \({K}_{c}\) is the correction factor considered for the high values of overburden stress–corrected shear wave velocity (\({V}_{s1}\)) due to cementation and aging, ranging between 0.6 to 0.8 (Ohta and Goto 1978; Robertson et al. 1992; Rollins et al. 1998). The value of \({V}_{s1}^{*}\) varies between 200 and 215 m/s, based on a percentage of fine content (FC). The \({V}_{s1}^{*}\) and \(FC\) are related as follows:

$$\begin{aligned}&V_{s1}^\ast=215\;\mathrm m/\mathrm s,\;\mathrm{for}\;\mathrm{sands}\;\mathrm{with}\;FC\leq5\%\\&V_{s1}^\ast=215-0.5(FC-5)\;\mathrm m/\mathrm s,\;\mathrm{for}\;\mathrm{sands}\;\mathrm{with}\;5\% < FC <{35}\%\\&V_{s1}^\ast=200\;\mathrm m/\mathrm s,\;\mathrm{for}\;\mathrm{sands}\;\mathrm{with}\;FC{\geq}35\%\end{aligned}$$
(7)

For older soil deposits, \({V}_{s1}\) is corrected for an equivalent young, clean soil as follows:

$${{(V}_{s1})}_{csa1}= {K}_{a1} . {K}_{cs} . {V}_{s1}$$
(8)

where \({K}_{a1}\) is an age factor and is considered equal to 1.0 for all Holocene soil deposits (Andrus and Stokoe 2000). The \({K}_{cs}\) is \(FC\) correction factor which can be determined as follows (Juang and Jiang 2000; Juang et al. 2002),

$$\begin{aligned}&K_{cs}=1.0,\;\mathrm{for}\;FC\;\leq5\%\\&K_{cs}=1+(FC-5)T,\;\mathrm{for}\;5\% < FC < 35\%\\&K_{cs}=1+30T,\;\mathrm{for}\;FC\geq35\%\end{aligned}$$
(9)

where,

$$T=0.009-0.0109\; ({V}_{s1} /100)+0.0038\;{({V}_{s1} / 100)}^{2}$$
(10)

  

Liquefaction potential index (LPI)

The \({F}_{s}\) is insufficient for evaluating liquefaction hazards (Toprak and Holzer 2003). Hence, the depth and thickness of layers susceptible to liquefaction are considered to understand the overall damage of any borehole. The liquefaction potential index (LPI) considers both the geometrical parameter of the liquefiable layer and the associated \({F}_{s}\) (Iwasaki 1978; Iwasaki et al. 1982). The \(\mathrm{LPI}\) map for the Kashmir valley was proposed based on SPT-based method, which suggests that maximum boreholes in Baramulla and Kupwara demonstrate very high LPI (Dar and Dubey 2015; Khan and Shah 2016; Sana and Nath 2016; Zahoor et al. 2019). The site condition for different values of \(\mathrm{LPI}\) and, the required plan of action for liquefaction mitigation is mentioned in Table 3.

Table 3 Likelihood of liquefaction at different probability (Iwasaki 1978; Iwasaki et al. 1982)
Full size table

The LPI can be determined as follows:

$$\mathrm{LPI }=\underset{0}{\overset{20}{\int }} F\left(z\right) . W\left(z\right) dz$$
(11)
$$\begin{aligned}&F\left(z\right)=1-F_s,\;\mathrm{for}\;F_s<1.0\\&F\left(z\right)=0,\;\mathrm{for}\;F_s\geq1.0\\&W\left(z\right)=10-0.5z,\;\mathrm{for}\;\mathrm z<20\;\mathrm m.\\&W(z)=0,\;\mathrm{for}\;\mathrm z>20\;\mathrm m\end{aligned}$$
(12)

Probability of liquefaction (\({P}_{L}\))

Based on the concept of LPI, the probability of liquefaction (\({P}_{L}\)) for any site can be estimated using Eq. (13), proposed by Juang et al. 2002. Table 4 presents the different states for the site under investigation based on the calculated probability of liquefaction.

Table 4 Liquefaction risk assessment based on probability of liquefaction (Chen and Juang 2000; Juang et al. 2002)
Full size table
$${P}_{L}= \frac{1}{1+ {\left(\frac{{F}_{s}}{0.96}\right)}^{4.5}}$$
(13)

Earthquake-induced liquefaction hazard investigation for the study area

Based on SPT (Boulanger and Idriss 2014), data collected from geotechnical consultancies and shear wave velocity (Andrus and Stokoe 2000; Andrus et al. 2004) estimated from multichannel analysis of surface waves (MASW) during field testing, the comprehensive liquefaction hazard assessment of the JR was carried out. Based on the calculated values of \({F}_{s}\), three distinct categories are fixed to define the potential liquefaction hazard for a particular site. The sites with \({F}_{s}\) < 1.0 are vulnerable to severe liquefaction, whereas sites with \({F}_{s}\) between 1.0 and 1.2 are likely to liquefy. Sites with \({F}_{s}\) > 1.2 are safe and not vulnerable to liquefaction. Table 5 demonstrates typical analytical findings for determining the \({F}_{s}\) against liquefaction at benchmark sites using both methodologies. The areas on the northeastern side like Rinaie, Shareef Bagh, Sonder in Kishtwar, and Gattigali and Punaja in Doda are not susceptible to liquefaction, as shown in Fig. 9a. In Gandhi Nagar, Talab Tillo, and Hakal in Jammu, along with several localities near the Tawi River’s bank including Simbal, Ram Bagh, and Jammu University, and eastern sections of Rajouri such as Laal Haveli and Jamia Masjid, the liquefaction hazard is severe. Samba and Kathua have a thick layer of young sedimentary formations with uniformly graded soils showing highly susceptible to liquefaction. Central portions of Udhampur, such as Jinghanu and Panchari, have a sand layer up to 4 m thick, indicating that they are expected to liquefy. Baflaiz, Morha near Mughal Road in Poonch, and Kauri and Bakkal village near Chenab bridge in Reasi exhibit exceptionally high \({F}_{s}\) against liquefaction, due to gravelly sand with high SPT (N) values. The increase in liquefaction resistance over time might be due to the deformation or compression of soil particles into more stable arrangements, as well as a higher degree of consolidation. The maps presented in Fig. 9b, c show the liquefaction hazard zonation based on the \({F}_{s}\) values determined using the methods suggested by Andrus and Stokoe (2000) and Andrus et al. 2004, respectively. Due to rock near the surface or at shallow depth, liquefaction threat is low in most portions of the eastern and central regions, according to these two maps: Sangaldan, Chachihal, and Banihal in Ramban, which are located in the northern part of the study area, have a \({F}_{s}\) against liquefaction more than 1.2.

Table 5 Liquefaction analysis for benchmark locations in the study area
Full size table
Fig. 9
figure 9

Liquefaction hazard assessment of Jammu Region using a SPT-based method (Boulanger and Idriss 2014), b shear wave velocity (\({V}_{s}\))-based method proposed by Andrus and Stokoe 2000, and c shear wave velocity (\({V}_{s}\))-based method proposed by Andrus et al. 2004

Full size image

Using the findings of both \({V}_{s}\)- and SPT-based techniques, an integrated liquefaction hazard map of the study area is also presented and illustrated in Fig. 10. As per this map, Banhore, Sawan Chak, Balour, and Talhar like areas near the bank of the Ravi River in the southern region have extremely low values of SPT (N) and \({V}_{s}\), making them especially vulnerable to liquefaction. Liquefaction is not occurring in the eastern and western parts of the JR which has high shear wave velocities and also not occurring in the central portion because of dense sand. The areas like Babliana, Mehmoodpur, Batna, Khan Pora, and Sawalkot is likely to liquefy as \({F}_{s}\) values vary between 1.0 and 1.2.

Fig. 10
figure 10

SPT and shear wave velocity–based integrated liquefaction hazard map of Jammu Region

Full size image

The \(\mathrm{LPI}\) of any site offers a better understanding of the borehole related to liquefaction. The determined \({F}_{s}\) values at different depths and thickness of layering for each borehole used to calculate the LPI for the sites under study. According to Iwasaki et al.’s (1982) method, the whole study area is classified into four different zones based on LPI values, as shown in Fig. 11. Based on LPI, boreholes near the bank of Tawi River and Ravi River in southern districts like Jammu, Samba, and Kathua have a very high liquefaction potential due to high water tables. The Simbal village in Jammu and SIDCO industrial area in Kathua show LPI > 15. The borehole at Jatah in the Samba district where sand blows were observed during the field survey fall in a high liquefaction potential zone. The maximum of the boreholes in Ramban, Rajouri, and Poonch have a 3–5-m gravel layer exhibit very low LPI, and are not vulnerable to liquefaction. The LPI for boreholes in Reasi and Udhampur ranges from 4.15 to 12.34, indicating a low to high liquefaction potential. The probability of liquefaction (\({P}_{L}\)) suggests the likelihood of occurrence of liquefaction features for any site in a seismically active region based on the \({F}_{s}\) against liquefaction (Juang et al. 2002). As shown in Fig. 12, the maximum portion in the northern part of the JR demonstrate \({P}_{L}\) < 0.42. The western portion of Jammu, Udhampur, and Kathua are likely to liquefy where \({P}_{L}\) varies between 0.65 and 0.85. As \({P}_{L}\) is less than 0.1 in Poonch and Doda, all sites are either unlikely to liquefy or will not liquefy. The \({P}_{L}\) is very high in regions like Chak Lalushah and Babliana in Jammu, as well as areas near the Ravi River in Kathua due to uniformly graded soil at shallow depth.

Fig. 11
figure 11

Liquefaction index (LPI) map of the Jammu region to identify the liquefaction vulnerability of sites subjected to strong ground motion during seismic events, with black stars indicating sites where liquefaction features were documented

Full size image
Fig. 12
figure 12

Spatial distribution of the probability of liquefaction for the Jammu region to check the liquefaction hazard associated with an earthquake

Full size image

There is no published study on liquefaction hazard assessment for the JR. Sana and Nath (2016) evaluated earthquake-induced liquefaction hazards across the whole Kashmir valley if we were looking for a study in and around Jammu. Except for the hills and places to the north of Dal Lake (Illahibagh and Zakura), which have minimal liquefaction susceptibility due to the presence of stiff soils, the majority of the lands have a moderate to extremely high liquefaction potential in Srinagar city (Zahoor et al., 2019). The probability of liquefaction estimated for different cities in the Himalayan region are calculated by different researchers are compared in Table 6.

Table 6 Probability of liquefaction (\({P}_{L}\)) for the JR compared to previous studies done for various cities in the Himalayas
Full size table

Conclusions

The Jammu Region (JR) in the northwestern Himalayas is jolted by frequent low to high-magnitude earthquakes including the deadliest 2005 Muzaffarabad earthquake and most recent 2019 Mirpur earthquake. Due to ongoing infrastructure development and urban expansion in such a seismically active area drives to perform liquefaction hazard assessment. In this study, the factor of safety (\({F}_{s}\)) against liquefaction is calculated using both SPT- and shear wave velocity (\({V}_{s}\))-based field methods and results are superimposed to generate an integrated map of liquefaction hazards in the JR. The subsurface borehole data was provided by local geotechnical consultancies, and the shear wave velocity was recorded at 243 places during multichannel analysis of surface waves (MASW)-based, geophysical testing using Tromino (Micromed s.p.a. 2009). During the field survey, a 20-m long ground rupture at Simbal near Jammu airport and sand blows near Degh Nala in Jatah were observed as liquefaction evidence of far-field 2019 Mirpur earthquake in Pakistan. The analysis of results showed that some areas are prone to liquefaction using the SPT approach but are not susceptible to liquefaction using the \({V}_{s}\) approach. The integrated liquefaction hazard map suggests that liquefaction is severe in the southern portion of the JR where CRR is less than CSR. Maximum parts of Udhampur and Rajouri are anticipated to liquefaction due to saturated loose to medium dense cohesionless soils. The northern part of Jammu, Reasi, and Ramban together with Kishtwar and Poonch have rocky strata at shallow depth and very high \({V}_{s}\), not susceptible to liquefaction.

The findings from both approaches to determine \({F}_{s}\) at different depths and thicknesses of each layer in a borehole used to prepare the map for liquefaction potential index (LPI) for the whole study area. This map illustrates that the central part has low to high, the northwestern part low to very low, and the southern part has a very high \(\mathrm{LPI}\). The southern portion of Jammu, Samba, and Kathua located on the Indo-Gangetic Plain have young alluvium and soft sedimentary deposits that reveal a very high probability of liquefaction. The sites in Doda, Poonch, and Kishtwar have very high \({V}_{s}\) and rocky strata exhibit a very low probability of liquefaction. The liquefaction hazard maps of the present study can be used for future design and analysis of infrastructures subjected to earthquake risk.