Spatiotemporal heterogeneity of b values revealed by a data-driven approach for the 17 June 2019 M_S 6.0 Changning earthquake sequence, Sichuan, China

The spatiotemporal heterogeneity of b values has great potential for helping in understanding the seismogenic process and assessing seismic hazard. However, there is still much controversy about whether it exists or not, and an important reason is that the choice of subjective parameters has eroded the foundations of much research. To overcome this problem, we used a recently developed non-parametric method based on a data-driven concept to calculate b values. The major steps of this method include (1) performing a large number of Voronoi tessellations and Bayesian information criterion (BIC) value calculation, selecting the optimal models for the study area, and (2) using the ensemble median (Q₂) and median absolute deviation (MAD) value to represent the final b value and its uncertainty. We investigated spatiotemporal variations in b values before and after the 2019 Changning M_S=6.0 earthquake in the Sichuan Basin, China. The results reveal a spatial volume with low pre-mainshock b values near the mainshock source region, and its size corresponds roughly with the rupture area of the mainshock. The anomalously high pre-mainshock b values distributed in the NW direction of the epicenter were interpreted to be related to fluid invasion. The decreases in b values during the aftershock sequence along with the occurrences of several strong aftershocks imply that b values could be an indicator of the stress state. In addition, we found that although the distribution characteristics of b values obtained from different methods of investigation are qualitatively consistent, they differ significantly in terms of their specific values, suggesting that the best way to study the heterogeneous pattern of b values is in the joint dimension of space-time rather than separately in time and space. Overall, our study emphasizes the importance of b-value studies in assessing earthquake hazards.

1 Introduction

The Gutenberg–Richter b value describes the corresponding frequency–magnitude distribution (FMD) characteristics by reflecting the relative proportion of the frequency of large and small earthquakes within a given space-time range. It is considered to be related to the stress conditions in the Earth's crust (e.g., Wyss, 1973; Urbancic et al., 1992; Mori and Abercrombie, 1997; Toda et al., 1998), complexity of the fault trace (Stirling et al., 1996), and the extent of creep (Amelung and King, 1997) and other factors. Experimental studies in the laboratory have shown that a weak and less resistant environment under stress would produce a high b value, while materials that are more compact and more resistant under pressure do not fail, which leads to a reasonably low b value (Aktar et al., 2004). In the case where the material and structure are clarified, a decreasing b value is considered to be related to increasing stress (Scholz, 1968) or pore pressure diffusion (Hainzl and Fischer, 2002; Lei and Satoh, 2007). For the above reasons, the b value has been widely considered in seismogenic environment analysis and seismic hazard research.

Spatial and temporal heterogeneity is an important topic in b-value research, especially under the assumptions that the local b values are inversely dependent on the applied shear stress and that low b values (b<0.7) can reflect the existence of locked faults or asperities. Therefore, the spatial and temporal heterogeneity of b values is considered an important clue for forecasting the location and size of potentially large earthquakes (Wiemer and Wyss, 1997; Schorlemmer and Wiemer, 2005; Murru et al., 2007). Using the spatial heterogeneity of the b value to identify possible asperities has been performed in some cases, such as for the San Jacinto–Elsinore fault system in southern California (Wyss et al., 2000), the Parkfield segment of the San Andreas fault (Wiemer and Wyss, 1997), and the case study of the 2014 Parkfield M=6.0 earthquake (Schorlemmer and Wiemer, 2005).

A model named the Asperity Likelihood Model (ALM) based on the above assumptions has been developed and used to forecast future earthquakes (Wiemer and Schorlemmer, 2007; Gulia et al., 2010). Research on the temporal heterogeneity of b values mainly includes using b-value time variation in the early aftershock sequence and the constructed system of the foreshock traffic light system (FTLS) to evaluate the risk of larger subsequent aftershocks (Gulia and Wiemer, 2019).

However, some research results show that the apparent variability in b values is not significant in some cases (Del Pezzo et al., 2003). For example, Amorèse et al. (2010) systematically examined the variation in b values in southern California to the depth of the crust and found that the hypothesis was not statistically significant. By using a data-driven approach, Kamer and Hiemer (2015) show that the spatial b values in most locations in California are distributed within a very limited range (0.94±0.04–1.15±0.06), and the previously reported spatial b-value variation is overestimated and mainly due to the subjective choice of parameters. Besides, the spatial and temporal heterogeneity of b values is also considered to be due to the subjective arbitrariness of the calculation rules and the lack of statistical robustness (Kagan, 1999).

