Layout 6 A statistical study of aftershock sequences Giorgio Ranalli Annali di Geofisica, Vol. 22, n. 4, 1969. ANNALS OF GEOPHYSICS, VOL. 53, N. 1, FEBRUARY 2010 ABSTRACT A comprehensive statistical study of the phenomenology of aftershock sequences is made in this paper. The spatial distribution of aftershocks indicates that they are mainly crustal events; however, deeper sequences also take place. The analysis of the distribution of aftershocks in 15 sequences with respect to time and magnitude leads to the statistical confirmation of a set of phenomenological laws describing the process, namely, the time-frequency law of hyperbolic decay of aftershock activity with time, the magnitude stability law, and the exponential magnitude- frequency distribution. The hypotheses involved are checked. The grouping of data and the statistical methods employed are chosen according to some basic well·confirmed assumptions regarding the nature of the process. Introduction A comprehensive study of the phenomenology of aftershock sequences must include space, time, and magnitude distributions. Sequences which occurred in many parts of the world have been studied by various authors and the results are scattered in the geophysical literature. The determination of the phenomenological laws describing the aftershock process is a statistical problem, and it is therefore desirable to employ a consistent statistical procedure and to perform appropriate tests. We feel that some of the methods employed in the past were not rigorous; consequently, this paper is an attempt to present a unified procedure for the statistical study of aftershock sequences. For this purpose, 15 sequences have been analyzed in detail, even if some of them had already been studied according to different methods. The study of the space distribution of aftershocks within a given sequence does not present particular problems. Its accuracy depends on the precision with which epicentral coordinates and focal depths are computed. It has been maintained [Page 1968a] that aftershocks are essentially a shallow phenomenon. A review of available information partly confirms this view, but allowance must be made for notable exceptions. The situation is more complicated with respect to time and magnitude distribution. The problem is basically that of finding a statistical relationship between the various quantities involved, and of estimating the parameters appearing in the statistical laws. First of all, when examining the data (which consist of the origin times ti and of the individual magnitudes Mi of the aftershocks), one must be reasonably sure to be dealing with a complete set; that is, ideally no aftershock with M ⩾ M*, where M* is the minimum magnitude detected in a sequence, should be missing in the time interval considered. Moreover, it has been shown empirically by Suzuki [1958] that the mode of grouping the data influences the results: therefore the mode of grouping should be as uniform as possible. Finally, it is desirable to employ a statistical procedure that does not contradict the underlying characteristics of the phenomenon observed, it is necessary to apply some statistical test to check the hypothesis being entertained, and confidence limits on the results must be given. More often than not, one or several of these conditions are not met in the study of aftershock sequences, and the significance of the results is therefore debatable. Among the several sequences for which origin times and magnitude of individual aftershocks have been published, 15 have been selected for detailed study. As far as possible we have tried to include sequences from different geographic regions, even if this implied considering a few sequences whose completeness may be in doubt. Since the determination of the statistical laws is more reliable when data are more abundant, no aftershock sequence consisting of less than 46 shocks has been included. Table 1 lists the aftershock sequences whose time and magnitude distributions are studied in this paper; it includes region of occurrence, literature reference, main shock parameters (t0 origin time; φ0, λ0 geographic coordinates; h0 focal depth; M0 magnitude), minimum magnitude of aftershock included (M*), focal depth of aftershocks (h), and total number K of aftershocks with M ⩾ M* recorded in the first 100 days. (In sequence (1) the first day after the main shock is excluded from the count.) Focal depths are in kilometers; the term «shallow» is taken to mean «crustal», and often (especially in California) «upper crustal». Magnitude are given in the M-scale or ML-scale, with the exception of sequence (1), where the m-scale has been employed. It will be seen that the results of the analysis confirm 43 A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES 44 Sequ ence R eg ion R eference t0 φ 0 λ 0 h 0 M 0 M * h K 1. A laska 1964 2. A leu tian I. 1957 3. L ong B each 1933 4. D esert H ot Sp. 1948 5. K ern C ou nty 1952 6. San Francisco 1957 7. Salinas 1963 8. P arkfield 1966 9. C halk idike 1932 10. W . T hessaly 1954 11. A m orgos 1956 12. M agnesia 1957 13. Z ante 1962 14. C rem asta 1966 15. H aw ke's B ay 1931 A laska A leu tian I. C alifornia ˝ ˝ ˝ ˝ ˝ G reece ˝ ˝ ˝ ˝ ˝ N ew Z ealand P age 1968b D u da 1962 B enioff 1951 R ichter et al. 1958 R ichter 1955 T ocher 1959 U dias 1965 M cE villy et al. 1967 P apazachos et al. 1967 ˝ ˝ ˝ ˝ C om ninak is et al. 1968 B enioff 1951 M ar 28.03:36:13 M ar 9.14:22:28 M ar 10.17:54:08 D ec 4.15:43:17 Ju l 21.11:52:14 M ar 22.19:44:21 Sept 14.19:46:17 Ju n 28.04:26:13 Sept 26.19:20:37 A pr 30.13:02:36 Ju l 9.03:11:40 M ar 8.12:21:13 A pr 10.21:37:13 Feb 5.02:01:43 Feb 3.10:15:00 61˚.0 N 51˚.3 N 33˚34´ N 33˚56´ N 35˚.0´ N 37˚40´ N 36˚52´ N 35˚57´ N 40˚.5 N 39˚3 N 36˚.7 N 39˚.3 N 37˚.6 N 39˚.1 N 39˚.5 S 147˚.8 W 175˚.8 W 117˚59´ W 116˚23´ W 119˚02´ W 122˚29´ W 121˚38´ W 120˚29´ W 23˚.7 E 22˚.2 E 25˚.8 E 22˚.6 E 20˚.1 E 21˚.6 E 177˚.0 E 20 km– shallow shallow 16 km shallow 3.3 km shallow 28 km 19 km 21 km 18 km 20 km 20 km– 8.5 8.3 6.3 6.5 7.7 5.3 5.4 5.5 6.9 7.0 7.5 6.8 6.3 5.9 7.6 4.5 5.9 3.9 3.0 4.0 2.0 1.0 2.0 3.4 3.2 3.5 3.0 3.6 3.4 4.1 < 35 km < 150 km shallow < 35 km shallow < 10 km < 14 km < 12 km shallow ˝ ˝ ˝ ˝ ˝ – 294 2057868 184 16046 17385 299 400 291 139 10371 T able 1.L ist of aftershock sequ ences selected for detailed stu dy. A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES the 3 basic laws describing the phenomenology of aftershock sequences, namely, the time-frequency law (Omori's law) as formulated by Mogi [1962], the magnitude stability law [Lomnitz 1966], and the magnitude-frequency law [Gutenberg and Richter 1954]. Spatial distribution of aftershocks The spatial distribution of the shocks in an aftershock sequence is naturally related to the location of the main shock. The following considerations are based on a comprehensive survey of available data and are not limited to the sequences listed in Table 1. If one traces on a map the boundary of the area in which the epicenters are located, the epicenter of the main shock is usually close to this boundary. Such is the case for all aftershock sequences of large earthquakes wich occurred in Japan from 1923 to 1963 [Matuzawa 1964]. When aftershock activity takes place along a fault segment, as is frequently the case in California, the domain of the epicenters is approximately elliptical with the long axis parallel to the active fault segment; often the main shock occupies, roughly speaking, one focus of the ellipse, and the aftershocks are concentrated toward the two ends. In some sequenze the aftershock epicenters are clustered in a very small area (only a few kilometers in length and width), but this is rather exceptional; usually they are spread out over a much larger area. In the aftershock sequence of the Aleutian Islands earthquake of March 9, 1957 [Duda 1962] the distribution of the epicenters follows closely the trend of tectonic activity along. the Aleutian arc. As to the focal depths of aftershocks, a review of available data by Page [1968a] indicates that, when hypocenter determinations are accurate, aftershocks are shallow events following a main shock which is itself shallow. Aftershock sequences therefore appear to be crustal phenomena, with the majority of shocks clustering in the upper layer of the crust (h ⩽ 20 km). 45 Table 2. Observed frequency of aftershocks. 46 There are, however, some exceptions. In the region of Greece, the Southern Sporades earthquake of March 18, 1926 had a reported depth of 50 km and was followed by 18 aftershocks with M ⩾ 3.9 which were recorded in Athens; the sequence of the Anatolia earthquake of March 18, 1953 (h = 50 km) comprised 21 aftershocks in the first 13 days; the sequence of the Zante earthquake of November 15, 1959 (h = 55 km) consisted of 18 shocks; the focal depths of individual aftershocks, however, were not determined [Papazachos et al. 1967]. In the aftershock sequence of the Kamchatka earthquake of November 4, 1952 the majority of shocks were located near the Mohorovicic discontinuity, but some of them had foci as deep as 60 km [Båth and Benioff 1958]. The sequence following the Aleutian Islands earthquake of March 9, 1957 had an average focal depth of 74 km, and individual shocks were as much as 150 km deep [Duda 1962]. Two more notable exceptions have occurred in Romania and in Central Asia. Iosif and Radu [1961] have studied the aftershock sequence following an earthquake with M = 7.4 and h = 150 km that took place in the region of Vrancea, Romania, on November 10, 1940. The focal depths of aftershocks (3.3 ⩽ M ⩽ 5.5) were of the same order. Lukk [1968] has studied the aftershock sequence of the Dzhurm earthquake of March 14, 1965, which occurred in the Pamir- Hindu Kush region and had a focal depth of 210 km. The observation period lasted for about 22 days, during which 390 aftershocks were recorded; their focal depths increased in time from 200 to more than 240 km. The examples of recorded subcrustal aftershock sequences, however, form a very small part of the total number of sequences known to date; on the other hand, subcrustal earthquakes are themselves much less numerous than crustal ones. Therefore the comparative index of aftershock activity should be given by the ratio of the percentages of crustal and subcrustal earthquakes which are followed by a sequence. At present there is a bias due to instrumentation which favours the detection of shallow aftershock sequences whereas deep ones may go undetected. Thus, the conclusion that aftershocks are generally a shallow phenomenon has to be accepted and at the same time it must be realized that exceptions exist and that the data are far from complete. Time distribution of aftershocks It is customary to regard aftershocks as random events in time, whose frequency is governed by some time-decay law. Jeffreys [1938], in a study of the aftershocks of the Tango, Japan, earthquake of March 7, 1927, found no sign of mutual dependence between aftershocks. That is, there was A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES Table 4. Fluctuations of the observed frequencies in the decay of after- shock activity with time. Table 3. Estimates of the parameters in the time-frequency law. Figure 1. Time distribution of aftershocks: Alaska 1964. 47 no indication that the chance of an aftershock in a given interval of time depended on anything but the time since the main shock, the aftershock frequency falling off with time according to Omori's law. The observed frequency showed only random departures from the law. It has since become a commonly accepted fact that aftershocks can be regarded as random independent events. It follows that any mathematical relationship relating time and frequency must not be interpreted as a physical "law" giving an exact correspondence, but as a statistical law of chance which is followed "on the average", observed frequencies showing random fluctuations from the theoretically expected values. The fact that aftershock sequences consist of independent random events does not imply that they are a simple Poisson process. In a simple Poisson process the probability of occurrence of one event in a given time interval is constant for all t; this is obviously not the case for aftershocks, where the probability of occurrence depends on the time elapsed since the main shock. But, as Jeffreys [1938] and many others have established, apart from the common dependance upon the main shock, no further mutual relation is found within the sequence. In this section the statistical decay law of aftershock activity is estimated for the 15 sequences listed in Table 1. The data have been grouped according to a procedure suggested by Utsu [1962]. The origin time t0 of the main A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES Figure 2. Time distribution of aftershocks: Aleutian Islands 1957. Figure 3. Time distribution of aftershocks: Kern County 1952. Figure 4. Time distribution of aftershocks: San Francisco 1957. Figure 5. Time distribution of aftershocks: Parkfield 1966. 48 shock has been taken as origin of the time axis, t0 = 0. The origin times ti of the aftershocks, obtained from the reference listed in Table 1, have been expressed in terms of days after the main shock. The first day has been excluded from the analysis because of its possible incompleteness with respect to the number of shocks counted, due to the high frequency of aftershocks. Usually, aftershocks occurring in the time interval 1 ⩽ t ⩽ 100 have been considered, unless a sequence comes to an end in a period of time shorter than 100 days. The time axis has been divided into logarithmically uniform intervals, such that their boundaries , satisfy the relation Now, if Ni is the number of aftershocks occurring in the time interval , the quantity represents the observed frequency per unit time interval. This observed frequency is associated with the centered value of the time interval concerned so that one obtains a set of points (ti, ni) in the (t, n)-plane. The data arranged in this fashion are shown in Table 2, the first column representing the centered time, the second the number of shocks in the time interval concerned, and the third the observed frequency. In the sequences (4), (7), (10), (12), (13) and (15), in which the number of shocks in some of the original time intervals was zero, the time intervals have been grouped two by two and ni and ti have been calculated accordingly. The ( ti, ni)-points usually show an approximately linear trend on doubly logarithmic paper. Consequently, it is reasonable to assume that the frequency of aftershocks per unit time n and the time t are related by an equation of the form (1) that is, Omori's law. The commonest procedure for estimating the parameters a and β is the least squares method, which has been applied to the great majority of aftershock sequences whose time distribution has been investigated so far. Accordingly, relation (1) is linearized by taking logarithms on both sides (2) Then, setting the following model is obtained A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES Table 5. Comparison of hyperbolic and exponential decay of aftershock activity. log 0.1 , 0.1, ..., 20.t i ii = = ) n t N i i i= D ) 2t t t 1 i i i= + + ) ) ( )n t t= a b- log ( ) log log .n t t= a b- log ( ) , log , log ,n t y t x= = = =a a b b-) ) ;y x= +a b ) ) ti ) t t t1i i i=D -+ ) ) ) 49 that is, the expected value of y is a linear function of x. Therefore i.e., the observation yi, consists of the expected value at the given xi plus a random fluctuation εi. In the method of least squares, the parameters a* and β* are chosen in such a way that the sum of squares of the vertical distances of the points from the regression line is minimized. For the linear model, this sum of squares is and the necessary conditions for ξ (a*, β*) to be minimum are from which one obtains (3) According to the Gauss-Markov theorem [Graybill 1961], the estimates of the parameters calculated according to the least squares method will be unbiased and of maximum efficiency if, and only if, the linear hypothesis is such that the random fluctuation ε has zero mean and constant variance (independent of x). In other words, if the calculated regression line is to give the expected value of y for each x, the observed values must be uncorrelated, and the probability distribution of y for each x must be symmetric. In many cases, especially when the fluctuation ε can be considered to be the sum of many independent factors, the conditions of the Gauss-Markov theorem are satisfied and the probability distribution of y may be regarded as approximately normal for every x. In other cases, however, the matter is very debatable, particularly when a transformation of coordinates is performed in order to linearize the least square model. It turns out that, under reasonable assumptions, neither the original model (1) nor the linearized model (2) satisfy the conditions of the Gauss-Markov theorem. This has been noted, for instance, by Page [1968b]. An appropriate method must take into account the probability distribution of n for each interval of time considered. On the other hand, it has been suggested [Yule and Kendall 1940] that the least squares method is approximately correct also when the conditions for its theoretical validity are not realized in practice. In order to clarify these matters, in the sequel we shall estimate the parameters in equation (1) by means of both the least squares method for the linearized model (2), and the maximum likelihood method [Cramér 1946], which takes into account the distribution of n for each time interval. First of all, we analyze this distribution. In the following discussion it is assumed that the deviations of the observed values ni from the expected value n(t) reflect actual random fluctuations and that other contributions to such deviations are negligible. Since the only other sources of deviations are errors of measurement, and these have been reduced to a minimum by excluding the first day after the main shock and by counting only welldefined aftershocks with M ⩾ M*, the assumption is most probably A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES Figure 6. Magnitude stability in time: Aleutin Islands. Table 7. Oscillations of mean magnitude. Table 6. Confidence limits on the decay parameter. ,y xi i i= + +a b f ) ) x,( ) y 1 1i k i i k i i 2 2p a b f a bR R= = - - = = ) ) ) )^ h 0 , 0= = da dp db dp ) ) , 1 , 1 k x x k x y x y y x y k y x k x 1 2 1 2 1 1 1 1 1 i k i i k i i k i i i k i i k i i k i i k i = = = = b a b R R R R R R R - - - = = = = = = = ) ) ) ` j Z [ \ ] ]] ] ] 50 correct. It is also assumed that, if the number of shocks expected in the i-th interval is E (Ni), the probability that the observed number is Ni , is given by (4) i.e., the number of shocks in each time interval is given by a Poisson distribution. The Poisson distribution is the most fundamental distribution for such discrete variates as the number of shocks in a given time interval, and it has been postulated for the case of aftershocks by a number of authors, e.g., recently, by Utsu [1962] and Page [1968b]. The expected value E (Ni) in the i-th interval is given by (5) where the approximation is introduced in order to avoid using the integral in equation (4); this is necessary because β is unknown and could be unity. The approximation has been checked numerically for some randomly selected samples and the error was found to be negligible. Then relation (4) becomes (6) The principle of the method of maximum likelihood is to take estimates of the unknown parameters that maximize the probability of obtaining the observed sample. Considering a sample of k independent values, each with a probability distribution p ( Ni; a, β), the probability that the sample consists precisely of these k values is (7) The function L (a, β) is called the likelihood function. The necessary condition for L (a, β) to have a maximum is Since 1n L (a, β) (where 1n stands for the natural logarithm) attains its maximum for the same values of aand β as L (a, β) itself, it is the function 1n L (a, β) which is commonly maximized. It follows from equations (6) and (7) A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES Table 8. Magnitude-frequency distribution (continues on next page). ; !p N E N N E N ei i i i N E N i i= -^^ ^ ^hh h h6 @ ( )E N n t dt t dt t ti t t i t t i 1 1 i i i i = = -a a D )b b- - + + ) ) ) ) ^ h # # ; , !p N N t t e ti i i i ti i=a b a D )b a D - - ) b-^ ^h h ( , ) ( ; , )L p N 1i k i=a b a bP= 0 , 0L L= = da d db d 51 that the likelihood function in the present problem is i.e., from which one obtains The maximum likelihood estimates of a and β are obtained by solving the equations i.e., (8) Equations (8) are the normal equations in the maximum likelihood method. From them one obtains, (9) i.e., (10) Equation (10) must be solved for β, and then a can be obtained from (9). A program has been written to solve equation (10) by the secant method [Ralston 1965]. The first approximation to the root of (10) and to the estimate of ahas been obtained by means of the least squares method for the linearized model, according to formulas (3). The estimates a1, β1 obtained by the least squares method, when compared with the maximum likelihood estimates, give an idea on how statistically reliable the least squares method is when the underlying assumptions are not met. The numerical results for the 15 sequences are summarized in Table 3. From left to right, the columns indicate the sequence involved, the least square estimates a1, β1 and the maximum likelihood estimates a, β. The parameter which characterizes a sequence is the decay parameter β, which measures the rate of decay in time of the frequency of aftershocks. It can be seen that differences between β1 and β are present but not very large. The decay parameter is usually around unity. Now, we proceed to check the validity of the time- frequency law. If the expected value of the frequency n varies in time according to equation (1); and using the approximation expressed by (5), the mean and the variance of the Poisson- distributed number of shocks are in each interval and therefore the standard deviation of n is A general theorem which holds for an arbitrary distribution with a second moment is Tchebychev's theorem [Cramér 1946]. It states that, if X is a random variable with mean E (X) and standard deviation D ( X), then the following inequality holds A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES ( , ) !L N t t e 1 t i k i i i ti i=a b a P D = )b a D - - ) b-^ h 1 ( , ) 1 !n L n N t t e 1 t i k i i i ti i=a b a R D = )b a D - - ) b-^ h; E 1 ( , ) 1 1 1 1 ! n L n N N n t N n t t t n N 1 1 1 1 1 i k i i k i i k i i i k i i i k i i=a b a b a R R R D R D R - + - - - = = = = = ) )b- 1 ( , ) 1 0 1 ( , ) 1 1 0 n L N t t n L N n t t t n t 1 1 1 1 i k i i k i i i k i i k i ii i = = = = da d a b a da d a b a R R D R R D - - + = = = = ) ) b b - - * 0 1 1 0 N t t N n t t t n t 1 1 1 1 i k i i k i i i k i i k i ii i = = a a R R D R R D - - = = = = ) ) b b - - * 1 1 t t N t t n t N n t 1 1 1 1 i k i i i k i i k i i i k i i i = =a R D R R D R = = = = ) )b b- - ( ) 1 1 0 F N t t n t N n t t t 1 1 1 1 i k i i k i i i k i i k i i i i = = b R R D R R D - - = = = = ) ) b b - - ( )E N N n t t t tV= = = aD D) )b-^ ^h h t t .D n = a D ) b-^ h where k is an arbitrary positive number. In other words, the probability that X assumes values outside the interval E (X) ± k D ( X) is less than 1/k². Conversely, from the viewpoint of sampling, in the long run less than (100/k²)% of the values assumed by X will fall outside the interval. The application of Tchebychev's inequality to check the validity of the time-frequency law has two limits, namely, it yields rather weak conditions, and the sample size is small. Nevertheless it gives at reasonably safe criterion for rejecting the validity of the law for sequences that show too wide fluctuations. Regarding the observed values of n as the result of random sampling from a population whose expected value varies with time according to (1), and choosing k = 2, Tchebychev's inequality takes the form that is, in the long run less than 25% of the observed values ni should fall outside the interval E (n) ± 2 D (n). If this condition is not satisfied, the assumption regarding the variation with time of the expected value n(t) must be rejected. In order to determine a confidence band according to Tchebychev's inequality, therefore, the quantities E (n) + 2 D (n), E (n) – 2 D (n) have been computed at selected ti, i = 1, 2, ..., k; and compared with the observed ni. Table 4 gives the results, showing in the columns from left to right the sequence, the total number of data points, the number of points outside the confidence band, and their percentage. Consequently, according to the selected criterion, sequenze (10), (12), (13) and (14) do not follow the assumed timefrequency law, inasmuch as the observed frequency cannot be explained only in terms of random fluctuations from the law. The other 11 sequences appear to follow the law and the fit is generally fairly good. .All the sequences which show considerable departures from the assumed law have occurred in the region of Greece; this fact might have some geotectonic significance. However, there exists the possibility that relatively poor instrumentation plays a part in some apparent irregularities. Sequence (14) originated under peculiar conditions; the frequency of aftershocks in it appears to be correlated to the variations in the water loading of a nearby artificial lake [Comninakis et al. 1968]. Figures 1 to 5 display on doubly logarithmic paper the results for some of the 11 sequences which appear to decay according to the postulated law. The dots represent the data points, the full line the fitted n(t), and the broken lines the confidence limits. It is interesting to note that most decay phenomena in physics are exponential, whereas the decay of aftershock activity with time appears to be hyperbolic. In order to check rapidly the possibility of exponential decay, the model has been linearized, and the parameters a2 and β2 have been estimated by the least squares method. Then the mean square deviation (standard error of estimate) of the data points from the fitted curve has been calculated by the formula (where k is the number of data points) and compared with A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES 52 Table 9. Estimation of the parameter b in the magnitude-frequency law. Table 10. Validity of the magnitude-frequency law. 1P X E X k D X k2; ; H G- ^ ^h h" , 2 4 1P n t t t ; ; H Ga a D - ) b b - -' 1 ( )n t e2 t2a= b- 2 1 log log logs k n t e 2 1 2e i k i i2 2= -a bR - - = l ^ h6 @ 53 the mean square deviation from the fitted time-frequency law that is, (where a1 and β1 are the least square estimates). If the assumed time-frequency law represents a better fit, in the least square sense, than the exponential decay. The results are summarized in Table 5. It can be seen that in all but three cases, namely, sequences (7), (10), and (15). It can be concluded that the frequency of aftershocks usually decreases in time hyperbolically and not exponentially. Confidence limits on the decay parameter The parameter β appearing in equation (1) is related to the rate of decay of aftershock activity with time and is therefore an important characteristic of the sequence under consideration. Consequently, it would be interesting to see whether the differences in the computed decay parameters for different sequences are significant or not. Confidence limits on β, however, cannot be calculated according to the usual least square procedure, which assumes the fluctuations to have a normal distribution. An approximate procedure for calculating confidence limits on β can be based on the addition theorem for the Poisson distribution [Cramér 1946]. Since in this procedure the parameter a is assumed to be known exactly, it will yield no more than an indication of the range in which the real value of β is likely to fall. The addition theorem for the Poisson distribution states that, if Ni, i = 1, 2, ..., k, is a sequence of stochastically independent and Poisson-distributed random variables with expected values Ei = E (Ni), then the sum will be Poisson-distributed with expected value In the present case, Ni being the number of shocks in the i-th time interval, the expected values are and The value of E (X) for each of the 11 sequences which follow the time-frequency law is such that the distribution of X can be approximated by the normal distribution. Therefore the standardized variable is approximately normally distributed with E (X*) = 0, V (X*) = 1. Therefore, the probability that X* assumes a value in the interval (λ1, λ2) is where Φ (X*) is the normal distribution function. In particolar Then the approximate 95% confidence limits on the decay parameter are obtained by solving for β the equations that is, recalling the definition of X*, A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES Figure 7. Magnitude-frequency distribution: Aleutian Islands 1957. Figure 8. Magnitude-frequency distribution: Kern County 1952. ( )n t t1 1a= b- 2 1 log log logs k n t t 2 1 1e i k i i i1 2= a bR - - - = ^ h6 @ E t ti i i= a D )b- E X t t . 1i k i i= a R D= )b-^ h X D X X E X = -) ^ ^ h h ( ) ( )P X1 2 2 1=G Gm m m mU U- )" , 1.96 1.96 .9 .P X 0 5=G G- )" , 1.96X = !) s s2 2e e1 l s s2 2e e1 l X N 1i k i= R= E X E 1i k i= R=^ h Naturally, the central value of X* is obtained when a, β are given by the maximum likelihood estimates. Assuming the value of a to be known exactly, the above equations take the form (11) Equations (11) have been solved by the secant method. The results for the 11 sequences with time-frequency law of the form given by (1) are shown in Table 6. It can be seen that the 95% confidence limits on the maximum likelihood estimate always contain the least square estimate of β, which therefore appears to be a good approximation to the real value of the decay parameter. The values listed in Table 6 must not be considered as exact. The limits rounded to the second decimal digit are probably fairly reliable. Most β-values cluster around the 0.9- 1.2 range. Sequence (11) stands by itself, showing a very rapid decay in activity. Also sequence (9) has a decay coefficient somewhat higher than usual. This fact hints at the possibility that sequences in Greece decay more rapidly; but no conclusion can be reached with such a small sample size. Magnitude stability in time When the characteristics of aftershock sequences with respect to magnitude are examined, two sources of error are added to the possibility of the incompleteness of data; namely, lack of accuracy in magnitude determination, and confusion between different magnitude scales. Unfortunately, authors sometimes do not specify which scale they are using. When one single sequence is being examined, no problems arise, because the data are consistent within the sequence; if, however, results for different sequences are to be compared, the use of different magnitude scales may affect the conclusions. In the sequel, the various "local" magnitude scales, for the purposes of comparison of results among sequences, have been assimilated to M. Thus the only distinction left is between M and m; the latter scale has been used only in sequence (1). Two aspects of the sequences have been examined in detail, namely, the variation of aftershock magnitudes with time, and the magnitude- frequency distribution. For all sequences except (1), where 1 ⩽ t ⩽100 days, the data for t ⩽100 days have been included in the analysis. We first consider the distribution of magnitude with respect to time. The overall mean magnitude, M, has been calculated for each sequences as where K is the total number of aftershocks in the sequence. Then the mean magnitude, M´, of each group of 10 successive aftershocks is computed, thereby eliminating large individual fluctuations. In almost all the sequences considered the mean magnitude M´ oscillates about during the whole length of the sequence and no appreciable decrease with time is detectable for t ⩽ 100 days. When a sequence lasts less than 100 days, the mean magnitude is stable throughout the sequence. Sometimes higher values of M´ are observed in the first few hours after the main shock, A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES 54 Figure 9. Magnitude-frequency distribution: San Francisco 1957. Figure 10. Magnitude-frequency distribution: Parkfield 1966. 1.96 t t N t t 1 1 1 i k i i i k i i k i i = ! a a R D R R D- = = = ) ) b b - - 1.96 0 1.96 0 . N t t t t N t t t t 1 1 1 1 1 1 i k i i k i i i k i i i k i i k i i i k i i = = a a a a R R D R D R R D R D - - - + = = = = = = ) ) ) ) b b b b - - - - * 1 , 1, 2, ...,M K M i K1i k i= =R= M 55 but this can almost certainly be attributed to the fact that some shocks with M ⩾ M* are not detected when the frequency is very large. Figure 6 gives an example. The dotted line represents . Thus, the law of magnitude stability in aftershock sequences (first proposed by Lomnitz 1966) is fully confirmed; during an aftershock sequence the mean magnitude of the shocks is constant in time. Consequently, the decrease of seismic activity with time is solely due to the decrease in aftershock frequency. Moreover, the fluctuations of M´ about the overall mean magnitude are not very wide. Table 7 summarizes the results for the sequences with a larger number of shocks. From left to right, the first column indicates the sequence, the second the overall mean magnitude, the third the number of calculated M í -points, and the fourth the percentage of such points which fall within the interval ± 0.20. It is to be noted that also the sequences with a lesser numher of aftershocks and not included in the table show remarkable magnitude stability. The only sequence in which M´ shows a decreasing trend with time is sequence (11). In all the others, M´ shows only random fluctuations from and the law of magnitude stability in time is satisfied. Magnitude-frequency distribution The frequency data for the 15 sequences of Table 1 are listed in Table 8. The first column gives the centered value of magnitude in the interval concerned (M ± 0.05); n(M) is the frequency; and N (M) the cumulative frequency. For brevity, intervals in which the frequency was zero have been omitted from the table. The magnitudes in sequence (1) are in the m-scale. The most commonly accepted form for the magnitude-frequency distribution, in case of both independent seismic events and aftershock sequences, is (12) where log is the logarithm to the base 10 and n(M) is the number of shocks with M ± dM [Gutenberg and Richter 1954]. Equation (12) is to be regarded as expressing a statistical relationship. Usually, the coefficients a and b have been calculated according to the least squares method. Suzuki [1958] has pointed out that this is not rigorous, because log n(M) is not symmetrically distributed with uniform variance for each magnitude interval ΔM. In this section we shall define n(M) is such a way that log n(M) = anot when M= 0, as in equation (12), but when M= M*, where M* is the minimum detectable magnitude in the sequence. Then the magnitude-frequency law takes the form (13) Converting to natural logarithms one obtains (14) where (15) Therefore Normalizing, i.e., imposing the condition that one has i.e., γ = b´. Thus we assume that the probability distribution of M takes the form (16) A procedure for estimating the parameter b´ in (16) can be derived as follows. The mean of the distribution is Approximating the population mean by the sample mean , given by where K is the total number of shocks, one has A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES Figure 11. Magnitude-frequency distribution: Cremasta 1966. log ( )n M a b M= - log ( ) ( ) ,n M a b M M M M M .H= - - ) ) 1n ( ) ( )n M a b M M M= - - )l l log , loga e a b e b = =l l ( )n M e e e , e .( ) ( )M M b M M aa b c c= = =- - - - ) )l l l l ( ) 1n M dM M = 3 ) # 1e dM b ( )b M M M c c = = 3 - - ) ) l l# ( ) ,n M b e M M .( )b M M H= )- - )l l 1M K M1i k i= R= ( ) 1E M M n M dM M b . M = = +) 3 ) l ^ h # M M M M M i.e., (17) The above procedure is but a particular case of the time- honored method of moments, introduced by K. Pearson and his school [Cramér 1946]. Formula (17) was also proposed by Utsu at a meeting of the Seismological Society of Japan (reported by Aki 1965). The estimation of the parameter b´ given by (17) has been proved by Aki [1965] to be equivalent to the maximum likelihood estimate, and therefore has several desirable statistical properties. Table 9 gives the results of calculations. From left to right, the columns indicate sequence, sample mean , minimum magnitude M*, and the estimations of the parameters b´, b appearing in relations (16) and (13) respectively . The estimate of b´ is given by (17), and b is given by (15). The minimum magnitude M* has been taken to be 0.05 units less than the values given in Table 8 because the value of M approximated to one decimal could actually come from anywhere in the interval M ± 0.05. It can be seen that b is usually slightly less than unity. The basic idea for examining the observed fluctuations from the magnitude-frequency law is as in the case of the time distribution of aftershocks. Suzuki [1958], among others, has argued that the number of shocks with M ± dM must follow a Poisson distribution. Accordingly, the problem is that of checking whether the observed deviation can be explained in terms of random fluctuations from the law. It is more convenient to consider the cumulative distribution because individual large fluctuations in small intervals, possibly due to inaccurate magnitude determination, are smoothed out in this fashion, and moreover the normal approximation can be employed. The basic idea of the method is due to Suzuki [1958]. Some minor modifications have been introduced, and the normal approximation to the Poisson distribution, instead of the Poisson distribution itself, has been used. From equation (16), it follows that the cumulative distribution function of magnitude has the form Assuming that the total number of aftershocks in a sequence, K, coincides with the theoretical value for M = M*, the expected value of the cumulative frequency at various M > M* is given by which, when transformed by taking logarithms on both sides, is a straight line on semi-logarithmic paper, with slope equal to –b if the logarithms are to the base 10. Now, to each magnitude range there corresponds a Poisson-distributed number of shocks. The cumulative number of shocks at a given magnitude M, therefore, is the summation of independent samples taken from each of the Poisson distributions corresponding to magnitudes greater or equal to M. According to the addition theorem for the Poisson distribution, such a cumulative number will also be Poisson-distributed. If the expected value is large enough, say, N (M) ⩾10 for all intervals, the Poisson distribution can be approximated by the normal distribution with mean N(M) and standard deviation . It is then possible to calculate the fiducial interval beyond which fluctuations are expected with a probability smaller than 5% (2.5% on each side). The limits of the interval such that are, for the normal distribution, ; i.e., in the present case, . Accordingly, the above quantities have been calculated at Then, by joining all points of ordinates respectively, one obtains a confidence band which should contain approximately 95% of the data points if they come from a population whose expected value is given by the magnitude-frequency law as expressed by (18). The rightmost interval (say, from Mj to ∞) has always be chosen in such a way that N ( M) ⩾ 10. Table 10 gives the results. It follows that sequences (1), (3), (9), (10), (12) and (13) do not appear to be governed by the assumed magnitude-frequency law. When these sequences are examined one by one, however, it is seen that several circumstances tend to decrease the weight that must A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES 56 Table 11. Confidence limits on the frequency parameter. 1M M b = +) l 1b M M = - ) l ( )M n M dM eF ( )b M M M = = 3 - - )l^ h # M K eN ( )b M M= - - )l^ h 0.95 ,P N1 2G G -m m" , 0.1 , 1, 2, ..., .M M i i Ki = =+ ) 2 and 2 ,M N M N M N MN i i i i+ -^ ^ ^ ^h h h h M N M^ h 2 NE N D!^ ^h h 2N M N M!^ ^h h 57 be assigned to their apparent irregularity. For sequence (1) the body-wave magnitudes as given by the U.S. Coast and Geodetic Survey have been employed. Page [1968b] found them inaccurate and modified them on the basis of records of five selected stations, thereby obtaining a better fit; however, no test on the fluctuations was performed in his study. Of the sequences occurring in Greece, it is worth noting that sequence (14), which was studied with a network improved with respect to the others, appears to follow the assumed law; this fact supports indirectly the conclusion that the irregular behavior of its aftershock frequency is due to changing local stress conditions. Sequences (10), (12) and (13) are also irregular with respect to aftershock frequency in time. This fact points to the likelihood that such irregularities are due to poor data, although no definite statement can be made. It can therefore be concluded that the magnitude- frequency law (18) is statistically followed by the large majority of the sequenze under consideration. Figures 7 to 11 give some examples; the logarithmic ordinate represent the cumulative frequency, the linear abscissa the magnitude. Confidence limits on the frequency parameter Differences between b-values may be significant or not. In order to decide the question, confidence limits must be set on the estimated value of the parameter. We shall follow in this matter a procedure suggested by Aki [1965]. Given a sample of K shocks with magnitudes Mi, i = 1, ..., K, let yi and Y be defined by where n(M) is given by formula (16). Clearly for all i, and therefore By the central limit theorem [Cramér 1946] the distribution of Y will be approximately normal if K is sufficiently large. Since it follows that the variable is approximately normally distributed with mean 0 and standard deviation 1; therefore Accordingly, the 95% confidence limits on b´ are obtained by solving the inequality which gives i.e., Table 11 summarizes the results for the 9 sequences where the assumed magnitude-frequency law appears to hold. According to usage, the parameter b, instead of b´, has been employed. It can be seen that the b-values cluster around the interval 0.8-1.0. Although sequence (2), which occurred along an active Island arc, shows an anomalously high value of the frequency parameter, the sample size is too small to support the contention that b has some geotectonic significance. Conclusion The basic statistical laws describing the phenomenology of aftershock sequences are confirmed by the present study. These laws are as follows: (1) Aftershock sequences are generally crustal events, although deeper ones also occur; (2) The frequency of aftershock occurrence within the same sequence decays in time according to the law where the decay parameter β is approximately equal to, or slightly greater than, unity; (3) The aftershock magnitudes, apart from individual fluctuations, show stability in time to the end of the sequence; (4) The frequency-distribution of magnitude in a sequence is of exponential form where the frequency parameter b = b´ log e is usually slightly less than unity The importance of the mode of grouping the data in a statistical analysis makes it desirable to introduce a standardized procedure. Furthermore, an appropriate statistical method must be employed, and the hypotheses A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES 1n ( ) b n Myi i= d d l Y y 1i k i= R= 1y b M M= - + ) l ( ) 0E y n M dMy M = = 3 ) ^ h # ( ) 1V y y n M dM b 2 2 M = = 3 ) l ^ h # 0 , V Y b KE Y 2= = l ^ ^h h Y Y E Y K b YY D == -) l^ ^ h h 1.96 1.96 0.95 .P Y =G G- )" , 1.96 1.96 K b Y G G- l 1.96 1 1.96 K b Y b M M 1i k iG GR- - -= )l l ` j 1 1.96 / 1 1.96 / M M K b M M K G G - - - + ) ) l ( ) , 100n t tt Ha= b- M K eN ( )b M M= - - )l^ h involved must be checked. A systematic treatment of data greatly increases the reliability of the results. It also turns out that the least squares method, when the observed values are uncorrelated, yields rather satisfactory results even if the conditions for its theoretical validity are not met. Acknowledgements. The author wishes to express his thanks to Professor A. E. Scheidegger of the University of Illinois at Urbana, who kindly discussed some questions which arose during the preparation of this paper and reviewed the manuscript. Numerical calculations were performed on the University of Illinois IBM System/360 computer, operating under a grant from the United States National Science Foundation. The support received from this agency is gratefully acknowledged. References Aki, K. (1965). Maximum Likelihood Estimate of b in the Formula log N = a – bM and its Confidence Limits, Bull. Earthq. Res. Inst., 43, 237-239. Båth, M. and H. Benioff (1958). The Aftershock Sequence of the Kamchatka Earthquake of November 4, 1952, Bull. Seism. Soc. Am., 48, 1-15. Benioff, H. (1951). Earthquakes and Rock Creep, I, Bull. Seism. Soc. Am., 41, 31-62. Comninakis, P., J. Drakopoulos, G. Moumoulidis and B. Papazachos (1968). Foreshock and Aftershock Sequences of the Cremasta Earthquake and Their Relation to the Waterloading of the Cremasta Artificial Lake, Ann. Geofisica, 21, 39-71. Cramér, H. (1946). Mathematical Methods of Statistics, Princeton. Duda, S. J. (1962). Phänomenologische Untersuchung einer Nachbebenserie aus dern Gebiet der Alëuteninseln, Freib. Forschungshefte, C-132, 5-90. Graybill, F. A. (1961). An Introduction to Linear Statistical Models, I, New York. Gutenberg, B. and C. F. Richter (1954). Seismicity of the Earth and Associated Phenomena, 2nd ed., Princeton. Iosif, T. and C. Radu (1961). Elastic Strain Characteristics in the Deep-Focus Earthquakes of the Vrancea Region, Stud. Cerc. Astr. Seism, 2, 269-285 (in Romanian). Jeffreys, Sir H. (1938). Aftershocks and Periodicity in Earthquakes, Gerl. Beitr Geophys., 53, 111-139. Lomnitz, C. (1966). Magnitude Stability in Earthquake Sequences, Bull. Seism. Soc. Am., 56, 247-249. Lukk, A. A. (1968). The Aftershock Sequence of the Dzhurm Deep-Focus Earthquake of March 14, 1965, Izv. Earth. Phys., 5, 83-85 (in Russian). Matuzawa, T. (1964). Study of Earthquakes, Uno Shoten, Tokio. McEvilly, T. V., W. H. Bakun and K. B. Casaday (1967). The Parlefield, California, Earthquakes of 1966, Bull. Seism. Soc. Am., 57, 1221-1244. Mogi, K. (1962). On the Time Distribution of Aftershocks Accompanying the Recent Major Earthquakes in and near Japan, Bull. Earthq. Res. Inst., 40, 107-124. Page, R. (1968a). Focal Depths of Aftershocks, J. Geophys. Res., 73, 3897-3903. Page, R. (1968b). Aftershocks and Microaftershocks of the Great Alaska Earthquake of 1964, Bull. Seism. Soc. Am., 58, 1131-1168. Papazachos, B., N. Delibasis, N. Liapis, G. Moumoulidis and G. Purcaru (1967). Aftershock Sequences of Some Large Earthquakes in the Region of Greece, Ann. Geofisica, 20, 1-93. Ralston, A. (1965). A First Course in Numerical Analysis, New York. Richter, C. F. (1955). Foreshocks and Aftershocks, in Earthquakes in Kern County, California, during 1952, Bull. Cal. Div. Mines., 171, 177-197. Richter, C. F., C. R. Allen and J. M. Nordquist (1958). The Desert Hot Springs Earthquakes and Their Tectonic Environment, Bull. Seism. Soc. Am., 48, 315-337. Suzuki, Z. (1958). A Statistical Study of the Occurrence of Small Earthquakes, III, Sci. Rep. Tohoku Univ., 5th Series, Geophys., 10, 15-27. Tocher, D. (1959). Seismographic Results from the 1957 San Francisco Earthquake, in San Francisco Earthquakes of March 1957, Cal. Div. Mines Sp. Rep., 57, 59-71. Udias, A. (1965). A Study of the Aftershocks and Focal Mechanism of the Salinas-Watsonville Earthquakes of August 31 and Septernber 14, 1963, Bull. Seism. Soc. Am., 55, 85-106. Utsu, T. (1962). On the Nature of Three Alaskan Aftershock Sequences of 1957 and 1958, Bull. Seism. Soc. Am., 52, 279-297. Yule, G. U. and M. G Kendall. (1940). An Introduction to the Theory of Statistics, 12th ed., London. © 2010 by the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved. A STATISTICAL STUDY OF AFTERSHOCK SEQUENCES 58