Significance of fractal correlation dimension and seismic b-value variation due to 15th July 2009, New Zealand earthquake of Mw 7.8

New Zealand earthquake that occurred on 15 th July 2009 (Mw 7.8) was analysed using fractal correlation dimension (Dc) and seismic b-value. We have analysed the earthquakes catalog of thirty-five years with a magnitude (mb ≥3.7), in order to observe a crucial information in terms of Dc value fluctuation for the event. The event is preceded by fall and anomalous change in Dc value in the year 2007 about two years prior to the mainshock. A sudden decrease in Dc value with highly clustered events is observed before the mainshock. The l ow value of Dc is an indicator of clustering and it shows how intermediate size events correlate with one another in the preparation process of this event. Here the low Dc value may be the indicator for high stress developer along the fault to produce large size earthquake. Moreover, we also observed abnormal fluctuation in b-value from 2003. The fractal clustering and scaling of earthquakes are indicated by b-value change prior to strong earthquake as a harbinger of stress correlation in various scales. The event is also marked for that occurred in the periphery of the positive Coulomb stress development, as obtained from three low Dc time windows' events. The drop in Dc value is not a single observation prior to this large event, but such pattern is also seen for other strong events in the study zone. One such well identified strong event is Mw 7.2 (2003) along with low Dc value prior to the event. Thus, stress correlation measured along with these indirect statistical tools gives the clue of self-organization of long wavelength of stress, which was not measured earlier with classical approaches. This type of study may provide a very useful information for hazard mitigation.