Based on the above viewpoints, the calculation reliability for research on the spatiotemporal heterogeneity of b values still needs to be solved, and the relationship between the spatiotemporal variation process of b values and the occurrence of strong earthquakes needs to be investigated for more earthquake cases. In this study, we will utilize data-driven b-value calculation methods that have been developed in recent years (Kamer and Hiemer, 2015; Nandan et al., 2017; Si and Jiang, 2019) for case studies of the 2019 Changning M_S=6.0 earthquake in Sichuan, China.

2 Method

In the traditional calculation of the Guttenberg–Richter magnitude–frequency b value, a fixed number of earthquakes (Hutton et al., 2010; Ogata, 2011) or a fixed minimum and maximum selection radius (Woessner and Wiemer, 2005) are generally used to select data and the maximum likelihood estimation is used to obtain b values. Because such calculations have strong subjectivity in their calculating rules, they have caused widespread controversy. The data-driven approaches to seismicity parameter calculation have been gradually developed in recent years (Sambridge et al., 2013; Kamer and Hiemer, 2015; Nandan et al., 2017; Si and Jiang, 2019) by using the Voronoi tessellation to create a large number of spatially random grids and covering the possibility of the segmentation of spatial regions, relying on the Bayesian information criterion (BIC) to select a part of the optimal models with the smallest BIC value and representing the final result of seismic activity parameters through the ensemble median value. Because the data-driven approach uses an automatic parametric calculation, it provides a possibility of solving the subjective problem of earthquake data selection.

Among those data-driven approaches, Si and Jiang (2019) developed a method using the continuous distribution function (hereafter referred to as the OK1993 model) given by Ogata and Katsura (1993), which has the advantage of simultaneously determining the minimum magnitude of completeness and obtaining b values. In this paper, we will use this approach to study the spatiotemporal heterogeneity of b values for the 2019 Changning M_S=6.0 earthquake.

The OK1993 model uses the seismic detection rate function q(M) to describe the complete detection degree of earthquake events with different magnitudes in the magnitude–frequency distribution:

$$\begin{array}{}\text{(1)}& q(M\mathrm{|}\mathit{\mu },\mathit{\sigma })={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}}}\underset{-\mathrm{\infty }}{\overset{M}{\int }}{e}^{-\frac{(x-\mathit{\mu }{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\mathrm{d}x},\end{array}$$

where M is the magnitude, the parameter μ represents the corresponding magnitude to the detection rate of 50 %, and σ indicates the corresponding magnitude range. The actual earthquake probability density function and the log-likelihood function of the OK1993 model can be expressed as

$$\begin{array}{}\text{(2)}& \begin{array}{rl}P(M\mathrm{|}\mathit{\beta },\mathit{\mu },\mathit{\sigma })& ={\displaystyle \frac{e^{-\mathit{\beta }M}q(M\mathrm{|}\mathit{\mu },\mathit{\sigma })}{\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}{e}^{-\mathit{\beta }M}q(M\mathrm{|}\mathit{\mu },\mathit{\sigma })\mathrm{d}M}}\\ & =\mathit{\beta }{e}^{-\mathit{\beta }(M-\mathit{\mu })+{\mathit{\beta }}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}/\mathrm{2}}q(M\mathrm{|}\mathit{\mu },\mathit{\sigma }),\end{array}\end{array}$$

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\ln L\left(\mathit{\theta }\right)& =n\ln \mathit{\beta }-\sum _{i=\mathrm{1}}^n\left[\mathit{\beta }{M}_i-\ln q\left(M_i\mathrm{|}\mathit{\mu },\mathit{\sigma }\right)\right]+n\mathit{\beta }\mathit{\mu }\\ & -{\displaystyle \frac{n}{\mathrm{2}}}{\mathit{\beta }}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}.\end{array}\end{array}$$

{M₁, M₂, …, M_n} in the above formula is the magnitude of a given series of observational events and the power exponent β=bln 10, The parameter [β, μ, σ] can be obtained by fitting the above formula using the maximum likelihood method. The Bayesian information criterion $\mathrm{BIC}=-\ln L\left(\mathit{\theta }\right)+k/\mathrm{2}\ln \left(n\right)$ is adopted to calculate the corresponding BIC value and select the optimal models. Since each grid node is composed of spatial coordinates [x, y] and three parameters [β, μ, σ] in the OK1993 model, the total number of freedom degrees is $k=\mathrm{5}\times \mathrm{no}.\mathrm{of}\mathrm{nodes}$ in the entire study region.

