¹Department of Geophysical Engineering, Division of Seismology, Istanbul University, 34850 Avcilar-Istanbul, Turrkey.

2Department of Geology and Geophysics, University of Edinburgh, Grant Institute, West Mains Road, Edinburgh EH9 3JW, UK.
 
 
 
 

Abstract. Seismically-active fault zones are complex natural systems exhibiting scale-invariant or fractal correlation between earthquakes in space and time, and a power-law scaling of fault length or earthquake source dimension consistent with the exponent b of the Gutenberg-Richter frequency-magnitude relation. The fractal dimension of seismicity is a measure of the degree of both the heterogeneity of the process (whether fixed or self-generated) and the clustering of seismic activity. Temporal variations of the b-value and the two-point fractal (correlation) dimension Dc have been related to the preparation process for natural earthquakes and rock fracture in the laboratory. These statistical scaling properties of seismicity may therefore have the potential at least to be sensitive short-term predictors of major earthquakes.

The North Anatolian Fault Zone (NAFZ) is a seismically-active dextral strike slip fault zone which forms the northern boundary of the westward moving Anatolian plate. It is splayed into three branches at about 31°E and continues westward toward the northern Aegean sea. In this study, we investigate the temporal variation of Dc and the Gutenberg-Richter b-value for seismicity in the western part of the NAFZ (including the northern Aegean sea) for earthquakes of Ms > 4.5 occurring in the period between 1900 and 1992. b ranges from 0.6-1.6 and Dc from 0.6 to 1.4.

The b-value is found to be weakly negatively correlated with Dc (r=-0.56). However the (log of) event rate N is positively correlated with b, with a similar degree of statistical significance (r=0.42), and negatively correlated with Dc (r=-0.48). Since N increases dramatically with improved station coverage since 1970, the observed negative correlation between b and Dc is therefore more likely to be due to this effect than any underlying physical process in this case. We present this as an example of how man-made artifacts of recording can have similar statistical effects to underlying processes.
 
 

1 Introduction

Since the development of the fractal description of the geometry of natural objects (Mandelbrot, 1982), it has been recognized that many complex space-time phenomena like seismicity may be described and interpreted in terms of fractal sets with power-law scaling. For example the b-value in the Gutenberg-Richter (1954) relation

log N = a - bm,                                                                                                                             (1)

implies a power-law relation between frequency of occurrence and the radiated energy, seismic moment or fault length, and is one of the most widely used statistical parameter to describe the scaling properties of seismicity. In this relation, N is the number of events with magnitude greater than m in a region within the studied time period, and a is a constant. The Gutenberg-Richter b-value has even been suggested as a generalized 'fractal' dimension of earthquake magnitude (e.g. Aki, 1981; Turcotte, 1992), even though it is not strictly a Haussdorff dimension. It has been pointed out that this interpretation rests on the assumption of a dislocation model for the seismic source and requires a scale-invariant recurrence interval (Sornette, Vanneste & Sornette, 1991). However it has nevertheless been proven to apply on empirical grounds for tensile fracture in the laboratory (Hatton, Main & Meredith, 1993). Nevertheless the question of application at a larger scale remains open, although it is an implicit assumption in much of the recent literature on the non-linear dynamics of earthquakes. However, the b-value has shown apparent systematic variations before and during major earthquakes (Smith, 1986) and in laboratory tests under controlled conditions (Main et. al., 1990).

Many studies have shown that the distribution of earthquakes are fractal in both time and space (e.g. Kagan and Knopoff, 1980; Sadovsky et al., 1984). Fractal structure also exists in the seismicity of microfracturing in rocks, notably for the b-value (Mogi, 1962; Scholz, 1968; Main et al., 1990) and the fractal correlation dimension (Hirata et al., 1987). A number of methods for calculating different 'fractal' dimensions have been developed and applied to investigate various spatial and temporal aspects of seismicity, fault structure and morphology (Turcotte, 1992). It is clear that a complete description of the fractal character of seismological data requires more than one dimension or scaling exponent, and these may also change at different locations or scales (Sornette, Vanneste & Sornette, 1991).

