Seismic Response of Sri Lanka using PSHA Technique

This paper describes a research study carried out to find the seismic response of Sri Lanka using Probabilistic Seismic Hazard Assessment (PSHA) considering seismic activities surrounding the region. Three attenuation relationships, namely Abrahamson and Silva, Raghukanth and Iyengar, and Campbell and Bozorgnia, verified to be reliable according to existing South Indian seismic data have been used in the analysis. The completeness check for the earthquake catalogue for Sri Lanka seismic scenario has been performed and 200 years of data were demonstrated to be satisfactory. PSHA technique was used with logic tree approach to minimize the effects of epidemic uncertainty. Then, Peak Ground Accelerations and Spectral Accelerations at bedrock level ignoring the overburden effects were computed for important cites in Sri Lanka and they were found to be highest in the west coastal zone within a range of 0.05 to 0.1 g and an approximate PGA value of 0.1 g in Colombo.


Introduction
Sri Lanka is located within a tectonic plate known as "Indo-Australia plate". According to the historical records, very few number of earthquakes were recorded within the country. The seismic scenario of Sri Lanka in reference to oceanic earthquakes in NW to SW region has been described in detail by Seneviratne et al. (2019) [22]. At present there is little evidence to identify seismic hazards originating from the oceanic region from SW to NE. With the Sri Lankan landmass being solidly connected to South Indian plate, the impact of earthquakes originating from South India may be of importance. However, there has been little seismic activities in the southern oceanic region though some scholars [8] have pointed out the possibility of an oceanic plate boundary forming to the deep south of Sri Lanka. Impact of seismic activities around Sri Lanka may be predicted using Probabilistic Seismic Hazard Assessment (PSHA) technique. This paper describes a study of seismic response of Sri Lanka considering all surrounding seismic scenarios. The overall results of PSHA are compared herewith those obtained from Deterministic Seismic Hazard Assessment (DSHA) [22].

Probabilistic Seismic Hazard Assessment (PSHA) Technique
DSHA and PSHA are the commonly used techniques for the determination of seismic hazards. The procedure of DSHA technique has already been described by Reiter (1990) [20].
The basic procedure of PSHA was initially developed by Cornell (1968) [7] and its computerized form was implemented by Mc Guire (1976 and 1978) [10]. Modern PSHA has gradually evolved by incorporating additional terms and computational tools in order to represent seismic hazards more accurately. The basic methodology involves calculating how often a suite of specified levels of ground motion will be exceeded at the site. The general procedure for a Cornell-McGuire PSHA comprises four fundamental steps.
The first step involves the identification and delineation of all potential sources of seismicity that may affect the site of interest.
These sources of seismicity may be represented as area sources, fault sources, or point sources, depending upon the geological nature of the source and available data.
In the second step, the temporal behaviour of earthquakes is assumed to follow a Poissonian process and it is determined for each source by establishing a magnitude recurrence relationship over the range of magnitudes that are likely to be generated by each seismic source.
The third step involves the use of Ground Motion Prediction Equations (GMPEs) which are same as predictive relationships in DSHA, to establish the conditional probability of exceedance of a pre-specified ground motion value for each site given the occurrence of an earthquake of a particular magnitude and location. The final step of the analysis computes the annual number of events that produces a ground motion parameter, e.g. Spectral Accelerations (SA), that exceeds a specified threshold level, z. This number of events per year is also called the "annual frequency of exceedance", and the inverse of it is called the "return period". Several probability distributions for each seismic source defined in the previous steps are determined from the past earthquake records.
This catalogue is shown in Figure 1 and a filtered version of this catalogue was used to describe the seismic scenario for oceanic earthquakes in the NW-SW in the region of Sri Lanka by Senevitatne et al. (2019) [22]. The completeness of the original catalogue has been verified by the authors following the procedure recommended by Stepp (1973) [23]. The above results shown in Figure 2 illustrate the completeness of the catalogue beyond 200 years.

