Site specific seismic hazard analysis of monumental site Dharahara, Kathmandu, Nepal

Abstract This study was carried out to estimate the seismic requirement for analyzing and designing the monumental structure, Dharahara, in Kathmandu, Nepal. A probabilistic seismic hazard analysis of the study region was performed to develop the elastic response spectrum for the bedrock site. The obtained spectrum was compared with the existing seismic codes, NBC 105:2020 and IS 1893:2016, to estimate the target response spectrum. One-dimensional equivalent linear (EQL) and non-linear (NL) ground response analyses were performed using DEEPSOIL v7, incorporating a suite of ground motions. The three borehole logs exhibited very soft soil in the region. Due to the presence of soft soil strata in the region under study, strength-corrected modulus reduction and damping curves were adopted in the current study to overcome the effect of larger shear strain. The average response spectra from EQL analyses predicted higher values than NL analyses. Finally, a comparison of the proposed design response spectrum with the code spectra was performed. Comparisons revealed that the mean EQL and NL outputs fall in between the design spectra of similar soil classes in NBC 105:2020 and IS 1893:2016 seismic codes.


Background
Dharahara was built by Bhimsen Thapa, the first Prime Minister of Nepal, in 1824 B.S and was the tallest tower (72.96 m) at the time of construction and is located in Sundhara, Kathmandu, Nepal. Dharahara is said to have been constructed for Queen Lalit Tripurasundari, the niece of Bhimsen Thapa. It was partially damaged during the 1833 major earthquake in Nepal and was restored afterward (Hutt 2019). It was significantly damaged by the Nepal-Bihar earthquake (M w ¼ 8.1) in 1934. After the 1934 earthquake, this historical lime masonry monument was raised to 61.88 m in height by Juddha Shamsher Jung Bahadur Rana, the Prime minister of Nepal at that time. In the 2015 Gorkha earthquake, this tower crumbled and collapsed, resulting in more than 180 casualties (Wang et al. 2016). Dharahara comprises the central core column with a diameter of 900 mm and the outer wall with varying thicknesses of 1.828 m at the bottom and 0.736 m at the top. The structure was a nine-story with a total height of 61.88 m from the base and 54.102 m from the top of the exterior stone footing steps. Also, the slab of cantilever verandah was projected at the top with a width of 900 mm and a thickness of 75 mm. Recently, after the Gorkha earthquake, in 2015, the Dharahara was reconstructed and will look similar to the old one; however, it would be equipped with modern facilities. As per the reconstruction plan, the new Dharahara would be 245 ft (72.96 m) tall with 11 stories when looking at the outside. However, it would comprise 21 stories from the inside. Since the damaged Dharahara was built with brick masonry long ago with traditional technology, the tower has historical value. Moreover, several cultural heritage structures in Nepal are in a need of seismic vulnerability and risk assessment as the country lies in a highly seismic region. Therefore, this study performs a site-specific seismic hazard assessment of the site comprising monumental structures located in the region of high seismicity and soft soil conditions.
The recurrence rate of a greater earthquake is high in Nepal as it is narrow and elongated along Greater Himalayas and lies towards the southern collisional boundary of the Indian and the Eurasian Plate (Ambraseys and Douglas 2004). Within the narrow width of Nepal three fault systems-Main Central Thrust (MCT), Main Boundary Thrust (MBT), and Himalayan Frontal Thrust (HFT) pass from the east to the west throughout the length with many other smaller 92 faults (BECA 1993). The Geological map including major faults with MFT, MBT, MCT, and several small faults is illustrated in Figure 1 based on the study of existing literature Pandey et al. (2002), Parajuli et al. (2021), Stevens et al. (2018), Rahman and Bai (2018), and Bothara et al. (2002). It is worth mentioning that Nepal accounts for about 17% of the largest earthquakes in the world (Ulak 2015). Along with the faults extended over the country, the geography is too complex. The higher Himalayas with an altitude of 5000 m and above consist of mainly vertical and steep rocky slopes with an excessive amount of quartzite, gneiss, marble, and schist. Similar to the higher Himalayan region the lesser Himalayan zone also consists of rock deposits and this zone is very active as well as uplifting at a high rate. In the midlands or Mahabharat range, soil formation consists of deeply weathered rocks and limestone, dolomite marble, and granites. Siwaliks which lie in the lowermost part of Nepal (Terai) comprise river deposits and coarse sediment in the north which is near the base of the Siwaliks range (Amatya and Jnawali 1994;Dahal 2010). Terai comprises several Dune valleys like in eastern Trijuga, western Dang, and central Chitwan and is filled with coarse to fine alluvium soil. Moreover, the Jhiku Khola fault, Chitlang fault, Kulekhani fault, Malagiri fault, and Kolphu Khola fault situated around the Kathmandu valley can generate potential earthquakes (Ulak 2015). Basic Exchange and Cooperation Agreement for Geo-spatial Cooperation (BECA) (BECA 1993) conducted a probabilistic seismic hazard analysis (PSHA) and risk assessment of Nepal. This was the early study performed for developing Nepal Building Code (NBC-105105, 1994; hereafter referred to as the NBC building code), however, many researchers (Pandey et al. 2002;Parajuli et al. 2021;Stevens et al. 2018;Rahman and Bai 2018;Ram and Wang 2013) have carried out PSHA of Nepal afterward. In the study, the researchers have adopted different earthquake sources and ground motion predictive equations (GMPEs). As a result, different levels of peak ground accelerations (PGA) were observed at the same locations. The PGA variation of the Kathmandu valley was estimated to be in the range of 0.4-0.55 g. However, in the revised NBC building code of Nepal (NBC 105 105: 2020), a PGA of 0.35 g (for bedrock ground motion) is recommended for Kathmandu valley.
The goal of earthquake resistance design is to design a fully functional structure with minimal damage during earthquake shaking which is achieved with the proper estimation of ground motion. The two hazard analysis approaches widely adopted includes deterministic and probabilistic seismic hazard analysis, that is, DSHA and PSHA, respectively. DSHA was the early developed hazard analysis and involves expert decisions regarding the identification of the earthquake source and its potential (Reiter 1991). The difference in these two approaches lies in the consideration of sources, that is, the DSHA accounts for scenario earthquake (considers worst-case ground motion) while PSHA accounts for the overall seismicity (likelihood occurrence of controlling earthquake, level of shaking, etc.). Hence, PSHA was adopted to evaluate the overall hazard of the current studied region. The study of seismic hazard analysis begins with the identification and delineation of earthquake sources (Kramer 1996).
Even though PSHA depicts the overall seismicity of the region with consideration of spatial uncertainty, size uncertainty, and temporal uncertainty for the estimation of the ground motion, the influence of local site conditions might also result in the variation of ground motion. The non-uniformity in soil type and the complex geology cause different levels of shaking and destruction during the earthquake. Hence, the major objective of the current study is to assess the seismicity of the study area and to perform the 1D site response analysis for site-specific seismic hazard assessment of Dharahara, Sundhara Kathmandu (27 42 0 2.52 00 , 85 18 0 42.84 00 ).