The construction of the data-driven approach can be achieved by Voronoi tessellation with limited boundaries. Voronoi tessellation refers to a unique set of continuous polygon partitioning schemes {P_i, i=1, 2, …, n} given by a set of spatial nodes S={s₁, s₂, …, s_n} in two-dimensional or three-dimensional space. The polygon $P_i=\mathit{\left\{}x\mathrm{|}\mathrm{dist}\right(x,s_i)\le \mathrm{dist}(x,s_j)$, i≠j}, where dist(a, b) denotes the Euclidean distance between two points. Voronoi tessellation also benefits from the uniqueness of its spatial division, so it is widely used in computer science, political elections, and many other studies (Rubner et al., 2000; Svec et al., 2007). The calculation steps of the data-driven approach include the following: (1) randomly throwing a certain number of nodes into the study area and performing Voronoi meshing, with the number of grid nodes gradually increasing from 2 to 40, ensuring that the Voronoi tessellation covers the possibility of various spatial region segmentation by each number of grid nodes being randomly thrown 100 times; (2) calculating OK1993 model parameters and BIC values for $(\mathrm{2}+\mathrm{3}+$ … $\mathrm{40})\times \mathrm{100}=\mathrm{81}\mathrm{900}$ Voronoi cells obtained from 3900 tessellations (or spatial calculation models) and summing the BIC values of all the Voronoi cells obtained from each tessellation and using this as the basis for judging whether this spatial calculation model is the optimal model; and (3) selecting, among the 3900 spatial calculation models, 100 models (marked as best-100) with smaller BIC values as the optimal models and using the parameters [β, μ, σ] of the ensemble median (Q₂) and median absolute deviation (MAD) as the final calculation results. The b value can be obtained by $b=\mathit{\beta }/\ln \mathrm{10}$.

The maximum likelihood calculation of the OK1993 model parameter is not performed when the number of earthquakes contained in a Voronoi cell is N₁<5, so the actual number of effective cells N_v obtained by each tessellation is used to distinguish the number of randomly thrown nodes. Although the value of N₁ may affect the parameter fitting error in some polygons with a small number of events, considering that the OK1993 model in the form of a continuous distribution function has the advantage of obvious fit adaptability compared to the traditional linear frequency–magnitude distribution (FMD) function in a small number of data cases, this setting also ensures that the spatial division can obtain more polygon calculation results, and the final result of the parameters is expressed by the ensemble median value, so the effect of this method of value taking on the final result is minimal.

In the above calculation steps, the setting of the maximum number of nodes, the number of random throws, etc. has obvious subjectivity. However, due to the fact that the data-driven approach actually obtains a very stable final result when the number of divisions and the number of grid nodes are sufficient (Si and Jiang, 2019) – for example, when the maximum number of nodes is 100, each type of node is randomly thrown 1000 times – the final result obtained when 1000 optimal models are selected is almost the same as the result of this paper.

3 Study region and data used

The 2019 Changning M_S=6.0 earthquake sequence occurred in the basin–mountain junction in the southern margin of the Sichuan Basin, where the tectonic activity is relatively weak. The seismicity in the area is mainly controlled by folds and associated faults. The intensity of historically destructive earthquakes is low in the area where aftershocks extend. No earthquake with a magnitude above 5.0 had been recorded in this area before the Changning M_S=6.0 earthquake. According to Yi et al. (2019), it is inferred that the occurrence of the Changning M_S=6.0 earthquake sequence may be related to the Baixiangyan–Shizitan anticline and the Shuanghechang anticline and their associated fault activities. Figure 1 shows the study area of this paper. We will focus on the rectangular area A^′B^′C^′D^′ where the aftershock sequence mainly occurred and the rectangular area ABCD where the surrounding earthquakes are active.

Figure 1 Distribution of seismicity in the Changning area. The red dots show the aftershocks of the Changning M_S=6.0 earthquake, and the green dots indicate the earthquakes that occurred before the Changning M_S=6.0 earthquake. Hexagonal stars mark the position of the mainshock and four aftershocks with magnitude no less than 5.0, and the corresponding focal mechanisms are marked. The dotted rectangular ABCD and A^′B^′C^′D^′ show the two spatial regions for calculating the b value and rotating the coordinate system, and the blue cross symbol gives the origin where the coordinate system is rotated. The blue triangles show the location of seismic stations that recorded these earthquakes, and the solid black lines represent active faults (He et al., 2019). The study region is shown in the location figure in the bottom left by a green rectangle.

