Relation of Mean Time Interval with Magnitude for Earthquakes in Northeast India and Its Surrounding Region

Based on active tectonics, structural orientation, focal mechanism of earthquakes and seismicity level, Northeast India and its surrounding region has been divided into five seismogenic zones namely, the Eastern Himalaya as Zone 1, the Eastern Himalaya Syntaxis as Zone 2, the Shillong plateau and Mikir hill as Zone 3, the Naga-Disang as Zone 4 and the Eastern Boundary as Zone 5. In the present study we proposed relation between earthquake mean time interval and the magnitude for earthquakes in these five seismogenic zones. We model the variation of the mean time interval with four different distribution models viz. the Poisson’s, the Weibull’s, the log-Weibull and the lognormal distributions. The mean time interval is found to be following the log-normal model the most amongst the four models tested in this study. The coefficients of the relation between logarithm of mean time interval and the magnitude as estimated from the complete catalog selected by two time windows exhibit mild differences limited to ± one standard deviation uncertainty in the coefficient and decreasing with an increase in the catalog time window length.


Introduction
Seismicity analysis assesses the behavior of earthquake occurrences with respect to previous events.This is one of the important components in seismic hazard assessment.In the time independent probabilistic seismic hazard assessment the basic assumption is that seismicity distribution follows the Poisson's distribution in space and time according to which the next event from a seismic source is independent of the previous event i.e. it is memory-less.Kagan and Jackson [1976], Kagan and Jackson [1991], Knopoff et al. [1996] and others observed clustering in time for large earthquakes, and therefore, proposed other types of distributions.Moreover for the seismic hazard and risk studies conducted for purposes relating to a shorter time interval (e.g. insurance renewal periods), temporal variations in the estimated hazard level may have significant effects [Musson et al. 2002].In the time dependent probabilistic seismic hazard analysis the activity rate of a region is considered to be dependent on time delay from the last earthquake and time interval between two events of the magnitude in the same source zone.
There are various methods to estimate the time interval between two events viz.Paleoseismicity study [McCalpin 2009], from the slip rate on the fault assuming characteristic earthquake behavior and relatively constant recurrence intervals [Wallace 1970] and historical catalog analysis.The slip rates on major faults in northeast India and its surrounding region have been investigated by various workers [viz.Angelier and Baruah 2009, Kundu and Gahalaut 2013, Gahalaut et al. 2013] though the precision is not sufficient enough to deduce the slip rate for each potential fault in the region.Northeast India has been the source of a number of great earthquakes viz.Assam Earthquake of magnitude M w 8.7 in 1950, Shillong Earthquake of M w 8.1 in 1897, Srimangal Earthquake of M w 7.6 in 1918 and Dhubri Earthquake of M w 7.1 in 1930.The return period for the great earthquakes like the 1897 Shillong Earthquake of M w 8.1 is expected to be of the order of thousand years [Bilham and England 2001].However, before the present the earthquake records are close to complete only for earthquakes in the most recent 200 years [Bilham 2004].Therefore, the relationship between mean time interval and magnitude has become important for computation of time dependent activity rate in the region.
A number of different methods have been proposed by various researchers to include time dependence in activity computation.A simple model proposed by Gere and Shah [1984] is that the longer the lapse time since the last event occurred, the sooner will the next event be nucleating.However, Davis et al. [1989] suggested that the reverse might also be true.Ward and Goes [1993] and Goes and Ward [1994] numerically proved that in the case of Weibull's distribution the correlation between the time lapse after the last event and the probability of occurrence of the next event may be positive or negative depending on the exponent of the distribution.Sornette and Knopoff [1997] analyzed this for other distribution models with memory and found that for any distribution that falls off at a faster rate than an exponential function at large time intervals, the correlation between the time lapse after the last event and the probability of occurrence of the next event is negative while for any distribution model that falls at a slower rate than an exponential one at large time intervals, the correlation between the time lapse after the last event and the probability of occurrence of the next event is positive.The negative correlation between the time lapse after the last event and the probability of occurrence of the next event has resemblance with seismic gap hypothesis proposed by Fedotov [1965].In the lognormal distribution the mean time interval is the most probable time lapse since the last event to the next.If the earthquake has not occurred near the mean time interval, it becomes increasingly likely that the occurrence time will be in the tail of the lognormal distribution [Musson et al. 2002].Various other distribution models have been tested by several workers viz.Weibull [Abaimov et. al. 2008], Brownian Passage time [Ellsworth 1999], lognormal [Musson 2002], Poisson, and gamma [Utsu 1984].Time independent seismicity analysis of Northeast India and its surrounding region has been performed by Thingabijam et al. [2008].Yadav et al. [2010] calculated the cumulative probability for earthquakes of magnitude greater than or equal to M w 7.0 using 20 events of magnitude ranging from M w 7.0 to M w 8.6 from northeast India and the surrounding region.

