## Focal mechanism estimation

We utilized automatic picking and applied focal mechanism determination processes to continuous seismic records after checking the quality of the extracted 100 event data by comparing the manually picked and automated results29.

Similar to a previous study30, we adopted a smoothed, depth-dependent velocity structure for hypocentre determination. Only events within the area covered by the 2017 network were considered. The focal mechanisms of the selected data are determined using the HASH algorithm31. The completeness of the detection process reached an M of approximately − 0.5 for 3,641 events within the 2017 network29. High-quality focal mechanism solutions were also obtained. For each event, the optimal solution was derived as the tensor average from all possible focal mechanism tensors. We also generated 1000 focal mechanisms for an event based on the solution distribution using HASH. The datasets are utilized in the subsequent processing. The number of focal mechanism solutions satisfying the criteria was 2,828 for the 2017 observations (Fig. A1). We created 1000 datasets of focal mechanism solutions by random sampling from every 1000 focal mechanism solutions. This dataset is used to calculate the confidence range of the stress field.

## Stress field

The normalised deviatoric stress tensors were estimated for the spatial bins dividing the target area. The tensor was derived from the seismic moment tensor data of the aftershocks, utilizing a method18 based on the equivalent relation between the inelastic strain and the released seismic moment tensor32. To stabilize the solution, we employed the stress field estimation using data in the 2017 observation and seismic observation performed in 2000 as a previous study22. The focal mechanisms of 2,285 events for 2000 observations were obtained using the processing sequence described above (Fig. A1). Subsequently, we created a moment tensor dataset based on the focal mechanism and earthquake magnitude (Ma). The tensor shape and scalar moment (M0) were obtained from the estimated focal mechanism and transformed value from Ma through an empirical relationship used in a previous study24: (log (M0) = 1.15Ma + 10.548). Notably, stress tensor estimation was independent of the fault plane selection, as it utilized moment tensor data. The stress tensor estimated from the seismic data provides a normalised tensor shape that is characterised by the principal directions and ratio among the eigenvalues. The deviatoric stress tensors were estimated in spatial blocks, set at intervals of 3 km horizontally and at 2.5 km depth (Fig. A2). In addition, we performed the same procedure for the block distribution shifted by half a block to the north and east from the origin to smoothen the results. The principal directions and relative values of the eigenvectors that revealed the normalised deviatoric stress tensor were estimated from the events in a spatial bin. The parameters in each spatial bin within a confidence range were estimated (Fig. A3). This result aligns with those of previous studies22–24. However, the parameters with substantial accuracy were mainly limited to spatial bins around the aftershock area with relatively high seismic activity. Consequently, in the subsequent discussion, we focused solely on the data within the spatial bins exhibiting significant accuracy, as shown in Fig. A2. Finally, we analysed 2,155 (1,412 for Mc = -0.3) focal mechanism data points to investigate the strength dependency of the b-value.

A comparison of the number of focal mechanisms with hypocentres, including spatial bins, is shown in Fig. A4 as the frequency-magnitude distribution. The number of focal-mechanism data points was slightly lower than that of the hypocentre data. This difference can be attributed to a less favourable geometrical condition of the station distribution or a lower signal-to-noise ratio of the seismograms used to determine the focal mechanism. The estimated b-values exhibit a discrepancy of approximately 0.07 due to the larger difference in the small magnitude range. However, no specific focal mechanism was identified in this dataset. Therefore, the relative b-value changes among the parts of the Mohr diagram are discussed.

## Shear and normal stress estimation for stress tensor and choice of fault planes

The normalised shear and normal stresses were calculated from the fault geometry and stress tensor (\(\sigma\) ). According to previous research14, normal stress vector \({\sigma }_{N}\) is a product of traction vector **T** and unit fault normal vector **n [**\({\sigma }_{N}=(T\bullet n)n\)**]**, where the traction vector is \(\sigma\) applied to **n**; \(T=\sigma n\). The shear stress vector \(\tau\) is expressed as \(\tau =T-{\sigma }_{N}\). From the seismic data (focal mechanism or moment tensor), we can only determine the normalised deviatoric stress tensor \(\widehat{\sigma }\), characterised by the principal directions of the eigenvalues and the ratio between the eigenvalues. This means that the magnitude of \(\sigma\) is arbitrary; therefore, the relative relationship between the shear and normal stresses was obtained. In this case, the relationship provides \({\varvec{\sigma }}_{\text{N}}/\left|T\right| \text{a}\text{n}\text{d} \varvec{\tau }/\left|T\right|\) from \(\widehat{\sigma }\) and **n**. Therefore, we can define a Mohr circle with a unit diameter, with its centre point at the coordinates’ origin. For instance, the maximum shear stress under the stress tensor acting on a plane with a normal vector that rotates 45° from the maximum principal direction to the minimum principal direction of the stress tensor is plotted at (normal stress, shear stress) = (0, 1) on the normalised (unit) Mohr circle. In this study, we used coordinates to express the location on the Mohr circle by r (radial length from the centre, i.e., ranging from 0 to 1) and θ (clockwise angle from the horizontal).

