Abstract

This paper presents use of semiempirical method for seismic hazard zonation. The seismotectonically important region of Uttarakhand Himalaya has been considered in this work. Ruptures along the lineaments in the area identified from tectonic map are modeled deterministically using semi empirical approach given by Midorikawa (1993). This approach makes use of attenuation relation of peak ground acceleration for simulating strong ground motion at any site. Strong motion data collected over a span of three years in this region have been used to develop attenuation relation of peak ground acceleration of limited magnitude and distance applicability. The developed attenuation relation is used in the semi empirical method to predict peak ground acceleration from the modeled rupture planes in the area. A set of values of peak ground acceleration from possible ruptures in the area at the point of investigation is further used to compute probability of exceedance of peak ground acceleration of values 100 and 200 gals. The prepared map shows that regions like Tehri, Chamoli, Almora, Srinagar, Devprayag, Bageshwar, and Pauri fall in a zone of 10% probability of exceedence of peak ground acceleration of value 200 gals.

1. Introduction

Estimation of seismic hazard is an important task for planners, engineers, and scientists. Peak ground acceleration serves as a basic tool for preparation of seismic hazard map in a tectonically active area. Attenuation relations are commonly used for computing peak ground acceleration at any site. Although attenuation relationship provides a good estimate of peak ground acceleration, these relations are mainly controlled by the data used for its derivation. With the advancement of seismology several techniques are used to simulate various strong motion parameters. These techniques are (i) stochastic simulation technique [1], (ii) empirical Green’s Function technique (EGF) [2, 3], (iii) composite fault modeling technique [4], and (iv) semi empirical technique [5]. Recently semi empirical technique has been used for simulating earthquake ground motion due to a buried rectangular rupture plane [5–13].

One of the advantages of this method is its use of attenuation relation of limited magnitude range for forecasting large earthquake. Seismic hazard in an area can be estimated by either the probabilistic seismic hazard assessment approach (PSHA) or the deterministic seismic hazard assessment approach (DSHA) [14]. Both approaches use the same datasets, which include earthquake sources, occurrence frequencies, and ground motion attenuation relationships. Himalayas are one of the most seismically active regions, which is capable of producing large damage due to probable earthquake. However, due to the limited recording of seismic activity in this region, there is no reliable catalogue of strong motion data and suitable attenuation relations which can be used for any seismic hazard study. For such situation the only option which is left is to use the attenuation relation prepared from data set from other regions. Such attempt has been made by different researchers like Khattri [15] and Bhatia et al. [16] for preparing the seismic zonation map for the Indian subcontinent.

The concept of point source approximation of finite earthquake source used in PSHA and DSHA has been one of the major constraints in representation of ground motion due to finite fault in the near-source region. Joshi and Patel [17] have formulated a method of seismic zonation, which is based on the deterministic modeling of finite ruptures along identified probable fault in an area using semi empirical approach. This method of seismic zonation based on semi empirical modeling of rupture plane has been used for preparing seismic hazard map of Doon valley [17], Assam valley [18], and the Uttarakhand Himalaya, India [19]. The technique of hazard zonation is based on the semi empirical simulation technique which in turn is dependent on the attenuation relation applicable in the region under study. In preparation of seismic zonation maps for the Uttarakhand Himalaya by Joshi and Mohan [19] the attenuation relation given by Abrahamson and Litehiser [20] has been used. Although the relation by Abrahamson and Litehiser [20] is based on worldwide data, it does not use any data from Himalayan earthquakes and clearly shows deviation from normality when this relation is used to predict Himalayan data [21]. The method of zonation presented by Joshi and Patel [17] also does not discuss the probability of exceedence of peak ground acceleration at a specific site and offers only a map of expected peak ground acceleration. This paper presents seismic hazard map of the Uttarakhand Himalaya based on modified method of seismic hazard zonation given by Joshi and Patel [17]. The modified method makes use of attenuation relations of limited applicability derived from regional strong motion data from the study area and uses it for preparing seismic hazard zonation map based on probability of exceedence of specific peak ground acceleration.