Seismotectonism and Source Zonation
Delineation of seismogenic source zones requires several parameters viz.homogeneous and complete historical catalog, neo-tectonic fault mapping and the geology of the region.Dutta [1964] and Gupta et al. [1986] divided Northeast India and the surrounding region into four seismogenic source zones, namely, the Eastern Himalayan Thrust zones (in present study seismogenic Zone 1), the Eastern Himalaya Syntaxis zone (in present study seismogenic Zone 2), the Shillong plateau (in present study seismogenic Zone 3) and the Arakon-Yoma subduction (in present study seismogenic Zone 4 and seismogenic Zone 5).Yadav et. al. [2009] divided the region into four source zones based on active tectonics, focal mechanism of earthquakes and the seismicity level.Due to non-availability of detailed map of neo-tectonism of the region Yadav et. al. [2011] divided the region into four seismic zones with overlapping boundaries.Studies from GPS data along the Arakan-Yoma subduction zone indicate significant differences in slip distribution from the northern to the southern part of the fault system [Vigny et al. 2003, Maurin et al. 2010].Also the structural orientation in this zone varies from NNE-SSW in the north to N-S in the south [Angelier and Baruah 2009].Therefore, the Arakan-Yoma subduction zone has been split into two seismogenic zones in this study, viz. the Naga-Disang (seismogenic Zone 4) and the Eastern Boundary (seismogenic Zone 5).The entire Northeast India and the surrounding region has been divided into five seismogenic zones namely, the Eastern Himalaya as seismogenic Zone 1, the Eastern Himalaya Syntaxis as seismogenic Zone 2, the Shillong Plateau and Mikir hill as seismogenic Zone 3, the Naga-Disang as seismogenic Zone 4 and the Eastern Boundary as seismogenic Zone 5 as depicted in Figure 1.The main tectonic features in the Eastern Himalayan zone are the Main Boundary Thrust (MBT) and the Main Central Thrust (MCT).In the Eastern Himalaya Syntaxis, located in the eastern end of the Himalaya, the East-West directed Eastern Himalaya structure intersects with the NNE-SSW Naga-Disang Thrust [Menahem et al. 1974].The dominant focal mechanism in this zone is of thrust fault type (Figure 2).The seismogenic Zone 3 comprises the Shillong Plateau and Mikir Hill.Shillong Plateau, the source of the 1897 earthquake of M w 8.1 has a complex tectonics.Bilham and England [2001] suggested a pop-up tectonics for the Shillong Plateau between the Oldham fault in the north and the Dauki fault in the south.The Mikir Hill is separated from the Shillong plateau in the East by the Kopili fault.The earthquakes in this zone have focal mechanism varying from thrust to strike-slip types.Kayal [1987Kayal [ , 1996] ] and Dian et al. [1984] suggested subduction tectonics and dragging of the dipping Indian lithosphere below the Indo-Burma ranges.Chen and Molnar [1990] obtained focal mechanism solutions for 10 earthquakes which show pure thrust to a mixture of reverse in the seismogenic Zone 4. The earthquake focal mechanism has dominantly strike-slip faulting type in the Eastern Boundary zone [Angelier and Baruah 2009] as shown in Figure 2.

