Driving stress and seismotectonic implications of the 2013 Mw5.8 Awaji Island earthquake, southwestern Japan, based on earthquake focal mechanisms before and after the mainshock

A driving stress of the Mw5.8 reverse-faulting Awaji Island earthquake (2013), southwest Japan, was investigated using focal mechanism solutions of earthquakes before and after the mainshock. The seismic records from regional high-sensitivity seismic stations were used. Further, the stress tensor inversion method was applied to infer the stress fields in the source region. The results of the stress tensor inversion and the slip tendency analysis revealed that the stress field within the source region deviates from the surrounding area, in which the stress field locally contains a reverse-faulting component with ENE–WSW compression. This local fluctuation in the stress field is key to producing reverse-faulting earthquakes. The existing knowledge on regional-scale stress (tens to hundreds of km) cannot predict the occurrence of the Awaji Island earthquake, emphasizing the importance of estimating local-scale (< tens of km) stress information. It is possible that the local-scale stress heterogeneity has been formed by local tectonic movement, i.e., the formation of flexures in combination with recurring deep aseismic slips. The coseismic Coulomb stress change, induced by the disastrous 1995 Mw6.9 Kobe earthquake, increased along the fault plane of the Awaji Island earthquake; however, the postseismic stress change was negative. We concluded that the gradual stress build-up, due to the interseismic plate locking along the Nankai trough, overcame the postseismic stress reduction in a few years, pushing the Awaji Island earthquake fault over its failure threshold in 2013. The observation that the earthquake occurred in response to the interseismic plate locking has an important implication in terms of seismotectonics in southwest Japan, facilitating further research on the causal relationship between the inland earthquake activity and the Nankai trough earthquake. Furthermore, this study highlighted that the dataset before the mainshock may not have sufficient information to reflect the stress field in the source region due to the lack of earthquakes in that region. This is because the earthquake fault is generally locked prior to the mainshock. Further research is needed for estimating the stress field in the vicinity of an earthquake fault via seismicity before the mainshock alone.


Introduction
On April 13, 2013 at 05:33 (local time), an M w 5.8 (M j 6.3) shallow inland earthquake occurred on Awaji Island in southwestern Japan (hereafter referred to as ' Awaji Island earthquake' or mainshock) (Fig. 1). Here, M j is the Open Access *Correspondence: imani@ni.aist.go.jp 1 Geological Survey of Japan, AIST, Tsukuba Central 7, 1-1-1 Higashi, Tsukuba, Ibaraki 305-8567, Japan Full list of author information is available at the end of the article Imanishi et al. Earth, Planets and Space (2020) 72:158 magnitude of the earthquake determined by the Japan Meteorological Agency (JMA). The focal mechanism of the mainshock (Fig. 1) and its aftershock distribution ( Fig. 2) indicate that it was caused by the movement of a western-dipping reverse fault, which was previously unrecognized. Strong ground motions during the mainshock were widely observed and the JMA seismic intensity of 6-(difficult to keep standing) was recorded at Awaji Island. The earthquake damaged more than 8000 houses and injured 35 people (Fire and Disaster Management Agency 2013), showing that even for a moderatesize event, shallow inland earthquakes have the potential to cause significant damage. If earthquakes occur in a densely populated area, they will create significant damage, such as the 2011 M w 6.2 (M L 6.3) Christchurch earthquake in New Zealand (Potter et al. 2015). Therefore, in addition to large earthquakes, it is important to investigate the generation mechanism of moderate earthquakes through detailed analyses. The Awaji Island earthquake occurred adjacent to the southwestern end of the disastrous 1995 M w 6.9 (M j 7.3) Kobe earthquake (Fig. 1). The focal mechanism of the Awaji Island earthquake is quite different from the Northeast-Southwest (NE-SW) trending right-lateral strike-slip Kobe earthquake, suggesting the existence of a local-scale driving stress variation. The epicenter of the mainshock was slightly shifted southwest from the extension of the fault zone of the Kobe earthquake ( Fig. 1). This observation indicates that the Awaji Island earthquake did not correspond to a direct aftershock of the Kobe earthquake. However, its occurrence was influenced by the Kobe earthquake via stress transfer; thus, it is necessary to quantitatively elucidate this effect as one of the factors causing the Awaji Island earthquake.
The occurrence of the Awaji Island earthquake has important implications in terms of seismotectonics in southwestern Japan, which constitutes a part of the Philippine Sea plate. The number of shallow inland earthquakes in the Kinki region ( Fig. 1), including the Awaji Island, increased from about 50 year before to 10 year after the occurrence of large interplate earthquakes along the Nankai trough (Utsu 1974;Hori and Oike 1999). Since the Awaji Island earthquake was captured by a dense seismic observation network (Okada 2004), this earthquake is an unusual case of inland earthquakes during the cycle of the Nankai trough earthquake. A detailed study of the earthquake will explain the causal relationship between the inland earthquake activity and the Nankai trough earthquake.
In this study, we focused on the driving stress of the Awaji Island earthquake to understand its generation mechanism. For this reason, we first determined the focal mechanisms  Nakata and Imaizumi (2002). Circles represent M0 + seismicity shallower than 20 km determined by the Japan Meteorological Agency (JMA) during the period from 1 January 2003 to 12 April 2013 (i.e., the day before the Awaji Island earthquake). The 'beach balls' represent the centroid moment tensor solutions by JMA. The surface projection of the fault model of the Kobe earthquake (Ide and Takeo 1997) is indicated by a rectangle. Topography is based on a 250 m grid digital elevation model published by the Geospatial Information Authority of Japan of the aftershocks of the Awaji Island earthquake, as well as earthquakes before the mainshock, by using the seismic records from regional high-sensitivity seismic stations. The stress tensor inversion method was then applied to infer the stress fields in the source region. Furthermore, we computed the Coulomb stress change and its evolution from the Kobe earthquake using a viscoelastic structure as an additional loading stress on the mainshock fault. We also discuss the seismotectonic implications of the Awaji Island earthquake in terms of the earthquake cycle in the Nankai trough and the origin of local stress heterogeneity. Figure 3a shows the distribution of the permanent seismic stations used in this study. The stations are operated by the National Research Institute for Earth Science and Disaster Resilience (NIED), Japan Meteorological Agency (JMA), Geological Survey of Japan, AIST (GSJ), Earthquake Research Institute, the University of Tokyo (ERI), and Disaster Prevention Research Institute, Kyoto University (DPRI). Each station is equipped with a set of three-component velocity transducers, which have a natural frequency of 1 or 2 Hz. Seismometers operated by the NIED and GSJ are installed at the bottom of boreholes (Okada et al. 2004;Imanishi et al. 2011b). We analyzed 188 earthquakes which had a JMA magnitude (M j ) ≥ 1.5 and focal depth ≤ 20 km between January 2003 and August 2013. The dataset consists of 31 earthquakes before the mainshock, the mainshock, and 156 aftershocks, whose locations are shown as circles with different colors in Fig. 3b. Out of all the data, the largest aftershock occurred 8 min after the mainshock and had a a b Hypocenter and focal mechanism of 2013 Awaji Island earthquake sequence

