Temporal and spatial changes in the b-value prior to the 2021 Luxian M S 6.0 earthquake in Sichuan, China

Abstract Using the seismic data collected since 2010 by the China Earthquake Networks Center and based on the maximum likelihood method, we analyzed the temporal and spatial anomalies of the b-value prior to the M S 6.0 Luxian earthquake on September 16, 2021. The results are as follows. (1) The Luxian earthquake broke the historical record of no M ≥ 6 earthquake occurring in the Huayingshan fault zone, revealing that the Huayingshan fault zone still has the seismogenic ability of medium and strong earthquakes. (2) Before the M S 6.0 Luxian earthquake, the b-values around the focal areas exhibited increasing–peaking–decreasing anomalous characteristics, indicating that this is an effective method of predicting the occurrence of moderate to strong earthquakes in the region. (3) In the past 10 years, the b-values of the middle segment of the Huayingshan fault zone have been relatively low, but the b-value of this section may be in an abnormal stage of slow increase, indicating that this segment may be in a stage of stress accumulation and concentration and is preparing for moderate to strong earthquakes. The middle segment of the Huayingshan fault zone may become the site of moderate to strong earthquakes in the future. We should continue to pay attention to it.


Introduction
The Gutenberg-Richter (GR) law expresses the earthquake magnitude-frequency relationship (Gutenberg and Richter 1944): LgN ¼ aÀbM, where the parameter b (known as the b-value) represents the factor proportional to the number of earthquakes of varying magnitudes in a region. Early laboratory rock-fracturing experiments (Scholz 1968), fluid extraction-related seismic activity experiments (Wyss 1973), and underground mine rock fracture experiments (Urbancic et al. 1992) revealed that the magnitude of the tectonic stress in a specific rock mass is inversely proportional to the magnitude of the b-value. Many studies have also shown that region includes the southern (blue dashed box, zone A in Figure 1(a)) and middle (blue dashed box, zone B in Figure 1(a)) segments of the Huayingshan fault zone. The earthquake catalog for the period from 1 January 2010 to 27 September 2021 used in this study was obtained from the China Earthquake Networks Center 1 . A total of 53,367 earthquakes with M ! 0.1 were included, including one magnitude 6.0 earthquake, three earthquakes with magnitudes of 5.0-5.9, 21 earthquakes with magnitudes of 4.0-4.9, 400 earthquakes with magnitudes of 3.0-3.9, and 5996 earthquakes with magnitudes of 2.0-2.9. There were 29,822 earthquakes with magnitudes of 1.0 to 1.9, and 17,124 earthquakes with magnitudes of less than 1.0 (The nation-wide and regional relationships between M S and M L in China are empirically obtained by  Wang and Yu (2009) derived by 6577 shallow earthquake magnitude data from the 'Annual Bulletin of Chinese Earthquakes' for 1990-2007, they obtained the empirical relationship between M S and M L is close to M S ＝ M L. Therefore, the conversion between M S and M L was not conducted in this paper, We use M to express the magnitude scale of an earthquake). The distributions of the earthquakes' epicenters and faults in the study area are shown in Figure 1(a).
In order to avoid the influence of data selection and obtain a reasonable and reliable earthquake catalog for further study, the analysis of the completeness of the earthquake catalog is very important. At present, the commonly used methods of estimating the minimum magnitude of completeness of a region include the entire-magnitude-range (EMR) method (Ogata and Katsura 1993), the maximum curvature (MAXC) method , the goodness of fit (GFT) method , and the M c estimated from the b-value stability (MBS) (Cao and Gao 2002). Among these methods, the MAXC and GFT methods are simple to use, but the resultant M c is relatively low and is often underestimated. The M c obtained using the MBS method is relatively high, and this method is greatly affected by the sample size. Thus, MBS method is only suitable for regional catalogues and has a large uncertainty. In contrast, the EMR method is able to yield an M c value that is stable, reliable, and closer to the true value. In this study, the Entire-Magnitude-Range (EMR) method (Woessner and Wiemer 2005) was used to obtain the magnitude of the completeness in the study area (M c ¼ 1.4; Figure 1(c)). The seismic data for the M ! 1.4 earthquakes were used for the calculation and analysis of the b-value in the study area.