2. Methodology: Determining the Coefficients of Attenuation Relation

Attenuation relation is one of the major requirements for estimations of seismic hazard of the area. Attenuation relations are generally derived from regression analysis. The first step in the regression analysis is the selection of regression model. In the present work model has been selected on the basis of dependency of peak ground acceleration on magnitude and hypocentral distance. Various functional dependencies have been checked and the model that gives the best correlation and minimum error is selected for final regression model. Most of the earthquakes in the Uttarakhand Himalaya originate at a shallow dipping plane of detachment at approximate depths of about 10–15 km. Therefore in order to have effect of earthquake originating from the plane of detachment a term “” is introduced in the present work which is equal to , which coincides with the approximate depth of plane of detachment. The following regression model has been selected in the present work:
where PGA is horizontal peak ground acceleration in gals, is moment magnitude, R is a hypocentral distance in km, accounts for the earthquake originating at an average basement depth of 15 km and is equal to (), and , and are regression coefficients, respectively. The following set of equations is obtained of different values of peak ground acceleration () obtained at different stations due to earthquake of magnitude () and hypocentral distance ():
For peak ground acceleration data of mth earthquake recorded at kth stations the following set of equations is obtaned:
The above set of equations can be written in a matrix form as follows:
The above matrix equation can be written as follows:
where “d” is the column matrix having peak ground acceleration data at various stations due to different earthquakes, “G” is a rectangular matrix defining the dependency of independent parameters of magnitude and hypocentral distance, and “m” is the unknown matrix having coefficient of regression analysis. The model matrix “m” is obtained by the following expression using least square inversion method:

In this expression, matrix defines the transpose of matrix . In actual case some of the eigen values of “” are very small so that the variance of solution becomes unacceptably large. To avoid this difficulty, we use the damped least squares method of Levenberg [23]. The inversion using damped least square method is expressed as follows:

In this equation, matrix “” is an identity matrix and “” is the damping factor. The criteria of selection of in the present work are based on those given by Dimri [24] and applied by Joshi [25]. Several iterations have been performed in this criteria to select the damping factor corresponding to the minimum root mean square error (RMSE). The root-mean-square error in each iteration has been computed by using the following formula:
where is an observed and is a calculated peak ground acceleration, respectively, and is a sample number.

Attenuation relation of limited applicability developed using method given in (7) can be used for deterministic modeling of rupture plane using semi empirical technique initially proposed by Midorikawa [5] and later modified by Joshi and Midorikawa [11]. This technique requires attenuation relation applicable in the study area. The method was initially proposed by Midorikawa [5] and is used to model rupture plane buried in a homogenous earth model. Modification in this technique has been made by Joshi and Midorikawa [11] to include the effect of the buried plane in the layered earth. This technique is used to simulate accelerogram due to rectangular rupture source and is given in stepwise manner as follows. (i)The rupture responsible for the target event is identified. (ii)Rupture is divided into several subfaults on the basis of self-similarity laws. (iii)Parameters of rupture plane are decided on the basis of earlier studies and source parameters of target events. (iv)One element is treated as a starting point of energy which acts as a nucleation point or the starting point of rupture. The rupture from nucleation point travels in all directions with an assigned rupture velocity. (v)On arrival of rupture front at the centre of an element energy is released in the form of envelope function defined by Kameda and Sugito [26] and modified by Joshi et al. [27] as follows:

In this expression, represents transmission factor [11] and the duration parameters, respectively. The parameter is determined using regional attenuation relation. In the present work PGA is calculated using developed attenuation relation using damped least square method given in (7). (vi)White Gaussian random noise of zero mean and one standard deviation (Figure 1(a)) is passed through several filters given by Boore [28] having basic shape of acceleration spectra given in Figure 1(b). (vii)This is windowed using the envelope function defined in (9) and given in Figure 1(d). (viii)In order to compensate the slip duration of target and small event the windowed filtered record is convolved with a correction function given by Irikura [29] and shown in Figure 1(f). (ix)The geometry of rupture propagation decides the time at which record is released from different elements. (x)The acceleration records released from each subfault reach the observation point depending upon the rupture propagation from the nucleation point to a particular subfault and the travel time of energy in the layered earth from that subfault to the observation point. Appropriate summation of all records gives resultant acceleration record at observation point (Figure 1(g)). The procedure of simulation of strong ground motion is illustrated in Figure 1.