The North Anatolian Fault Zone (NAFZ) is one of the most active dextral strike slip faults of the world (Fig. 1). It has a total length of about 1500 km, extending from the Karl?ova triple junction (labeled K in Fig. 1) in eastern Turkey to mainland Greece (Barka, 1992). Dextral motions on the NAFZ are known to have taken place since the late Miocene (?engör, 1979). About 900 km of surface rupture has been created along the fault zone from Erzincan to Adapazar? (between 31°E and 41°E) by a series of six large (Ms>7) westward migrating earthquakes between 1939 and 1967 (Ambraseys, 1970; Barka, 1992). During this suggested migration, some segments remained unbroken, and have been suggested as potential 'seismic gaps' to be broken in the future (Toksöz et al., 1979). Figure 2 shows the distribution of seismicity in the area, with the area of present interest outlined, and Figure 3 (after Öncel et al., 1994a) shows the temporal distribution of epicenters along strike for earthquakes of Ms>6 occurring along the trend of the whole NAFZ between 1900 and 1992. An examination of Fig. 3 indicates only a weak westward migration for events of this magnitude, with pronounced fluctuations about the secular trend. Beginning from Adapazar? at about 31°E, the NAFZ splays into three strands (southern, middle, and northern strands) in the Marmara and the north Aegean regions (Barka and Kadinsky-Cade, 1988). Dextral motion on the western part of the NAFZ is clearly expressed by an abundance of physiographic and geologic features (Allen, 1969; Barka, 1992). In the western part of the NAFZ, the northern and the middle strands are called the ?zmit-Sapanca and Iznik-

Figure 3: Temporal evolution of the longitude of earthquake epicenters along the North Anatolian fault zone for magnitudes Ms6 (after Öncel et al., 1994a). A weak westward drift is indicated by secular trend indicated by the dashed regression line (r=-0.31).

Mekece faults, respectively (Sipahio?lu and Matsuda,

the1986). Also, seismic studies indicate that the activity

along the northern strand has been much higher than the other strands (Crampin et. al, 1985; Üçer and Alptekin, 1988).

This paper examines temporal variations in the fractal characteristics of seismicity in the western part of the NAFZ, including events in the northern Aegean sea (longitudes between 24°-31°E) over a period from 1900 to 1992. In particular we examine the effect of improved station coverage on the calculation of appropriate scaling parameters. The seismicity data in the area as a whole show no strong overall separation into strands, and are therefore analyzed in a single exercise. In a separate paper we investigate a similar temporal evolution of seismicity for the more well-defined north Anatolian fault zone between longitudes 31°-41°E (Öncel et al., 1994a). Spatial variations of these properties for the whole Anatolian complex are examined in Öncel et al. (1994b). The data for this time period provide an opportunity to examine the temporal variability of the b-value and the fractal (correlation) dimension (Dc) of earthquake epicenter locations. Previous studies have shown varying degrees of negative correlation between these parameters (e.g. Hirata, 1989), which must be accounted for in any physical model of seismicity (e.g. Main, 1992). We show here that such an apparent correlation may also occur in an analysis of temporal variations of the fractal properties of seismicity due to improvements in station coverage.

2 Method of analysis

b and Dc were estimated for sliding windows of 100 consecutive events, the windows being advanced by 10 events between each calculation. In almost any regional tectonic zone this number of events is at present all the data will allow for studies over this time scale. Nevertheless, it might be questioned whether a set of 100 events for a fractal computation is sufficient, or valid, especially given that such estimates are commonly made on the basis of many thousands of data points (Kurths &Herzel, 1987). This question was adressed in an example of a non-stationary system (Havstad &Ehlers, 1989) in which estimates of the fractal dimension made with a data set of 100 events were within 10% of the estimate made with 12800 points. This corresponds roughly to the order of magnitude error estimated in the presented study, implying that the method used does

