A correlation between the b-value and the fractal dimension from the aftershock sequence of the 1999 Chi-Chi, Taiwan, earthquake

GJI

Seismology

FA S T T R A C K PA P E R

A correlation between the b-value and the fractal dimension from the

aftershock sequence of the 1999 Chi-Chi, Taiwan, earthquake

Chien-Chih Chen,

Wei-Chien Wang,

Young-Fo Chang,

Yih-Min Wu

and Yuan-Hsi Lee

1_{Department of Earth Sciences and Graduate Institute of Geophysics, National Central University, Jhongli, Taoyuan, Taiwan 320. E-mail: chencc@ncu.edu.tw} 2_{Institute of Seismology, National Chung-Cheng University, Ming-Hsiung, Chia-Yi, Taiwan 621}

3_{Department of Geosciences, National Taiwan University, Taipei, Taiwan 106}

Accepted 2006 September 21. Received 2006 September 21; in original form 2006 June 4

S U M M A R Y

We analyse large number of aftershock events from the 1999 Chi-Chi, Taiwan, earthquake (ML7.3) recorded around the epicentre area of the main shock in central Taiwan where events

can be precisely located, due to dense coverage of modern seismometers. The seismicity is characterized by the b-value of the Gutenberg–Richter relation and the fractal (correlation) dimension D of earthquake hypocenters calculated from sliding windows containing more than 100 events. Over a span of 6 months after the main shock, we find a positive correlation between b and D from the aftershock sequence. Our result thus suggests that Aki’s relation D= 3b/c could be possibly applied to one single fracturing process and its related aftermath.

Key words: Aki’s relation, b-value, Chi-Chi earthquake, fractal dimension.

1 R E L AT I O N S H I P B E T W E E N T H E b - VA L U E A N D T H E F R A C T A L D I M E N S I O N O F E A RT H Q UA K E S

The relationship between the b-value of the Gutenberg–Richter re-lation and the fractal dimension D of earthquakes is a widely dis-cussed topic during the last three decades (e.g. Aki 1981; King 1983; Turcotte 1986; Hirata 1989; Wang 1991; Volant & Grasso 1994; On-cel et al. 1996; Lapenna et al. 1998; Henderson et al. 1999; Legrand 2002; Wyss et al. 2004; Mandal & Rastogi 2005). Since Aki (1981) proposed a simple relationship between the b-value and the fractal dimension with a positive correlation D= 3b/c (where c is the slope of log moment versus magnitude relation and is about 1.5), both the positive (e.g. Guo & Ogata 1995; Legrand 2002; Pascua

et al. 2003; Oncel & Wilson 2004) and negative (e.g. Hirata 1989;

Henderson et al. 1994; Oncel et al. 1996; Wang & Lee 1996) cor-relations between those two scaling exponents have been reported and debated in the literature. In some cases (e.g. Henderson et al. 1999; Mandal & Rastogi 2005; Mandal et al. 2005), the correlation could even change from a negative one to positive.

The paper of Hirata (1989) was probably the first article demon-strating the derivation of Aki’s relation D= 3b/c. He in the paper then presented a negative, instead of positive, correlation between the b-value and fractal dimension of earthquakes in Tohoku region of Japan. Later, Wang & Lin (1993) also showed a negative cor-relation between those two exponents in western Taiwan and they explained a negative correlation for observed seismicity might be reasonable due to the isolation of asperities or barriers on the fault

planes, which indicates the calculation of those two scaling param-eters (Hirata 1989; Wang & Lin 1993) was actually produced from events associated with various fault planes/systems.

Yet another interpretation for the negative correlation was pro-vided in Henderson et al. (1999). They found a negative correlation between the b-value and the fractal dimension of induced seismic-ity during the initiation of water injection into a well in the Geysers geothermal area, California. They suggested such behaviour could occur under conditions of rapid stress loading resulted from water injection. When rapid injection rate overcomes the mechanism of dilatant hardening, the diffusion of pore pressure triggers numerous small earthquakes (high b-value) manifested by the spatial cluster-ing (low D). Interestcluster-ingly, a positive correlation between b and D was also observed in Henderson et al. (1999) over the timespans when well injection activity was fairly constant, indicating a slowly loaded system. Mandal & Rastogi (2005) presented a quite similar result from the natural aftershock sequence of the 2001 Bhuj earth-quake in India. It is not clear how the initiation and the induction of the Bhuj main shock and aftershock events could be related to a hydraulic fracturing process, although aftershocks were widely distributed within an ‘inferred’ fluid-filled rock mass with high-density fractures (Mandal & Rastogi 2005). Nevertheless, Mandal & Rastogi (2005) presented a negative correlation for the first 2 months after the main shock followed by a positive correlation for a later period of aftershock activity.