We used earthquake catalogs and bulletins provided by the Sichuan Regional Seismic Network from 1 January 2009 to 17 July 2019. To obtain relatively reliable parameters such as the epicenter location and focal depth, the double-difference algorithm hypoDD (Waldhauser and Ellsworth, 2000) was used to relocate the earthquakes. Among the data we used are a total of 21 246 seismic events that meet the requirements of the hypoDD method with not less than four arrivals, including 516 649 P-wave arrivals and 506 809 S-wave arrivals, and 59 permanent seismic stations and temporary seismic stations are used which are located in Sichuan and the surrounding provinces. We used a 12-layer one-dimensional crustal velocity model (Xie et al., 2012) during the relocation. The ratio of V_P to V_S is set to 1.730.

A total of 18 371 earthquake events were relocated (Fig. 1), of which the smallest event had a magnitude of −1.0. Among them were 13 728 and 4642 earthquakes before and after the M_S=6.0 mainshock, respectively. The horizontal and vertical uncertainties are 0.425, 0.457, and 0.654 km, respectively. The average root mean square (rms) of the travel-time residuals of the locations was reduced to 0.162 s. In addition, a total of 2875 events were discarded, accounting for 13.53 % of the number of earthquakes in the original catalog, and most of their magnitudes ranged between M_L 0.3 and M_L 1.4. Considering that the data-driven approach used in this paper is the selection and ensemble averaging of a large number of random space partitioning schemes and that the OK1993 model is a continuous function of the magnitude–frequency distribution, the effect of these excluded events on the calculation result of the b value can be ignored.

From the spatial distribution of the relocated earthquakes shown in Fig. 1, the aftershocks are mainly distributed in the northwest direction of the mainshock epicenter and extend along the Changning anticline with a length of about 27 km, which is much longer than the rupture scale of about 10 km for a M=6 earthquake in accordance with the empirical formula given by Wells and Coppersmith (1994). Besides, the shape of the aftershock distribution is not simply linear; there are obvious inflections in the middle segment, and in the northwest there is a branch approximately perpendicular to the direction of the aftershock distribution. There are relatively few aftershocks near the epicenter of the mainshock, and a large number of aftershocks occurred in the northwest.

In the aftershock sequence of the Changning M_S=6.0 earthquake, there are four aftershocks with magnitudes exceeding M_S=5.0, which are 17 June 2019 M_S=5.1, 18 June 2019 M_S=5.3, 22 June 2019 M_S=5.4, and 4 July 2019 M_S=5.6 earthquakes.

To facilitate the calculation of b values and the display of the results, we have selected only the events within the rectangular area A^′B^′C^′D^′ where almost all aftershocks are concentrated and the rectangular area ABCD where a large number of earthquakes had taken place before the mainshock occurred. The positions of these earthquakes were transformed by Cartesian coordinates and rotated according to the origin point (28.395^∘ N, 104.986^∘ E) of the coordinates so that the aftershock sequence can be spread horizontally in the new coordinate system. The epicenter distribution after coordinate transformation in Fig. 2a–c shows the spatiotemporal distribution on the distance versus rank of index 2-D map of the earthquake within the rectangular frame A^′B^′C^′D^′.

Figure 2 Distribution of seismicity for b-value calculations. (a) Rotation of the coordinate system to the seismic distribution along the direction of the aftershock distribution; (b) projection of the earthquakes in the rectangular frame A^′B^′C^′D^′ on the depth profile; (c) the temporal and spatial distribution on the distance versus rank of index 2-D map of the earthquakes within the rectangular frame A^′B^′C^′D^′. The meaning of the symbols is the same as in Fig. 1.

4 Spatial distributions of b values in surface and depth profiles

According to the technical process of the data-driven approach described above, after Voronoi tessellation, calculation of the BIC values, and selection of the optimal models, the ensemble median (Q₂) and ensemble median absolute deviation (MAD) of b values can be obtained. Figure 3 shows an example of calculating the parameters of the OK1993 model in terms of the frequency–magnitude distribution based on a data-driven approach. Figure 3a is the distribution of those BIC values corresponding to the number of effective cells N_V, and the red dots are the selected best-100 models. Figure 3b shows an example of the best-100 models, that is, in the case of N_V=20, the Voronoi tessellation in the rectangular study area ABCD and the distribution of b values obtained by its calculation. Figure 3c shows an example of the fitting result of the OK1993 model corresponding to a cell in Fig. 3b. The OK1993 model parameters obtained by the fitting are b=0.714, μ=0.376, and σ=0.247.