Methodologies
The methodology for the present analysis comprises of two steps viz. 1) pre-processing of the earthquake catalog and 2) computation of mean time interval and selection of optimal model for the variation in mean time interval with magnitude.We have used cumulative probability to compare the variation in mean time interval with magnitude in four different model types viz.Poisson's, Weibull's, Log-weibull and Lognormal distributions for the selection of the optimum distribution model.The Kolmogorov Smirnov (K-S) and Chi-Square methods have been used to evaluate the model fit [Lilliefors 1967].

Pre-Processing of Catalog
The earthquake catalog for the study region has been selected from South Asia earthquakes catalog prepared by Nath et al. [2010].The selected catalog has been updated up to 31 st December 2012.Earthqua- kes reported from various global earthquake catalogs viz.United States Geological Survey (USGS), Global Centroid Moment Tensor (GCMT), International Seismological Center (ISC), India Meteorological Department (IMD) and National Geophysical Research Institute (NGRI) served as data sources for the catalog.The earthquake catalog has been made homogeneous in moment magnitude (M w ) by using empirical relations between moment magnitude (M w ) with the body wave magnitude (m b ), the surface wave magnitude (M s ) and the local magnitude (M L ) proposed by Nath et al. [2010] using orthogonal regression on the earthquake events included from the same study region.
Declustering of the catalog is the selection of main-shocks from a cluster of foreshocks, mainshocks and aftershocks.The method of either Gardner & Knopoff [1974] or Reasenberg [1985] is mostly used because of their application simplicity [Stiphout et. al., 2012].Gardner and Knopoff [1974] defined the method to identify aftershocks and foreshocks in the earthquake catalog using inter event distance in time and space.But this method ignored the secondary and higher order aftershocks [Stiphout et al. 2012].The method proposed by Musson [2000] is similar to the method given by Reasenberg [1985].The main difference is that the Reasenberg's method considers the first event of a sequence as the mainshock, and a subsequent larger earthquake becomes a ''larger mainshock'' [Reasenberg and Jones 1989] while the Musson's method considers the largest event in a sequence to be the mainshock; in case two events of equal magnitude occur, the first event is the mainshock [Musson 2000].In the present analysis we applied declustering algorithm developed by Musson [2000] with a fixed time window of 80 days which has been found to be optimal for the study region.This value was determined through inspection of  several clusters of seismicity in the catalog to establish an appropriate value through direct observation.
The declustered catalogue comprises the area between 17°-33° N and 86°-100° E. The study region comprises of all the five seismogenic zones and a total of 27092 events of which 8549 events are of M w ≥ 5.0 and 963 events are of M w ≥ 6.0.In the declustered catalog the number of events with magnitude M w ≥ 5.0 and M w ≥ 6.0 are 3630 and 689 respectively while the total events are 6209.The declustured events within the five seismogenic zones along with faults have been depicted in Figure 2. The number of independent events in the time completeness windows A and B within the five seismogenic zones is provided in Table 1.The Frequency Magnitude Distribution (FMD) plot for the five seismogenic zones in Northeast India and its surrounding region is shown in Figure 3.
We used the method given by Zuniga and Wyss [1995] to infer the minimum magnitude of complete-  plots A and B for each seismogenic zone represent the result for the catalogs starting from 2012 onwards in the descending order of the year and the other starting from 2005 onwards in the descending order of the year, respectively.The completeness time window length for the catalog starting from 2012 and 2005 onwards in the descending order has been listed in Table 2 for all the five seismogenic zones.

