Stress relaxation arrested the mainshock rupture of the 2016 Central Tottori earthquake

After a large earthquake, many small earthquakes, called aftershocks, ensue. Additional large earthquakes typically do not occur, despite the fact that the large static stress near the edges of the fault is expected to trigger further large earthquakes at these locations. Here we analyse ~10,000 highly accurate focal mechanism solutions of aftershocks of the 2016 Mw 6.2 Central Tottori earthquake in Japan. We determine the location of the horizontal edges of the mainshock fault relative to the aftershock hypocentres, with an accuracy of approximately 200 m. We find that aftershocks rarely occur near the horizontal edges and extensions of the fault. We propose that the mainshock rupture was arrested within areas characterised by substantial stress relaxation prior to the main earthquake. This stress relaxation along fault edges could explain why mainshocks are rarely followed by further large earthquakes. Rupture during the 2016 Central Tottori earthquake, Japan stopped within regions characterised by stress relaxation before the earthquake, according to analyses of aftershock focal mechanisms. This could explain why large earthquakes typically do not follow each other.

F ollowing a large earthquake, shear stress increases along the edges and extensions of the fault where the slip of the earthquake did not occur 1 . This stress is added to the stress that accumulated before the earthquake; thus, if the fault and its surroundings were in a state of uniform stress before the earthquake, it is thought that the fault edges and extensions would subsequently be in a very unstable state, allowing earthquakes to be easily generated 2 . This is demonstrated in Fig. 1a, in which the stress profile along the fault strike is shown for a vertical strikeslip fault with a simplified slip model (see "Methods"). It is assumed that the same amount of shear stress as the averaged stress drop on the mainshock fault is applied uniformly on the fault, inside and outside the mainshock rupture area during an interseismic period. This stress level is thought to indicate the apparent strength of the fault, since the mainshock occurs when the shear stress increases to this level. If the strength is assumed to be uniform on the fault and its extension, it is clear that the stress after the mainshock becomes greater than the strength near the edges of its rupture area and beyond. Thus, large earthquakes will be generated there. However, in intraplate regions, large earthquakes are rarely followed by other large earthquakes in adjacent areas; typically, only smaller earthquakes occur and, then, decrease in frequency roughly with the reciprocal of time after the mainshock. However, occasionally seismic sequences are characterised by the occurrence of multiple large earthquakes [3][4][5] , and forecasting such phenomena remains a major issue in seismology 3 . At present, the cause of seismic activity following large earthquakes is controversial and includes static [6][7][8] and dynamic triggers [9][10][11] as well as the effect of pore fluid pressure 12 , and no clear physical model has been established to date. Therefore, aftershocks are forecasted statistically 13 , and it has been reported that aftershock activity is sometimes unusual before other large earthquakes 3 .
The question of whether or not a subsequent major earthquake will occur along the edges and extensions of a mainshock fault is compounded by the questions of why rupture propagation of a large earthquake stops at its edge and how the size of the earthquake is determined 14,15 . When rupture propagation stops at an ordinary high-strength barrier, it may be possible for a large event to occur 16 . On the contrary, in the case of low strength and/or low-stress barriers, such as volcanos [16][17][18][19] , this may not be common. Ke et al. 19 emphasised the importance of the low-stress barrier called the earthquake arrest zone, using laboratory experiments and dynamic rupture simulations. Endpoints of earthquake faults, large bends 20 , offsets 21 , or damage zones, such as horsetail splays or branch faults 22 , are recognised from the surface traces of active faults. However, a systematic study of strike-slip fault surface traces showed that at one endpoint of nearly half of all faults analysed, the reason for the cessation of rupture is unknown 23 . Moreover, the stress  Fig. 1e. c Case of stress concentration by the aseismic slip below the seismic belt 38 . The aseismic slip (B) occurred in the ductile fault zone below the seismic belt shown by pink in Fig. 1e. d Slip profiles of the coseismic slip and the aseismic slip around the mainshock fault for the case (b). e Geometries of the mainshock fault and the surrounding aseismic fault, and the ductile fault below the seismic belt. The slip distribution of the mainshock fault is shown by colour contours.
state at the edge of earthquake faults has yet to be estimated from observations.