Geology of the Kathmandu valley and relevance of site-specific hazard analysis
In the case of Kathmandu valley, the soil beneath the city is mainly formed by young soft fluvio-lacustrine sediments. The southern part of the valley is dominated by siltyclayey soil deposits whereas the northern part (around the Bagmati river) is dominated by sand and fluvial gravel deposit. Kathmandu consists of consolidated fluvio-lacustrine Quaternary sediments located above the bedrock with lacustrine and fluvial soil deposits up to 500 m in depth (Yoshida 1984;Dill et al. 2003;Gaha et al. 2022;Bhusal and Paudel 2021). The valley once was a lake of soft soil (lake) deposit up to 100 m thickness which leads to the possibility of amplification of ground motion resulting in excessive destruction of structures (Kumahara et al. 2016). Hence, for Kathmandu valley, NBC 105:2020 has specified to consider very soft soil named soil type D (hereafter referred to as NBC soft soil type D). This depicts the possibility of amplification of seismic waves as it passes through the soft soil layers, especially for long periods. The ground motion parameters like the amplitude of ground motion, frequency content, and duration are altered due to the local soil site conditions (Kramer 1996). The damage might even be excessive under minor bedrock shaking due to amplification of the ground motion. Therefore, due to the soil profiles in several regions of Nepal and the seismicity of the region, the seismic hazard analysis of the studied site location was deemed necessary. This is also because, during the development of the NBC building code, the average seismicity of the region is considered for seismic zonation of the whole country. However, the seismicity and damage during the seismic event might vary based on the local site conditions. These scenarios have already been experienced during past earthquakes, for example, the Michoacan earthquake of moment magnitude, M w ¼ 8 (Mexico, 1985) resulted in moderate damage in the nearby region of its epicenter but resulted in excessive damage some 350 km away in Mexico City (Celebi et al. 1987). Also, during the Loma Prieta earthquake, the Modified Mercalli Intensity (MMI) VII scale was experienced in the epicentral region, but the intensity of MMI IX was felt some 100 km away from the epicenter (Seed 1990). Moreover, amplification of seismic waves and severe destruction were observed at locations around 350 km away from the epicenter during Bhuj EQ (M w ¼ 7:7) in 2001 (Kumar et al. 2017). Further, it was also observed that during seismic events like Hyogo-Ken Nanbu (Kobe) earthquake in Japan (1995), Spitak earthquake in Armenia (1988), Chi-Chi earthquake in Taiwan (1999), etc., exceeded the provisions stipulated in the code (response spectrum; Stewart 2008; Griffiths et al. 2016) studied the challenges related to site response analysis for soft soil under the ground motion of high intensity. The use of equivalent linear (EQL) and non-linear (NL) analyses for soft soil sites resulted in the estimation of large shear strains (3-10%) which finally resulted in a typical spectral shape. It was recommended to modify the modulus reduction and damping curves for strains greater than approximately 0.1% to obtain a realistic soil model for site response analysis. Also, to have a realistic representation of the shear strength of soil, Stewart and Kwok (2008), Kwok et al. (2008), and Yee et al. (2013) have formulated several methods to modify the G/G max curve at a higher value of shear strain. These modifications can be adopted to produce G/G max at shear strain values greater than 1%. The work of Stewart and Kwok (2008) and Kwok et al. (2008) was further modified by Hashash et al. (2010) for a non-linear site response analysis. As mentioned in the study by Griffiths et al. (2016), Aaqib et al. (2018), and , the modification by Hashash et al. (2010) usually produces a reliable estimate of ground shaking and shear strain.
Several numerical tools are developed to perform the 1D ground response analysis. These tools adopt analysis methodologies like linear, equivalent linear (EQL), and non-linear (NL) analysis for ground response analysis. The selection of the analysis method depends on the assumption adopted during numerical modeling, availability of input parameters, objective, and significance of the research problem (Kramer 1996). In the current study, equivalent linear (EQL) and non-linear (NL) approaches were adopted for ground response analysis. The EQL analysis is known to estimate the soil response under seismic simulation with reasonably good accuracy at small strains. In this analysis method, the dynamic property of soil is adjusted based on the compatibility equation of state (strain compatibility). Even though the EQL approach is computationally efficient, it still provides an approximation of the actual inherent non-linearity during the seismic ground response. Kramer (1996) mentioned that the EQL is not capable to capture the stiffness degradation of the soil during cyclic loading. Also, the behavior of soil during seismic loading is highly non-linear. In the case of the NL approach, the actual non-linear response of the soil deposit can be estimated based on the direct numerical integration in the time domain. The equation of motion can be integrated into small steps with the input of the advanced constitutive model depicting the non-linear response of the soil. Kim et al. (2016) have highlighted the relative difference between EQL and NL analysis during 1-D site response analysis. The authors observed that at high frequency, NL results are higher than EQL due to the overdamping. However, near the resonant site frequency, NL ordinates are lower than EQL due to the formation of resonant response. The abovementioned difference was observed to be less distinct for spectral acceleration ratios as compared to Fourier amplitude ratios. The analysis methodology and details of ground response analysis adopted in current study are explained in Section 3.