Two disadvantages involved in the calculation of b and D for the abovementioned researches are mainly the mixture of events from different fault systems (Hirata 1989; Wang & Lin 1993)

Figure 1. Map showing the epicentres of the Chi-Chi main shock (black

star) and its aftershocks (red dots) occurred within 6 months after the main shock. Thick black line is the Chelungpu thrust ruptured when the main shock occurred. We analysed the seismicity within the blue square in this study. The UTM coordinate system is used for this map.

and the limited numbers of events used (Henderson et al. 1999; Mandal & Rastogi 2005). Tackling the effects of mixture and lim-ited earthquake population, we therefore re-investigate the correla-tion between b and D from the abounding aftershock sequence of the 1999 Chi-Chi earthquake (Fig. 1). We focus our issue in this paper on whether the positive or negative correlation holds true.

2 C A L C U L AT I O N O F T H E b - VA L U E A N D T H E F R A C T A L D I M E N S I O N O F E A RT H Q UA K E S

We give brief reviews to the estimations of the b-value and the fractal dimension of earthquakes in this section.

2.1 The estimation of the b-value

The number of earthquakes N with magnitude greater than M is re-lated to the magnitude by log N= a − b M, which is widely known as the Gutenberg–Richter relation (Gutenberg & Richter 1944). In this study, the b-value of the Gutenberg–Richter relation was estimated by a weighted least-squares fitting method (Shi & Bolt 1982) for considering the different weighting of data points in the log-normal Gutenberg–Richter relation. The maximum likelihood estimation for the b-value is often adopted and claimed to be a better estima-tion than the general least-squares method (e.g. Hirata 1989; Mandal & Rastogi 2005). However, when simultaneously calculating both the b-value and the following fractal dimension, we prefer using the same regression technique for balancing the potential effect due to the regression technique on both estimations. Besides, we have also

confirmed that, in most occasions of evaluating the b-value, both the maximum likelihood and the weighted least-squares estimators give results within 10 per cent of each other.

2.2 The estimation of the fractal correlation dimension The analysis of correlation dimension (Grassberger & Procaccia 1983) is a powerful tool for quantifying the self-similarity of a ge-ometrical object, a point set in a vector space for example. Let us consider a set of 3-D vectors Xi(i = 1, 2, . . . , N) for, say, the

lo-cations of hypocenters. Grassberger & Procaccia (1983) define the correlation sum C(r) as C(r )= 1 N 1 N− 1 N  i=1 N  j=1, j=i (r −  Xi − Xj), (1)

where is the Heaviside step function, that is,  (x) = 0 for x ≤ 0 and (x) = 1 for x > 0, and || · || is the Euclidean norm. Given a large amount of data points and small correlation distance r, they showed that a power-law relation between C(r) and r holds true and defined the correlation dimension D by

C(r )∝ rD_{. (2)}

For considering the different weighting of data points (r, C) and bal-ancing the potential fitting effect caused by the regression method, as mentioned in the previous section, the weighted least-squares fit-ting method was again used to fit the slope of the log–log plot of C(r) versus r.

3 A N A L Y S I S O F T H E A F T E R S H O C K S E Q U E N C E F R O M T H E C H I - C H I M A I N S H O C K

An earthquake with moment magnitude 7.6 took place in central Taiwan on 1999 September 21 (or at UTC 17:47 20 September) and the epicentre was located near a small town of Chi-Chi (Fig. 1). The Chi-Chi earthquake ruptured an approximately 100-km segment of the Chelungpu fault (Fig. 1), which represents a geological bound-ary separating the foothills from the plain areas in the central west-ern Taiwan. The Chi-Chi event also inflicted severe damage in the central and northern Taiwan. The Chi-Chi earthquake is the largest earthquake to occur on the Taiwan Island over the past hundred years (Shin 2000). The earthquake catalogue (Fig. 1) for Taiwan region is routinely released from the Taiwan Central Weather Bureau (CWB). In this study, we used data of earthquakes occurred in Taiwan area within 6 months after the Chi-Chi main shock, because most of the aftershock events occurred before 2000 March (Wu & Chen 2006). The CWB seismic network was operated in a trigger-recording mode before the end of 1993. Since 1994, the operation has been changed to continuously recording mode and manual identifications of earth-quake events, thus has greatly enhanced the detection sensitivity of earthquakes.