Introduction
The occurrence of high frequency earthquake in any region reveals the development of stress in the crust that leads to deformation. Such areas must be focused for understanding the mechanism of earthquakes and hence to find out the probable potential zone for future large size earthquakes. Several workers have applied various methodologies to find out the precursory signature for future earthquakes, but a working model is yet to be achieved. Puzzling earth scientist has the quest to understand these complex tectonic patterns due to the heterogeneity of the crust and non-availability of any such model for precursory study. However, globally earthquakes rarely seem to be having surface outcrops of active faults. Most of the proposed models are of geological rather than actual physical base. Thus, with a view to have a more valuable perceptive of the physical model of earthquake occurrence, the fractal statistics have been used in which crustal deformation pattern may be studied using fractal dimension. Mandelbrot (1983) was the first to put the idea of fractal geometry and dimension to describe the scale invariance in natural phenomena. Latter, the concept of fractal structure has been applied by several researchers (Hirata., 1989;Feng and Seto., 1999;Dimitriu, et al., 2000;Kayal, et al., 2012;Roy and Mondal., 2012;Han, et al., 2015) to understand the complex mechanism of occurrence of earthquakes. Moreover, distribution of earthquakes and its pattern of evolution in different parts of the globe were also analysed (Legrand et al., 1996;Oncel and Wilson., 2002;Roy et al, 2012;Roy and Gupta., 2015;. Recently Michas, et al., (2015) have shown the significance of earthquake study using fractal analysis.
Scientific community (Hirata et al., 1987) has also reported the laboratory experiment of decrease in fractal dimension with the development of rock fracturing process. Similarly the temporal variation of fractal dimension (Dc) of seismicity gives the idea of preparation process for the large size earthquakes. Several large size earthquakes in different parts of the world have been studied statistically and have been found to be associated with some clustering of earthquakes before the main shock earthquake (Roy and Nath., 2007;Roy and Mondal., 2010). The scaling nature in distribution of earthquake occurrences in space and time domain are the general characteristics, which are required for understanding the physics of earthquake leading to seismic hazard assessment (Sobolev et al., 1991;Tomilin and Voinov., 1995;Knopoff et al., 1996;Telesca et al., 2003aTelesca et al., , 2003bMukhopadhyay et al., 2004;Zamani and Agh-Atabai., 2009;Teotia and Kumar., 2011;Özturk., 2018).
The seismicity distribution in New Zealand ( Fig.1) shows that the deformation is widespread, diffused and persistent. Here the scattered nature of seismicity is evident on both regional and local scales (Fig.1). New Zealand is located along a zone of contact between Pacific and Australia plates. The motion of these plates leads to the major seismicity in this country. The 15 th July 2009 earthquake was the largest in the past 80 years. This thrust-related earthquake occurred within the Fiordland subduction zone of southwestern South Island at 9:22 UTC. The Mw 7.8 earthquake ruptured the plate boundary interface between the subducting Australia plate and the overlying Pacific plate (Beavan et al., 2010). Mahesh, et.al (2011) reported that this quake has significantly increased the Coulomb stress on the overriding plate particularly on the offshore portion of the Alpine fault.
Here we have applied two statistical parameters Dc and seismic b-value fluctuation with respect to one another for more than thirty years to see the anomalous behaviour before the New Zealand earthquake (Mw 7.8). Further, Coulomb stress change prior to event has been studied to show the significance of Dc and b-value study.

Geo-tectonic setting of the study area:
Intense seismicity is experienced by the people of New Zealand, due to its tectonic plate motion (Fig.1). Australia and Pacific plate's boundary run along the Puysegur Fault, Puysegur Trench, the Alpine Fault and the Hikurangi Trench. The country is adjacent to the subducting oceanic crust of Pacific plate below the continental crust of Australia plate to the northeast at Hikurangi margin. On the other hand, Australia plate is subducting below the Pacific plate to the southwest of the country, at Puysegur margin. A very long Alpine fault of 650 km crosses the country and link the two subduction margin. During the early Miocene, displacement of about 480 km was estimated due to strike slip deformation along the Alpine Fault. New Zealand's Southern Alps is considered to be formed due to the estimated shortening of about 70-90 km across the Alpine Fault, during the last 6.4 Myr. The plate motion is varying between 35-40 km/Myr from northern Puysergur Trench to Southern Hikurangi margin (De Mets et al., 1994). The intense seismic activity in the Fiordland region was marked with both shallow and intermediate earthquakes. Shallow earthquakes of dip-slip and strike-slip focal mechanisms are found, distributed in a narrow zone along the Fiordland margin (Smith and Davery, 1984;Reyners, 1989;Anderson and Webb, 1994).Tectonically and geographically the region is very complex.

Data set
For the fractal analysis of the New Zealand earthquake, we prepared the catalog using the global catalog of USGS and catalog of Geo-Net, considering earthquakes of (mb≥ 3.7) during the period 1973-2009. In our analysis epicentral distance based fractal analysis was considered but not the hypocentral location due to the presence of higher error values in hypocentre. The completeness magnitude plays a vital role in finding out the significance of stress analysis. Considering the significance of statistical analysis, we merged the two catalogue to make our data set complete for this type of analysis. We estimate the completeness magnitude using classical Gutenberg-Richter (1944) relation as given in Fig.2a. As we know that the Earthquake scaling and frequency of occurrence relations requires that small earthquakes will be just as important as larger ones in redistributing the forces that drive relative displacements across active faults of any dimension, including plate boundaries (Hanks, 1992). So we have taken care of that in our analysis. While considering the catalogue of USGS and local catalogue of Geo-Net, for the entire periods from 1973-2009, we felt that the earthquakes which are missing in one agencies are filled up by the other agency which are very much necessary to give better idea about the stress pattern / seismicity distribution. Here we focussed in finding the pattern of Dc value only before the main event of 15 th July, 2009. So the data set after 2009 has not been taken for the analysis. The temporal variation of Dc and b-value are estimated by taking successive 100 events' windows in the latitude range 40°S -50°S and longitude range 162° E -178° E. (Fig.1). We have taken all parameters from Global Centroid Moment Tensor (GCMT) catalog for Coulomb stress analysis.
Correlation function is defined as: where N = total number of vector pairs, rij = distances between the points of a set H = Heaviside step function. r = length scale Fractal structure of the epicentral distribution is obtained by the relation given below:- The plot of C(r) versus r in logarithmic coordinates is used to get correlation dimension (Dc).Two plots for log C(r) versus log r are shown in Fig 2b for low Dc value before the main event and Fig.3c for high Dc value for the window of main event.

Seismic b-value
The b-value has been computed using the classical Guttenberg Richter relation (Gutenberg & Richter., 1944), and the equation is given as :- Where N(m) is the number of earthquakes with magnitude larger than or equal to m, 'b' is scaling parameter and 'a' is a constant. The plot for data set is shown in Fig.2a and the power law is valid. Self-similar property of earthquakes occurrence is considered in this law. The constant 'a' depends on the overall seismicity of the region. The 'b'value is normally 1, but it varies from 0.5 to 1.5 depending on the seismicity and tectonic setting of a seismically active region Jones, 1989 Kisslinger andJones, 1991).
Here we have considered sampling of 100 events for estimation of b-value with standard deviation of ~ 0.1b. The error value increases as the number of sampling earthquakes decreases. We can estimate b-value with minimum of 25 events having standard deviation as low as ~ 0.25b (Bernice., 1983).The estimation of b-value is more reliable if the sampling subset contains minimum of 100 earthquakes. The higher the size of subset more reliable it is (Utsu, 1965).