Figure 4: Cumulative (solid line) and discrete (dots) frequency-magnitude statistics for the area and time period interest. The roll-off in the frequency of small events indicates that the catalogue is substantially complete for Ms4.5.

not significantly underestimate the errors involved. However the relative changes reported here are probably more significant than absolute values as a result. A window containing a fixed number of events n, rather than one pertaining to a fixed period of time, was chosen since it permits one to control the reliability in estimating b and Dc.

The b-value is estimated using the maximum likelihood method (Aki, 1965):

,                                                                                                                                       (2)

where  is the average magnitude and mo is the threshold magnitude for complete reporting of the earthquake magnitudes. The 95% confidence limits for this estimate are ± . The b-value varies between 0.5 and 1.6 in the present work, implying a typical error at this confidence level of ± 0.1-0.2 for n=100. Figure 4 shows that log-linear scaling is appropriate for events in the catalogue above magnitude 4.5 Ms, implying that the catalogue is substantially complete in this magnitude range. The threshold magnitude mo is taken as 4.45 in this study, because magnitudes are usually quoted to one decimal place.

The temporal distribution of seismicity is characterized by the fractal dimension of the distribution of earthquake epicenters. In principle, the fractal dimension of hypocentres could be calculated to avoid the possible introduction of bias due to the projection errors. However the depths of earthquakes are much less well determined than their horizontal positions, so it was considered preferable to use epicenters. The fractal correlation dimension of the distribution of earthquake epicenters for the same 100-event windows was estimated using the correlation integral method (Grassberger and Procaccia, 1983). There are other approaches to the quantitative description of the temporal distribution of the fractal properties of earthquakes (e.g. the box counting algorithm - Sadovsky et al., 1984), but the correlation technique has been more extensively used previously due to its simplicity and reliability (e.g. Kagan and Knopoff, 1980; Hirata, 1989; Henderson et al., 1992). The correlation dimension depends on the distance between epicenters, so that fluctuations in this parameter can be more directly related to physical models based on fluctuations in the elastic interactions between individual earthquake events (e.g. Main, 1992).

The correlation dimension, Dc, is found from:

where r is the distance between two epicenters and C(r) is the correlation function:

                                                                                                                         (4)

where N is the number of points separated by a distance R less than r. The distance r (in degrees) between two events is calculated from:

                                                                      (5)

where(qi, fi) and (qj, fj) are the latitudes and longitudes of the ith and jth events respectively (Hirata, 1989). The correlation dimension is estimated by fitting a straight line to a plot of log C(r) against log r (using 1°=111km) over a data range where a power law can be justified over 1.5 orders of magnitude. The lower bound is usually determined by the epicentral resolution, and the upper bound by the influence of the finite size of the study zone (Kagan and Knopoff, 1980). The observed lower and upper bound in this study for power-law scaling were found to be in the range 5km<r<160km, consistent with this suggestion. We estimated the value and standard error in Dc by linear regression over this range.
 
 

3 The data

The earthquake catalogue for the NAFZ used in this study was compiled from the International Seismological Center Data File, augmented by the Turkish national network and by many local seismic networks, with good station coverage particularly since 1970 (Öncel et al., 1994a).
Figure 5: Example calculations for the correlation dimension and the seismic b-value at three different time periods. In (a), (b) and (c) the left-hand diagram gives the whole correlation plot, and middle diagram gives the actual regression fit to the 1.5 orders of magnitude of data for which a linear fit could be obtained (r here in km). The right-hand diagrams show the frequency-magnitude distributions for the same two time periods.

Events in the time period 1900-1992 are available, quoting surface-wave magnitude Ms. The station coverage in recent times is particularly good for the western Anatolian zone. (More reliable seismological data in Turkey have been obtained after WWSSN was set up in 1964). In this study, we divided our data set into 25 time slices for the western part of the NAFZ.
 
 

4 Results

