American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) ISSN (Print) 2313-4410, ISSN (Online) 2313-4402 © Global Society of Scientific Research and Researchers http://asrjetsjournal.org/ Analysis of India-Nepal Border Region Earthquake Wave for Studying the Layering Information of the Crust Md. Omar Faruka*, Md. Abdullah Al Asad b , A. T. M. Saiful Islamc, Md. Saiful Islamd aAssistant Professor, Dept. of ICE, Pabna University of Science and Technology, Pabna-6600, Bangladesh b,cLecturer, Dept of APECE, Bangabandhu Sheikh Mujibur Rahman Science & Technology University, Gopalgonge-8100, Bangladesh dAssociate Professor, Dept. of ETE, Pabna University of Science and Technology, Pabna-6600, Bangladesh aEmail: fom_06@yahoo.com bEmail: asad_135uits@yahoo.com cEmail: saifulslm6@gmail.com dEmail: msislam0701@gmail.com Abstract Seismic surface wave group velocity dispersion has computed for India-Nepal border region earthquake of magnitude 5.0 occurred on 28 March 2012 of 23:40:14 UTC by graphical method. A model taking subsurface layer parameters is also constructed to compute the group velocity dispersion by modified Haskell matrix method. Group velocity dispersion by graphical method is then interpreted from model parameters. Sensitivity and the statistical errors of the model are studied and presented in this research. Interpreted crustal structure of the India-Nepal region shows that there are four major subsurface layers of thickness 4.0 km, 8.0 km, 11.0 km and 20.0 km. Keywords: Layering information; period; group velocity; earthquake data; seismic wave; model parameter. 1. Introduction There are different techniques for earthquake seismic surface wave analysis. Group velocity dispersion analysis is one of the most useful techniques. ------------------------------------------------------------------------ * Corresponding author. 290 http://asrjetsjournal.org/ mailto:fom_06@yahoo.com mailto:saifulslm6@gmail.com American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 Seismic surface wave dispersion analysis of the local earthquakes can be used to study the layering information of the crust using simple models of continental or oceanic crust. Direct and indirect modeling techniques are commonly used in determination of the earth's interior from seismic surface wave dispersion. Direct modeling determines the layering information from observed surface wave dispersion [1, 2, 3, 4]. On other hand, the most widely used indirect modeling techniques deal with trial-and-error procedures. Dispersion of seismic surface wave is computed for different model parameters to see how the computed dispersion matches with observed dispersion [5]. Ewing in [6] first introduced such model for the oceanic crust using Rayleigh wave dispersion. In this research, the group velocity dispersion has been computed and analyzed using graphical method [6] for the up-down component of the ground accelerated earthquake seismic wave of India-Nepal border region recorded at Dhaka Meteorological Department on 28 March 2012 of 23:40:14 UTC. 2. Material and Methods 2.1Earthquake Data Magnitude 5.0 earthquake occurred on 28 March 2012 at India-Nepal border region. This region is 110 km far from Dinajpur city, Bangladesh and also not so far from the Himalayan frontal arch. Table 1 lists the source parameters of the selected event. The event was recorded at Dhaka Meteorological Department seismic station located at 23.780 N and 90.380E and the station is equipped with a three component digital broad–band sensor which can record up-down, north-south and east-west components. The recorded earthquake seismic wave is shown in Figure 1. Table 1: Earthquake Source parameters Date Origin Time Location Depth (Km) Distance of Epicenter (Km) Mw 28th March 2012 23:40:14(UTC) 05:40:14 (BST) 26.093440N, 87.75130E 40 369 4.1 2.2 Methods Seismic surface wave group velocity dispersion is considered as the factor, which has relationship with structure of the crust. Group velocity from recorded earthquake wave and multilayered crustal model can be obtained respectively by graphical method and modified Haskell matrix method as explained below. • Graphical Method Graphical method is basically a technique of group velocity dispersion determination. In this method the travel times ( )t of some chosen phases along the surface wave train are measured and plotted on a graph versus the 291 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 order number ( )n of the chosen phases. Usually the travel times of the wave crests and troughs are read. The ( )tn, curve built by these points is then approximated with linear segments. The period is determined by the slope of these lines and the corresponding travel times are read from the midpoints of the segments [6]. The group velocity, gU of seismic surface wave can be obtained as: 0 500 1000 1500 2000 2500 3000 3500 4000 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x 10 5 Time in Second A m pl itu de Figure 1: Up-down ground accelerated earthquake seismic wave recorded at Dhaka Meteorological Department of 5.0M earthquake on 28 March 2012 at 23:40:14 UTC that occurred at India-Nepal border region t U g ∆ = (1) Where ∆ is epicenteral distance and t is the travel time. • Modified Haskell Matrix Method Modified Haskell matrix method for the case of 1−n homogeneous, isotopic elastic layers over a half-space matrix can be written as [7]: 121 ........ AAAAEJ mnn −=  (2) Where mA is the 44Χ Haskell matrix for the m’th layer and nE  is the half-space inversion matrix. Then the secular function (dispersion relation) can be written as: 292 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 0|),( 21122211 12 12 =−==∆ JJJJJCT (3) Hence two columns or rows of J are necessary for this result and it requires a 42Χ matrix to store the product as each layer. It has seen that Haskell matrix poses a loss of significant figures in the secular function [8, 9, 10]. In order to minimize the losses a 66Χ matrix is employed where the elements are second order sub- determinant of the Haskell matrix. The matrix can be written as: 1 6 1 − = ∑= m j j m ij m i RBR (4) Where m ijB = ij kl mA | Hence the secular function: ( ) 0, ==∆ nn BRcT (5) nB is the matrix of sub-determinants of the half-space nE  . However, B matrix shows that 1 1 1 ii BR = (6) and mm RR 43 = [Using Eqn. (4)] This phenomenon leads to define a 55Χ matrix B  rather than 66Χ matrix and the modified B  matrix can be expressed as: B11 B12 2B13 B15 B16 B21 B22 2B23 B25 B26 B  = B31 B32 (2B33-1) B35 B36 (7) B51 B52 2B53 B55 B56 B61 B62 2B63 B65 B66 and similarly m iR is thus reduced from a six dimensional vector to five. Dispersion relation (Equation. 5) can be solved numerically according to the model parameters (Vp, Vs, ρ and 293 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 thickness) in the form of group velocity versus time period plot. On other hand same plot can also be obtained from the recorded earthquake data using graphical method. Hence crustal interpretations are now possible in an indirect way by matching both dispersion relations. 3. Sensitivity of Earth Model Parameter Rayleigh wave phase velocity and group velocity (dispersion data) are the function of four parameters: S-wave velocity, P-wave velocity, density, and layer thickness [11]. Each of the parameters contributes to the dispersion curve. Using above parameters an initial earth model is constructed (Table 2). The group velocity with period is also computed and shown in Figure 2. It is observed that the variations in S-wave velocities have a dramatic effect on Rayleigh wave group velocities (Figure 2). S-wave velocities are changed by 1% in the model (Table 2), an average change of 2.81% in group velocity. After changing the S-wave velocity, the group velocity is represented by the diamond. Effects of the 1% changes in S-wave velocity are quite dramatic in comparison to similar changes in P-wave velocity and density (Figure 2).Group velocities are influenced much less by changes in density than P-wave velocity. A 1% increase in P-wave velocities (Table 2) represents an average group velocities change 0.48% and average group velocity changes 0.0067% for changing the density. This significant change in density has very subtle effect on group velocity (Figure 2).The effect of layer thicknesses on Rayleigh wave group velocities can be minimized by dividing the sub-surface into thinner layers within each unique and constant S-wave interval velocity. When the model (Table 2) defines a thickness increase of 1% the average change in group velocities is approximately 0.012% (Figure 2). According to inspection (Figure 2) it can be said that the group velocity increases with S-wave, P-wave velocity and density, but decreases with increasing thickness. The S-wave velocity is the dominant parameter influencing the changes in group velocity for this particular type model. Table1 2: Initial Earth model parameters. Layer number Vp (km/s) Vs (km/s) ρ (gm/cc) h (km) 1 5.54 3.20 2.54 5.0 2 5.63 3.25 2.57 12.0 3 5.89 3.40 2.65 15.0 Half-space 6.10 3.53 2.72 Infinite 4. Model Error Estimation The aim of the current research is to study the layering information of the crust / shallow depth hence the depth of 46.0 km is being considered here. The S-wave velocities of the layers are free to change during the inversion. Consequently, the P-wave velocities are estimated using the Vp/Vs ratio 1.732. Poisson’s ratio (σ) in each layer was assumed to be .25 and the densities (ρ) are calculated from the P-wave velocities (Vp) using the relation 0.32Vp+0.77 [12]. Starting from initial estimates, the model parameters are iteratively improved until a good fit between the theoretical and 294 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 observed dispersion curves is obtained. During the inversion, a number of criteria were adapted to calculate the goodness of fit. These criteria are the standard error of estimate (SE), mean residual (MR), average absolute residual (AR), weighted root mean square error (RMS) and the percent of signal power fit (SPF). These criteria are computed by [13]: 0 10 20 30 40 50 60 70 80 90 100 3.17 3.18 3.19 3.2 3.21 3.22 3.23 3.24 3.25 3.26 Period in Second G ro up V el oc ity in K m /S Change in Vs Initial Model Change in Vp Change in Density Change in Thickness Figure 2: Sensitivity of group velocity dispersion obtained by 1% changes in each model parameter of the initial earth model as shown in Table 1. ( ) 1 1 2 − − = ∑ = N meanobsi SE N i (8) MR = (9) ( )∑ = − = N i SE prediobsi N AR 1 1 (10) ( ) ( )1 1 2 − − = ∑ = NN mean prediobsi RMS N i (11) Where ‘obs’ is the observed group velocity at each period, ‘mean’ is the mean of the observed group velocities, N is the number of observations at each period and ‘pred’ is the predicted group velocity of the current model. Estimated errors of the models are shown in table 3. 295 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 ( ) ( ) 1001 2 1 2 1 X obs predobs SPF N i i n i ii             − −= ∑ ∑ = = (12) 5. Crustal Thickness Measurement Group velocity dispersions are estimated in this section using graphical method (Eqn.1) and Haskell modified matrix method (Eqn. 5) as discussed below: 5.1. Group Velocity Estimation from Earthquake data Group velocity is computed for the earthquake data recorded at Dhaka Meteorological Department seismic station, Bangladesh (located at 23.780 N and 90380 E) equipped with a three component digital broad–band sensor (Figure 1) and earthquake source parameters are shown in table 1. Figure 3 shows order number (n) versus travel time (t) plot and Figure 4 shows group velocity variation with time period. This dispersion relation is computed by graphical method (Eqn. 1). 0 200 400 600 800 1000 1200 0 2000 4000 6000 8000 10000 12000 14000 16000 18000 Travel Time in Second O rd er N o. (n ) Figure 3: order number versus travel time curve for India-Nepal border region earthquake data 5.2 Group Velocity Estimation from Model There are eight models (D1-D8) are considered in this work. Model based group velocity is computed using modified Haskell matrix method (Eqns. 2-7). The computed group velocity according to model parameters are shown in Figs. 5-12 also show the group velocity computed by graphical method. 296 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 0 20 40 60 80 100 120 3.29 3.295 3.3 3.305 3.31 3.315 Period in Second G ro up V el oc ity in K m /S Figure 4: Group velocity dispersion curve for India-Nepal border region earthquake data 0 10 20 30 40 50 60 70 80 90 100 3.29 3.295 3.3 3.305 3.31 3.315 Period in Second G ro up V el oc ity in K m /S Observed data Measured data Model D1 Vp (Km/s) Vs (Km/s) Density (gm/cc) Thickness (Km) 5.726 3.3065 2.602 3.0 5.802 3.350 2.626 18.0 6.148 3.550 2.737 8.0 5.767 3.330 2.615 19.0 6.283 3.628 2.780 Inf. Figure 5: group velocity dispersion obtained from India-Nepal border region earthquake data and from modeling D1. Rectangular box contained the model parameters. 0 10 20 30 40 50 60 70 80 90 100 3.29 3.295 3.3 3.305 3.31 3.315 Period in Second G ro up V el oc ity in K m /S Observed data Measured data Model D2 Vp (Km/S) Vs(Km/S) Density (gm/cc) Thickness (Km) 5.725 3.306 2.602 4.0 5.923 3.420 2.665 20.0 5.975 3.450 2.628 9.0 6.062 3.500 2.709 13.0 6.283 3.628 2.780 Inf. Figure 6: Group velocity dispersion obtained from India-Nepal border region data and from modeling D2. Rectangular box contained the model parameters. 297 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 0 10 20 30 40 50 60 70 80 90 100 3.29 3.295 3.3 3.305 3.31 3.315 Period in Second G ro up V el oc ity in K m /S Observed data Measured data Model D3 Vp (Km/S) Vs (Km/S) Density (gm/cc) Thickness (Km) 5.725 3.306 2.602 4.0 5.923 3.42 2.665 20.0 5.975 3.45 2.682 9.0 6.062 3.50 2.709 12.0 6.283 3.628 2.780 Inf. Figure 7: Group velocity dispersion obtained from India-Nepal border region earthquake data and from modeling D3. Rectangular box contained the model parameters. 0 10 20 30 40 50 60 70 80 90 100 3.29 3.295 3.3 3.305 3.31 3.315 Period in Second G ro up V el oc ity in K m /S Observed data Measured data Model D4 Vp (Km/S) Vs (Km/S) Density (gm/cc) Thickness (Km) 5.725 3.306 2.602 3.0 5.923 3.420 2.665 20.0 5.975 3.450 2.682 9.0 6.062 3.500 2.709 13.0 6.283 3.628 2.780 Inf. Figure 8: Group velocity dispersion obtained from India-Nepal border region earthquake data and from modeling D4. Rectangular box contained the model parameters. 0 10 20 30 40 50 60 70 80 90 100 3.29 3.295 3.3 3.305 3.31 3.315 Period in Second G ro up V el oc ity in K m /S Observed data Measured data Model D5 Vp (Km/S) Vs (Km/S) Density (gm/cc) Thickness (Km) 5.725 3.306 2.602 4.0 5.767 3.330 2.615 10.0 5.802 3.350 2.626 11.0 5.888 3.400 2.654 20.0 6.283 3.628 2.780 Inf. Figure 9: Group velocity dispersion obtained from India-Nepal border region earthquake data and from modeling D5. Rectangular box contained the model parameters. 298 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 0 10 20 30 40 50 60 70 80 90 100 3.29 3.295 3.3 3.305 3.31 3.315 Period in Second G ro up V el oc ity in K m /S Observed data Measured data Model D6 Vp (Km/S) Vs (Km/S) Density (gm/cc) Thickness(Km) 5.725 3.306 2.602 4.0 5.923 3.420 2.665 20.0 5.975 3.450 2.682 8.0 6.062 3.500 2.709 12.0 6.283 3.628 2.780 Inf. Figure10: Group velocity dispersion obtained from India-Nepal border region earthquake data and from modeling D6. Rectangular box contained the model parameters. 0 10 20 30 40 50 60 70 80 90 100 3.29 3.295 3.3 3.305 3.31 3.315 Period in Second G ro up V el oc ity in K m /S Observed data Measured data Model D7 Vp (Km/S) Vs (Km/S) Density (gm/cc) Thickness (Km) 5.726 3.3065 2.602 3.0 5.802 3.3500 2.626 18.0 6.148 3.5500 2.737 8.0 5.767 3.3300 2.615 15.0 6.283 3.6280 2.780 Inf. Figure 11: Group velocity dispersion obtained from India-Nepal border region earthquake data and from modeling D7. Rectangular box contained the model parameters. 0 10 20 30 40 50 60 70 80 90 100 3.29 3.295 3.3 3.305 3.31 3.315 Period in Second G ro up V el oc ity in K m /S Observed data Measured data Model D8 Vp (Km/S) Vs (Km/S) Density (gm/cc) Thickness (Km) 5.726 3.306 2.602 4.0 5.767 3.330 2.615 8.0 5.802 3.350 2.626 11.0 5.888 3.400 2.654 20.0 6.283 3.628 2.780 Inf. Figure 12: Group velocity dispersion obtained from India-Nepal border region earthquake data and from modeling D8. Rectangular box contained the model parameters. 299 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 6. Interpretation From the modeling subsurface layers are estimated I th as seen that (Figs. 5-12) group velocity obtained from earthquake data and from models have the similar characteristics as both are varying with period, and to a reasonable maximum velocity. Therefore, interpretations are made from model parameters as shown in the rectangular box in Figs. 5-12. None of the plots is found matched exactly and it should not match as the models consider only four variables, in fact there should be few more variables. Hence, statistical errors are analyzed as explained in section (4). Using Eqns. 8-12 the computed errors are shown in Table 3. Table 3: Data fit criteria Model SE MR AR RMS SPF D1 0.0173109 0.0002477 0.0143114 0.0003304 99.99997829 D2 0.0173109 0.0011323 0.0654135 0.0003305 99.99997000 D3 0.0173109 0.0009765 0.0564137 0.0003252 99.99997234 D4 0.0173109 0.0011079 0.0640002 0.0003299 99.99997057 D5 0.0173109 0.0004663 0.0269418 0.0003499 99.99997487 D6 0.0173109 0.0009186 0.0530661 0.0003240 99.99997317 D7 0.0173109 -0.0002366 0.0136731 0.0002824 99.99997831 D8 0.0173109 0.0001649 0.0095266 0.0003144 99.99998440 According to estimated statistical errors (Table 3) the model E is found more acceptable. Hence it can be said that the India-Nepal earthquake wave is indicated that there are four major subsurface layers and layer thicknesses are shown in Figure 12. 7. Conclusion There are few challenges to set up the model parameters. Most critical constraint is to consider the Poisson’s ratio of 0.25. In real cases the ratio might be different for different subsurface layers and hence the interpreted crustal structure from model might not be appropriate. However for the computational advantages Vp/Vs ratio or Poisson’s ratio were kept fixed as it has seen in many contributions to use the value of 1.732 or 0.25 respectively On other hand from the sensitivity of the models, it has shown that the thickness of the layers is a vital factor therefore, thickness setting in the model is also found to be difficult [12]. However, from the investigations it is revealed that the setting of total depth rather than individual thicknesses of the subsurface layers can provide better interpretations that are more acceptable. Hence, total depth of 48.0 km is considered in our models. Instead of above limitations interpretation made from the four models are seemed good enough with the group velocities obtained by graphical methods as shown in Figs 5-12. Group velocity dispersion from the eight models (D1-D8) (Figs. 5-12) and considering statistical error analysis (Table 3), it can be said that all the models are very nearer to an acceptable matching level though the statistical confidence level SPF should be 91.5% but our results are around 99.99998440%. Considering all errors studying in this research (Table. 3) 300 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 model D8 is seemed more acceptable of India-Nepal border region. Hence the interpreted subsurface layers of the studied this earthquake data shows that there are four major subsurface layers having respectively the thickness and density of 4.0 km, 2.602 gm/cc; 8.0 km, 2.615 gm/cc; 11.0 km, 2.626 gm/cc; 20.0 km, 2.654 gm/cc. Acknowledgements Authors acknowledge Bangladesh Meteorological Department for providing seismological data. References [1]. J. Dorman and M. Ewing, “Numerical Inversion of Seismic Surface Wave Dispersion Data and Crust- Mantle Structure in the New York-Pennsylvania Area,” Journal of Geophysical Research, vol. 67, pp. 5227-5241, Dec, 1962. [2]. T. V. McEvilly, “Central U. S. Crust-Upper Mantle Structure from Love and Rayleigh Wave Phase Velocity Inversion,” Bulletin of the Seismological Society of America, vol. 54, pp. 1997-2015, Dec, 1964. [3]. L. W. Braile and G. R. Keller, “Fine Structure of the Crust Inferred from Linear Inversion of Rayleigh- Wave Dispersion,” Bulletin of the Seismological Society of America, vol. 65, pp. 71-83, Feb, 1975. [4]. S. Bloch, A. L. Hales and M. Landisman, “Velocities in the Crust and Upper Mantle of Southern Africa from Multi-Mode Surface-Wave Dispersion,” Bulletin of the Seismological Society of America, vol. 59, pp. 1599-1630, Aug, 1969. [5]. J. Dorman, M. Ewing and J. Oliver, “Study of Shear-Velocity Distribution in the Upper Mantle by Mantle Rayleigh Waves,” Bulletin of the Seismological Society of America, vol. 50, pp. 87-115, Jan, 1960. [6]. M. Ewing and F. Press, “Crustal Structure and Surface Wave Dispersion; Part II, Solomon Islands Earthquake of July 29, 1950,” Bulletin of the Seismological Society of America, vol. 42, pp. 315-325, Oct, 1952. [7]. T. H. Watson, “A Note on Fast Computation of Rayleigh Wave Dispersion in the Multilayered Elastic Half-Space,” Bulletin of the Seismological Society of America, vol. 60, pp. 161-166, Feb, 1970. [8]. L. Knopoff, “A Matrix Method for Elastic Wave Problems,” Bulletin of the Seismological Society of America, vol. 54, pp. 431-438, Feb, 1964. [9]. J. Dunkin, “Computation of Modal Solutions in Layered, Elastic Media at High Frequencies,” Bulletin of the Seismological Society of America, vol. 55, pp. 335-358, Apr, 1965. [10]. E. N. Thrower, “The Computation of the Dispersion of Elastic Waves in Layered Media,” Journal of Sound and Vibration, vol. 2, pp. 210-226, July, 1965. [11]. J. Xia, R. D. Miller and C. B. Park, “Estimation of Near-Surface Shear-Wave Velocity by Inversion of Rayleigh Waves,” Geophysics, vol. 64, pp. 691–700, May-June, 1999. [12]. C. J. Ammon, G. E. Randall and G. Zandt, “On the Non-uniqueness of Receiver Function Inversions,” Journal of Geophysical Research, vol. 95, pp. 15303-15318, Sep, 1990. [13]. K. M. A. Elenean, K. S. Aldamegh, H. M. Zharan, and H. M. Hussein, “Regional Waveform 301 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 15, No 1, pp 290-302 Inversion of 2004 February 11 and 2007 February 09 Dead Sea Earthquakes,” Geophysical Journal International, vol. 176, pp.185–199, Jan, 2009. 302