Estimation of peak ground acceleration and spectral acceleration for South India with local site effects: probabilistic approach

In this work an attempt has been made to evaluate the seismic hazard of South India (8.0 ◦ N–20 N; 72 E–88 E) based on the probabilistic seismic hazard analysis (PSHA). The earthquake data obtained from different sources were declustered to remove the dependent events. A total of 598 earthquakes of moment magnitude 4 and above were obtained from the study area after declustering, and were considered for further hazard analysis. The seismotectonic map of the study area was prepared by considering the faults, lineaments and the shear zones in the study area which are associated with earthquakes of magnitude 4 and above. For assessing the seismic hazard, the study area was divided into small grids of size 0.1 ×0.1, and the hazard parameters were calculated at the centre of each of these grid cells by considering all the seismic sources with in a radius of 300 km. Rock level peak horizontal acceleration (PHA) and spectral acceleration (SA) values at 1 s corresponding to 10% and 2% probability of exceedance in 50 years have been calculated for all the grid points. The contour maps showing the spatial variation of these values are presented here. Uniform hazard response spectrum (UHRS) at rock level for 5% damping and 10% and 2% probability of exceedance in 50 years were also developed for all the grid points. The peak ground acceleration (PGA) at surface level was calculated for the entire South India for four different site classes. These values can be used to find the PGA values at any site in South India based on site class at that location. Thus, this method can be viewed as a simplified method to evaluate the PGA values at any site in the study area.