Figure 5 displays examples of the calculation of b and Dc for the western part of the NAFZ (24°-31°E). The correlation plots for three different time windows (Figs. 5a, b & c) exhibit broad scale invariance both in the discrete and cumulative statistics between 0.7°<log r<2.2°, or 5km<r<160km, giving regression coefficients of the order r=+0.98 in the cumulative statistics, and r=0.74 to 0.98 in the discrete data. It is important to establish significant correlation in the discrete data since cumulative data have the effect of artificially smoothing any real scatter, and would for example produce a positive regression coefficient  even for a random discrete data set. The frequency-magnitude distributions for the same time windows are also shown in Fig. 5, and exhibit a regression coefficient in the cumulative statistics of r=-0.98. Figures 5 b and c show second-order fluctuations about the scale-invariant trend at a length scale of 20 km. This trend is not matched in Fig 5a for the early part of the catalogue. Bumps on the frequency-magnitude curves in Figs 5 b and c can also be seen at a magnitude 6 or so, corresponding to a fault length of 10-20 km. This is also less marked in Figure 5a. One possibility is that the better-determined, more recent data are capable of resolving these second-order effects. The most likely explanation of why this particular length scale shows some deviation from strict scale-invariance is the finite width of the predominantly seismogenic zone in the area (about 10 km). Figure 6 shows the temporal variation in b and Dc together with their standard errors. The fluctuations are statistically significant at this level (66% confidence) but also at the 95% confidence level, since in both cases the fluctuations are also significant at 95% confidence. The arrows on Figure 6 mark the time periods used in Fig. 5. Figure 6 also shows that the negative correlation apparent in the temporal variation of b and Dc is confirmed by regression (r=-0.56).
Figure 6: The top diagram shows the temporal evolution of the b-value and correlation dimension, plotted with error bounds at one standard deviation. The three arrows correspond to the time periods shown in Figure 5. The bottom diagram shows the regression plot of b and Dc. Both diagrams imply a negative correlation between b and Dc.
 
 

Figure 7(a) shows the temporal variation in seismic event rate N (i. e. n=100 events divided by the size T of the variable time window). The event rate shows a gradual Figure 7(a) shows the temporal variation in seismic event rate N (i. e. n=100 events divided by the size T of the variable time window). The event rate shows a gradual

Figure 7: Seismic event rate N as a function of time, b-value and correlation dimension. The event rate is defined by N=n/T, where T is the time period for n=100 events.
 
 

improvement since 1940 (the first available point), with a marked peak in 1980-1990, and a smaller peak in 1955-1965. These results are strongly suggestive of an improvement in station coverage since 1960, and particularly since 1970. In order to examine this possible systematic effect on the trend of Figure 6, Figure 7(b) and 7(c) show individual regression plots of log N vs. b and Dc respectively. log N is used rather than N because it is directly related to the parameter a of the Gutenberg-Richter
relation. log N is positively correlated with b (r=0.42) and negatively correlated with Dc (r=-0.48). Both of these are statistically comparable with the negative correlation between b and Dc (r=-0.56) observed in Figure 6.
 
 

5 Discussion

The results presented here show an apparent negative correlation between b and Dc during their temporal evolution, but with a regression coefficient comparable to that between (log) event rate and these individual parameters. Since event rate increases dramatically with improved station coverage, particularly with the advent of more dense local seismic networks since 1970, we must conclude that this apparent correlation is more likely to be an artifact of improved coverage rather than a result of the underlying dynamics. This is in direct contrast to our previous results on the central and eastern Anatolian fault zone (Öncel et al., 1994a), where the event rate showed no dramatic increase since 1970, and where the correlation between event rate and b or Dc had the opposite sign to those found here, and with regression coefficients of a much smaller magnitude than the much stronger correlation found between b and Dc (r=-0.85). The latter results (Öncel et al., 1994a) were interpreted with a model of increased clustering (higher Dc) and an elevated probability of larger events (lower b) being associated with high event rate (even during times of relatively poor local station coverage), consistent with an increase in stress concentration (Main, 1992). Here we also have higher Dc being associated with lower b, but with a low event rate at a time when station coverage was poor, particularly for events in the northern Aegean sea.