Method
The methods commonly used to calculate the b-value include the maximum likelihood method (MLE) (Aki 1965;Utsu 1965;Shi and Bolt 1982) and the least squares method (LSM) (Okal and Kirby 1995;Main 2000;Z€ oller et al. 2002). The MLE regards the contribution of each seismic event as the same and averages the magnitudes of all of the earthquakes with the same weight. Compared with the LSM, the MLE weights small and medium earthquakes more. A large number of studies (Goldstein et al. 2004;Newman 2005;Bengoubou-Val erius and Gibert 2013;Nava et al. 2017) have shown that the standard deviation of the LSM is more than twice the standard deviation of the MLE. The MLE can obtain an accurate estimate when the sample size is small. Therefore, the MLE was used in this study to calculate bvalue: In Equation (3-1), N is the total number of earthquake samples, M 0 is the initial magnitude, lge ¼ 0.4343, and M is the average magnitude of a group of earthquakes, namely, M ¼ 1 N P N i¼1 M i : According to the results of Shi and Bolt (1982), the standard deviation of the b MLE , db, is: (3-2) In Equation (3-2), N is the total number of earthquake samples, M i is the magnitude of each earthquake, and M is the average magnitude of a group of earthquakes, namely, M ¼ 1

Temporal distribution of b-values
To reflect the change in the b-value in the study area over the past 10 years, the seismic data for the M ! M c earthquakes in the study area since 2010 were used. We used the MLE method to calculate the temporal changes in the b-values. The sample window size was 500 earthquakes. The minimum number of earthquakes was 50. The number of bootstraps was 20. The calculation results are shown in Figure 2(a). As is shown in Figure 2 the differences in the b-values before and after 2015 were significant. Based on the analysis of the frequency-time distribution map of the M ! Mc earthquakes ( Figure  2(b)) and the seismic monitoring capabilities of the study area, we believe that this may be related to the greatly improved monitoring capacity, monitoring level, and quality of the observation data of the Sichuan Earthquake Monitoring Network after 2015. As a result, more low-magnitude seismic data became suddenly available for calculating the b-value, thereby increasing the b-value. From 2018 to the present, the b-values in the study area fluctuated greatly, with many sudden increases and decreases. Based on the solid red lines in Figure 2(a) and the regional earthquake catalog, we found that since 2018, there have been 11 moderate to strong earthquakes of magnitude 4.5 or above in the study area. Before many moderate to strong earthquakes, the b-value exhibits anomalous characteristics, that is, it increases suddenly, peaks, and then gradually decreases to a low value. Earthquakes occurred in the time periods when the b-value changed from a peak to a low value, such as the earthquakes 1-6 shown in Figure 2(a), including the magnitude 6.0 Luxian earthquake in Sichuan on 16 September 2021. Such anomalous characteristics appeared before theses earthquakes occurred. The change in the b-value reflects the change in the local stress state, and its anomalous characteristics can also be used as one of the effective methods of predicting the occurrence of strong earthquakes in this region. The above analysis of the characteristics of the temporal changes in the b-values reflects the stress changes and b-value changes in the entire study area. It can be seen from Figure 1(a) that since 2010, three earthquakes with magnitudes of 5.0 to 5.9 (solid purple circles in Figure 1(a)) occurred in Weiyuan County outside of the Huayingshan fault zone, and in the Qingjiangbai District, Chengdu, Sichuan, in the northwestern corner of the region. Within range of 30 km around the Huayingshan fault zone, there were no earthquakes with magnitudes of 5.0 or higher from 2010 until the Luxian earthquake. The statistics for a large region may conceal local regional anomalies. In order to reveal the anomalous characteristics of the b-value prior to the Luxian earthquake and the potential seismic hazards for the other areas of the Huayingshan fault zone, we analyzed the temporal changes in the b-values in the Huayingshan fault zone in more detail, and used the MLE method to calculate the temporal changes in the b-values in two small areas within the Huayingshan fault zone: the southern segment (zone A) and the middle segment (zone B). The parameters for zones A and B were selected as follows. The sample window sizes for zones A and B were 100 and 200 earthquakes, respectively. The minimum numbers for zones A and B were 50 earthquakes, respectively, and the number of bootstraps for zones A and B was 20. The calculation results for Zone A are shown in Figure 3  Peak segment and low-value segment 1.0 ÃÃ 1.932 Note: the KS test was used to test the difference between the high b-value segment and the low b-value value segment before several earthquakes. ÃÃ Indicates at significance levels of p < 0.05 and p < 0.01, the differences in the b-values between these two segments are significant. The D value is the absolute difference between the cumulative probability of the sample size distribution of the two groups. The larger the D value is, the greater the gap in the distribution between the two groups of data is. The Z value represents the test statistic; and # indicates that the difference is not significant.
It can be seen from Figure 3(a) that the b A -value of the southern segment of the Huayingshan fault zone varied between 0.584 and 1.333, with an average of 0.853. Before several earthquakes with magnitudes of 4.5 or above occurred in zone A (earthquakes 1 to 5 in Figure 3(a)), the b-value suddenly increased and then gradually decreased to a low value after peaking. We conducted the KS test on the curve segments with high and low b-values prior to these earthquakes in zone A. The test results are shown in Table 1. The D values in Table 1 are the absolute differences between the cumulative probability of the distribution of the sample sizes of the two groups. The larger the D value is, the larger the gap between the distribution of two groups of data is, and the more significant the difference is. Z is the test statistic. The analysis results show that at significance levels of p < 0.05 and p < 0.01, except for the M4.8 earthquake that occurred on 28 December 2016 (earthquake #2), the differences in the pre-earthquake b-values were significant for the other earthquakes. The magnitude 4.6 earthquake on 23 July 2021 occurred two months before the Luxian earthquake, and its epicenter was about 18 km northeast of the epicenter of the Luxian earthquake. Figure 3(b) clearly shows that 14 months before this earthquake, the bvalue of the southern segment of the Huayingshan fault zone began to increase slowly and reached a peak (about 1.13) about 8 months before the earthquake. Then, it slowly decreased and the magnitude 4.6 earthquake occurred. After the magnitude 4.6 earthquake occurred, the b-value in zone A increased rapidly and decreased rapidly within a short period of time (about 20 days). Then, the magnitude 6.0 Luxian earthquake occurred. These phenomena reflect the stress states in different periods of time during the preparation for the moderate to strong earthquakes in the southern segment of the Huayingshan fault zone, and they also reflect the stress accumulation during the earthquake preparation periods and the stress release process after the earthquakes occurred.  Table 1) that at significance levels of p < 0.05 and p < 0.01, the differences in the b-values between these two segments in zone B are significant. Based on the earthquake catalog for zone B, we found that zone B is currently dominated by small earthquakes, and there have not been any M ! 4 earthquakes. However, the fact that the b-value of zone B has been in a slowly increasing stage should not be ignored since this indicates that zone B may be in the increasing stage of the increasing-peaking-decreasing process before the occurrence of a moderate to strong earthquake. In preparation for future moderate to strong earthquakes, continuous attention should be paid to the temporal evolution of the future b-values in this area.

Spatial distribution of the b-value
The b-values reflect the regional crustal stress rupture strength and state. The spatial scanning results of the b-value can be used to identify areas with low b-values, and the low b-values can be used as high effective shear stress indicators to reveal and infer the relative current accumulated stress levels of the different segments of the active fault zone, to identify possible asperities, and to further determine the possibility of strong earthquakes occurring in the active fault zone. In order to explore the spatial differences in the b-values during different periods and in different regions of the focal areas and neighboring areas before the M6.0 Luxian earthquake occurred, we divided the b-value spatial scanning into four time periods according to the temporal distribution characteristics of the b-values in the above-mentioned study area. We obtained the spatial distributions of the b-values in the different time periods and the different regions, which reflect the different stress accumulation levels and seismic activity levels in these time periods and regions.
The four periods in which the b-value spatial scanning was conducted are as follows. The first period was from 1 January 2010 to 15 September 2021 (prior to the magnitude 6.0 Luxian earthquake). The second period was from 11 November 2012 to 31 December 2016. The third period was from 1 January 2017 to 30 April 2020. The fourth period was from 1 May 2020 to 15 September 2021. We established an appropriate scheme for the spatial scanning of the b-value based on the collected seismic data. For a fixed grid cell size, as the scanning radius increases, the blank area on the b-value spatial distribution map decreases. Therefore, the scanning radius r should be appropriately increased for areas with a smaller number of seismic records. In addition, for the same scanning radius r, although the overall patterns of the bvalue spatial distribution maps calculated using different mesh sizes are roughly the same, the b-value spatial distribution map calculated using a smaller mesh size (e.g. 0.02 Â0.02 ) is more refined than those calculated using a larger mesh size (such as 0.1 Â0.1 or 0.05 Â0.05 ). After many tests and comparisons, we chose a mesh size of 0.02 Â 0.02 to set up the grid. With each grid cell's node as the center of a circle, the earthquakes with M ! M c were selected within the circular statistical unit with a radius of 50 km (to ensure that there were enough samples, grid cells with less than 50 earthquakes were not included in the calculation and are shown as blank), and the MLE was used to perform the spatial scanning of the b-values. We obtained a spatial distribution map of the b-values before the M6.0 Luxian earthquake ( Figure  4), a spatial distribution map of the standard deviations of the b-values ( Figure 5), a spatial distribution map of the percentage fit of the G-R power law equation (Figure  6), as well as the mean b values and the mean standard deviation values db in the different time periods and different regions (Table 2). Figure 4 shows that the spatial distributions of the b-values in different time periods and regions in the study area were extremely uneven and exhibit considerable differences, reflecting the differences in the seismic activity levels in the different regions and time periods, and thus reflecting the spatial differences in the seismic hazard. Figure 5 shows the spatial distributions of the standard deviations of the bvalues during the different time periods. The low standard deviation values clearly indicate that the estimated values of the spatial b-values are statistically significant. Figure 6 shows the spatial distribution of the percentage fit of G-R power law equation during the different time periods, which was used to estimate the goodness of fit. We first computed the absolute difference between the observed and synthetic distribution of the number of events in each magnitude bin. Then, we divided it by the total number of observed events to normalize the distribution (Singh et al. 2012). Figure 6 shows that the overall spatial estimate of the percentage fit of the G-R power law is >80% in most of the regions, indicating that the estimation of the b-value is statically significant.
As can be seen from Table 2, the b-values of the middle and southern segments of the Huayingshan fault zone were lower in the first and second periods than in the entire study area. In particular, the average b value (0.742 ± 0.041 and 0.726 ± 0.050, respectively) of the middle segment was much lower than the average value (0.95 ± 0.05 and 0.971 ± 0.064, respectively) for the study area. We used the KS test to  Table 2. At significance levels of p < 0.05 and p < 0.01, the b-values of the middle and southern segments of the Huayingshan fault zone significantly differed from that of the study area, and the difference between the middle segment and the study area was larger than that between the southern segment and the study area. During the first and second periods, in the Huayingshan fault zone, the b-value of the middle segment has been lower than that of the southern segment, indicating that the Huayingshan fault zone has a relatively high level of stress accumulation. In the third period, the b-value of the southern part of the Huayingshan fault zone was still lower than the average level of the entire region. The b-value (1.050 ± 0.069) of the middle segment of the Huayingshan fault zone increased. This is consistent with the temporal distribution characteristics of the gradual increase in the b-value in zone B after 2017 and indicates that the area is currently dominated by small and medium earthquakes. We should continue to pay attention to the temporal and spatial changes in the b-value in this segment, which can provide information for predicting the occurrence of moderate to strong earthquakes in the future. The fourth period was the temporal anomaly period of the bvalue before the M6.0 Luxian earthquake (16 months) occurred. As can be seen from the spatial distribution map of the b-value for this period, the b-value around the focal areas of the Luxian earthquake was relatively low, much lower than the average b-value of the study area and lower than the b-values in the previous three periods. The low b-value reflected the fact that the southern segment of the Huayingshan fault zone was in a relatively high stress state before the M6.0 Luxian earthquake, confirming that the spatial anomalies of low b-value can be used as one of the effective methods of determining the hazard of strong earthquakes in the active fault zone. The KS test was used to analyze the significance of the differences between the two sample sizes of the calculated b-values for the southern and middle segments of the Huayingshan fault zone and the b-value of the study area, ÃÃ Indicates at significance levels of p < 0.05 and p < 0.01, the differences in the b-values between these two sample are significant. The D value is the absolute difference between the cumulative probability of the distribution of the two sample sizes of the two groups. The larger the D value is, the greater the gap in the distribution between the two groups of data is.

Conclusions
In this study, we analyzed the temporal and spatial distribution characteristics of the b-values in the focal areas and adjacent areas before the M S 6.0 Luxian earthquake occurred on 16 September 2021. The results of this study follow.
1. The M S 6.0 Luxian earthquake was the tectonic response of the Bayan Har Block in western Sichuan pushing into the Sichuan Basin. It was the first historical earthquake with a magnitude of 6 or higher to occur in the Huayingshan fault zone. The Huayingshan fault zone does not have late Quaternary activity, but it is a seismogenic structure. This earthquake revealed that the Huayingshan fault zone still has the seismogenic ability to produce moderate to strong earthquakes, and there is still a hazard of earthquakes of magnitude 6 or higher occurring. 2. Before the M S 6.0 Luxian earthquake, the b-values around the focal areas (the southern segment of the Huayingshan fault zone) exhibited increasing-peaking-decreasing anomalous characteristics, and the anomalous characteristics of the low b-values appeared in the focal areas and their vicinities. This reflects the stress changes during the build-up to moderate to strong earthquakes in the southern segment of the Huayingshan fault zone and suggests that temporal and spatial b-value anomalies can be used as one of the effective method to predict the occurrence of moderate to strong earthquakes in this region. 3. In the past 10 years, the b-values of the middle segment of the Huayingshan fault zone have been relatively low, but the b-value of this section may be in an abnormal stage of slow increase, indicating that this segment may be in a stage of stress accumulation and concentration and is building up to moderate to strong earthquakes. The middle segment of the Huayingshan fault zone may become the site of moderate to strong earthquakes in the future. We should continue to pay attention to the temporal and spatial development of the b-values in this segment. Note 1. http://10.5.160.18/