INTRODUCTION

The most common statistical tool for seismic hazard analysis is the Gutenberg - Richter (G-R) relation (^{Gutenberg and Richter, 1942}, ^{Ishimoto and Ida, 1939}, ^{Richter, 1958}), that describes the distribution of earthquake magnitudes as

where N is the number of earthquakes with magnitude greater than or equal to M. Parameter α₁ is the logarithm of the total number of earthquakes with M ≥ M ₁ , and parameter b, the slope of the relation commonly referred to as the b-value, is a measure of the relative quantities of small to large earthquakes and is usually ~1. M₁ (often denoted M _c ) is the minimum magnitude for which coverage is complete, so that for smaller magnitudes log₁₀N and M are not linearly related.

While α₁ depends on the sample time and overall seismicity rate, b is related to the local geology and to the level of ambient stress (^{Scholz, 1968}; ^{Ghosh et al., 2008}) and is related to the fractality of fractures (^{Aki, 1981}; ^{Öncel et al., 2001}), so that it varies in time and space (^{Enescu and Ito, 2002}). Several authors have found precursory changes in the b-value before the occurrence of major earthquakes (e.g. ^{Shaw et al., 1992}; ^{Wyss and Wiemer, 2000}; ^{Enescu and Ito, 2001}; ^{Márquez-Ramírez, 2012}). Hence, b is an extremely important parameter in seismic hazard studies, both for estimating seismicity occurrence rates and as a precursor to large earthquakes.

There is no explicit upper limit to the magnitudes used in (1) and, since there must be a physical limit to how large an earthquake can be, many authors have proposed ways of truncating or modifying the G-R relation to account for a maximum possible magnitude M_max (^{Utsu, 1999}). However, for most studies, the G-R relation ceases to be linear for magnitudes way below a maximum possible magnitude; obviously, magnitudes below that corresponding to log₁₀(M)<0 are either under- or over-sampled, but under- or over-sampling are common for magnitudes corresponding to log₁₀(M)<0.5 or 0.6. Let us denote by M₂ the magnitude above which over- or under-sampling occur, according to the completeness of the seismic catalog; then the G-R histogram will behave linearly only within the [M₁,M₂] range.

Often the b-value is determined by fitting a straight line to the linear part of the G-R histogram that, since magnitudes are usually rounded to one decimal place, commonly has classes ΔM = 0.1 wide.

Another common way to estimate the b-value is through the Aki-Utsu relation, which is based on the fact that the G-R relation (1) is a reverse cumulative distribution according to which magnitudes are distributed exponentially as

where β = bln(10) . (c.f. ^{Lomnitz, 1974}), and β is related to the mean of the distribution, μ, as β=1(μ-M₁) (c.f. ^{Parzen, 1960}), so that

^{Aki (1965)} showed that the maximum-likelihood estimate of b, b_m , is given by

where M- is the sample mean, and ^{Utsu (1965)} pointed out that, since magnitudes are rounded to ΔM, the actual minimum magnitude is M1U=M1-∆M/2, so that

Formula (5), which we will refer to as the Aki-Utsu estimate, has been widely used as a simple and straightforward way of estimating b directly from the magnitude sample mean, with no explicit need for a G-R histogram. However, in too many cases, people do not realize the need for correctly determining M1U, and of having a representative sample from the linear part of the histogram that is large enough so that M-≅μ, which is an absolute requirement for b_m ≅ b. In too many cases M- is estimated from too small samples, either because data are scarce or because no attention is paid to the matter of representativity, and small samples result in estimates with too much scatter to be useful (^{Kramer, 2014}; ^{Nava et al., 2017a}).

Whichever the method used to estimate the b-value, the data has to fulfill two requirements: first, the number of data in the[M₁,M₂] range and the range itself need to be large enough so that a straight line can be adequately fitted or so that the observed mean magnitude is representative of the distribution mean; second, the data should be homogeneous.

The first requirement is not particularly difficult to meet when considering a large area or a long-time history, but when trying to have a good definition in time and/or space, which requires short time and/or space windows, then having a representative sample may be difficult.

The second requirement can have three aspects. Homogeneity in time: network coverage usually changes in time, but there are sophisticated techniques to deal, at least partially, with this aspect (e.g. ^{Kijko, 2004}; ^{Kijko and Smit, 2014}); also, b can change in time, in which case the measured value will be an average over time. Homogeneity in space: b does change from place to place so that a measure using data from a large region will yield a space average of the local values. Homogeneity in magnitude: it is not uncommon to report moment magnitudes for large earthquakes and use some other scale, such as coda or duration, for small and very small earthquakes; unless the small earthquake scale is correctly calibrated so that it measures like the large magnitude scale, then b-value determinations will be erroneous.