Coulomb stress change
Coulomb Stress change may give the idea of stress change in the earth crust due to occurrence of the earthquakes prior to large size event and hence the pattern of seismicity in the study region. Observation of increased stress due to an event may explain how other events may be induced by the previous one in the fault zone. The earthquakes that occur in the increased static stress are supposed to have been advanced toward failure by the positive stress increment, where the negative stress region are considered to be more relaxed. Static stress change may be the cause of rupturing processes due to previous events. It has been seen that the small, medium and large aftershocks generally do occur in the region of 'stressed-up' by a mainshock or damaging earthquakes and that generally do not occur in the 'stress-shadowed' region (Ellsworth et al.,1981;Simpsons and Reasenberg, 1994;Jaume and Sykes,1996). Researchers have also shown the relation between stress change and seismicity rate using some studies on the Coulomb stress relative to some global large earthquakes (Reasenberg and Simpon;Stein et al 1992;King, et .al, 1994;Stein, 1999). Here, we have applied in this study to see how the stress is correlated with the low Dc value patches which are considered as the possible highly stressed area.
Here, we have computed the stress changes for entire study region by considering the 2009 earthquake (MW 7.8) as mainshock earthquake source, whose strike, dip and rake values are 26°, 24° and 140° respectively (Mahesh et al., 2011). Coulomb failure stress, CFS has been computed in the elastic half-space associated by mainshock on the assumption of Poisson's ratio of 0.25 and shear modulus values of 32 GPa (Okada, 1992). Furthermore, we have utilized the equation (4) to examine the coulomb stress changes given by King et al., (1994), which is related to the shear stress changes and normal stress changes (King et al., 1994;Stein, 1999;Shan et al., 2013;Liu et al., 2017).
CFS 0 fault plane loaded CFS 0 fault plane relaxed where, is shear stress changes, ՛ is the friction coefficient and is normal stress changes (positive for unclamping). For CFS calculation, the constant friction coefficient value has been assumed for fault ՛ = 0.40 and also uniform slip model has been considered to estimate the CFS for the entire study region.
We checked the result through two different tests for entire region by using different values of coefficient of friction 0.5 and 0.6 to check the reliability of results and it has been found that very small changes are existing and that can be neglected. Previously, several authors (Ellis et al., 2006;Boulton et al., 2012;Barth et al., 2013;) utilized the low coefficient of friction in the range of 0.12-0.37 but it was found that these values show the weak results. Moreover, the value of coefficient of friction ( ՛ = 0.40) has been considered as standard value by several authors Toda and Enescu 2011;Xiong et al., 2017;Pandey et al., 2019) for different tectonic setting to compute the Coulomb stress changes that's why we have used the coefficient of friction value 0.4 for this study region. Generally, the coefficient of friction is dependent on the category of rocks (Byerlee and Brace 1968;Byerlee, 1978) therefore entire result is based on rock type.
We have used this parameter to observe the change in stress in two different plots, one with three time windows events that produced low Dc value and second with the entire data set used in the study region to see the pattern of stress transfer. We have estimated Coulomb stress changes by utilizing all data available from GCMT catalogue for three low Dc time windows and the entire Dc time windows from GCMT catalog.

Results and discussion
Fault systems for a seismically active region are multifaceted natural arrangement unveiling scale-invariant or fractal deformation. The manifestations of earthquakes follow the fractal distribution property both temporally and spatially. The distribution of fault systems and earthquakes manifestation may be quantified by a single value called fractal dimension (Dc) which follow power law distribution. The deformation in general obey power-law scaling of fault length or earthquake distribution. We can measure both the heterogeneity and the clustering of earthquakes evolution using fractal dimension of earthquakes' scattering (Oncel et al., 1996). A highly stressed region is associated with dense clustering of smaller and medium size earthquakes which indicates strain release of the accumulated energy due to tectonic plate motion. These highly clustered region of earthquakes may be the most possible weak zones for potential site for future large size earthquake.
High clustering of events in a region has been observed with low fractal dimension value leading to information of seismicity before large earthquakes (Telesca et al., 2001;Roy and Padhi., 2007;Roy and Mondal., 2010). Significance of occurrences of low fractal dimension associated with high clustering before the large earthquake have also been reported by many researchers (Hirata et al., 1987;Roy and Ram., 2006;Roy and Nath., 2007;Roy and Mondal., 2010 ). The degree of clustering of earthquakes is inversely proportional to the Dc value. So the lower value is associated with high clustering and vice-versa for a region (Hirata et al., 1987;Kagan et al., 2006). Fractal correlation dimension (Dc) of seismicity of a region is utilised to find out the estimated clustering of earthquakes.