The location error of the CWB catalogue has not been systemat-ically estimated. However, based on previous studies (Cheng 2000; Wu et al. 2003), the error in hypocentral location is within 1 km in western Taiwan. Particularly, Cheng (2000) relocated more than 2000 aftershocks of the Chi-Chi event with his 3-D velocity model and found an average shift of 0.6 km in focal depth when compar-ing with the results in the CWB catalogue. In this study, we used earthquake data with focal depth less than 30 km, which almost involves all the focal depths of aftershocks from the Chi-Chi main shock. For the spatial selection of events, we selected earthquakes

mainly distributed around the surface rupturing fault, that is, the Chelungpu thrust, of the Chi-Chi main shock (Fig. 1) to avoid data mixture from different fault systems as mentioned above. The se-lected region for our analyses was confined to an area with diagonal points of (200 000/2610 000) and (270 000/2700 000) in the UTM coordinate system (Fig. 1), which is with the dimensions of 70 km× 90 km. There were thus more than 26 000 events with magnitude larger than 1 selected for the prescribed space–time window in our analysis.

For evaluating the spatial variations in the b-value and the hypocentral correlation dimension of the aftershock sequence, the sliding-window technique has been used on the earthquake data. If the system is approximately stationary over many relatively local ar-eas, it is possible to determine the fractal dimension by using many small data sets covering a local area (Havstad & Ehlers 1989). The volume for the b and D analyses is then defined as a column with the dimensions of 20 km× 20 km × 30 km, almost cubic in shape. The centres of calculated volumes are evenly spaced out 2 km apart, thus giving the overlap of 90 per cent between two neighbouring columns and producing smooth patterns dense enough for both variations in the b and D values.

Before doing the b and D calculations, for each calculated column, we have determined the magnitude of completeness (Mc) which is

Figure 2. Maps of: (a) the b-value; (b) the fractal correlation dimension D; (c) the magnitude of completeness Mcand (d) the earthquake intensity, that is, the

earthquake number, obtained from the aftershock sequence occurred within 6 months after the Chi-Chi main shock. The UTM coordinate system is used for

each map. Black open circles denote the main shock and large aftershocks with ML≥ 6. Thick black line denotes the Chelungpu thrust. Grid points marked

by solid squares represent the areas with less than 100 earthquakes, and crosses for the areas with the D/b ratio≥ 2.5.

estimated by taking the maximum value of the derivative of the frequency-magnitude distribution (Wyss et al. 1997). We show the obtained Mcmap in Fig. 2(c). Only earthquake data with magnitude

larger than Mcwas used to calculate the b and D values. The b and

D calculations were conducted for the column containing more than 100 events to ensure the reliable values of b and D, and actually most of the columns calculated contain more than 200 events (Fig. 2d). Figs 2(a) and (b) are the obtained maps of the b and D values, respectively. We have blanked out those grid points having less than 100 events with white solid squares on each map in Fig. 2 and plotted the main shock and large aftershocks (ML≥ 6) on the maps. The

standard deviation for b or D is in general less than 10 per cent of the calculated b or D value. It is obvious from Fig. 2(a) that the b-values are low in the areas where large events occurred. It seems reasonable that a large aftershock produced a few secondary aftershocks with moderate sizes and, therefore, lowered the b-value. On the other hand, those large aftershocks look like enveloping the region with D values≥ ∼2.2 in the eastern part of the studied area (Fig. 2b), thus indicating an inward fracturing process from the rupturing surfaces confined by those large aftershocks.

We then have a scatter plot of b versus D for each grid point (Fig. 3). For comparison, in Fig. 3, we also plotted two straight lines, they are D= 2b and 3b (solid lines in Fig. 3). Quite interestingly,

Figure 3. Scatter plot of the b-value and the fractal dimension D calculated

from the aftershock sequence. Also shown in this plot are two straight lines

of D= 2b and 3b. There are more than 65 per cent of points (b, D) lying in

between D= 2b ± 0.6, which are represented by two dash lines.

we could find a fairly well positive correlation between b and D with a slope of∼2. Although the scatter of points (b, D) from D = 2b is quite large and such scattering prevents an effective regression from obtaining a slope of∼2, most points obviously follow the trend of D= 2b. There are more than 65 per cent of points (b, D) lying in between D= 2b ± 0.6 (two dash lines in Fig. 3) and more than 90 per cent lying in between D= 2b ± 0.8. Furthermore, it seems that a cluster of points (b, D) in Fig. 3 has been distributed towards the direction of high D and/or low b, to the upper-left corner. This cluster of points seems following another trend of D= 3b. When we plot a histogram for the ratios of D to b (Fig. 4), it turns out that the D/b ratios have a quite striking bimodal distribution with two modes of∼2 and ∼3. It is quite possible that our obtained data points of (b, D) actually consist of two different relations between

