## 3 Estimation of \(b\)-Value and Fractal Dimension  
Following the definition of the body-wave earthquake magnitude, the magnitude of an AE event was determined according to the log of the maximum amplitude (\(A_{max}\)) of the AE signal (roughly say, the vibration velocity of the elastic wave [21]):  
\[\mathrm{M}=\log(A_{max}/A_{0})+M_{0} \tag{1}\]  
where \(A_{0}\) was the threshold for detection. Such a definition results in magnitudes consistent with earthquake magnitudes. Since wave attenuation is a function of ray path, some events having insufficient energy could be recorded by only the detector nearest the AE hypocenter. As a result, the number of events counted generally differed between two detectors. Due to the unknown calibration parameter of the PZT sensors, relative magnitudes were normally used. However, \(M_{0}=-11\) was used to obtain a rough estimation of the absolute AE magnitude. As mentioned in the introduction, it is commonly accepted that AE activity in rocks follows the Gutenberg and Richter (GR) relation for frequency-magnitude distribution [22]:  
\[\log_{10}N=a-bm \tag{2}\]  
where \(N\) is the number of events of magnitude \(m\) or greater, and \(a\) and \(b\) are constants. Data recorded by the peak detectors were used to calculate seismic \(b\)-values. In order to investigate the temporal evolution of the \(b\)-value, \(b\)-values were calculated using the maximum likelihood method for a running window of a constant number of AE events (\(n\)) with a running step of \(n\)/4 events (75% overlap). Previous studies have highlighted that the number of AE events influences \(b\)-value analysis [23]. Thus, a sensitivity analysis was performed for a different number of windows. Figure 2 shows \(b\)-values vs. time for \(n=500\), 1000, and 2000. It is clearly demonstrated that a larger \(n\) results in a smoother \(b\)-value curve, whereas a smaller \(n\) results in a higher resolution of \(b\)-value evolution. As shown in Figure 2, the results for \(n=2000\) were sufficiently smooth while maintaining the most basic characteristics of \(b\)-value evolution, and thus, they were used in the present study. Waveform data were used for hypocenter determination. First, the method used for the automatic measurement of the arrival times of seismic signals [24] was used in the pickup of P-wave 1 st motion for AE. Acoustic emission hypocenters were determined automatically using the first arrival time data of P-waves. Time-varied mean anisotropic velocities averaged from measured P-wave velocities (Vp) along multiple paths were used in hypocenter determination. For an event, clear first motions from at least four sensors were required for estimating four unknown parameters (three for source location and one for the original time). Due to heterogeneities in velocity structure and uncertainty in arrival times, more sensors generally resulted in better precision. In the present study, AE hypocenters determined by at least eight P-wave arrival times were used. The estimated location error was less than 2 mm for most hypocenters. In order to examine the fractal behaviors of AE hypocenters, the correlation integral function was used for the estimation of multi-fractal dimensions [25]:  
\[C_{q}(r)=\frac{1}{n}\left[\sum_{j=1}^{n}\left(\frac{n_{j}(R\leq r)}{n-1} \right)^{q-1}\right]^{1/(q-1)}(q=-\infty,\ldots,-1,0,1,\ldots,\infty) \tag{3}\]  
where \(n_{j}\) (\(R<r\)) is the number of hypocenter pairs separated by a distance \(R\) less than or equal to \(r\), \(q\) is an integer, and \(n\) is the total number of AE events analyzed. The AE hypocenter distribution generally exhibits a power law for any \(q\), \(C_{q}(r)\propto r^{D_{q}}\), and \(D_{q}\) defines the fractal dimension. The present study focused on \(D_{2}\), which coincides with the correlation dimension. Spatial correlation length (SCL) was also used to characterize the spatial distribution of AE hypocenters, which does not require a power-law distribution.