Figure 3 An example of calculating the parameters of the OK1993 model in terms of the frequency–magnitude distribution based on a data-driven approach. (a) Distribution of BIC values versus the number of effective cells N_v in the Voronoi tessellation. The black dots and error bars are commensurate with the mean value and 1 standard deviation of BIC values under the corresponding N_v, respectively. The top horizontal dashed line marks the BIC values of the entire spatial region without mesh generation (BIC₀, N_v=1). The red dots show the BIC values with the best-100 solutions selected, while the gray dots are the other BIC results according to N_v. (b) Example of Voronoi tessellation of N_v=20 and one of the best-100 models selected. The white circles are the positions of the Voronoi nodes, and the resulting partitions are color coded by their estimated b values (obtained from the β value in the OK1993 model). (c) Example of fitting result for the frequency–magnitude distribution (FMD) of the OK1993 model in the Voronoi cell indicated by a thick line in graph (b).

We calculated the distribution of the ensemble median b value in the rectangular region ABCD and the depth profile of the rectangular region A^′B^′C^′D^′. The results are shown in Fig. 4. Figure 4a and b are the result before the Changning M_S=6.0 earthquake and the entire study period, respectively. The results show that the b values exhibit a strong heterogeneous spatial distribution in the rectangular region ABCD before the Changning M_S=6.0 earthquake. Low b values are mainly distributed in the eastern half of the area, with its lowest value being b=0.732 and located near the epicenter of the mainshock. Low-b-value contours are mainly distributed in the NE–SW direction and are consistent with the direction of the Shuanghechang anticline and its associated faults passing through the main epicenter. In the western part of the rectangular region ABCD, high b values are distributed, with a largest value of b=2.200. This indicates that before the Changning M_S=6.0 earthquake, the differential stress near the epicenter of the mainshock was high, but the spatial scale of this larger differential stress was much smaller than the scale of the aftershock spatial distribution. The spatial distribution of b values calculated using all seismic events (see Fig. 4b) shows that the area with low b values in the region ABCD is significantly enlarged, and the b values in the rectangular region A^′B^′C^′D^′ are almost less than 1.0 and further reduced to 0.698 near the epicenter of the mainshock. This phenomenon of a significant decrease in the b value of the aftershock sequence after the mainshock widely exists in many earthquake cases (El-Isaa and Eatonb, 2014; Gulia and Wiemer, 2019).

Figure 4 The spatial distribution of the ensemble median b values of the best-100 solutions for N_v=2–N_v=40 in the Changning area. (a) The ensemble median b values before the Changning M_S=6.0 earthquake is distributed on the horizontal plane after the rotation; (b) the ensemble median b values obtained by calculation of the whole earthquake including the aftershocks of the Changning M_S=6.0 earthquake distributed on the horizontal plane after the rotation; (c) distribution of the ensemble median b values before the occurrence of the Changning M_S=6.0 earthquake in the rectangular frame A^′B^′C^′D^′ on the depth profile; (d) distribution of ensemble median b values obtained by calculation of all earthquakes including aftershocks of the Changning M_S=6.0 earthquake in the rectangular frame A^′B^′C^′D^′ on the depth profile. The black dots on each graph mark the seismic events used in the calculation.

Figure 4c and d show the distribution of the ensemble median b value on the depth profile of the rectangular area A^′B^′C^′D^′ and correspond to the results before the Changning M_S=6.0 earthquake and all study periods, respectively. The calculation results after considering the depth information of the earthquake show that b values also have strong heterogeneity at different depths. Among them, in Fig. 4c, low b values are mainly distributed at a depth of 4–15 km, which contains the source of the Changning M_S=6.0 earthquake and the 17 June 2019 M_S=5.1 earthquake. The lowest b value is about 0.493, which is much smaller than the minimum value in Fig. 4a. In Fig. 4d, considering the occurrence of the Changning M_S=6.0 earthquake sequence, the distribution area of low b values expands in the NW direction, and the lowest b value is about 0.501, which is close to that in Fig. 4c. Compared with the results obtained by ignoring the depth information of the earthquake in Fig. 4a and b, the results obtained by Fig. 4c and d reveal more significant heterogeneity of b values when we extend the investigation to different depths of the crust. Lower b values may indicate that there should be greater differential stress at the depth where the source area of the mainshock is located, and it is easily ignored by b-value calculations that usually do not consider the depth information of earthquake events.