Need for the assessment of seismic hazard for cultural heritage site
The major purpose of Dharahara during its construction period was to summon the people and provide information on government announcements. Its construction style also depicts the harmony between all religious beliefs as the minaret was constructed in Mughal style with a European-style railing and the top floor comprising a Shiva shrine (Hutt 2019). It is a symbol of religious unity, a cultural icon of Kathmandu, and is documented as a national heritage site. Dharahara has historical and cultural significance, therefore, during the design of such a structure, the seismicity and ground response analysis of the site should be performed.
The impact of the seismic event on world heritage cultural sites is beyond the economic loss. The consequences of such seismic events on cultural sites might result in loss of employment, and disruption of tourism to name a few among others. To mitigate the earthquake risk for European regions and towns RISK-UE project funded by the European Commission was started in 1999. The project proposed a macro seismic model and mechanical model for the assessment of vulnerability and the evaluation of the seismic risk scenarios. The macroseismic model uses macro seismic hazard maps and vulnerability curves while the mechanical model uses capacity curves to evaluate the vulnerability of the existing building stock. Lagomarsino and Giovinazzi (2006) and Despotaki et al. (2018) performed the simplified assessment for the seismic risk evaluation of the heritage site in Europe. The hazard curves were developed from the hazard model SHARE and were later combined with the fragility function for the calculation of the probability of damage or collapse of the structure. Melissianos et al. (2021) conducted the risk assessment of monolithic columns of the monumental structure of Greece. PSHA followed by site-specific hazard analysis was carried out for the site comprising Aphaia temple in Aegina. Finally, the fragility analysis was performed to estimate the seismic risk of the column for seismic excitation. The study suggested various limit states which define the damage and collapse of the column. Hassan et al. (2020) performed the modeling of site-specific ground motion for a historical site and further calculated the required seismic demand for the studied site. The local site effect was considered to perform the site-specific seismic hazard analysis and to compute the response spectra. The research work suggests a suitable seismic assessment methodology for the preservation and risk reduction of the heritage site under earthquake shaking.
Such studies are helpful to study the nearby monumental structures and regions where similar hazards and identical ground motions are expected. The performance of such monumental structures during extreme seismic events can also be used as a knowledge transfer of the constructional aspect of such structures to the next generation. Hence, reducing the seismic risk associated with the damage to cultural heritage has encouraged conducting several research works on collaboration in several countries. One of its kind is the PROTHEGO (PROTection of European Cultural Heritage from Geo-hazards) (Margottini et al. 2017). The project usually adopts InSAR (radar interferometry) techniques to monitor the monumental structures located in Europe which were damaged during seismic events. This research work helps to enhance the national level management of the cultural heritage practices of the country, identify the potential damages and risk associated with the failure of such structures, enhance the institutional support through the transfer of knowledge and innovation, and improve the performance of such structures by adopting suitable strengthening approaches before (disaster preparedness) or after major seismic events. Moreover, it should be mentioned herein that the assessment of vulnerability and associated risk of cultural heritage structures is not only defined by the geological characteristic alone but also by the structural characteristics and the type of monument. The monuments are usually clustered into different typological categories like churches, temples, towers, palaces, etc., and are associated with the belief of local people residing in the area.
In the current study, the site response analysis for the city, Kathmandu comprising cultural heritage structures is performed. Kathmandu is the historically enriched capital city comprising several structures of historical importance. The past seismic events have always resulted in severe damage to the structures and are mainly due to local amplification of soil as it has a significant impact on increasing the scale of destruction. Hence, the estimation of seismic demand of Kathmandu valley is deemed required to carry out the planning for the preservation of historic structures. Moreover, the site response analysis of the study region is also required as the region comprises very soft soil which might result in the amplification of seismic waves. In the current study, PSHA followed by 1D site response analysis was carried out and elastic design spectra were developed for location comprising monumental structure (location comprising Dharahara). This study is unique as there are no such studies carried out for such regions comprising soft soil.