We next describe the problems encountered with this third requirement, while trying to determine whether changes in the b-value were observable before and after large earthquakes in the Mexican subduction zone.

The Mexican Subduction Zone

The tectonic activity in the south and southeast of México is governed by the subduction of the Rivera and Cocos plates under the North American plate (Figure 1). It is along this subduction zone where the largest earthquakes in Mexico have occurred. The subduction dip angle of the Cocos plates changes from sub-horizontal below central México to major dip [25⁰-30⁰] near Chiapas to the southeast (^{Pérez-Campos et al., 2008}). The convergence velocity also changes, the Cocos plate has a velocity of 4 to 5 cm/yr in its western part, and the eastern part has velocity around 6 to 7 cm/yr relative to the North American plate (^{Dañobeitia et al., 2016}; ^{Nuñez-Cornú et al., 2016}; ^{Gutiérrez et al., 2015}; ^{Kostoglodov and Bandy, 1995}; ^{Kostoglodov and Pacheco, 1999}).

Figure 1 Tectonic plates interaction in México. Numbers indicate the velocity in cm/yr, and arrows show the local direction of the subducting plates relative to the North American plate. Red for the Cocos plate, Green for the Rivera plate, and Yellow for the Pacific plate.  

International catalogs do not list magnitudes small enough for reliable b-value determinations in this region, so that we decided to use data from the Servicio Sismológico Nacional (SSN, Mexico’s National Seismological Service) seismic catalog, from 1988 (^{Zúñiga et al., 2000}) to 2019.

Event Selection

We considered all earthquakes with M ≥ 7.0 along the subduction zone, and selected a region around each mainshock, according to the spatial distribution of its aftershocks, to sample regions subject to the stresses that would cause the mainshock. Some of the regions were adjusted to avoid including events associated with another mainshock. Only mainshocks with enough data for pre- and post-event time windows were considered for b determination; Figure 2 shows the chosen areas, and Table 1 lists their magnitude, occurrence time, location, and the total number of events in window.

Figure 2 The squares indicate regions used for b value determination for large earthquakes along the Mexican subduction zone; from West to East: Guerrero 2014, Oaxaca 2012, Oaxaca 2018, Oaxaca-Chiapas 2017, Chiapas-Guatemala 2012. Dots represent epicenters: red M < 4.0, green 4.0 ≤ M < 5.0, blue 5.0 ≤ M < 6.0, grey 6.0 ≤ M < 7.0; Yellow diamond with 7.0 ≤ M < 8.0, and Yellow stars with M ≥ 8.0.  

Table 1 Earthquakes with enough data for b determination. 

M	Time	Lon	Lat	Depth	State Year	Total
7.50	2012.2179	-98.4570	16.2640	18.00	Oaxaca 2012	8927
7.30	2012.8516	-92.3160	14.0270	17.10	Chiapas/Guatemala 2012	6637
7.20	2014.2948	-101.4600	17.0110	18.00	Guerrero 2014	3857
8.20	2017.6855	-94.1030	14.7610	45.90	Oaxaca/Chiapas 2017	30143
7.20	2018.1287	-98.0140	16.2180	16.00	Oaxaca 2018	16833

The cumulative number of earthquakes curve for each region was plotted to determine times and threshold magnitudes for homogeneity; then, the cumulative curves were used to determine time windows. One of the time windows was chosen before the mainshock when stresses are expected to be high, and since the catalog does not extend backwards in time long enough to sample b-values before the stress build-up, in order to have a low-stress reference value, a second window was chosen after the mainshock liberated the stored stress and after the significant part of the aftershock activity, when seismicity was back to background level, so that we could sample low-stress b-values without aftershock noise.

