Determination of Response Spectra for Sri Lankan Cities using Finite Difference Method

Sri Lanka has been considered as a seismically safe country in the past considering the large distance from the island to active plate boundaries. However, with the increased degree of urbanization, the possible impact of intraplate earthquakes on population centres within the island has become important. In this context, deterministic and probabilistic seismic hazard assessments have been carried out by few researchers so far. However, these studies have not considered the influence of variation of bedrock profile on the seismic wave propagation. In the study presented here, a numerical simulation is carried out to investigate the effect of variation of the bedrock profile on the seismic wave propagation in Sri Lanka. The acceleration time histories of seven real time earthquake records, selected from the PEER database, were used as input in a numerical finite difference model simulating the two-dimensional bedrock response. The resultant spectra at the bedrock surface thus obtained are compared with those derived from the early studies. The comparison clearly shows the effect of ground elevation profile on seismic wave propagation, which should be taken into account especially when designing earthquake resistant structures in higher elevation areas.


Introduction
Sri Lanka, an island situated in the Indian ocean close to South India, has been considered as a safe country, especially with related to earthquakes and associated hazards, by the general public and most of the authorities. As a result, there are no proper guidelines for earthquake resistant structures. Even though there were few seismographs installed across the country, the primary focus has been to contribute to global networks assisting the prediction of major seismic events within the region. The instruments required to monitor the seismic response required for engineering purposes are not installed even now. This paper describes an analysis carried out to determine the bedrock response at major cities in Sri Lanka due to seismic activity in the region of Comorin ridge and failed Mannar rift zone [6] (Seneviratne et al. 2019) [6]. The absence of real data from earthquakes in proximity to the study area led to the use of international data bases by the researchers.

Seismic Zone Considered in the Analysis
There are no records of major earthquakes in or around Sri Lanka in the recorded history. The most significant seismic event to affect Sri Lanka had occurred in 1615 in close proximity to Colombo causing nearly 2000 causalities. Since then, except for ground shaking, no causalities have been recorded though few tremors have been experienced from time to time.
However, a careful examination of the past records available indicates that there is concentration of seismicity in proximity to Comorin Ridge and failed Mannar rift zone on the west coast of Sri Lanka. The seismic activities in the SW to NE arc appear to be significant. This scenario has been examined in depth by Seneviratne et al. (2019) [6]. Figure 1 illustrates the seismicity of Comorin ridge and failed Mannar rift zone since 1900. Gutenberg and Richter (1954) [3] relationship indicates a magnitude of 6.9 on the Richter scale at a return period of 475 years. This is described in detailed by Seneviratne et al. (2019) [6]. The analyses described in this paper are based on M 6.9 earthquake occurring at an epicentral distance of 90 km from the western coast line as suggested by Seneviratne et al. (2019) [6].

Figure 1 -Epicentres of Earthquakes Occurred on the Comorin Ridge and the Surrounding Area
This position of the epicenter lies along the Comorin ridge and failed Mannar rift zone which also contain the actual epicenter of the 1938 earthquake, considered to be the most significant earthquake to seismically affect Sri Lanka since 1900.

Formulation of the Problem
The seismic waves generated at the hypocentre of an earthquake propagate outward along the bedrock and reaches the continental shelf of the west coast of Sri Lanka. According to Sreejith et al. (2008) [7]., the bedrock level rises from Comorin ridge towards the Continental shelf by about 3 km with most of the rise occurring close to the continental shelf as shown in Figure 2.
In the analyses conducted here, the seismic activity from the focus to the continental shelf is not modelled. By using a design earthquake, the time history records of seismic waves arriving at the continental shelf were assumed.
Five cross sections running in SW to NE direction passing through Colombo-Batticaloa, Kalutara-Kalmunai, Galle-Arugambay, Chilaw-Vakarei and Puttulum-Trincomalee were modelled in the analyses. Figure 3 illustrates the cross sections. These cross sections were selected covering the main cities in Sri Lanka to analyse the effect of topography along the ground profile on the seismic wave propagation.