Probabilistic seismic hazard analysis
Seismic ground motion is uncertain, unpredictable, and occurs with spatial and temporal uncertainty without early warning. PSHA involves the consideration of spatial uncertainty, size uncertainty, and temporal uncertainty for the estimation of ground motion. The identification of sources is very complex and is one of the most important steps for the estimation of the likelihood of occurrence of ground shaking for a particular site. Several researchers (Parajuli et al. 2010;Pandey et al. 2002;Ram and Wang 2013) have performed the PSHA of Nepal by using a different source. The researchers have estimated the different hazard levels. In the current study, mixed types of sources adopted in the study by Parajuli et al. (2021) were considered. A total of 1228 (M w ! 4.5) historic seismic events were collected within the rectangular coordinates bounded by coordinates 78 30 0 00 00 E, 25 30 0 00 00 N) and 89 30 0 00 00 E, 31 30 0 00 00 N around the study area. To address the size uncertainty, the minimum and maximum magnitude of the earthquake were adopted. The minimum magnitude of Mw 4.5 and maximum magnitude based on the equation proposed by Wells and Coppersmith (1994) were adopted for the sources adopted in the current study. The temporal and spatial biasness of the seismic events needs to be filtered out so that the undesired aftershocks can be removed from the recorded seismic database. The removal of aftershocks was performed using Gardner and Knopoff (1974). Then, a completeness analysis of each seismic source was carried out following the study Stepp (1972). Using the completeness period, Gutenberg Richter relationship (recurrence law) was developed for the best fit of the frequency magnitude relationship. The completeness of earthquake and estimation of recurrence parameters for the central part of the study area are plotted in Figure 2. Moreover, to estimate the seismic hazard, the proper selection of GMPEs also has the major role. Hence, the subsequent step of seismic hazard analysis is to identify the proper attenuation relationship for the study region. To the best of the author's knowledge, not a single attenuation relationship has been developed for the studied region to date. In such a scenario, the suitable attenuation relationships developed for a similar tectonic environment from a global database should be considered. The ground motion predictive equation adopted in the current study was adopted from Parajuli et al. (2021). Therefore, since the study comprises several sources and predictive equations, a logic tree approach (see Parajuli et al. 2021) with equal weights was adopted. Several attenuation laws developed for the subduction interface and active shallow crustal earthquakes were considered by the logic tree approach (Parajuli et al. 2021) while performing PSHA. In this study, the bedrock shear wave velocity is assumed to be 760 m/s, as per the recommendations of Asimaki et al. (2017) for the Kathmandu valley. The return period of 475 years (10% probability of exceedance in 50 years) is typically used as a design basis earthquake and was provided as an input in the Open Quake Engine developed by global earthquake model (GEM) , .
The distribution of PGA for the different locations of Nepal and the studied region are shown in Figure 3(a) and (b), respectively for a 10% probability of exceedance in 50 years. It can be observed that the variation of PGA is very less due to the narrow studied area (Sundhara, Kathmandu; see Figure 3(b)). Also, the rock Uniform hazard spectra (UHS) for 475 years of return period for the Sundhara region having PGA 0.424 g was estimated and illustrated in Figure 3(c). The maximum spectral acceleration was calculated to be 1.07 g at 0.2 s. Baker (2011) has provided the methodology to estimate the target spectrum based on conditional probability. The estimation of conditional probability depends on the fundamental period of vibration and GMPEs. The estimation of the exact fundamental period is a bit complex for the current studied monumental structure. Moreover, the single GMPE is not developed for the studied region so multiple GMPEs were adopted. The consideration of all the GMPEs to estimate the conditional mean spectrum (target spectrum) is more complex.
Hence, in the current study, the uniform hazard spectrum was adopted to estimate the elastic design spectrum. Newmark and Hall (1982) developed the procedure to estimate the elastic design spectra based on PGA, soil type, damping ratio, and percentile spectrum amplification factor. Spectral amplification factors for rock and competent soil conditions as the ratio of v/a ¼ 36 (in/s)/g and ad/v 2 ¼ 6 were developed where v, a, and d represent velocity, acceleration, and displacement, respectively. Also, the spectral bounds equal to 2.71, 2.3, and 2.01 were used for acceleration, velocity, and displacement components for unit PGA (Newmark and Hall 1982). These bounds were developed for an 84.1 percentile spectrum amplification factor, and a 5% damping ratio, and were used for the construction of design spectra. A typical elastic design spectra is shown in Figure 4(a) which was developed by Newmark and Hall (1982) to construct smooth elastic design spectra. Figure 4 illustrates that distinct periods were fixed, and the spectrum was divided into acceleration, velocity, and displacement-sensitive regions. Elastic design spectra were constructed for bedrock using spectral bounds. A similar methodology was used by Baruwal et al. (2020) for developing the design spectra for Pokhara valley, Nepal, and therefore, is adopted in the current study.