Figure 1: Method of simulation of strong ground motion showing (a) white Gaussian noise, (b) filter representing actual earthquake process, (c) filtered white noise, (d) envelope of accelerogram released by i, j element within the rupture plane shown in (f), (e) multiplication of the envelope with filtered white noise, (f) rupture plane divided into elements in the layered earth with i, j, element releasing record which is convolved with the correction factor F(t), (g) summation of all accelerograms from various elements according to their arrival time at the observation point, and (h) the simulated acceleration record and its acceleration spectrum. Figure modified after Joshi [10].

2.2. Preparation of Seismic Hazard Map Using Semi Empirical Technique

The semi empirical modeling technique gives a chance to model a finite rupture plane for simulation of strong ground motion. This technique has been used by Joshi and Patel [17], Joshi and Mohan [19] for preparing of seismic hazard map. The technique of preparation of seismic hazard map given by Joshi and Patel [17] and Joshi and Mohan [19] lacks representation of probability associated with the observed peak ground acceleration at a different location. In the present work modifications in the technique given by Joshi and Patel [17] have been made to include the role of probability of occurrence of specific of peak ground acceleration value at a given point. The method of preparation of seismic hazard zonation map is given as follows. (i)The first step is identification of active lineaments in the region. The length of a possible rupture along these lineaments is measured from the same map. The length and width of possible ruptures along these lineaments are calculated using the empirical relation given by Wells and Coppersmith [30] and Kanamori and Anderson [31], respectively. (ii)The entire region is divided into a grid consisting of several observation points at which acceleration record is simulated using semi empirical technique. The peak ground acceleration parameter from the simulated record at the observation point is used for preparation of seismic hazard map. (iii)At each observation point peak ground acceleration is stored by modeling, each rupture along selected lineament in the area using semi empirical modeling technique given by Midorikawa [5]. For, “” number of lineaments “” values of peak ground accelerations, that is, , are obtained at that observation point. In this process the peak ground accelerations are also obtained for various possibilities of nucleation point. For a rupture divided into subfaults of size there are possibilities of nucleation point. Therefore the process of simulation generates a data set of peak ground acceleration which consists of all possibilities of ruptures. The database includes contributions from ruptures within 100 km radius from the observation point. The probability of exceedence of peak ground acceleration of 100 gal and 200 gal is computed using the obtained database of PGA values from several models at the observation site. (iv)The frequency magnitude relation for this region shown in Figure 2 is calculated on the basis of available data from the United State Geological Survey (USGS) and is given as follows:
where is the magnitude of earthquake and is number of earthquakes equal to or more than . (v)This process is repeated for all observation points and the probability of exceedence of peak ground acceleration at each point is computed. Contours of the expected acceleration have been used for defining various zones. This parameter is finally used in the preparation of a seismic hazard zonation map.

Figure 2: Frequency magnitude relation for the study area.

For the purpose of preparation of the seismic hazard zonation map, the computer software “EQHAZ” has been developed in the FORTRAN. Required inputs are parameters of the rupture plane and velocity model of the region.

3. Geology and Tectonics of Study Area

