## 2 \(b\)-Value Calculation Methods  
In our estimation process of the \(b\)-value, we employed the widely recognized Gutenberg-Richter relationship [10] in statistical seismology, which is expressed as \(logN=a-bM\), where \(N\) denotes the cumulative number of earthquakes with magnitude equal to or greater than \(M\), and \(a\) and \(b\) represent constants. Here, the \(a\)-value indicates the level of seismic activity within the region, while the \(b\)-value, representing the slope of the magnitude-frequency distribution, is typically associated with the region's geological structure and prevailing stress state. Given the linear nature of this relationship, the Least Squares Method (LSM) has been traditionally utilized to estimate \(a\)- and \(b\)-values. However, LSM tends to assign equal weight across all magnitude bins, which can lead to disproportionate influences from smaller sample sizes in higher-magnitude bins on the estimated \(b\)-value. Such disproportionate effects can result in unstable \(b\)-value estimates, which might lead to unreliable predictions in seismic risk assessments. To address these limitations, we recommend and have utilized the maximum likelihood estimation (MLE) method in this study. The MLE method treats the magnitude \(M\) as a continuous random variable and assumes that it follows a power-law distribution, as outlined by [PERSON] [11] and [PERSON] [12]. They provided foundational formulas for the maximum likelihood estimation of the \(b\)-value, enhancing the precision of \(b\)-value calculations by effectively accounting for different sample sizes across magnitude bins. This approach significantly reduces the biases associated with LSM, thereby providing more stable and reliable estimates of \(b\)-values. Based on this, [PERSON] and [PERSON] separately provided expressions for the maximum likelihood estimation of the \(b\)-value:  
\[\hat{b}=\frac{\log(e)}{\overline{M}-M_{C}} \tag{1}\]  
\(\overline{M}\) represents the average magnitude, and \(M_{C}\) is the magnitude of completeness. [PERSON] [13] further refined this method, improving the accuracy of the \(b\)-value estimations under varied seismic conditions:\[\hat{b}=\frac{\log(e)}{\overline{M}-(M_{C}-\Delta M/2)} \tag{2}\]  
The MLE method calculates the \(b\)-value based on the probability density values across different magnitude segments, avoiding the disproportionate weight issues characteristic of the LSM. This enhancement in methodological approach leads to more accurate and stable \(b\)-value estimates, which are crucial for understanding the underlying stress states in seismic regions. Recent research emphasizes the critical role of accurate \(b\)-value estimation in understanding regional stress levels and seismic risk. Studies have demonstrated that lower \(b\)-values are often associated with high seismic coupling and geological asperities, particularly in regions like the Chilean megathrust, highlighting areas of significant stress accumulation that are prone to large earthquakes [14]. For instance, [PERSON] et al. (2015) observed rapid changes in \(b\)-values before and after the Tohoku-oki earthquake, providing insights into the randomness of megathrust earthquakes and the dynamics of stress recovery [15]. This study further illustrates how sudden shifts in \(b\)-values can signify underlying stress adjustments that might not be apparent from surface observations alone. Other investigations have linked rapid changes in reservoir water levels to lower \(b\)-values, suggesting that such environmental modifications can induce seismic activity through increased pore pressure and resultant stress [16]. Moreover, a trend of decreasing \(b\)-values preceding major earthquakes has been observed, indicating rising shear stress within fault zones. These findings collectively validate the use of \(b\)-values as reliable indicators of stress variations, serving as essential tools for seismic risk assessment and the development of mitigation strategies [17].