Introduction
The Peninsular India is one of the oldest Archaean shield regions in the world, whose continental interiors are generally Correspondence to: T. G. Sitharam (sitharam@civil.iisc.ernet.in) considered as seismically stable. A recent study by Kumar et al. (2007) indicate a lithospheric thickness of 80-100 km for Indian plate. This lithospheric depth is very small when compared to the depth of lithosphere below other continental masses like Africa, Antarctica and Australia. At present Indian plate is moving with a high velocity (Jade, 2004), which is attributed to the low lithospeheric depth as one of the reasons (Kumar et al., 2007). Gangrade and Arora (2000) have also reported that a slow and steady accumulation of seismic energy is occurring in this region which may lead to earthquakes of moderate to significantly high magnitudes. Some of the major earthquakes reported in Peninsular India since 18th Century include Mahabaleshwar (1764), Kutch (1819), Damooh hill (near Jabalpur, 1846), Mount Abu (1848), Coimbatore (1900), Son-Valley (1927), Satpura (1938), Anjar (1956), Koyna (1967), Killari (1993), Jabalpur (1997) and the recent Bhuj earthquake (2001). Mandal et al. (2000) highlighted that during the last four decades stable continental shield region (SCR) of India has experienced moderate seismic activity.
Earthquake hazard is controlled by three factors -source properties, path characteristics, and local site effects. For assessing the seismic hazard the important factors to be considered are: past earthquake data, earthquake source characteristics in the region and the attenuation relationships. Seismic hazard may be analyzed deterministically by considering a particular earthquake scenario, or probabilistically, by considering the uncertainties involved in earthquake size, location, and time of occurrence (Kramer, 1996). In this paper an attempt has been made to estimate seismic hazard at rock level in terms of peak horizontal acceleration (PHA) and spectral acceleration (SA) using probabilistic seismic hazard analysis (PSHA). A new seismotectonic map of the study area has been prepared for this purpose, which consist of the past earthquake data, collected from various sources, and the seismic sources associated with these earthquakes. After declustering the earthquake data, their completeness was analyzed using Stepp's (1972) method and the values of the seismicity parameters, "a" and "b", were estimated using G-R method. Using the seismic hazard parameters and the seismic sources of the region, seismic hazard curves (mean annual rate of exceedance versus PHA (g) and SA (g) for 1 s), rock level response spectrum and contours maps of PHA (g) with 10% and 2% probability of exceedance in 50 years and spectral acceleration corresponding to 1 s and 5% damping with 10% and 2% probability of exceedance in 50 years were prepared.
The peak ground acceleration (PGA) and SA values at ground surface may vary significantly from the values at bed rock level. These variations, either amplification or deamplification, will depend upon the site conditions. For the evaluation of PGA values, the site characterization has to be done and this will require analysis and study of a large amount of geological, seismological and geotechnical data.
The PGA values at the ground surface for the entire study area were evaluated by considering four different site classes viz. site class A to D (BSCC, 2001).
A design response spectrum, which will be extremely useful for designers, has been developed for Mumbai for site class D, based on Eurocode 8 (2005) and IS (BIS-1893(BIS- , 2002 code recommendation. The comparison of design values obtained are also presented in this paper.
The seismic hazard estimation for Peninsular India was previously attempted by Jaiswal and Sinha (2007) and their study area starts from 10 • N. In the present study, the study area selected is from 8 • N and the earthquake data and details of seismic sources were collected for an area of about 300 km from the boundary of study area, as per Regulatory Guide 1. 165 (1997). Apart from this there are major differences in the methodology adopted for the evaluation of seismic hazard. The previous study adopted a zone less approach and the seismic sources were not identified. However in this work an extensive study was carried out to identify and to map the seismic sources from SEISAT (2000) and from remote sensing images (Ganesha Raj and Nijagunappa, 2004). About 189 active seismic sources, which were associated with earthquakes of moment magnitude 4 and above, were identified in the study area. The earthquake catalogue for this study contains 598 events with moment magnitude 4 and above against 184 events considered by Jaiswal and Sinha (2007). In the current work data until 2006 was considered where as Jaiswal and Sinha (2007) had considered the data till 2002. Moreover the present study used a regional attenuation relationship developed by Iyengar and Raghukanth (2004), whereas, Jaiswal and Sinha (2007) gave 50% weightage to this relationship and the rest to other attenuation relations. The variation of the spectral acceleration at 1 s, with 10% and 2% probability of exceedance in 50 years, for the entire study area is also provided in the present study. The variation of peak horizontal acceleration and spectral acceleration has been evaluated for return periods of 475 years and 2500 years. Site response effects were also considered and the PGA values at ground surface for the entire study area were developed for four different site classes. The site response study has not been attempted previously for this region. A design response spectrum for Mumbai, based on site class D, has been developed based on Eurocode 8 (2005) and Indian standards code (BIS-1893(BIS- , 2002.

Study area
South India (8.0 • N-20 • N; 72 • E-88 • E) is considered as a stable continental shield region with moderate seismicity. The seismic regionalization studies started in 1959 by Tandon (1956) and Krishna (1959) to demarcate areas of potential earthquake damage in the Indian subcontinent. The intensity-based mapping of Indian Subcontinent was presented by Guha (1962) and Gubin (1968). The probabilistic seismic hazard studies for this area was done by several researchers including Basu and Nigam (1977), Kaila and Rao (1979), Khatri et al. (1984) and Sinha (2006, 2007). Even though researchers like Reddy (2003), Purnachandra Rao (1999), Gangrade and Arora (2000), etc. have highlighted the need for seismic study of southern Peninsular India, no detailed seismic hazard analysis has been done for the entire South India. The area selected in this study is between latitude 8 • N-20 • N and longitude 72 • E-88 • E, which is shown as the hatched portion in Fig. 1. Seismic study area covers about 300 km from the boundary of study area, as per Regulatory Guide 1. 165 (1997), is also shown in Fig. 1. The study area consists of five states of India, viz. Kerala, Tamil Nadu, Karnataka, Andhra Pradesh, and Goa. The union territory of Pondicherry and some parts of the states of Maharashtra, Orissa and Chhattisgarh also fall in the study area. Total population in South India is about 300 million and there are about eight cities with population exceeding one million in this region. Mumbai, the financial capital of India and the fifth most populous city in the world, alone accounts for a population of 19 million. These facts clearly demonstrate the need to identify the active seismic sources and to assess the seismic hazard in this region considering the local site effects. The study area includes two mountain ranges of Western Ghats and Eastern Ghats and a plateau formed by major rivers like Tungabhadra, Kaveri, Krishna and Godavari. A map showing various geotectonic features of South India (Indian peninsula) like, major fault zones, thrusts, structural trends, granitic intrusions and the various lithotectonic provinces was presented by Valdiya (1973).

Seismicity characteristics
According to the seismic zonation map of BIS-1893BIS- (2002, study area falls in zones II and III. Ramalingeswara Rao (2000) has calculated the strain rate for Indian shield region and the values obtained were found to be the second highest for SCR of the world. Seeber et al. (1999) has reported that the between 1960 and 1990 the seismicity of peninsular India showed a three fold increase. Estimation of seismic hazard values in any region needs the complete details of past earthquakes with a uniform magnitude scale. The earthquake catalogue for this region was prepared by consolidating the available information from different sources and covers the time period 1820-2006. The earthquake data were collected from different sources like Coast and Geodetic Survey, earthquakeinfo (2005), Guaribidanur Array, Geological Survey of India, Indian Meteorological Department, International Seismological Centre, International Seismological Summary, Kalpakkam Atomic Reactor, National Geophysical Research Institute, Seismological Observatory, and United States Geological Survey. In addition to this, a few more data were collected from the catalogues published by different researchers like Oldham (1883), Milne (1911), Turner (1911), Gubin (1968), Kelkar (1968), Gosavi et al. (1977) and intensities of the earthquakes as assigned by Srivastava and Ramachandran (1983). The details of data collected from different sources like, duration of data and the magnitude of data are given in Table 1. The data collected were carefully analyzed to remove the duplicate events. In some cases the magnitude for the same event reported by different agencies were slightly different and in those cases we have selected the higher magnitude along with its epicentral co-ordinates in the analysis. The data obtained were in different scales such as body wave magnitude (m b ) surface wave magnitude (M s ), local magnitude (M L ) and the earthquake intensity scale (I ). Hence these data were converted to a common scale of moment magnitude (M w ). The local magnitudes were converted  Heaton et al. (1986), the body wave and surface wave magnitude were converted using the relations suggested by Scordilis (2006) and the intensity values were converted using the empirical relation (M=(2/3) I +1), where I is the earthquake intensity value. A declustering algorithm was used to remove the dependent events from this catalogue. Most of the declustering algorithms are developed for active tectonic conditions; hence it may not be suitable in the present study. Here we have considered a criterion based on uniform time (>30 days) and space (>30 km) window between successive events. After declustering there were 1588 events in the catalogue, out of which 10 events of M w >=6.0, 135 events of M w >=5.0 and 598 events of M w >=4.0. All the events of M w >=4.0 were used for further analysis. The information on earthquakes with magnitude greater than or equal to 6.0 is given in Table 2.

Regional recurrence relation
The Seismic activity of a region is characterized by the Gutenberg and Richter (1944) earthquake recurrence law. According to this law, Where N is the total number of earthquakes of magnitude M and above in an year and a and b are the seismic parameters of the region. The completeness of the catalogue was evaluated using the Stepp's (1972) method. From the analysis it was found that the magnitude range 4-5 was complete for 70 years, magnitude range of 5-6 was complete for 80 years and magnitudes above 6 were complete for 100 years (Fig. 2). The frequency magnitude relationship obtained from the data for a magnitude interval of 0.1 for the study area is shown in Fig. 3. From this analysis, the earthquake recurrence relation obtained is The hazard parameters are a=4.58 and b=0.891. The "b" value obtained for the study area matches well with the previous studies of Ram and Rathor (1970), Kaila et al. (1972), Ramalingeswara Rao and Sitapathi Rao (1984), Anbazhagan et al. (2009) and Sinha (2006, 2007).

Identification of seismic sources
The seismic sources in the study area were identified from the seismotectonic atlas (SEISAT, 2000), published by the Geological Survey of India. The SEISAT maps are available in A 0 size sheets and each map covers an area of 3 • ×3 • . The study area was covered in 19 sheets of SEISAT. The required pages of SEISAT were scanned and after georeferencing these images, the earthquake data was superimposed on this. The sources, which are associated with earthquake events of magnitude 4 and above, were identified and were used in the hazard analysis. Apart from this some additional lineaments and faults were obtained from the remote sensing image analysis of Ganesha Raj and Nijagunappa (2004). The seismic study area covers about 300 km from the boundary of study area. The seismotectonic map along with the details of earthquakes is shown in Fig. 4. The details regarding location, length, number of earthquakes of magnitude 4 and above occurred along a source and the maximum reported magnitude were collected for all the sources. A total of 189 seismic sources were identified in the study area, which includes 39 faults, 139 lineaments, 2 ridges and 9 shear zones. The length of seismic sources varied form 3 km to 592 km. The details of the source parameters are given in Table 3.

Seismic hazard analysis
Probabilistic seismic hazard analysis (PSHA) was initially developed by Cornell (1968). Many researchers have adopted this methodology for evaluating hazard and recently this method has been adopted by Iyengar and Ghosh (2004), Raghu Kanth and Iyengar (2006) and Anbazhagan et al. (2009) for the probabilistic seismic hazard analysis of Delhi, Mumbai and Bangalore respectively. In the present study also similar methods are used to evaluate the seismic hazard of South India. The hazard curves obtained from PSHA show the variation of peak horizontal acceleration (PHA) or spectral acceleration (SA) against mean annual rate of exceedance. For calculating the seismic hazard values, the entire study area was divided into grids of size 0.1 • ×0.1 • (about 9500 cells) and the hazard values at the centre of each grid cell were calculated by considering all the seismic sources with in a radius of 300 km. For this purpose a new program was developed in MATLAB.
The occurrence of an earthquake in a seismic source is assumed to follow a Poisson's distribution. The probability of ground motion parameter at a given site, Z, will exceed a specified level, z, during a specified time, T is represented by the expression: Where ν (z) is the mean annual rate of exceedance of ground motion parameter, Z, with respect to z. The function ν (z) incorporates the uncertainty in time, size and location of future earthquakes and uncertainty in the level of ground motion they produce at the site. It is given by: Where N n (m 0 ) is the frequency of earthquakes on a seismic source n, having a magnitude higher than a minimum magnitude m 0 (in this study it is taken as 4.0). f n (m) is the probability density function for a minimum magnitude of m 0 and a maximum magnitude of m u ; f n (r|m) is the conditional probability density function (probability of occurrence of an earthquake of magnitude m at a distance r from the site for a seismic source n); P (Z>z|m, r) is the probability at which the ground motion parameter Z exceeds a predefined value of z, when an earthquake of magnitude m occurs at a distance of r from the site. The integral in Eq. (4) can be replaced by summation and the density functions f n (m) and f n (r|m) can be replaced by discrete mass functions. The resulting expression for ν (z) is given by: Where λ n (m i ) is the frequency of occurrence of magnitude m i at the source n obtained by discretizing the earthquake recurrence relationship for the source n. The magnitude recurrence model for a seismic source specifies the frequency of seismic events of various sizes per year. The seismic parameters were determined using Gutenberg-Richter (G-R) magnitude-frequency relationship which is given in Eq. (2). The maximum magnitude (m u ) for each source was estimated based on maximum reported magnitude at the source plus 0.5 (Iyengar and Ghosh, 2004;Raghu Kanth and Iyengar, 2006). The recurrence relation for each fault, capable of producing earthquake magnitudes in the range m 0 to m u was calculated using the truncated exponential recurrence model developed by McGuire and Arabasz (1990), and it is given by the following expression: Where β=bln (10) and N i (m 0 ) weightage factor for a particular source based on deaggregation. Deaggregation procedure followed here is similar to the one followed by Iyengar and Ghosh (2004), Raghu Kanth and Iyengar (2006) and The weigthing factor for length, L i is the length of the fault i. The earthquake event weighting factor (χ i ) has been taken as the ratio of the past events associated with source i to the total number of events in the region as given below: χ i = Number of earthquakes close to the source Total number of earthquakes in the region (9) The recurrence relation of source i has been calculated by averaging both weighting factors as given below: The value of m 0 for all the faults were taken as 4 and the value of m u for each of the faults were calculated as explained above. The magnitude range for each fault was divided into small intervals and the recurrence rates of each of these magnitudes were calculated using the above equations. In PSHA, the probability of occurrence of an earthquake anywhere in the fault is assumed to be the same. The uncertainty involved in the source to site distance is described by a probability density function. The shortest and longest distance from each source to the centre of the grid is calculated. The hypocentral distance was calculated by considering a focal depth of 15 km, similar to the one used for DSHA by Sitharam et al. (2006) and Sitharam and Anbazhagan (2007). A probable source zone depth of 10 km has been considered by Bhatia et al. (1997) in an exercise to develop seismic hazard map of the shield region of India in GSHAP. The conditional probability distribution function of the hypocentral distance R for an earthquake of magnitude M=m is assumed to be uniformly distributed along a fault. It is given by Kiureghian and Ang (1977) as: Where X(m) is the rupture length in kilometers, for an event of magnitude m, estimated using the Wells and Coppersmith (1994) equation, which is as given below: X(m) = min 10 (−2.44+0.59(m i )) , fault length (14) The notations used in the Eqs. (11)(12)(13)(14) are explained in Fig. 5. Since the study area is located in peninsular India, the attenuation relation (for peak horizontal acceleration and spectral acceleration) for the rock site in Peninsular India developed by Raghu Kanth and Iyengar (2007) was used. The attenuation relation suggested for the region is: Where y, M, R and ε refer to PHA/spectral acceleration (g) at the bed rock level, moment magnitude, hypocentral distance and error associated with the regression respectively. The coefficients in Eq. (15), c 1 , c 2 , c 3 , and c 4 are obtained from Raghu Kanth and Iyengar (2007). The normal cumulative distribution function has a value which is most efficiently expressed in terms of the standard normal variables (z) which can be computed for any random variables using transformation as given below (Kramer, 1996): Where PHA is the various targeted peak acceleration levels which will be exceeded. ln PHA is the value calculated using attenuation relationship equation and σ ln PHA is the uncertainty in the attenuation relation expressed by the standard deviation.

Evaluation of seismic hazard
For generating seismic hazard curves, a newly developed MATLAB program was used. The input for this program include the parameters of the faults, viz. latitude and longitude of the starting and end point of the faults, number of earthquakes along the fault, the maximum magnitude, the seismicity parameters and the grid size. The program will divide the study area into grids of specified size and the seismic hazard values at the centre of each grid cell will be calculated. The hazard curves for mean annual rate of exceedance for PHA and Spectral acceleration for 1 s at rock level were calculated for each of the grid cells. Uniform hazard response spectrum (UHRS) at rock level with 2% and 10% probability of exceedance in 50 years were also calculated for all the grid points. UHRS with 2% and 10% probability of exceedance in 50 years for important cities in South India are shown in Figs. 6 and 7.

Local site effects
In the previous sections the evaluation of PHA and SA values at rock levels were discussed in detail. Due to the variations in the local site conditions, the surface level peak ground acceleration (PGA) and spectral acceleration values can be different from the rock level values. The method which is widely followed for classifying the site is the average shear wave velocity in the top 30 m, which is commonly referred to as V 30 s . This method has been adopted by many codes for estimating the design ground motion at a site, after incorporating the local site conditions. In the present study four site classes were considered, viz. site class A to D. The shear wave velocity ranges for each site class, as suggested by National Earth- ln F s = a 1 y br + a 2 + ln δ s Where a 1 and a 2 are regression coefficients, y br is the spectral acceleration at rock level and δ s is the error term. The value of spectral acceleration for different site classes can be obtained using Eq. (18).
Where y s is the spectral acceleration at the ground surface for a given site class.

Results and discussions
In this study an attempt has been made to evaluate the seismic hazard based on PHA and SA at rock level and PGA at ground level, based on different site classes, for South India based on probabilistic seismic hazard analysis. Using a newly developed MATLAB program, the hazard curves and UHRS for 2% and 10% probability exceedance in 50 years are calculated for 9500 grid points in the study area. The program was validated by comparing with results obtained by Raghu Kanth and Iyengar (2006) for Mumbai. The values obtained in this study matches well with the previous results. To compare the seismic hazard for important cities in South India, hazard curves for mean annual rate of exceedance versus PHA and mean annual rate of exceedance versus SA for 1 s (both at rock level) are given in Figs. 8 and 9. The rock level PHA distribution for the entire study area with 10% probability of exceedance in 50 years (return period of 475 years) is shown in Fig. 10. The maximum PHA value of 0.37 g was estimated at the Koyna region. A PHA value of 0.37 g corresponds to zone V as per the current seismic hazard map of India (BIS-1893(BIS- , 2002. However, the Koyna region is placed in the seismic zone IV, where the PHA value varies for 0.16 g to 0.24 g (Fig. 11). The PHA values in the central region of South India, around Bangalore, are in the range of 0.11 g to 0.16 g. This region also shows higher seismic hazard than the code (BIS-1893(BIS- , 2002) specification. The PHA value for the Bangalore region reported by Sitharam and Anbazhagan (2007) is 0.146 g, which also matches with the present study. Similar results were obtained in the studies by Jaiswal and Sinha (2007)   and Anbazhagan et al. (2009) for Bangalore region. The variation of spectral acceleration at 1 s for 10-percentage probability of exceedance in 50 years is shown in Fig. 12. Here also the maximum SA at rock level is obtained at Koyna region.
In the recent years the building codes are revised to consider the 2% probability of exceedance in 50 years, which corresponds to a return period of 2475 years. The seismic hazard curves for PHA at rock level for 2% probability of exceedance in 50 years is given in Fig. 13. The increase in PHA values corresponding to a decrease in the probability of exceedance is very small for most of the regions in South India. In these hazard curves also the Koyna region is having the highest PHA value of 0.66 g. The spectral acceleration at 1 s with 2% probability of exceedance for 50 years at rock level is given in Fig. 14.
Uniform hazard response spectra (UHRS) were developed for bedrock and A, B, C and D type soil classes for Mumbai and Bangalore. For both the cities the probabilities of ex-  ceedance considered were 2% and 10% in 50 years and the results are shown in Figs. 15 to 18. These results clearly show the variation of predominant frequency with change in soil types. For the cities of Mumbai and Bangalore, the period of oscillation corresponding to maximum SA varies from 0.05 s at bedrock level to 0.2 s at ground surface for site class D. The knowledge of soil condition is essential for determining the PGA at ground surface because the amplification of spectral acceleration values varies with site conditions. However for a large area like South India, which covers more than one million square kilometers, getting the     Fig. 24. This result shows that the Eurocode method is giving higher values when compared to the IS code values. In Eurocode the soil classification is done based on V 30 s , where as the IS code classification is given as dense, medium and soft soil based on SPT-N values.

2002) for site class D and is shown in
Even though we have followed an entirely different methodology for the seismic hazard analysis, the pattern of variation of PHA values at rock level for a return period of 475 years obtained from this study matches well with the results obtained by Jaiswal and Sinha (2007). The slightly higher value of PHA obtained in the present study is mainly due to the regional attenuation relation used in this study (Iyengar and Raghu Kanth, 2004), the different methodology followed in identifying seismic sources and the modified earthquake catalogue used. More over this study estimates  the SA for 1s and the PGA values at the ground surface level for the entire study area. Response spectrum for different site classes for Mumbai and Bangalore are also presented here.

Conclusions
In the present study an attempt has been made to evaluate the seismic hazard of South India using probabilistic approach. Seismic hazard was evaluated using the earthquake data obtained till December 2006 and the details of seismic sources were obtained form SEISAT-2000 and current literature. The variation of PHA as well as SA for 1 s were estimated for 2% and 10% probability of exceedance in 50 years. It is observed that the difference in PHA values and the SA values for the probability of exceedance of 2% and 10% is very less  for major portion of the study area. Further, local site effects have been considered and a simple methodology is presented to estimate ground level response spectrum for any local site based on PSHA and site class. Based on the present investigation, the following conclusions have been reached.
1. The PHA values obtained for some of the regions are much higher than what is specified in the seismic zoning map of India (BIS-1893(BIS- , 2002. This indicates that more detailed study of the area is required and the seismic zoning map of IS (BIS-1893(BIS- , 2002) code has to be revised.
2. The PHA, SA and UHRS values are developed by incorporating the uncertainties in magnitude, hypo central distance and attenuation of seismic waves.

Site class -A Site class -B
Site class -C Site class -D 3. The seismotectonic map was prepared by considering all the known seismic sources which are associated with earthquakes of magnitude of 4 and above. This can be used for the microzontation studies of different parts of South India.
4. The period of oscillation for maximum spectral acceleration for the study area (at bedrock level) was found to be less than 0.1 s.

The UHRS for various site conditions at Mumbai and
Bangalore clearly shows the variation in predominant frequency and spectral acceleration, with change in soil conditions.
6. This study considers published seismic sources from SEISAT (2000) and remote sensing studies by Ganesha Raj and Nijagunappa (2004). The identification of new seismic sources will have to be carried out using recent satellite data and this work is under progress. These new sources need to be further used for revising the hazard.
7. The PGA value for the entire South India has been evaluated for different site classes in this study. By knowing the site class at any location in the study area, the ground level PGA values can be obtained based on the simplified methodology adopted in this study.