The study area lies in the seismically active region of Uttarakhand Himalaya. The northeastern part of the study area is occupied by the Tibetan Plateau while the southwestern sector is the Indo-Gangetic plains. The northeastern part of the study area exposes rocks of the Trans Himalayan tectogen along with late-to posttectonic granitoid and ophiolite bodies [22]. South of this, the belt between the Karakoram Fault and the Indian Suture zone is occupied by cover rocks, affected by Himalayan orogeny, late to post tectonic granitoid ophiolite and lithopackets of the accretionary complex [22]. The foothill belt constitutes foredeep sediments affected by the terminal phases of the Himalayan orogeny. South of foothills the tract is covered by the alluvial fill of the Gangetic foredeep. The northernmost tectonic element is the dextral Karakoram Fault, which is subparallel to the Indus Suture Zone (ISZ), occurring to the south. The Indus Suture Zone (ISZ) is an important tectonic zone, which demarcates the northern boundary of the Indian Plate [22]. Within the main Himalayan belt the high-grade complex of the central crystalline is bound to the north and south by the Martoli Thrust and the Main Central Thrust (MCT), respectively. A similar high grade complex, the Almora crystalline, is delimited on either side by the North Almora Thrust (NAT) and South Almora Thrust (SAT). The main Himalayan belt is represented from the Tertiary frontal fold belt by the Main Boundary Thrust (MBT) and the southern limit of the frontal belt is marked by the Main Frontal Thrust (MFT), which has its surface manifestation at places [22]. A number of N-S to NNE-SSW faults of limited spatial extension have also been mapped from the area. The geology and tectonics of the region are after GSI [22] as shown in Figure 3.

Figure 3: Location of strong motion recorders in the Kumaon and Garhwal Himalaya. The geology and tectonics of the region are after Geological Survey of India [22]. Empty triangles denote the stations maintained by the National Geophysical Research Institute and Department of Earth Sciences, Indian Institute Technology, Roorkee, and half-filled triangles denote stations maintained by the Department of Earthquake Engineering, Indian Institute Technology, Roorkee.

3.1. Case Study: Attenuation Relation for Kumaon and Garhwal Himalaya

The Kumaon region has been selected as an area between latitude 29° and 33° and longitude 80° and 81° which covers 40% of total area of Uttarakhand Himalaya. It was observed that most earthquakes in this part of Himalaya originate from the basement thrust at an average depth of 15 km. Since our database contains very few earthquakes having hypocentral distance less than 15 km, we have rejected that data. Remainder data covers almost 87% of the total data. The hypocentral distance and magnitude range for this data set is between 15 to 100 km and 3.5 ≤ , respectively. Assuming the logarithmic dependency of the magnitude scale on amplitude of ground motion we have assumed three different possible models for regression analysis. These models explain linear, power, and logarithmic dependency of dependent variables on various independent variables like , , and . Using the damped least square inversion scheme regression coefficients have been computed by using (7). The following attenuation relation is selected on the basis of minimum root mean square error and maximum correlation:

Comparison of observed and calculated peak ground acceleration using developed attenuation relation given in (11) is shown in Figure 4.

In the present work the Garhwal Himalaya has been selected as an area between latitude 29° and 33° and longitude 78° and 80° which covers 60% area of the Uttarakhand Himalaya. The data for developing attenuation relation for the Garhwal Himalaya is selected from the events recorded in the strong motion network operating in the Garhwal Himalaya. Peak ground acceleration from 29 records of horizontal component has been used as input data. The distance and magnitude ranges for this data set are between 20 and 210 km and , respectively. The coefficient of attenuation relation is calculated using damped least square method given in (7) and the relation is given as follows:

Comparison of observed and calculated peak ground acceleration using developed attenuation relation for Garhwal Himalaya region given in (12) is shown in Figure 5.

3.2. Simulation Using Semi Empirical Approach: A Discussion

Attenuation relations developed in present work have been used to simulate strong motion records of the Uttarkashi () and the Chamoli () earthquakes. Strong motion records of these two earthquakes have been simulated by Joshi [10], Kumar [13] using attenuation relation developed by Abrahamson and Litihiser [20] in the semi empirical approach. In an attempt to check the applicability of the attenuation relation developed for Garhwal Himalaya in the present work simulation of strong motion record has been made at two near-field stations that recorded these earthquakes.