Finite Difference Analysis
Finite difference analysis was performed along the above five cross sections which are parallel to each other and aligned SW-NE with an interval of 100 km. Especially in FLAC, the explicit Lagrangian calculation scheme and the mixeddiscretization zoning technique is used to ensure that plastic collapse and flow are modelled accurately. Another advantage of FLAC is that, as no matrices are formed during the analysis without excessive memory requirements of the computer, large twodimensional calculations can be made. As the computational time of a large finite difference model of this nature is very high, a personal computer having 16 GB RAM, with a processor of core i7 was used for the analysis.

Model Development
In the literature, there are many examples of analytical model studies based on wave propagation carried out as part of micro seismic studies. However, the authors are not aware of any analysis carried out on macro seismic scale such as in the present study. In the absence of proper guidelines available in the literature, the development of the present model was carried out on an explorative basis. Checks were conducted wherever possible to ensure the reliability of the results obtained.
Modelling of seismic wave propagation due to a distant seismic event, starting from the epicentre to any point of interest, would involve a large computing time and cost. Therefore, to simplify the problem to a realistic level, the arrival wave at the coastal line of Sri Lanka was assumed as input to the model, and the wave propagation behaviour from west coast to east coast was analysed. Figure 4 illustrates the model cross section developed for the section passing through Colombo-Batticaloa. The ground elevation profile was obtained using elevation from the mean sea level at 100 m intervals along the longitudinal direction. The levels were obtained from 1:50000 maps of the Survey Department of Sri Lanka. It was assumed that the depth of overburden is negligible, and therefore the bedrock level and the ground level were taken as the same. This is justified as the depth of overburden is rarely more than 50 m and the effect of it on the bedrock response to a seismic event is likely to be relatively insignificant.

Figure 4 -Model Cross Section Developed for the Section Passing through Colombo-Batticaloa
The cross section was discretised using a step size less than one tenth of the shortest wave length of the significant seismic wave spectrum (0.1 Hz to 25 Hz). This was fine tuned after a series of trial runs minimizing numerical difficulties at the shortest running time. The model was developed with the west and east coastal lines at either end. According to the past earthquake records, the epicentral depth is in the range of 10-15 km. Therefore, a depth of 15 km from the mean sea level was selected as the height of the model.
The length of the model sections varied from 200 km to 250 km depending on the cross section. Each cross section passes through two cities on either end. The cities were decided based on the importance of the city from administrative point of view and the elevation. As one of the objectives of this study is to have a look at the effect of ground profile on wave propagation, the height of a city from the mean sea level was taken into consideration.

Material Properties
According to geology of Sri Lanka described by Cooray (1984) [1], over 80% of land area of Sri Lanka is covered by ancient crystalline rocks composed of metamorphose sediments of Precambrian era intruded by granite and pegmatite. The north and north west coastal belts are underlain by sedimentary rocks, mainly limestone. Most of the rock types exhibit joint fracture folding and faulting.
The effect of differences in bedrock type and structure are not considered in the present analysis. The bedrock is assumed to be homogenous and isotropic. The bedrock properties used in the analysis are given in Table 1. Dynamic elastic modulus, dynamic Poisson's ratio, bulk density and shear wave velocity are taken as the average values of those given by Jayawardena (2001) [4] based on testing of a large number of borehole cores and irregular rock of different rock types taken from various locations in Sri Lanka. Shear modulus (G) and shear wave velocity (Vs) were calculated from Elastic modulus (E) and Poisson's ratio (ν) using the relations: G = E/2(1+ν), Vs = √(G/p). The calculated values are compared with those reported by Jayawardena (2001) [4]. The difference in the parameters observed and calculated are insignificant.
Though Jayawardena (2001) [4] in his study has provided some mechanical properties of the most available rocks in Sri Lanka, there is hardly any evidence available on the exact location and depth of these rock types across the country. This lack of information has led to a major problem in developing the Numerical Model especially because the Finite Difference programme has been developed in a way where the material properties at a particular location should be known.

