Get 20M+ Full-Text Papers For Less Than $1.50/day. Subscribe now for You or Your Team.

Learn More →

Estimation of ground motion during the 18th September 2011 Sikkim Earthquake

Estimation of ground motion during the 18th September 2011 Sikkim Earthquake Geomatics, Natural Hazards and Risk Vol. 3, No. 1, February 2012, 9–34 Original article Estimation of ground motion during the 18th September 2011 Sikkim Earthquake S.T.G. RAGHUKANTH*, K. LAKSHMI KUMARI, and B. KAVITHA Department of Civil Engineering, IIT Madras, Chennai, India (Received 30 October 2011; in final form 30 November 2011) This article simulates the ground motion during the devastating 18th September 2011 Sikkim earthquake by analytical approaches. Surface level peak ground acceleration (PGA) values at several stations in the epicentral region have been estimated by the finite fault seismological model. A PGA contour map in the epicentral region is provided. It is found that very near to the epicentre, PGA would have reached more than 0.3g. The SPECFEM3D_GLOBE package which models the complete three-dimensional wave velocity structure and topography is used to simulate the displacement time histories on a global scale. A contour map of peak ground displacement (PGD) and permanent ground residual displace- ments (PGRD) in the near source region is provided. The estimated ground motions are validated with that obtained from damage values. 1. Introduction The state of Sikkim with Gangtok (27.338N, 88.68E) as the capital city lies in the Himalayan region (Figure 1). There has been an increase in population and infrastructure developments in the state of Sikkim in last few decades. As per census 2011, the total population of Sikkim is 0.6 million which is about 12.36% growth from the last census of 2001. The state of Sikkim comprising of four districts covers an area of about 7096 sq km with the population density registered at 86 persons per sq km. The Indian standard code of practice for earthquake resistant design of structures, IS 1893 (2002), has identified Sikkim state in zone of high seismic hazard (Zone IV) with the design peak ground acceleration as 0.24g. The 18th September 2011 event (M 6.9) which occurred on the Indo–Nepal border caused severe damages in the state of Sikkim. According to the India Meteorological Department (IMD) the epicentre of this quake was at 27.718N, 88.28E. This earthquake was widely felt in major cities, including Kolkata, Dhaka, Kathmandu and Chittagong. Several strong earthquakes of magnitudes greater than M 6 have occurred in Sikkim Himalaya in the past, but none of them have lead to a large scale damages compared to the 18th September 2011 event. This can be attributed to the increase in the growth of population and infrastructure development in this region over the last few decades. Although damages were reported in the entire Sikkim, towns and villages in north Sikkim district were worst affected due to proximity of the epicentre. This event triggered several landslides and rolling boulders in the *Corresponding author. Email: raghukanth@iitm.ac.in Geomatics, Natural Hazards and Risk ISSN 1947-5705 Print/ISSN 1947-5713 Online ª 2012 Taylor & Francis http://www.tandf.co.uk/journals http://dx.doi.org/10.1080/19475705.2011.646323 10 S.T.G. Raghukanth et al. Figure 1. Sikkim state. Red dots indicate location of landslides. Available in colour online. epicentral region. A total of 354 landslides were mapped in Teesta valley by National Remote Sensing Centre (http://www.nrsc.gov.in/) from satellite images after the earthquake. The strong shaking triggered landslides which caused damage to lifeline structures and significant building collapse. The official reports indicate that about 97 people were killed in Sikkim state. Casualties have been reported in neighbouring Nepal, Bhutan and Tibet. The financial losses due to this event are estimated to be around US$ 20 billion. The damage during the earthquake is mainly attributed to the strong near-source ground motion, surface ruptures and widespread slope failures in the epicentral region. For engineering purposes, the most sought-after information during an earthquake is the ground motion time history in the epicentral region. In this study, analytical and stochastic seismological approaches are used to estimate the ground motion during this event. These models are viable alternatives and popular among engineers and seismologists for computing the ground motion during earthquakes in regions lacking Strong Motion Accelerogram (SMA) data. These models have been successfully used to estimate ground motion during several past earthquakes (Singh et al. 1999, Iyengar and Raghukanth 2006, Raghukanth 2008, Raghukanth and Bhanuteja 2011). The regional velocity model, quality factor and stress drop are available for the Himalayan region (Parvez et al. 2002, Singh et al. 2004). With this limited information, one can use seismological approaches to estimate ground motion during this event. In this study, estimates of peak ground acceleration (PGA) values for the Sikkim event have been carried out by three different approaches. In the first approach, the finite fault seismological model of Boore (2009) is used for simulating acceleration time histories. Acceleration time histories are computed on grid size of 1 sq km covering a region of 200 km6 200 km in Sikkim region corresponding to hard rock type site condition. The site classification of Sikkim region based on topography as reported in www.earthquake.usgs.gov/hazards/apps/ vs30/ website combined with the correction factors available in IBC (2003) is used for estimating the surface level PGA. In the second approach, spectral finite element method (SPECFEM) is used to simulate ground motion. SPECFEM is based upon a Estimation of ground motion 11 weak formulation of the equations of motion and combines the flexibility of a finite- element method with the accuracy of a global pseudospectral method. The software package, SPECFEM3D_GLOBE which is a collection of Fortran subroutines for modelling the global seismic wave propagation is used to compute the displacement time histories in the epicentral region. Apart from these two analytical approaches, an estimate of PGA is obtained from the reported MMI values. 2. Seismotectonic setting and source of Sikkim earthquake The seismotectonic map of Sikkim Himalaya with superimposed seismicity is shown in Figure 2. It can be observed that two main Himalayan faults namely main boundary thrust (MBT) and main central thrust (MCT) are crossing the Sikkim state (Geological Survey of India (GSI) 2000). The MCT is not parallel to the MBT and arches in the form of a mushroom (Ray 2000). Apart from this, MCT also takes closed forms in Sikkim Himalaya. A large number of lineaments which extend for several kilometers lie transverse to the Himalayan faults in this region. The Tista Figure 2. Seismotectonic map of Sikkim Himalaya (GSI 2000). Available in colour online. 12 S.T.G. Raghukanth et al. lineament is thought to separate central and eastern Himalaya. The east Patna fault extends in to the Himalaya in the form of Arun lineament. The Kanchenjunga lineament which lies parallel to Arun lineament also extends from Indo-Gangetic plain to well inside the Himalayan region. The Tista and Purnia-Everest lineaments exhibit NW-SE trending where as Kanchenjunga, Arun and Everest lineament exhibit NE-SW trending. Detailed analysis of the seismicity in the Sikkim Himalaya indicates that movements along faults like the Gangtok lineament and the Tista lineament are responsible for many earthquakes in this region. De and Kayal (2003) analysed microearthquake data and observed that Sikkim Himalaya exhibit complex tectonic activity. The focal depth of the earthquakes in Sikkim Himalaya lies in between 0 and 45 km and occur by both thrust and strike-slip mechanism. Monsalve et al (2006) studied the seismicity in eastern Nepal and reported a bimodal distribution of focal depth in the Himalaya. They concluded that earthquakes in Himalaya occur both in the upper crust and around the crust-mantle boundary in the 50–100 km depth range. Recently, Hazarika et al. (2010) studied 356 local earthquakes in Sikkim and reported that seismic activity in this region can be attributed to transverse tectonics rather than underthrusting of India plate. The focal depth of the events in Sikkim Himalaya could be as high as 70 km and moment tensors indicate predominantly strike slip faulting. They attributed the strike slip mechanism in Himalayan thrust environment to an oblique convergence of Indian- Asian collision. Sikkim has experienced several damaging earthquakes in the past. The MMI in Sikkim Himalaya has been reported to be VIII for the great Bihar-Nepal earthquake (Mw 8.3) of 1934 (GSI 2000). The 1988 event (Mw 6.5) which occurred in the epicentral region of 1934 great earthquake caused landslides and subsidence in Sikkim (De and Kayal 2003). Damages were observed in East Sikkim for the recent 14th February 2006 (Mw 5.3) moderate event. The focal mechanism solution and field investigations indicate that the Tista lineament caused 18th September 2011 event. The USGS (http://neic.usgs.gov/) website reported the strike, dip and rake angle of the rupture plane as 1308,908 and 1688, respectively. The focal depth of this event is obtained as 35 km. The Global Centroid Moment Tensor database (http:// www.globalcmt.org/) obtained the focal depth as 47 km with strike, dip and rake angle as 1338,738 and 1638. The focal mechanism solutions reported by the above two agencies indicate nearly strike slip faulting. On the other hand, the India Meteorological Department (IMD) estimated the focal depth of this event as 10 km. 3. Simulation of SMA using stochastic finite-fault simulation method Due to unavailability of strong motion records for the Sikkim earthquake, ground motions are simulated by an analytical model. The stochastic finite fault approach of Boore (2009) which is an improved version of the methodology proposed by Motazedian and Atkinson (2005) is used to simulate the rock level ground motion for the seven geological regions in India. The theory and application of seismological models for estimating ground motion has been discussed in detail by Boore (1983, 2009) and Motazedian and Atkinson (2005). The application of seismological model for estimating ground motion in northeastern India has been discussed in detail by Raghukanth and Somala (2009) and Raghukanth and Dash (2009). A brief description of the method starting with the frequency domain representation of the ground acceleration will bring out the advantages of this approach. In this method, Estimation of ground motion 13 the rectangular fault plane is divided into N number of subfaults and each subfault is represented as a point source. The Fourier amplitude spectrum of ground motion [A(r,f)] due to the jth subfault at a site is derived from the point source seismological model, expressed as: A ðr; fÞ¼ CH S ðfÞFðfÞDðr; fÞPðfÞ: ð1Þ j j j Here, C is a scaling factor, S (f) is the source spectral function, D(r,f) is the diminution function characterizing the quality of the region, P(f) is a filter function, F(f) is the site dependent function that modifies the bed rock motion in the vertical direction and H is a scaling factor used for conserving the energy of high-frequency spectral level of subfaults. In this study, following Brune (1970), the principal source model S (f) is taken as: 2 0j hi S ðfÞ¼ð2pf Þ : ð2Þ 1 þ f f 0j Here, f is the corner frequency and M is the seismic moment of the jth subfault. 0j 0j The three important seismic source parameters M , f and the stress drop (Ds) are 0j 0j related through: 1=3 Ds 1=3 6 1=3 f ¼ 4:9  10 ðN Þ N V : ð3Þ 0j s Rj 0j Here, V stands for the shear wave velocity in the source region, corresponding to bedrock conditions. N is the cumulative number of ruptured subfaults by the time Rj rupture reaches the jth subfault. The spatial modifying function D(r, f) is given by: pfr Dðr; fÞ¼ Gexp ; ð4Þ V QðfÞ where G is the geometric attenuation factor. The other term denotes anelastic attenuation with hypocentral distance r and the quality factor as Q. The spatial spread of the ground motion depends sensitively on the quality factor of the local region. The constant C of Equation (1) is: pffiffiffi R 2 yf C ¼ ð5Þ 4prV where hR i is the radiation coefficient averaged over an appropriate range of yf azimuths and take-off angles and r is the material density at the focal depth. The coefficient 2 in the above equation arises as the product of the free surface amplification and partitioning of energy in orthogonal directions. The scaling factor for jth sub-fault, H based on the squared acceleration spectrum is taken as (Boore 2009) pffiffiffiffi H ¼ N ; ð6Þ 0j 14 S.T.G. Raghukanth et al. where f is the corner frequency at the end of the rupture, which can be obtained by substituting N (t) ¼ N in Equation (3). The filter function P (f) is taken as: R j pffiffiffiffi 1 þðf f Þ 0j P ðfÞ¼ : ð7Þ j 1 þðf=f Þ The moment of jth subfault is computed from the slip distribution as follows: M D 0 j M ¼ ð8Þ 0j j¼1 Here, D is the average final slip acting on the jth subfault. M is the total seismic j 0 moment on the fault. To further account for earthquake rupture, Motazedian and Atkinson (2005) introduced the concept of pulsing area where the cumulative number of active subfaults, N increases with time at the initiation of rupture and Rj becomes constant at some fixed percentage of the total rupture area. This parameter determines the number of active subfaults during the rupture of jth subfault. These many subfaults are used in computing the corner frequency in Equation 3. The above is a general finite source model expressed in the frequency domain, valid for any region if only the various controlling parameters can be selected suitably. The effects that are specific to Sikkim earthquake appear in the quality factor Q(f), stress drop value (Ds), focal depth, dip angle of the fault plane and the site amplification function F(f). Thus, these five parameters have to be selected to simulate ground motion for this event. Here lies strength of this approach since almost all required parameters for Himalayan region are available. The spatial modifier of Equation 4 consists of a general geometrical attenuation term G. This is taken following Singh et al. (1999) as: ; for r  100 G ¼ ð9Þ pffiffi ; for r > 100 10 r The shear wave velocity and density at the focal depth are fixed at 3.5 km/s and 2900 kg/m , respectively corresponding to bedrock (Singh et al. 1999, Hazarakia et al. 2010). In the past, several seismologists have analysed instrumental data to arrive at the value of stress drop that can happen in Himalayan region for a given magnitude of earthquake. Singh et al. (2004) have studied the 1991 Uttarkashi earthquake and 1999 Chamoli earthquake and derived the frequency dependent 0.8 quality factor for Himalayan region as 253f . The other important regional parameter is the stress drop which can vary from 50 to 200 bars in Himalayan region (Singh et al. 2004). Since there are discrepancies regarding the focal depth for this 18th September event, it is varied from 25 km to 50 km. The regional velocity model used is the one reported by Parvez et al. (2003). Details of the regional model are presented in Table 1. Modification between bedrock and A-type site is a linear problem in one-dimension and hence for such sites, amplification can be directly found by the quarter-wavelength method of Boore and Joyner (1997). The estimated amplification function is presented in Figure 3. Apart from the above parameters, the S-wave radiation coefficient (hR i) varies randomly within particular intervals. yf Here, following Boore and Boatwright (1984), the S-wave radiation coefficient is taken to be in the interval 0.48–0.64. In addition to the above parameters the slip Estimation of ground motion 15 Table 1. Regional Velocity Model (Parvez et al 2003). Depth (km) V (km/s) V (km/s) Density (kg/m ) Q Q p s p S 0–4.9 4.5 2.5 1800 3600 1600 4.9–19 5.9 3.1 2400 2200 1100 19–32 6.1 3.5 2600 1400 500 32–50 7.1 4.1 2900 800 200 4 50 7.6 4.2 3300 800 200 Figure 3. Amplification function for A-type sites in Sikkim region. distribution and pulsing percentage area is also required in the simulation. The pulsing percentage area is varied from 25% to 75% (Atkinson and Boore 2006). 3.1. Slip model The stochastic finite fault approach requires slip distribution on the rupture plane which is highly complex. Mai and Beroza (2002) have proposed stochastic characterization of slip, in which slip distribution is modelled as a anisotropic Von Karman random field. They modelled slip distributions of several past earthquakes and developed empirical equations for estimating correlation lengths along the two dimensions of the fault plane from the magnitude of the earthquake. In this approach, once the dimensions and orientation of rupture plane are fixed, the slip distribution on the fault plane is generated as a spatial random field with correlation lengths depending on the magnitude of the earthquake. Accordingly, in the present analysis, the correlation lengths for these four events in strike and dip direction are determined from the scaling relations of Mai and Beroza (2002). From the derived correlation lengths, slip is simulated as a random field by the spectral representation method of Shinozuka and Deodatis (1996). The hypocentre in each case is taken nearer to the regions of maximum slip release (Mai et al. 2005). This sample slip field on the fault plane is discretized into subfaults of size 1 km6 1 km and seismic moment of each point source is computed using Equation 8. Since the moment 16 S.T.G. Raghukanth et al. distribution on the rupture plane is modelled as random field, one can generate an ensemble of time histories at the surface by numerical simulation. Based on the scaling relations of Mai and Beroza (2002), the dimensions of the rupture plane are estimated from magnitude as 35 km6 20 km. In Figure 4, two samples of the slip field are shown. To cover all possible peculiarities in the slip distribution, fifty samples of slip field are generated for each event. 3.2. Sample ground motion From the above discussion it is seen that once the stress drop, slip field, site amplification and Q factors are known, the Fourier amplitude spectrum of ground acceleration for any combination of magnitude (M ) and hypocentral distance (r) can be expressed within the limitations of a finite source seismological source mechanism model. This is a random process and hence in the time domain this represents an ensemble of accelerograms. These samples can be retrieved from Equation 1 in three steps. First, a Gaussian stationary white noise sample of length equal to the strong motion duration (Boore and Atkinson 1987) T ¼ðÞ 1=f þ 0:05 r ; ð10Þ Figure 4. Sample slip realizations for the Sikkim earthquake (M 6.9). Available in colour online. Estimation of ground motion 17 is simulated for each subfault. This sample is multiplied by a suitable non-stationary modulating function suggested by Saragoni and Hart (1974). The simulated sample is Fourier transformed into the frequency domain. This Fourier spectrum is normalized by its root mean-square value and multiplied by the seismological source-path-site function A(f) of Equation 1. The resulting function is transformed back into the time domain to get a sample of the ground motion accelerogram for each subfault. The simulated acceleration time histories for all the subfaults are summed up with time delay of Dt to account for the rupture velocity to obtain the final ground acceleration, aðtÞ¼ a t þ Dt : ð11Þ j j j¼1 Any number of accelerograms can be simulated in this fashion to compute samples of response spectra corresponding to 5% damping. Since the stress drop, focal depth, dip, amplification function, radiation coefficient, pulsing percentage area are random variables and slip as a random field, the uncertainty arising out of these parameters is included in the simulation. The input parameters used in the seismological model for the Sikkim event are reported in Table 2. Accordingly, 50 samples of these seven input parameters are generated from uniform distribution and these are combined using the Latin Hypercube sampling technique (Iman and Conover 1980) to simulate 50 sets of random seismic parameters for the Sikkim event. By the above procedure, acceleration time histories are computed on grid size of 1 sq km covering a region of 200 km6 200 km in Sikkim region. In Figure 5(a,b) mean and standard deviation of rock level peak ground acceleration are shown. It can be observed that the contours are axisymmetric and the PGA is of the order of 0.22g+ 0.06g near the epicentre. The standard deviation of PGA decreases with increase in epicentral distance. Table 2. Input Parameters in the stochastic finite fault seismological model. Parameter Sikkim earthquake Mag., M 6.9 Location Latitude 27.71 Longitude 88.20 Depth to top of the fault (km) 25–50 Stress drop (bars) 50–200 0.8 Q( f ) 253f Site effects Quarter-wavelength method Fault dimensions (km) 356 19 Fault orientation Strike(8) 120–140 Dip (8) 70–90 Pulsing area percentage 25–60 Rupture velocity 0.8b Slip distribution Random field 18 S.T.G. Raghukanth et al. Figure 5. (a) Contour map of mean PGA (g) at A-type rock. Black dots indicate location of landslides; (b) Standard deviation of PGA (g) at A-type rock. Available in colour online. Estimation of ground motion 19 3.3. Surface level PGA The local site property is expressed in terms of the average shear wave velocity in the top 30 m of the soil. This is the recognized IBC (2003) approach, wherein sites are classified as A: (V 4 1.5 km/s); B: (0.76 km/s5V  1.5 km/s); C: (0.36 km/s 5 30 30 V  0.76 km/s); D: (0.18 km/s 5 V  0.36 km/s). E- and F-type sites with V 30 30 30 0.18 km/s are susceptible for liquefaction and failure. The above simulated PGA values are valid for A-type rock site condition. The stochastic finite fault simulation described above can be extended to the surface with the help of soil profiles and V values sampled from the region. Thus, for a specific site, precise correction of the bedrock results will be possible only when the soil section data are available, including the variation of V with depth. Unfortunately, such information on lateral variation of shear wave velocity profiles in Sikkim is not available. Wald and Allen (2007) developed empirical equations relating slope of the topography variations with shear wave velocity from global database. In this study, the approach proposed by Wald and Allen (2007) is used for classifying the local site condition in the epicentral region. The variation of shear wave velocity in the top 30 m (V ) obtained from global V map server (http://earthquake. usgs.gov/hazards/apps/vs30/) is shown in Figure 6. It can be observed that most of the sites correspond to C-type site Figure 6. Shear wave velocity in top 30 m (V , m/s), (http://earthquake.usgs.gov/hazards/ apps/vs30/). Available in colour online. 20 S.T.G. Raghukanth et al. condition. Once V is known, the amplification factors given in IBC (2003) listed in Tables 3 and 4 are used to estimate the surface level PGA and 1 s spectral acceleration. In Figure 7 the resulting surface level PGA contour is shown along with the location of landslides. It can be observed that the contours are axisymmetric which can be attributed to the focal depth of the earthquake. The effect of the local site condition on PGA can also be observed in the figure. The PGA in the near source region is as high as 0.35g. At Lachen, Ravangla, Teesta, Chungthang and Mangan, PGA lies in between 0.25g and 0.29g. PGA estimates obtained from the stochastic finite fault simulation at important villages and towns which suffered damage are shown in Table 5. The contours of spectral acceleration (S )at 1 s natural time period are shown in Figure 8. The spectral acceleration in the epicentral region is as high as 0.2g. At places like Chungthang and Mangan which suffered heavy damage it is of the order of 0.14g. 4. PGA from MMI data The PGA values estimated above are based on an analytical model. The above results are now correlated with field observation of damage after the event. Although strong motion data are not available, information on damages at several villages and towns in epicentral region are available from filed investigations and media reports. The MMI values are fixed up based on these damages. These intensities are reported in Table 5. It can be observed that maximum intensity of VIII has been felt at Chungthang.To validate the estimates obtained from the seismological model, the MMI values have been converted into PGA values using an empirical relation. Previously, Iyengar and Raghukanth (2003) compared Indian MMI data with instrumental PGA data with the help of 43 sample values and obtained a linear relation between ln(PGA) and MMI as: lnðÞ PGA=g¼ 0:6782 MMI  6:8163; s lnðÞ e ¼ 0:7311: ð12Þ Table 3. Short period correction factor for various site conditions (IBC 2003). Site category 0.1g 0.2g 0.3g 0.4g 0.5g A 0.8 0.8 0.8 0.8 0.8 B 1.0 1.0 1.0 1.0 1.0 C 1.2 1.2 1.1 1.0 1.0 D 1.6 1.4 1.2 1.1 1.0 Table 4. Long period correction factor for various site conditions (IBC 2003). Site category 0.1g 0.2g 0.3g 0.4g 0.5g A 0.8 0.8 0.8 0.8 0.8 B 1.0 1.0 1.0 1.0 1.0 C 1.7 1.6 1.5 1.4 1.3 D 2.4 2.0 1.8 1.6 1.5 Estimation of ground motion 21 Figure 7. Surface level peak ground acceleration (g). Black dots indicate location of landslides. Available in colour online. The above equation is for the expected PGA values only. In practice, PGA will be distributed as a random variable and hence recorded values will show considerable variation in a region with nearly constant MMI. In Table 5, the MMI-PGA values are compared to the corresponding PGA values estimated from stochastic approaches. It can be observed that the comparison between the two estimates is consistent. 5. Displacement time history The above one-dimensional stochastic finite fault seismological model is efficient in simulating accelerations including high frequencies. However, this model do not incorporate seismic phases that partition the energy into high frequency body waves and low frequency surface waves and hence is not valid in simulating displacement time histories, which are controlled by low frequencies. Moreover, the displacements are sensitive to fault rupture process. The three-component displacement time histories including the near-source features like directivity and fling are required to understand the evolution of nonlinear structural deformations in large bridges, dams and other important structures. Since many hydro-power projects are under construction in Sikkim state, it would be of interest to engineers to know the surface displacements caused by this moderate event. Analytical techniques based on kinematic models and elastic wave propagation approaches are very efficient in 22 S.T.G. Raghukanth et al. Table 5. Comparison among various PGA estimates. Latitude Longitude MMI–PGA PGA (g) Station (8N) (8E) MMI Reference (g) Stoch. fault Lachen 27.73 88.55 VII India today 0.13 0.297 Lachung 27.7 88.75 VII India today 0.13 0.226 Chungthang 27.61 88.63 VIII The hindu; earthquake-report.com 0.25 0.268 Mangan 27.52 88.54 VIII news.oneindia.in 0.25 0.281 Ravangla 27.30 88.36 VI asc-india.org 0.06 0.248 Teesta 27.03 88.26 VII The Hindu 0.13 0.150 Dikchu 27.38 88.52 VII asc-india.org 0.13 0.248 Rangpo 27.18 88.53 VI The Economic Times 0.06 0.181 Taplejung 27.35 87.66 V earthquake-report.com 0.03 0.180 Yuksom 27.37 88.22 VI asc-india.org 0.06 0.279 Singtam 27.23 88.50 VII news.oneindia.in 0.13 0.223 Gangtok 27.20 88.40 VII www.reuters.com 0.13 0.199 Darjeeling 27.05 88.26 VII Hindustan Times 0.13 0.156 Kalimpong 27.07 88.48 VII NDTV 0.13 0.156 jorethang 27.13 88.28 VII haalkhabar.in 0.13 0.186 Pentong 27.51 88.53 VI The Times of India 0.06 0.278 Yadong 27.51 88.97 V earthquake-report.com 0.03 0.150 Sankhuwasabha 27.58 87.33 VII Nepal Mountain News 0.13 0.136 Sakyong 27.33 88.60 VII The Times of India 0.13 0.209 Nayabazaar 27.13 88.26 VI haalkhabar.in 0.06 0.185 Pegong 27.59 88.65 VIII The Economic Times 0.25 0.255 Lingzya 27.55 88.45 VIII The Times of India 0.25 0.301 Estimation of ground motion 23 Figure 8. Surface level Spectral acceleration (g) Contours at T ¼ 1 s. Available in colour online. explaining source directivity effects as well as effects of underground medium on ground motion. The Sikkim event also triggered extensive landslides in the epicentral region which is mainly due to varying Himalayan topography (Figure 1). This demands simulation of ground motion by including the topography of the region. In this study, spectral finite element method (SPECFEM) is used to simulate ground motion for the Sikkim earthquake. SPECFEM is based upon a weak formulation of the equations of motion and combines the flexibility of a finite-element method with the accuracy of a global pseudospectral method. The software package, SPEC- FEM3D_GLOBE which is a collection of Fortran subroutines for modelling the global seismic wave propagation is used to compute the displacement time histories in the epicentral region. The effect of lateral variations in material properties, topography and bathymetry, the oceans, have been incorporated in simulating ground motion. 5.1. Spectral finite element method In this study, the spectral finite element method introduced by Patera (1984) and developed by Komatitsch and Tromp (2002a,b) for global seismic wave propagation is used to model the Earth for simulating the ground motion during the Sikkim event. The Earth’s interior can be broadly divided into four layers namely crust, mantle, outer core and inner core. The crust, mantle and inner core are solid where as the outer core is composed of liquid material. 24 S.T.G. Raghukanth et al. 5.1.1. Crust and mantle. The governing equation for the Earth medium in spherical domain (r,y,f) subjected to a unit impulsive double couple applied at location (r )in the crust is given by (Dahlen and Tromp 1998) @ uðr; tÞ @uðr; tÞ rðrÞ þ 2o  ¼ H uðr; tÞþ fðr; tÞð13Þ @t @t where u(r,t) is the displacement field, r(r) is the material density, o is the Earth’s angular rotation vector and f(r,t) denotes the body force. The elasto-gravity operator, H is given by: H uðr; tÞ¼r  sðuÞþrðru  gÞ  rðruÞg: ð14Þ In the above equation, g is the gravitational acceleration. The body force components for a unit impulsive double couple are: fðr; tÞ¼M rdðr  r ÞSðtÞð15Þ where d denotes the Dirac delta function, S(t) denotes the source time function and M contains the components of the symmetric seismic moment tensor related to the orientation of the fault plane as (Aki and Richards 1980) M ¼M sin ðdÞcos ðlÞsin ð2fÞþ sin ð2dÞsin ðlÞsin ðfÞ yy 0 M ¼MðÞ sin ðdÞcos ðlÞcos ð2fÞþ 0:5sin ð2dÞsin ðlÞsin ð2fÞ yj 0 M ¼MðÞ cos ðdÞcos ðlÞcos ðfÞþ cos ð2dÞsin ðlÞsin ðfÞ ry 0 ð16Þ M ¼ M sin ðdÞcos ðlÞsin ð2fÞ sin ð2dÞsin ðlÞcos ðfÞ jj 0 M ¼ MðÞ cos ðdÞcos ðlÞsin ðfÞ cos ð2dÞsin ðlÞcos ðfÞ rj 0 M ¼ MðÞ sin ð2dÞsin ðlÞ rr 0 where f, d and l are the strike, dip and rake angles of the fault. The seismic moment (M ) is estimated from the average slip (D), fault area (A) and rigidity modulus (m) as: M ¼ mAD ð17Þ 5.1.2. Outer core. Assuming the fluid to be stably stratified and isentropic, the equation of motion in the fluid outer core can be written as: @ uðr; tÞ @uðr; tÞ 1 rðrÞ þ 2o  ¼r kr u þ u  g ð18Þ @t @t r where k is the fluid’s bulk modulus. 5.1.3. Inner core. The equation of motion for inner core, which is composed of solid material similar to crust and mantle region can be written as: @ uðr; tÞ @uðr; tÞ rðrÞ þ 2o  ¼ H uðr; tÞð19Þ @t @t Estimation of ground motion 25 5.1.4. Boundary conditions. The vanishing of stresses at the free surface on continental crust can be expressed as: sðÞ r; t nrðÞ ¼ 0 8r 2  ; ð20Þ where  is the free surface of the continental crust and n being the unit normal to the interface. To model the oceans and large lakes on Earth’s surface, the stress free boundary condition at the ocean crust boundary (OCB) is modified to include the fluid pressure. Neglecting the variations of gravity and angular rotation in the oceans, if h is the height of the water column, the fluid pressure at the ocean crust boundary can be written as: @ uðr; tÞ sðr; tÞ nðrÞ¼ r h nðrÞ 8r 2  ; ð21Þ OCB @t with r being the density of the sea water. At crust-mantle boundary (CMB) and inner core boundary (ICB) normal component of velocity and normal stress have to be continuous. This can be expressed as: @u nrðÞ ¼ 0 8r 2 ICB; CMB @t ð22Þ ½ jsðÞ r; t  nðrÞj¼ 0 8r 2 ICB; CMB The equations of motion in crust and mantle, outer core and inner core (Equations 13, 18 and 19) has to be solved subjected to the above boundary conditions (Equations 21–22). The spectral finite element method is based on the integral or weak formulation of the problem. Earth is divided into a number of hexahedral non- overlapping elements. In SPECFEM high-degree Lagrange interpolant is used to represent functions on the element. In the hexahedral volume element, all the functions are interpolated by triple products of Lagrange polynomials of degree n at the classical n þ1 nodes called as Gauss-Lobatto-Legendre (GLL) quadrature points. The numerical integration of these functions over volume elements are approximated using the GLL integration rule. The contributions from all the elements that share a common global grid point are summed before the time marching. The method is implemented in parallel computing using the message passing technique. (Komatitsch and Tromp 2001). The SPECFEM3D Globe package (Bozdag et al. 2010) which is a collection of Fortran subroutines is used to simulate three-dimensional global wave propagation for the Gujarat earthquake. The SPECFEM3D package has been implemented on the VEGA super cluster available at IIT Madras using job scheduling software ‘portable batch system’ (Raghukanth and Bhanuteja 2011). The VEGA super cluster is based on HP Proliant DL160 G5 servers with Quad-core Intel Xeon processors. This Cluster consists of 128 nodes. Each node consists of eight processors. The first step in SPECFEM is to generate the spectral-element mesh. This is shown in Figure 9. The Earth is divided into six chunks based on the cubed sphere mapping. In this study, to allow for a denser mesh, among the six chunks that constitute the globe, single chunk mesh is used for estimating the ground motion. The chunk is centered at 308N and 808E. The size of the chunk is 6086 608 with absorbing boundary conditions at the artificial edges of the mesh (Figure 9). This chunk is divided into 256 slices with each slice 26 S.T.G. Raghukanth et al. Figure 9. Spectral element mesh of the globe. Red line indicates the size of the chunk used in this study. Available in colour online. assigned to one processor. The total number of processors needed to handle the mesh would be 256. Each slice is further subdivided into 406 40 spectral elements at the earth’s surface. The total number of spectral elements in the one chunk mesh is 13.2 million. Each spectral element consists of 125 grid points. Since grid points are shared among the elements, the entire globe is represented by 872.2 million grid points. The total number of degrees of freedom in entire mesh is 2.31 billion. The average grid spacing at the Earth’s surface is about 2.61 km. The simulated displacement time histories are valid up to a shortest period of 4.5 s. Once all the three regions in Earth’s interior have been discretized, the global system of equations to be solved can be written as: M€u þ Wu_ þ Ku þ Bu ¼ F ð23Þ where M and K are the global mass matrix and stiffness matrix. W contains terms related to angular rotation vector and B is related to the boundary interactions at Estimation of ground motion 27 CMB and ICB and F is the source term. Explicit expressions for the mass, stiffness, W and B matrices at elemental level and further construction of these matrices at global level are available in the article of Komatitsch and Tromp (2002a,b). A explicit second-order finite difference scheme is used to march the above equation in time. The finite source rupture model reported in Figure 4 is used for simulating ground motion. The 3D velocity model for the mantle and crust are taken from the work of Kustowski et al. (2008) and Bassin et al. (2000). The topography and bathymetry model available at 5-min interval grid (ETOPO5) is from the US National Oceanic and Atmospheric Administration (http://www.ngdc.noaa.gov/). The 35 km6 19 km fault plane is divided into 665 subfaults and slip acting on each subfault is represented as a point double-couple varying in time as a smooth ramp function. The moment tensor at each subfault is computed from slip and orientation of the fault plane (Equations 15 and 16). The rise time in the ramp function for each subfault is estimated from the empirical relation of Somerville et al (1999) as: 1=3 T ¼ 2:03  10 M ð24Þ The triggering time of each point source is the estimated from the ratio of the distance between the point source and the hypocentre to rupture velocity. The three components of displacements on the surface of the Earth, for each double couple source are obtained by solving the Equation 23. The ground displacements due to the complete fault rupture, is obtained by summing-up the contribution of all the point double-couples. A typical ground motion simulation including the mesher and solver requires about 19 h using 256 processors on an HP Proliant DL160 G5 servers operating at 3 Ghz speed. In this study, displacement time histories are simulated for two sample realization of slip field. The focal depth of the two samples are 35 km and 26 km. In Figure 10(a,b), the simulated displacement time histories at several other stations which suffered heavy damage in the Sikkim for the both the slip models shown in Figure 4. Permanent displacements can be observed in the simulated time histories (Figure 10) at several stations in the near-field region. The displacement field has been computed at a spacing of 2.6 km covering a region of dimension 350 km6 350 km around the fault. The contours of peak ground displacement in horizontal and vertical direction corresponding to slip model shown in Figure 4(a) is shown in Figure 11(a,b) along with location of landslides. It can be observed that PGD values are high at Mangan, Chungthang, Dikchu and Gangtok. The comparison between the PGD contours and location of landslides is good. It is well known that, at low frequencies, near-field ground motions will be strongly influenced by the dimensions of the fault plane, slip distribution on the rupture plane, direction of rupture propagation and orientation of the station with respect to the fault. The ground motions also follow certain radiation patterns and exhibit rupture-directivity effects. These low frequency characteristics can be observed from the above figures. Figure 12(a–c) presents the contours of permanent surface displacements in the north, east and vertical directions corresponding to the slip model shown in Figure 4(a). The near-source features such as permanent ground displacement can be observed at sites near the fault. These static displacements decrease with an increase in epicentral distance and are valid at bedrock level. Since the top soil would behave in a nonlinear fashion, the large displacements in the linear model would hint at ground failures due to the high level of strains. The field 28 S.T.G. Raghukanth et al. Figure 10. Analytical surface displacement time histories (Cm.) o-PGD, Blue line-slip model of Figure 4(a). Red line simulated from the slip model shown in Figure 4(b). Available in colour online. Estimation of ground motion 29 Figure 11. (a) Peak ground horizontal displacement (cm) in the for the slip model shown in Figure 4(a) (Focal depth ¼ 35 km). Black dots indicate location of landslides. (b) Peak ground displacement (cm) in vertical direction the for the slip model shown in Figure 4(a) (Focal depth ¼ 35 km). Black dots indicate location of landslides. Available in colour online. 30 S.T.G. Raghukanth et al. Figure 12. (a) Permanent ground displacement (cm) simulated from slip model shown in Figure 4(a) (Focal depth ¼ 35 km). East component. Black dots indicate location of landslides. (b) Permanent displacement (cm) simulated from slip model shown in Figure 4(a) (Focal depth ¼ 35 km). North component. Black dots indicate location of landslides. (c) Permanent ground displacement (cm) simulated from slip model shown in Figure 4(a) (Focal depth ¼ 35 km). Vertical component. Black dots indicate location of landslides. Available in colour online. Estimation of ground motion 31 investigations carried out after the earthquake indicated ground cracks and openings on the order of a few centimetres in the near field region which matches with the permanent displacements simulated from SPECFEM (Figure 12). 5. Summary and conclusions This article estimates the ground motion during the 18th September 2011 Sikkim earthquake (M 6.8) by analytical approaches. In absence of strong motion data, it is more scientific to use the stochastic finite fault seismological model (Boore 2009) for generating large samples of artificial ground motion time histories. With the help of such a sample where in the uncertain source, path and site parameters can be randomly varied, one can arrive at reliable ground motion relations for a given earthquake. Stochastic seismological model for simulation of ground motion depends on region specific earthquake model parameters such as stress drop, slip distribution, dip of the fault plane, focal depth, site amplification function and the quality factor. These model parameters are available for Himalayan region. This helps in fixation of the limit of uncertainties in the input parameters for Sikkim earthquake (Table 1). Accelerations at A-type rock level have been simulated in a 200km6 200km region around the epicentre. The approach of Wald and Allen (2007) has been used to obtain the V values in the epicentral region. The soil has been classified into B- type, C-type and D-type sites depending on the V and correction factors given in IBC (2003) are used to arrive at the surface level PGA. The contour map of mean PGA is shown in Figure 7. To validate the simulated ground motions, further estimates of PGA have been obtained from the MMI values. The comparisons at specific stations are shown in Table 5. It is seen that the PGA values obtained from both the approaches are consistent. The differences at some stations are attributed to local site effects and the complex topography of the region. The spectral finite element method is used to simulate the ground motion in the entire Earth medium including the Himalayan topography. The computational difficulties put limits on the high frequency end of the spectrum and hence numerical results are limited to displacement time histories only. Since the required computation for SPECFEM is large, the ground motion time histories have been obtained for two sample realizations of the slip distribution (Figure 4). Contours of a PGD and PGRD over a region of 350 km6 350 km, in and around the epicentre have been obtained for one sample realization of the slip distribution. The simulated time histories contains permanent ground residual displacement indicating ground failures in the near source region which matches with reported in the field investigations. The simulated ground motions in this study bring out the directivity effects due to fault orientation and rupture direction. Figure 12 reflects the permanent movement at the end of the earthquake shaking due to the effect of slip across the fault. These static deformations are observed in the near neighbourhood of faults and decrease sharply with distance as one moves away from the faults. The displacements will be of use in design of infrastructural facilities such as tunnels, pipelines, transmission lines that run across faults in Sikkim. The damages to structures during the Sikkim earthquake have been mainly attri- buted to poor construction practices in this region which increases the vulnerability of this region. In order to reduce the vulnerability, one requires ground motion of damaging earthquakes. The results obtained from this study can be of varied use in 32 S.T.G. Raghukanth et al. infrastructural developments in Sikkim Himalaya. For example, the obtained PGA and displacement time histories will be of use in the aseismic design of new structures and also for assessing the performance of existing structures. The contour maps of PGA, PGD and PGRD presented in Figures 7, 11 and 12 will be useful to engineers in selecting the amplitude of motions to be considered in design to achieve the objective of no collapse and preservation of life safety for events of comparable magnitude in this region. The estimated ground motion contours obtained from this study can be directly used to understand the pattern of damage and also the behaviour of structures which did not suffer appreciable damage during the main event. Since such moderate earthquakes will inevitably occur in the future, the results obtained in this study can be an important aid for emergency planning. Since near-field strong motion records are sparse in the Himalayan region, the estimated time histories in this study will be useful in understanding ground motion during a future earthquake of comparable magnitude in the Himalayan region. Acknowledgements Thanks to Prof. C.V.R. Murthy for his encouragement and help in fixing the MMI values at several towns in Sikkim state. References AKI, K. and RICHARDS, P.G., 1980, Quantitative Seismology: Theory and Methods, Vol. I. pp. 1–700. (New York: WH Freeman). ATKINSON, G.M. and BOORE, D.M., 2006, Earthquake ground-motion prediction equations for eastern North America. Bulletin of the Seismological Society of America, 96, pp. 2181– BASSIN, C., LASKE, G. and MASTERS, G., 2000, The current limits of resolution for surface wave tomography in North America. EOS Trans AGU, 81, pp. F897. BOORE, D.M., 1983, Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra. Bulletin of the Seismological Society of America, 73, pp. 1865–1894. BOORE, D.M., 2009, Comparing stochastic point-source and finite-source ground-motion simulations: SMSIM and EXSIM. Bulletin of the Seismological Society of America, 99, pp. 3202–3216. BOORE, D.M. and ATKINSON, G.M., 1987, Stochastic prediction of ground motion in eastern North America. Bulletin of the Seismological Society of America, 4, pp. 460–477. BOORE, D.M. and BOATWRIGHT, J., 1984, Average body-wave radiation coefficients. Bulletin of the Seismological Society of America, 74, pp. 1615–1621. BOORE, D.M. and JOYNER, W.B., 1997, Site amplifications for generic rock sites. Bulletin of the Seismological Society of America, 87, pp. 327–341. BOZDAG, E., CHEN, M., HJORLEIFSDOTTIR, V., KOMATITSCH, D., LIU, Q., MESCHEDE,M., MICHEA, D., PETER, D., SAVAGE, B., SCHUBERTH, B., SIEMINSKI, A., TAPE, C., TROMP,J. and ZHU, H., 2010, SPECFEM3D_GLOBE User Manual. Computational infra- structure for geodynamics, pp. 1–67. (France: Princeton University, California Institute of Technology and University of Pau/CNRS/INRIA). BRUNE, J.N., 1970, Tectonic stress and the spectra of seismic shear waves from earthquakes. Journal of Geophysical Research, 75, pp. 4997–5009. DAHLEN, F.A. and TROMP, J., 1998, Theoretical Global Seismology, pp. 1–1027. (Princeton, NJ: Princeton University Press). Estimation of ground motion 33 DE, R. and DE KAYAL, J.R., 2003, Seismotectonic model of the Sikkim Himalaya: constraint from microearthquake surveys. Bulletin of the Seismological Society of America, 933, pp. 1395–1400. Geological Survey of India (GSI), 2000, Seismotectonic Atlas of India and Its Environs, pp. 1– 87. (Kolkata: Geological Survey of India). HAZARIKA, P., RAVI KUMAR, M., SRIJAYANTHI, G., SOLOMON RAJU, P., PURNACHANDRA RAO,N. and SRINAGESH, D., 2010, Transverse tectonics in the Sikkim Himalaya: evidence from seismicity and focal-mechanism Data. Bulletin of the Seismological Society of America, 100, pp. 1816–1822. IBC, 2003, International Building Code, pp. 1–673. (Washington, DC: International Code Council, ICC). IMAN, R.L. and CONOVER, W.J., 1980, Small sample sensitivity analysis techniques for computer models, with an application to risk assessment. Communications in Statistics, A9, pp. 1749–1842. IS 1893, 2002, The Indian Standard Code, Criteria for Earthquake Resistant Design of Structures, pp. 1–40. (New Delhi: Bureau of Indian Standards). IYENGAR, R.N. and RAGHUKANTH S.T.G., 2003, Attenuation of strong ground motion and site specific seismicity in Peninsular India. In Proceedings of National Seminar on Seismic Design of Nuclear Power Plants, 21–22 February 2003, SERC, Chennai (Delhi: Allied Publishers Pvt. Ltd.), pp. 269–291. IYENGAR, R.N. and RAGHUKANTH, S.T.G., 2006, Strong ground motion estimation during the Kutch, India earthquake. Pure Applied Geophysics, 163, pp. 153–173. KOMATITSCH, D. and TROMP, J., 2001, Modeling of seismic wave propagation at the scale of the earth on a large Beowulf. In Proceedings of the ACM/IEEE supercomputing SC’2001 conference, CD-ROM. KOMATITSCH, D. and TROMP, J., 2002a, Spectral-element simulations of global seismic wave propagation–I. Validation. Geophysical Journal International, 149, pp. 390–412. KOMATITSCH, D. and TROMP, J., 2002b, Spectral-element simulations of global seismic wave propagation–II. 3-D models, oceans, rotation, and self-gravitation. Geophysical Journal International, 150, pp. 303–318. KUSTOWSKI, B., EKSTROM, G. and DZIEWONSKI, A.M., 2008, Anisotropic shear-wave velocity structure of the Earth’s mantle: a global model. Journal of Geophysical Research, 113, pp. B06306, doi:10.1029/ 2007JB005169. MAI, P.M. and BEROZA, G.C., 2002, A spatial random field model to characterize complexity in earthquake slip. Journal of Geophysical Research, 107, pp. 2308–2329, doi:10.1029/ 2001JB000588. MAI, P.M., SPUDICH, P. and BOATWRIGHT, J., 2005, Hypocenter locations in finite-source rupture models. Bulletin of the Seismological Society of America, 95, pp. 965–980. MONSALVE, G., SHEEHAN, V., SCHULTE-PELKUM, V., RAJAURE, S., PANDEY, M.R. and WU,F., 2006, Seismicity and one-dimensional velocity structure of the Himalayan collision zone: earthquakes in the upper crust and upper mantle. Journal of Geophysical Research, 111, pp. B10301, doi:10.1029/2005JB004062. MOTAZEDIAN,D.and ATKINSON, G.M., 2005, Stochastic finite-fault modeling based on a dynamic corner frequency. Bulletin of the Seismological Society of America, 95, pp. 995–1010. PARVEZ, I.A., VACCARI, F. and PANZA, G.F., 2003, Deterministic seismic hazard map of India and adjacent areas. Geophysical Journal International, 55, pp. 489–508. PATERA, A.T., 1984, A spectral element method for fluid dynamics: laminar flow in a channel expansion. Journal of Computational Physics, 54, pp. 468–488. RAGHUKANTH, S.T.G., 2008, Ground motion estimation during the Kashmir earthquake of 8th October 2005. Natural Hazards, 46, pp. 1–13. RAGHUKANTH, S.T.G. and BHANUTEJA, B., 2011, Ground motion simulation for 26th January 2001 Gujarat earthquake by spectral finite element method. Journal of Earthquake Engineering, (In Press). 34 S.T.G. Raghukanth et al. RAGHUKANTH, S.T.G. and DASH, S.K., 2009, Deterministic seismic scenarios for North East India. Journal of Seismology, 14, pp. 143–167. RAGHUKANTH, S.T.G. and SOMALA, S.N., 2009, Modeling of strong motion data in Northeastern India: Q, Stress Drop and Site Amplification. Bulletin of the Seismological Society of America, 99, pp. 705–725. RAY, S.K., 2000, Culmination zones in Eastern Himalaya. Geological Survey of India (Special Publication), 55, pp. 85–94. SARAGONI, G.R. and HART, G.C., 1974, Simulation of artificial earthquakes. Earthquake Engineering and Structural Dynamics, 2, pp. 249–267. SHINOZUKA, M. and DEODATIS, G., 1996, Simulation of multi-dimensional gaussian stochastic fields by spectral representation. Applied Mechanics Reviews – ASME, 49, pp. 29–53. SINGH, S.K., GARCIA, D., PACHECO, J.F., VALENZUELA, R., BANSAL,B.K.and DATTATRAYAM, R.S., 2004, Q of the Indian shield. Bulletin of the Seismological Society of America, 94, pp. 1564–1570. SINGH, S.K., ORDAZ, M., DATTATRAYAM, R.S. and GUPTA, H.K., 1999, A spectral analysis of the 21 May 1997, Jabalpur, India, earthquake (Mw ¼ 5.8) and estimation of ground motion from future earthquakes in the Indian shield region. Bulletin of the Seismological Society of America, 89, pp. 1620–1630. SOMERVILLE, P.G., IRIKURA, K., GRAVES, R., SAWADA, S., WALD, D., ABRAHAMSON,N., IWASAKI, Y., KAGAWA, T., SMITH, N. and KOWADA, A., 1999, Characterizing crustal earthquake slip models for the prediction of strong ground motion. Seismological Research Letters, 70, pp. 59–80. WALD, D.J. and ALLEN, T.I., 2007, Topographic slope as a proxy for seismic site conditions and amplification. Bulletin of the Seismological Society of America, 97, pp. 1379–1395. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png "Geomatics, Natural Hazards and Risk" Taylor & Francis