The rupture models considered for the Uttarkashi and the Chamoli earthquakes are assumed to be similar to those modeled by Joshi [10]. Actual record of the Chamoli earthquake recorded at the Gopeshwar station is shown in Figure 6(a). Figures 6(b) and 6(c) show the records simulated using present relation and Abrahamson and Litehiser’s [20] respectively, in the semi empirical approach. The comparison shows present relation simulates records having peak ground acceleration closer to the actual value. Further comparison of response spectra at 5% damping shown in Figure 6(d) indicates that the obtained curve from simulated record using attenuation relation developed in the present work matches closely with that prepared from observed records. Comparison of actual and simulated acceleration record for the Uttarkashi earthquake at Bhatwari station is shown in Figures 6(e), 6(f), and 6(g), respectively. Simulated records using present relation and that given by Abrahamson and Litehiser’s [20], respectively, in the semi empirical approach is shown in Figures 6(e) and 6(f), respectively. Comparison of response spectra at 5% damping in Figure 6(h) indicates that simulated records using present attenuation relation give a good match in high frequencies as compared to that using attenuation relation given by Abrahamson and Litehiser [20] in semi empirical approach.The comparison of simulated and observed records using semi empirical approach establishes the suitability of developed attenuation relation in modeling earthquake source in the area.

Figure 6: (a) Actual acceleration record, (b) simulated acceleration record using developed attenuation relation for Garhwal Himalaya, (c) simulated acceleration record using the attenuation relation given by Abrahmson and Litehiser [20], and (d) comparison of response spectra at 5% damping from actual and simulated records at Gopeshwar station. (e) Actual acceleration record, (f) simulated acceleration record using developed attenuation relation, (g) simulated acceleration record using the attenuation relation given by Abrahmson and Litehiser [20], and (h) comparison of response spectra at 5% damping from actual and simulated records at Bhatwari station. The black, blue and red lines represent the actual record, simulated record using developed attenuation relation, and that using the attenuation relation given by Abrahamson and Litehiser [20], respectively.

Simulation from semi empirical technique is dependent on rupture dimension and other controlling parameters like starting point of rupture, velocity of rupture, and dip of rupture plane. The PSHA technique uses the attenuation relation to check all possibilities of peak ground accelerations due to an earthquake from identified fault. The attenuation relation gives peak ground acceleration at a particular distance due to an earthquake originating at a point. It is not possible to include the effect of finite fault on peak ground acceleration at that point using only attenuation relation. However such difficulties can be overcome by using modified semi empirical technique. This technique can give several estimates of peak ground acceleration at a particular point by changing various parameters of finite fault. In an attempt to check dependence of peak ground acceleration on location of nucleation point, rupture velocity, and dip of rupture plane, the rupture plane of the Chamoli earthquake on, 29 March, 1999, modeled by Joshi [10] has been used for simulation acceleration record at Gopeshwar station. The nucleation point within this rupture plane is changed iteratively and accelerograms are simulated for each case and are shown in Figure 7. It is seen that peak ground acceleration varies from 176 to 417 gals. Plot of peak ground acceleration obtained from simulated records using different nucleation point in Figure 8 clearly shows drastic variation of peak ground acceleration with respect to mean value. Therefore for estimation of seismic hazard at a particular point all possibilities of nucleation point within rupture plane have been considered in the present work.

Figure 7: Simulated acceleration records at Gopeshwar station for parameters of the Chamoli earthquake on 29 March, 1999, by considering all possibilities of nucleation point within the rupture plane. The shaded rectangle shows the nucleation point which gives maximum ground acceleration.

3.3. Zonation Map of Uttarakhand Himalaya

The area of study forms part of the Kumaon and the Garhwal regions of the Uttarakhand Himalaya within latitude 29°–33°N and longitude 78°–81°E. The main seismic zone extending from Uttarkashi in the west to Dharchula in the east is subparallel to the Himalayan trend and lies in close proximity to Main Central Thrust (MCT). This zone is dominated by shallow focus (0–40 km) events, though some deeper events are also recorded mainly from the Kaurik Fault Zone [22]. Within a period of 181 years from 1816 to 1997, 297 seismic events occurred in this region, out of which 32 events had magnitude [32]. In the present work tectonic map prepared by GSI [22] has been used. A total of 67 lineaments from this map were identified as active lineaments [18]. The possibility of entire or portion of rupture along these lineaments that can generate earthquake of magnitude greater than or equal to 6.0 is visualized in the present work. A total of 111 such possibilities were modeled and are shown in Figure 9. The magnitude of an earthquake generated by rupture along these lineaments is calculated by using the regression relation given by Wells and Coppersmith [30].