Macrozonation studies
Various researchers have studied the seismicity and tectonics within and outside of Sri Lanka in order to highlight the potential seismic hazard. Many Indian scientists have carried out studies to define seismic characteristics of peninsular India. Menon et al (2010) [12] had carried out a probabilistic seismic hazard macro-zonation of Tamil Nadu in Southern India in 2010. They have scientifically identified eleven seismogenic zones related to the Southern part of India around Tamil Nadu including Sri Lanka.
Based on above studies, a macro-zonation has been proposed by the authors to study seismic response of Sri Lanka as shown in Figure 3. The same seismogenic zones used by Menon et al. (2010) [12] has been used in this study except the Sri Lankan zone which is replaced by the line source of Comorin Ridge and Mannar lineament, the dominant seismic features of the zone. This is justified as very few earthquake records are available in the sea region from SW to NE; also according to Vitanage (1995) [24] and Fernando & Kulasinghe (1986) [9] only micro seismic activities had originated within the country.

Gutenberg & Richter Relationships
Following the steps of PSHA technique the probability of occurrence of earthquakes within each area source (macro zone) and the line source were determined from the earthquake catalogue [14]. As an example, Figure 4 defines the above relationship for the Comorin Ridge and Mannar lineament line source. The complete relationships (gradient and intercept) for all the sources determined from semi logarithmic plots of magnitude verses inverse of return period are given in Table 1.

Selection of Attenuation Relationships
Attenuation relationships are used to determine the seismic response at a site within one area source due to an earthquake in another area source. It is assumed that an earthquake in another area source will occur at the point on the boundary between the two areas along the line which connects the area source. Therefore, the shortest distance to the area source boundary was used as epicentral distance in these attenuation relationships to calculate the Peak Ground Accelerations (PGAs) and Spectral Accelerations (SAs).
Three ground motion attenuation relationships were chosen in this study to derive the PGA and the SA values. The relationship by Abrahamson and Silva (1997) [3] was developed from earthquakes recorded worldwide, especially in the West and North America.
The above relationship is recommended for use in the analysis of shallow crustal seismic events worldwide. Raghukanth and Iyengar (2007) [18] developed an empirical relationship for Peninsular India based on a stochastic seismological model and compared its predictions with instrumental data from Koyona (1967) and Bhuj (2001) earthquakes in India. In addition to the above two attenuation relationships, the attenuation relationship developed by Campbell and Bozorgnia (2008) [5] which is validated for western United States and similar tectonically active regions of shallow crustal faulting has been chosen for the current study. The selected attenuation relationships to predict the median response are: (i). Abrahamson & Silva (1997) [3]: Base form of the magnitude and distance dependence of the attenuation relationship is only considered in this study ignoring the effects of the fault type and sites located on the hanging wall.  12 Where c0, c1, c2, c3, c4, c5 and c6 are the coefficients given in Campbell and Bozorgnia (2008) for the period range of 0.01 second to 7.5 second. M and RRUP are the moment magnitude and the closest distance to the rupture plane, respectively.
In the above equations, Sa is the Spectral Acceleration and g, the gravitational acceleration.
Attenuation of ground motion predicted by these three relationships is compared with strong motion records of two earthquakes from Peninsular India. Four strong motion records from the Jabalpur earthquake (1997) and Bhuj aftershock (2001), both of Magnitude Mw 5.7, have been used in the comparison of predicted and observed attenuation characteristics ( Figure 5). As seen from the figure, predicted and observed behavior is similar and therefore, these attenuation relationships are used for subsequent analyses with equal weightage. Figure 6 illustrates the spectral accelerations predicted by the three selected attenuation relationships at 1 second. As seen from the figure, the predicted ground accelerations from the three attenuation relationships lie very close to each other at all epicentral distances of interest (10 -500 km). Similar agreement was also obtained at 0.5 seconds. where c1, c2, c3 and c4 are coefficients of attenuation equation for Southern India given in Raghukanth and Iyengar (2007). M and r refer to moment magnitude and hypocentral distance respectively. ln(br) is the error term. It is zero for the median response.