Results and discussion
Estimation of the true mainshock fault plane and its horizontal edges. In this study, we estimated the stress state at the horizontal fault edges of a simple, moderately sized, strike-slip intraplate earthquake by precisely estimating the locations of the fault plane of the mainshock, including their horizontal edges, with an accuracy of~200 m, from approximately 10,000 accurate focal mechanism solutions (see "Methods"). These solutions were determined from the data of both 44 permanent stations and 69 temporary seismic stations installed immediately after the occurrence of the 2016 Central Tottori earthquake (Mw6.2) 24 that occurred in the seismic zone near the coast of the Sea of Japan in the San'in district 25 . The distributions of the earthquake's aftershocks and seismic stations, and examples of focal mechanism solutions are shown in Fig. 2 and Supplementary Fig. 1, respectively. A foreshock of magnitude 4.2 occurred about 2 h before the mainshock at almost the same location, as the mainshock and small earthquakes continued 26 .
In recent years, it has become clear that many aftershocks occur outside the fault planes of mainshocks 27,28 . The changes in stress around the edges of faults caused by the mainshocks are tensional on one side and compressional on the other. Furthermore, the direction of the maximum compressive stress generated by the mainshocks is rotated 90°, with the fault as a boundary. It has previously been reported that aftershock focal mechanism solutions are in consistent with the 90°rotation 29 . For the 2016 Central Tottori earthquake, it was confirmed that nearly all aftershocks occurred outside the fault plane, since the orientations of the tension (T) axis of the optimal focal mechanism solutions for the Central Tottori earthquake vary by 20-40°near the fault around the southern edge, as shown in Fig. 3a. Similar off-fault aftershocks were detected for the 2000 Western Tottori Prefecture earthquake, and it is estimated from precise aftershock distributions that they occurred outside the damage zone 28 . It is also shown in Fig. 3a that with respect to the pressure (P) axes, those with small plunges show similar differences. These are results from a depth of 5 ± 0.5 km, but similar results are seen at other depths, as shown in Supplementary Fig. 2. The locations of the mainshock fault planes and their horizontal edges were estimated from focal mechanism solutions and GNSS data 30 (see "Methods"). In the northern part, the fault planes are complicated and were estimated to be separated into two, the western fault (WF) and the eastern fault (EF), as shown in Fig. 3. Thus, the results obtained in the southern part are mainly discussed in this paper. However, similar shifts of T and P axes were observed around the northern edge of the WF.
It is worth noting that normal faulting aftershocks were concentrated on the eastern side of the southern fault edge, namely in the tensional region (Fig. 3c). The regions around the fault edges where the normal faulting aftershocks were concentrated matched well with the distribution of changes in the Coulomb failure stress (ΔCFS) for a normal fault calculated from the slip distribution of the mainshock, as shown in Fig. 3d.
Since normal faulting aftershocks are considered to be more easily triggered where the ΔCFS is large 31 , we assumed that horizontal fault edges were located at the positions where the spatial average of the ΔCFS value for the normal faulting aftershocks was at maximum, and estimated the locations of the edges, as shown in Supplementary Figs. 3 and 4 (see "Methods").
A lot of aftershocks occurred around the fault edges, due to the stress concentration. In contrast, aftershocks also occurred in the central part, where the ΔCFS values were small for both strike-slip and normal faulting aftershocks. It is difficult to explain this by using a simple fault model as that employed in this study. The source model inferred from strong-motion waveforms has two large slip regions in the northern and southern parts 32 , and this model can explain those aftershocks. The purpose of this study is to investigate the properties of the horizontal fault edges, and the GNSS data used in this study do not have a resolution high enough for solving such a small region 30 , therefore, this problem requires further research.
Stress state near the fault edge. Figure 4 shows the distribution of the P axis plunges of the aftershocks near the estimated southern fault edge at depths of 5 ± 0.5 km, together with ΔCFS for normal (Fig. 4a) and strike-slip ( Fig. 4b) faulting aftershocks. These are results from a depth of 5 km, while those from all the depths are shown in Supplementary Figs. 5 and 6. The normal faulting aftershocks were notably concentrated in the region where the ΔCFS for the normal fault was larger. Meanwhile, in the region where the ΔCFS for the strike-slip fault with the same focal mechanism solution as the mainshock was larger, there were almost no strike-slip faulting aftershocks, which are denoted by blue colours. The strike-slip faulting aftershocks also appear to avoid the large ΔCFS region. At other depths in the southern part, a concentration of normal faulting aftershocks was observed, in particular at depths of 4 and 6 km, and strike-slip faulting aftershocks are absent, except at~8 km ( Supplementary Fig. 5).
The ΔCFS values in this study were calculated for stress changes due to the mainshock. It is known that before the mainshock, the stress field in this district was prone to strike-slip faulting earthquakes 25,29 (Fig. 2). If this stress field was spatially uniform, the total shear stress after the mainshock (expressed as the sum of the stress field before the mainshock and the change in stress due to the mainshock) would be expected to greatly increase along the fault edge. Many historical, large earthquakes in this area 25 have suggested that the stress state is critical (i.e., close to failure). If the stress concentration along the fault edge is further added to this critical uniform stress field, it is expected to be extremely prone to earthquakes; however, at least during the four years after the mainshock, no larger earthquakes have occurred.
What happened near the fault edges? There are four possible hypotheses (1-4) for the observations concerning the fault edges: (1) The stress field before the mainshock was homogeneous, and the strength along the edge and extension of the fault may have been very large, and in extreme cases, there were no preexisting planes of weakness in the extension; (2) the stress field before the earthquake may not have been homogeneous, such that the concentration of stress only occurred in the area of the mainshock slip; (3) before the earthquake, stress relaxation occurred along the edges and extensions; or (4) after the earthquake, inelastic deformation may have occurred along the edges and the extension.
Hypothesis (1) corresponds to an ordinary high-strength barrier, as stated above. This is hypothesised, since large earthquakes are inevitably generated at the edge if the strength of the fault is uniform under the assumed homogeneous stress field before the mainshock, as shown in Fig. 1a. However, it is known that in the case of this barrier, aftershocks are concentrated around the fault edges 16,33 . Furthermore, hypothesis (1) is unlikely because, as shown at the depth of 5 km in Fig. 4b, several normal faulting aftershocks occurred, even in the region where there were almost no strike-slip faulting aftershocks despite of the large ΔCFS for the strike-slip fault. Meanwhile, hypothesis (4) is inconsistent with the results of a geodetic study, which showed that post-seismic slip was not detected along the fault edges and extensions, but instead in the co-seismic slip area and in the near-surface shallower region 30 .
Aseismic slip along the edges and extensions is a strong candidate that realises hypothesis (2), as shown in Fig. 1b This is an extreme case, in which aseismic slip (A) mainly occurred on the fault shown by grey around the mainshock fault and the magnitude of the slip is the uniform block slip minus the coseismic slip, as shown in Fig. 1d. The stress due to the aseismic slip is shown as a dashed line in Fig. 1b, which is the reverse of the stress change due to the coseismic slip 34 . This would have occurred near the region of strong coupling at the plate boundary 35 . After the earthquake, the stress before the earthquake is cancelled by the coseismic change, then the stress state becomes flat and there is no stress concentration. In this case, no large earthquakes occur at the edge and in the extension. However, this may not have been the case, for the same reason why we rejected hypothesis (4), since an aseismic slip area is generally characterised by post-seismic slip 36 . Furthermore, postseismic slip is usually accompanied by the migration of aftershocks 37 , which was not observed for the Central Tottori earthquake 30 .
As the above case is an extreme one, we examined another possible candidate for stress accumulation only near the earthquake fault: the aseismic slip in the ductile fault zone below the seismic belt in the San'in district 25,38 . Nishimura and Takada (2017) explained GNSS velocity data by the aseismic slip below a depth of 16 km shown in Fig. 1e. The magnitude of the aseismic slip is calculated as the magnitude of the accumulated shear stress averaged on the earthquake fault is equal to that of the stress change by the mainshock. The stress states before and after the mainshock are shown in Fig. 1c as the dashed and dotted lines, respectively. Although unlike the case of the uniform stress shown in Fig. 1a, the region where the stress after the mainshock is greater than the strength is limited to the vicinity of the edges, the region is approximately a half of the length of the mainshock fault. In this case, it is unlikely that an earthquake as large as that expected in the case of the uniform stress will occur at the edge and in the extension. However an earthquake of a considerable magnitude, the fault of which is longer than half the mainshock fault, could occur there, and it is thought that the observed results cannot be explained by this model.
Hypothesis (3) corresponds to another barrier 17 , as stated above, and is schematically visualised in Fig. 5, in which the mainshock slip is confined by the regions of stress relaxation, by the near-surface layer, and the deeper ductile fault zone. The mainshock fault should have slightly penetrated into these regions, but presently, we have no clear evidence for this. In this case, a large earthquake is unlikely to occur in this extension because the shear stress decreased due to stress relaxation. The reason why post-seismic deformation was not detected along the fault edges and extensions may be a large relaxation time due to low temperatures. Additionally, a large earthquake might not occur in the extension of the fault but in a conjugate direction 5 , for which the evaluation of the stress state requires further study.
We finally discuss the northern edge. There were few aftershocks along the EF in the northern part, and very few aftershocks occurred along the edge and in the northern extension, as shown in Fig. 3 and Supplementary Fig. 2. This is in harmony with the results in the southern part. Concerning the WF, the northern edge was not well determined and estimated to be at Y = 6.1−6.9 km, based on focal mechanism data (see "Methods", Supplementary Fig. 12). If the edge was at Y = 6.9 km, in one extreme case, this position almost coincides with the limit of the aftershock distribution, as shown in Fig. 3 and Supplementary Fig. 2. As shown in Fig. 2, a quaternary volcano is located near this northern edge 39 . In the case of the other extreme, Y = 6.1 km, as shown at the depths of 5-6 km ( Supplementary Fig. 6b), there were also few earthquakes in the region where the ΔCFS for the strike-slip faulting aftershocks were larger. However, this tendency is less remarkable than that around the southern edge, possibly because another minor fault other than the WF and EF may exist due to branching or horsetail structures 22 .
Since we were able to precisely determine the positions of an earthquake fault and its horizontal edges, by utilising a large number of well-determined focal mechanism solutions obtained from the rapidly deployed high-density seismic network 24 , we believe that our results are reliable. However, only one moderate earthquake was studied and whether or not the observed trends are common to other earthquakes is an important issue requiring further study. A similar pattern of the P axis azimuth was recognised around the southern edge of the 2000 Western Tottori earthquake (Mw6.6), and it was estimated that the deviatoric stress in the 2 × 2 km 2 area around the edge was as small as 5 MPa at a depth of 7 km 29 . Moreover, according to the results of seismic velocity tomography in the San'in district, including aftershock regions of moderate to large earthquakes, low-velocity anomalies were systematically estimated to be located near both ends of the aftershock regions 40,41 . These results suggest that stress relaxation at the fault edge may be a common phenomenon for large intraplate earthquakes that are not followed by larger earthquakes.