Table 1 -Summary of Mechanical Properties of Different Rock Types in Sri Lanka
Therefore, as an approximation in this study, the average value of mechanical properties of the mostly available rock types were used.

Boundary Conditions
There are three types of boundary conditions (quiet, free or fixed) which can be specified in FLAC software. After several trial runs on models with different boundary conditions, having quiet boundaries on all three sides on the model (as shown in Figure 5), was considered to be the most appropriate. Quiet boundaries allow wave dispersion through them without reflection. The top surface (bedrock profile) acts as a free boundary which allows wave reflection and surface wave generation.
A trial run with fixed boundary at the model base led to wave trapping totally masking the effect of attenuation, which was not supported by past observations from literature.

Selection of Input Ground Motion
The selected seismic scenario, which is based on the previous studies, comprises an earthquake of magnitude 6.9 having a 475 year return period with epicenter at 90 km distance from the west coast and a depth of 15 km below the ground surface (Seneviratne, 2016) [6]. There are no suitable seismic events recorded in and around Sri Lanka which fit the above requirement. Therefore, a set of earthquake records which suits the above conditions was selected from the PEER (Pacific Earthquake Engineering Research Center) database (PEER, 2006) for the analysis. The details of the selected earthquakes are given in Table 2.
The above procedure gives the input motion to the model at the bedrock surface. Due to dispersion, the input motion varies with depth below the bedrock surface. To select the most appropriate input motion variation, three input motion distributions, namely uniform, triangular and parabolic (Figure 6), were considered. In the cases of uniform and triangular ground motions, the profile development is based on PGA at the bedrock surface. The parabolic input motion profile was developed in such a way that it has PGA at the ground surface and two third (2/3) of the PGA at mid height. In both parabolic and triangular    Figure  7 illustrates the attenuation characteristics of Kalutara-Kalmunai profile due to seismic event simulating "Friuli" earthquake obtained for the three input motion profiles (uniform, triangular and parabolic).
In the case of uniform input motion profile, the attenuation characteristics show an unusual peak at about 40 km from Kalutara which is unusual considering the ground profile. In the case of triangular input motion, very low values of PGA were shown from Kalutara (close to input boundary) to a distance of 40 km inland which is not acceptable. In contrast, the attenuation characteristics shown by the case of parabolic motion profile was typical of what is expected. Therefore, the parabolic input motion was selected as the most suitable way of applying the input acceleration in this study.

Results and Discussion
The variation of PGA with distance for each selected east-west profile was obtained by averaging the response given by FLAC analysis for the seven earthquakes selected from PEER database ( Table 2). As an average, the time to run one model is around three hours with the post processing time of around one hour to obtain the final response spectrum. This makes four hours on average to complete one model.
Figures 8a-e illustrate the variation of PGA and ground elevation above mean sea level with distance along the profiles. Generally PGA decreases with distance. This is evident from the PGA distance relationships for Chilaw-Vakarai and Puttalam-Trincomalee sections throughout. A wider plateau seems to enhance the PGA as seen from the Galle-Arugambay and Kalutara-Kalmunai sections. This may be due to wave trapping. If the elevation decreases again the above effect seems to die out with the distance, at a rate lower than on the level ground.
The attenuation characteristics presented here are discussed in great detail in Dananjaya (2015) [8]. PGA at the bedrock level in 52 cities covered in this study under the assumed seismic scenario are given in Table 3.

Conclusions
When the results of the numerical simulation of ground response attenuation are studied, it can be noted that there is a significant amplification of the PGA where there is a highland for a broader width, but not amplified when there are sudden peaks. Since the width of the plateau seems to be the demarcating factor in the amplification of the PGA, this could be due to the fact that waves can be trapped in the plateau. When there are sudden peaks, the waves could be just passing by without entering inside the peak because of direct propagation of reflection. When the results are studied carefully it can be noted that the elevations beyond 500 m distributed for a distance beyond 100 km has amplified the PGA values. This is a crucial issue as most the attenuation relationships does not interpret this scenario and hence will be predicting lower PGA values for such locations.