Figure 5 shows the spatial distribution of the median absolute deviation (MAD) of b values by the data-driven approach according to Fig. 4. The ensemble MAD b value is smaller in most regions of Fig. 5a–d, especially in the rectangular region A^′B^′C^′D^′, which implies that these regions have relatively stable distribution and reliable ensemble median b values. As a comparison with Figs. 4 and 5, we also used the Changning M_S=6.0 earthquake and aftershocks to calculate the ensemble median b values and MAD b values. For the corresponding results, please see Fig. S2.

Figure 5 The spatial distribution of the median absolute deviation (MAD) of the b values by the data-driven approach according to Fig. 4. (a) The ensemble MAD b values before the Changning M_S=6.0 earthquake is distributed on the horizontal plane after the rotation; (b) the ensemble MAD b values obtained by calculation of the whole earthquake including the aftershocks of the Changning M_S=6.0 earthquake distributed on the horizontal plane after the rotation; (c) distribution of the ensemble MAD b values before the occurrence of the Changning M_S=6.0 earthquake in the rectangular frame A^′B^′C^′D^′ on the depth profile; (d) distribution of ensemble MAD b values obtained by calculation of all earthquakes including aftershocks of the Changning M_S=6.0 earthquake in the rectangular frame A^′B^′C^′D^′ on the depth profile.

5 Spatiotemporal heterogeneity of b values

Considering that the b value usually changes over time before and after a strong earthquake, this paper not only examines the spatial distribution of b values in the surface and depth profiles but also discusses the spatiotemporal distribution of b values for earthquake events in the rectangular area A^′B^′C^′D^′ where the Changning M_S=6.0 sequence is located. Due to the strong temporal and spatial inhomogeneity of seismic activity, especially clustering in time, there are great difficulties in obtaining a stable and reliable b value and clearly showing the temporal and spatial variation in the b value. In order to reduce this difficulty to a certain extent, here we use the index of earthquake occurrence instead of time; that is, the earthquake is projected on a pseudo-time axis of the index number of the occurrence time sequence. Using the same calculation method as in Figs. 4 and 5, the distributions of ensemble median b values and ensemble MAD b values on the distance–index map are obtained. The corresponding results are shown in Figs. 6 and 7. Considering the possible abrupt change in the regional stress field due to strong earthquakes such as the Changning M_S=6.0 earthquake, we adopt two schemes to study the spatiotemporal distribution of b values. One is to study the seismicity before and after the mainshock as a whole, and the other is to study the seismicity before and after the mainshock as two independent periods. The calculation results under the two schemes are shown in Fig. 6a and b, respectively.

Figure 6 Spatiotemporal distribution of the ensemble median b values of the best-100 solutions for N_v=2–40 in a 2-D space consisting of the distance along strike and rank of index. (a) The ensemble median b values obtained from all data before and after the Changning M_S=6.0 earthquake; (b) the ensemble median b values obtained from the data before and after the Changning M_S=6.0 earthquake considered separately. The vertical dotted line shows where the M_S=6.0 earthquake occurred. The timescale is marked on the x axis (above the axis), including the time of whole years marked by long ticks and the half-year time marked by short ticks.

Figure 7 Spatiotemporal distribution of the median absolute deviation (MAD) of the b values of the best-100 solutions for N_v=2–40 in a 2-D space consisting of distance along strike and rank of index. (a) The ensemble MAD b values obtained from all data before and after the Changning M_S=6.0 earthquake; (b) the ensemble MAD b values obtained from the data before and after the Changning M_S=6.0 earthquake considered separately. The vertical dotted line shows where the M_S=6.0 earthquake occurred. The timescale is marked on the x axis (above the axis), including the time of whole years marked by long ticks and the half-year time marked by short ticks.

It can be seen in Fig. 6a that before the Changning M_S=6.0 earthquake, in the segment between −5 and −10 km near the A^′–B^′ end and with a length of about 10 km (NW direction of the aftershocks in Fig. 1), there were relatively stable high b values, with the maximum value exceeding 2.0. In the segment between −5 and 12 km near the C^′–D^′ end and with a length of about 17 km (the SE direction of the aftershocks in Fig. 1, including the nucleation point of the mainshock), there were relatively stable low values before the Changning M_S=6.0 earthquake, and the range of the low b values gradually narrowed and became concentrated towards the nucleation point of the mainshock. After the Changning M_S=6.0 earthquake occurred, the b values in the entire spatial range from A^′–B^′ to C^′–D^′ decreased significantly. Among them, the b values in the 0–12 km segment where the nucleation point of the mainshock is located recovered rapidly, while the b values in the 0 to −15 km segment increased at a slower rate.