Computation of Mean Time Interval
Mean Time Interval (MTI) for magnitude M is the average time interval between two earthquakes of magnitude greater than or equal to M. M varies from M w 5.0 to a maximum magnitude for which there are at least two earthquakes in the complete catalog.Two complete catalogs have been given by completeness time windows A and B selected using the Stepp's method as explained above.The catalogs within the completeness time windows A and B are the complete catalogs with latest years 2012 and 2005.MTI is varying with magnitude and exhibits positive correlation with the increasing magnitude.

Selection of Model for MTI Variation
MTI computed from the complete catalog has been fitted with four distribution models viz. the Poisson's, the Weibull's, the log-weibull and the lognormal models.We applied the method of maximum likelihood [Harrish and Stocker 1998] to infer the model parameters of the four distribution functions as depicted in Table 3. Thereafter the variation in mean time interval has been fitted with the four models mentioned above and given in the Appendix.
The method of cumulative probability has been used to compare the fitting with four distribution functions as depicted in Figure 5 for all the five seismogenic zones.The subplots A and B in Figure 5, show MTI computed with the catalog windows A and B, respectively, as discussed earlier.As shown in Figure 5 the cumulative probability with the lognormal model is found to be in good agreement with cumulative probability calculated from empirical cumulative distribution function (ECDF), discussed in details in the result and discussion sections to follow.We also tested for model fit by applying appropriate method like Kolmogorov-Smirnov (K-S) technique which is the most appropriate approach for testing for the continuous distribution [Miller, 1956, Marsaglia et al. 2003] and is, therefore, considered for lognormal, Weibull's and the log-Weibull models while for the Poisson's model the Chi-Square test [Cochran 1952] has been applied which is preferred for testing the discrete distributions.The lognormal model is found to be the optimum model for checking the variation of MIT with magnitude.Therefore, the MTI variations in all the five seismogenic zones have been fitted with the lognormal model as given in equation ( 1) using the method of least square minimization; Where, MTI is the mean time interval for magnitude greater than or equal to M, and α and β are the model coefficients which remain constant in each seismogenic zone.

Results
Based on active tectonics, structural orientation, focal mechanism of earthquakes and seismicity level, the Northeast India and its surrounding region have been divided into five seismogenic zones modified from Yadav et al. [2011].In each seismogenic zone the minimum magnitude of completeness is depicted in Figure 3 as the frequency magnitude distribution (FMD) plot.The time completeness test proposed by Stepp [1973] has been applied to the entire declustered catalog with a fixed magnitude range varying from M w 5.0 to M w 7.0 and the results presented in Figure 4. We computed MTI for magnitude M varying from M w 5.0 to M w 7.0, a magnitude that occurred at least twice in the complete catalog.The catalog with the completeness time windows A and B has been considered complete for the magnitude range M ≥ M w 4.0 to M ≤ M w 7.0.The MTI has been fitted in four different distribution models viz. the Poisson's, the Weibull's, the log-Weibull and the lognormal models.Figure 5 depicts the comparison of the cumulative probability computed from the four different distribution models with that from ECPD.The real data represent the cumulative probability computed using ECPD.
In order to judge the effect of different catalog time windows on the model fit, we performed the model fit test for MTI computed from catalog with time windows A and B. The subplot A and the subplot B in Figure 5 depict the comparison of the cumulative probability computed from ECPD with four different models for the two catalog time windows A and B. The cumulative probability estimated by lognormal model is found to be in good agreement with the real data for both the time windows in all the five seismogenic zones.The MTI variation with magnitude has been fitted in equation (1) using the least square minimization method as depicted in Figure 6 and Figure 7 for windows A and B, respectively.
The coefficients of the relation given by equation (1) for the two complete catalogs in all the five seismogenic zones have been listed in Table 4.The difference in the model coefficients β between the two time windows A and B in the seismogenic Zones 1 and 5 is found to be 0.24 and 0.20, respectively, while in the seismogenic Zone 3 the minimum difference is found to be of the order of 0.01.Moreover for the seismogenic Zones 1 and 5 the standard deviation in the β value is 0.08 and 0.05, respectively, for the time window A.