Design spectra
For the bedrock designation, IS 1893:2016 have adopted soil type I, and NBC 105:2020 have adopted soil type A. A comparison of Elastic design spectra with UHS and response spectra (RS) from IS 1893:2016 and NBC 105:2020 for the bedrock site is shown in Figure 4(b). From Figure 4(b), it can be observed that the elastic design spectra have exceeded the stipulated response spectra of the existing code provisions so were adopted as a target spectrum in the current study.

1D ground response analysis
The study on the effect of local soil conditions to alter the induced ground motion in bedrock for Dharahara was carried out by DEEP SOIL V7 . DEEP SOIL was developed at the University of Illinois and adopts Fast Fourier transform to evaluate the response of soil layers to the input ground motion. The layers and properties of soil were explicitly modelled/provided in the numerical tool with a soil model idealized as lumped MDOF system. The equations of dynamic equilibrium representing the soil model were then solved by Newmark b method. The local soil condition influences the attenuation of the seismic waves and hence modifies the ground motion. The behavior of soil is non-linear and hysteretic during seismic events. The cyclic shear strain induced during seismic events increases the damping ratio (D) and reduces the shear modulus (G) as shown in Figure 5. The hysteresis stress-strain behavior with the initial shear modulus (G max ) and secant shear modulus (G) is shown in Figure 5(a). The generalized shear modulus degradation curve and damping curve (D) are shown in Figure 5(b) and expressed in Equation (1) and (2), respectively.
where W D is the area of the complete hysteresis loop and W S is the equivalent elastic stored energy. Several input parameters required for the numerical tool include the thickness of each layer of soil, unit weight, and shear modulus or shear wave velocity of the material. A suite of site response analyses was performed in the frequency domain and 5% damped response spectral accelerations were calculated. A description of input data consisting of the soil profile and input motion for 1D site response analysis is discussed in the subsequent section.