Estimation of ground motion during the 18th September 2011 Sikkim Earthquake

Loading next page...
 
/lp/taylor-francis/estimation-of-ground-motion-during-the-18th-september-2011-sikkim-Hb0heSJxzb

References (47)

Publisher
Taylor & Francis
Copyright
Copyright Taylor & Francis Group, LLC
ISSN
1947-5713
eISSN
1947-5705
DOI
10.1080/19475705.2011.646323
Publisher site
See Article on Publisher Site

Abstract

Geomatics, Natural Hazards and Risk Vol. 3, No. 1, February 2012, 9–34 Original article Estimation of ground motion during the 18th September 2011 Sikkim Earthquake S.T.G. RAGHUKANTH*, K. LAKSHMI KUMARI, and B. KAVITHA Department of Civil Engineering, IIT Madras, Chennai, India (Received 30 October 2011; in final form 30 November 2011) This article simulates the ground motion during the devastating 18th September 2011 Sikkim earthquake by analytical approaches. Surface level peak ground acceleration (PGA) values at several stations in the epicentral region have been estimated by the finite fault seismological model. A PGA contour map in the epicentral region is provided. It is found that very near to the epicentre, PGA would have reached more than 0.3g. The SPECFEM3D_GLOBE package which models the complete three-dimensional wave velocity structure and topography is used to simulate the displacement time histories on a global scale. A contour map of peak ground displacement (PGD) and permanent ground residual displace- ments (PGRD) in the near source region is provided. The estimated ground motions are validated with that obtained from damage values. 1. Introduction The state of Sikkim with Gangtok (27.338N, 88.68E) as the capital city lies in the Himalayan region (Figure 1). There has been an increase in population and infrastructure developments in the state of Sikkim in last few decades. As per census 2011, the total population of Sikkim is 0.6 million which is about 12.36% growth from the last census of 2001. The state of Sikkim comprising of four districts covers an area of about 7096 sq km with the population density registered at 86 persons per sq km. The Indian standard code of practice for earthquake resistant design of structures, IS 1893 (2002), has identified Sikkim state in zone of high seismic hazard (Zone IV) with the design peak ground acceleration as 0.24g. The 18th September 2011 event (M 6.9) which occurred on the Indo–Nepal border caused severe damages in the state of Sikkim. According to the India Meteorological Department (IMD) the epicentre of this quake was at 27.718N, 88.28E. This earthquake was widely felt in major cities, including Kolkata, Dhaka, Kathmandu and Chittagong. Several strong earthquakes of magnitudes greater than M 6 have occurred in Sikkim Himalaya in the past, but none of them have lead to a large scale damages compared to the 18th September 2011 event. This can be attributed to the increase in the growth of population and infrastructure development in this region over the last few decades. Although damages were reported in the entire Sikkim, towns and villages in north Sikkim district were worst affected due to proximity of the epicentre. This event triggered several landslides and rolling boulders in the *Corresponding author. Email: raghukanth@iitm.ac.in Geomatics, Natural Hazards and Risk ISSN 1947-5705 Print/ISSN 1947-5713 Online ª 2012 Taylor & Francis http://www.tandf.co.uk/journals http://dx.doi.org/10.1080/19475705.2011.646323 10 S.T.G. Raghukanth et al. Figure 1. Sikkim state. Red dots indicate location of landslides. Available in colour online. epicentral region. A total of 354 landslides were mapped in Teesta valley by National Remote Sensing Centre (http://www.nrsc.gov.in/) from satellite images after the earthquake. The strong shaking triggered landslides which caused damage to lifeline structures and significant building collapse. The official reports indicate that about 97 people were killed in Sikkim state. Casualties have been reported in neighbouring Nepal, Bhutan and Tibet. The financial losses due to this event are estimated to be around US$ 20 billion. The damage during the earthquake is mainly attributed to the strong near-source ground motion, surface ruptures and widespread slope failures in the epicentral region. For engineering purposes, the most sought-after information during an earthquake is the ground motion time history in the epicentral region. In this study, analytical and stochastic seismological approaches are used to estimate the ground motion during this event. These models are viable alternatives and popular among engineers and seismologists for computing the ground motion during earthquakes in regions lacking Strong Motion Accelerogram (SMA) data. These models have been successfully used to estimate ground motion during several past earthquakes (Singh et al. 1999, Iyengar and Raghukanth 2006, Raghukanth 2008, Raghukanth and Bhanuteja 2011). The regional velocity model, quality factor and stress drop are available for the Himalayan region (Parvez et al. 2002, Singh et al. 2004). With this limited information, one can use seismological approaches to estimate ground motion during this event. In this study, estimates of peak ground acceleration (PGA) values for the Sikkim event have been carried out by three different approaches. In the first approach, the finite fault seismological model of Boore (2009) is used for simulating acceleration time histories. Acceleration time histories are computed on grid size of 1 sq km covering a region of 200 km6 200 km in Sikkim region corresponding to hard rock type site condition. The site classification of Sikkim region based on topography as reported in www.earthquake.usgs.gov/hazards/apps/ vs30/ website combined with the correction factors available in IBC (2003) is used for estimating the surface level PGA. In the second approach, spectral finite element method (SPECFEM) is used to simulate ground motion. SPECFEM is based upon a Estimation of ground motion 11 weak formulation of the equations of motion and combines the flexibility of a finite- element method with the accuracy of a global pseudospectral method. The software package, SPECFEM3D_GLOBE which is a collection of Fortran subroutines for modelling the global seismic wave propagation is used to compute the displacement time histories in the epicentral region. Apart from these two analytical approaches, an estimate of PGA is obtained from the reported MMI values. 2. Seismotectonic setting and source of Sikkim earthquake The seismotectonic map of Sikkim Himalaya with superimposed seismicity is shown in Figure 2. It can be observed that two main Himalayan faults namely main boundary thrust (MBT) and main central thrust (MCT) are crossing the Sikkim state (Geological Survey of India (GSI) 2000). The MCT is not parallel to the MBT and arches in the form of a mushroom (Ray 2000). Apart from this, MCT also takes closed forms in Sikkim Himalaya. A large number of lineaments which extend for several kilometers lie transverse to the Himalayan faults in this region. The Tista Figure 2. Seismotectonic map of Sikkim Himalaya (GSI 2000). Available in colour online. 12 S.T.G. Raghukanth et al. lineament is thought to separate central and eastern Himalaya. The east Patna fault extends in to the Himalaya in the form of Arun lineament. The Kanchenjunga lineament which lies parallel to Arun lineament also extends from Indo-Gangetic plain to well inside the Himalayan region. The Tista and Purnia-Everest lineaments exhibit NW-SE trending where as Kanchenjunga, Arun and Everest lineament exhibit NE-SW trending. Detailed analysis of the seismicity in the Sikkim Himalaya indicates that movements along faults like the Gangtok lineament and the Tista lineament are responsible for many earthquakes in this region. De and Kayal (2003) analysed microearthquake data and observed that Sikkim Himalaya exhibit complex tectonic activity. The focal depth of the earthquakes in Sikkim Himalaya lies in between 0 and 45 km and occur by both thrust and strike-slip mechanism. Monsalve et al (2006) studied the seismicity in eastern Nepal and reported a bimodal distribution of focal depth in the Himalaya. They concluded that earthquakes in Himalaya occur both in the upper crust and around the crust-mantle boundary in the 50–100 km depth range. Recently, Hazarika et al. (2010) studied 356 local earthquakes in Sikkim and reported that seismic activity in this region can be attributed to transverse tectonics rather than underthrusting of India plate. The focal depth of the events in Sikkim Himalaya could be as high as 70 km and moment tensors indicate predominantly strike slip faulting. They attributed the strike slip mechanism in Himalayan thrust environment to an oblique convergence of Indian- Asian collision. Sikkim has experienced several damaging earthquakes in the past. The MMI in Sikkim Himalaya has been reported to be VIII for the great Bihar-Nepal earthquake (Mw 8.3) of 1934 (GSI 2000). The 1988 event (Mw 6.5) which occurred in the epicentral region of 1934 great earthquake caused landslides and subsidence in Sikkim (De and Kayal 2003). Damages were observed in East Sikkim for the recent 14th February 2006 (Mw 5.3) moderate event. The focal mechanism solution and field investigations indicate that the Tista lineament caused 18th September 2011 event. The USGS (http://neic.usgs.gov/) website reported the strike, dip and rake angle of the rupture plane as 1308,908 and 1688, respectively. The focal depth of this event is obtained as 35 km. The Global Centroid Moment Tensor database (http:// www.globalcmt.org/) obtained the focal depth as 47 km with strike, dip and rake angle as 1338,738 and 1638. The focal mechanism solutions reported by the above two agencies indicate nearly strike slip faulting. On the other hand, the India Meteorological Department (IMD) estimated the focal depth of this event as 10 km. 3. Simulation of SMA using stochastic finite-fault simulation method Due to unavailability of strong motion records for the Sikkim earthquake, ground motions are simulated by an analytical model. The stochastic finite fault approach of Boore (2009) which is an improved version of the methodology proposed by Motazedian and Atkinson (2005) is used to simulate the rock level ground motion for the seven geological regions in India. The theory and application of seismological models for estimating ground motion has been discussed in detail by Boore (1983, 2009) and Motazedian and Atkinson (2005). The application of seismological model for estimating ground motion in northeastern India has been discussed in detail by Raghukanth and Somala (2009) and Raghukanth and Dash (2009). A brief description of the method starting with the frequency domain representation of the ground acceleration will bring out the advantages of this approach. In this method, Estimation of ground motion 13 the rectangular fault plane is divided into N number of subfaults and each subfault is represented as a point source. The Fourier amplitude spectrum of ground motion [A(r,f)] due to the jth subfault at a site is derived from the point source seismological model, expressed as: A ðr; fÞ¼ CH S ðfÞFðfÞDðr; fÞPðfÞ: ð1Þ j j j Here, C is a scaling factor, S (f) is the source spectral function, D(r,f) is the diminution function characterizing the quality of the region, P(f) is a filter function, F(f) is the site dependent function that modifies the bed rock motion in the vertical direction and H is a scaling factor used for conserving the energy of high-frequency spectral level of subfaults. In this study, following Brune (1970), the principal source model S (f) is taken as: 2 0j hi S ðfÞ¼ð2pf Þ : ð2Þ 1 þ f f 0j Here, f is the corner frequency and M is the seismic moment of the jth subfault. 0j 0j The three important seismic source parameters M , f and the stress drop (Ds) are 0j 0j related through: 1=3 Ds 1=3 6 1=3 f ¼ 4:9  10 ðN Þ N V : ð3Þ 0j s Rj 0j Here, V stands for the shear wave velocity in the source region, corresponding to bedrock conditions. N is the cumulative number of ruptured subfaults by the time Rj rupture reaches the jth subfault. The spatial modifying function D(r, f) is given by: pfr Dðr; fÞ¼ Gexp ; ð4Þ V QðfÞ where G is the geometric attenuation factor. The other term denotes anelastic attenuation with hypocentral distance r and the quality factor as Q. The spatial spread of the ground motion depends sensitively on the quality factor of the local region. The constant C of Equation (1) is: pffiffiffi R 2 yf C ¼ ð5Þ 4prV where hR i is the radiation coefficient averaged over an appropriate range of yf azimuths and take-off angles and r is the material density at the focal depth. The coefficient 2 in the above equation arises as the product of the free surface amplification and partitioning of energy in orthogonal directions. The scaling factor for jth sub-fault, H based on the squared acceleration spectrum is taken as (Boore 2009) pffiffiffiffi H ¼ N ; ð6Þ 0j 14 S.T.G. Raghukanth et al. where f is the corner frequency at the end of the rupture, which can be obtained by substituting N (t) ¼ N in Equation (3). The filter function P (f) is taken as: R j pffiffiffiffi 1 þðf f Þ 0j P ðfÞ¼ : ð7Þ j 1 þðf=f Þ The moment of jth subfault is computed from the slip distribution as follows: M D 0 j M ¼ ð8Þ 0j j¼1 Here, D is the average final slip acting on the jth subfault. M is the total seismic j 0 moment on the fault. To further account for earthquake rupture, Motazedian and Atkinson (2005) introduced the concept of pulsing area where the cumulative number of active subfaults, N increases with time at the initiation of rupture and Rj becomes constant at some fixed percentage of the total rupture area. This parameter determines the number of active subfaults during the rupture of jth subfault. These many subfaults are used in computing the corner frequency in Equation 3. The above is a general finite source model expressed in the frequency domain, valid for any region if only the various controlling parameters can be selected suitably. The effects that are specific to Sikkim earthquake appear in the quality factor Q(f), stress drop value (Ds), focal depth, dip angle of the fault plane and the site amplification function F(f). Thus, these five parameters have to be selected to simulate ground motion for this event. Here lies strength of this approach since almost all required parameters for Himalayan region are available. The spatial modifier of Equation 4 consists of a general geometrical attenuation term G. This is taken following Singh et al. (1999) as: ; for r  100 G ¼ ð9Þ pffiffi ; for r > 100 10 r The shear wave velocity and density at the focal depth are fixed at 3.5 km/s and 2900 kg/m , respectively corresponding to bedrock (Singh et al. 1999, Hazarakia et al. 2010). In the past, several seismologists have analysed instrumental data to arrive at the value of stress drop that can happen in Himalayan region for a given magnitude of earthquake. Singh et al. (2004) have studied the 1991 Uttarkashi earthquake and 1999 Chamoli earthquake and derived the frequency dependent 0.8 quality factor for Himalayan region as 253f . The other important regional parameter is the stress drop which can vary from 50 to 200 bars in Himalayan region (Singh et al. 2004). Since there are discrepancies regarding the focal depth for this 18th September event, it is varied from 25 km to 50 km. The regional velocity model used is the one reported by Parvez et al. (2003). Details of the regional model are presented in Table 1. Modification between bedrock and A-type site is a linear problem in one-dimension and hence for such sites, amplification can be directly found by the quarter-wavelength method of Boore and Joyner (1997). The estimated amplification function is presented in Figure 3. Apart from the above parameters, the S-wave radiation coefficient (hR i) varies randomly within particular intervals. yf Here, following Boore and Boatwright (1984), the S-wave radiation coefficient is taken to be in the interval 0.48–0.64. In addition to the above parameters the slip Estimation of ground motion 15 Table 1. Regional Velocity Model (Parvez et al 2003). Depth (km) V (km/s) V (km/s) Density (kg/m ) Q Q p s p S 0–4.9 4.5 2.5 1800 3600 1600 4.9–19 5.9 3.1 2400 2200 1100 19–32 6.1 3.5 2600 1400 500 32–50 7.1 4.1 2900 800 200 4 50 7.6 4.2 3300 800 200 Figure 3. Amplification function for A-type sites in Sikkim region. distribution and pulsing percentage area is also required in the simulation. The pulsing percentage area is varied from 25% to 75% (Atkinson and Boore 2006). 3.1. Slip model The stochastic finite fault approach requires slip distribution on the rupture plane which is highly complex. Mai and Beroza (2002) have proposed stochastic characterization of slip, in which slip distribution is modelled as a anisotropic Von Karman random field. They modelled slip distributions of several past earthquakes and developed empirical equations for estimating correlation lengths along the two dimensions of the fault plane from the magnitude of the earthquake. In this approach, once the dimensions and orientation of rupture plane are fixed, the slip distribution on the fault plane is generated as a spatial random field with correlation lengths depending on the magnitude of the earthquake. Accordingly, in the present analysis, the correlation lengths for these four events in strike and dip direction are determined from the scaling relations of Mai and Beroza (2002). From the derived correlation lengths, slip is simulated as a random field by the spectral representation method of Shinozuka and Deodatis (1996). The hypocentre in each case is taken nearer to the regions of maximum slip release (Mai et al. 2005). This sample slip field on the fault plane is discretized into subfaults of size 1 km6 1 km and seismic moment of each point source is computed using Equation 8. Since the moment 16 S.T.G. Raghukanth et al. distribution on the rupture plane is modelled as random field, one can generate an ensemble of time histories at the surface by numerical simulation. Based on the scaling relations of Mai and Beroza (2002), the dimensions of the rupture plane are estimated from magnitude as 35 km6 20 km. In Figure 4, two samples of the slip field are shown. To cover all possible peculiarities in the slip distribution, fifty samples of slip field are generated for each event. 3.2. Sample ground motion From the above discussion it is seen that once the stress drop, slip field, site amplification and Q factors are known, the Fourier amplitude spectrum of ground acceleration for any combination of magnitude (M ) and hypocentral distance (r) can be expressed within the limitations of a finite source seismological source mechanism model. This is a random process and hence in the time domain this represents an ensemble of accelerograms. These samples can be retrieved from Equation 1 in three steps. First, a Gaussian stationary white noise sample of length equal to the strong motion duration (Boore and Atkinson 1987) T ¼ðÞ 1=f þ 0:05 r ; ð10Þ Figure 4. Sample slip realizations for the Sikkim earthquake (M 6.9). Available in colour online. Estimation of ground motion 17 is simulated for each subfault. This sample is multiplied by a suitable non-stationary modulating function suggested by Saragoni and Hart (1974). The simulated sample is Fourier transformed into the frequency domain. This Fourier spectrum is normalized by its root mean-square value and multiplied by the seismological source-path-site function A(f) of Equation 1. The resulting function is transformed back into the time domain to get a sample of the ground motion accelerogram for each subfault. The simulated acceleration time histories for all the subfaults are summed up with time delay of Dt to account for the rupture velocity to obtain the final ground acceleration, aðtÞ¼ a t þ Dt : ð11Þ j j j¼1 Any number of accelerograms can be simulated in this fashion to compute samples of response spectra corresponding to 5% damping. Since the stress drop, focal depth, dip, amplification function, radiation coefficient, pulsing percentage area are random variables and slip as a random field, the uncertainty arising out of these parameters is included in the simulation. The input parameters used in the seismological model for the Sikkim event are reported in Table 2. Accordingly, 50 samples of these seven input parameters are generated from uniform distribution and these are combined using the Latin Hypercube sampling technique (Iman and Conover 1980) to simulate 50 sets of random seismic parameters for the Sikkim event. By the above procedure, acceleration time histories are computed on grid size of 1 sq km covering a region of 200 km6 200 km in Sikkim region. In Figure 5(a,b) mean and standard deviation of rock level peak ground acceleration are shown. It can be observed that the contours are axisymmetric and the PGA is of the order of 0.22g+ 0.06g near the epicentre. The standard deviation of PGA decreases with increase in epicentral distance. Table 2. Input Parameters in the stochastic finite fault seismological model. Parameter Sikkim earthquake Mag., M 6.9 Location Latitude 27.71 Longitude 88.20 Depth to top of the fault (km) 25–50 Stress drop (bars) 50–200 0.8 Q( f ) 253f Site effects Quarter-wavelength method Fault dimensions (km) 356 19 Fault orientation Strike(8) 120–140 Dip (8) 70–90 Pulsing area percentage 25–60 Rupture velocity 0.8b Slip distribution Random field 18 S.T.G. Raghukanth et al. Figure 5. (a) Contour map of mean PGA (g) at A-type rock. Black dots indicate location of landslides; (b) Standard deviation of PGA (g) at A-type rock. Available in colour online. Estimation of ground motion 19 3.3. Surface level PGA The local site property is expressed in terms of the average shear wave velocity in the top 30 m of the soil. This is the recognized IBC (2003) approach, wherein sites are classified as A: (V 4 1.5 km/s); B: (0.76 km/s5V  1.5 km/s); C: (0.36 km/s 5 30 30 V  0.76 km/s); D: (0.18 km/s 5 V  0.36 km/s). E- and F-type sites with V 30 30 30 0.18 km/s are susceptible for liquefaction and failure. The above simulated PGA values are valid for A-type rock site condition. The stochastic finite fault simulation described above can be extended to the surface with the help of soil profiles and V values sampled from the region. Thus, for a specific site, precise correction of the bedrock results will be possible only when the soil section data are available, including the variation of V with depth. Unfortunately, such information on lateral variation of shear wave velocity profiles in Sikkim is not available. Wald and Allen (2007) developed empirical equations relating slope of the topography variations with shear wave velocity from global database. In this study, the approach proposed by Wald and Allen (2007) is used for classifying the local site condition in the epicentral region. The variation of shear wave velocity in the top 30 m (V ) obtained from global V map server (http://earthquake. usgs.gov/hazards/apps/vs30/) is shown in Figure 6. It can be observed that most of the sites correspond to C-type site Figure 6. Shear wave velocity in top 30 m (V , m/s), (http://earthquake.usgs.gov/hazards/ apps/vs30/). Available in colour online. 20 S.T.G. Raghukanth et al. condition. Once V is known, the amplification factors given in IBC (2003) listed in Tables 3 and 4 are used to estimate the surface level PGA and 1 s spectral acceleration. In Figure 7 the resulting surface level PGA contour is shown along with the location of landslides. It can be observed that the contours are axisymmetric which can be attributed to the focal depth of the earthquake. The effect of the local site condition on PGA can also be observed in the figure. The PGA in the near source region is as high as 0.35g. At Lachen, Ravangla, Teesta, Chungthang and Mangan, PGA lies in between 0.25g and 0.29g. PGA estimates obtained from the stochastic finite fault simulation at important villages and towns which suffered damage are shown in Table 5. The contours of spectral acceleration (S )at 1 s natural time period are shown in Figure 8. The spectral acceleration in the epicentral region is as high as 0.2g. At places like Chungthang and Mangan which suffered heavy damage it is of the order of 0.14g. 4. PGA from MMI data The PGA values estimated above are based on an analytical model. The above results are now correlated with field observation of damage after the event. Although strong motion data are not available, information on damages at several villages and towns in epicentral region are available from filed investigations and media reports. The MMI values are fixed up based on these damages. These intensities are reported in Table 5. It can be observed that maximum intensity of VIII has been felt at Chungthang.To validate the estimates obtained from the seismological model, the MMI values have been converted into PGA values using an empirical relation. Previously, Iyengar and Raghukanth (2003) compared Indian MMI data with instrumental PGA data with the help of 43 sample values and obtained a linear relation between ln(PGA) and MMI as: lnðÞ PGA=g¼ 0:6782 MMI  6:8163; s lnðÞ e ¼ 0:7311: ð12Þ Table 3. Short period correction factor for various site conditions (IBC 2003). Site category 0.1g 0.2g 0.3g 0.4g 0.5g A 0.8 0.8 0.8 0.8 0.8 B 1.0 1.0 1.0 1.0 1.0 C 1.2 1.2 1.1 1.0 1.0 D 1.6 1.4 1.2 1.1 1.0 Table 4. Long period correction factor for various site conditions (IBC 2003). Site category 0.1g 0.2g 0.3g 0.4g 0.5g A 0.8 0.8 0.8 0.8 0.8 B 1.0 1.0 1.0 1.0 1.0 C 1.7 1.6 1.5 1.4 1.3 D 2.4 2.0 1.8 1.6 1.5 Estimation of ground motion 21 Figure 7. Surface level peak ground acceleration (g). Black dots indicate location of landslides. Available in colour online. The above equation is for the expected PGA values only. In practice, PGA will be distributed as a random variable and hence recorded values will show considerable variation in a region with nearly constant MMI. In Table 5, the MMI-PGA values are compared to the corresponding PGA values estimated from stochastic approaches. It can be observed that the comparison between the two estimates is consistent. 5. Displacement time history The above one-dimensional stochastic finite fault seismological model is efficient in simulating accelerations including high frequencies. However, this model do not incorporate seismic phases that partition the energy into high frequency body waves and low frequency surface waves and hence is not valid in simulating displacement time histories, which are controlled by low frequencies. Moreover, the displacements are sensitive to fault rupture process. The three-component displacement time histories including the near-source features like directivity and fling are required to understand the evolution of nonlinear structural deformations in large bridges, dams and other important structures. Since many hydro-power projects are under construction in Sikkim state, it would be of interest to engineers to know the surface displacements caused by this moderate event. Analytical techniques based on kinematic models and elastic wave propagation approaches are very efficient in 22 S.T.G. Raghukanth et al. Table 5. Comparison among various PGA estimates. Latitude Longitude MMI–PGA PGA (g) Station (8N) (8E) MMI Reference (g) Stoch. fault Lachen 27.73 88.55 VII India today 0.13 0.297 Lachung 27.7 88.75 VII India today 0.13 0.226 Chungthang 27.61 88.63 VIII The hindu; earthquake-report.com 0.25 0.268 Mangan 27.52 88.54 VIII news.oneindia.in 0.25 0.281 Ravangla 27.30 88.36 VI asc-india.org 0.06 0.248 Teesta 27.03 88.26 VII The Hindu 0.13 0.150 Dikchu 27.38 88.52 VII asc-india.org 0.13 0.248 Rangpo 27.18 88.53 VI The Economic Times 0.06 0.181 Taplejung 27.35 87.66 V earthquake-report.com 0.03 0.180 Yuksom 27.37 88.22 VI asc-india.org 0.06 0.279 Singtam 27.23 88.50 VII news.oneindia.in 0.13 0.223 Gangtok 27.20 88.40 VII www.reuters.com 0.13 0.199 Darjeeling 27.05 88.26 VII Hindustan Times 0.13 0.156 Kalimpong 27.07 88.48 VII NDTV 0.13 0.156 jorethang 27.13 88.28 VII haalkhabar.in 0.13 0.186 Pentong 27.51 88.53 VI The Times of India 0.06 0.278 Yadong 27.51 88.97 V earthquake-report.com 0.03 0.150 Sankhuwasabha 27.58 87.33 VII Nepal Mountain News 0.13 0.136 Sakyong 27.33 88.60 VII The Times of India 0.13 0.209 Nayabazaar 27.13 88.26 VI haalkhabar.in 0.06 0.185 Pegong 27.59 88.65 VIII The Economic Times 0.25 0.255 Lingzya 27.55 88.45 VIII The Times of India 0.25 0.301 Estimation of ground motion 23 Figure 8. Surface level Spectral acceleration (g) Contours at T ¼ 1 s. Available in colour online. explaining source directivity effects as well as effects of underground medium on ground motion. The Sikkim event also triggered extensive landslides in the epicentral region which is mainly due to varying Himalayan topography (Figure 1). This demands simulation of ground motion by including the topography of the region. In this study, spectral finite element method (SPECFEM) is used to simulate ground motion for the Sikkim earthquake. SPECFEM is based upon a weak formulation of the equations of motion and combines the flexibility of a finite-element method with the accuracy of a global pseudospectral method. The software package, SPEC- FEM3D_GLOBE which is a collection of Fortran subroutines for modelling the global seismic wave propagation is used to compute the displacement time histories in the epicentral region. The effect of lateral variations in material properties, topography and bathymetry, the oceans, have been incorporated in simulating ground motion. 5.1. Spectral finite element method In this study, the spectral finite element method introduced by Patera (1984) and developed by Komatitsch and Tromp (2002a,b) for global seismic wave propagation is used to model the Earth for simulating the ground motion during the Sikkim event. The Earth’s interior can be broadly divided into four layers namely crust, mantle, outer core and inner core. The crust, mantle and inner core are solid where as the outer core is composed of liquid material. 24 S.T.G. Raghukanth et al. 5.1.1. Crust and mantle. The governing equation for the Earth medium in spherical domain (r,y,f) subjected to a unit impulsive double couple applied at location (r )in the crust is given by (Dahlen and Tromp 1998) @ uðr; tÞ @uðr; tÞ rðrÞ þ 2o  ¼ H uðr; tÞþ fðr; tÞð13Þ @t @t where u(r,t) is the displacement field, r(r) is the material density, o is the Earth’s angular rotation vector and f(r,t) denotes the body force. The elasto-gravity operator, H is given by: H uðr; tÞ¼r  sðuÞþrðru  gÞ  rðruÞg: ð14Þ In the above equation, g is the gravitational acceleration. The body force components for a unit impulsive double couple are: fðr; tÞ¼M rdðr  r ÞSðtÞð15Þ where d denotes the Dirac delta function, S(t) denotes the source time function and M contains the components of the symmetric seismic moment tensor related to the orientation of the fault plane as (Aki and Richards 1980) M ¼M sin ðdÞcos ðlÞsin ð2fÞþ sin ð2dÞsin ðlÞsin ðfÞ yy 0 M ¼MðÞ sin ðdÞcos ðlÞcos ð2fÞþ 0:5sin ð2dÞsin ðlÞsin ð2fÞ yj 0 M ¼MðÞ cos ðdÞcos ðlÞcos ðfÞþ cos ð2dÞsin ðlÞsin ðfÞ ry 0 ð16Þ M ¼ M sin ðdÞcos ðlÞsin ð2fÞ sin ð2dÞsin ðlÞcos ðfÞ jj 0 M ¼ MðÞ cos ðdÞcos ðlÞsin ðfÞ cos ð2dÞsin ðlÞcos ðfÞ rj 0 M ¼ MðÞ sin ð2dÞsin ðlÞ rr 0 where f, d and l are the strike, dip and rake angles of the fault. The seismic moment (M ) is estimated from the average slip (D), fault area (A) and rigidity modulus (m) as: M ¼ mAD ð17Þ 5.1.2. Outer core. Assuming the fluid to be stably stratified and isentropic, the equation of motion in the fluid outer core can be written as: @ uðr; tÞ @uðr; tÞ 1 rðrÞ þ 2o  ¼r kr u þ u  g ð18Þ @t @t r where k is the fluid’s bulk modulus. 5.1.3. Inner core. The equation of motion for inner core, which is composed of solid material similar to crust and mantle region can be written as: @ uðr; tÞ @uðr; tÞ rðrÞ þ 2o  ¼ H uðr; tÞð19Þ @t @t Estimation of ground motion 25 5.1.4. Boundary conditions. The vanishing of stresses at the free surface on continental crust can be expressed as: sðÞ r; t nrðÞ ¼ 0 8r 2  ; ð20Þ where  is the free surface of the continental crust and n being the unit normal to the interface. To model the oceans and large lakes on Earth’s surface, the stress free boundary condition at the ocean crust boundary (OCB) is modified to include the fluid pressure. Neglecting the variations of gravity and angular rotation in the oceans, if h is the height of the water column, the fluid pressure at the ocean crust boundary can be written as: @ uðr; tÞ sðr; tÞ nðrÞ¼ r h nðrÞ 8r 2  ; ð21Þ OCB @t with r being the density of the sea water. At crust-mantle boundary (CMB) and inner core boundary (ICB) normal component of velocity and normal stress have to be continuous. This can be expressed as: @u nrðÞ ¼ 0 8r 2 ICB; CMB @t ð22Þ ½ jsðÞ r; t  nðrÞj¼ 0 8r 2 ICB; CMB The equations of motion in crust and mantle, outer core and inner core (Equations 13, 18 and 19) has to be solved subjected to the above boundary conditions (Equations 21–22). The spectral finite element method is based on the integral or weak formulation of the problem. Earth is divided into a number of hexahedral non- overlapping elements. In SPECFEM high-degree Lagrange interpolant is used to represent functions on the element. In the hexahedral volume element, all the functions are interpolated by triple products of Lagrange polynomials of degree n at the classical n þ1 nodes called as Gauss-Lobatto-Legendre (GLL) quadrature points. The numerical integration of these functions over volume elements are approximated using the GLL integration rule. The contributions from all the elements that share a common global grid point are summed before the time marching. The method is implemented in parallel computing using the message passing technique. (Komatitsch and Tromp 2001). The SPECFEM3D Globe package (Bozdag et al. 2010) which is a collection of Fortran subroutines is used to simulate three-dimensional global wave propagation for the Gujarat earthquake. The SPECFEM3D package has been implemented on the VEGA super cluster available at IIT Madras using job scheduling software ‘portable batch system’ (Raghukanth and Bhanuteja 2011). The VEGA super cluster is based on HP Proliant DL160 G5 servers with Quad-core Intel Xeon processors. This Cluster consists of 128 nodes. Each node consists of eight processors. The first step in SPECFEM is to generate the spectral-element mesh. This is shown in Figure 9. The Earth is divided into six chunks based on the cubed sphere mapping. In this study, to allow for a denser mesh, among the six chunks that constitute the globe, single chunk mesh is used for estimating the ground motion. The chunk is centered at 308N and 808E. The size of the chunk is 6086 608 with absorbing boundary conditions at the artificial edges of the mesh (Figure 9). This chunk is divided into 256 slices with each slice 26 S.T.G. Raghukanth et al. Figure 9. Spectral element mesh of the globe. Red line indicates the size of the chunk used in this study. Available in colour online. assigned to one processor. The total number of processors needed to handle the mesh would be 256. Each slice is further subdivided into 406 40 spectral elements at the earth’s surface. The total number of spectral elements in the one chunk mesh is 13.2 million. Each spectral element consists of 125 grid points. Since grid points are shared among the elements, the entire globe is represented by 872.2 million grid points. The total number of degrees of freedom in entire mesh is 2.31 billion. The average grid spacing at the Earth’s surface is about 2.61 km. The simulated displacement time histories are valid up to a shortest period of 4.5 s. Once all the three regions in Earth’s interior have been discretized, the global system of equations to be solved can be written as: M€u þ Wu_ þ Ku þ Bu ¼ F ð23Þ where M and K are the global mass matrix and stiffness matrix. W contains terms related to angular rotation vector and B is related to the boundary interactions at Estimation of ground motion 27 CMB and ICB and F is the source term. Explicit expressions for the mass, stiffness, W and B matrices at elemental level and further construction of these matrices at global level are available in the article of Komatitsch and Tromp (2002a,b). A explicit second-order finite difference scheme is used to march the above equation in time. The finite source rupture model reported in Figure 4 is used for simulating ground motion. The 3D velocity model for the mantle and crust are taken from the work of Kustowski et al. (2008) and Bassin et al. (2000). The topography and bathymetry model available at 5-min interval grid (ETOPO5) is from the US National Oceanic and Atmospheric Administration (http://www.ngdc.noaa.gov/). The 35 km6 19 km fault plane is divided into 665 subfaults and slip acting on each subfault is represented as a point double-couple varying in time as a smooth ramp function. The moment tensor at each subfault is computed from slip and orientation of the fault plane (Equations 15 and 16). The rise time in the ramp function for each subfault is estimated from the empirical relation of Somerville et al (1999) as: 1=3 T ¼ 2:03  10 M ð24Þ The triggering time of each point source is the estimated from the ratio of the distance between the point source and the hypocentre to rupture velocity. The three components of displacements on the surface of the Earth, for each double couple source are obtained by solving the Equation 23. The ground displacements due to the complete fault rupture, is obtained by summing-up the contribution of all the point double-couples. A typical ground motion simulation including the mesher and solver requires about 19 h using 256 processors on an HP Proliant DL160 G5 servers operating at 3 Ghz speed. In this study, displacement time histories are simulated for two sample realization of slip field. The focal depth of the two samples are 35 km and 26 km. In Figure 10(a,b), the simulated displacement time histories at several other stations which suffered heavy damage in the Sikkim for the both the slip models shown in Figure 4. Permanent displacements can be observed in the simulated time histories (Figure 10) at several stations in the near-field region. The displacement field has been computed at a spacing of 2.6 km covering a region of dimension 350 km6 350 km around the fault. The contours of peak ground displacement in horizontal and vertical direction corresponding to slip model shown in Figure 4(a) is shown in Figure 11(a,b) along with location of landslides. It can be observed that PGD values are high at Mangan, Chungthang, Dikchu and Gangtok. The comparison between the PGD contours and location of landslides is good. It is well known that, at low frequencies, near-field ground motions will be strongly influenced by the dimensions of the fault plane, slip distribution on the rupture plane, direction of rupture propagation and orientation of the station with respect to the fault. The ground motions also follow certain radiation patterns and exhibit rupture-directivity effects. These low frequency characteristics can be observed from the above figures. Figure 12(a–c) presents the contours of permanent surface displacements in the north, east and vertical directions corresponding to the slip model shown in Figure 4(a). The near-source features such as permanent ground displacement can be observed at sites near the fault. These static displacements decrease with an increase in epicentral distance and are valid at bedrock level. Since the top soil would behave in a nonlinear fashion, the large displacements in the linear model would hint at ground failures due to the high level of strains. The field 28 S.T.G. Raghukanth et al. Figure 10. Analytical surface displacement time histories (Cm.) o-PGD, Blue line-slip model of Figure 4(a). Red line simulated from the slip model shown in Figure 4(b). Available in colour online. Estimation of ground motion 29 Figure 11. (a) Peak ground horizontal displacement (cm) in the for the slip model shown in Figure 4(a) (Focal depth ¼ 35 km). Black dots indicate location of landslides. (b) Peak ground displacement (cm) in vertical direction the for the slip model shown in Figure 4(a) (Focal depth ¼ 35 km). Black dots indicate location of landslides. Available in colour online. 30 S.T.G. Raghukanth et al. Figure 12. (a) Permanent ground displacement (cm) simulated from slip model shown in Figure 4(a) (Focal depth ¼ 35 km). East component. Black dots indicate location of landslides. (b) Permanent displacement (cm) simulated from slip model shown in Figure 4(a) (Focal depth ¼ 35 km). North component. Black dots indicate location of landslides. (c) Permanent ground displacement (cm) simulated from slip model shown in Figure 4(a) (Focal depth ¼ 35 km). Vertical component. Black dots indicate location of landslides. Available in colour online. Estimation of ground motion 31 investigations carried out after the earthquake indicated ground cracks and openings on the order of a few centimetres in the near field region which matches with the permanent displacements simulated from SPECFEM (Figure 12). 5. Summary and conclusions This article estimates the ground motion during the 18th September 2011 Sikkim earthquake (M 6.8) by analytical approaches. In absence of strong motion data, it is more scientific to use the stochastic finite fault seismological model (Boore 2009) for generating large samples of artificial ground motion time histories. With the help of such a sample where in the uncertain source, path and site parameters can be randomly varied, one can arrive at reliable ground motion relations for a given earthquake. Stochastic seismological model for simulation of ground motion depends on region specific earthquake model parameters such as stress drop, slip distribution, dip of the fault plane, focal depth, site amplification function and the quality factor. These model parameters are available for Himalayan region. This helps in fixation of the limit of uncertainties in the input parameters for Sikkim earthquake (Table 1). Accelerations at A-type rock level have been simulated in a 200km6 200km region around the epicentre. The approach of Wald and Allen (2007) has been used to obtain the V values in the epicentral region. The soil has been classified into B- type, C-type and D-type sites depending on the V and correction factors given in IBC (2003) are used to arrive at the surface level PGA. The contour map of mean PGA is shown in Figure 7. To validate the simulated ground motions, further estimates of PGA have been obtained from the MMI values. The comparisons at specific stations are shown in Table 5. It is seen that the PGA values obtained from both the approaches are consistent. The differences at some stations are attributed to local site effects and the complex topography of the region. The spectral finite element method is used to simulate the ground motion in the entire Earth medium including the Himalayan topography. The computational difficulties put limits on the high frequency end of the spectrum and hence numerical results are limited to displacement time histories only. Since the required computation for SPECFEM is large, the ground motion time histories have been obtained for two sample realizations of the slip distribution (Figure 4). Contours of a PGD and PGRD over a region of 350 km6 350 km, in and around the epicentre have been obtained for one sample realization of the slip distribution. The simulated time histories contains permanent ground residual displacement indicating ground failures in the near source region which matches with reported in the field investigations. The simulated ground motions in this study bring out the directivity effects due to fault orientation and rupture direction. Figure 12 reflects the permanent movement at the end of the earthquake shaking due to the effect of slip across the fault. These static deformations are observed in the near neighbourhood of faults and decrease sharply with distance as one moves away from the faults. The displacements will be of use in design of infrastructural facilities such as tunnels, pipelines, transmission lines that run across faults in Sikkim. The damages to structures during the Sikkim earthquake have been mainly attri- buted to poor construction practices in this region which increases the vulnerability of this region. In order to reduce the vulnerability, one requires ground motion of damaging earthquakes. The results obtained from this study can be of varied use in 32 S.T.G. Raghukanth et al. infrastructural developments in Sikkim Himalaya. For example, the obtained PGA and displacement time histories will be of use in the aseismic design of new structures and also for assessing the performance of existing structures. The contour maps of PGA, PGD and PGRD presented in Figures 7, 11 and 12 will be useful to engineers in selecting the amplitude of motions to be considered in design to achieve the objective of no collapse and preservation of life safety for events of comparable magnitude in this region. The estimated ground motion contours obtained from this study can be directly used to understand the pattern of damage and also the behaviour of structures which did not suffer appreciable damage during the main event. Since such moderate earthquakes will inevitably occur in the future, the results obtained in this study can be an important aid for emergency planning. Since near-field strong motion records are sparse in the Himalayan region, the estimated time histories in this study will be useful in understanding ground motion during a future earthquake of comparable magnitude in the Himalayan region. Acknowledgements Thanks to Prof. C.V.R. Murthy for his encouragement and help in fixing the MMI values at several towns in Sikkim state. References AKI, K. and RICHARDS, P.G., 1980, Quantitative Seismology: Theory and Methods, Vol. I. pp. 1–700. (New York: WH Freeman). ATKINSON, G.M. and BOORE, D.M., 2006, Earthquake ground-motion prediction equations for eastern North America. Bulletin of the Seismological Society of America, 96, pp. 2181– BASSIN, C., LASKE, G. and MASTERS, G., 2000, The current limits of resolution for surface wave tomography in North America. EOS Trans AGU, 81, pp. F897. BOORE, D.M., 1983, Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra. Bulletin of the Seismological Society of America, 73, pp. 1865–1894. BOORE, D.M., 2009, Comparing stochastic point-source and finite-source ground-motion simulations: SMSIM and EXSIM. Bulletin of the Seismological Society of America, 99, pp. 3202–3216. BOORE, D.M. and ATKINSON, G.M., 1987, Stochastic prediction of ground motion in eastern North America. Bulletin of the Seismological Society of America, 4, pp. 460–477. BOORE, D.M. and BOATWRIGHT, J., 1984, Average body-wave radiation coefficients. Bulletin of the Seismological Society of America, 74, pp. 1615–1621. BOORE, D.M. and JOYNER, W.B., 1997, Site amplifications for generic rock sites. Bulletin of the Seismological Society of America, 87, pp. 327–341. BOZDAG, E., CHEN, M., HJORLEIFSDOTTIR, V., KOMATITSCH, D., LIU, Q., MESCHEDE,M., MICHEA, D., PETER, D., SAVAGE, B., SCHUBERTH, B., SIEMINSKI, A., TAPE, C., TROMP,J. and ZHU, H., 2010, SPECFEM3D_GLOBE User Manual. Computational infra- structure for geodynamics, pp. 1–67. (France: Princeton University, California Institute of Technology and University of Pau/CNRS/INRIA). BRUNE, J.N., 1970, Tectonic stress and the spectra of seismic shear waves from earthquakes. Journal of Geophysical Research, 75, pp. 4997–5009. DAHLEN, F.A. and TROMP, J., 1998, Theoretical Global Seismology, pp. 1–1027. (Princeton, NJ: Princeton University Press). Estimation of ground motion 33 DE, R. and DE KAYAL, J.R., 2003, Seismotectonic model of the Sikkim Himalaya: constraint from microearthquake surveys. Bulletin of the Seismological Society of America, 933, pp. 1395–1400. Geological Survey of India (GSI), 2000, Seismotectonic Atlas of India and Its Environs, pp. 1– 87. (Kolkata: Geological Survey of India). HAZARIKA, P., RAVI KUMAR, M., SRIJAYANTHI, G., SOLOMON RAJU, P., PURNACHANDRA RAO,N. and SRINAGESH, D., 2010, Transverse tectonics in the Sikkim Himalaya: evidence from seismicity and focal-mechanism Data. Bulletin of the Seismological Society of America, 100, pp. 1816–1822. IBC, 2003, International Building Code, pp. 1–673. (Washington, DC: International Code Council, ICC). IMAN, R.L. and CONOVER, W.J., 1980, Small sample sensitivity analysis techniques for computer models, with an application to risk assessment. Communications in Statistics, A9, pp. 1749–1842. IS 1893, 2002, The Indian Standard Code, Criteria for Earthquake Resistant Design of Structures, pp. 1–40. (New Delhi: Bureau of Indian Standards). IYENGAR, R.N. and RAGHUKANTH S.T.G., 2003, Attenuation of strong ground motion and site specific seismicity in Peninsular India. In Proceedings of National Seminar on Seismic Design of Nuclear Power Plants, 21–22 February 2003, SERC, Chennai (Delhi: Allied Publishers Pvt. Ltd.), pp. 269–291. IYENGAR, R.N. and RAGHUKANTH, S.T.G., 2006, Strong ground motion estimation during the Kutch, India earthquake. Pure Applied Geophysics, 163, pp. 153–173. KOMATITSCH, D. and TROMP, J., 2001, Modeling of seismic wave propagation at the scale of the earth on a large Beowulf. In Proceedings of the ACM/IEEE supercomputing SC’2001 conference, CD-ROM. KOMATITSCH, D. and TROMP, J., 2002a, Spectral-element simulations of global seismic wave propagation–I. Validation. Geophysical Journal International, 149, pp. 390–412. KOMATITSCH, D. and TROMP, J., 2002b, Spectral-element simulations of global seismic wave propagation–II. 3-D models, oceans, rotation, and self-gravitation. Geophysical Journal International, 150, pp. 303–318. KUSTOWSKI, B., EKSTROM, G. and DZIEWONSKI, A.M., 2008, Anisotropic shear-wave velocity structure of the Earth’s mantle: a global model. Journal of Geophysical Research, 113, pp. B06306, doi:10.1029/ 2007JB005169. MAI, P.M. and BEROZA, G.C., 2002, A spatial random field model to characterize complexity in earthquake slip. Journal of Geophysical Research, 107, pp. 2308–2329, doi:10.1029/ 2001JB000588. MAI, P.M., SPUDICH, P. and BOATWRIGHT, J., 2005, Hypocenter locations in finite-source rupture models. Bulletin of the Seismological Society of America, 95, pp. 965–980. MONSALVE, G., SHEEHAN, V., SCHULTE-PELKUM, V., RAJAURE, S., PANDEY, M.R. and WU,F., 2006, Seismicity and one-dimensional velocity structure of the Himalayan collision zone: earthquakes in the upper crust and upper mantle. Journal of Geophysical Research, 111, pp. B10301, doi:10.1029/2005JB004062. MOTAZEDIAN,D.and ATKINSON, G.M., 2005, Stochastic finite-fault modeling based on a dynamic corner frequency. Bulletin of the Seismological Society of America, 95, pp. 995–1010. PARVEZ, I.A., VACCARI, F. and PANZA, G.F., 2003, Deterministic seismic hazard map of India and adjacent areas. Geophysical Journal International, 55, pp. 489–508. PATERA, A.T., 1984, A spectral element method for fluid dynamics: laminar flow in a channel expansion. Journal of Computational Physics, 54, pp. 468–488. RAGHUKANTH, S.T.G., 2008, Ground motion estimation during the Kashmir earthquake of 8th October 2005. Natural Hazards, 46, pp. 1–13. RAGHUKANTH, S.T.G. and BHANUTEJA, B., 2011, Ground motion simulation for 26th January 2001 Gujarat earthquake by spectral finite element method. Journal of Earthquake Engineering, (In Press). 34 S.T.G. Raghukanth et al. RAGHUKANTH, S.T.G. and DASH, S.K., 2009, Deterministic seismic scenarios for North East India. Journal of Seismology, 14, pp. 143–167. RAGHUKANTH, S.T.G. and SOMALA, S.N., 2009, Modeling of strong motion data in Northeastern India: Q, Stress Drop and Site Amplification. Bulletin of the Seismological Society of America, 99, pp. 705–725. RAY, S.K., 2000, Culmination zones in Eastern Himalaya. Geological Survey of India (Special Publication), 55, pp. 85–94. SARAGONI, G.R. and HART, G.C., 1974, Simulation of artificial earthquakes. Earthquake Engineering and Structural Dynamics, 2, pp. 249–267. SHINOZUKA, M. and DEODATIS, G., 1996, Simulation of multi-dimensional gaussian stochastic fields by spectral representation. Applied Mechanics Reviews – ASME, 49, pp. 29–53. SINGH, S.K., GARCIA, D., PACHECO, J.F., VALENZUELA, R., BANSAL,B.K.and DATTATRAYAM, R.S., 2004, Q of the Indian shield. Bulletin of the Seismological Society of America, 94, pp. 1564–1570. SINGH, S.K., ORDAZ, M., DATTATRAYAM, R.S. and GUPTA, H.K., 1999, A spectral analysis of the 21 May 1997, Jabalpur, India, earthquake (Mw ¼ 5.8) and estimation of ground motion from future earthquakes in the Indian shield region. Bulletin of the Seismological Society of America, 89, pp. 1620–1630. SOMERVILLE, P.G., IRIKURA, K., GRAVES, R., SAWADA, S., WALD, D., ABRAHAMSON,N., IWASAKI, Y., KAGAWA, T., SMITH, N. and KOWADA, A., 1999, Characterizing crustal earthquake slip models for the prediction of strong ground motion. Seismological Research Letters, 70, pp. 59–80. WALD, D.J. and ALLEN, T.I., 2007, Topographic slope as a proxy for seismic site conditions and amplification. Bulletin of the Seismological Society of America, 97, pp. 1379–1395.

Journal

"Geomatics, Natural Hazards and Risk"Taylor & Francis

Published: Feb 1, 2012

There are no references for this article.