Figure 4. Histogram of the D/b ratio. The D/b ratio has an interesting

bimodal distribution with two modes of∼2 and ∼3.

D and b, that is, D= 2b and 3b, plus some noises from the locating procedure of the catalogue and our calculating process for b and D values. For instance, the D value is possibly overestimated due to the noise involved in the data (e.g. Kantz & Schreiber 1997). An obvious example of noise is the location error of hypocenters, particularly in focal depth. The nature of infinite dimension for stochastic process biases the calculated D towards to a higher value.

We have marked in Fig. 2 the areas with the D/b ratio larger than 2.5, the crosses on each map. As it can be found in Fig. 2(a), the major part of the areas having the D/b ratio greater than 2.5 is related to the areas having low b-values and large aftershocks as well. It is possible that the major body of the earthquake data in those areas is consti-tuted by a large amount of smaller earthquakes, which were mainly the offspring of moderate-sized aftershocks as above-mentioned. Numerous smaller earthquakes dominate our calculations of fractal dimension and are distributed within crustal volumes, thus gener-ating large values in D. We have noted that Legrand (2002) exactly proposed a relation of D= 3b for small earthquakes.

4 C O N C L U D I N G R E M A R K S O N A K I ’ S R E L A T I O N O F D = 3 b/c

Abounding aftershock sequence of the Chi-Chi earthquake gives us a very good opportunity to reinvestigate the correlation between the b-value and the fractal dimension D, which was proposed al-most three decades ago and has been continuously debated (Aki 1981; King 1983; Turcotte 1986; Hirata 1989; Wang 1991; Volant & Grasso 1994; Oncel et al. 1996; Lapenna et al. 1998; Henderson et al. 1999; Legrand 2002; Wyss et al. 2004; Mandal & Rastogi 2005). In this study, we have found positive correlation between b and D, which is much consistent with Aki’s relation D= 3b/c in cases of c∼1.5 and, probably, c ∼1.0 as well (Aki 1981; Legrand 2002). There is no observation on the c value particularly for the Taiwan region and we consider that 1.5 (for intermediate events) and/or 1.0 (for small events) might be reasonable values or at least appropriate to the main shock and aftershock sequence of the Chi-Chi event.

It is worth to note that our analyses of b and D were basically di-rected to one single fracturing process, since we had concentrated the analyzed spatiotemporal scale to the immediate aftershocks nearby the surface rupture of the Chi-Chi main shock. Therefore we provide an evidence for Aki’s relation (Aki 1981) and Hirata’s interpretation (Hirata 1989). It is definitely true that a single fracturing process is not necessary equivalent to a single fault plane. Hirata (1989) speculated the connectivity of the asperity or barrier distribution is a crucial factor affecting the validness of Aki’s relation. Thus, based on the present study, we suggest that Aki’s relation could be possibly applied to a connective fault system. Similar result was also obtained in the locked segment of the San Andres fault near Parkfield, where Wyss et al. (2004) shows that the relationship D approximate to 2b holds. We are now interested in whether the widely reported nega-tive correlations are actually related to the effect of the mixture of different fault systems. Much more work needs to be done to resolve this issue more fully.

A C K N O W L E D G M E N T S

CCC is grateful for research support from the National Science Council (ROC) and the Department of Earth Sciences (NCU, ROC). The effort of Taiwan CWB to maintain the CWBSN is highly

appreciated. The authors are also grateful to Drs D. Legrand and J. Trampert and one anonymous reviewer for their helpful sugges-tions, which have largely improved our manuscript.

R E F E R E N C E S

Aki, K., 1981. A probabilistic synthesis of precursory phenomena, in

Earth-quake Prediction: An International Review, pp. 566–574, eds Simpson,

D.W. & Richards P. G., AGU, Washington, DC.

Cheng, W.B., 2000. Three-dimensional crustal structure around the source area of the 1999 Chi-Chi earthquake in Taiwan and its relation to the aftershock locations, Terrest. Atmos. Oceanic Sci., 11, 643– 660.

Grassberger, P. & Procaccia, I., 1983. Characterization of strange attractors,

Phys. Rev. Lett., 50, 346–349.

Guo, Z. & Ogata, Y., 1995. Correlation between characteristic parameters of aftershock distributions in time, space, and magnitude, Geophys. Res.

Lett., 22, 993–996.

Gutenberg, B. & Richter, C.F., 1944. Frequency of earthquakes in California,

Bull. seism. Soc. Am., 34, 185–188.

Havstad, J.W. & Ehlers, C.L., 1989. Attractor dimension of nonstation-ary dynamical systems from small data sets, Phys. Rev. A, 39, 845– 853.