Discussion
The FMD plot for the earthquakes exhibits that the completeness magnitude for the earthquake catalog in the five seismogenic zones varies from M w 5.0 in the seismogenic Zones 1, 2 and 3 to M w 5.2 in the seismogenic Zone 5 as depicted in Figure 3.In order to make equal completeness magnitude in all the five seismogrnic zones we rounded it off to M w 5.0 for all seismogenic zones.In each of the five seismogenic zones the largest magnitude M for which there is at least two events with magnitude ≥ M is smaller than M w 7.0.Therefore, the magnitude range from M w 5.0 to  M w 7.0 has been considered for the catalogs time completeness test using the Stepp's method.The time completeness test with two windowed catalogs -one up to 2012 and the other up to 2005 exhibits that the former has a larger completeness time length as compared to the later as depicted in subplot A and subplot B respectively in Figure 4. Thus time span of the complete catalog time window A is larger than the time span of the complete catalog time window B as given in Table 1.
In each of the five seismogenic zones the higher magnitude earthquake has larger value of MTI.As depicted in Figure 5, for small values of MTI the real data have higher deviations from the model derived values, while for large MTI the deviation is lesser, which can be attributed to the saturation effect in the distribution model.The real data represent the cumulative probability computed from MTI using empirical cumulative distribution function (ECDF).As shown in Figure 5, however, among the four models the cumulative probability computed with the Weibull's model and the lognormal model is found to be comparable with the real data.Moreover among the two comparable models, for the lognormal model the difference is smaller except in the seismogenic Zone 5 for the time window B. Nevertheless, for the time window B in the seismogenic Zone 5 the real data is not monotonic thereby exhibiting a misfit with both the comparable models i.e. the lognormal and the Weibull's models.For both the time windows A and B cumulative probability computed from ECDF and that computed by the lognormal model are in good agreement.
Thus lognormal model yields cumulative probability comparable with that computed from ECDF for both the catalog windows A and B in all the five seismogenic zones, therefore, it is considered to be the best fit model for testing the MTI variation with magnitude.The linear fit of logarithm of MTI with magnitude in the five seismogenic zones with two catalog windows has been depicted in Figures 6 and 7 respectively.As given in Table 4, in all the seismogenic zones the estimated coefficients from the catalog time window A has smaller value of standard deviation as compared to that for the catalog time window B, which can be attributed to the time window length.With an increase in the time span of the complete catalog, the standard deviation in the estimated coefficients is found to be decreasing.Thus the coefficients of equation ( 1) derived from the complete catalog time window A is preferred for the estimation of MTI at various magnitudes in each seismogenic zone of Northeast India and its surrounding region.
In the time dependent probabilistic seismic hazard assessment the activity rate for earthquake of a magnitude has a functional dependence on time lapse since the last event occurred, the time interval (inter-arrival time) and the uncertainty in time interval for the magnitude.The time lapse since the last event of an earthquake magnitude can be computed from an earthquake catalog.The computation of MTI for a magnitude from the catalog need time span of complete catalog equal to at least twice the MTI for the earthquake of that magnitude.Northeast India and its surrounding region have experienced earthquakes of great magnitudes having return period of the order of thousand years [Bilham and England 2001].Therefore, the relation between MTI and the magnitude becomes an important proposition in the estimation of time interval at various magnitudes in the seismogenic source zones of the study region.The coefficients have been estimated for five different seismogenic zones that comprises number of faults.

Conclusions
In the present study a linear relationship between logarithm of mean time interval and magnitude has been established for five seismogenic zones in Northeast India and its surrounding region.Lognormal model is found to be the best fit model for defining the variation of mean time interval with magnitude after testing four distribution models.It is observed that model coefficients vary within the catalog time windows.However the difference observed are limited to ± one standard deviation associated with the coefficients.The larger the time span for the complete catalog time window, the smaller will be the standard deviation in the model coefficients.In the absence of fault specific slip rates for all the major active faults in Northeast India and its surrounding region the relation between mean time interval and magnitude becomes an integral part for the computation of time-dependent activity rate in the region.The coefficients have been estimated for five different seismogenic zones that comprises of numerous faults.The present relation estimates the mean time interval between earthquake of a magnitude considering all earthquakes greater than or equal to the magnitude within the seismogenic zone.