To estimate the normal and shear stresses on the unit Mohr circle, a fault-normal vector is required. We have two candidates for the normal vector for each event from the focal mechanism solution owing to the existence of two nodal planes. The following scheme was used to select one of the two planes. Initially, we eliminated the fault plane from the dataset if it had a misfit angle (dθ), which is the discrepancy between the slip vector of the focal mechanism and the maximum shear stress direction expected from the stress field exceeding 30°. We then set two conditions that are used in the plane selection: (1) misfit angle (dθ) and (2) comparison lengths of segments (dp) to the line defined by Coulomb failure function (CFF) for a certain friction coefficient \(\left(\mu \right)\) from two points calculated by the normal vectors of the candidates. A schematic of this process is shown in Fig. A5. Shear strength in CFF is expressed as \({{\tau }}_{0}+\mu ({s}_{n}-p)\), where t0, *s**n*, and \(p\propto \text{d}p\) represent the cohesive strength, normal stress on a plane, and pore fluid pressure proportion to dp. A point touching the CFF line on the unit Mohr circle with angle θc can be expressed by \({{\theta }}_{\text{c}}={\text{tan}}^{-1}\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mu $}\right.\right)\) as shown in Fig. A5b and reported by Jaeger et al.14. Unfavourable fault planes for the stress field require a high value of p to generate rupture, which is more challenging to achieve naturally than in situations with a low *p* owing to difficulties sustaining a high *p*.

The dθ range for each plane was calculated from the stress fields within a 95% confidence interval in the spatial bin containing the event. Initially, we selected a plane with a smaller dθ range than the other plane (condition 1). However, if the dθ ranges overlapped, we chose a plane with a smaller dp (condition 2). This criterion is known as selection based on fault instability, as defined in a previous study33. To calculate dp, we set *µ* = 0.6 based on the typical internal friction coefficient27. However, despite assuming *µ* = 0.2 or 1.0, the result did not change (Fig. A6).

After completing the above process for each event, we compared the selected fault plane for an estimated stress tensor with that for another one. As described in “stress field”, spatial block distribution shifted half the block size horizontally, meaning two candidates of the stress field exist for most events. We used a fault plane and stress field combination to provide a relatively small misfit dθ.

## b-value estimation

The b-value, which represents the power-law decay of the frequency-magnitude distribution, has been estimated for various regions and exhibits spatial and temporal variations. We followed the standard method to estimate the b-value34–36. We set a consistent magnitude range (minimum and maximum magnitudes) for each part of the Mohr circle, as comparing the b-values across different parts of the circle is crucial for this study. We used the magnitude range from − 0.3 to 2, with the lower bound determined using ZMAP software26. The upper bound is the set position at which the segment of the frequency-magnitude relation is slightly deviates from the linear relation between log N and M (i.e. off from the Gutenberg-Richter relation). Additionally, the number of earthquakes above the upper bound was too small for the analysis of the b-value variation within the Mohr circle.

We estimated the b-values for the parts of the Mohr circle that contained over 50 events, excluding events close to the co-seismic fault of the 2000 W-Tottori EQ. The heterogeneity of the stress field may lead to improper estimation of the fault direction relative to the stress field. Based on the co-seismic fault model21, we calculated the distance to the fault planes and selected events that were located farther than a certain distance (R). In this study, we set R = 0.5 km. The results for different distance settings are shown in Fig. A7a). Increasing R significantly decreased the number of usable events due to the concentration of hypocentres around the co-seismic fault planes (Fig. A7b). The number of excluded events with larger misfits, as described in the previous section, decreased with increasing R, indicating the influence of the nearby faults on events. We selected R = 0.5 km, which led to the utilization of 65% of events from the initial set, while 8% were excluded..

The unit Mohr circle is divided into parts expressed by range (r minimum, r maximum; θ° minimum, θ° maximum) in the coordinate shown in the inset of Fig. 2(c): those ranges are (0.8, 1: 10, 30), (0.8, 1; 20, 40), (0.8, 1; 40, 60), (0.8, 1; 60, 80), (0.8, 1; 70, 90), (0.8, 1; 80, 100), (0.8, 1; 90, 110), (0.5, 0.8; 22.5, 67.5), (0.5, 0.8; 45, 90), (0.5, 0.8; 67.7, 112.5), (0, 0.8; 0, 90), and (0, 0.8; 90, 180).