Henderson, J., Main, I.G., Pearce, R.G. & Takeya, M., 1994. Seismicity in north-eastern Brazil—fractal clustering and the evolution of the b-value,

Geophys. J. Int., 116, 217–226.

Henderson, J.R., Barton, D.J. & Foulger, G.R., 1999. Fractal clustering of induced seismicity in the Geysers geothermal area, California, Geophys.

J. Int., 139, 317–324.

Hirata, T., 1989. A correlation between the b-value and the fractal dimension of earthquakes, J. geophys. Res., 94, 7507–7514.

Kantz, H. & Schreiber, T., 1997. Nonlinear Time Series Analysis, Cambridge University Press, Cambridge, UK.

King, G., 1983. The accommodation of large strains in the upper lithosphere of the earth and other solids by self-similar fault systems: The geometrical origin of b-value, Pure appl. Geophys., 121, 761–815.

Lapenna, V., Macchiato, M. & Telesca, L., 1998. 1/ fβfluctuations and

self-similarity in earthquake dynamics—observational evidences in southern Italy, Phys. Earth planet. Inter., 106, 115–127.

Legrand, D., 2002. Fractal dimensions of small, intermediate, and large earthquakes, Bull. seism. Soc. Am., 92, 3318–3320.

Mandal, P. & Rastogi, B.K., 2005. Self-organized fractal seismicity and b

value of aftershocks of the 2001 Bhuj earthquake in Kutch (India), Pure

appl. Geophys., 162, 53–72.

Mandal, P., Mabawonku, A.O. & Dimri, V.P., 2005. Self-organized fractal seismicity of reservoir triggered earthquakes in the Koyna-Warna seismic zone, western India, Pure appl. Geophys., 162, 73–90.

Oncel, A.O. & Wilson, T., 2004. Correlation of seismotectonic variables and GPS strain – measurements in western Turkey, J. geophys. Res., 109(B11), B11306.

Oncel, A.O., Main, I., Alptekin, O. & Cowie, P., 1996. Temporal variations in the fractal properties of seismicity in the north Anatolian fault zone

between 31◦E and 41◦E, Pure appl. Geophys., 147, 147–159.

Pascua, M.A.R., Devicente, G., Calvo, J.P. & Perezlopez, R., 2003. Simi-larities between recent seismic activity and paleoseismites during the late Miocene in the rxternal Betic Chain (Spain)—relationship by b value and the fractal dimension, J. Struct. Geol., 25, 749–763.

Shi, Y. & Bolt, B.A., 1982. The standard error of the magnitude-frequency

b value, Bull. seism. Soc. Am., 72, 1677–1687.

Shin, T.C., 2000, Some seismological aspects of the 1999 Chi-Chi earthquake in Taiwan, Terrest., Atmos. Oceanic Sci., 11, 555–566.

Turcotte, D.L., 1986. A fractal model for crustal deformation,

Tectono-physics, 132, 261–269.

Volant, P. & Grasso, J.R., 1994. The finite extension of fractal geometry and power-law distribution of shallow earthquakes–a geomechanical effect,

J. geophys. Res., 99, 21 879–21 889.

Wang, J.H. & Lee, C.W., 1996. Multifractal measures of earthquakes in west Taiwan, Pure appl. Geophys., 146, 131–145.

Wang, J.H. & Lin, W.H., 1993. A fractal analysis of earthquakes in west Taiwan, Terrest., Atmos. Oceanic Sci., 4, 457–462.

Wang, J.H., 1991. A note on the correlation between b-value and fractal dimension from synthetic seismicity, Terrest., Atmos. Oceanic Sci., 2, 317–329.

Wu, Y.M., Chang, C.H., Hsiao, N.C. & Wu., F.T., 2003. Relocation of the 1998 Rueyli, Taiwan, earthquake sequence using three-dimensions veloc-ity structure with stations corrections, Terrest., Atmos. Oceanic Sci., 14, 421–430.

Wu, Y.M. & Chen, C.C., 2006. Seismic reversal pattern for the 1999 Chi-Chi, Taiwan, Mw 7.6 earthquake, Tectonophysics, in press.

Wyss, M., Sammis, C.G., Nadeau, R.M. & Wiemer, S., 2004. Fractal dimen-sion and b-value on creeping and locked patches of the San-Andreas fault near Parkfield, California, Bull. seism. Soc. Am., 94, 410–421.

Wyss, M., Shimazaki, K. & Wiemer, S., 1997. Mapping active magma cham-bers by b values beneath the off-Ito volcano, Japan, J. geophys. Res., 102, 20 413–20 422.