A New b Value Estimation Method in Rock Acoustic Emission Testing

The Gutenberg‐Richter (G‐R) relationship has been regarded as the fundamental description of the size distribution from large‐scale‐earthquake to small‐scale‐laboratory rock ruptures. However, because deviation from the G‐R relationship has often been reported, especially when the amplitude distribution is used for b value estimation in rock acoustic emission testing, the effect of attenuation should be considered. Here, we perform a detailed analysis on the deviation of the size distribution from a G‐R law and discuss modification of the attenuation on the doubly truncated size distribution from a more general perspective. We find that the attenuation will not modify the size distribution, and the b value is theoretically verified to be unchanged within an interval specified by the minimum and maximum amount of attenuation. Based on these discussions, we propose a new b value estimation procedure for rock acoustic emission testing and apply it to a dilation rupturing test; the results confirm that b value depends on material heterogeneity and stress.


Introduction
In seismology, source energy E and seismic moment M 0 are the fundamental physical measure of seismic event; M 0 can be transformed to moment magnitude M w by a logarithmic relation to get uniform and consistent magnitudes for earthquake catalog. As a widely used single parameter to characterize the earthquake, the local magnitude M can also be scaled by M w , and the proper scaling between M and M w is an important prerequisite for any seismic hazard assessment (Deichmann, 2006(Deichmann, , 2017. The source energy E and seismic moment M 0 of earthquakes commonly follows a power law size distribution which is an intrinsic characteristic of the frequency-size distribution for earthquake, while the exponential Gutenberg-Richter (G-R) law in magnitude M as a proxy source parameter is equivalent to the power law in a true source parameter. The G-R law is expressed on a logarithmic scale given by where N is the number of earthquakes with a magnitude ≥M, a and b are constants (Gutenberg & Richter, 1944), and magnitude M is a logarithmic measure of the radiated energy or the seismic moment of the source. Here, the parameter b describes the size distribution scaling, which is often referred to as the b value, and the spatial and temporal variation of the b value is an important parameter for seismic hazard assessment.
As b value is a statistical parameter, its inferred value is influenced by many factors, and the errors in estimation can result in potentially misleading b value variations (Roberts et al., 2015). One of the most important reasons for errors in b value estimation is the fact that the size distribution is not usually log linear in the whole range of magnitudes, and sometimes it follows neither the logarithmic linear nor any smoothly nonlogarithmic linear law but can be more complex (Lasocki & Papadimitriou, 2006;Vorobieva et al., 2016). In fact, it is often observed that changes in the size distribution for larger earthquakes can be explained by finite size effects (Casertano, 1982;Clauset, 2007) and statistical artifacts (Main, 2000), and those for smaller events are caused by incomplete reporting in earthquake catalogs (Cao & Gao, 2002;Kwiatek et al., 2014;Mignan, 2012;Raub et al., 2017;Wiemer & Wyss, 2000). Meanwhile, errors in b value estimation are also related to the data volume. It is considered in some researches that a size of 50 events can be adopted as the minimum number of events for stable b value estimation (Kurz et al., 2006;Murru et al., 1999;Schorlemmer et al., 2004); however, in some other studies, a size of 200 events was accepted as the minimum data volume when considering the effect of estimation uncertainties (Amorese et al., 2010;Roberts et al., 2015). Other factors such as bin width, the cumulative or incremental data counting, and the estimation algorithm will all influence the b value result (Aki, 1965;Amorese et al., 2010;Bender, 1983;Greenhough & Main, 2008;Nava et al., 2016;Weichert, 1980). Acoustic emission (AE) generated by fracturing of rock has been found to obey the same form of size distribution, and there have been numerous studies indicating that in many ways the mechanism of earthquake foreshock sequences could therefore be reproduced through AE from the statistical behavior of microfracture activity observed in laboratory fracture experiments Lei, 2003;Lockner, 1993;Mogi, 1962;Scholz, 1968Scholz, , 2015Vorobieva et al., 2016). Moreover, the variation of b value obtained in rock AE deformation tests and earthquakes has always been compared to establish analogies in the damage process and precursory analysis (Goebel et al., 2013;Lennartz-Sassinek et al., 2014;Lockner et al., 1991;Main et al., 1989;Scholz, 1968), and specifically the effects of pore fluid pressure on b value and precursors have also been focused through laboratory experiments (Sammonds et al., 1992). Because of the acquisition threshold and the limitation of the maximum output voltage in some AE equipment, a similar deviation from the G-R relationship of the size distribution will also appear in laboratory rock AE tests, and the influence of data volume, cumulative or incremental data counting, bin width, and estimation algorithm on the b value result is obvious (Cosentino et al., 1977;Liakopoulou et al., 1994;Lockner, 1993).
Furthermore, it is worth noting that unlike the b value in earthquakes which uses the frequency-magnitude distribution for estimation, the apparent frequency-amplitude distribution is mostly used for b value estimation in rock AE tests (Cox & Meredith, 1993;Liakopoulou et al., 1994;Weiss, 1997), and this apparent amplitude is directly derived from the sensor-collected signal that is attenuated from the source. Consequently, the distribution of these apparent amplitudes is not the physical size distribution of the sources, so attenuation could also modify the AE b value (Lavrov, 2005;Lockner et al., 1991;Unander, 1993;Weiss, 1997). Although some researchers have made provisions when using apparent amplitude to estimate b value by adjusting the AE amplitude for 1 r geometric spreading and averaging over transducers to give an equivalent event amplitude (Lockner et al., 1991) or by using the root-mean-square principle to obtain a relative AE magnitude (Kwiatek et al., 2014), the amplitude distribution actually cannot fully represent the crack size distribution. The attenuation effect was first proposed by Lockner et al. (1991), and then the attenuation effect on b value has been theoretically discussed by Unander (1993) and Weiss (1997). They assumed a constant attenuation coefficient and uniformly distributed sources in the medium, and they considered that the attenuation is composed of a geometric term and a "dissipative" exponential term. In fact, the attenuation coefficient is not constant in the rock material. Not only is it a frequency-dependent parameter but it also varies with the direction of elastic wave propagation because of the anisotropy of the rock, making it extremely difficult to theoretically describe the attenuation in the general case.
Since estimation result of b value can be influenced by many factors especially the deviation of frequency-magnitude distribution from G-R relationship. Therefore, a reliable b value estimation procedure is critical in rock AE testing which can well ensure the robustness of comparative precursory analysis between laboratory AE tests and earthquakes. Here we perform a detailed analysis of the deviation of the size distribution from G-R relationship and intend to examine the influence of statistical method on size distribution measurement. Meanwhile, we investigate the modification of attenuation on doubly truncated size distribution from a more general perspective, where we assume that sources in any distribution are attenuated to a certain extent and amplitudes are measured on the sample boundary by acoustic sensors, then location of the sources and sensors and the amplitudes and locations of the sources are independent. Based on the deviation discussion, we propose a new b value estimation procedure which specify the minimum data volume and statistical method and employ Fisher optimal split and Global Search algorithm (here we named this procedure as FGS) to determine the logarithmic linear segment in the frequency-amplitude distribution. In order to verify the reliability of FGS method, we design a static dilation rock rupturing AE test by injecting nonexplosive cracking agent into three predrilled boreholes to form a specific fracture surface in cubic rock specimen (as shown in Figure 1). This experimental design is to ensure that the sensor-collected AE signals are all generated by expansion rupturing in rock specimen and does not rely on source location to identify valid rupturing data. The results of this paper are of great help for deep understanding of attenuation effect on b value which show that the inferred b value for AE apparent amplitude is equal to that at source in a finite interval, and the data truncation actually is the fundamental cause of attenuation modification of the size distribution. Furthermore, the proposed FGS method in this paper also provides a new way for b value estimation through apparent amplitude in rock AE tests.

Deviation From the G-R Law of Frequency-Size Distribution
The G-R relationship of the size distribution is an intrinsic characteristic of earthquakes and rock AE events that indicates that the number of small-scale ruptures is much greater than the number of large-scale ones. The G-R relationship can be interpreted as a manifestation of the self-organized critical behavior of earth dynamics (Bak & Tang, 1989;Christensen et al., 2002;Clauset, 2007;Newman, 2005). However, deviation from the G-R relationship has always been observed in earthquake and rock AE experiments (Cox & Meredith, 1993;Liakopoulou et al., 1994;Lockner, 1993;Main, 2000;Weiss, 1997). One of the reasons for this deviation is inadequate data acquisition. At the lower magnitude end, this is apparent from the difficulty in determining the magnitude of completeness, Mc (Cao & Gao, 2002;Godano et al., 2014;Mignan, 2012;Radziminovich et al., 2019;Raub et al., 2017). Meanwhile, small-scale ruptures can be masked by larger ones, which is also a very significant factor for inadequate data acquisition in rock AE experiments when avalanche destruction occurs during the final loading stage or in high strain rate loading tests (Cox & Meredith, 1993;Liakopoulou et al., 1994;Liu et al., 2015Liu et al., , 2020Meredith & Main, 1990).
Another reason for this deviation is due to the statistical methods used in b value estimation. In fact, because of the accuracy of equipment acquisition and the finite size of earthquake and AE amplitudes, the b value is estimated in a doubly truncated exponential distribution (Cosentino et al., 1977;Page, 1968). Generally, the density function of a doubly truncated distribution generated by an underlying exponential one can be expressed as follows: where α and β are constants, β= bln (10), and "other" is any other distribution that is excluded from consideration. Then, the cumulative frequency distribution in the interval where N total is the total number of events in the catalog. If the cumulative frequency distribution is expressed again on a logarithmic scale, the probability function becomes Accordingly, it is found from this equation that when x is closer to the upper limit of the magnitude a x , x is not linearly related to log 10 N(x), and the logarithmic cumulative frequency distribution behavior can be significantly altered. However, if the incremental statistical method is used, the incremental frequency distribution in the interval is And again the probability function on a logarithmic scale will be which shows that, for a certain bin width, x i − 1 is linearly related to log 10 N(x i − 1 ). The above analysis means that, theoretically, the cumulative frequency distribution used for b value estimation will inevitably result in deviation from an exponential distribution while the incremental frequency distribution will not.
In addition, owing to the natural smoothing effect of plotting cumulative frequency data, a regression analysis based on the cumulative frequency distribution will systematically increase the goodness of fit and affect the computation of the magnitude of completeness, Mc (Amorese et al., 2010;Main, 2000;Schorlemmer et al., 2004;Wiemer & Wyss, 2000). Moreover, the inaccurate magnitude determinations and breaks in a logarithmic frequency-magnitude distribution slope will possibly be smoothed out in cumulative frequency distribution, which cause the estimated b value to incorrectly describe the size distribution scaling in practice.

Effect of Attenuation on AE Frequency-Amplitude Distribution
The attenuation of elastic waves in rock material is hard to describe theoretically, so it is difficult to analyze the modification of attenuation on size distribution. Here we perform an analysis from another more general perspective to investigate the effect of attenuation on b value by assuming that sources in any distribution are attenuated to a certain extent and amplitudes are measured on the sample boundary by acoustic sensors independent of the composition of the attenuation terms.

Continuous Conditional Probability Density Function of the Attenuated Source Amplitude
Suppose the AE source obeys any distribution in a bounded region Ω (Ω maybe a 3-D solid, surface, or curve) and assume that the b value and the region remain unchanged in some given time interval discussed. Let the source decibel amplitude (Z) probability density function in exponential distribution which truncated at an interval [a 0 , a x ] be as follows: where α and β are constants and "other" means any other distribution that is excluded from consideration.
We define (x 1 , y 1 , z 1 ) in Ω and the origin (0, 0, 0) as the coordinates of the source and sensor, respectively. In addition, the attenuation amplitude can be summarized by the relation A =A 0 g(x 1 , y 1 , z 1 ) , where 0 < g(x 1 , y 1 , z 1 ) < 1, g(x 1 , y 1 , z 1 ) is the percentage attenuation in some given manner and g can be an any expression, A and A 0 are the voltage amplitude of the AE source after and before attenuation respectively, and the attenuation relation on a logarithmic scale is Letting X = 20log 10 (A/10 −6 ) and Z = 20log 10 (A 0 /10 −6 ) be the decibel amplitudes, we can simplify Equation 8 as X = Z+20log 10 g(x 1 , y 1 , z 1 ).
The density function of the coordinate of the source is f (x 1 , y 1 , z 1 ), which is in any distribution. Suppose that the source amplitude Z and its location (x 1 , y 1 , z 1 ) are independent of each other, so that f Z (z)f (x 1 , y 1 , z 1 ) is the joint density function. The probability of the attenuated decibel amplitude X can be given by We assume the region Ω be a subspace of three-dimension space and Ω could be any given 1 to 3 dimensional bounded subspace. In order to facilitate the discussion, the integral of Ω is denoted as ∭ Ω dω. Then the probability of the attenuated decibel amplitude X can be calculated as follows: where γ is a constant.
Þas the maximum and minimum attenuation amplitudes from the source to sensor, respectively, so −20 log 10 g (x 1 , y 1 , z 1 ) − m ≥ 0 and −20 log 10 g(x 1 , y 1 , z 1 ) − M ≤ 0. When a 0 − m ≤ x ≤ a x − M, we have a 0 ≤ a 0 -− m − 20 log 10 g(x 1 , y 1 , z 1 ) ≤ x − 20log 10 g(x 1 , y 1 , z 1 ) ≤ a x − M − 20log 10 g(x 1 , y 1 , z 1 ) ≤ a x , so a 0 ≤ x − 20log 10 - Then we can continue to calculate the F X (x) as follows: The probability of the attenuated amplitude X in interval [a 0 − m, a x − M] can be given by ln 10 dω; Now we only consider the interval [a 0 − m, a x − M] of X; it means that we will calculate the conditional probability distribution of X known that X is in the interval [a 0 − m, a x − M]. In order to facilitate the discussion, X in the interval [a 0 − m, a x − M] is replaced by X 0 . Then, the conditional probability distribution function for the attenuated amplitude X 0 can be given by 10.1029/2020JB019658

Journal of Geophysical Research: Solid Earth
Then the conditional probability density function of X 0 is

Discrete Probability Density Function of the Attenuated Source Amplitude
Similarly, suppose AE source Z 2 obeys any distribution in a region Ω with exponentially distributed amplitude in decibels. If Y is the attenuation amplitude in decibels from source to sensor, then its probability can be given as where a i is a discrete attenuation value and a 1 < a 2 < ⋯ < a n . Then, the amplitude X after attenuation will be Suppose that the source amplitude Z 2 and Y are independent of each other, then the probability of the attenuated amplitude X can be given by The density function is Then, the conditional probability distribution function for the attenuated amplitude X 1 can be given by 10.1029/2020JB019658

Journal of Geophysical Research: Solid Earth
And the conditional probability density function of X 1 is As shown in Equations 7 and 14 or Equations 16 and 20, when the exponential source distribution is doubly truncated at a finite interval, the distribution function of the attenuated amplitudes obeys the same exponential distribution at a certain interval, and the β value is theoretically verified to be unchanged. This means that, if the amplitudes of AE sources obey an exponential distribution and the AE sources are in any distribution in a bounded region, attenuation will not affect the b value of the sources inferred from the AE amplitude data within a certain interval. That is, if the attenuated amplitude distribution measured on the sample boundary obeys the G-R relationship in a certain interval, its b value is the same as that of the AE source amplitude distribution.

A New b Value Estimation Procedure for Rock AE and Its Application
Based on the discussion above, we propose a new procedure when using apparent amplitude for rock AE b value estimation with the following steps: Step 1: Minimum number of events for b value estimation.
From the statistical point of view, the minimum amount of data for estimation does not have a universal set value. It can be selected after the level of significance is specified in practical applications. A commonly recognized axiom is that the power of a statistical test increases with the increase in sample size (Siegel, 1956). For logarithmic linear fitting of a G-R law distribution, the average error of the estimated scaling parameter becomes <1% when the sample size exceeds 50 (Clauset, 2007), and 50 events were also adopted as the minimum number of events for stable b value estimation in various studies (Kurz et al., 2006;Murru et al., 1999;Schorlemmer et al., 2004), while the stability estimation of b value needs to include the space-time window of 100 earthquakes, as illustrated by Shi and Blot (1982). Amorese et al. (2010) and Robert et al. (2015) also stated that the computation of b value includes at least~200 events. In comparison to the cumulative frequency distribution, larger data dispersion will occur if the incremental frequency distribution is used for b value estimation. Therefore, it is necessary to increase the data volume to ensure estimation robustness, so a minimum of 200 events should be utilized in the incremental frequency distribution.
Step 2: Perform incremental data counting and choose the bin width.
As discussed in section 2, the cumulative frequency distribution can inevitably deviate from a G-R relationship compared with the incremental frequency distribution, and the cumulative frequency distribution will also smooth the distribution breaks during the statistics process. Therefore, to display frequency distribution characteristics of AE data more realistically, the AE data should be counted incrementally.
Another parameter used for AE frequency-amplitude distribution statistics is bin width. In seismology, the selection of bin width is critical to grouping magnitudes. Improper choice of the bin width may introduce a significant bias in the magnitude of completeness and b value estimation (Marzocchi & Sandri, 2009;Schorlemmer et al., 2004;Wiemer & Wyss, 2000). Because earthquake magnitudes are given with one digit after the decimal point, many researches have demonstrated that it is reasonable to set the bin width as the magnitude round-off interval of 0.1 (Bender, 1983;Lasocki & Papadimitriou, 2006;Main, 2000). As AE amplitudes were measured in decibels, and the amplitudes are usually divided by 20 to produce a b value comparable to that reported in the seismic literature (Cox & Meredith, 1993;Liakopoulou et al., 1994;Sagar et al., 2012;Weiss, 1997), the divided value will have two digits after the decimal point with a 0.05 round-off interval. Therefore, the better bin width for amplitude grouping needs to be set to 0.05; this means that the number of amplitudes should be counted in steps of 1 dB. In fact, owing to the large amount of AE data collected in rock deformation tests, the coarser bin width will smooth the distribution breaks, while the finer bin width can demonstrate the frequency-amplitude distribution characteristics in more detail. Therefore, no matter what kind of equivalent amplitude or AE magnitude is used for b value estimation in the rock AE test, the bin width should be set to the size of the round-off interval at most.

Journal of Geophysical Research: Solid Earth
Step 3: Determine the logarithmic linear segment in the frequency-amplitude distribution.
As discussed earlier, the deviation of the size distribution from the G-R law will appear at both lower and upper amplitude ends. However, the analysis in section 3 demonstrated that, if the attenuated amplitude distribution is exponential in a certain interval, its b value is the same as that of the AE source amplitude distribution. Therefore, a specific method is needed for determining the logarithmic linear segment in the incremental frequency-amplitude distribution. Here, we propose a method based on the Fisher optimal split and global search algorithm; the main idea of this method is as follows: First, the slopes of any two successive points are calculated, and we want to find a benchmark slope (S(A i0 )) in the linear segment and an interval width centered on S(A i0 ). By referencing the idea of the 3σ principle of mean value in statistics, the relation between that interval width and the standard deviation will be established.
Since the points at both ends of the logarithmic frequency-amplitude distribution deviate significantly from the linearity, the standard deviation (std0) of all slopes is large, so only the standard deviation (std1) with slopes <0 is considered. We believe that as long as S(A i0 ) is found to be substantially accurate, the slopes in the linear segment should fall in the interval [S(A i0 ) − std1, S(A i0 )+std1]. In order to make this interval more accurate, the scale parameter r 0 is introduced. Then the slopes in the linear segment will fall in the interval [S(A i0 )−r 0 × std1, S(A i0 )+r 0 × std1], and S(A i0 ) is the slope corresponding to the maximum number of slopes that fall in the interval. In order to correct S(A i0 ) and the corresponding interval width, another scale parameter u is introduced, and a step size h is obtained from u and std0. The interval is adjusted to [S(A i0 ) − kh, S(A i0 )+kh], k = 1, 2, 3, 4, …, m. According to the change in the number of slopes falling in each interval at various k, Fisher optimal split algorithm (Fisher, 1958) is used to divide the number into two categories to obtain the interval  (Greenhough & Main, 2008); Fisher optimal split algorithm is used again to correct the interval, and the corresponding left and right endpoints of the reasonable linear segment can be obtained.
According to the instructions above, scale parameters r 0 and u which used to adjust internal width are set between 0.1-1 and 1-100, respectively. The fitting error variance of the linear segment calculated under the given initial values of r 0 and u is taken as the target value. Because any values of r 0 and u can get a target value, a global search algorithm is used to run 1,000 times to determine the optimal value of r 0 and u by finding the minimum error variance of generalized linear regression. Finally, the linear segment of the logarithmic frequency-amplitude distribution is screened out and the corresponding b value can be obtained by generalized linear regression. The specific algorithm of this method is demonstrated in the following seven steps: (1) Suppose that there are n points in the logarithmic frequency-amplitude distribution, A i − 1 and A i are the amplitudes of two successive points, and the slope for A = A i is defined as Then, the n−1 slopes are computed for each amplitude increment and the corresponding standard deviation is expressed by std0. Similarly, the standard deviation of all S(A i ) < 0 is expressed by std1.
(2) We define an interval [S(A i ) − r 0 std1, S(A i )+r 0 std1], where r 0 is a scaling parameter that is initially set to 0.1. For each slope where S(A i ) < 0, if the maximum number of slopes fall within that defined interval at a certain point i, then S(A i ) is replaced by S(A i0 ), which is defined as the benchmark slope, and the number of slopes falling within that defined interval at point i is marked as s. A step size h is defined as , and u is another scaling parameter that is initially set to 10.

Journal of Geophysical Research: Solid Earth
(3) Then, the interval defined in (2)  (4) The Fisher optimal split algorithm is used to divide {temp2 m − 1 } into two categories, each with a minimum sum of squared deviation. If k = k i is the optimal split point, the corresponding interval [S (A i0 ) − k i h, S(A i0 )+k i h] will be obtained. Because slopes less than zero are the ones to be considered, the interval can be adjusted to [S(A i0 ) − k i h, 0).
(5) To further optimize the benchmark slope S(A i0 ), a finer interval [S(A i0 ) − std2, 0), where std2 is the standard deviation of slopes falling in the interval [S(A i0 ) − k i h, 0), is used to remove the abnormal slopes at both ends. Then, a new benchmark slope for points within interval [S(A i0 ) − std2, 0) can be obtained using generalized linear regression.
(6) The new benchmark slope obtained in (5) is assigned to S(A i0 ), and the calculations in (3) and (4) are repeated to get the final optimal split point. Then the smallest amplitude and the largest amplitude in the interval [S(A i0 ) − k i h, 0) are the left and right endpoints, respectively, of the logarithmic linear segment.
(7) A global search algorithm is used to run 1,000 calculations to determine the optimal value of r 0 and u by finding the minimum error variance of the generalized linear regression, and finally, the logarithmic linear segment of the amplitude distribution is screened out.
Here we call this newly proposed b value estimation procedure the FGS method, and the complete workflow of FGS is shown in Figure 2.
To investigate the performance of FGS, a rupturing scheme using a nonexplosive expansion agent was designed (as shown in Figure 1) to conduct an AE test on granite, limestone, and red sandstone, which is intended to ensure that the sensor-collected AE signals are all generated by expansion rupturing in rock specimens. The AE activities were recorded by six  sensors with a resonant frequency of 140 kHz at a 5 MHz sampling frequency, and the signals collected by the sensor with largest data volume were used for b value estimation. The estimated b value using FGS method and corresponding standard error of each rock specimen are given in Table 1 and shown in Figure 3a. The temporal variations of the b values of granite, limestone, and red sandstone are shown in Figures 3b-3d. These were estimated using a running window of 4,916, 200, and 280 events and a running step of 2,458, 100, and 100 events, respectively. Figure 3a shows that the FGS method visually identifies the linear segment in the logarithmic incremental frequency-amplitude distribution, and the b values of the three rock specimens as listed in Table 1 are distinct. As these three kinds of rock specimens are bonded by mineral grains of different sizes and different types of discontinuities, their rupture scaling during deformation will also be different. Red sandstone is composed of fine-grained particles and seldom has large discontinuities, which accounts for a high proportion of small-scale ruptures. Granite has various large mineral grains and defects or voids. Limestone contains numerous of joints created during deposition, and these generate more large-scale ruptures. Therefore, the b value of red sandstone is the largest, followed by granite, and the b value of limestone is the smallest. Furthermore, Figures 3b-3d show that the changes in b value are seen consistently if we use different metrics of source size, that is, AE energy or amplitude. The appearance of a high-energy and high-amplitude signal is usually triggered by large-scale ruptures, which in turn leads to a decrease in the b value. The experimental results in this research show that the b value estimated by using the FGS method can well demonstrate the size distribution characteristic in rock AE tests.

Conclusions
1. The G-R relationship is an intrinsic characteristic of the frequency-size distribution from earthquake to laboratory rock ruptures, indicating that the number of small-scale ruptures is much greater than the number of large-scale ones. Nevertheless, deviation from G-R relationship is inevitable in practice. In fact, whether small-scale events cannot be fully recorded owing to limited equipment accuracy or small-scale ruptures are masked by larger ones that mainly occur in the avalanche destruction stage, the cause of the deviation of the size distribution from G-R law can be ultimately attributed to inadequate data acquisition. Furthermore, from our discussion of attenuation effect on the amplitude distribution, if the source amplitude distribution interval in Equations 7 and 16 is set as infinite (−∞, +∞), a finite attenuation actually cannot cause deviation on both ends of the amplitude distribution; however, once the size distribution is truncated at a certain interval, deviation will occur at both ends. This indicates that data truncation actually is the fundamental cause of attenuation modification of the size distribution. Inadequate data acquisition and data truncation will both cause deviation at both ends of the size distribution; therefore, the data catalog is critical, and the b value should be estimated in a finite interval. The proposed b value estimation procedure (FGS) in this study is designed to search that finite interval, and the application also shows its reliability. 2. Because the b value is a statistical value, its estimation result is significantly influenced by the data volume. Generally, a sufficiently large number of events can improve the robustness of the b value. However, when spatial or temporal variations in b value are investigated, the pursuit of high resolution of the variations in b values will reduce the number of events, which correspondingly leads to a decrease in reliability and robustness of the estimation. In other words, the more we want to increase the resolution of the variations in b values, the less precisely we are able to simultaneously estimate the b value because fewer events are used. Nevertheless, a reliable estimation is crucial in b value analysis. Therefore, the robustness of b value estimation should be guaranteed in priority even at the sacrifice of resolution of the spatial or temporal variation.

Data Availability Statement
The original AE data of dilation rupturing tests used in this work will be uploaded to a public repository called Digital Rocks Portal. The data archiving is underway now; we temporarily provide a copy of these AE data as supporting information named Data Sets S1-S3 for review purposes.