Data and Resources
In present analysis we have used earthquake catalog for a period of 1902-2012.The online available sources of the earthquake catalog used in this study are: National Geophysical Data Center (NGDC), United State Geological Survey (USGS), Global Centroid Moment Tensor (GCMT), International Seismological Center (ISC) and India Meteorological Department (IMD).Link: http://www.earthqhaz.net/sacat/

Appendix
Weibull's and Log-weibull Distribution Functions: Weibull's distribution function is a two parametric distribution model.The probability density function for Weibull's distribution is given by, (2) where, t is the variable which in our case is MTI and γ and α are the scale and shape parameters respectively.
The cumulative distribution function for the Weibull's model is given as, for the log-Weibull distribution function the cumulative distribution function is given by, Poisson's Model: Poisson's distribution is a discrete probability distribution and is given as, where, t = 1, 2, 3…inter arrival time in month and λ is the model parameter.Cumulative distribution function for Poisson's model is given by, Lognormal Model: The probability density function for the lognormal distribution is given as, (7) where, t is a variable and µ and σ are the model parameters.The cumulative distribution function of the lognormal model is given as, (8) where, ϕ is the cumulative distribution function of the normal distribution.In the present analysis we applied maximum likelihood method to estimate the model parameters that has been given in Table 3 for all four different distribution models.

Figure 3 .
Figure 3. Frequency Magnitude Distribution (FMD) plots for the five seismogenic zones in Northeast India and its surrounding region.

Figure 4 .
Figure 4. Plots between standard deviation and mean of occurrence rate with catalog time interval to extract complete catalog time window length in the five seismogenic zones, (A) catalog time window starting from 2012 and (B) catalog time window starting from 2005.

Figure 5 .
Figure 5.Comparison of the observed Mean Time Interval (MTI) with four different distribution models in the five seismogenic zones, subplot (A) for the complete catalog with the time window starting from 2012, and subplot (B) for the complete catalog with time window starting from 2005.

Figure 6 .
Figure 6.Least square fitting of logarithms of MTI with magnitude for the catalog time window starting from 2012.Two subplots (A) and (B) are used to avoid intersection of lines.

Figure 7 .
Figure 7. Least square fitting of logarithms of MTI with magnitude for the catalog time window starting from 2005.Two subplots (A) and (B) are used to avoid intersection of lines.

Table 1 .
Number of events after declustering within windows A and B in five seismogenic zones.

Table 2 .
Complete catalog time length for windows A and B for the catalog starting from 2012 and 2005 respectively in the five seismogenic zones.
[1973].The Stepp's method determines the fraction of the catalog time span in which the mean rate of occurrence is stable for a fixed magnitude range (magnitude needs to be greater than Mc) using the Stepp's plot.The time upto which the slope in Stepp's plot is constant gives completeness time span of the catalog.In the pre-ferent time windowed catalogs on the model fit the two completeness time windows have been inferred from the two catalogs -one up to 2012 and the other up to 2005.The plot of standard deviation of mean occurrence rate with time length for time completeness for all the five seismogenic zones are shown in Figure4.The

Table 3 .
Estimated coefficients of the four distribution models by two complete catalog time windows A and B using the method of maximum likelihood.α* represent γ, µ and λ of the Weibull, lognormal and Poisson' model β* is for α, σ l for Weibull and lognormal models.The parameters for Weibull and log-Weibull are same.

Table 4 .
Coefficients of the relation between log MTI and magnitude given by equation (1), in the five seismogenic zones for two complete catalog time windows A and B given in Table2.