Figure 9: Modified tectonic map after GSI [22]. The shaded rectangle shows the rupture dimension and its name used for modeling for seismic hazard map.

The average depth of basement thrust in this region is 12 km. Recent study of Uttarkashi earthquake (1991) and Chamoli earthquake (1999) has shown that the rupture causing these earthquakes originated at a basement at depth of around 12 km. This is one of the reasons to place the rupture responsible for the hypothetical earthquakes at a depth of 12 km. Parameters like strike of the rupture plane are measured from the tectonic map whereas the dip of modeled rupture plane is assumed as 15°. In the present work we are using the attenuation relation having magnitude constraints of ; therefore the magnitude of elementary earthquake is selected as 4.8 which lies well within the constraint limit of the data set used for developing attenuation relation. The number of subfaults required for modeling target events is decided on the basis of self-similarity laws given by Kanamori and Anderson [31]. For all ruptures the rupture velocity is assumed as 2.6 km/sec which is similar to that used for modeling of the Uttarkashi and Chamoli earthquakes by Joshi [6, 10]. The velocity model given by Yu et al. [32] is used in the present work and is selected on the basis of its applicability for modeling of the Uttarkahsi and Chamoli earthquakes by Joshi [6, 8, 10] and Yu et al. [32]. This velocity model is given in Table 1.

The entire region is divided into 100 square grids of length 10 km and each corner of grid is assumed as an observation point. All elements within the modeled rupture plane are treated as nucleation point and peak ground acceleration value from all possibilities of nucleation point is computed from simulated record. A set of peak ground accelerations from different rupture sources has been prepared. This set is used to prepare the probability of exceedence of peak ground acceleration of value 100 and 200 gals, respectively. The prepared zonation map for two different cases is shown in Figures 10 and 11, respectively. It is seen that for the map of peak ground acceleration of 100 gal a large part of area falls in the probability of 10% of exceedence of this peak ground acceleration value. The area becomes smaller when the probability map for exceedence of 200 gals is prepared. It is seen that many locations like Tehri, Chamoli, Almora, Srinagar, Devprayag, Bageshwar, and Pauri fall in a zone of 10% probability of exceedence of peak ground acceleration of value 200 gals. The peak ground acceleration recorded during the Uttarkashi and the Chamoli earthquake at near-source region is between 253 and 352 gals. The destruction during the Uttarkahsi and the Chamoli earthquakes was severe in the near-source region. This indicates that the area which falls under 10% probability of exceedence for peak ground acceleration value of 200 gals in the present study can expect similar destruction in future earthquake and needs to be looked at more carefully.

4. Conclusions

A technique of seismic zonation map based on deterministic modeling of rupture plane is presented in this paper and is applied to the Uttarakhand Himalaya. The prepared seismic zonation map for Uttarakhand Himalaya show that major area falls in 10% and 20% probability of exceedence of peak ground acceleration 100 gals and 200 gals, respectively. This indicates that Uttarakhand Himalaya and the surrounding region are highly vulnerable to seismic hazard. The Tehri, Chamoli, Almora, Srinagar, Devprayag, Bageshwar, and Pauri regions experience peak ground acceleration more than 200 gal. This study explains the utility of the proposed scheme for seismic zonation in any seismically active region where scarcity of digital seismic data poses hindrance to any seismic hazard studies.

Acknowledgments

This work has been done under Indo-Mexico Project no. DST/INT/MEX/RPO-05/2007 sponsored by the Department of Science and Technology, Government of India and Government of Mexico. The authors also thank Indian Institute of Technology, Roorkee, India, for supporting this research work.