Role of multifractal analysis in understanding the preparation zone for large size earthquake in the North-Western Himalaya region

Seismicity has power law in space, time and magnitude distributions and same is expressed by the fractal dimensionD, Omori’s exponentp andb-value. The spatiotemporal patterns of epicenters have heterogeneous characteristics. As the crust gets self-organised into critical state, the spatio-temporal clustering of epicenters emerges to heterogeneous nature of seismicity. To understand the heterogeneous characteristics of seismicity in a region, multifractal studies hold promise to characterise the dynamics of region. Multifractal study is done on seismicity data of the NorthWestern Himalaya region which mainly involve seismogenic region of 1905 Kangra great earthquake in the North-Western Himalaya region. The seismicity data obtained from USGS catalogue for time period 1973–2009 has been analysed for the region which includes the October 2005 MuzafrabadKashmir earthquake ( Mw = 7.6). Significant changes have been observed in generalised dimension Dq , Dq spectra and b-value. The significant temporal changes in generalised dimensionDq , b-value andDq − q spectra prior to occurrence of Muzaffrabad-Kashmir earthquake relates to distribution of epicenters in the region. The decrease in generalised dimension andb-value observed in our study show the relationship with the clustering of seismicity as is expected in self-organised criticality behaviour of earthquake occurrences. Such study may become important in understanding the preparation zone of large and great size earthquake in various tectonic regions. Correspondence to: S. S. Teotia (teotiass@rediffmail.com)


Introduction
Fracture exhibits a fractal structure over a wide range of fracture scale, i.e. from the scales of micro fracture to mega faults (Aki, 1981;Allègre et al., 1982;Brown and Scholz, 1985;Turcotte, 1986a, b;Scholz and Aviles, 1986;Hirata, 1987).Recent studies have shown that many natural phenomena such as the spatial distribution of earthquakes (Telesca et al., 2003a(Telesca et al., , b, 2004)), fluid turbulence are heterogeneous fractals (Mandelbrot, 1989).For heterogeneous fractals, a unique fractal dimension is not sufficient to characterise them and it differs by the method used to calculate it.Such fractal are called multifractals and they are characterised by generalised D q or the f (α) spectrum (Hentschel and Procraccia, 1983;Halsey et al., 1986).A wide variety of heterogeneous phenomena such as diffusion -limited aggregation, dendritic crystallization, dielectric breakdown, various fingering and river-flow, are multifractals and multifractal analysis must be used to characterise those complex phenomena (Stanley and Meakin, 1988).If the spatiotemporal distribution of earthquakes is multifractal, the fractal property of earthquakes in different areas and different time spans must be compared using whole spectrum of D q or f (α) (Hirabayashi et al., 1992).Several studies have been made to investigate the temporal variation of heterogeneity in seismicity using multifractal analysis in various seismic regions (Hirata and Imoto, 1991;Hirabayashi et al., 1992;Li et al., 1994;Dimitriu et al., 2000;Teotia et al., 1997;Telesca et al., 2005;Teotia and Kumar, 2007).Recently it has been recognised that the parameters which describe seismicity in a region show a spatial and temporal evolution which may be associated with the process of generation of large size/great size earthquakes i.e. evolving of fracture from micro fracture to mega fault (Telesca et al., 2003a, b).The change in seismicity pattern of earthquakes is reflected in generalised Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.dimension D q and D q spectra of the seismicity.Therefore, the study of temporal variation D q and D q spectra may be used to study the changes in the seismicity structure before the occurrence of large earthquake, which may prove to be useful in understanding the preparation zone of large to great size earthquakes.Accordingly (Hirata and Imoto, 1991) have performed multifractal analysis of micro earthquake data of Kanto region using correlation integral method.Hirabayashi et al. (1992) have done multifractal analysis of seismicity in regions of Japan, California and Greece using various methods such as fixed mass and fixed radius methods.Li et al. (1994) have performed the multifractal analysis of spatial distribution of earthquakes of Tanshan region (M l > 1.8) using extended Grassberger-Procraccia method of dimension D q estimation.Teotia et al. (1997), Teotia and Kumar (2007), Teotia (2000), Sunmonu et al. (2001) have studied the multifractal characteristics of seismicity of the Himalaya region and Sumatra region (m b ≥ 4.5).The 8 October 2005, Muzaffrabad-Kashmir earthquake (M w = 7.6) killed more than 80 000 people and was the deadliest in the history of the Indian subcontinent.The earthquake occurred on a rupture plane 75 km long and 35 km wide with strike of 331 • and a dip angle of 29 • .The epicenter of the event was located north of Muzaffrabad (34.493 • N and 73.629 • E) in Kashmir Himalaya (USGS1 ).The fractal nature of earthquake occurrence of the Northwest Himalayan region is also reported by Roy and Mandal (2009).However this study deals with the multifractal analysis of seismicity of region which resulted in Muzaffrabad-Kashmir earthquake of (M w = 7.6) to understand the intrinsic nature of seismicity controlling multifractal characteristics i.e. generalised dimension D q and D q spectra in the region.