From the results before and after the Changning M_S=6.0 earthquake shown in Fig. 6b, it can be seen that the occurrence of the mainshock has a greater impact on the continuity of time-variant b values. This means that the spatiotemporal evolution image of the b values given in Fig. 6a over the entire study period is not physically valid. Correspondingly, the decrease in pre-mainshock b values and the sudden expansion of the low b values may be a kind of artifact caused by the subsequent aftershocks brought into the calculation (Lei et al., 2019).

Compared with Fig. 6a, the results in Fig. 6b show that before the Changning M_S=6.0 earthquake, the shortening and concentration changes for the low-b-value segment near the C^′–D^′ end and the expansion process of the high-b-value segment near the A^′–B^′ end are performed simultaneously. This implies that a significantly higher differential stress area is concentrated towards the nucleation point of the mainshock. Figure 7a and b show the distribution of ensemble MAD b values according to Fig. 6a and b, where higher ensemble MAD b values mainly appear in some areas with higher b values in Fig. 6a and b.

6 Discussion

In the pattern of b-value spatial heterogeneity before strong earthquakes, the locations of rupture nucleation points, sliding distributions, and aftershock distributions of some strong earthquakes were observed to correspond to areas with lower b values, such as the Parkfield M=6.0 earthquake on 28 September 2004 (Wiemer and Wyss, 1997; Schorlemmer et al., 2004; Schorlemmer and Wiemer, 2005). However, the significant spatial heterogeneity of b values obtained from the studies of these earthquakes is suspected to be related to the subjective arbitrariness of the calculation rules (Kamer and Hiemer, 2015). The calculation results based on the data-driven method (Si and Jiang, 2019) in this paper show that significant spatial heterogeneity of b values can still be observed before the Changning M_S=6.0 earthquake, especially on the depth profile of the fault. Moreover, according to the empirical relationship between the magnitude and rupture scale of Wells and Coppersmith (1994), the low-value spatial scale of b<0.75 in Fig. 4c is also close to the rupture length of about 10 km for the M=6.0 mainshock. This also means that it is still feasible to use the spatial heterogeneity of the b values to identify the locked asperities and determine the location of future strong earthquakes if more cases are verified.

There is still much controversy over the temporal variation pattern of b values in the source area before a strong earthquake. Although the decrease in b values prior to failure was found in laboratory fracturing experiments on relatively complete rock samples (e.g., Thompson et al., 2006; Lei, 2019) and in the case studies of strong earthquakes (Nanjo et al., 2012; Schurr et al., 2014; Bayrak et al., 2017; Huang et al., 2020), a large number of reported temporal variations in b values before actual strong earthquakes are still considered to have no statistically significant predictive power (Parsons, 2007). Some studies have found that the temporal variation in b values corresponding to asperities is synchronized with the loading rate and shear stress (Tormann et al., 2013). Schorlemmer et al. (2004) and Wiemer and Wyss (2002) studied some earthquake cases and concluded that the b value is quite stable over time and it is difficult to observe a significant change. The study of the relationship between acoustic emission events and stress in the stick-slip experiment shows that the complexity of the temporal variations in b values observed when sliding on rough fault planes may be due to fault-structure heterogeneity (Goebel et al., 2013). In this study of the Changning M_S=6.0 earthquake, we did not simply examine the temporal variations in b values in a fixed spatial range but investigated the migration pattern of the b value in a 2-D spatiotemporal dimension. We found that as the time approaches the occurrence of the mainshock, the spatial range of the low b values gradually shrinks and focuses on the vicinity of the rupture nucleation point, and the b values do not decrease significantly. Under the assumption that the fault-structural heterogeneity will not change in the short term and based on previous understandings of the correlation between high b values and fluid-induced seismicity, the migration pattern in this paper may be explained by the erosion of fluid in the high-differential-stress area where the nucleation point is located.

