The oscillatory behaviour of the aftershocks rate of the 2001 Bhuj earthquake, India: observation and interpretation

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 oscillatory 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 superimposed oscillations present in the aftershock activity are extracted by differencing the observed aftershock activity from the statistical fits. The frequencies of these oscillations are found to be 0.062 h, 0.078 h, 0.102 h, 0.118 h, 0.141 h, 0.164 h, 0.233 h and 0.476 h. Some of the plausible causes for this kind of oscillations present in the aftershock activity are also discussed in this paper.


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 property 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 earthquake was similar to that of an earlier earthquake, 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 earthquake were recorded at GBA, India, at distance 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.

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, Electronics and Instrumentation Group, Seismology Division, Bhabha Atomic Research Centre, Trombay, 400085 Mumbai, India; e-mail: nairgj@apsara.barc.ernet.inshocks 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 relation, the number of aftershocks, N, occurring per unit time is proportional to 1/t n , when t >> t 0 .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 aftershock activity in the rest of the world.A similar power law for the rate of aftershocks is also 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 defined as (Stein, 1999) where ∆τ is the shear stress on a fault and µ is the frictional co-efficient.The rate of aftershocks N, for such a process at a time t, is obtained from the solution of the following differential equation, for the decay of the stress s (Shaw, 1993), 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 exponential relation for the rate of decay of aftershocks will not arise as a solution of eq.(2.3).The scaling and translation of time ordinate is permitted in equation (2.2).Thus N 0 ⋅t 0 -n denotes the rate of aftershocks at time t = 0.The generalized Omori's Law is here in after referred as the hyperbolic law in this text.

Instrument and ground motion data
The ground motion data of the aftershock sequence of this Bhuj earthquake were recorded at GBA, (fig. 1) India, in digital format and used for the analysis.The recording instruments were vertical component seismometers of natural frequency 1 Hz.The sampling frequency 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 only considered for the analysis.The most prominent arrival from these aftershocks was the Lg phase and it was used to quantify the Lg magnitude.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 observed 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 aftershocks 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 tuning was also done for some of the selected events to further confirm whether their origin coordinates matched the expected area.

Statistical analysis
The plot of rate of aftershocks against time for 869 h, following the main shock, is shown in fig. 3     standard deviation (SD) and Levenberg-Marquardt χ 2 values (Levenberg, 1944;Marquardt 1963;More, 1977) are 0.37 and 0.16 respectively for the exponential fit and 0.58 and 0.37 respectively for the hyperbolic fit.The Levenberg-Marquardt χ 2 value for the exponential fit is denoted as χe 2 while it is denoted as χh 2 for the hyperbolic fit.For 5χe 2 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 decay 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 hyperbolic law of the form N = N0⋅(t) -n as shown in fig. 4 as traces 'b' e 'c' respectively.For the exponential fit the amplitude N 0, the time constant Tc, SD and χe 2 values are found to be 11.05, 35.91 h, 0.73 and 0.54 respectively.For the hyperbolic fit, the amplitude N0, the exponent n, the SD and χh 2 and are found to be 31.93,0.65, 1.18 and 1.44 respectively.For 4 χe 2 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 statistical 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 period is shown in fig. 5 as trace 'p' and its linear 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 oscillatory behavior of aftershocks over hyperbolic trend was observed for cumulative Benioff strain by Sornette and Sammis (1995).

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    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 oscillations correspond to random time series or deterministic 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 Appendix.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 oscillations present in the aftershock activity are not random.

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 spectrum is observed to be random.In order to give an   8.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.185 h -1 , 0.217 h -1 , 0.233 h -1 and 0.476 h -1 (rounded off to three decimal places).Thus by observing either the fourier amplitude or the phase spectrum one can ascertain whether a given signal under consideration is deterministic 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 aftershock 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 differencing the observed rate of decay and the statistical 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 correspond to random time series.A few plausible causes which can generate this kind of oscillatory 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., In order to relate the above causes to the oscillatory behavior of aftershock activity of this Bhuj earthquake, a correlation should exist between the rate of aftershocks and measurements from tidal gauges, tilt meters and gravity meters.However, during the occurrence of aftershock sequence of this Bhuj earthquake no such measurements were available.

Fig. 1 .
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 inset box shows the lay-out of GBA seismic array, at a distance of ~1300 km from the epicentre, India.

Fig. 2 .
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.

Fig. 3 .
Fig. 3. Trace 'a' shows the plot of rate of aftershocks of Bhuj earthquake, 2001, against time, observed 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 deviation (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.

Fig. 4 .
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.et/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.

Function
, from the observed decay curve, trace 'a' of fig. 4.These subtracted curves are smoothed separately 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 filtered curves, figs.6 and 8, are fourier transformed to obtain the amplitude and phase spectrum which are shown in figs.7 and 9 respectively.The amplitude and phase spectra for the Statistical errors (rounded off to two decimal 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 .
Fig.5.The log-linear plot of observed rate of aftershocks 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 respectively for the trace 's'.All parameters and errors are rounded off to two decimal places.

Fig. 6 .
Fig. 6.The plot of the difference between the observed 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.

Fig. 8 .
Fig. 8.The plot of the difference between the observed 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 present in the rate of aftershocks.
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 frequencies 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 example, a synthetic signal is generated, fig. 10 trace 'a', for a deterministic case.The fourier amplitude and phase spectrum of this deterministic signal are shown in fig. 10 as traces 'b' and 'c' respectively.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 figure.Runs Test is performed on the data of phase spectra of figs. 10 and 11 to check for randomness.The z scores (test statistic values) of the figures 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 spectra 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 deterministic, while the phase spectrum of fig.11 is random.

Fig. 10 .Fig. 11 .
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'.