Methodology
For the analysis of seismicity of the North-Western Himalaya, the earthquake data set of USGS for the grid (36 • N, 72 • E), (30 • N, 72 • E), (28 • N, 85 • E), (35 • N, 85 • E) was analysed.The tectonic map along with grid is shown in Fig. 1.The most of the seismic activity is controlled by two lineaments i.e.Main Central Thrust (MCT) and Main Boundary Fault (MBF).The time window for analysis was chosen to be 1973-2009.To gain an unbiased and homogeneous data set, we restricted the USGS data by setting a lower limit of magnitude for period 1973-2009.The magnitude of completeness is shown in the frequency magnitude plot and it is found to be complete for magnitude threshold in m b ≥ 4.5 (see Fig. 2a).The distribution of earthquakes occurring in the selected region for time period  is shown in Fig. 2b.
The available methods for calculating D q are the box counting method, the fixed-radius method and the fixed-mass method (Greenside et al., 1982;Grassberger et al., 1988;Baddi and Broggi, 1988).These methods work well provided the numbers of data points are very large.The extended Grassberger and Procaccia method (Grassberger and Procaccia, 1983) was used to calculate the generalised dimension D q and D q spectrum.This method has been applied by many authors on earthquake data of different regions for estimation of generalised dimension D q even for small data points.It is described as: logC q (r) = D q log(r)(r → 0) (1) where r is the scaling radius, N is the total number of data points within a search region in a certain time interval (also called the sample volume): X i is the epicentral location (given in latitude and longitude) of the i-th event, X j is the epicenter (given in latitude and longitude) of the j -th event, C q (r) is the q-th integral and H (.) is the Heaviside step function.
For estimating D q spectra as a function of time, a time series of earthquake epicenters has to be formed and divided into subseries (subsets).Let set {X i ,M i }, i = 1, M be a complete set of earthquakes occurring in time period analysed, and M i the magnitude of an earthquake occurring at time t i .Thus the earthquake constitutes a time series of N elements.The time series consists of 902 events in the region.We consider this time series as the original data set for multifractal analysis of this region.In this study the original data set for the North-Western Himalaya is divided into 27 subsets (i.e.S1-S27).Each subset in the region consists of 100 events with an overlap 70 events.The shift of 30 events has been used in the analysis.The same number of events in each subset will avoid any non stationarity effect in data set.There is no limitation of number of events in the formation of subsets and shift used to move from one subset to another subset.For studying the spatio-temporal variation of fractal dimensions, i.e. generalised dimensions (D q ) and the number of events i.e. 100 is large enough to provide the reliable estimate of fractal dimensions.The correlation integral is calculated using Eq. ( 2) for the epicentral distribution X i of the subset.The distance r between subset is calculated by using spherical triangle (Bullen and Bolt, 1985;Hirata, 1989).For -2 0 2 Fig. 3. (a) logC q (r) vs. logr for q = 2, -2 for subset S1; (b) logC q (r) vs. logr for q = 2, -2 for subset S10; (c) K vs. logr for subset S1; (d) K vs. logr for subset S10.
epicentral distribution having a fractal structure, the following power law relationship is obtained in the scaling region C q (r) ∼ r D q .
An appropriate scaling region has been estimated before the computation of generalised dimension D q .For this we have used Li et al. (1994) method to determine the scaling region in the western Himalaya.Three point curves for q = −2.0,0.0, 2.0 estimated in selected region for two different subsets i.e. subsets 1 and 10 are shown in Fig. 3a  and b.There may be more linear segments in each curve.It is known that slopes of adjacent points in the graph will be constant when these points lie on a linear segment.The graph of logr vs. K for two subsets is shown in Fig. 3c and d.K is the slope of adjacent points of the graph of logr vs. logC q (r).The scaling region has been selected for the range in r for which the variation in K slope is minimum for all three values of q.The scaling region has been found to be consistent irrespective of subsets.Based on this method K have been found to be minimum for range of logr (from -2.0 to 0.7) i.e. first seven points in logC q (r) vs. logr plot for all three values of q.The generalised dimensions for all 27 subsets have been estimated by least square fit from first seven points (Table 1).The errors shown in Table 1 of fractal dimensions variation in best fit line in different subsets will be reflected in least square errors.
The b-value in the Gutenberg-Richter relation (logN (m > M) = a − bM, where M is the magnitude and N (m > M) is the number of events of magnitude larger than M) was estimated for {M i} 100+30j i=1+30j of 27 subsets generated from the original data set by changing j from 0 to 26.Magnitude distribution obeys the Gutenberg-Richter relation, the b-value may be estimated by the maximum likelihood method, which gives better estimation of the b-value than the least square method (Utsu, 1965;Aki, 1965).When the number of events exceeds about 50, the maximum likelihood method provides us a stable estimation of b-value (Utsu, 1965).The number of events in each subset, N = 100 is also large enough to provide reliable estimates of b-value using maximum likelihood method (Utsu, 1965).The uncertainty is estimated by (Aki, 1965) which is given by σ b = b/ √ N, where N is number of earthquakes in the subset.The 95% confidence error in the b-values is 0.08 in our estimated b-values of 27 subsets.