As an example, Figure 3 shows the selection of pre and post-events for the 2012.85, M 7.3 earthquake located at the Mexico (State of Chiapas) ( Guatemala border. From the cumulative curve for the whole region (Figure 3-Top), a clear change in the slope, possibly due to changes in the seismic network, is apparent around 2011, so that we used events from this point on to the end of the catalog. The pre- and post-mainshock time windows, which will be referred to as W1 and W2 from now on, are shown in Figure 3 (Bottom).

Figure 3 Cumulative curve for the 2012.85 earthquake located in the border of Chiapas/Guatemala. Top: Complete catalog. Bottom: Time windows, pre-earthquake W1 with 1109 events, and post-earthquake W2 with 2463. 

RESULTS

Figure 3 shows the time window for the Chiapas-Guatemala 2012 earthquake; Figure 4 shows the windows selected for the rest of the earthquakes listed in Table 1. Because the aftershock sequence of the M _W 8.2 Oaxaca-Chiapas earthquake is not finished yet, it was not possible to have a post-quake window for this event.

Figure 4 Pre and post-mainshock time windows for the earthquakes listed in Table 1. 

Figures 5 to 9 show the fit of b_m to the G-R distribution (left), and the likelihood distribution and b_x choice (right) for events listed in Table 1; the b_m and b_x values are listed in Table 2.

Figure 5 Oaxaca 2012. Left: G-R fit for b_m ; (blue) circles are the G-R histogram, and (red) diamonds show the corresponding non-cumulative distribution; N is the number of data in the linear range from M₁ to M₂ ; M.L. is the straight line corresponding to b_m the maximum likelihood fit; L.S. is the dashed line for the least-squares fit (shown for comparison only). Right: b likelihood distribution; Δb is the class width, N_r is the number of Monte Carlo realizations, s is the pseudo-random number generator seed; the thin dashed red line indicates b_m , the black dash-dot line indicates b_x . Short vertical lines with crosses and asterisks indicate 75% and 90% confidence ranges, respectively. Top: Pre-event W1. Bottom: Post-event W2.  

Figure 6 Chiapas-Guatemala 2012. Same conventions as in Figure 5. 

Figure 7 Guerrero 2014. Same conventions as in Figure 5. 

Figure 8 Oaxaca-Chiapas 2017. Same conventions as in Figure 5. 

Figure 9 Oaxaca 2018. Same conventions as in Figure 5. 

From these figures, it is clear that the G-R distributions do not have a single slope; all of them show a large slope for small magnitudes and a smaller one for larger magnitudes, meeting around magnitude 4.5. This feature is explained when topic Mitos y Realidades, in Sismos y Volcanes CDMX, mobile application software (2019), is consulted; the app. states that the SSN employs M_w for events larger than 4.5 and coda magnitude for events smaller than 4.5 (which scale is used for 4.5 magnitude events remains a mystery).

Because of this change in slope, it was not possible to obtain linear ranges over a wide enough magnitude interval; due to space and time limitations of the windows, only the small magnitude range was (barely) adequate for b_m estimation. As evidenced by the non-cumulative histograms shown in the G-R plots, there were not enough data in the larger magnitudes range to obtain adequate estimates. Thus, we had to base our estimates on short ranges over only the smaller magnitudes.

A characteristic of the most likely source b-value method is that when the linear range is wide, some two or more magnitude units, say, and the number of data within the range is larger than about 2,000; the source b distributions are narrow and b_m and b_x are equal or differ by little. For the determinations presented here, with narrow ranges and few data, the source b distributions are wide.

We will now proceed to describe our results and defer their interpretation and assessment to the next section. The results are summarized in Table 2.

DISCUSSION AND CONCLUSIONS

The change in slope around M4.5 made it necessary to estimate b values using the smaller magnitudes only, and the measured b-values are larger than the semi-theoretical 1.5 maximum value (^{Olsson, 1999}). These large values and the change in slope indicate that the M _c scale used by the SSN is not consistent with the M _W scale, so that the above mentioned maximum value does not apply to SSN small magnitudes. There were not enough data to obtain Aki-Utsu b-value estimates from the larger magnitudes, but least-squares fits (where possible for large magnitudes) yield values smaller than 1.5.

Hence, the b-values obtained here are useful only for comparisons among themselves, and cannot be used for comparisons with values obtained from data sets with true moment magnitudes.

The b-values for the Oaxaca 2012 earthquake are definitely smaller for the high stress regime before the mainshock than for the low-stress regime after it. Results for other earthquakes consistently show b-values to be smaller before than after the main events, but the spreads in source b distributions make it impossible to discard the corresponding null hypotheses with any significant degree of confidence. We consider, however, that these results do strongly suggest smaller b values before the mainshocks than in low-stress regimes.

The question whether b-values are a useful precursor tool for the Mexican subduction zone remains an open question and will remain so until the SSN scales for small and large magnitudes agree (hopefully both scaling as M _W ), so that reliable b-value determinations based on an appropriately wide magnitude range are possible. Meanwhile, let these observations be a caveat for researchers planning to work with b-values from the Mexican subduction zone.