Hypocenter determination
We determined the hypocenters of the earthquakes by applying a maximum-likelihood estimation algorithm (Hirata and Matsu'ura 1987) to the picked arrival times. The P-and S-wave arrival times were identified manually. We used a one-dimensional velocity model (Ide and Takeo 1997), which consists of four layers (Table 1). We first located the hypocenters of all events without station corrections. We then relocated them by introducing the station corrections, which were obtained using the average of the differences between the observed and the theoretical travel times at each station. To obtain the final locations, this procedure was repeated until the reduction in the root mean square of the arrival time residuals (RMS) became less than 0.01 s. After three iterations, the RMS decreased from 0.20 to 0.06 s for the P-wave and from 0.52 to 0.12 s for the S-wave. We then applied the double-difference earthquake relocation algorithm (hypoDD) (Waldhauser and Ellsworth 2000) to the traveltime differences formed from the data using the same velocity structure. Each event is linked to their neighbors through commonly observed phases with a maximum hypocentral separation between the linked events (5 km). This resulted in 22,296 P-wave pairs and 28,807 S-wave pairs. Figure 4 shows the hypocenter distribution determined in this study. In comparison with the distribution determined by JMA, our DD location is clustered and shows a sharper view of seismicity. The aftershocks define the western-dipping plane with a dip angle of ~ 60° that extends between ~ 12 and 17 km in depth and 10 km along the strike of the plane (Fig. 4b), which is consistent with the focal mechanism solution of the mainshock. The mainshock hypocenter is located almost in the middle of  the aftershock distribution. Although this distribution can also be observed in the JMA catalog, our DD location makes its feature clearer. The pre-shocks were dispersed and did not overlap with the aftershock distribution, suggesting that the pre-shocks alone cannot illuminate a fault plane associated with the mainshock. Based on the field surveys and the structural analysis of the fault rocks, Lin et al. (2015) reported that there is an active fault in the source region (Yamada fault) that has been active since the late Pliocene. They inferred that the Awaji Island earthquake occurred along the fault. However, Fig. 4b shows that the extension of a western-dipping plane to the surface does not match the surface trace of the Yamada fault. We conclude that the Yamada fault did not rupture in the Awaji Island earthquake.

Focal mechanism solutions using body wave amplitudes
We determined the focal mechanism solutions of the pre-shocks and aftershocks using P-wave polarity with absolute P-and SH-wave amplitudes, following the a b Fig. 4 Hypocenter distributions determined in the present study. Blue, hypomh (Hirata and Matsu'ura 1987); Red and yellow star, hypoDD (Waldhauser and Ellsworth 2000); black, JMA catalogue. a Pre-shock events. b The mainshock and aftershocks. Yamada fault reported by Lin et al. (2015) is shown by an orange line on map procedure of Imanishi et al. (2011a). We first determined the focal mechanisms and seismic moments for earthquakes with at least eight P-wave polarities. We then obtained the amplitude correction for each station by calculating the logarithmic average of the ratios between the synthetic and observed amplitudes of all analyzed events. Using these station corrections, the focal mechanism solutions and seismic moments were re-determined. We tested the stability of the solution by plotting all focal mechanisms whose residuals were ≤ 1.1× the minimum residual value. We rejected the ambiguous solutions where multiple solutions were possible. We defined focal mechanism uncertainties for each event based on the average of the Kagan angles (Kagan 1991) between the best-fitting solution and all the solutions that had residual values less than 1.1 times the minimum residual value. Here, the Kagan angle becomes 0° when the two mechanisms are the same and 120° when they differ the most. The average uncertainty of all the focal mechanisms was 11°.
In total, we obtained 28 pre-shock and 134 aftershock focal mechanisms whose moment magnitudes ranged from 1.7 to 3.1 and from 1.6 to 3.4, respectively. The relationship between M w and M j is shown in Fig. 14.
The focal mechanisms of pre-shocks and aftershocks are plotted on the maps in Fig. 5a and b, respectively. We classified these mechanisms as reverse (green), strike-slip (red), or normal (blue) faulting earthquakes according to the definitions by Frohlich (1992). For the pre-shocks, strike-slip earthquakes are prominent. In contrast, reverse-faulting type aftershocks occur throughout the aftershock area, which is consistent with the mainshock focal mechanism. A considerable number of aftershocks with strike-slip components also occur. Katao et al. (2014) determined the P-wave first motion focal mechanism solutions of the aftershocks by adding the data from 10 temporary seismic stations, installed after the mainshock. The spatial pattern of the focal mechanisms determined in this study agreed with their result, suggesting that even if there are few seismic stations above the source, reliable solutions can be obtained by including the amplitude data in the focal mechanism determination.
The directions of the P-axes are shown in Fig. 6. Most P-axes are sub-horizontal and oriented in the E-W direction, which conforms to the general tectonic trend around this area (Townend and Zoback 2006;Terakawa and Matsu'ura 2010;Matsushita and Imanishi 2015;Yukutake et al. 2015). Looking closely, however, the P-axis orientation locally changes in space, specifically, the WSW-ENE direction in the main part of the aftershock area, and WNW-ESE in the surrounding region. This is further discussed quantitatively via stress tensor inversion analysis in "Stress field in and around the Awaji Island earthquake" section.

Accurate fault geometry of mainshock by incorporating aftershock focal mechanisms
The aftershock distribution shows the mainshock fault plane western-dipping with a dip angle of ~ 60° (Fig. 4b).
In order to define the fault geometry more accurately, which is necessary for the following analyses, we incorporated the aftershock focal mechanisms into the determination. We focused on Kagan's angles (Kagan 1991), which are the minimum rotation angles of the aftershock focal mechanisms relative to a reference mechanism corresponding to the double-couple component derived from the centroid moment tensor (CMT) solution of the mainshock (strike = 175°, dip = 60°, rake = 95°) (Fig. 1). In Fig. 7, Kagan's angles of the focal mechanisms ≤ 30° (> 30°) are plotted as blue (red) circles. Except for one aftershock, the blue circles define a clear alignment, illuminating the mainshock fault plane. We infer that the blue circles, except for the easternmost one, correspond to the on-fault aftershocks and the others (the red circles and the easternmost blue circle) correspond to the off-fault ones. The geometry of the on-fault aftershocks defines a mainshock rupture plane that strikes N175° E and dips at 62°. Based on this fault geometry, we obtained the mainshock fault slip model (see Appendix B). The observed seismograms can be explained by almost pure reverse-faulting slips near the hypocenter (Fig. 15).

Stress field in and around the Awaji Island earthquake
We determined the stress field in and around the source region from the population of focal mechanisms by applying the stress tensor inversion method from Michael (1984). The parameters solved for in the inversion are: (1) the orientation of the three principal stress axes and (2) the relative magnitude of the principal stresses defined by φ = (σ 2 − σ 3 )/(σ 1 − σ 3 ) , where σ 1 , σ 2 , and σ 3 are the maximum, intermediate, and minimum compressive principal stresses, respectively. One plane must be chosen from each focal mechanism as the actual fault plane, because the inversion uses the direction of the tangential traction on the fault plane. We select the fault plane in each focal mechanism that minimizes the misfit angle, the difference in the angle between observed and theoretical slip directions, while inverting for the stress tensor based on the grid search algorithm from Michael (1987). This algorithm is considered a rational approach in situations where the choice of fault plane cannot be made based on relative hypocentral locations or geological Focal mechanism solutions determined in this study (lower hemisphere of equal-area projection). a Pre-shock focal mechanism solutions, where different colors are used to differentiate reverse (green), strike-slip (red), and normal (blue) faulting mechanisms. A triangle diagram (Frohlich 1992) with the color scale is presented in the right, where each focal mechanism is plotted by open circles. The number within each set of parentheses shows the number of events for that faulting mechanism. Beach balls with an asterisk in the upper right indicate a reverse-faulting earthquake. b Distribution of aftershock focal mechanisms for each faulting mechanism. The map area corresponds to the red rectangle in a. The hypoDD location of the mainshock is shown by a yellow star information (Kastrup et al. 2004). To calculate confidence regions for the stress tensor, we applied the bootstrap resampling technique by assuming that each nodal plane had the same probability of being selected during the resampling. Following Michael (1987), we used 2000 bootstrap samples to obtain the 95% confidence region.
As for the pre-shocks, we divided the analyzed earthquakes into two regions: within the main part of the source area and other (Fig. 8a). This division is based on the observation of spatial change in the P-axis orientation described in "Focal mechanism solutions using body wave amplitudes" section. In the following, we call the former earthquakes as 'Pre-shock2′ and the latter as 'Pre-shock1′. The stress tensor inversion results are shown in Fig. 8b. The ratio of earthquakes with large misfit angles is high for aftershocks, which may reflect the presence of small-scale stress heterogeneity caused by the mainshock. According to Michael (1991), the quality of the results can be examined based on the average misfit angle ( β ) for a given amount of focal mechanism errors. In the case of the present study (focal mechanism errors of 11°), we can conclude that the uniform stress component is adequately retrieved within a 95% level, when β does not exceed ~ 36° [see Fig. 8 of Michael (1991)]. The values of β are 10° (Pre-shock1), 7° (Pre-shock2), and 24° (aftershocks), suggesting a valid fit to the single spatially uniform stress tensor even in aftershocks.
Stress fields estimated from the pre-shock focal mechanisms (i.e., Pre-shock1 and Preshock2) indicate that both are characterized by a strike-slip faulting stress regime but by different horizontal stress orientations, where the maximum compressive stress ( σ 1 ) is oriented in the ENE-WSW direction within the source area (Pre-shock2) and in the ESE-WNW direction outside the source area (Pre-shock1) (Fig. 8c). The stress ratio φ is also different in two datasets, suggesting that there is a notable difference in the stress field inside and outside the source area.
Regarding the stress field derived from aftershock focal mechanisms, σ 1 is oriented in the ENE-WSW direction, which is consistent with that derived from Pre-shock2 (Fig. 8c). In addition, the observation of a low stress ratio ( φ = 0−0.3 ) is much more similar to the case of Pre-shock2 than Pre-shock1. The 95% confidence regions of σ 2 and σ 3 form a girdle and trend roughly in the NNW-SSE direction. Together with the observation of a low Kagan's angle (Kagan 1991) of aftershock focal mechanisms relative to the mainshock. a Map view and b cross-section. Blue and red circles indicate Kagan's angle below and above 30°, respectively. The hypoDD location of the mainshock is shown by a yellow star stress ratio, these confidence regions indicate that principal stresses σ 2 and σ 3 differ only slightly in magnitude and that the stress field is characterized by a mixture of reverse and strike-slip faulting. Except for the faulting style, the stress fields within the source area, estimated from the two datasets (before and after the mainshock), have common characteristics. Although the temporal stress change is often attributed to coseismic stress perturbation (Hardebeck and Okada 2018), as shown later ("Driving stress inferred from slip tendency analysis" section), we infer that the difference in the faulting style is an artifact caused by incomplete data sampling before the mainshock.
For the stress tensor inversion, we also tried the approach from Vavryčuk (2014) to check the robustness of the result in Fig. 8. This method is a modification of Michael (1984Michael ( , 1987, in which the fault planes are chosen by applying the fault instability constraint and the stress tensor is determined in iterations. The uncertainties are calculated as the maximum differences between the results of the inversion for noise-free and noisy data with 1000 noise realizations, requiring assumptions with regard to focal mechanism error. Although Vavryčuk (2014) suggests that Michael's method yields a biased value of the stress ratio, we confirmed that, for the present datasets and focal mechanism error of 11°, both methods return almost identical results.

Driving stress inferred from slip tendency analysis
We evaluated whether the stress fields estimated from Pre-shock2 and aftershocks, both of which are located within the source area, could cause an almost pure reverse-faulting Awaji Island earthquake. We used the normalized slip tendency defined as T s ′ = T s / max(T s ) = T s /μ, where T s is the ratio of the shear stress to the effective normal stress on the fault plane and μ is the static frictional coefficient (Lisle and Srivastava 2004). The value of T s ′ ranges from 0 to 1, where larger values correspond to a greater slip potential. We supposed that the frictional sliding envelope is tangential to the Pre-shock2 aftershocks Pre-shock1 (15) Pre-shock2 (13) aftershocks (   The left-hand panels show the principal stress axes with their 95% confidence regions plotted on the lower hemisphere stereonets. The middle panels show the misfit angle for the data with respect to the optimal stress tensor determined by the stress tensor inversion. Here, the misfit angle represents the angle between the tangential traction predicted by the optimal solution and the observed slip direction on each plane, as determined from the focal mechanism. The right-hand panels show the frequency of the stress ratio φ, which belongs to the 95% confidence region. c Frequency of the azimuth of σ 1 from the 95% confidence region of the stress tensor inversion. Here, 180° is subtracted from azimuth, if necessary, to ensure that all orientations can be plotted in the range of 0° to 180° σ 1 -σ 3 Mohr circle, which makes it possible to compute T s ′ without knowing the absolute stress value (Lisle and Srivastava 2004). Collettini and Trippetta (2007) defined the favorably oriented planes of the slip by 0.5 < T s ′ ≤ 1 and the unfavorably oriented planes by 0 ≤ T s ′ ≤ 0.5. We assumed the fault geometry estimated by the on-fault aftershocks (Fig. 7) and set the static friction coefficient to 0.6. The slip angle was estimated based on the Wallace-Bott hypothesis (Bott 1959;Wallace 1951) that slip on faults occurs in the direction of the maximum resolved shear stress. Figure 9a presents a histogram of the T s ′ values computed from the 95% confidence regions of the assumed stress tensor, where the green and black bars correspond to the results of the stress fields derived from Pre-shock2 and aftershocks, respectively. The results indicate that the stress field estimated from the aftershocks is favorably oriented to the mainshock fault. Figure 9b presents a histogram of slip angles inferred from the stress tensors included in the 95% confidence region. Note that the inferred slip angles have a peak at 100°-110° (101° in the case of the optimal stress tensor), which is consistent with the result of the slip inversion (see Fig. 15b). In contrast, the stress field estimated from the Pre-shock2 is unfavorably oriented to the mainshock fault. This is a natural consequence because it is difficult for faults to be reactivated under a strike-slip faulting regime in which the maximum principal stress is perpendicular to the strike of the fault. Naturally, faults can be reactivated under an extremely low friction coefficient or an extremely high fluid pressure as is shown with the San Andreas Fault (Zoback et al. 1987). However, they do not always cause reverse faulting. In fact, the inferred slip angle has a broad distribution. We suggest that the dataset before the mainshock does not have sufficient information to reflect the stress field in the source region due to the lack of earthquakes in the region, specifically because the earthquake fault was locked prior to the mainshock.
Based on the above considerations, we conclude that the stress field estimated from the aftershocks is appropriate for the driving stress of the Awaji Island earthquake where the source fault originally had a high slip potential of reverse faulting. To further confirm this, we investigate whether the inferred driving stress is compatible with the occurrence of pre-shocks in the source area, i.e., Pre-shock2. We calculated the misfit angle of each focal mechanism of Pre-shock2 with respect to the optimal stress tensor derived from aftershocks, where we selected the lowest misfit angle between the two nodal planes. Figure 10 presents a histogram and spatial distribution of misfit angles. All earthquakes have Pre-shock2 Fig. 9 Driving stress of the Awaji Island earthquake. a Slip tendency analysis for the Awaji Island earthquake fault. The frequency of the normalized slip tendency (T s ') is computed from the 95% confidence region of the assumed stress tensors. The black, blue, and green colors show the results obtained assuming the stress tensors derived from aftershocks, Pre-shock1, and Pre-shock2, respectively. b Slip angle expected by the stress tensor derived from aftershocks, Pre-shock1, and Pre-shock2. The meaning of colors is the same as in a. The frequency is computed from the stress tensors included in the 95% confidence region. Inverted triangle shows the average rake angle of mainshock (Appendix B) Misfit angles of Pre-shock2 with respect to the optimal stress tensor derived from aftershocks (i.e., a driving stress of Awaji Island earthquake). a Frequency. b Spatial distribution small misfit angles less than 30°, which indicates that the inferred driving stress is also adequate for explaining the data in the source area before the mainshock.

Stress transferred by 1995 M W 6.9 Kobe earthquake
Regarding an additional loading stress of the Awaji Island earthquake, we evaluated the Coulomb stress change ( CFF ) and its evolution after the Kobe earthquake, following the approach from Ohtani and Imanishi (2019). We assumed a one-dimensional viscoelastic structure used in Shikakura et al. (2014), in which the thickness of the elastic layer is 35 km and the Maxwell viscosity of the viscoelastic layer is 5.0 × 10 18 Pa·s. For the fault slip of the Kobe earthquake, we used a finite fault source model derived from the inversion of the near-source strong-motion data (Ide and Takeo 1997).
Most of the slip exists beneath Awaji Island, which is a feature commonly found in other slip models (Wald 1996;Koketsu et al. 1998;Sekiguchi et al. 2000). We assumed the receiver fault of the Awaji Island earthquake based on the on-fault aftershocks and the slip model; i.e., (strike, dip, slip) = (175°, 62°, 101°). An apparent friction coefficient of 0.4 was used for the CFF calculation. Figure 11a shows the computed coseismic CFF at the depth of 14.4 km (relocated hypocenter of the mainshock), which suggests that the Awaji Island earthquake lies in a positive Coulomb stress change area. The coseismic CFF value at the mainshock hypocenter is 0.09 MPa, which exceeds the earthquake triggering threshold of 0.01 MPa (Stein 1999). The 18-year delay suggests that the coseismic positive CFF was insufficient to trigger the slip on the fault plane of the Awaji Island earthquake at the time. Figure 11b shows the postseismic change in CFF at the hypocenter, which decreased over time. This indicates that a gradual stress change by viscoelastic relaxation cannot explain the occurrence of the Awaji Island earthquake. The magnitude of the stress decrease is 3 kPa in 18 year (3% of the coseismic CFF ), which is insufficient to return the stress level to that before the Kobe earthquake. In other words, the earthquake was still likely to occur. Saito et al. (2018) computed the stress change at the depth of 10 km due to the interseismic plate locking along the Nankai trough (slip deficit) estimated using the GNSS data . This showed that around the source region, a reverse-faulting stress change with the σ 1 -direction, similar to the driving stress of the Awaji Island earthquake, is induced in the order of a few kPa/year [see Fig. 5 of Saito et al. (2018)]. With this in mind, we inferred that the stress change caused by the slip deficit contributed to the occurrence of the Awaji Island earthquake, pushing the earthquake fault over its failure threshold in 2013.

Existence of local stress heterogeneity in the source area
The stress tensor inversion revealed that the region outside the source area is characterized by a strikeslip faulting stress regime with σ 1 oriented in the ESE-WNW direction sub-horizontally. This stress field is widely observed throughout the region (Townend and Zoback 2006;Terakawa and Matsu'ura 2010;Matsushita and Imanishi 2015;Yukutake et al. 2015). In contrast, σ 1 rotates in a counter-clockwise direction about 30° within the source area and becomes perpendicular to the strike of the mainshock fault plane (Fig. 8). This is an important factor in the activation of reverse-faulting Awaji Island earthquake. Strictly speaking, the σ 1 -direction is a little oblique to the strike of the fault, which explains why the slip direction of the mainshock slightly exceeds 90° (Fig. 15). To demonstrate how the stress field in the source area locally deviates from the surrounding region, we show the distribution of the maximum principal compressive stress direction (S Hmax ) over a wider area. We computed the maximum horizontal compressive stress directions (S Hmax ) from the optimal stress parameters, determined via stress tensor inversion, using the method of Lund and Townend (2007). In the present case, the orientation of S Hmax is almost consistent with that of σ 1 . Figure 12a compiles S Hmax orientations determined in the present study and those reported by Townend and Zoback (2006) and Matsushita and Imanishi (2015). Figure 12b illustrates that the orientation observed throughout the region is prominent from E-W to ESE-WNW. Although local fluctuations can be identified, we can confirm that the stress field in the source area locally deviates from the surrounding region.
Here, we investigate whether a regional stress could cause the Awaji Island earthquake. Assuming that the stress field estimated from the Pre-shock1 is representative of regional stress, we computed T s ′ and the slip angle for the mainshock fault plane via the same approach as in "Driving stress inferred from slip tendency analysis" section. The result shows that the regional stress field is unfavorably oriented to the mainshock fault and the expected slip motion is left-lateral (blue bars in Fig. 10). This indicates that the regional stress field cannot move the mainshock fault as a reverse fault, supporting the existence of local stress heterogeneity in the source area. This result demonstrates that the regional-scale   Townend and Zoback (2006) and Matsushita and Imanishi (2015). Circles in the right map correspond to earthquakes used in the stress tensor inversion (see Fig. 8). b Rose diagram of S Hmax orientation shown in a. The red and blue lines indicate S Hmax orientation within and outside the source area determined in the present study, respectively. Regarding the orientation within the source area, the average value is shown stress field (tens to hundreds of km) is not sufficient for investigating the generation mechanism of M6-class earthquakes, highlighting the importance of detailed local-scale (< tens of km) stress information.

Origin of local stress heterogeneity
What are the factors that created the local stress heterogeneity in the source area of the Awaji Island earthquake? One possible factor is a stress anomaly caused by the Kobe earthquake. As shown above, however, the coseismic stress change at the mainshock hypocenter was at an order of 0.1 MPa, which is only 3% of the coseismic stress drop of the Awaji Island earthquake (Appendix B). It is considered that there was almost no contribution to the local stress field by the Kobe earthquake. According to the geological structure of this region, there is a south-eastern uplift structure (the Ichinomiya flexure) on the north-eastern side of the source region and a northwestern uplift structure (the Ayuhara flexure) on the southwestern side (Takahashi et al. 1992) (Fig. 13). The local tectonic movement that created these two opposite flexures would have generated local stress heterogeneity, similar to the case of stress anomaly produced by the interactions of discontinuous parallel fault segments (Crider and Pollard 1998;Vavryčuk and Adamova 2018). Since there is no information on the geometry of faults that created the flexures nor on the amount of accumulated slip, a quantitative discussion is not possible. Depending on the setting of unknown parameters, however, it would be possible to locally rotate the maximum principal stress axis counter-clockwise with respect to the surroundings. It is important to note that the mechanism alone does not always produce a local reversefaulting stress field. Ogata (2013) detected anomalies in seismicity and crustal deformation around the source during the period leading up to the Awaji Island earthquake. The research shows that these anomalies are well explained by an aseismic reverse fault slip at the deeper extension of the fault. Based on the analysis of focal mechanism solutions, Matsumoto et al. (2015) found that the stress concentration occurred around the fault before the mainshock, suggesting that the possible cause was a deep aseismic slip. We expect that a transient aseismic slip repeatedly occurs at the deeper extension of the fault, such as at the plate boundary (Obara and Kato 2016). As a result, a reverse-faulting stress field can be locally formed over a long period of time because a deeper slow slip transfers stress to the shallower part. The summary  (Nakata and Imaizumi 2002), and the black lines show Ichinomiya and Ayukawa flexures (Takahashi et al. 1992). The red and black arrows show the orientation of the maximum principal stress in the source region and in the surrounding area, respectively. Blue circles represent the on-fault aftershocks. The line X-Y is perpendicular to the strike of the mainshock fault plane. Topography is based on a 10 m grid digital elevation model published by the Geospatial Information Authority of Japan of local seismotectonics causing stress heterogeneity is shown in Fig. 13.

Seismotectonic implications of the Awaji Island earthquake
We explore the seismotectonic implications of the Awaji Island earthquake, together with the implications of the 2018 M w 5.8 northern Osaka earthquake (Fig. 1). As mentioned in "Introduction" section, the activity of shallow inland earthquakes in the Kinki region increased from about 50 year before to 10 year after the occurrence of large interplate earthquakes along the Nankai trough (Utsu 1974;Hori and Oike 1999). Mitogawa and Nishimura (2020) showed that a model incorporating the viscoelasticity successfully reproduces the activation of shallow inland earthquakes, not only before but also after the interplate earthquakes. Under a more realistic model setting, Shikakura et al. (2014) explained this feature by computing the long-term viscoelastic Coulomb stress change on inland earthquake faults. They incorporated all possible deformation sources: the E-W compression observed as the Niigata-Kobe Tectonic Zone (Sagiya et al. 2000), historical interplate earthquakes along the Nankai trough, interseismic plate locking, and historical inland earthquakes in southwestern Japan. They found that the frequency of N-S trending reverse-faulting earthquakes increases before the occurrence of interplate earthquakes. Contrastingly, the frequency of NW-SE trending left lateral and NE-SW trending right-lateral strike-slip earthquakes increases after interplate earthquakes occur. Considering the recurrence interval of 100-200 year along the Nankai trough (Headquarters for Earthquake Research Promotion 2013) and the elapsed time since the last earthquakes (1944 Showa-Tonankai and1946 Showa-Nankai earthquakes), the occurrence of the 2013 reverse-faulting Awaji Island and the 2018 northern Osaka earthquakes (see Appendix C) may indicate that the Nankai trough is now in the later stage of the earthquake cycle and the Kinki region encounters the active period before the next Nankai event. The present study suggests that stress build-up, due to the interseismic plate locking, largely contributed to the occurrence of the two earthquakes. However, these two earthquakes alone are not sufficient to reach this conclusion. It is important to investigate whether the N-S trending reverse-faulting earthquakes occur more frequently than other types of earthquakes, while also considering the microearthquakes.

Caution for a stress field inferred from pre-shocks alone
In the present study, we demonstrate that the stress field estimated from the pre-shocks is not suitable for the driving stress of the Awaji Island earthquake. We infer that the earthquake fault was locked before the mainshock; therefore, the earthquake data around the most important location were incomplete. The seismic activity in the source fault is generally low before the occurrence of the earthquake; therefore, the same may be true for other faults. This raises a serious problem when evaluating the slip potential of future earthquakes and estimating the rupture scenario from dynamic simulations, based on the stress field estimated from earthquakes before the mainshock. This also indicates that caution must be exercised when estimating absolute stress level based on the temporal stress change caused by earthquakes (see Hardebeck and Okada 2018, and references therein), which requires exact sampling of earthquakes at the same location before and after the earthquakes. Based on a close examination of Fig. 5a, we can identify three reverse-faulting earthquakes within the source area, which suggests that an accurate stress field can be estimated using only pre-shocks. In the present study, we specifically analyzed earthquakes of M j ≥ 1.5. We did not consider most of the earthquakes for stress field estimation (compare blue circles in Figs. 2 and 3) because it is difficult to determine stable focal mechanism solutions of relatively small earthquakes. There is a sophisticated approach that determines the stress tensor, not from previously obtained focal mechanisms, rather from original datasets such as P-wave polarities and amplitude data (Abers and Gephart 2001;Horiuchi et al. 1995;Imanishi 2018;Iwata 2018). In this approach, the stress field can be determined even if only a small number of P-wave polarities and/or amplitude data are available from a single event. The approach is one potential strategy of addressing the challenge of incomplete datasets in the source fault before the earthquake, and will contribute to the development of more accurate seismic-hazard assessment before future earthquakes occur.

Conclusion
In the present study, we investigated the driving stress of the 2013 M w 5.8 Awaji Island earthquake using preshock and aftershock focal mechanism solutions. The stress tensor inversion using two datasets (pre-shock or aftershock focal mechanisms) revealed that the direction of the maximum principal stress within the source area locally deviates from the surrounding area, which as a result becomes perpendicular to the strike of the mainshock fault plane. The faulting style differs depending on whether the pre-shock or aftershock focal mechanisms are used; the former indicates a strike-slip faulting and the latter contains a reverse-faulting component. Since it shows a high slip potential in the slip direction of the mainshock slip model, the stress field estimated from the aftershocks can be considered as the driving stress of the Awaji Island earthquake where the stress field locally contains a reverse-faulting component with ENE-WSW compression. We suggest a possibility that local tectonic movement, i.e., the formation of flexures in combination with recurring deep aseismic slips, is one of the factors that have formed the local-scale stress heterogeneity (Fig. 13). The existing knowledge on regional-scale stress (tens to hundreds of km) could not predict the occurrence of the reverse-faulting Awaji Island earthquake, emphasizing the importance of the stress information at a local scale (< tens of km). The coseismic Coulomb stress change, caused by the 1995 Kobe earthquake, increased along the fault plane of the Awaji Island earthquake (0.09 MPa), and promoted the occurrence of the Awaji Island earthquake. Contrastingly, the postseismic Coulomb stress change was negative (3 kPa in 18 year). We conclude that the gradual stress build-up, due to interseismic plate locking along the Nankai trough, overcame the postseismic stress reduction in a few years, pushing the Awaji Island earthquake fault over its failure threshold in 2013. The observation that the earthquake occurred in response to the interseismic plate locking will facilitate further research on the relationship between the inland earthquake activity and the Nankai trough earthquake.
The present study highlights the problem that the dataset before the mainshock may not have sufficient information to reflect the stress field in the source region. This is because of a lack of earthquakes in the region because the earthquake fault is generally locked prior to the mainshock. Caution should be exercised when estimating the stress field in the vicinity of the earthquake fault based on seismicity before the mainshock alone.

Authors' contributions
KI and TU designed the present study. MO calculated the Coulomb stress change due to the Kobe earthquake considering the viscoelastic relaxation. TU determined the slip model of the Awaji Island earthquake. KI performed the remaining parts of the analysis and drafted the manuscript. All authors discussed the results and contributed to the final manuscript. All authors read and approved the final manuscript.

Availability of data and materials
The Japan Meteorological Agency's earthquake catalog is available on the webpage of the Japan Meteorological Agency (https ://www.data.jma.go.jp/ svd/eqev/data/bulle tin/hypo.html). Waveform data are available on the webpage of the National Research Institute for Earth Science and Disaster Prevention (National Research Institute for Earth Science and Disaster Resilience 2019a, b). S Hmax orientations of Townend and Zoback (2006) are extracted from the World Stress Map database (https ://www.world -stres s-map.org/) (Heidbach et al. 2016). The slip model of the Kobe earthquake is taken from Earthquake Source Model Database (https ://equak e-rc.info/srcmo d/) (Mai and Thingbaijam 2014).

Competing interests
The authors declare that there are no competing interests.
Author details 1 Geological Survey of Japan, AIST, Tsukuba Central 7, 1-1-1 Higashi, Tsukuba, Ibaraki 305-8567, Japan. 2 Earthquake Research Institute, The University of Tokyo, 1-1-1 Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan. Figure 14 shows the comparison of M w determined in the present study, and M j . We also plotted the M w values of the 1995 Kobe earthquake and the Awaji Island earthquake. The M j -M w relationship has a slope of 1/2 for M j < 3-4 and 1 for larger M j . The offset in M j relative to M w seen in larger earthquakes is typically observed for shallow inland earthquakes in Japan (Uchide and Imanishi 2018). The red curve in Fig. 14 shows the M j -M w scaling relation determined in another area in Japan (Uchide and Imanishi 2018), which has a slope similar to that of the Awaji Island case, but with a certain offset. Uchide and Imanishi (2018) suggested that the lack of high-frequency components in observed seismic data is due to the anti-aliasing filter and the anelastic attenuation effect during wave propagation. Since the sampling rate and The red curve is the M j -M w scaling determined in another area in Japan (Uchide and Imanishi 2018) the anti-aliasing filter are identical in this study and in Uchide and Imanishi (2018), the differences in the anelastic effect and high-frequency radiation from the source may play a role.

Appendix B: Slip model of the mainshock
The slip model of the mainshock was estimated using the slip inversion analysis method (Uchide and Ide 2007). We assumed a single fault plane with the same strike and dip as those estimated from on-fault aftershocks. The fault model spans 12 km in the strike direction and 16.5 km in the dip direction. We set the rupture initiation point as the relocated mainshock hypocenter: 34.426° N, 134.831° E, at the depth of 14.4 km. The spatial distribution of the fault slip-rate is represented by linear spline functions over grids on the fault plane. The nodal points of the linear spline slip-rate functions were set every 1.5 km in both the strike and dip directions and every 0.5 s in time. Since this is a reverse-faulting event, we allowed  the rake angle to vary between 45° and 135° by preparing two non-negative slip functions with 45° and 135° of rake angles for each grid. The time window (3 s) at each grid point starts at time r/v hr , where r is the distance from the rupture initiation point to the grid, and v hr is the hypothetical rupture velocity. We determined v hr = 3.0 km/s by trial and error. The unknown model parameters were peaks of the slip-rate function (total = 700). We applied a temporal smoothing constraint and searched the intensity of the smoothing constraint to minimize the Akaike's Bayesian Information Criterion (ABIC) (Akaike 1980;Yabuki and Matsu'ura 1992). We used three-component strong-motion data observed by KiK-net (Fig. 15a), which is operated by NIED (Okada et al. 2004). Original acceleration data were numerically integrated into velocity, resampled every 0.25 s, bandpass-filtered between 0.1 and 0.5 Hz, cut from 1.0 s before to 19.0 s after the S arrivals, and normalized by the maximum absolute value in each component. The lower and upper limits of the passband were restricted by source duration and the quality of Green's functions, respectively. We calculated Green's functions using the reflection and transmission method (Kennett and Kerry 1979) and wavelength integral (Bouchon 1981) based on the one-dimensional structures of seismic velocity and quality (Q) factors (Table 1). The anelastic attenuation was considered by introducing complex velocity (Takeo 1985).
The obtained slip model is illustrated in Fig. 15b. The variance reduction, (1 -v m /v d ) × 100 [%] (v m : variance of misfit; v d : variance of data), is 43.3%. Although the variance reduction is not always high, this model captures the main pulse of the observed seismograms (Fig. 15c). The slip is concentrated in the area deeper than the hypocenter. The total seismic moment ( M o ) is 3.4 × 10 17 Nm (M w 5.6). The on-fault aftershocks are located around the high slip area, whose slip is 15 cm or more. Based on this observation, we inferred that the dimension of the coseismic slip ( S ) is approximately 6 km× 6km . The coseismic stress drop ( σ ) of 4 MPa was calculated using the formula: �σ = 2.5M o /S 1.5 (Udías et al. 2014). The average rake angle of grids with a slip of 15 cm or more was 101°.
Appendix C: Stress transferred by the 1995 Kobe earthquake: the case of the 2018 M w 5.6 (M j 6.1) northern Osaka earthquake As a reference, we also calculated the CFF on the fault plane for the 2018 M w 5.6 (M j 6.1) northern Osaka earthquake, which occurred about 25 km east of the northeast end of the Kobe earthquake fault (Fig. 1). Previous studies reveal that the northern Osaka earthquake consists of two fault segments. Here, the mainshock rupture was initiated on a NNW-SSE-striking reverse fault and then propagated to an adjacent ENE-WSW-striking strikeslip fault (Hallo et al. 2019;Kato and Ueda 2019;Li et al. 2019). A large non-double-couple component of the CMT solution (Fig. 1) exhibits this complex rupture process. Here, we evaluated the CFF on the reverse fault, i.e., the first ruptured plane, whose fault parameter is (strike, dip, slip) = (345°, 50°, 85°) (Kato and Ueda 2019). Different from the Awaji Island earthquake, coseismic CFF at the hypocenter (34.8443° N, 135.6217° E, 13.0 km) was negative (− 0.01 MPa), which delays the occurrence of the earthquake. Although the postseismic CFF is positive (0.002 kPa in 23 year), the impact of the Kobe earthquake on the occurrence of the northern Osaka earthquake was negligible because of the low absolute value. The stress change around the source region of the northern Osaka earthquake, caused by the interseismic plate locking, was in the order of a few kPa/year as a mixture of the reverse and strike-slip faulting with E-W compression, which is consistent with a stress field causing the northern Osaka earthquake (Hallo et al. 2019). Here, the stress field is mixed because of the direction of minimum principal stress induced by the interseismic plate locking deviates from the vertical direction [see Fig. 5 of Saito et al. (2018)]. We inferred that the stress build-up due to the interseismic plate locking offset the coseismic stress reduction in a few years, and plays an important role in the occurrence of the northern Osaka earthquake.