The aim of dynamic rupture inversion is to gain a better understanding the causality of the rupture process not only in spatio-temporal evolution but also in stress condition. Let us suppose that the rupture process is governed by linear slip-weakening friction (Fig. 4). The fault strength is defined as a function of on-going slip \({\Delta }u\) as follows,

$$\tau \left({\Delta }u\right)={\tau }_{r}+{\Delta }{\tau }_{b}H\left(1-\frac{{\Delta }u}{{D}_{c}}\right)$$

where the breakdown strength drop \({\Delta }{\tau }_{b}\) is the difference between the peak strength \({\tau }_{p}\) and the residual strength \({\tau }_{r}\), namely \({\Delta }{\tau }_{b}\equiv {\tau }_{p}-{\tau }_{r}\), and \(H(\bullet )\) is the Heaviside function. We consider a planar fault in an infinite homogeneous medium, so that normal stress play no explicit role and \({\tau }_{r}=0\) can be assumed without loss of generality. The parameter \({D}_{c}\) is called the critical slip displacement, characterizing the slip-weakening friction, and generally considered to be scale-dependent (Ohnaka, 2003; Aochi & Ide, 2004). The initial shear stress \({\tau }_{0}\) is set such that \({\tau }_{0}<{\tau }_{p}\). The rupture stability is often examined using non-dimensional parameters such as the S-value (Das & Aki, 1977) or \(\kappa\) (Madariaga & Olsen, 2000):

$$S=\frac{{\tau }_{p}-{\tau }_{0}}{{\tau }_{0}-{\tau }_{r}} \text{a}\text{n}\text{d} \kappa =\frac{{\left({\tau }_{0}-{\tau }_{r}\right)}^{2}}{G\bullet \left({\tau }_{p}-{\tau }_{r}\right)}\frac{L}{{D}_{c}}$$

where \(G\) is the rigidity of the medium and \(L\) is a characteristic length of the rupture process on a causal fault.

We prepare a small nucleation patch and a larger target patch to describe fault heterogeneity, similar to previous works (Aochi & Twardzik, 2020; Aochi & Ruiz, 2021). They are circular with radii \(R={r}_{0}\) and \(R={r}_{1}\) respectively for simplicity. We seek to determine the relative position and frictional parameters of the target patch (Fig. 4), on which uniform frictional parameters (\({D}_{c}, {\tau }_{p},{\tau }_{r}\)) are assumed. On the other hand, for the nucleation patch, we assume that its size is equal to half of the target, \({r}_{1}=2 {r}_{0}\) and that \({D}_{c}\) is proportional to the distance from the center (Aochi & Ide, 2004; Uemura et al., 2023). Consequently, the value of \({D}_{c}\) on the target patch is twice that at the edge of the nucleation patch. Outside of the patches, no weakening, \({\tau }_{p}={\tau }_{r}\), is supposed to gradually halt the rupture propagation. We then limit the number of simulations by discretizing key parameters in advance. We test three pairs of (\({\tau }_{p},{\tau }_{r}\)) and nine relative locations of the target patch with respect to the given nucleation patch, while following the \({D}_{c}\)-scaling adopted in Ide & Aochi (2005) and Aochi & Twardzik (2020). Numerical simulations of the dynamic rupture process are carried out using boundary integral equation method (Aochi et al., 2000) and renormalization technique (Aochi & Ide, 2004). The simulation starts with a grid size of \({\Delta }s\)=32 m for a fault size of 2048 m \(\times\) 2048 m, and when the rupture reaches at the edge of the model dimension, it is normalized to a larger scale by a factor of 4, i.e. to \({\Delta }s\)=128 m and again to \({\Delta }s\)=512 m. The radius of the target patch \({r}_{1}\) is adjusted to obtain the magnitude expected to Mw6, i.e. \({r}_{1}\) = 4.92 km, 3.93 km and 2.95 km for the high, medium and low stress cases, respectively. Synthetic seismograms are calculated on the basis of dynamic rupture source models using a finite difference method (Aochi & Madariaga, 2003) supposing a regional 1D structure (Nishimura et al., 2023).

Figure 5 shows the comparison of the simulated slip area for the nine different target patch positions in the case of (\({\tau }_{p},{\tau }_{r}\)) = (10 MPa, 6 MPa). The magnitude ranges from Mw 5.98 (case 00020) to Mw6.04 (case 00022, 00024, 00026 and 00028). Seismograms from the six stations are compared in Fig. 6. As in the focal mechanism inversion, we use only the beginning of the waveforms for 10 seconds in the frequency range between 0.05 and 0.5 Hz. Here we test two possible fault planes (southeast or northwest dipping), three stress levels and nine target patch locations. We run the procedure twice, i.e., using the two different locations (#4 and #5 in Fig. 2) and two slightly different focal mechanisms. The focal mechanism of the second run is the one obtained at position #7, where the maximum slip is obtained in the model 00021 of the first run. The trend in solution convergence is very similar. In general, convergence is better when the stress level is high, namely the ruptured area is smaller, and when the rupture directivity is the same as the rake direction (numbered as model 11, 21, 31). Overall, models 00021 and 00031 of the southeast dipping fault are the two best, indicating a relatively higher stress drop and a localized ruptured area.

Figure 6 compares the vertical displacement at station ISK001 at which a temporal GPS station has been installed nearby showing uplift of about a few tenths’ centimeters before the arrival of the mainshock deformation (Nishimura, 2024). It is briefly confirmed by the double integration of the acceleration data. However, the acceleration data may be subject to drift and do not show the final stable displacement after the arrival of the strong ground shaking from the M7.6 mainshock. For this reason, static displacement is not used in the previous inversion. However, from a qualitative point of view, the high-stress model does not retain sufficient final displacement, because the ruptured area is small, namely far from station ISK001. Only the two models of low and medium stress values with the appropriate directivity can leave the final displacement of 20 cm, which is consistent with the insight from GPS measurements.