Results and discussion
Observation of seismicity suggests a relationship between the distribution of earthquake magnitudes and the distribution of earthquake epicenters (Hirata, 1989;Henderson et al., 1994).The fractal dimension D 2 and D −2 for all 27 subsets are shown in Table 1. Figure 4 shows the spatio-temporal variation of generalised fractal dimension D q (q = 2,−2) along with b-values in all 27 subsets analysed in the study.It is evident that there is a significant increase in clustering prior to occurrence of Muzaffrabad-Kashmir earthquake.The same is evident in D q − q spectra.The flat D q − q spectra in subset 21 shows the convergence of epicenters to a point i.e. confinement of seismicity in the zone of preparation of large size Muzaffrabad-Kashmir earthquake (M w = 7.6) leading to very small values of generalised dimension.Subset 21 includes very few epicenters of aftershock events of Muzaffrabad-Kashmir earthquake of 8 October 2005.The steep slope in D q spectra also gives better indication for preparation zone for large size earthquake in the Himalayan region which is considered to be store house of elastic energy and have potential to generate moderate, large size and great size earthquakes.The same has been observed for the North-Western Himalaya region which shows the steep slope in D q spectra from subset S15 (1995.544-2001.759)to subset S20 (2003.527-2004.773)prior to occurrence of large size (M w = 7.6) Muzaffrabad-Kashmir earthquake which belong to subset S21 (see Fig. 5).Steep slope in D q −q spectra relates to clustering in the zone of preparation or nucleation.Therefore, in the regions having zone of preparation for large size earthquake events M > 7, the completeness of earthquake catalogue for m b ≥ 4.5 may give substantial indication in D q − q spectra as is found in this study and also in other studies of the Himalayan and Sumatra region (Teotia et al., 1997;Teotia and Kumar, 2007;Sunmonu et al., 2001).
The b-value varies from 0.94 ± .09 to 1.8 ± .18 in the subsets (S1-S27) for the time period 1973-2009 in the region.The consistent and significant decreases in D q and b-value have been observed prior to occurrence of Muzaffrabad-Kashmir earthquake.We note that the changes in b-values as well as D 2 starts from subset S15 (1995.544-2001.759)with lowest value in subset S21 (2004.773-2004.773).The Muzaffrabad-Kashmir large size earthquake lies in subset 21.The decrease in correlation dimension Dc is also reported prior to Muzaffrabad earthquake (Roy and Mandal, 2009) however the completeness of data is not carried out in their study which used merged catalogue of microseismic data of Wadia Institute of the Himalaya Geology and USGS catalogue.Therefore, the decrease in correlation dimension in study by Roy and Mandal (2009) may not entirely be associated with intrinsic nature of seismicity in the region.However, in our study the completeness of USGS catalogue is ascertained, therefore, the decrease in generalised dimension D q and b-value may be completely associated with the intrinsic nature of seismicity of the northwestern Himalaya.
The consistent decrease in spatial dimension (D s ) may correspond to introduction of clustering of seismicity in the region shows the region preparedness for the occurrence of major size events (De Rubeis et al., 1993).De Rubeis et al. (1993) have also found the temporal evolution for three seismic zones of Italy by means of the correlation integral fractal method and the decrease in fractal dimension prior of occurrence of major events also supports the results in our study.The decrease in the fractal dimension with the evolution of rock fracture was also reported by Hirata et al. (1987).These studies and our study for the Himalayan region show similar behaviour in spatio-temporal variation of fractal dimension as the seismicity in the region evolves from extended distribution of epicenters to clustered distribution and again back to distributed seismicity.This is evident from the decrease in D q from subset 15 to subset 21 and increase in D q from subset 21 to subset 27 (see Fig. 4).De Rubeis et al. (1993) have also found the same pattern of distributed Nonlin.Processes Geophys., 18,[111][112][113][114][115][116][117][118]2011 www.nonlin-processes-geophys.net/18/111/2011/ seismicity and increase in fractal dimension D s from the minimum value during the major sequence, indicating the tendency of seismic activity to spread in space.Therefore, the spatio-temporal evolution of seismicity can be best understood in the framework of self-organised criticality.The concept of Self-Organised Criticality (SOC) can be applied to distribution of seismicity (Bak andTang, 1988, 1989).The definition of Self-Organised criticality is that a natural system is in a marginally stable state when perturbed from this state, will evolve back to the state of marginal stability.
A simple Cellular-automata model illustrates self-organised criticality.The results in this study supports the evolution of epicenters distribution from a random spatial distribution to an organised fractal structure as is given in modified SOC model (Ito and Matsuzaki, 1990).Thus the application of nonlinear geophysics (NG) methods is important for extreme phenomena and new hazard assessment techniques (Rundle et al., 2003).