Correlation fractal dimension (Dc) analysis
The present study for the New Zealand earthquake (Mw 7.8) shows a very good and clear indication of highly clustered events before the mainshock. We observed densely clustered events both in the case of using entire data set for this study (Fig. 1), as well as more prominently in the case of three low Dc value contributing events (Fig.3). These low Dc values are observed before the mainshock (Fig.4). Diagnostic pattern is revealed with highest Dc value in the year 2003 and one lowest Dc value in the year 2007 prior to the 15 th July 2009 earthquake of Mw 7.8 as shown in Fig. 4 and Table 1. This anomalous variation in Dc value occurred just few years before the event. In the process of the mainshock three consecutive low Dc value were observed, which are the most clustered pattern of seismicity (Fig.3) (Table 1). Dc value of 100 events windows clearly depicts the fluctuating nature of statistical value called Dc with time (Fig.4) and Table 1. However it may be noted that the low Dc value characteristics is unique with space and time. The reason behind is particularly due to the possible approach of algorithm to quantify actual relative clustering spatiotemporally. Further it is known to us that the physical or geological facets change with space and time. Especially the geochemical, thermodynamic etc. processes controlling the elasto-dynamics of the entire seismogenic study zone. In fact in support of these we can see the how the earthquakes cluster may easily be observed with the data set as shown in Fig.3 and Fig.1. Similarly the Dc patches are plotted using the above data set and we find the lowest Dc patch happened to be close to the mainshock (Fig.5) and (Fig. 6). The low Dc contour as observed in Fig.5 gives clear indication of highly stressed zones (Roy and Padhi., 2007;Oncel et al., 1996;Roy and Mondal 2010) and the main event occurred close to it. Essentially the Coulomb stress study also shows that the main event has occurred in the 'stressed-up' region ( Fig.7). Here the clustering pattern of intermediate size events prior to main event is an example of strain liberation in self-organized fashion. Stress wavelength especially the longer one i.e regional one tries to get correlated with the smaller one in local scale. But this correlation can be observed only by seeing the clustering of earthquakes in fractal scale invariance or statistically self-similar fashion. We observed highly anomalous variation in Dc value for this large size earthquake.
The study area has also experienced other strong events. The strong events in different time windows is listed in table 1. Relative low Dc value for the earthquake Mw 7.2 (2003) is also observed in Fig. 4. The clustering of earthquakes prior to this event is also seen (Fig.8). The clustering of events prior to other strong earthquakes are also noticed (Fig. 9). Other two strong earthquakes of magnitude Mw 7.1 (2004) and Mw 7.4 (2007) are observed with clustering pattern Fig.9 & Fig.10. Mw 7.2 (2003) earthquake occurred after two decades of event Mw 7.2 (1981).On the other hand events Mw 7.1 (2004) and Mw 7.4 (2007) occurred during the first and third years respectively after the event of Mw 7.2 (2003). The occurrence of these two events of magnitude 7.1 and 7.4 within a very short period indicate the probable reason for the fracturing process of the main event of Mw 7.8 in 2009. We do not have the sufficient data to see the Dc value variation for the event of Mw 7.2 (1981). Low Dc value variation prior to strong events can be observed in Fig.4.