For the spatiotemporal heterogeneity of the b value of the aftershocks of the 2019 Changning M_S=6.0 earthquake, we noticed that the aftershocks expanded spatially to areas with high pre-mainshock b values in the northwest direction, and the length of the aftershock area was significantly longer than the rupture scale of the earthquake (see Fig. 6b). Since the aftershocks do not exhibit relatively slow spatiotemporal migration behavior, the physical mechanism that drives the aftershocks of this earthquake cannot be explained by either the traditional stress corrosion model (Das and Scholz, 1981) or the frictional afterslip model (Perfettini et al., 2018; Koper et al., 2018). Some views suggest that aftershock activity in high-b-value regions may be related to the reactivation of highly fractured fault zones, the redistribution of stress fields, and the role of fluids trapped in microfractures (Aktar et al., 2004). Long et al. (2020), imaging the velocity structure of the area where the Changning M_S=6.0 earthquake was located, showed that there is an obvious S-wave low-velocity anomaly at the depth of 3 to 8 km in the northwestern segment of the aftershock. In this paper, this S-wave low-velocity anomaly region also corresponds to the distribution of high b values, which may be related to fluid intrusion. Therefore, we deduce that the abundant aftershocks produced by this mainshock and the active area that exceeds the rupture scale of the mainshock are more likely to be caused by the mainshock which triggered a series of complex structural aftershocks northwest of the nucleation point. The dynamic expansion of the region of high pre-mainshock b values to the nucleation point also creates conditions for the triggering of a large number of aftershocks that are widespread spatially.

In addition, b values of the aftershocks first dropped rapidly to about 0.5, then gradually recovered, and returned to the pre-seismic level after the fourth aftershock of magnitude 5 (excluding high-b-value areas). The phenomenon that the b values of the aftershock sequence decrease immediately after the mainshock and then rapidly recover has been observed in many earthquake cases (El-Isaa and Eatonb, 2014; Tormann et al., 2015). Unlike most aftershock sequences, where the b value generally increases by 20 % after the mainshock, this sudden decrease in the b value is considered to be related to the occurrence of subsequent strong aftershocks or larger earthquakes (Gulia and Wiemer, 2019). In the aftershock sequence of the Changning M_S=6.0 earthquake, the rapidly decreasing b value of the aftershocks was accompanied by four strong aftershocks with magnitudes greater than 5.0, which is consistent with the phenomenon revealed by previous studies. This may also support the idea of discrimination between foreshocks and aftershocks by real-time monitoring of the b value in aftershock sequences (Gulia and Wiemer, 2019). However, it needs to be pointed out that similarly to the problem of sudden changes in the spatiotemporal distribution of b values before and after the mainshock, it cannot be ruled out that four strong aftershocks with M>5 will affect the continuity of the b values to a certain extent.

7 Conclusions

To reveal whether there is spatiotemporal heterogeneity of b values before and after the 2019 Changning M_S=6.0 earthquake and to overcome the subjectivity of the choice of data used for calculation, we applied the OK1993 model of magnitude–frequency distribution according to the data-driven idea to calculate b values. We also investigated the distribution characteristics of b values in three different ways: horizontal surface distribution, depth profile distribution, and the distance–rank of index map. The main conclusions are as follows:

1. The b values before and after the Changning M_S=6.0 earthquake showed strong spatiotemporal heterogeneity in the horizontal surface distribution, in the depth profile distribution, and on the distance–rank of index map. Among them, before the Changning M_S=6.0 earthquake, there were obvious low-b-value distributions near the epicenter of the mainshock and within the depth range of 3 to 12 km. The correlation shows that there may be significantly higher differential stress in the source area before the Changning M_S=6.0 earthquake. The northwestern segment of the aftershocks has a b-value distribution that is distinctly high, which coincides with the S-wave low-velocity anomaly region shown by the velocity structure imaging.

2. The b-value spatiotemporal distribution results show that before the Changning M_S=6.0 earthquake, the high-b-value region of the NW segment spread by aftershocks gradually expanded and approached the nucleation point as the time approached the failure time of the mainshock. This may be related to the fluid intrusion in the rock. A large number of aftershocks were produced, and the area where the aftershocks were spread was significantly larger than the rupture scale of the mainshock. The mainshock may have triggered seismicity in the NW direction where the fluid intrudes.

3. The b values of the aftershocks of the Changning M_S=6.0 earthquake decreased rapidly and gradually recovered after the mainshock, indicating a higher differential stress level in the aftershock area. The time variation in low b values is synchronized with the occurrence of strong aftershocks with M≥5.0, showing the application potential that can be used to distinguish between foreshocks and aftershocks.

4. Although the distribution characteristics of b values before and after the Changning M_S=6.0 earthquake were qualitatively consistent when they were studied in different space-time dimensions, there were significant differences in specific b values. For example, the minimum b value of the Changning M_S=6.0 earthquake in the depth profile distribution is about 0.493, but it is about 0.732 when the seismic depth information is ignored and the value is only calculated on the surface. This inconsistency needs special attention when studying the spatiotemporal heterogeneity of b values.