Conclusions
The seismicity data of studied region show multifractal behaviour.Multifractal analysis of seismicity data holds promise in understanding the zone of preparation which may result in damaging earthquakes.The spatio-temporal variation of D q and D q −q spectra provides some lessons to learn from nature.This study supports the use of spatio-temporal variation of D q at small q (i.e.q = −2, 2) to detect the seismicity change in regional area.The steep slope in D q − q spectra is detected before the occurrence of Muzaffrabad-Kashmir earthquake (M w = 7.6) of 8 October 2005.However, the steep slope of D q spectra may not necessarily be due to increase in D q for negative value of q and decrease in D q for positive values of q when compared from one subset to consecutive subset in case seismicity is converging to very localise distribution, i.e. from randomised extended distribution to self-organised clustered distribution of epicenters.The changes in D q are consistent with the tectonics of the region which evolves from plane filling of epicenters to line filling of epicenters, i.e. large size/great size earthquake occurrence along Main Central Thrust (MCT) in the Himalaya region.The same has been observed in this study which shows significant decrease in generalised dimension, i.e.D 2 and D −2 as seismicity evolves from plane filling distribution to line filling distribution.Our study supports the spatio-temporal variation of b-value, D q and D q spectra for assessing the zone of preparation for different size earthquakes.The completeness for lower magnitude threshold by ensuring better monitoring and location capabilities in the Himalaya region may help in further constraining the results and will provide better understanding of zone of preparation for major, large and great size earthquake.

Fig. 4 .
Fig. 4. The temporal variation of generalised dimension D q (q = 2, −2) along with the temporal variation of b-value.The D q (q = 2, −2) and b-values are shown with error bars.

Fig. 5 .
Fig.5.The seismicity of subsets S15, S20 and S21.The Muzaffrabad-Kashmir earthquake of October 2005 lies in subset S21.The corresponding D q − q spectra of these subsets are also shown.

Table 1 .
Showing the subsets description along with time periods and fractal dimension.