Seismic b-value analysis
Similar to the Dc value anomaly, we are able to see such significant changes in b-value before the mainshock. Changes in the b-value are seen in the plot of frequency versus magnitude as shown by Fig.11. The plot of b-value clearly depicts the highly abnormal changes in the b-value since August 2003 Fig. 11 and Table 1.The event is observed to have occurred after the abnormal increase and decrease of b-value. There were rapid variations in Dc and bvalue almost six years before the mainshock. The b-value of each 100 events windows were determined and listed in Table.1. The lowest b-value is 0.37 and highest is 1.32. Highest Dc value of 1.65 corroborates with the b -value of 1.29. This may be supported with the power law validity for both the methods with the dispersed seismicity. Three low b-value patches of Fig.12 show the anomalous changes near the main events, which is not reflected by b-value of all the data set used (Fig.13). Further, Coulomb stress change for three low Dc time windows clearly indicate the event has occurred in the positive CFS region (Fig.7), which is again not reflected in the entire data set (Fig.14) that used for Coulomb stress study.
The change or variation in b-value from usual value of 1 is dependent on other related factors of the study region. High b-value was reported by Mogi (1962) due to increase in the material heterogeneity or crack density. Similarly decrease in b-value was noticed due to increase in applied shear stress (Scholz., 1968a;1968b),or an increase in effective stress (Wyss., 1973 ). The b-value gradually decreases before a major fracture, even under a constant load (Mogi 1981). So, there is a possibility of gradual decrease in the b-value before a large size earthquake. We observed increase and decrease in b-value prior to this event as in Fig.11.No matter whatever, the mechanism may be, it is generally recognized that there is a tendency for the b-value to decrease prior to a large earthquake. So undoubtedly this provides a clue in foreshadowing large size earthquakes. But this does not necessarily mean that a large earthquake will occur when there is a decrease in the b-value. Essentially this observation will help us to understand the crustal activities for future study.

Coulomb Stress analysis
We have calculated the Coulomb stress changes of 126 earthquakes dataset for the entire region of the study area, as depicted in the Fig.14. From this figure, it has been noticed that the highest stress changes are obtained in the northern part of the study region i.e. towards the Australia plate. On the other hand Pacific Plate region is dominated by low-stress changes. Low or average stress changes have been observed in the rest of the study area. The Coulomb stress change ranges lie between 0.05 to -0.05 bar. Fig.8 depicts the stress change due to three low Dc time windows. Low to high stress variation is noticed in different parts of the study region.
Distributions of Coulomb stress estimated using the entire data set, which shows that the event Mw 7.8 has occurred in the negative values of CFS region Fig.14.Then the area should be considered as stable due to low stress and hence no large earthquake should be expected. On the other hand, Coulomb stress obtained by using the three low Dc time windows clearly indicates that the event has occurred in the positive values of CFS region. (Fig.7). Thus the low Dc patches (Fig.5) which is considered the possible highly stressed region corroborates well with stressed region those obtained by Coulomb stress change, as depicted by Fig.7 for three low Dc time windows. Not only this, distribution of low b-value patches showed the abnormal pattern of the events that has occurred in that region as given by Fig.7 and Fig.12. We observed that the low Dc patches which are possibly highly stressed region found to correlate well with other seismotectonic parameters.

Conclusion
The aim of this study is to observe how the variation of Dc and b-value for New Zealand earthquake of Mw 7.8 event behave prior to mainshock and its possible correlation with Coulomb stress. We observed highly anomalous fluctuation of these two statistical parameters prior to the main shock. On the other hand, we don't find such high fluctuation in the values of these parameters in the region with other strong events, though we have observed drop in Dc value. This may be the possible reason for instability in the area for large size earthquake. Correlation of Dc/b -value with Coulomb stress shows that the Mw 7.8 event has occurred in the positive value of CFS region as observed from three low Dc time windows, as shown by Fig.5, Fig.7 & Fig. 12. This result also show that the possible highly stressed region, as seen in Fig. 5 corroborated well with the Coulomb stress change Fig.7 for the considered Mw7.8 event.
These two seismotectonic parameters that have an important role to play in providing information of pattern of seismic activity of the region. We may come up with very good relation with GPS based crustal deformation observation for longer period, and hence segregate the highly hazardous zone for future large size earthquake. But the important thing is to be noted that for this type of analysis, the data set should be of large quantity and should be complete even for lower magnitude earthquakes. Fig.1 Tectonic map of study area with significant fault system and seismicity in and around the main event that occurred during the study period. Fig.2b Plot for log C(r) versus log r is shown for low Dc value before the main event. Fig.2c Plot for log C(r) versus log r is shown for high Dc value for the window of main event. Fig.3 Spatial distribution of seismicity that contribute to three low Dc value.These events have induced low Dc value windows prior to main shock.