Soil profile
Standard Penetration Test (SPT) N-value was measured with intervals of 1.5 m up to a depth of 65 m for the three boreholes. All of these three locations consist of sand, clay, dark grey medium dirty clayey sand, silty sandy, clay dark grey silty clay types of soils. The groundwater table was encountered at a depth of approximately 6-8 m from the ground surface. SPT values ranging from 6 to 23 were observed with the majority of test values showing SPT values less than 10. Due to a lack of shear wave velocity (V s ) measurements, in situ standard penetration test blow count measurements were used to determine the shear wave velocity. Gautam (2017) andJICA (2002) developed the empirical relationship between V s and N value for Kathmandu valley (all types of soil). The average V s values given by (3) and (4) were adopted at different depths.
Also, time-average shear wave velocity within 30 m depth V s 30 average shear wave velocity above bedrock V s soil and site natural time period (TG) for different boreholes were calculated based on Equation 5.
where h i and V si are the thickness and shear wave velocity for each considered layer. The calculation of borehole parameters is shown in Table 1. NBC 105:2020 have classified the soil type as type A to D ranging from hard to very soft soil. According to the NBC building code (NBC 105:2020), the soil type in Kathmandu valley is very soft soil, that is, type D. Also, based on the borehole log the soil type observed in the studied location is type D similar to that recommended in NBC 105:2020 for the location under study. According to NBC 105:2020, the sites comprising of low amplitude natural period greater than 1 s can be referred to as very soft soil. From Table 1 it can be observed that site natural time period for all of the investigated borehole logs is greater than 1 s. This investigation shows that there is a high probability of soil amplification around the studied area, which requires a 1D site response analysis to be carried out. The stratigraphy of the boreholes under investigation is shown in Figure 6.
The different soil parameters obtained at the site were found in the following range: specific gravity (G) 2.61 À 2.69 g/cm 3 , plasticity index (PI) 4.4 À 36.1, cohesion (C) 0.75 À 7.5 kPa, and angle of friction 18.5 -29.5 . The index property relationships developed by Darendeli (2001) were used to evaluate the dynamic soil parameter required for analysis. For site response analyses, different variables are required such as overconsolidated ratio (OCR), plastic index (PI), earth pressure coefficient (K 0 ), frequency cycle (N), shear strength, and loading frequency (f). PI, cohesion, and angle of friction were measured in the field. The angle of friction was used to calculate K 0 using the Jaky (1944), shear strength was calculated based on the Mohr-Coulomb criterion (Kramer 1996) while the remaining parameters, that is, OCR, N, and f were taken as 1, 10, and 1, respectively, as recommended by Darendeli (2001). In DEEPSOIL, two models are present to characterize the non-linear dynamic behavior of soil named the General Quadratic Hyperbolic (GQ/H) and Modified Kodner Zelasko (MKZ) model. The GQ/H model requires shear strength as an input parameter while the MKZ model estimates it from empirical equations which result in either overestimation or underestimation of shear strength ). Due to the presence of soft soil deposits in the location under investigation, this study uses the GQ/H model to account for the shear strength because it was reported that soft soil sites may result in an estimation of large shear strains (3-10%) and generate unusual characteristics in the predicted surface ground motion (Griffiths et al. 2016). The bedrock was modelled as an elastic half-space with a V s equal to 760 m/s as this was used in the estimation of uniform hazard spectra for the study site. Figure 7(a) and (b) illustrates the modulus reduction and damping curves incorporated in DEEPSOIL for the site response analyses, respectively.

Input motions
The selection of suitable input ground motions is crucial for the 1D site response analysis. Mostly the ground motion recorded in the studied area comprises records in soft soil, so the ground motions occurring in the seismic environment similar to the region under consideration were selected from the peer ground motion database. The ground motions with similar frequency characteristics to the study region were selected from the PEER database because of a lack of data availability. This selection and scaling criteria are consistent with the studies of , Nguyen et al. (2020), Tran et al. (2021), Kim et al. (2018), and Kiani and Khanmohammadi (2015). Ground motions with shear wave velocity ranging from 760 to 1500 m/s, Moment magnitude (Mw) in the range of 5.5 À 8.5 Mw, and distance ranging from 0 to 100 km were selected for the 1D site response analyses (Niroula & Chamlagain,). Furthermore, the ground motions with Mw > 5 and distance < 100 km are usually investigated for engineering significance (Douglas and Aochi 2008), which motivated these record selection criteria. The selected ground motions having different PGA values ranging from 0.03 to 0.31 g were plotted in Figure 8(a). The required PGA for the studied area was estimated to be 0.424 g. Hence, the time histories were scaled to