( )
Where c0, c1, c2, c3, c4, c5 and c6 are the coefficients given in Campbell and Bozorgnia (2008) for the period range of 0.01 second to 7.5 second. M and RRUP are the moment magnitude and the closest distance to the rupture plane, respectively.
In the above equations, Sa is the Spectral Acceleration and g, the gravitational acceleration.
Attenuation of ground motion predicted by these three relationships is compared with strong motion records of two earthquakes from Peninsular India. Four strong motion records from the Jabalpur earthquake (1997) and Bhuj aftershock (2001), both of Magnitude Mw 5.7, have been used in the comparison of predicted and observed attenuation characteristics ( Figure 5). As seen from the figure predicted and observed behavior is similar and therefore, these attenuation relationships are used for subsequent analyses with equal weightage.  Figure 6 illustrates the spectral accelerations predicted by the three selected attenuation relationships at 1 second. As seen from the figure, the predicted ground accelerations from the three attenuation relationships lie very close to each other at all epicentral distances of interest (10 -500 km). Similar agreement was also obtained at 0.5 seconds.

Logic Trees
Logic tree methodology was used to address the epistemic uncertainty of various parameters and relationships in the PSHA calculation. The maximum cut off magnitude is based on the maximum historical earthquake (MHE) in each source zone of the earthquake catalogue. As an alternative, the MHE increased by 0.3 units has been considered. 60% of weight has been ( ) been considered. 60% of weight has been assigned to the maximum magnitude and balance 40% of weight was assigned to MHE increased by 0.3 units. For the selected attenuation relationships already described in previous section, it was difficult to assign a higher weight to one equation over the other due to closeness of prediction. Hence equal weightages have been assigned to all three relationships. Figure 7 shows this logic tree probability distribution.

Calculation Procedure
For each of the point of interest in the site, say point A (e.g. Colombo), the effect of the seismic activity in the surrounding source zones was modeled by placing the epicenters within the source zone but at the point closest to the point A. The magnitude of the earthquake was taken according to the return period considered (for e.g. 475 years) and the resulting spectral acceleration at A was calculated using the attenuation relationships and the logic trees described in Sections 6 and 7. The resultant bedrock response spectrum at point A was taken as the envelope of all the response spectra obtained in the above manner. PGA at point A is the intercept of response spectrum at zero period.  Figure 11 illustrates the calculated 5% elastic damping Response Spectra (RS) for 50 year, 475 year and 2475 year return periods at bed rock level in Colombo. Similarly, the response spectra at other major cities in Sri Lanka were calculated.

Conclusions
In this study, an improved seismic hazard assessment is performed for Sri Lanka incorporating new earthquake catalogue and zonation approach, ground motion prediction equations, and treatment of the inherent uncertainties. The new catalogue is compiled for area bounded by latitudes 0 0 N to 20 0 N and longitudes 70 0 E to 90 0 E based on historical and instrumental earthquake data from 1063 to 2012. For the new zonation approach, eleven seismogenic area source zones and one line source were defined. The GMPEs are validated with the past earthquake data for their predictions in the South Indian region. Epistemic uncertainty in the computed seismic hazard is accounted for within a logic tree framework with alternatives for time completeness analyses, maximum cut-off magnitude and GMPEs.
Peak Ground Acceleration (PGA) for rocky or hard soil sites in Sri Lanka varies from 0.1g in West coastal line to 0.05g in the East coastal line for the 475-year return period event while it is fairly constant in the North-South direction. This is mainly due to the fact that effect of line source of Comorin Ridge and Mannar lineament, the dominant seismic features of the zone, becomes dominant in the PGA map (as well as spectral acceleration maps) for the 475 return period event. Also, 5% elastic damping response spectrum was proposed for rocky or hard soil sites based on the results obtained.
Furthermore, the present results indicate that the hazard distribution in Sri Lanka is significantly higher than that specified previously by Abayakoon (1996) and Peiris (2007). This may due to the use of GMPEs as appropriate for different seismotectonic regimes for which multiple models including a line source representing epicentres of past earthquakes in Comorin Ridge and Mannar lineament were considered.