A review and new insights on the estimation of the b-value and its uncertainty

The estimation of the b-value of the Gutenberg-Richter Law and its uncertainty is crucial in seismic hazard studies, as well as in verifying theoretical assertions, such as, for example, the universality of the Gutenberg-Richter Law. In spite of the importance of this issue, many scientific papers still adopt formulas that lead to different estimations. The aim of this paper is to review the main concepts relative to the estimation of the b-value and its uncertainty, and to provide some new analytical and numerical insights on the biases introduced by the unavoidable use of binned magnitudes, and by the measurement errors on the magnitude. It is remarked that, although corrections for binned magnitudes were suggested in the past, they are still very often neglected in the estimation of the b-value, implicitly by assuming that the magnitude is a continuous random variable. In particular, we show that: i) the assumption of continuous magnitude can lead to strong bias in the b-value estimation, and to a significant underestimation of its uncertainty, also for binning of ∆M = 0.1; ii) a simple correction applied to the continuous formula causes a drastic reduction of both biases; iii) very simple formulas, until now mostly ignored, provide estimations without significant biases; iv) the effect on the bias due to the measurement errors is negligible compared to the use of binned magnitudes. Mailing address: Dr. Warner Marzocchi, Istituto Nazionale di Geofisica e Vulcanologia, Via D. Creti 12, 40128 Bologna, Italy; e-mail: marzocchi@bo.ingv.it 1272 Warner Marzocchi and Laura Sandri In this frame, the fractal geometry distribution and the earthquake dynamics are the spatial and temporal signatures of the same phenomenon. On the other hand, several lines of empirical evidence dispute the universality of the GR Law. In fact, in spite of a «first order» validity of the GR Law with a constant b-value ≈ 1, significant spatial and temporal variations in the bvalue seem to take place (Miyamura, 1962; Healy et al., 1968; Pacheco and Sykes, 1992; Kagan, 1997; Wiemer and McNutt, 1997; Wiemer and Wyss, 1997; Wiemer et al., 1998; Wiemer and Katsumata, 1999). These variations are usually attributed to many processes, such as fault heterogeneity (Mogi, 1962), the stress level imposed on rocks (Scholz, 1968), and pore pressure variations. Another important discrepancy of the hypothesis of the universality of the GR Law is that the latter seems to hold only for a finite range of magnitudes (e.g., Pacheco and Sykes, 1992; Pacheco et al., 1992; Scholz, 1997; Triep and Sykes, 1997; Knopoff, 2000). For example, Kagan (1993, 1994) found that a gamma distribution better matches the requirement of a maximum seismic energy released by an earthquake (cf. Wyss, 1973). In all the cases reported above the b-value should not be a constant in the earthquake catalog. The empirical validation of the «universality» hypothesis, as well as the identification of spatial and temporal changes, passes through the estimation of the b-value and its uncertainty. Different formulas were proposed in the past (e.g., Aki, 1965; Utsu, 1965, 1966; Shi and Bolt, 1982; Bender, 1983; Tinti and Mulargia, 1987) which take into account the unavoidable binning of the magnitudes in different ways. As a matter of fact, Aki’s formulas, which consider magnitude a continuous random variable, are still the most used. Implicitly, this means that, at least for instrumental measurements, the effect of the binning is considered negligible, and that all the formulas mentioned above give almost comparable estimations. In this paper, we review the most diffuse formulas so far proposed and investigate analytically and numerically on the main properties of such estimations. In particular, we estimate the reliability and the possible bias of the formulas as a function of the number of data and of the measurement errors in the magnitude. Here, we do not consider the estimations provided through the least squares technique. In fact, even though this approach was still employed in many recent scientific papers (Pacheco and Sykes, 1992; Pacheco et al., 1992; Karnik and Klima, 1993; Okal and Kirby, 1995; Scholz, 1997; Triep and Sykes, 1997; Main, 2000), the use of the least squares technique does not have any statistical foundation (e.g., Page, 1968; Bender, 1983). 2. Estimation of b and b v t through the Maximum Likelihood technique In the first papers that described the Maximum Likelihood (ML) estimation (Aki, 1965; Utsu, 1965), the magnitude M was considered a continuous Random Variable (RV). If eq. (1.1) holds, the probability density function (pdf) of M is ln f M b 10 10 10 10 bM bM bM min max = ^ ^ h h (2.1) where Mmin and Mmax are, respectively, the minimum and the maximum magnitude allowed. If Mmax >> Mmin, eq. (2.1) becomes M 10 b M M min . ln f b 10 = ^ ^ ] h h g (2.2) Note that the passage from eq. (2.1) to eq. (2.2) requires that, in practice, the GR Law holds for a range of magnitudes Mmax – Mmin ≥ 3. This assumption might be questionable in studying variations of the b-value in smaller ranges (e.g., Knopoff, 2000). The ML estimation of eq. (2.2) consists of choosing the b-value which maximizes the likelihood function (Fisher, 1950), that is ln b M 10 1 thresh = n t t ^ ` h j (2.3) where nt is the sampling average of the magnitudes, and Mthresh is the threshold magnitude 1273 A review and new insights on the estimation of the b-value and its uncertainty which usually corresponds to the minimum magnitude for the completeness of the seismic catalog. Hereinafter, the symbol ^ distinguishes the estimation by the true value of the parameter. We have also attached an asterisk to b in order to distinguish it from another estimation which we shall introduce later. The uncertainty is estimated by (Aki, 1965) N b b = v ) ) t t t (2.4) where N is the number of earthquakes. A major contribution was provided later on by Shi and Bolt (1982) who provided a new formula to estimate the error of the b-value . b N N M 2 30 1 ( ) b SB i i N


Introduction
The Gutenberg-Richter Law (from now on GR Law) (Gutenberg and Richter, 1954) is certainly one of the most remarkable and ubiquitous features of worldwide seismicity.In the most common form it reads where N is the number of events with magnitude M, and a and b are two constant coefficients.
The scientific importance of the GR Law is linked to its apparent ubiquity.It emerges in a variety of tectonic settings and depth ranges, in seismic catalogs ranging from a few months to centuries, and in natural as well as induced seismicity.This ubiquity makes the GR law very useful in seismic hazard studies, and very attractive from a theoretical point of view.In past years, many efforts were devoted to understanding the physical meaning of such a power law, but the conclusions are diverging.Some authors (e.g., Bak and Tang, 1989;Ito and Matsuzaki, 1990), by assuming a presumed universality of the parameter b and the implicit scale independence in eq.(1.1), propose a model of Self-Organized Criticality (SOC) to explain earthquake occurrence.

A review and new insights on the estimation of the b-value and its uncertainty
Warner Marzocchi and Laura Sandri Istituto Nazionale di Geofisica e Vulcanologia, Bologna, Italy

Abstract
The estimation of the b-value of the Gutenberg-Richter Law and its uncertainty is crucial in seismic hazard studies, as well as in verifying theoretical assertions, such as, for example, the universality of the Gutenberg-Richter Law.In spite of the importance of this issue, many scientific papers still adopt formulas that lead to different estimations.The aim of this paper is to review the main concepts relative to the estimation of the b-value and its uncertainty, and to provide some new analytical and numerical insights on the biases introduced by the unavoidable use of binned magnitudes, and by the measurement errors on the magnitude.It is remarked that, although corrections for binned magnitudes were suggested in the past, they are still very often neglected in the estimation of the b-value, implicitly by assuming that the magnitude is a continuous random variable.In particular, we show that: i) the assumption of continuous magnitude can lead to strong bias in the b-value estimation, and to a significant underestimation of its uncertainty, also for binning of ∆M = 0.1; ii) a simple correction applied to the continuous formula causes a drastic reduction of both biases; iii) very simple formulas, until now mostly ignored, provide estimations without significant biases; iv) the effect on the bias due to the measurement errors is negligible compared to the use of binned magnitudes.
In this frame, the fractal geometry distribution and the earthquake dynamics are the spatial and temporal signatures of the same phenomenon.On the other hand, several lines of empirical evidence dispute the universality of the GR Law.In fact, in spite of a «first order» validity of the GR Law with a constant b-value ≈ 1, significant spatial and temporal variations in the bvalue seem to take place (Miyamura, 1962;Healy et al., 1968;Pacheco and Sykes, 1992;Kagan, 1997;Wiemer and McNutt, 1997;Wiemer and Wyss, 1997;Wiemer et al., 1998;Wiemer and Katsumata, 1999).These variations are usually attributed to many processes, such as fault heterogeneity (Mogi, 1962), the stress level imposed on rocks (Scholz, 1968), and pore pressure variations.Another important discrepancy of the hypothesis of the universality of the GR Law is that the latter seems to hold only for a finite range of magnitudes (e.g., Pacheco and Sykes, 1992;Pacheco et al., 1992;Scholz, 1997;Triep and Sykes, 1997;Knopoff, 2000).For example, Kagan (1993Kagan ( , 1994) ) found that a gamma distribution better matches the requirement of a maximum seismic energy released by an earthquake (cf.Wyss, 1973).In all the cases reported above the b-value should not be a constant in the earthquake catalog.
The empirical validation of the «universality» hypothesis, as well as the identification of spatial and temporal changes, passes through the estimation of the b-value and its uncertainty.Different formulas were proposed in the past (e.g., Aki, 1965;Utsu, 1965Utsu, , 1966;;Shi and Bolt, 1982;Bender, 1983;Tinti and Mulargia, 1987) which take into account the unavoidable binning of the magnitudes in different ways.As a matter of fact, Aki's formulas, which consider magnitude a continuous random variable, are still the most used.Implicitly, this means that, at least for instrumental measurements, the effect of the binning is considered negligible, and that all the formulas mentioned above give almost comparable estimations.In this paper, we review the most diffuse formulas so far proposed and investigate analytically and numerically on the main properties of such estimations.In particular, we estimate the reliability and the possible bias of the formulas as a function of the number of data and of the measurement errors in the magnitude.
Here, we do not consider the estimations provided through the least squares technique.In fact, even though this approach was still employed in many recent scientific papers (Pacheco and Sykes, 1992;Pacheco et al., 1992;Karnik and Klima, 1993;Okal and Kirby, 1995;Scholz, 1997;Triep and Sykes, 1997;Main, 2000), the use of the least squares technique does not have any statistical foundation (e.g., Page, 1968;Bender, 1983).

Estimation of b and b v t through the Maximum Likelihood technique
In the first papers that described the Maximum Likelihood (ML) estimation (Aki, 1965;Utsu, 1965), the magnitude M was considered a continuous Random Variable (RV).If eq.(1.1) holds, the probability density function (pdf) of M is where Mmin and Mmax are, respectively, the minimum and the maximum magnitude allowed.If M max >> Mmin, eq.(2.1) becomes Note that the passage from eq. (2.1) to eq. (2.2) requires that, in practice, the GR Law holds for a range of magnitudes Mmax -Mmin ≥ 3.This assumption might be questionable in studying variations of the b-value in smaller ranges (e.g., Knopoff, 2000).
The ML estimation of eq. ( 2.2) consists of choosing the b-value which maximizes the likelihood function (Fisher, 1950), that is where n t is the sampling average of the magnitudes, and Mthresh is the threshold magnitude which usually corresponds to the minimum magnitude for the completeness of the seismic catalog.Hereinafter, the symbol ^ distinguishes the estimation by the true value of the parameter.We have also attached an asterisk to b * in order to distinguish it from another estimation which we shall introduce later.The uncertainty is estimated by (Aki, 1965) where N is the number of earthquakes.
A major contribution was provided later on by Shi and Bolt (1982) who provided a new formula to estimate the error of the b-value where N is the number of earthquakes.Compared to eq. ( 2.4), eq. ( 2.5) provides a reliable estimation also in the presence of possible (time and/or spatial) variations of the b-value (Shi and Bolt, 1982).
As a matter of fact, the magnitude is not a continuous variable and it is not devoid of measurement errors.In practical cases, the uncertainties on the measured magnitudes lead to the use of «binned» magnitudes, i.e. the magnitudes are grouped by using a selected interval ∆M.For instance, for instrumental measurements, the magnitude interval used for the grouping is ∆M = 0.1; for magnitude estimation of historical events, the grouping can even be ∆M ≥ 0.5.
In the following sections, we consider separately the effects due to the binning of the magnitudes and the measurement errors on the estimations of the b-value and its uncertainty.We remark that a wrong choice of the completeness threshold magnitude Mthresh can also cause a significant bias, as shown by Wiemer and Wyss (2000).Here, we do not deal with this issue, but it should be kept in mind in any practical analysis.

Influence of the binned magnitudes
In this section, we first study the biases introduced by the binning of the magnitudes in the estimation of the b-value and its uncertainty made through eqs.(2.3) and (2.4).Then, we check through numerical simulations the reliability of some formulas reported in scientific literature which either assume the magnitudes as a continuous RV or take properly into account the binning of the magnitudes.

Effects of the binned magnitudes on Aki's (continuous) formulas
The use of eq. ( 2.3) produces a biased estimation of the real b-value, because of two factors: i) the average µ of a continuous RV with a power law distribution is different from the average of the same «binned» RV; ii) Mthresh ≠ Mmin.Hereinafter we will call these biases θ1, and θ2, respectively.
The influence of θ1 was estimated by Bender (1983).She found that the sample average n t computed from binned data is systematically higher than the true value µ.This is due to the fact that the real (continuous) magnitudes in the interval M i − ∆Μ/2 ≤ M < Mi + ∆Μ/2 are not symmetrically distributed around the central value Mi.Moreover, she showed that the bias θ 1 is negligible for ∆M = 0.1, while it is very important for larger ∆M (for example, ∆M = 0.6), that might be necessary to evaluate the b-value for historical catalogs.
In any case, we argue that this correction could have not been largely employed also because the influence of the bias θ 2 is erroneously considered negligible.Indeed, the difference between b ) t and b t , respectively the «non-corrected» and «corrected» estimation, is certainly considerable because in a power law distribution the average µ is very close to the minimum value of the distribution Another crucial aspect, until now completely neglected, concerns the effect of the bias θ 2 on the estimation of the uncertainty given by eq.(2.4) and the «corrected» form which reads The positive bias θ2 has two competitive effects.It leads to an increase in the estimated uncertainty; in fact, comparing eqs.(3.4) and (2.4) we obtain On the other hand, the bias θ 2 leads to an increase of the dispersion of the RV b ) t around its expected value.By neglecting the bias θ1, the true variance of b ) t is where var ( n t ) is the variance of the RV n t .In the same way, we have .
The validity of eqs.(3.6) and (3.7) deserves further explanations.In particular, these equations assume that E (b ) t ) = b ) and E(b t ) = b, respectively.If we take the expected value of Taylor's expansion around the true value µ of eqs.(2.3) and (3.1), we see that these assumptions hold only for small deviations of n t , i.e. for large datasets.Numerical investigations have shown that the biases are negligible for datasets with 50 or more earthquakes.
From eqs. (3.5) and (3.8), we can conclude that the true dispersion of the RV b b v ^h increases more than the increase in the estimation of the uncertainty b v ) t t .In other words, eq. ( 2.4) provides an underestimation of the true dispersion.

Binned formulas
After the correction suggested by Utsu (1966), Bender (1983), Tinti and Mulargia (1987) provided formulas to estimate the bvalue, by properly taking into account the grouping of the magnitudes.Remarkably, besides very few exceptions (e.g., Frohlich and Davis, 1983), these formulas were almost ignored in subsequent applications.We argue that the reasons are mainly of a technical nature.Bender's (1983) formula, for example, can be solved only numerically.Moreover, in her analysis she gave more emphasis to the bias θ1 introduced by the use of the continuous approximation (eq.( 2.3)), concluding that the latter provides almost unbiased estimations of the b-value if the magnitude interval for the grouping is ∆M = 0.1.
A definite improvement to the estimation of the b-value was provided by Tinti and Mulargia (1987).Their formula reads where and the associated asymptotic error is where N is the number of earthquakes.In this case, we think the very scarce use of these formulas was probably due to some kind of crypticism of the paper.

Numerical check
In order to check the reliability of the formulas described above, we simulate 1000 seismic catalogs, for different catalog sizes.The magnitudes M are obtained by binning, with ∆M = 0.1 (as for the instrumental magnitudes), a continuous RV distributed with a pdf given by eq.(2.2); in other words, Mi is the magnitude attached to all the synthetic seismic events with real continuous magnitude in the range M i -0.05 ≤ M < Mi + 0.05.
In fig.1a,b we report the medians of b ) t , b t and bTM t calculated in 1000 synthetic catalogs as a function of the number of data, for the case b = 1 and b = 2.To each median is attached the 95% confidence interval, given by the interval between the 2.5 and 97.5 percentile.From fig. 1a,b, we can see that the estimation bTM t (Tinti and Mulargia, 1987) is bias free, also for a small dataset.As regards the continuous formulas, with and without correction (respectively eqs.(3.1) and (2.3)), we can see that the bias θ2 reported in fig.1a,b is comparable to the theoretical expectation given by eq.(3.3).The corrected estimation b t is undoubtedly much closer to the real b-value.The slight underestimation of b t (much less than 1% of the real b-value) is due to the bias θ1 previously discussed (Bender, 1983).Therefore, at least for ∆M = 0.1, θ 1 can be neglected (e.g., Bender, 1983), but θ2 is certainly relevant.
In order to evaluate the reliability of the estimations of the uncertainty, it is necessary to compare each estimation with the true dis- (3.12) The null hypothesis is that the two variances are the same.
The results are reported in fig.2a,b.The most interesting result is relative to the systematic rejection of the F test for the non-corrected case, given by eq.(2.4).This means that the use of the non-corrected estimation b ) t (eq.( 2.3)) leads to an underestimation of its real dispersion.This fact is very important because it might suggest variations in the b-value whereas they do not exist.Note that the amplitude of this bias is correctly estimated by eq.(3.8).On the contrary, all the other estimations of the uncertainty do not have significant biases.It is important to note that the error estimated through eq.(2.5) (Shi and Bolt, 1982) is unbiased if the b-value is estimated through the corrected (eq.(2.3)) and non-corrected (eq.(3.1)) formulas.

Effects of the measurement errors
In order to take into account the measurement errors, we consider the «real» magnitude M + as a sum of two independent RVs where M is the earthquake magnitude devoid of measurement errors distributed with a pdf given by eq.(2.1), and ε simulates the measurement errors distributed as a Gaussian noise.The choice of the normal distribution has been discussed in detail by Tinti and Mulargia (1985), and Rhoades (1996).
In the case the two RVs, M and ε, are independent, the statistical cumulative distribution of M + can be obtained as follows (Ventsel, 1983;Tinti and Mulargia, 1985;Rhoades, 1996) where G( M + ) is the cumulative distribution of M + , f 1 (M) is given by eq.(2.1), and  where + n / is the sample average of M + .The first three addenda represent the estimation of the bvalue without measurement errors (cf.eq. ( 2.3)), while the term inside the braces is the additional part which takes into account the measurement errors.

Numerical check
The contribution of the term inside the bracket in eq.(4.5) to the b-value estimation for different σε can be studied by adding measurement errors to the synthetic catalogs generated as described before.
In fig.3a,b we report the results for different σε.Compared to the case without measurement errors (fig.1a,b), the only observable difference is a global negative bias present in all the estimations for b = 2 and σε = 0.5.This can be due to the improper binning of ∆M = 0.1 in the presence of larger measurement errors (σε = 0.5).
As regards the estimation of the uncertainties b v t t and b v ) t t , we operate as in the previous case, by plotting the results of the F test (fig.4a,b).The results are almost the same as those obtained for the case without measurement errors (fig.2a,b).Note that also in this case the use of eq.(2.4) leads to a significant underestimation of the uncertainty.

Final remarks
The purpose of this paper was to review the main concepts of the estimation of the b-value and its uncertainty, and to provide new insights on the reliability of different formulas reported in scientific literature.The detection and quantification of possible significant biases in such estimations is a crucial issue in many scientific applications, such as hazard studies and any attempt to test the constancy or the universality of the b-value.Here we studied analytically and numerically the possible biases introduced by the use of the binned magnitudes and measurement errors on different formulas.
We found that, in general, the influence of the measurement errors appears negligible com-pared to the effects of the binned magnitudes if improperly dealt with.In particular, we found that the most commonly used formulas (Aki, 1965;see eqs. (2.3) and (2.4)), which assume the magnitude as a continuous random variable, produce a strong bias in the estimation of the b-value, and a significant underestimation of its uncertainty.This means that any spatial or temporal variation of the b-value obtained by eqs.(2.3) and (2.4) has to be regarded with a strong skepticism.On the contrary, we show how the same continuous formula with a small correction to take into account the binned magnitudes (eqs.(3.1) and (3.4)) drastically reduces the biases of the b-value and of its uncertainty.We also verified that other very simple formulas provided in the past (e.g., Tinti and Mulargia, 1987) do not have significant biases and they work very well also with a small dataset contaminated by realistic measurement errors.

Fig
Fig. 2a,b.F test values (see eq. (3.12)) for the case of fig.1a,b (see text).Open and solid squares are, respectively, the F test values relative to the uncertainties given by eqs.(2.4) and (2.5).Open and solid circles are the F test values relative to the same uncertainties calculated for b t instead of b ) t (i.e.eq.(3.4)).Asterisks are the F test values relative to the uncertainty calculated through eq.(3.11).The dotted line represents the critical value to reject the null hypothesis at a significance level of 0.05.
Fig. 3a,b.As for fig.1a,b, but relative to the case with added measurement error to the magnitudes.The standard deviations of the measurement errors are σε = 0.1, 0.3, 0.5 (from left to right).