Results
This section presents the results of equivalent linear and nonlinear 1D site response analyses. The spectral accelerations, PGA's distribution along with the depth, mobilized shear stress along with the depth, and variation of maximum shear strains for borehole 1 (BH-1) for all the input ground motions are shown in Figures 9, 10, 11, and 12, respectively. In Figure 9, the obtained response spectra (RS) from EQL and NL analysis with code provision (NBC 105:2022 and IS 1893:2016) and target response spectrum for BH-1 are compared. For both types of analyses, deamplification of the input ground motions was observed for short periods. In the case of EQL, the average RS is deamplified up to periods 1.0 s. Similarly, in the case of NL analysis, deamplified zone is observed up to periods of 1.5 s, followed by the amplification of input ground motion, as shown in Figure 9. The zone of deamplification is mostly observed at spectral periods less than the site period. In the case of soft soil, a higher strain level was observed which causes a relatively overdamping value for high frequency. Even though the modifications were performed for large strain values, the results exhibit deamplification for the shorter time period. Moreover, the comparison of response spectra from EQL and NL analysis was also performed with NBC 105:2020 (NBC soil Type D) and IS 1893:2016 (IS soil Type III/soft soil) which are commonly adopted for the seismic design of structures in the study area. It should be noted that there is a considerable difference in the response spectra of NBC soil type D and IS soil type III (IS soft soil). It is also worth mentioning that the response spectra calculated in this study lie in between the NBC soil type D and IS soil type III response spectra. This demonstrates that the results from this study are well in between the envelopes used in the region for soft soil conditions and depict a reasonable response.
The profile of PGA with depth is shown in Figure 10. The PGA is demonstrated to be de-amplified along the depth for both linear and non-linear analyses. For the depth of 10-20 m, larger deamplification was observed from both analyses which are due to the shear wave velocity reversal in the soil profile. The PGA profile from the EQL analysis is slightly smoother than from the NL analysis. This is attributed to the use of non-linear stress-strain response in NL analysis over time during earthquake shaking which minimizes the effect of overdamping. As shown in Figure 10, the PGA value at the surface was obtained to be equal to 0.32 and 0.212 g from EQL and NL analysis, respectively.
The comparison of mobilized shear and maximum strain with depth from EQL and NL analysis for BH-1 is shown in Figure 11 and Figure 12, respectively. Mobilized shear refers to yield stress developed by soil to resist deformation under    applied load. Generally, the mobilized shear increases with an increase in depth which is shown in Figure 11. The strain level significantly affects the soil amplification. Generally, the higher the shear strain higher the damping as a result there are chances of soil deamplification. Also, it can be observed from Figure 10 and Figure 12, that the variation of strain along the depth results in the variation of PGA.
Similarly, spectral acceleration, PGA distribution along with the depth, mobilized shear stress, and variation of maximum strain along the depth for the borehole 2 (BH-2) and borehole 3 (BH-3) are shown in Figure 13, Figure 14, Figure 15, and Figure 16, respectively. In Figure 13, the comparison of obtained response spectra (RS) from EQL and NL analysis with code provision and target response spectrum for BH-2 and BH-3 is compared. For both analyses, PGA was observed to be deamplified like BH-1. The three-bore hole tests comprise of similar soil profile as shown in Figure 6; therefore, the results were observed to be similar. On one hand, the larger amplification in the case of NL analyses at moderate to long periods is attributed to the over-damping of short periods and high frequencies in the EQL analyses. On the other hand, the smaller amplification in the case of NL analyses at short periods is attributed to the incoherence in the input ground motion because of the instantaneous variation of stiffness which is associated with the stress reversal produced in the course of nonlinear stress-strain Cox et al. (2011). As shown in Figure 14 the PGA value at the surface was obtained to be equal to 0.348 and 0.21 g from EQL and NL analysis, respectively for BH-2. Similarly, the PGA value at the surface was obtained to be equal to 0.352 and 0.208 g from EQL and NL analysis, respectively for BH-3.
The variation of mobilized shear strength along the depth is shown in Figure 15. It can be observed that at greater depth shear strength was predicted to be higher from NL analysis as compared to EQL analysis. Also, it can be observed from Figures 14  and 16 that the variation of strain along the depth results in a variation of PGA. It should be mentioned herein that as the seismic waves propagate through soil layers below the ground surface, the local soil and site conditions can affect the amplitude, frequency content, and duration of seismic waves.
Amplification factor versus time was plotted as shown in Figure 17. Figure 17 was plotted for the average amplification of all three borehole logs. There is a zone of deamplifications at periods less than 1.0 s for both EQL and NL analyses. The mean Figure 13. Comparison of obtained response spectra from EQL and NL analysis with code provision and target response spectrum: (a) Borehole-2 (BH-2); (b) Borehole-3 (BH-3).
amplification peaks at 2.8 and 3.0 s in the EQL and NL analyses, respectively. The peaking of amplification beyond site periods is attributed to period elongation characteristics of soft soils in site response analysis.   The short period (F a ) and long period (F v ) seismic site coefficients are expressed in terms of the ratio of response spectral accelerations in the short period and the long period where the interval for F a ranges from 0.1 to 0.4 s and F v ranges from 0.4 to 2.0 s (Borcherdt 2002). These seismic site coefficients were calculated based on (8) and (9) as recommended by Borcherdt (2002). A similar methodology was used by  for the development of Korean site classification and is adopted in the current study.  The calculated mean value of F a is 0.68 and 0.56 for EQL and NL, respectively, whereas F v is 1.22 and 1.0 for EQL and NL, respectively. This demonstrates a lack of short-period amplification in the study area where the mean F a is lower than 1.0 in both sets of analyses. F v also results in nominal values because of long site periods. This is also attributed to the fact that the amplification peaks outside the F v interval, that is, beyond 2.0 s, as illustrated in Figure 17.
After the estimation of seismic site coefficients for the location under investigation, the design spectra were calculated based on NEHRP (FEMA 273 273, 1997) two-point methodology. The comparison of the proposed design response spectrum with the existing response spectrum was performed as shown in Figure 18. From Figure 18 it can be observed that the spectral acceleration adopted by NBC 105: 2020 (soil type D) is higher as compared to the proposed design spectra from the current study for the studied region. Also, from Figure 18 it can be observed that the spectral acceleration adopted by IS 1893:2016 (soft soil) is lower as compared to the proposed design spectra from the current study for the studied region. This trend is consistent throughout the whole period range. This study discusses the site response analysis of the considered site based on EQL and NL analysis methods. Both of the methods have their advantages and disadvantages in terms of estimation of the site response. From the current study, the spectral accelerations were observed to vary from EQL and NL analysis as shown in Figure 18. The result showed that the short period amplification predicted from the NL analysis is smaller which is similar to the result obtained in the research of Griffiths et al. (2016) for soft soil.