Increased clustering may occur because of (a) alignment of epicenters in 1D along a linear fault strike or (b) 2D spatial clustering of epicenters during fault nucleation or in aftershock sequences (Fig. 1 of Main, 1992). Fig. 2 of this paper shows a strong 1D alignment of epicenters, particularly in the northern Aegean sea, and Figure 5 shows that increased clustering is associated with a rise in the probability of closely-spaced events, consistent with (a) according to Fig. 1a of Main (1992). Thus increased clustering in this paper is associated with the development of more apparent 1D aligned seismicity (lower Dc), higher detection rate and higher b-value (Figures 6 & 7). A higher b-value would result from more complete reporting of smaller magnitude events. Therefore the correlation observed in this paper are more likely to be associated with improved epicentral resolution and better detection of smaller magnitudes in the northern Aegean sea than any underlying dynamics.
 
 

6 Conclusion

Statistics of earthquakes suggests a weak negative correlation between the magnitude-frequency distribution b-value and the fractal distribution of earthquake epicenters (Hirata, 1989; Henderson et al., 1992). In this study, we found a weak negative correlation (r=-0.56) for temporal variations of the b-value and the correlation dimension Dc for earthquakes along the western part of the NAFZ and the northern Aegean sea. However these results are more likely to be an artifact of improved station coverage than the underlying dynamics, particularly for events showing strong epicentral alignment along individual fault zones in the northern Aegean sea.
 
 

Acknowledgements. The work described above was carried out while AOÖ was a visiting researcher at the University of Edinburgh. We are grateful to referees for their constructive comments during the revision of this paper.

References

Aki, K., Maximum likelihood estimate of b in the formula log N = a-bm and its confidence limits, Bull. Earthquake Res. Inst. Tokyo Univ. 43, 237-239, 1965.

Aki, K., A probabilistic synthesis of precursory phenomena, in Simpson, D. W. & Richards, P. G. (eds.) Earthquake Prediction, American Geophysical Union, Maurice Ewing Series 4, 566-574, 1981.

Allen, C., Active faulting in northern Turkey, Contrib. 1577, 32, Div. Geol. Sci., Calif. Inst. Technol., Pasadena, 1969.

Ambraseys, N., Some characteristic features of the North Anatolian Fault Zone, Tectonophysics 9, 143-165, 1970.

Barka, A., The North Anatolian fault zone, Annales Tectonicae 6, 164-195, 1992.

Barka, A. & Kadinsky-Cade, K., Strike-slip fault geometry in Turkey and its influence on earthquake activity. Tectonics 7, 663-684, 1988.

Crampin, S., R. Evans, and Ücer, S. B., Analysis of local earthquakes the Turkish Dilatancy Projects (TDP1 and TDP2), Geophys. J. R. Astr. Soc. 83, 1-16, 1985.

Gutenberg, B. & Richter, C.F., Seismicity of the Earth and Associated Phenomena, Princeton University Press, Princeton, N.J., 1954.

Grassberger, P., and Procaccia, I., Measuring the strangeness of strange attractors, Physica 9D, 189-208, 1983.

Hatton, C.G., Main, I.G. & Meredith, P.G., A comparison of seismic and structural measurements of scaling exponents during tensile subcritical crack growth, J. Struct. Geol. 15, 1485-1495, 1993.

Havstad, J.W.& Ehlers, C.L., Attractor dimension of non-stationary dynamical systems from small data sets, Phys. Rev. A39, 845-853.

Henderson, J., Main, I. G., Meredith, P.G. and Sammonds, P. R., The evolution of seismicity-observation, experiment and a fracture-mechanical interpretation, J. Struct. Geol. 14, 905-913, 1992.

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

Hirata, T., Satoh, T. & Ito, K., Fractal structure of spatial distribution of microfracturing in rock, Geophys. J. R. astr. Soc., 90, 369-374, 1987.