Methods
Aftershock locations. We only used manually picked data from stations within 50 km of the epicentre. We located hypocentres using the hypomh_ps programme 25 , assuming a one-dimensional velocity structure, estimated in the aftershock area of the 2000 Western Tottori Prefecture earthquake 42 , which was adjacent to the current study area. In the hypomh_ps programme, the origin time was eliminated from the location problem. Thus, we could estimate reliable hypocentres without a trade-off between the origin time and hypocentre 43,44 . We first determined hypocentres without station corrections. Then, we calculated the averages of the observed minus calculated (O − C) values for the P-and S-waves at each station for earthquakes that had ≥ 30 manually picked P-wave arrival time data points and used them as station correction information. While calculating the averages of the O − C values, values > 0.5 s for P-waves and > 1.5 s for S-waves were omitted. Most errors in the hypocentre calculation were 30-40 m in the horizontal direction and approximately twice that in the vertical direction. For more details, we refer to Iio et al. 24 . We confirmed that the estimated aftershock distribution was very similar to that estimated with tomoDD 45 ; however, the distribution was remarkably different from that estimated with the GrowCrust algorithm 46,47 , possibly because it allows iterative changes in origin times.
In this study, we analysed aftershocks that occurred from 06:33 Japan Standard Time (JST) on October 22 to December 15 in 2016. Although aftershocks that occurred from 14:07 JST on October 21 (when the mainshock occurred) to 06:33 JST on October 22 (when the temporal network began operating) were not included in our dataset, we confirmed from the unified catalogue of the Japan Meteorological Agency (JMA) that the hypocentral distribution during this period was very similar to that in the subsequent period.
Focal mechanism solutions. Focal mechanisms for aftershocks were determined using a grid search algorithm described in detail by Iio et al. 48 . In this study, we only used focal mechanisms with P-wave polarities of 15 or more, 10 or fewer faultplane solutions, and an uncertainty of < 60°(estimated from the Kagan angle 49 between the two most different solutions). The number of focal mechanisms that remained after the above screening process was 10,075, and the average uncertainty of the focal mechanisms was 13.2°. In this study, we selected one optimal faultplane solution from multiple solutions for each event and analysed the P and T axes of the solution. If all the best solutions for an event had no inconsistent polarities, we selected the solution with the nodal planes that were the farthest from the polarity data on the focal sphere; however, if the best focal mechanisms had inconsistent polarities, we selected a solution with a nodal plane closest to the inconsistent polarities 50 . Examples of focal mechanism solutions are shown in Supplementary Fig. 1. Those for aftershocks around the southern edge in Fig. 3 are shown separately for the eastern and western sides of the fault. Magnitudes ranged from −0.6 to 2.4. Each focal mechanism solution was well-determined, and many normal faulting aftershocks occurred on the eastern side.
Estimation of the true mainshock fault plane-T axis azimuth profiles. It appears that the azimuths of the P and T axes were shifted near the mainshock fault plane at all depths ( Fig. 3 and Supplementary Fig. 2). To confirm that these data reflected the stress state in the aftershock area, we performed stress inversion analysis according to the ordinary method, details of which are presented in Iio et al. 48 . We assumed that the stress state was uniform in each sub-region, which was 1 km cuboid, and the length of the sub-region in the X direction was reduced to include focal mechanisms less than 50. Each sub-region overlapped with the adjacent sub-region by half its length. Further, each sub-region was set to exclude the WF, at which clear T axis shifts were observed. As shown in Supplementary  Fig. 7, we found that the distributions of the azimuths of the P and T axes were in accordance with those of σ 1 and σ 3 , respectively; thus, it is reasonable to explain them based on the changes in stress caused by the mainshock.
Changes in the azimuths of the T axis along the fault-normal direction were investigated for every cuboid region, including the 1 × 1 km 2 subfault over the entire fault, and the results at a depth of 5 km are shown in Supplementary Fig. 8 (all results are shown in Movie S1). We only used the data with T axis plunges < 30°, and 140°< T axis azimuths < 260°, from −6 to −3 km in the X-direction. On the boundary of the fault, step-wise changes were observed for many of the regions, with changes showing opposite patterns in the north and south that were consistent with the changes in stress due to the mainshock slip. Assuming that the fault was located at the positions of these steps, we approximated the T axis azimuth data with a step-wise function (tanh(aX)) and attempted to detect the three-dimensional positions of the fault using a nonlinear least-squares method. The value of parameter a of the stepwise function was set as 20, considering errors in hypocentre locations. We weighted the data such that it was inversely proportional to the error of focal mechanism solutions measured by the Kagan angle 49 , as shown in Supplementary Fig. 8. The minimum error was set as 4°. Initial values were set by a grid search. Fault positions were estimated at grids at the centre of each subfault and the grids were distributed as each subfault overlapped the adjacent one by half its length. When data were not fitted well, results of the initial grid search are shown, but 95% confidence intervals (CIs) of estimated fault positions and T axis azimuths are not shown. For many grids, the estimated fault positions are different from ordinary estimations that were obtained by minimising the sum of the squared distance between the fault positions and aftershock hypocentres 28 , beyond the 95% CIs, suggesting that it may not be possible to accurately estimate the mainshock fault plane only from aftershock hypocentre data.
Estimation of the true mainshock fault plane-single fault model. We only used data estimated at the grids where the 95% CIs of T axis azimuths on both sides of the fault did not overlap, the AIC was significant when compared to the model without the step, the amount of data on one side of the estimated fault was ≥ 5, Z > −11 km, and −3 km < Y < 7 km; the last condition was used to exclude data far from the expected mainshock fault plane. From these data, we estimated the mainshock fault plane by local linear approximation using the 'fit' function in MATLAB v.9.3.0 (MathWorks, USA) with the "lowess" (locally weighted scatter plot smooth) option. We used fault position data weighted inversely proportional to the 95% CIs and set a parameter for the extent of data used by the local regression to be 0.1, which corresponded to a cut-off length of~2 km, that is, twice the subfault length. If smaller values for this parameter had been used, a bumpy fault plane would have been estimated. The estimated fault plane is shown in Supplementary Fig. 9, along with the data used. The data are also shown in Supplementary Fig. 3b with the distributions of T axes. In this approximation, finer grids were estimated at the intervals of~200 m in the fault strike direction and 500 m in the dip direction, for subsequent calculations of stress. However, it appears that several fault position data are distributed on the eastern side, apart from the estimated fault, in the region of Y > 2 km. Furthermore, we found that slip on the estimated fault cannot fully explain the coseismic displacements deduced from the GNSS measurements shown by red arrows in Supplementary Fig. 10 30 .
Estimation of the true mainshock fault plane-two-fault model. Then, we estimated two fault planes from the same dataset obtained above. First, we estimated the WF by using the data where X < −4.2 km if Y > 3 km. If Y ≤ 3 km, all the data were used. Next, we set the weights of the fault positions of Y < 2 km estimated in the first step, to 1, namely the fault positions of Y < 2 km were almost fixed, and estimated the EF by using the data where X ≥ −4.2 km if Y > 3 km. The estimated WF and EF are shown in Supplementary Fig. 11. It was found that the data from the eastern side, apart from the single fault, were well fitted by EF.
In these steps, the edges of the fault planes were not determined, since extrapolation was also performed outside the aftershock region. Supplementary Figs. 9 and 11 were drawn considering the results for the fault edges described below.
Estimation of locations of the southern fault edge. Here, we first estimated the southern fault edge. Tensional stress is generally concentrated on one side of the horizontal edges of strike-slip faults. As the regional stress state across the islands of Japan is essentially compressional 51 , normal faulting earthquakes are the most likely to occur in regions where tensional stress is concentrated. Thus, it can be expected that the horizontal edges of the mainshock fault of the Central Tottori earthquake would be located where the ΔCFS for normal faulting aftershocks was maximised. In the calculation of ΔCFS, μ was set to 0.4, following previous aftershock studies 6 . We then assumed several fault models with fault edges in different locations, and calculated ΔCFS for normal faulting aftershocks with a Paxis plunge > 60°. For simplicity, horizontal edges were not allowed to vary with depth. The upper and lower edges were set to depths of 1.6 and 11.8 km, respectively, following the aftershock distribution. The distributions of slip magnitudes were assumed to be tapered linearly by widths of 4 km from the horizontal, 3.5 km from the upper, and 5 km from the lower edges of the fault, following the fault model with varying slip magnitudes estimated from InSAR and GNSS data 30 . Here, we did not refer to the fault slip models estimated from strong-motion data 32,46 , because these models cannot explain the coseismic displacements 30 . We confirmed that the minor changes in the slip distribution, such as different taper lengths, rarely affect the results. The absolute amount of slip was normalised so that the seismic moment matched that determined by the fault model 30 .
We calculated the contribution from all of the triangular subfaults obtained by the Delaunay separation of the grid data of the curved fault plane estimated above using the Meade programme 52 , to solve the problem of singularities 53 by shifting points to calculate. The distributions of ΔCFS for normal faulting aftershocks are shown in Supplementary Fig. 3 for different positions of the southern fault edge. The moment ratio of the western and eastern faults was assumed to be 0.4:0.6, and both northern edges were set at Y = 6.1 km, which did not affect the following result. Aftershocks within 100 m from the fault plane were excluded, because errors in hypocentres for those events could cause substantial adverse effects. The ΔCFS values were computed for both the nodal planes and larger values were adopted. We averaged ΔCFS within a cuboid area with a height of 13 km and a horizontal cross-section of 2 km 2 , centred on assumed fault edges. The relationship between the mean ΔCFS and the positions of horizontal fault edges is shown in Supplementary Fig. 12a. The position of the southern edge was well-estimated as the position that maximises the mean ΔCFS value (Y = −2.6 km). This is within 1σ CIs of those of the uniform slip model 54 .
Estimation of locations of the northern fault edges. To estimate the locations of the northern fault edges, we utilised the coseismic displacements deduced from the GNSS measurements 30 (which are shown by red arrows in Supplementary Fig. 10), because few aftershocks occurred around the northern edge of the EF, as seen in Fig. 3 and Supplementary Fig. 2, and several GNSS stations are located near the northern edges ( Supplementary Fig. 10). We assumed several variations for the locations of both the fault edges and the moment ratio of the WF and EF and calculated the differences between the observed and calculated coseismic displacement vectors. For simplicity, the locations of the horizontal edges were not allowed to vary with depth. As in the previous section, the distributions of slip magnitudes were assumed to be tapered linearly by widths of 4 km from the horizontal, 3.5 km from the upper, and 5 km from the lower edges, for both the faults.
The results are shown in Supplementary Fig. 13. We assumed the locations of the northern fault edge of the WF from 5.7 to 7.3 km in 0.2 km steps. We calculated the RMS misfits (mm) between the observed and calculated vectors for locations of the northern fault edge of the EF from 5.7 to 7.3 km in steps of 0.2 km and the moment ratios of the WF from 0 to 1 in steps of 0.2. For many cases of the WF, the minimum value was obtained for the location of the northern fault edge of the EF of 6.1 km and the moment ratios of the WF of 0.4. A ratio of 0 or 1 meant only the EF or WF, respectively, and misfits were exceptionally large in these cases. Using these values, changes in the misfits are shown for different locations of the northern fault edge of the WF in Supplementary Fig. 12c. A broad minimum with a gentle peak was observed around Y = 6.5 km.
By the same method used in the previous section, we estimated the location of the northern fault edge of the WF from focal mechanism data. However, as shown in Supplementary Fig. S12b, a broad peak can be observed from Y = 6.1 to 6.9 km, and the northern edge could not be clearly determined using the focal mechanism data or the GNSS data.
The estimated faults are compared with P and T axes, and stress states in Supplementary Figs. 2 and 7, respectively. The extent of the estimated northern edges of the WF is indicated by pink and blue dashed lines from Y = 6.1 to 6.9 km in Supplementary Figs. 2 and 7, respectively. In the case of Y = 6.9 km, it can be observed that the northern edge is located at the northern end of the aftershock distribution, and that very few aftershocks occur north of the edge in Supplementary Fig. 2. It appears that in this case, the WF almost penetrates the aftershock distribution. However, in Fig. 3 and Supplementary Fig. 6, we used the northern edge of Y = 6.1 km and calculated the stress field as the opposite extreme case. This is because, in the case of the northern edge of Y = 6.9 km, it is not effective to show ΔCFS around the edge, due to the absence of aftershocks north of the edge. The estimated fault models with the distributions of slips are shown in Fig. 4. In the southern part, the WF and EF are unified because of the same geometries.
Calculation of shear stress on a simple planar vertical strike-slip fault and in the extension. To evaluate the stress concentration near the horizontal edges of a mainshock fault, we assumed a simpler fault slip model and calculated the stress state after the mainshock in Fig. 1. We assumed that the fault plane was a simple planar vertical strike-slip fault with the same horizontal and lower edges as the EF. For simplicity, the upper edge was set at the surface and the slip distribution was assumed to be uniform in the vertical direction. This implicitly assumes that postseismic deformation occurs near the surface as was partly observed by the GNSS and InSAR measurements 30 . The slip distribution was tapered horizontally by the square root taper of the width of 4 km, to evaluate the stress state near the horizontal edges. The average stress drop on the mainshock fault was computed to be 3.4 MPa. The shear stresses on the fault and extension due to the slip distributions were calculated at depths of 5 km along the fault strike on a plane 100 m orthogonal to the fault by the programme of Nikkhoo and Walter (2015) 53 . We assumed a Poisson's ratio of 0.25 and a rigidity of 35 GPa. The reason for calculating this on a plane 100 m orthogonal to the fault was practically to avoid singularities on faults and the stress concentration at points very close to the fault edge.

Data availability
All seismic data needed to evaluate the conclusions in the paper are available in the Zenodo repository, https://zenodo.org/deposit/5036124#, and available from the Data