Conclusion
This study incorporates seismic hazard analysis followed by site-specific response analysis of the region comprising of monumental structure Dharahara. To preserve the cultural heritage structure and to carry out the risk management plan the vulnerability associated with such structures needs to be assessed. The classical definition of risk addresses human or monetary loss however the value of heritage is intangible and is beyond the direct economic loss. The cultural heritage links the new and old generations of people residing in a specific area in the evolution of human consciousness. Hence, preserving our cultural heritage structures results in preserving our identity. The current research work suggests a seismic assessment methodology for the preservation and risk reduction of the heritage site under earthquake shaking. Such studies are helpful to study the nearby monumental structures and regions where similar hazards and identical ground motions are expected. The performance of such monumental structures during extreme seismic events can also be used as a knowledge transfer of the constructional aspect of such structures to the next generation.
1D site-specific response analysis was carried out in DEEPSOIL v7 using three boreholes. The soil profile was modelled and a suite of equivalent linear and non-linear analyses were performed. To estimate the target response spectrum, PSHA was carried out for the study region. A total of ten different time histories were selected based on target response spectrum, fault mechanism, shear wave velocity, etc., and were provided as input in the numerical tool DEEPSOIL v7. The obtained target response spectrum was compared with existing code provisions, that is, NBC 105:2020 and IS 1893:2016. From the study following conclusions are obtained.
For the study region, response spectrum for bedrock was obtained to be higher from PSHA as compared to NBC 105:2020, and IS 1893:2016. Therefore, during the construction and rehabilitation of structures with historical significance, site response analysis should be carried out. For the simulation results, PGA was observed to vary along with the depth; therefore, soil amplification and its characteristics should be considered during the selection of an appropriate location and type of footing. It is worth mentioning that for a site consisting of very soft soil, spectral acceleration value was demonstrated to be deamplified for short periods up to 1.0 s from EQL analysis and 1.5 s from NL analysis and was amplified afterward. This result shows that the structures having a long period have a severe effect when it is located in a very soft soil region. In the case of soft soil, a higher strain level was observed which causes a relatively over damping value for high frequency and effect lesser for EQL than NL analysis as a result EQL analysis predicted a higher response than NL analysis. Seismic site coefficients F a was demonstrated to have values lower than 1.0 whereas the mean value for F v was calculated to be 1.22 from EQL analysis. The lower values of F a are attributed to the lack of short-period amplification in the region due to long short periods. The lower values of F v are attributed to the intervals used for its calculation. It is recommended to revisit the F v intervals for the study region in the future. It was observed that the design spectrum calculated in this study lies in between NBC 105: 2020 and IS 1893:2016 (soil type D) for the studied region, throughout the whole period range. The design spectra from two seismic codes are demonstrated to be upper and lower bounds in the region under study, which could result in either overestimation or underestimation of the spectral acceleration. Therefore, a site-specific response is recommended.