Completo_DEF.qxd 1301 ANNALS OF GEOPHYSICS, VOL. 47, N. 4, August 2004 Key words aftershocks – Bhuj earthquake – Omori’s Law – GBA array – Lg magnitude 1. Introduction The destructive Bhuj earthquake of Mw 7.7 on January 26, 2001, which ravaged through the Kachchh and adjoining regions of India, caused a death toll of 30 000 and loss of prop- erty worth more than 10 million US dollars (Gupta et al., 2001). The India Metrological Department (IMD) gave the epicentre of this earthquake as 23.40N and 70.28E with focal depth of 25 km. The magnitude of this earth- quake was similar to that of an earlier earth- quake, which had hit the northern part of Kachchh, India, in 1819 (Oldham, 1926). The epicentre of this Bhuj earthquake and the area of occurrence of aftershocks are shown in fig. 1. The aftershocks following this major earth- quake were recorded at GBA, India, at dis- tance of around 1300 km from the epicentre. The layout of seismometers at GBA station is also given in the in-set of fig. 1. In this paper the temporal behavior of the aftershocks in the initial 144 h is analysed. 2. Theory The stress energy following a main shock is slowly released as a sequence of aftershocks. According to Omori (1894), the rate of after- Mailing address: Dr. G. Jayachandran Nair, Electro- nics and Instrumentation Group, Seismology Division, Bhabha Atomic Research Centre, Trombay, 400085 Mum- bai, India; e-mail: nairgj@apsara.barc.ernet.in The oscillatory behaviour of the aftershocks rate of the 2001 Bhuj earthquake, India: observation and interpretation G. Jayachandran Nair, M.M. Abdul Razak, A.G. Venkatesh Prasad and Erulankatil Unnikrishnan Electronics and Instrumentation Group, Seismology Division, Bhabha Atomic Research Centre, Trombay, Mumbai, India Abstract A damaging earthquake of Mw 7.7, which struck the Bhuj region of India on January 26, 2001, was followed by a large number of aftershocks. The aftershock data available at Gauribidanur Seismic Array Station (GBA), India, till 869 h following the main shock were compiled. The plot of the aftershocks rate with time was found to be os- cillatory decay. There was a sharp decrease of the aftershocks rate in the initial 144 h from the main shock and this paper presents the analysis of the temporal characteristics of aftershock activity during this period. A statistical best fit for the rate of aftershocks is performed using the generalised Omori’s law and the exponential decay law. The statistical errors for the exponential fit are found to be lower than that of the generalised Omori’s fit. The super- imposed oscillations present in the aftershock activity are extracted by differencing the observed aftershock activ- ity from the statistical fits. The frequencies of these oscillations are found to be 0.062 h–1, 0.078 h–1, 0.102 h–1, 0.118 h–1, 0.141 h–1, 0.164 h–1, 0.233 h–1 and 0.476 h–1. Some of the plausible causes for this kind of oscillations present in the aftershock activity are also discussed in this paper. 1302 G. Jayachandran Nair, M.M. Abdul Razak, A.G. Venkatesh Prasad and Erulankatil Unnikrishnan shocks decays with time as 1/t, where t is the time from the main shock. A generalization of Omori’s Law was proposed by Utsu (1961), in which the rate of aftershocks decays with time as 1/(t + t0)n, where t0 is a constant and n is the exponent. From this generalized Omori’s rela- tion, the number of aftershocks, N, occurring per unit time is proportional to 1/tn, when t >> t0. The aftershock activity of the California earthquakes gave the value of n lying between 0.5 and 1.5 (Reasenberg and Jones, 1989), which is also observed in most of the after- shock activity in the rest of the world. A sim- ilar power law for the rate of aftershocks is al- so obtained due to viscoelastic creep (Schaff et al., 1998) and stress transfer mechanism (Stein, 1999). In the above mechanism there is an implicit assumption that the coulomb stress change is proportional to the sum of normal stress (∆σd) and pore pressure change (∆P). The coulomb stress change (∆σf) is de- fined as (Stein, 1999) ∆σf = ∆τ + µ(∆σd + ∆P) (2.1) where ∆τ is the shear stress on a fault and µ is the frictional co-efficient. The rate of after- shocks N, for such a process at a time t, N = N0⋅(t + t0)–n (2.2) is obtained from the solution of the following differential equation, for the decay of the stress s (Shaw, 1993), ds/dt = –A⋅s (n+1)/ n. (2.3) The nature of the solution of eq. (2.3) takes the form as eq. (2.2), under the assumption that s and N are linearly related. As the exponent n is empirically observed to be less than 1.5, an ex- ponential relation for the rate of decay of after- shocks will not arise as a solution of eq. (2.3). The scaling and translation of time ordinate is permitted in equation (2.2). Thus N0 ⋅t0–n denotes Fig. 1. The geology and neotectonics of Kachchh region, India, affected by the Bhuj earthquake. The star indicates the IMD epicenter. The circle around the star shows the area of occurrence of aftershocks. The in- set box shows the lay-out of GBA seismic array, at a distance of ~1300 km from the epicentre, India. 1303 The oscillatory behaviour of the aftershocks rate of the 2001 Bhuj earthquake, India: observation and interpretation the rate of aftershocks at time t = 0. The general- ized Omori’s Law is here in after referred as the hyperbolic law in this text. 3. Instrument and ground motion data The ground motion data of the aftershock sequence of this Bhuj earthquake were record- ed at GBA, (fig. 1) India, in digital format and used for the analysis. The recording instru- ments were vertical component seismometers of natural frequency 1 Hz. The sampling fre- quency was 20 samples per second. There were a large number of overlapping events in the first three hours following the main shock and it was very difficult to separate them. Hence the aftershocks which originated later than three hours following the main shock were on- ly considered for the analysis. The most promi- nent arrival from these aftershocks was the Lg phase and it was used to quantify the Lg mag- nitude. All the aftershocks used in this analysis had Lg magnitude varying from 5.9 to 3.5. The minimum detectable Lg magnitude at GBA from Bhuj was ~ 3.2. A few examples of ob- served digital waveforms of aftershocks corre- sponding to USGS mb 5.5 and 4.7 are shown in fig. 2. The nature of the rate of aftershocks recorded till 869 h from the main shock is shown in fig. 3 as trace ‘a’. A total of 320 af- tershocks were recorded during this period. From fig. 3, it is observed that there is a sharp decrease of rate of aftershocks in the initial 144 h, an increase in activity from 144 h to 300 h and a slowly decreasing activity from 300 h onwards. In order to confirm that the recorded aftershocks originated from the Bhuj region, the difference between the time of arrival of S and P phases and Lg- and P-phases were checked and found to match the expected times from the Bhuj region. Digital array tun- ing was also done for some of the selected events to further confirm whether their origin coordinates matched the expected area. 4. Statistical analysis The plot of rate of aftershocks against time for 869 h, following the main shock, is shown in fig. 3 as trace ‘a’. Traces ‘b’ and ‘c’ of fig. 3 indicate the exponential fit, 10.95.e–t/36.53, and the hyperbolic fit, 39.72.t–0.76, respectively. The Fig. 2. Typical aftershock digital waveforms of Bhuj earthquake, 2001, observed at GBA for USGS mb 5.5 and 4.7. The date and time of arrival at GBA are shown in the figure. 1304 G. Jayachandran Nair, M.M. Abdul Razak, A.G. Venkatesh Prasad and Erulankatil Unnikrishnan Fig. 3. Trace ‘a’ shows the plot of rate of after- shocks of Bhuj earthquake, 2001, against time, ob- served at GBA, for 869 hour following the main shock, The trace ‘b’ shows the exponentially fitted curve, 10.95e–t/36.53, for the trace ‘a’. The standard de- viation (SD) and χ 2 are 0.37 and 0.16 respectively for the trace ‘b’. Trace ‘c’ shows the hyperbolically fitted curve, 39.72.t–0.76, for the trace ‘a’. The SD and χ 2 are 0.58 and 0.37 respectively for the trace ‘c’. All parameters and errors are rounded off to two decimal places. Function Constants Standard Deviation (SD) χ 2 N0.exp (–t/Tc) N0 = 11.05 Tc = 35.91 0.73 0.54 N0.(t) –n N0 = 31.93 n = 0.65 1.18 1.44 Table I. Fitted parameters and statistical errors (rounded off to two decimal places) in the initial 144 h of rate of decay of aftershocks of the Bhuj earth- quake, 2001. Fig. 4. Trace ‘a’ shows the plot of rate of aftershocks of Bhuj earthquake, 2001, against time in the initial 144 hour. Trace ‘b’ shows the exponentially fitted curve, 11.05.e – t/35.91, for the trace ‘a’. The SD and χ 2 are 0.73 and 0.54 respectively for the trace ‘b’. Trace ‘c’ shows the hyperbolically fitted curve, 31.93.t – 0.65 for the trace ‘a’. The SD and χ 2 are 1.18 and 1.44 respectively for the trace ‘c’. All parameters and errors are rounded off to two decimal places. standard deviation (SD) and Levenberg-Mar- quardt χ 2 values (Levenberg, 1944; Marquardt 1963; More, 1977) are 0.37 and 0.16 respec- tively for the exponential fit and 0.58 and 0.37 respectively for the hyperbolic fit. The Leven- berg-Marquardt χ 2 value for the exponential fit is denoted as χe2 while it is denoted as χh2 for the hyperbolic fit. For 5χe2 deviation value of the exponential fit, a confidence level of 95% is obtained and the corresponding hyperbolic fit gives 94% confidence level. The rate of de- cay of aftershocks in the initial 144 h is smoothed by a three point adjacent average method and is shown in fig. 4 as trace ‘a’. This curve is analytically fitted by an exponential law of the form N = N0.exp (–t/Tc) and hyper- bolic law of the form N = N0⋅(t)–n as shown in fig. 4 as traces ‘b’ e ‘c’ respectively. For the ex- ponential fit the amplitude N0, the time con- stant Tc, SD and χe2 values are found to be 11.05, 35.91 h, 0.73 and 0.54 respectively. For the hyperbolic fit, the amplitude N0, the expo- nent n, the SD and χh2 and are found to be 31.93, 0.65, 1.18 and 1.44 respectively. For 4 χe2 deviation value of the exponential fit, a confidence level of 95% is obtained while the hyperbolic fit gives 89% confidence level. The exponential fit gives greater confidence level than the hyperbolic fit. The results of these sta- tistical fittings are given in table I. The plot of logarithm of the aftershocks rate against time (here-in-after called log-linear plot) during this 1305 The oscillatory behaviour of the aftershocks rate of the 2001 Bhuj earthquake, India: observation and interpretation period is shown in fig. 5 as trace ‘p’ and its lin- ear fit is given as trace ‘r’ in the same figure. The plot of logarithm of rate of aftershocks against logarithm of time (here in after called log-log plot) is shown in fig. 5 as trace ‘q’ and its linear regression fit is given as trace ‘s’ in the same figure. The statistical errors of these log-linear and log-log fits are given in table II. The linearity of the log-linear curve is more prominent than that of the log-log curve as found from the statistical errors given in table II. It can be concluded from the trends of these curves shown in fig. 5 that the instantaneous slope of trace ‘q’ is time dependent. The large SD for the exponential and hyperbolic fits can be attributed to the superimposed oscillations present in the aftershock activity. Such oscilla- tory behavior of aftershocks over hyperbolic trend was observed for cumulative Benioff strain by Sornette and Sammis (1995). 5. Extraction of oscillatory periods The superimposed oscillations present in the rate of decay of aftershocks are analytically extracted by subtracting the exponential and hyperbolic fits, traces ‘b’ e ‘c’ of fig. 4, from the observed decay curve, trace ‘a’ of fig. 4. These subtracted curves are smoothed separate- ly by a three point adjacent average method to avoid sharp discontinuities. These smoothed curves are filtered through a high pass filter with a cut-off frequency at 0.05 h–1. The filtered curves for the exponential and hyperbolic fits are shown in figs. 6 and 8 respectively. The fil- tered curves, figs. 6 and 8, are fourier trans- formed to obtain the amplitude and phase spec- trum which are shown in figs. 7 and 9 respec- tively. The amplitude and phase spectra for the Function Standard Deviation (SD) χ 2 Log-linear 0.39 0.15 Log-log 0.47 0.22 Table II. Statistical errors (rounded off to two dec- imal places) for linear fits of log-linear and log-log plots of rate of decay of aftershocks (initial 144 h) of the Bhuj earthquake, 2001. Fig. 5. The log-linear plot of observed rate of after- shocks of Bhuj earthquake 2001 is shown as trace ‘p’ and trace ‘r’ represents the linear fit, –0.02.t + 2.16, for the trace ‘p’. The SD and χ 2 for are 0.39 and 0.15 respectively for the trace ‘r’. The log-log plot of the rate of aftershocks is shown as trace ‘q’ and trace ‘s’ represents the linear fit, –1.07.ln(t)+4.84, for the trace ‘q’. The SD and χ 2 are 0.47 and 0.22 respec- tively for the trace ‘s’. All parameters and errors are rounded off to two decimal places. Fig. 6. The plot of the difference between the ob- served rate of aftershocks, trace ‘a’ of fig. 4, and its exponential fit, trace ‘b’ of fig. 4, after filtering through a High Pass Filter with a cut off frequency at 0.05 h–1. This curve shows the oscillatory part present in the rate of aftershocks. 1306 G. Jayachandran Nair, M.M. Abdul Razak, A.G. Venkatesh Prasad and Erulankatil Unnikrishnan Fig. 7. The fourier amplitude and phase spectra of the oscillatory part shown in fig. 6. The amplitude spectrum shows prominent peaks at 0.062 h–1, 0.078 h–1, 0.102 h–1, 0.118 h–1, 0.141 h–1, 0.164 h–1, 0.182 h–1, 0.213 h–1, 0.233 h–1 and 0.476 h–1 (rounded off to three decimal places). h–1 and 0.476 h–1. In order to find the nature of aftershocks it is necessary to examine the phase spectra of figs. 7 and 9, whether these oscilla- tions correspond to random time series or de- terministic series. Runs test (Bradley, 1968) is performed to check whether a given sequence of data is random or not. The result of this test is called as z score or test statistic value. If z score is greater than 1.96, then the data under consideration are not random at a confidence level of 95%. This Runs test is explained in Ap- pendix. Using this test, the z score of the phase spectra of figs. 7 and 9 are found to be – 7.487 and – 9.635 respectively. Since the z score  for the phase spectra are greater than 1.96, the phase spectra of figs. 7 and 9 are not random at a confidence level of 95%. Hence the oscilla- tions present in the aftershock activity are not random. 6. Discussion and conclusions Generally the fourier spectra obtained from random signals will show either the amplitude and/or phase of oscillations is random. In most of the geophysical signals, the fourier phase spec- trum is observed to be random. In order to give an Fig. 8. The plot of the difference between the ob- served rate of aftershocks, trace ‘a’ of fig. 4, and its hyperbolic fit, trace ‘c’ of fig. 4, after filtering through a High Pass Filter with a cut off frequency at 0.05 h–1. This curve shows the oscillatory part pres- ent in the rate of aftershocks. Fig. 9. The fourier amplitude and phase spectra of oscillatory part shown in fig. 8. The amplitude spec- trum shows prominent peaks at 0.062 h–1, 0.078 h–1, 0.102 h–1, 0.118 h–1, 0.141 h–1, 0.164 h–1, 0.185 h–1, 0.217 h–1, 0.233 h–1 and 0.476 h–1 (rounded off to three decimal places). exponential fit are given in fig. 7. Similarly the amplitude and phase spectra for the hyperbolic fit are shown in fig. 9. From the amplitude spectra of figs. 7 and 9 the periods of prominent common spectral peaks are picked and the fre- quencies are found to be 0.062 h–1, 0.078 h–1, 0.102 h–1, 0.118 h–1, 0.141 h–1, 0.164 h–1, 0.233 1307 The oscillatory behaviour of the aftershocks rate of the 2001 Bhuj earthquake, India: observation and interpretation example, a synthetic signal is generated, fig. 10 trace ‘a’, for a deterministic case. The fourier am- plitude and phase spectrum of this deterministic signal are shown in fig. 10 as traces ‘b’ and ‘c’ re- spectively. Using the fourier amplitude spectrum of fig. 10 trace ‘b’, another synthetic signal is generated using random phases and is shown in fig. 11, trace ‘a’. The fourier amplitude and phase spectra of this signal, fig. 11 trace ‘a’, are given as traces ‘b’ and ‘c’ respectively in the same fig- ure. Runs Test is performed on the data of phase spectra of figs. 10 and 11 to check for random- ness. The z scores (test statistic values) of the fig- ures are found to be –3.34 and 0.54 respectively. As the z score for the phase spectrum of fig. 10 is greater than 1.96, this phase spectrum is not random at 95% confidence level while that of fig. 11 is random. Though the fourier amplitude spec- tra of figs. 10 and 11 are the same, their phase spectra are different since they have different z scores. The phase spectrum of fig. 10 is deter- ministic, while the phase spectrum of fig. 11 is random. Thus by observing either the fourier am- plitude or the phase spectrum one can ascertain whether a given signal under consideration is de- terministic in time or random. In this analysis of the aftershock activity of Bhuj earthquake, the fourier phase spectra of figs. 7 and 9 do not show random characteristics. Hence it is concluded that the oscillations present in the aftershock activity do not correspond to random time series. From the analysis of the initial 144 h of af- tershock activity, the exponential decay fit gives lower statistical errors than that of the generalized Omori’s fit. The superimposed oscillations present in the aftershock activity during this period are extracted by differenc- ing the observed rate of decay and the statisti- cal fits. The prominent common frequencies of these oscillations are found to be 0.062 h–1, 0.078 h–1, 0.102 h–1, 0.118 h–1, 0.141 h–1, 0.164 h–1, 0.233 h–1 and 0.476 h–1. It is also shown, using synthetic signals, that the oscillations present in the aftershock activity do not corre- spond to random time series. A few plausible causes which can generate this kind of oscil- latory decay of aftershocks are (a) the elastic or visco-elastic oscillations of the crustal blocks following the main shock and/or (b) gravitational perturbations, tidal waves, etc., Fig. 10. Trace ‘a’ shows a synthetic signal obtained by superposition of two decaying sinusoidal curves. Traces ‘b’ and ‘c’ show the fourier amplitude and phase spectra respectively for the trace ‘a’. a b c Fig. 11. The synthetic signal generated from the fourier amplitude spectrum of fig. 10, trace ‘b’, and random phases is shown as trace ‘a’. Traces ‘b’ and ‘c’ show the fourier amplitude and phase spectra re- spectively for the trace ‘a’. Compare the phase spec- tra shown in figs. 10 and 11. a b c causing modulation in the energy release through the aftershocks. In order to relate the above causes to the os- cillatory behavior of aftershock activity of this Bhuj earthquake, a correlation should exist be- tween the rate of aftershocks and measurements from tidal gauges, tilt meters and gravity me- ters. However, during the occurrence of after- shock sequence of this Bhuj earthquake no such measurements were available. Acknowledgements The authors express their sincere thanks to the staff at GBA for their efficient operation of the GBA seismic array during the aftershock se- quence of Bhuj earthquake and the excellent quality of digital data provided by them. The authors express their sincere gratitude to Shri B.V. Dinesh for his valuable assistance in ob- taining high quality of data. 1308 G. Jayachandran Nair, M.M. Abdul Razak, A.G. Venkatesh Prasad and Erulankatil Unnikrishnan Appendix. Runs test. Runs test (Bradley, 1968) is performed to check whether a given sequence of data is random or not. According to this test, the average value of the data is computed. Each member of the data whose value is above this average is assigned as +1, and those that are below are assigned as –1. The total number of +1 is taken as n1 and that of –1 is taken as n2. The number of consecutive changes of +1 to –1 and vice versa is taken as the length of the run R. Then compute n n n n1 1 2 1 2)= + +a ^ h (A.1) n n n n n n n n n n2 2 1 2 1 2 1 2 1 2 1 2 2 1 2) ) ) ) ) )= - - + + -v ^ ^ ^`h h hj (A.2) .z R= -a v^ h (Α.3) This z value is called the z score or test statistic value. It is compared with the Cumulative Standard Normal Distribution curve given by / /expZ P Z z u du1 2 2 z 2) )G= = -z r 3- #_ ` b _i i k i . (A.4) A z score with absolute value greater than 1.96 indicates non-randomness at the confidence level of 95%. REFERENCES BRADELY, J. (1968): Distribution-Free Statistical Tests, ch. 12 (on line: http://www.itl.nist.gov/div898/handbook/ eda/section3/eda35d.htm). GUPTA, H.K., N.P. RAO, B.K. RASTOGI and D. SARKAR (2001): The deadliest intraplate earthquake, Science, 291, 2101-2102. LEVENBERG, K. (1944): A Method for the solution of certain problems in Least Squares, Quart. Appl. Math., 2, 164-168. MARQUARDT, D. (1963): An Algorithm for Least-Squares Estimation of Non-linear parameters, SIAM J. Appl. Math., 11, 431-441. MORE, J.J. (1977): The Levenberg-Marquardt Algorithm: implementation and theory, in Numerical Analysis, ed- ited by G.A. WASTON (Springer Verleg), Lecture Notes in Mathematics, 630, 108-116. OLDHAM, R.D. (1926): The cutch (kachh) earthquake of 16th June 1819 with a revision of great earthquake of 12th June 1897, Mem. Geol. Surv. India, 46, 71-146. OMORI, F. (1894): Investigation of Aftershocks, Rep. Earth- quake Inv. Comm., 2103-39. REASENBERG, P.A. and L.M. JONES (1989): Earthquake hazard after a mainshock in California, Science, 243, 1173-1176. SCHAFF, D.P., G.C. BEROZA and B.E. SHAW (1998): Post- seismic response of repeating aftershocks, Geophys. Res. Lett., l25 (24), 4549-4552. SHAW, B.E. (1993): Generalised Omorìs Laws for after- 1309 The oscillatory behaviour of the aftershocks rate of the 2001 Bhuj earthquake, India: observation and interpretation shocks and foreshocks from simple dynamics, Geo- phys. Res. Lett., 20, 907-910. SORNETTE, D. and C.G. SAMMIS (1995): Complex critical exponents from renormalization group theory of earth- quakes: implications of earthquake predication, J. Phys. I, 5, 607-619. STIEIN, R.S. (1999): The role of stress transfer in earth- quakes occurrence, Nature, 402, 605-609. UTSU, T. (1961): A statistical study on the occurrence of af- tershocks, Geophys. Mag., 30, 521-605. (received September 19, 2003; accepted March 26, 2004)