Kagan, Y.Y. & Knopoff, L., Spatial distribution of earthquakes: the two point correlation function, Geophys. J. R. astr. Soc. 62, 303-320, 1980.

Kurths, J. & Herzel, H., An attractor in solar time series, Physica 25D, 165-172.

Main, I.G., Damage mechanics with long-range interactions: correlation between the seismic b-value and the fractal two-point correlation dimension, Geophys. J. Int. 111, 531-541, 1992.

Main, I.G., Meredith, P.G., Sammonds, P.R. & Jones, C., Influence of fractal flaw distributions on rock deformation in the brittle field, in Knipe R.J. & Rutter, E.H. (eds.), Deformation Mechanisms, Rheology and Tectonics, Geol. Soc. London Special Publication 54, 81-96, 1990.

Mandelbrot, B.B., The fractal Geometry of Nature, Freeman Press, San Francisco, 1982.

Mogi, K., Magnitude-frequency relation for elastic shocks accompanying fractures of various materials and some related problems in earthquakes, Bull. Earthquake Res. Inst. Tokyo Univ. 40, 831-853, 1962.

Öncel, A.O., Main, I.G. Alptekin, Ö. & Cowie, P.A., Temporal variations in the fractal properties of seismicity in the north Anatolian fault zone between 31°-42°E, in proc., conference on Deformations, Istanbul, (submitted 1994a).

Öncel, A.O., Main, I.G. Alptekin, Ö. & Cowie, P.A., Spatial variations in the fractal properties of seismicity in the north Anatolian fault zones, Tectonophysics, (submitted 1994b).

Sadovsky M.A., Golubeva, T.V., Pisarenko, V.F. & Shnirman, M.G., Characteristic dimensions of rock and hierarchic properties of seismicity, Izv. Ac. Sci. USSR Phys. Solid Earth. (English Translation) 20, 87-96, 1984.

Scholz, C.H., The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes, Bull. seism. Soc. Am. 58, 399-415, 1968.

Şengör, A. M. C., The North Anatolian transform fault: Its age, offset and tectonic significance, J. Geol. Soc. London 136, 269-282, 1979.

Sipahioğlu, S. and Matsuda, T., Geology and Quaternary fault in the Iznik-Mekece area, in Işıkara, A. M. & Honkura, Y. (Eds.)., Electric and Magnetic Research on Active Faults in the North Anatolian Fault Zone, 25-41, 1986.

Smith, W.D., Evidence for precursory changes in the frequency-magnitude b-value, Geophys. J. Roy. Astron. Soc. 86, 815-838, 1986.

Sornette, D., Vanneste, C. & Sornette, A., Dispersion of b-values in Gutenberg-Richter law as a consequence of a proposed fractal nature of continental faulting, Geophys. Res. Lett. 18, 897-900, 1991.

Toksöz M.N., Shakal A.F. & Michael A.J., Space-time migration of earthquakes along the N. Anatolian fault zone and seismic gaps, Pure Appl. Geophys. 117, 1258-1269, 1979.

Turcotte, D.L., Fractals and chaos in geology and geophysics, Cambridge University press, Cambridge, 1992.

Üçer, S.B. & Alptekin, Ö., Present and historical seismicity of the Marmara region, in Proceedings of the Workshop on Historical Seismicity and Seismotectonics of the Mediterranean Region, Istanbul, 6-36, 1988.
 
 

FIGURE CAPTIONS
Figure 1: Major tectonic features of Turkey. As a result of collision between the Arabian and the Eurasian plates the Anatolian and the Northeast Anatolian blocks are escaping sideways. The North Anatolian Fault and the East Anatolian Fault respectively defines the northern and southern east boundary of the westward moving Anatolian block.

Figure 2: Map showing the epicenter distribution of earthquakes which occurred between 1900 and 1992 in the NAFZ and in the EAFZ. The solid parallelogram encloses the area studied here.