Detecting the stress anomaly before occurrence the 2016 Kumamoto Earthquake, Kyushu Island, Japan

Detecting the stress anomaly before occurrence the 2016 Kumamoto Earthquake, Kyushu Island, Japan | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Detecting the stress anomaly before occurrence the 2016 Kumamoto Earthquake, Kyushu Island, Japan Yushi Nagayama, Satoshi Matsumoto, Kentaro Emoto, Takeshi Matsushima, and 1 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-6532782/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract To evaluate future seismic activity around fault zones, understanding the stress state and its buildup is crucial. There are two major factors of earthquake occurrence, the stress concentration and the fault strength weakening, but it is difficult to distinguish between the two. Therefore, when estimating the stress field around the fault zone, it is important to consider the stress concentration prior to the earthquake. In this study, we modeled the stress field incorporating pre-seismic stress concentration at the focal area of the 2016 Kumamoto earthquake (Mj 7.3), Kyushu, Japan, using focal mechanism data and slip distribution on the mainshock fault plane. The estimated stress field shows the presence of stress concentration around the earthquake fault, corresponding to approximately 29% of the co-seismic stress change before the earthquake sequence. This stress concentration may reflect the strength heterogeneity and crustal deformation in the upper and lower crust, as suggested in previous studies. In addition, we also found the stress concentration after the mainshock, and its stress concentration may contribute to aftershock and post-seismic deformation. The 2016 Kumamoto Earthquake Kyushu stress field stress concentration Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 1 Introduction Earthquakes occur when stress acting on a medium reaches its shear strength. Stress concentration and weakening of fault strength are primary factors by which the shear stress acting on the fault plane exceeds its strength. However, it is difficult to distinguish between the two. Thus, when evaluating the stress state around a fault zone, it is important to account for the presence of the stress concentration and the mechanism of its buildup. Various factors contribute to the stress loading and its concentration on the earthquake fault. For the case of stress loading on an over-riding plate in the plate subduction system, the interplate coupling between the continental and subducting oceanic plate causes stress loading and concentration. In addition, the anelastic deformation in the lower crust beneath a fault zone also causes the stress concentration in the fault zone (e.g., Iio et al., 2002 , 2004 ). Furthermore, the heterogeneity in the upper crust contributes to the stress build up and strain accumulation around the fault zone (e.g., Meneses-Gutierrez et al., 2018 ). Nevertheless, it is important to investigate the above factors that are responsible for stress build, identifying the origin of the stress field is difficult because of a lack of information about absolute stress levels. As pointed out by studies estimating absolute stress fields (e.g., Mitsuoka et al., 2020 ; Terakawa et al., 2025 ), it is required to detect spatial-temporal heterogeneous stress field with known absolute stress change such as fault slip of a large earthquake occurrence. For the 2016 Kumamoto Earthquake, the slip distribution on the mainshock fault plane (e.g., Asano and Iwata 2021 ) is available, and a large amount of seismic data can be utilized due to the high seismicity prior to the mainshock. Therefore, we focus on the focal area of the 2016 Kumamoto earthquake sequence to obtain detailed insights into the stress field around the fault. The 2016 Kumamoto earthquake sequence occurred in the Hinagu and Futagawa Fault Zones in Kyushu Island, southwest Japan (Fig. 1 ). The mainshock (Mj 7.3) occurred on April 16, 2016, following the largest foreshock (Mj 6.5) on April 14, 2016. The large fault slip and co-seismic stress change due to the mainshock were estimated from various kinds of the observation data, such as geodetic data (e.g. Fukahata and Hashimoto 2016 ) and seismic data (e.g. Asano and Iwata 2016 , 2021 ; Mitsuoka et al., 2020 ). Additionally, the presence of the post-seismic deformation was estimated from the geodetic data (e.g. Pollitz et al., 2017 ; Moore et al., 2017 ; Liu et al., 2024 ). Around the Hinagu and Futagawa Faults, the spatially heterogeneous stress field in strike-slip or normal fault type was estimated in both pre- and post-mainshock periods by some previous studies (e.g., Matsumoto et al., 2015b , 2018 ; Mitsuoka et al., 2020 ). Furthermore, in the focal area, the presence of a low b-value and its spatial variation have been estimated prior to the earthquake sequence (Nanjo et al., 2016 ). Therefore, these results might indicate stress concentration in the focal area of the 2016 Kumamoto earthquake. While some previous studies have examined the spatial-temporal changes and the absolute magnitude of the stress field in the target area (e.g. Mitsuoka et al., 2020 ; Terakawa et al., 2025 ), there is no study that has discussed the distribution and factors of the stress concentration before and after the earthquake sequence around the earthquake fault. In this study, we model the stress field in the focal area of the 2016 Kumamoto earthquake sequence by combining the regional stress field and co-seismic stress change due to the mainshock and detect the stress concentration in the fault zone prior to the mainshock occurrence. Additionally, we discuss the factor of stress concentration in the source region from the perspective of crustal strength heterogeneity. 2 Data We used the dataset of focal mechanism in the target area (longitude: 130.5° − 131.2°E, latitude: 32.5° − 33.1° N, depth: 0-20km). We analyzed data for three periods; Period I, II and III are from 1993 to the Mj6.5 largest foreshock occurrence on 14 April 2016, from the mainshock on 16 April 2016 to 19 November 2021 and from 20 November 2021 to 31 August 2023, respectively. The seismic stations were operated by Kyushu University, the National Research Institute for Earth Science and Disaster Prevention (NIED) and the Japan Meteorological Agency (JMA). The station distributions for three periods are shown in Fig. 2 . We note that high density seismic stations were installed by Kyushu University around the Hinagu Fault Zone in Period III. We determined focal mechanisms from P wave polarity data observed at 10 or more seismic stations by using the HASH algorithm (Hardebeck and Shearer 2002 ). The Institute of Seismology and Volcanology (SEVO), Kyushu University picked these P wave polarity data manually in Period I and II and automatically in Period III. Magnitude ranges for normalized stress tensor estimated using Matsumoto ( 2016 ) are 1.5–4.0 for Period I and II and 0.5–4.0 for Period III, respectively. Finally, we obtained 2681, 6251 and 591 good quality focal mechanism solutions for Period I, II and III, respectively (Fig. 3 ). In order to obtain co-seismic stress change, we adopted co-seismic slip distribution estimated by Asano and Iwata ( 2021 ). They estimated the slip distribution of the mainshock fault using the strong ground motion waveform data around the earthquake fault. 3 Method We modeled the stress field around the focal area of the 2016 Kumamoto earthquake sequence using the method proposed by Matsumoto et al. ( 2015a ). We assumed that the deviatoric stress tensor within the target region can be represented as sum of a uniform regional stress field and co-seismic stress change resulting from the mainshock. In order to express stress field taking into account the stress concentration around the focal area prior to the mainshock occurrence, we write the deviatoric stress tensor at position of the \(\:k\) -th eartquake \(\:{\varvec{x}}_{\varvec{k}}\) in the target region as $$\:\begin{array}{c}{S}_{ij}^{b}\left({\varvec{x}}_{\varvec{k}}^{\varvec{b}}\right)={S}_{ij}^{r}-\left(1-\alpha\:\right)\varDelta\:{S}_{ij}\left({\varvec{x}}_{\varvec{k}}^{\varvec{b}}\right)\#\left(1\right)\end{array}$$ $$\:\begin{array}{c}{S}_{ij}^{a}\left({\varvec{x}}_{\varvec{k}}^{\varvec{a}}\right)={S}_{ij}^{r}+\alpha\:\varDelta\:{S}_{ij}\left({\varvec{x}}_{\varvec{k}}^{\varvec{a}}\right)\#\left(2\right)\end{array}$$ $$\:\begin{array}{c}{S}_{ij}^{r}=\xi\:{s}_{ij}\:\#\left(3\right)\end{array}$$ where \(\:{S}_{ij}\left({\varvec{x}}_{\varvec{k}}\right)\) is the \(\:i\) and \(\:j\) components of the deviatoric stress tensor at \(\:\:{\varvec{x}}_{\varvec{k}}\) , \(\:{S}_{ij}^{r}\) is the uniform regional deviatoric stress tensor, \(\:{\Delta\:}{S}_{ij}\left({\varvec{x}}_{\varvec{k}}\right)\) is the stress tensor equivalent to the co-seismic stress change at \(\:{\varvec{x}}_{\varvec{k}}\) , and \(\:{s}_{ij}\) is the normalized regional deviatoric stress tensor. The superscripts \(\:b\) and \(\:a\) represent the periods before the largest foreshock and after the mainshock, respectively. \(\:\xi\:\) determines the magnitude of \(\:{S}^{r}\) , while \(\:\alpha\:\) represents the degree of the stress concentration around the earthquake fault prior to the earthquake sequence. When \(\:\alpha\:\) is close to 0, the stress concentration prior to the earthquake is observed around the fault. Conversely, for \(\:\alpha\:\) close to 1, the stress field before earthquake was relatively uniform compared to after the earthquake. Among these parameters, \(\:{s}_{ij}\) and \(\:{\Delta\:}{S}_{ij}\) can be estimated or calculated from the observed data and fault model. Therefore, we determined the unknown parameters \(\:\alpha\:\) and \(\:\xi\:\) by a grid search analysis. 3.1 Estimation of regional stress field We estimated the regional deviatoric stress field in the target region using the method proposed by Matsumoto ( 2016 ). In this method, the deviatoric stress tensor is estimated from the sum of seismic moment tensors of earthquakes in the target region based on the plasticity theory. The deviatoric stress tensor can be obtained by the following formula: $$\:\begin{array}{c}{s}_{ij}\propto\:\sum\:_{k=1}^{K}{M}_{ij}^{k}\:\#\left(4\right)\end{array}$$ where \(\:{s}_{ij}\) is the normalized deviatoric stress tensor, \(\:{M}_{ij}^{k}\) is the seismic moment tensor released by \(\:k\) -th earthquake, and \(\:K\) is the total number of earthquakes in the target area. \(\:{M}_{ij}^{k}\) is obtained as the product of \(\:{s}_{ij}\) and the seismic moment \(\:{M}_{O}\) . \(\:{M}_{O}\) is determined using an empirical relationship between the seismic moment and the magnitude of each earthquake \(\:M\) as (Matsumoto et al., 2016 ) $$\:\begin{array}{c}\text{log}{M}_{O}=1.151M+10.548\:\#\left(5\right)\end{array}$$ We estimated the regional deviatoric stress using data from both before (Period I) and after (Period II and III) the mainshock. However, there is difference in the number of events among the periods, so we estimated the normalized regional stress tensor by averaging the normalized stress tensors for Period I and Period II and III and re-normalizing. As a result of estimating regional stress field, the \(\:{\sigma\:}_{3}\) axis, the minimum principal compression stress axis of the stress field, is horizontal and oriented approximately N-S (Fig. 4 ). This result is consistent with the previous studies on the stress field around the Futagawa and Hinagu Fault Zones in both pre- and post-mainshock periods (e.g. Matsumoto et al., 2015b ; Matsumoto et al., 2018 ; Mitsuoka et al., 2020 ). In addition, the regional deviatoric stress tensor we estimated had a high stress ratio ( \(\:({\sigma\:}_{2}-{\sigma\:}_{3})/({\sigma\:}_{1}-{\sigma\:}_{3})\approx\:0.86\) ,where \(\:{\sigma\:}_{1},\:{\sigma\:}_{2}\) and \(\:{\sigma\:}_{3}\) are the maximum, intermediate and minimum principal compression stresses, respectively). This high stress ratio might be attributed to the events with a normal faulting in the target region, and it consists with the results of Mitsuoka et al. ( 2020 ) estimating high stress ratio around the earthquake fault. 3.2 Calculation of stress tensor due to co-seismic fault slip We calculated the stress change tensor \(\:{\Delta\:}{S}_{ij}\left({\varvec{x}}_{\varvec{k}}\right)\) due to the co-seismic fault slip at each earthquake’s hypocenter \(\:{\varvec{x}}_{\varvec{k}}\) using the formula developed by Okada ( 1992 ). We adopted the slip distribution of the mainshock obtained by Asano and Iwata ( 2021 ) from strong ground motion data. We assumed that the medium is a homogeneous half space with a rigidity of 33 GPa. Figure 5 shows the differential stress ( \(\:{\sigma\:}_{1}-{\sigma\:}_{3}\) ) of the \(\:{\Delta\:}{S}_{ij}\) around the focal area at a depth of 10 km. The differential stress exceeds 0.1 MPa in most of the target area and is significantly larger (> 10 MPa) near the mainshock fault, particularly around the Futagawa Fault, compared to the surrounding areas located farther from the earthquake fault. 3.3 Grid search analysis and definition of RMS residuals value We determined two unknown parameters \(\:\alpha\:\) and \(\:\xi\:\) in equations (1) and (2) using the grid search analysis. To determine the optimal pair of \(\:\alpha\:\) and \(\:\xi\:\) , we first calculated the residuals as the differences between the observed unit slip vector \(\:\widehat{\varvec{s}}=({\widehat{s}}_{1},\:{\widehat{s}}_{2},\:{\widehat{s}}_{3})\) from the focal mechanism solutions at located \(\:{\varvec{x}}_{\varvec{k}}\) and the unit vector \(\:\widehat{\varvec{\tau\:}}=({\widehat{\tau\:}}_{1},\:{\widehat{\tau\:}}_{2},\:{\widehat{\tau\:}}_{3})\) of the maximum shear stress direction on the fault plane. The unit vector \(\:\widehat{\varvec{\tau\:}}\) was calculated from the deviatoric stress tensor \(\:S\) at the same position \(\:{\varvec{x}}_{\varvec{k}}\) and matrix \(\:A\) formulated from the fault normal vector components through their combinations, represented by Michael ( 1984 ). Then, we calculated the RMS (Root Mean Square) residuals value between \(\:\widehat{\varvec{s}}\) and \(\:\widehat{\varvec{\tau\:}}\) from all focal mechanisms. Because there is a large difference in the number of the events between pre- and post-earthquake periods, to avoid bias in the results, we calculated the RMS residual value by weighting it with the inverse of the number of focal mechanisms each period. The \(\:RMS\) is defined as the following formula: $$\:\begin{array}{c}RMS=\sqrt{\frac{\left[\frac{1}{3{K}_{b}}RS{S}_{b}+\frac{1}{3{K}_{a}}RS{S}_{a}\right]}{2}}\:\#\left(6\right)\end{array}$$ $$\:\begin{array}{c}RS{S}_{b}=\sum\:_{k=1}^{{K}_{b}}\sum\:_{i=1}^{3}{\left({\widehat{\varvec{s}}}_{i,k}^{b}-{\widehat{\varvec{\tau\:}}}_{i,k}^{b}\right)}^{2}\:\#\left(7\right)\end{array}$$ $$\:\begin{array}{c}RS{S}_{a}=\sum\:_{k=1}^{{K}_{a}}\sum\:_{i=1}^{3}{\left({\widehat{\varvec{s}}}_{i,k}^{a}-{\widehat{\varvec{\tau\:}}}_{i,k}^{a}\right)}^{2}\:\#\left(8\right)\end{array}$$ where \(\:K\) is the total number of the focal mechanism data, and \(\:{\widehat{\varvec{s}}}_{i,k}\) and \(\:{\widehat{\varvec{\tau\:}}}_{i,k}\) are \(\:i\) - component of the unit vectors \(\:\widehat{\varvec{s}}\) and \(\:\widehat{\varvec{\tau\:}}\) for the \(\:k\) -th focal mechanism, respectively, \(\:RSS\) is the residual sum of squares between the unit vector \(\:\widehat{\varvec{s}}\) and \(\:\widehat{\varvec{\tau\:}}\) and the subscript \(\:b\) and \(\:a\) correspond to before the largest foreshock (Period I) and after the mainshock (Period II and III), respectively. In this approach, a focal mechanism generally has two nodal planes. Therefore, to determine the appropriate fault plane, we calculated \(\:\widehat{\varvec{\tau\:}}\) at both planes and identified the one with the smaller \(\:RSS\) value as the fault plane. The grid search ranges are 1 kPa to 100 MPa and 0 to 1 for \(\:\xi\:\) and \(\:\alpha\:\) , respectively. We determined the optimal pair of parameters by minimizing the \(\:RMS\) . In addition, in Periods II and III, we used only events at least 1 km away from the mainshock fault plane. 4 Result 4.1 Result of the grid search and the optimum stress field model As a result of the grid search, the optimal pair of the parameters was determined to be \(\:\alpha\:=0.71\) and \(\:\xi\:=5.62\) MPa. This obtained \(\:\alpha\:\) indicates that 29% of the co-seismic stress change was concentrated around the earthquake fault prior to the earthquake sequence. Furthermore, we estimated the confidence intervals for \(\:\alpha\:\) and \(\:\xi\:\) using the bootstrap resampling. In this approach, focal mechanism solutions for each period were resampled 1000 times, allowing for duplication, and the 95% confidence interval was calculated. During this process, the regional stress field was also re-estimated for each resampled dataset (Fig. 4 b). As a result, the confidence intervals of \(\:\alpha\:\) and \(\:\xi\:\) were estimated to be 0.55–1.00 and 3.16–7.94 MPa, respectively (Fig. 6 ). The confidence interval of \(\:\alpha\:\) suggests that a stress field with a stress concentration equivalent to 0–45% of the co-seismic stress change before the earthquake sequence better explains the observed data. We calculated the deviatoric stress field in the target area using the previously determined optimal parameters ( \(\:\alpha\:=0.71,\:\xi\:=5.62\) MPa) and Equations (1) and (2). Figure 7 shows the distribution of the differential stress at a depth of 12km before and after the mainshock, respectively. Before the earthquake sequence, there were two dominant high differential stress regions at the southern part of the Hinagu Fault and the central part of the Futagawa Fault at a depth of 12 km. In contrast, at the ends of the source fault plane (i.e. the southern edge of the Hinagu Fault, the northeastern of the Futagawa Fault and the junction of the Hinagu and Futagawa Faults), the differential stress was lower than that in the surrounding areas. On the other hand, after the mainshock, areas with high differential stress are located near the Hinagu and Futagawa Fault planes and the extension zones of both Faults, while differential stress is low in the central part of the Hinagu Fault and around the Futagawa Fault. 4.2 The relationship between stress concentration and hypocenters of each earthquake To verify the consistency of the stress concentration and observed seismic activity, we compared the distribution of stress concentration and hypocenters of earthquakes (Fig. 7 ). These hypocenters were determined using data manually picked by SEVO. As shown in Fig. 7 , the low seismicity region corresponds to relatively low differential stress, while the high seismicity regions are associated with high differential stress during both the pre- and post-earthquake periods. In addition, by comparing the pre- (Period I) and post-earthquake sequence (Periods II and III), we found that stress built up, particularly at the southwestern end of the rapture zone of the Hinagu Fault, the western side of the Hinagu Fault, and the northern end of the Futagawa Fault. These areas were characterized by low seismicity in Period I and high seismicity in Periods II and III. Therefore, both spatially and temporally, there was a strong correspondence between the distribution of the differential stress and the seismic activity around the earthquake fault. However, the correspondence was not observed in some regions around the fault. For example, in Period I, although the differential stress around the central part of Futagawa Fault and the southern part of the Hinagu Fault is significantly higher than in the surrounding areas, seismic activity is not particularly high. Furthermore, to the north of the central part of the Futagawa Fault and the southeast part of the Hinagu Fault, an aftershock cluster with high activity exists, while the differential stress is low in Periods II and III. On the other hand, in the central part of the Futagawa Fault, a few aftershocks were observed, despite high differential stress in the stress field in Periods II and III. These inconsistencies between the distribution of hypocenters and differential stress is further addressed in the Discussion section. 4.3 Misfits of the stress field model To evaluate the optimal stress field model, we calculated the misfit angles between the unit slip vector \(\:\widehat{\varvec{s}}\) of focal mechanisms and the expected unit vector \(\:\widehat{\varvec{\tau\:}}\) oriented to the direction of the maximum shear stress on the fault plane from the stress tensor. We defined the misfit angle \(\:{\theta\:}_{mis}\) as $$\:\begin{array}{c}{\theta\:}_{mis}=\text{arccos}\left(\widehat{\varvec{s}}\bullet\:\widehat{\varvec{\tau\:}}\right)\:\#\left(9\right)\end{array}$$ We calculated \(\:\widehat{\varvec{\tau\:}}\) and \(\:{\theta\:}_{mis}\) at the two nodal planes of the focal mechanism and selected the plane with smaller \(\:{\theta\:}_{mis}\) as the fault plane, in the same way as the previously defined RMS value. As a result of misfit angle \(\:{\theta\:}_{mis}\) , some areas had large \(\:{\theta\:}_{mis}\) values (Fig. 8 ). In the stress field prior to the earthquake sequence, the areas with large \(\:{\theta\:}_{mis}\) were located to the north, away from the Futagawa Fault, and in the shallow region (0-10km depth) at the northeastern end of the Futagawa Fault. On the other hand, in the stress field in Periods II and III, some regions with particularly large \(\:{\theta\:}_{mis}\) located in the southeast of the Hinagu Fault source area and near the junction of the Hinagu and Futagawa Faults. The potential causes of the large \(\:{\theta\:}_{mis}\) value described above are also examined in detail in the Discussion section. 5 Discussion As described in Result section, there are some inconsistencies between the distribution of hypocenters and differential stress (Fig. 7 ). The inconsistency, in Period I, characterized by low seismic activity despite high differential stress in the central part of the Futagawa Fault and the southern part of the Hinagu Fault, may be attributed to higher crustal strength in these regions relative to its surrounding areas. This suggests that greater stress could be sustained without being released by earthquakes in these regions. In addition, other inconsistencies also observed in Periods II and III, namely, high seismic activity exists while the differential stress is low around the north of the central part of the Futagawa Fault and the southeast part of the Hinagu Fault and low seismic activity exist while the difference stress is particularly high in the central part of the Futagawa Fault. A possible explanation for this exception is the difference in geometry between the assumed and the actual fault planes. In this study, we adopted the flat fault plane model proposed by Asano and Iwata ( 2021 ). However, the actual fault plane may not be perfectly flat. In addition, the after-slip pointed out by some previous studies (e.g., Pollitz et al., 2017 ; Moore et al., 2017 ; Liu et al., 2024 ) may have affected seismic activity by releasing stress. Alternatively, in this region, it is also possible that strain accumulated by the mainshock remains without releasing stress. Moreover, as mentioned in Result section, there are some regions which have large \(\:{\theta\:}_{mis}\) value in the target area. In Period I, large \(\:{\theta\:}_{mis}\) values were located to the north away from the Futagawa Fault, and in the shallow region at the northeastern end of the Futagawa Fault (Fig. 8 a and 8 b). The former region is consisting of multiple clusters of earthquakes. Since these earthquakes are away from the mainshock fault, they may have been triggered by factors not considered in this study’s stress field model, such as local stress heterogeneity, resulting in large \(\:{\theta\:}_{mis}\) . The large misfit in the latter region is located near Mt. Aso, which also could be attributed to the heterogeneity of the stress field due to the strong heterogeneous structure (e.g., Savage et al., 2016 ). On the other hand, in Periods II and III, some regions with particularly large \(\:{\theta\:}_{mis}\) located in the southeast of the Hinagu Fault source area and near the junction of the Hinagu and Futagawa Faults (Fig. 8 c and 8 d). These regions might be affected by stress heterogeneity and/or the fluid intrusion that existed before the earthquake sequence, as reported by Terakawa et al. ( 2025 ). They estimated the absolute background stress field in the focal area of the 2016 Kumamoto earthquake from the perspective of the effective friction coefficient and strain energy by combining the background stress field and co-seismic stress changes due to both the largest foreshock and the mainshock. They suggested that there were some earthquakes occurred due to factors not accounted for in their stress field model. Even though we took into account the co-seismic stress change due to the largest foreshock in addition to the mainshock, the average values of the \(\:{\theta\:}_{mis}\) around the Hinagu Fault decreased only slightly, and large \(\:{\theta\:}_{mis}\) values remained in this area compared to the surrounding regions (Figure S4). Therefore, the large \(\:{\theta\:}_{mis}\) in these regions might have been caused by factors other than the regional stress field or co-seismic stress changes during the largest foreshock and the mainshock. Some previous studies estimated the stress field in both pre- and post-mainshock periods around the Hinagu and Futagawa Fault zones. Matsumoto et al. ( 2015b , 2018 ) reported the stress field around the Hinagu and Futagawa Fault zones and their surrounding area during the pre-mainshock period. They estimated the spatially heterogeneous stress field in the target area. In particular, Matsumoto et al. ( 2018 ) reported the presence of the lateral heterogeneous and depth dependent stress field around the source fault, which implies that stress concentration factor \(\:\alpha\:\) is not completely 1. On the other hand, the spatial heterogeneity of the stress field around the focal area after the mainshock was also estimated in detail by some previous studies (e.g., Yoshida et al., 2016 ; Mitsuoka et al., 2020 ). Therefore, the stress field during the post-earthquake period also cannot be expressed by only a uniform regional stress field, namely the model with \(\:\alpha\:=0\) is contradictory to these findings. Following them, the optimal model with \(\:\alpha\:=0.71\) , which has the stress concentration equivalent to 29% of the co-seismic stress change, is qualitatively consistent with the stress field estimated in previous studies. However, there are differences in the absolute magnitude of the deviatoric stress tensor compared to previous studies (Figure S6). The magnitude of the regional deviatoric stress tensor representing the regional stress in the target area was slightly smaller than that estimated by Mitsuoka et al. ( 2020 ) for the pre-earthquake period (= 7.8 MPa) using focal mechanism data, but there was no significant difference exceeding the 95% confidence interval. However, Terakawa et al. ( 2025 ), based on the the effective friction coefficient and strain energy by combining the background stress field and co-seismic stress changes due to both the largest foreshock and the mainshock, estimated the absolute magnitudes of the deviatoric stress tensors are 37–65 MPa and 39–70 MPa around the Hinagu Fault and Futagawa Fault depth at 10 km, respectively. The difference in the differential stresses between this study and Terakawa et al. ( 2025 ) might be attributed to the evaluation of the stress field in the vicinity of the earthquake fault. We have estimated the absolute stress by excluding the focal mechanism data in this study located near the mainshock fault plane, within 1.0 km, during Periods II and III. When we use the events whose distances from the fault plane of the mainshock are 0.5 km or larger in the grid search analysis, \(\:\xi\:\) increases and \(\:\alpha\:\) decreases (Figure S7). This leads to a larger magnitude of the deviatoric stress tensor before the mainshock. However, in this case, the RMS residual and the confidence interval are larger than those in the original case. This implies the possibility that localized high stress existed close to the earthquake fault and radiated large seismic energy, which was not captured in this analysis. As mentioned in the introduction, one of the factors that can cause stress loading on the fault is the anelastic deformation in the lower crust beneath the fault zone. In Kyushu Island, the distribution of anelastic deformation rate in the lower crust prior to the Kumamoto earthquake sequence was estimated by Yuasa and Matsumoto ( 2023 ) based on GNSS observation data. They estimated a higher anelastic strain rate (> 0.4 µ/yr) in the lower crust from the BSG to the focal area of the Kumamoto Earthquake than its surrounding areas. Because they used large blocks (20 km × 20 km) to estimate the distribution of anelastic strain rate, it is difficult to examine the detailed distribution of anelastic strain rate around the earthquake fault. However, near the central part of the Futagawa Fault, where the stress concentration was found in this study before the mainshock, is located just above the block with the highest anelastic strain rate (~ 0.6 µ/yr) in the focal area. In addition, the stress change rate in the upper crust (at a depth of 7 km) was caused by the anelastic deformation in the lower crust, which was > 8 kPa/yr around the Futagawa Fault. Therefore, the stress concentration, especially large in the Futagawa Fault zone, during the pre-earthquake period might have been caused by the higher anelastic deformation in the lower crust beneath the source fault than it was in the rest area. In addition, the strength heterogeneity in the upper and lower crust is also the factor causing the localized stress concentration. Some previous studies have revealed the strength heterogeneity in the crust and proposed that it might have caused strain accumulation and generation of large earthquakes (Usui et al., 2024 ). In the focal area of the 2016 Kumamoto earthquake sequence, various crustal structure surveys, such as velocity structure and electrical resistivity structure, have been conducted, and the heterogeneity in crustal structure was pointed out (e.g., Shito et al., 2017 ; Aoyagi et al., 2020 ; Aizawa et al., 2021 ). In particular, the electrical resistivity structure, sensitive to the presence of fluid, provides insights into the distribution of the mechanical strength in the crust. Aizawa et al. ( 2021 ) revealed the low-resistivity zones in the deep part of the Hinagu Fault. They interpreted these low-resistivity zones as weak structures; the deep ones associated with high-temperature magmatic fluid and the shallow ones enriched with clay. In this study, the stress concentration regions prior to the mainshock were estimated to be in the southeastern part of the Hinagu and the central part of the Futagawa Faults, which are adjacent to the low-resistivity zone identified by Aizawa et al. ( 2021 ). Therefore, around the fault zone, such fluid-rich weak zones in the crust might have contributed to the strain accumulation in these regions and stress concentration around them. Based on the above, it can be inferred that the stress concentration prior to the mainshock around the focal area of the 2016 Kumamoto Earthquake was caused by strain accumulation in the weak zones of both the lower and upper crust (Fig. 9 ). 6 Conclusion We estimated the stress field model in the focal area of the 2016 Kumamoto earthquake, including the stress concentration before the earthquake sequence using high quality focal mechanism data and the slip distribution data of the mainshock fault plane. We applied the method proposed by Matsumoto et al. ( 2015a ) to the target region and determined the magnitudes of the stress concentration and regional stress field by the grid search analysis. As a result, the optimal model had the stress concentration equivalent to approximately 29% of the co-seismic stress change around the fault zone before the mainshock. This stress concentration distributed around the weak strength structure in the crust indicated by previous studies. This result indicates that the deformation could occurs in the low-viscosity structures in the lower crust and low-elasticity structures in the upper crust, and the stress concentrated around there and loaded on the earthquake fault. In addition, there were also the stress concentration around the fault zone after mainshock, and it might have driven aftershock activity and post-seismic deformation at the extended part of the earthquake fault. We revealed the relationship between the distribution of the stress concentration and the strength heterogeneity in the crust around the fault zone of the 2016 Kumamoto earthquake sequence. Therefore, in the future studies, estimating the heterogeneity of the strength structure and the stress field around the fault zone with high resolution may provide the presence of the stress concentration and the potential of the future seismic activity. Abbreviations BSG: Beppu-Shimabara graben MLT: the median tectonic kine GNSS: Global Navigation Satellite System Declarations Ethics approval and consent to participate Not applicable Consent for publication Not applicable List of abbreviations BSG: Beppu-Shimabara graben MLT: the median tectonic kine GNSS: Global Navigation Satellite System Availability of data and materials Dataset of focal mechanism data is available upon reasonable request. If any readers might be interested in accessing the data, please contact Satoshi Matsumoto ( [email protected] ). Competing interests The authors declare that they have no competing interests. Funding The Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, under its The Third Earthquake and Volcano Hazards Observation and Research Program (Earthquake and Volcano Hazard Reduction Research). Authors' contributions YN analyzed the data and designed the research with help from SM. All authors contributed to the discussion of the research and read the final manuscript. Acknowledgements We used seismic data from the JMA and Hi-net (NIED). This study was supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, under its The Third Earthquake and Volcano Hazards Observation and Research Program (Earthquake and Volcano Hazard Reduction Research). Authors' information Department of Earth and Planetary Sciences, Graduate School of Science, Kyushu University, Fukuoka 819-0395, Japan Yushi Nagayama & Kono Taiki Institute of Seismology and Volcanology, Faculty of Science, Kyushu University, Kyushu University, Shimabara 855-0843, Japan Satoshi Matsumoto, Takeshi Matsushima & Kentaro Emoto Endnotes References Aizawa K, Takakura S, Asaue H, Koike K, Yoshimura R, Yamazaki K, Komatsu S, Utsugi M, Inoue H, Tsukamoto K, Uyeshima M, Koyama T, Kanda W, Yoshinaga T, Matsushima N, Uchida K, Tukashima Y, Matsushima T, Ishigara H, Muramatsu D, Teguri Y, Shito A, Mtsumoto S, Shimizu H (2021) Electrical conductive fluid-rich zones and their influence on the earthquake initiation, growth, and arrest processes: observations from the 2016 Kumamoto earthquake sequence, Kyushu Island, Japan. Earth, Planets and Space 73:1-12. https://doi.org/10.1186/s40623-020-01340-w Aoyagi Y, Kimura H, Mizoguchi K (2020) Seismic velocity structure at the southern termination of the 2016 Kumamoto Earthquake rupture, Japan. Earth, Planets and Space 72:1-14. https://doi.org/10.1186/s40623-020-01276-1 Asano K, Iwata T (2016) Source rupture processes of the foreshock and mainshock in the 2016 Kumamoto earthquake sequence estimated from the kinematic waveform inversion of strong motion data. Earth Planets and Space 68:147. https://doi.org/10.1186/s40623-016-0519-9 Asano K, Iwata T (2021) Revisiting the source rupture process of the mainshock of the 2016 Kumamoto Earthquake and implications for the generation of near‐fault ground motions and forward‐directivity pulse. Bulletin of the Seismological Society of America 111.5: 2426-2440. https://doi.org/10.1785/0120210047 Fukahata Y, Hashimoto M (2016) Simultaneous estimation of the dip angles and slip distribution on the faults of the 2016 Kumamoto earthquake through a weak nonlinear inversion of InSAR data. Earth, Planets and Space 68: 1-10. https://doi.org/10.1186/s40623-016-0580-4 Hardebeck, J L, Shearer P M (2002) A new method for determining first-motion focal mechanisms. Bulletin of the Seismological Society of America 92.6: 2264-2276. https://doi.org/10.1785/0120010200 Iio Y, Sagiya T, Kobayashi Y, Shiozaki I (2002) Water-weakened lower crust and its role in the concentrated deformation in the Japanese Islands. Earth and Planetary Science Letters 203.1:245-253. https://doi.org/10.1016/S0012-821X(02)00879-8 Iio Y, Sagiya T, Kobayashi Y (2004) Origin of the concentrated deformation zone in the Japanese Islands and stress accumulation process of intraplate earthquakes. Earth, Planets and Space 56.8: 831-842. https://doi.org/10.1186/BF03353090 Liu Y, Cui X, Hu Y, Zhang J, Chen Y (2024) Integrated investigation on heterogeneous lower crust rheology in Kyushu and afterslip behavior following the 2016 Mw7. 1 Kumamoto earthquake. Geophysical Research Letters 51.7: e2023GL107606. https://doi.org/10.1029/2023GL107606 Matsumoto S, Katao H, Iio Y (2015a) Determining changes in the state of stress associated with an earthquake via combined focal mechanism and moment tensor analysis: Application to the 2013 Awaji Island earthquake, Japan. Tectonophysics 649: 58-67. https://doi.org/10.1016/j.tecto.2015.02.023 Matsumoto S, Nakao S, Ohkura T, Miyazaki M, Shimizu H, Abe Y, Inoue H, Nakamoto M, Yoshikawa S, Yamashita Y (2015b) Spatial heterogeneities in tectonic stress in Kyushu, Japan and their relation to a major shear zone. Earth, Planets and Space 67: 1-9. https://doi.org/10.1186/s40623-015-0342-8 Matsumoto S (2016) Method for estimating the stress field from seismic moment tensor data based on the flow rule in plasticity theory. Geophysical Research Letters 43.17: 8928-8935. https://doi.org/10.1002/2016GL070129 Matsumoto S, Nishimura T, Ohkura T (2016) Inelastic strain rate in the seismogenic layer of Kyushu Island, Japan. Earth, Planets and Space 68: 1-7. https://doi.org/10.1186/s40623-016-0584-0 Matsumoto S, Yamashita Y, Nakamoto M, Miyazaki M, Sakai S, Iio Y, Shimizu H, Goto K, Okada T, Ohzono M, Terakawa T, Kosuga M, Yoshimi M, Asano Y (2018) Prestate of stress and fault behavior during the 2016 Kumamoto earthquake (M7. 3). Geophysical Research Letters 45.2: 637-645. https://doi.org/10.1002/2017GL075725 Meneses‐Gutierrez A, Sagiya T, and Sekine S (2018) Crustal deformation process in the Mid‐Niigata region of the Niigata‐Kobe tectonic zone as observed by dense GPS network before, during, and after the Tohoku‐Oki earthquake. Journal of Geophysical Research: Solid Earth 123.7: 6072-6085. https://doi.org/10.1029/2018JB015567 Michael A J (1984) Determination of stress from slip data: faults and folds. Journal of Geophysical Research: Solid Earth 89.B13: 11517-11526. https://doi.org/10.1029/JB089iB13p11517 Mitsuoka A, Shito A, Matsumoto S, Yamashita Y, Nakamoto M, Sakai S, Iio Y, Shimizu H, Goto K, Okada T (2020) Spatiotemporal change in the stress state around the hypocentral area of the 2016 Kumamoto earthquake sequence. Journal of Geophysical Research: Solid Earth 125.9: e2019JB018515. https://doi.org/10.1029/2019JB018515 Moore J DP, Yu H, Tahg C, Wang T, Barbot S, Peng D, Masuti S, Dauwels J, Hsu Y, Lambert V, Nanjundiah P, Wei S, Lindsey E, Feng L, Shibazaki B (2017) Imaging the distribution of transient viscosity after the 2016 M w 7.1 Kumamoto earthquake. Science 356.6334: 163-167. https://doi.org/10.1126/science.aal3422 Nanjo K Z, Izutsu J, Orihara Y, Furuse N, Togo S, Nitta H, Okada T, Tanaka R, Kamogawa M, Nagao T (2016) Seismicity prior to the 2016 Kumamoto earthquakes. Earth, Planets and Space 68: 1-10. https://doi.org/10.1186/s40623-016-0558-2 Okada Y (1992) Internal deformation due to shear and tensile faults in a half-space. Bulletin of the seismological society of America 82.2: 1018-1040. https://doi.org/10.1785/BSSA0820021018 Pollitz F F, Kobayashi T, Yarai H, Shibazaki B, Matsumoto T (2017) Viscoelastic lower crust and mantle relaxation following the 14–16 April 2016 Kumamoto, Japan, earthquake sequence. Geophysical Research Letters 44.17: 8795-8803. https://doi.org/10.1002/2017GL074783 Savage M. K, Aoki Y, Unglert K, Ohkura T, Umakoshi K, Shimizu H, Iguchi M, Tameguri T, Ohminato T, Mori J (2016) Stress, strain rate and anisotropy in Kyushu, Japan. Earth and Planetary Science Letters 439: 129-142. https://doi.org/10.1016/j.epsl.2016.01.005 Shito A, Matsumoto S, Shimizu H, Ohkura T, Takahashi H, Sakai S, Okada T, Miyazaki H, Kosuga M, Maeda Y, Yoshimi M, Asano Y, Okubo M (2017) Seismic velocity structure in the source region of the 2016 Kumamoto earthquake sequence, Japan. Geophysical Research Letters 44.15 (2017): 7766-7772. https://doi.org/10.1002/2017GL074593 Shimizu H, Iio Y, Sakai S, Okada T, Takahashi H, Watanabe T, et al. (2016) Urgent joint seismic observation of the 2016 Kumamoto earthquake—Seismic activities and their background. Paper presented at the Japan Geoscience Union 2016 Meeting, Japan. 22-26 May 2016 https://confit.atlas.jp/guide/event/jpgu2016/subject/MIS34-02/detail Terakawa T, Asano K, Urata Y (2025) Constraint on the background stress in the source region of the 2016 Kumamoto earthquake sequence based on temporal changes in elastic strain energies and coseismic stress rotation. Geophysical Journal International 240.1: 174-188. https://doi.org/10.1093/gji/ggae376 Usui Y, Uyeshima M, Hase H, Ichihara H, Aizawa K, Koyama T, Sakanaka S, Ogawa T, Yamaya Y, Nishitani T, Asamori K, Ogawa Y, Yoshimura R, Takakura S, Mishina M Morita Y (2024) Three‐dimensional electrical resistivity structure beneath a strain concentration area in the back‐arc side of the northeastern Japan arc. Journal of Geophysical Research: Solid Earth 129.5: e2023JB028522. https://doi.org/10.1029/2023JB028522 Yoshida K, Hasegawa A, Saito T, Asano Y, Tanaka S, Sawazaki K, Urata Y, Fukuyama E (2016) Stress rotations due to the M6. 5 foreshock and M7. 3 main shock in the 2016 Kumamoto, SW Japan, earthquake sequence. Geophysical Research Letters 43.19: 10-097. https://doi.org/10.1002/2016GL070581 Yuasa Y, Matsumoto S (2023) Anelastic deformation in the lower crust and its implications for tectonic dynamics in Kyushu, Southwest Japan. Tectonophysics 846 : 229674. https://doi.org/10.1016/j.tecto.2022.229674 Supplementary Files AdditionalfileNagayamaetal.pdf GraphicalabstractNagayamaetal.png Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-6532782","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":456585160,"identity":"460fdc25-edb1-48aa-81ee-05b9f3a5d7e1","order_by":0,"name":"Yushi Nagayama","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABU0lEQVRIie2RMUsDMRTHEw/ulkgdAwf2K7wiWIvafJUcgXY5sCB06lAQMpW6pqjfoSIIblcO2qW0uB1kaRGcbhAKRbAUL6VCr9bBTfB+Q/Ly4Mf/JUEoI+OvYyE4WxW4mSzw1V6dflYqv1ZQuEPZougMnmek8ZhH1D+a1WpjlqMiwB0ZsiLlNkWNMrJuUjGlln/pkr4uNKl/7CrQXkdVOL6XofekjNIXCN8GmwoEPneJrbFRLAKaw3AIeCpDDtHFnCI7QFjxlDKOxQdZatY0gxEYMbZWGEQmZfldiap9d19qL1HAJRDgrtMCMxjuGgXLHUpsl+7aWkjyWk8U4amB5D01qnrd4cQ68dqCbN9lXH2J4rk+v3bEw4wsyix3ZfWmrfopgwHH0du8fFhIvxhCBPZIstkHG+kBttdV0iQFtfUzzgS/mz2XSl9s1HmKMjIyMv41n9qIfLgl/BQuAAAAAElFTkSuQmCC","orcid":"","institution":"Department of Earth and Planetary Sciences, Graduate School of Science, Kyushu University","correspondingAuthor":true,"prefix":"","firstName":"Yushi","middleName":"","lastName":"Nagayama","suffix":""},{"id":456585161,"identity":"a301043b-579a-4d3b-b91d-561667b9adbb","order_by":1,"name":"Satoshi Matsumoto","email":"","orcid":"","institution":"Institute of Seismology and Volcanology, Faculty of Science, Kyushu University","correspondingAuthor":false,"prefix":"","firstName":"Satoshi","middleName":"","lastName":"Matsumoto","suffix":""},{"id":456585162,"identity":"5782c865-46d4-4312-a649-304bb9a57bd2","order_by":2,"name":"Kentaro Emoto","email":"","orcid":"","institution":"Institute of Seismology and Volcanology, Faculty of Science, Kyushu University","correspondingAuthor":false,"prefix":"","firstName":"Kentaro","middleName":"","lastName":"Emoto","suffix":""},{"id":456585163,"identity":"5417d2f9-d50c-4f2c-a4af-1e7af8fce326","order_by":3,"name":"Takeshi Matsushima","email":"","orcid":"","institution":"Institute of Seismology and Volcanology, Faculty of Science, Kyushu University","correspondingAuthor":false,"prefix":"","firstName":"Takeshi","middleName":"","lastName":"Matsushima","suffix":""},{"id":456585164,"identity":"5b421617-77c5-4d1b-90b7-92e127e49dac","order_by":4,"name":"Taiki Kono","email":"","orcid":"","institution":"Department of Earth and Planetary Sciences, Graduate School of Science, Kyushu University","correspondingAuthor":false,"prefix":"","firstName":"Taiki","middleName":"","lastName":"Kono","suffix":""}],"badges":[],"createdAt":"2025-04-26 05:45:05","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-6532782/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-6532782/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":83014699,"identity":"0c8f5ce8-8e87-40c1-b08f-cc9bc0a7655a","added_by":"auto","created_at":"2025-05-19 06:05:07","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":5845060,"visible":true,"origin":"","legend":"\u003cp\u003eMap of Kyushu Island and the target region in this study.\u003c/p\u003e\n\u003cp\u003e(a) Location of Kyushu Island. The red ellipse and black line show the location of Beppu-Shimabara graben (BSG) and the median tectonic line (MTL), respectively. The red triangles show the location of active volcanoes in Kyushu Island. The black dashed rectangle indicates the target area of this study. (b) Map of the Kumamoto region. The purple lines indicate active faults. The blue and red stars show the hypocenter of the largest foreshock (Mj 6.5) and the mainshock (Mj 7.3), respectively. The red rectangles show the fault plane model of the mainshock (Asano and Iwata 2021). The blue line indicates the caldera rim of Mt. Aso. The black dots show the hypocenters picked and determined manually by the Institute of Seismology and Volcanology (SEVO), Kyushu University, from 16 April 2016 to 31 August 2023, at depths of 0 - 20 km.\u003c/p\u003e","description":"","filename":"Figure1Nagayamaetal.png","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/827a4875468cfcbd52df9b19.png"},{"id":83014704,"identity":"4f60c9a6-d02e-4f60-afac-584b3f4c2cd1","added_by":"auto","created_at":"2025-05-19 06:05:09","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":2248080,"visible":true,"origin":"","legend":"\u003cp\u003eMap of the seismic stations used in this study.\u003c/p\u003e\n\u003cp\u003e(a), (b) and (c) shows the distribution of the seismic stations in Period I, II and III, respectively. The blue inverted triangles indicate the seismic stations operated by Kyushu University, the National Research Institute for Earth Science and Disaster Prevention (NIED) and the Japan Meteorological Agency (JMA). The seismic stations installed by the Group for Urgent Joint Seismic Observation of the 2016 Kumamoto Earthquake (Shimizu et al., 2016) are shown by green inverted triangles respectively. The red circles indicate the high density seismic stations installed by Kyushu University around the Hinagu Fault Zone. The black rectangle and purple lines show the location of the target area and active faults.\u003c/p\u003e","description":"","filename":"Figure2Nagayamaetal.png","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/86d845b7634f5a7ed98999ee.png"},{"id":83014835,"identity":"3cc717f2-ae1d-4a08-892a-43e31664333b","added_by":"auto","created_at":"2025-05-19 06:05:36","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":4165014,"visible":true,"origin":"","legend":"\u003cp\u003eP and T axes projected horizontal planes in each period.\u003c/p\u003e\n\u003cp\u003e(a), (b) and (c) show P axes in Period I, II and III, respectively. The green lines indicate P axes with a dip angle more than 45°. (d), (e) and (f) show T axes in Period I, II and III, respectively. The yellow lines indicate T axes with a dip angle more than 45°.\u003c/p\u003e","description":"","filename":"Figure3Nagayamaetal.png","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/600f284d4084815b90159d15.png"},{"id":83012842,"identity":"957843e7-5e0c-41df-b7d6-92856c6b9ccf","added_by":"auto","created_at":"2025-05-19 05:39:45","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":1342359,"visible":true,"origin":"","legend":"\u003cp\u003eThe uniform regional stress field estimated in this study.\u003c/p\u003e\n\u003cp\u003e(a) The regional deviatoric stress tensor in the target area. The orange beach ball shows the regional stress tensor. The gray circles are the focal mechanisms used in estimation of the regional stress. (b) The principal compression stress axes of the regional deviatoric stress tensor. The direction of each stress axis is plotted on the lower hemisphere projection. The colored dots show the principal compression stresses of the regional stress field re-estimated when determining α and ξ within the 95% confidence interval in the bootstrap method; the blue, green and red dots indicate σ\u003csub\u003e1\u003c/sub\u003e, σ\u003csub\u003e2\u003c/sub\u003e and σ\u003csub\u003e3 \u003c/sub\u003erespectively. The black circle, cross and plus symbols show the principal compression stresses of the regional stress field estimated from original dataset\u003c/p\u003e","description":"","filename":"Figure4Nagayamaetal.png","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/da6a25d8b4ca5ab63ac2158b.png"},{"id":83014690,"identity":"9031c871-b906-49ad-93d5-35fb0c48cb94","added_by":"auto","created_at":"2025-05-19 06:04:16","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":2073941,"visible":true,"origin":"","legend":"\u003cp\u003eSee image above for figure legend\u003c/p\u003e","description":"","filename":"Figure5Nagayamaetal.png","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/c8c57472909b0acf2fcab99f.png"},{"id":83012841,"identity":"6bf05f9c-0fcc-4014-b9af-9dea02eac1b1","added_by":"auto","created_at":"2025-05-19 05:39:44","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":1053842,"visible":true,"origin":"","legend":"\u003cp\u003eSee image above for figure legend\u003c/p\u003e","description":"","filename":"Figure6Nagayamaetal.png","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/eb1181f444c44711b59de92d.png"},{"id":83012848,"identity":"b9891875-bbd3-4059-ae49-3091c2e4a6be","added_by":"auto","created_at":"2025-05-19 05:39:45","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":3132302,"visible":true,"origin":"","legend":"\u003cp\u003eStress field (a) before and (b) after the mainshock at a depth of 12 km.\u003c/p\u003e\n\u003cp\u003eThe color contour indicates the magnitude of the differential stress (σ\u003csub\u003e1\u003c/sub\u003e-σ\u003csub\u003e3\u003c/sub\u003e). The blue and red stars show the hypocenter of the largest foreshock (Mj 6.5) and the mainshock (Mj 7.3), respectively. The black dots show the hypocenters picked and determined manually by the Institute of Seismology and Volcanology (SEVO), Kyushu University, in each period. The black rectangles show the fault plane model of Asano and Iwata (2021), and the black dashed lines indicate the location of the fault plane at a depth of 12km.\u003c/p\u003e","description":"","filename":"Figure7Nagayamaetal.png","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/461447fabf6413c0a106b33b.png"},{"id":83012847,"identity":"32769cbc-fef2-4b3c-9745-99fb030aa43f","added_by":"auto","created_at":"2025-05-19 05:39:45","extension":"png","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":1622922,"visible":true,"origin":"","legend":"\u003cp\u003eDistributions of the average misfit angle.\u003c/p\u003e\n\u003cp\u003e(a) and (b) show the distributions in the shallow (\u0026lt; 10 km) and deep (\u0026gt; 10 km regions), respectively, for Period I. The size of the circle shows the number of the used seismic events in the grid (0.05°×0.05°) and the color of the circle indicates the average misfit angle in the grid. (c) and (d) are the same as (a) and (b) but for Periods II and III.\u003c/p\u003e","description":"","filename":"Figure8Nagayamaetal.png","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/6200b94b441ec04e0c7031bd.png"},{"id":83012854,"identity":"f1364e1e-28f2-4f0b-9629-5c43e293bf9f","added_by":"auto","created_at":"2025-05-19 05:39:45","extension":"png","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":775183,"visible":true,"origin":"","legend":"\u003cp\u003eConceptual illustration of the mechanism of stress concentration around the Futagawa and Hinagu Faults.\u003c/p\u003e\n\u003cp\u003eThe blue plane and red ellipses indicate the mainshock fault plane and the location of stress concentration. (a) The model before the earthquake sequence. The pink ellipses indicate the weak zone in the crust. Stress Concentration occurs around the weak zones. (b) The model after the mainshock. Stress concentration is induced by fault slip during the mainshock, and aftershocks, after slip or/and strain accumulation occur around the stress concentration regions.\u003c/p\u003e","description":"","filename":"Figure9Nagayamaetal.png","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/de551e3f3d2c1dcfacca65eb.png"},{"id":92600036,"identity":"0e47ebf7-5003-4f21-829c-8c75a3490bbb","added_by":"auto","created_at":"2025-10-01 13:58:15","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":24760372,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/a9546354-6aaf-4bbf-84e6-1cd35317789d.pdf"},{"id":83014694,"identity":"f09d605c-9ec1-4167-b115-ba91a278604e","added_by":"auto","created_at":"2025-05-19 06:04:22","extension":"pdf","order_by":13,"title":"","display":"","copyAsset":false,"role":"supplement","size":1811178,"visible":true,"origin":"","legend":"","description":"","filename":"AdditionalfileNagayamaetal.pdf","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/7ce124a0a07fb4e9865bca18.pdf"},{"id":83012856,"identity":"6768839c-bc08-45a5-ac08-c14852313eff","added_by":"auto","created_at":"2025-05-19 05:39:45","extension":"png","order_by":14,"title":"","display":"","copyAsset":false,"role":"supplement","size":295197,"visible":true,"origin":"","legend":"","description":"","filename":"GraphicalabstractNagayamaetal.png","url":"https://assets-eu.researchsquare.com/files/rs-6532782/v1/e304cf6f182060903f21b4be.png"}],"financialInterests":"","formattedTitle":"Detecting the stress anomaly before occurrence the 2016 Kumamoto Earthquake, Kyushu Island, Japan","fulltext":[{"header":"1 Introduction","content":"\u003cp\u003eEarthquakes occur when stress acting on a medium reaches its shear strength. Stress concentration and weakening of fault strength are primary factors by which the shear stress acting on the fault plane exceeds its strength. However, it is difficult to distinguish between the two. Thus, when evaluating the stress state around a fault zone, it is important to account for the presence of the stress concentration and the mechanism of its buildup.\u003c/p\u003e \u003cp\u003eVarious factors contribute to the stress loading and its concentration on the earthquake fault. For the case of stress loading on an over-riding plate in the plate subduction system, the interplate coupling between the continental and subducting oceanic plate causes stress loading and concentration. In addition, the anelastic deformation in the lower crust beneath a fault zone also causes the stress concentration in the fault zone (e.g., Iio et al., \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2002\u003c/span\u003e, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2004\u003c/span\u003e). Furthermore, the heterogeneity in the upper crust contributes to the stress build up and strain accumulation around the fault zone (e.g., Meneses-Gutierrez et al., \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2018\u003c/span\u003e). Nevertheless, it is important to investigate the above factors that are responsible for stress build, identifying the origin of the stress field is difficult because of a lack of information about absolute stress levels. As pointed out by studies estimating absolute stress fields (e.g., Mitsuoka et al., \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Terakawa et al., \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2025\u003c/span\u003e), it is required to detect spatial-temporal heterogeneous stress field with known absolute stress change such as fault slip of a large earthquake occurrence. For the 2016 Kumamoto Earthquake, the slip distribution on the mainshock fault plane (e.g., Asano and Iwata \u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e2021\u003c/span\u003e) is available, and a large amount of seismic data can be utilized due to the high seismicity prior to the mainshock. Therefore, we focus on the focal area of the 2016 Kumamoto earthquake sequence to obtain detailed insights into the stress field around the fault.\u003c/p\u003e \u003cp\u003eThe 2016 Kumamoto earthquake sequence occurred in the Hinagu and Futagawa Fault Zones in Kyushu Island, southwest Japan (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). The mainshock (Mj 7.3) occurred on April 16, 2016, following the largest foreshock (Mj 6.5) on April 14, 2016. The large fault slip and co-seismic stress change due to the mainshock were estimated from various kinds of the observation data, such as geodetic data (e.g. Fukahata and Hashimoto \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2016\u003c/span\u003e) and seismic data (e.g. Asano and Iwata \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2016\u003c/span\u003e, \u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Mitsuoka et al., \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Additionally, the presence of the post-seismic deformation was estimated from the geodetic data (e.g. Pollitz et al., \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Moore et al., \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Liu et al., \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). Around the Hinagu and Futagawa Faults, the spatially heterogeneous stress field in strike-slip or normal fault type was estimated in both pre- and post-mainshock periods by some previous studies (e.g., Matsumoto et al., \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2015b\u003c/span\u003e, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Mitsuoka et al., \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Furthermore, in the focal area, the presence of a low b-value and its spatial variation have been estimated prior to the earthquake sequence (Nanjo et al., \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). Therefore, these results might indicate stress concentration in the focal area of the 2016 Kumamoto earthquake. While some previous studies have examined the spatial-temporal changes and the absolute magnitude of the stress field in the target area (e.g. Mitsuoka et al., \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Terakawa et al., \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2025\u003c/span\u003e), there is no study that has discussed the distribution and factors of the stress concentration before and after the earthquake sequence around the earthquake fault. In this study, we model the stress field in the focal area of the 2016 Kumamoto earthquake sequence by combining the regional stress field and co-seismic stress change due to the mainshock and detect the stress concentration in the fault zone prior to the mainshock occurrence. Additionally, we discuss the factor of stress concentration in the source region from the perspective of crustal strength heterogeneity.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e"},{"header":"2 Data","content":"\u003cp\u003eWe used the dataset of focal mechanism in the target area (longitude: 130.5\u0026deg; \u0026minus;\u0026thinsp;131.2\u0026deg;E, latitude: 32.5\u0026deg; \u0026minus;\u0026thinsp;33.1\u0026deg; N, depth: 0-20km). We analyzed data for three periods; Period I, II and III are from 1993 to the Mj6.5 largest foreshock occurrence on 14 April 2016, from the mainshock on 16 April 2016 to 19 November 2021 and from 20 November 2021 to 31 August 2023, respectively. The seismic stations were operated by Kyushu University, the National Research Institute for Earth Science and Disaster Prevention (NIED) and the Japan Meteorological Agency (JMA). The station distributions for three periods are shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e. We note that high density seismic stations were installed by Kyushu University around the Hinagu Fault Zone in Period III. We determined focal mechanisms from P wave polarity data observed at 10 or more seismic stations by using the HASH algorithm (Hardebeck and Shearer \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e2002\u003c/span\u003e). The Institute of Seismology and Volcanology (SEVO), Kyushu University picked these P wave polarity data manually in Period I and II and automatically in Period III. Magnitude ranges for normalized stress tensor estimated using Matsumoto (\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2016\u003c/span\u003e) are 1.5\u0026ndash;4.0 for Period I and II and 0.5\u0026ndash;4.0 for Period III, respectively. Finally, we obtained 2681, 6251 and 591 good quality focal mechanism solutions for Period I, II and III, respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). In order to obtain co-seismic stress change, we adopted co-seismic slip distribution estimated by Asano and Iwata (\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). They estimated the slip distribution of the mainshock fault using the strong ground motion waveform data around the earthquake fault.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e"},{"header":"3 Method","content":"\u003cp\u003eWe modeled the stress field around the focal area of the 2016 Kumamoto earthquake sequence using the method proposed by Matsumoto et al. (\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e2015a\u003c/span\u003e). We assumed that the deviatoric stress tensor within the target region can be represented as sum of a uniform regional stress field and co-seismic stress change resulting from the mainshock. In order to express stress field taking into account the stress concentration around the focal area prior to the mainshock occurrence, we write the deviatoric stress tensor at position of the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:k\\)\u003c/span\u003e\u003c/span\u003e-th eartquake \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\varvec{x}}_{\\varvec{k}}\\)\u003c/span\u003e\u003c/span\u003e in the target region as\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$\\:\\begin{array}{c}{S}_{ij}^{b}\\left({\\varvec{x}}_{\\varvec{k}}^{\\varvec{b}}\\right)={S}_{ij}^{r}-\\left(1-\\alpha\\:\\right)\\varDelta\\:{S}_{ij}\\left({\\varvec{x}}_{\\varvec{k}}^{\\varvec{b}}\\right)\\#\\left(1\\right)\\end{array}$$\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Equb\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equb\" name=\"EquationSource\"\u003e\n$$\\:\\begin{array}{c}{S}_{ij}^{a}\\left({\\varvec{x}}_{\\varvec{k}}^{\\varvec{a}}\\right)={S}_{ij}^{r}+\\alpha\\:\\varDelta\\:{S}_{ij}\\left({\\varvec{x}}_{\\varvec{k}}^{\\varvec{a}}\\right)\\#\\left(2\\right)\\end{array}$$\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Equc\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equc\" name=\"EquationSource\"\u003e\n$$\\:\\begin{array}{c}{S}_{ij}^{r}=\\xi\\:{s}_{ij}\\:\\#\\left(3\\right)\\end{array}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003ewhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{S}_{ij}\\left({\\varvec{x}}_{\\varvec{k}}\\right)\\)\u003c/span\u003e\u003c/span\u003e is the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:i\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:j\\)\u003c/span\u003e\u003c/span\u003e components of the deviatoric stress tensor at\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\:{\\varvec{x}}_{\\varvec{k}}\\)\u003c/span\u003e\u003c/span\u003e, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{S}_{ij}^{r}\\)\u003c/span\u003e\u003c/span\u003e is the uniform regional deviatoric stress tensor, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\Delta\\:}{S}_{ij}\\left({\\varvec{x}}_{\\varvec{k}}\\right)\\)\u003c/span\u003e\u003c/span\u003e is the stress tensor equivalent to the co-seismic stress change at \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\varvec{x}}_{\\varvec{k}}\\)\u003c/span\u003e\u003c/span\u003e, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{s}_{ij}\\)\u003c/span\u003e\u003c/span\u003e is the normalized regional deviatoric stress tensor. The superscripts \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:b\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:a\\)\u003c/span\u003e\u003c/span\u003e represent the periods before the largest foreshock and after the mainshock, respectively. \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\xi\\:\\)\u003c/span\u003e\u003c/span\u003e determines the magnitude of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{S}^{r}\\)\u003c/span\u003e\u003c/span\u003e, while \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e represents the degree of the stress concentration around the earthquake fault prior to the earthquake sequence. When \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e is close to 0, the stress concentration prior to the earthquake is observed around the fault. Conversely, for \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e close to 1, the stress field before earthquake was relatively uniform compared to after the earthquake. Among these parameters, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{s}_{ij}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\Delta\\:}{S}_{ij}\\)\u003c/span\u003e\u003c/span\u003e can be estimated or calculated from the observed data and fault model. Therefore, we determined the unknown parameters \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\xi\\:\\)\u003c/span\u003e\u003c/span\u003e by a grid search analysis.\u003c/p\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003e3.1 Estimation of regional stress field\u003c/h2\u003e \u003cp\u003eWe estimated the regional deviatoric stress field in the target region using the method proposed by Matsumoto (\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). In this method, the deviatoric stress tensor is estimated from the sum of seismic moment tensors of earthquakes in the target region based on the plasticity theory. The deviatoric stress tensor can be obtained by the following formula:\u003cdiv id=\"Equd\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equd\" name=\"EquationSource\"\u003e\n$$\\:\\begin{array}{c}{s}_{ij}\\propto\\:\\sum\\:_{k=1}^{K}{M}_{ij}^{k}\\:\\#\\left(4\\right)\\end{array}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003ewhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{s}_{ij}\\)\u003c/span\u003e\u003c/span\u003e is the normalized deviatoric stress tensor, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{M}_{ij}^{k}\\)\u003c/span\u003e\u003c/span\u003e is the seismic moment tensor released by \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:k\\)\u003c/span\u003e\u003c/span\u003e-th earthquake, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:K\\)\u003c/span\u003e\u003c/span\u003e is the total number of earthquakes in the target area. \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{M}_{ij}^{k}\\)\u003c/span\u003e\u003c/span\u003e is obtained as the product of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{s}_{ij}\\)\u003c/span\u003e\u003c/span\u003e and the seismic moment \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{M}_{O}\\)\u003c/span\u003e\u003c/span\u003e. \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{M}_{O}\\)\u003c/span\u003e\u003c/span\u003e is determined using an empirical relationship between the seismic moment and the magnitude of each earthquake \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:M\\)\u003c/span\u003e\u003c/span\u003e as (Matsumoto et al., \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2016\u003c/span\u003e)\u003cdiv id=\"Eque\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Eque\" name=\"EquationSource\"\u003e\n$$\\:\\begin{array}{c}\\text{log}{M}_{O}=1.151M+10.548\\:\\#\\left(5\\right)\\end{array}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eWe estimated the regional deviatoric stress using data from both before (Period I) and after (Period II and III) the mainshock. However, there is difference in the number of events among the periods, so we estimated the normalized regional stress tensor by averaging the normalized stress tensors for Period I and Period II and III and re-normalizing.\u003c/p\u003e \u003cp\u003eAs a result of estimating regional stress field, the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\sigma\\:}_{3}\\)\u003c/span\u003e\u003c/span\u003e axis, the minimum principal compression stress axis of the stress field, is horizontal and oriented approximately N-S (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). This result is consistent with the previous studies on the stress field around the Futagawa and Hinagu Fault Zones in both pre- and post-mainshock periods (e.g. Matsumoto et al., \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2015b\u003c/span\u003e; Matsumoto et al., \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Mitsuoka et al., \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). In addition, the regional deviatoric stress tensor we estimated had a high stress ratio (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:({\\sigma\\:}_{2}-{\\sigma\\:}_{3})/({\\sigma\\:}_{1}-{\\sigma\\:}_{3})\\approx\\:0.86\\)\u003c/span\u003e\u003c/span\u003e,where \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\sigma\\:}_{1},\\:{\\sigma\\:}_{2}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\sigma\\:}_{3}\\)\u003c/span\u003e\u003c/span\u003e are the maximum, intermediate and minimum principal compression stresses, respectively). This high stress ratio might be attributed to the events with a normal faulting in the target region, and it consists with the results of Mitsuoka et al. (\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2020\u003c/span\u003e) estimating high stress ratio around the earthquake fault.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003e3.2 Calculation of stress tensor due to co-seismic fault slip\u003c/h2\u003e \u003cp\u003eWe calculated the stress change tensor \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\Delta\\:}{S}_{ij}\\left({\\varvec{x}}_{\\varvec{k}}\\right)\\)\u003c/span\u003e\u003c/span\u003e due to the co-seismic fault slip at each earthquake\u0026rsquo;s hypocenter \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\varvec{x}}_{\\varvec{k}}\\)\u003c/span\u003e\u003c/span\u003e using the formula developed by Okada (\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e1992\u003c/span\u003e). We adopted the slip distribution of the mainshock obtained by Asano and Iwata (\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e2021\u003c/span\u003e) from strong ground motion data. We assumed that the medium is a homogeneous half space with a rigidity of 33 GPa. Figure\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e shows the differential stress (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\sigma\\:}_{1}-{\\sigma\\:}_{3}\\)\u003c/span\u003e\u003c/span\u003e) of the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\Delta\\:}{S}_{ij}\\)\u003c/span\u003e\u003c/span\u003e around the focal area at a depth of 10 km. The differential stress exceeds 0.1 MPa in most of the target area and is significantly larger (\u0026gt;\u0026thinsp;10 MPa) near the mainshock fault, particularly around the Futagawa Fault, compared to the surrounding areas located farther from the earthquake fault.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003e3.3 Grid search analysis and definition of RMS residuals value\u003c/h2\u003e \u003cp\u003eWe determined two unknown parameters \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\xi\\:\\)\u003c/span\u003e\u003c/span\u003e in equations (1) and (2) using the grid search analysis. To determine the optimal pair of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\xi\\:\\)\u003c/span\u003e\u003c/span\u003e, we first calculated the residuals as the differences between the observed unit slip vector \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{s}}=({\\widehat{s}}_{1},\\:{\\widehat{s}}_{2},\\:{\\widehat{s}}_{3})\\)\u003c/span\u003e\u003c/span\u003e from the focal mechanism solutions at located \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\varvec{x}}_{\\varvec{k}}\\)\u003c/span\u003e\u003c/span\u003e and the unit vector \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{\\tau\\:}}=({\\widehat{\\tau\\:}}_{1},\\:{\\widehat{\\tau\\:}}_{2},\\:{\\widehat{\\tau\\:}}_{3})\\)\u003c/span\u003e\u003c/span\u003e of the maximum shear stress direction on the fault plane. The unit vector \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{\\tau\\:}}\\)\u003c/span\u003e\u003c/span\u003e was calculated from the deviatoric stress tensor \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:S\\)\u003c/span\u003e\u003c/span\u003e at the same position \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\varvec{x}}_{\\varvec{k}}\\)\u003c/span\u003e\u003c/span\u003e and matrix \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:A\\)\u003c/span\u003e\u003c/span\u003e formulated from the fault normal vector components through their combinations, represented by Michael (\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e1984\u003c/span\u003e). Then, we calculated the RMS (Root Mean Square) residuals value between \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{s}}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{\\tau\\:}}\\)\u003c/span\u003e\u003c/span\u003e from all focal mechanisms. Because there is a large difference in the number of the events between pre- and post-earthquake periods, to avoid bias in the results, we calculated the RMS residual value by weighting it with the inverse of the number of focal mechanisms each period. The \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:RMS\\)\u003c/span\u003e\u003c/span\u003e is defined as the following formula:\u003cdiv id=\"Equf\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equf\" name=\"EquationSource\"\u003e\n$$\\:\\begin{array}{c}RMS=\\sqrt{\\frac{\\left[\\frac{1}{3{K}_{b}}RS{S}_{b}+\\frac{1}{3{K}_{a}}RS{S}_{a}\\right]}{2}}\\:\\#\\left(6\\right)\\end{array}$$\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Equg\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equg\" name=\"EquationSource\"\u003e\n$$\\:\\begin{array}{c}RS{S}_{b}=\\sum\\:_{k=1}^{{K}_{b}}\\sum\\:_{i=1}^{3}{\\left({\\widehat{\\varvec{s}}}_{i,k}^{b}-{\\widehat{\\varvec{\\tau\\:}}}_{i,k}^{b}\\right)}^{2}\\:\\#\\left(7\\right)\\end{array}$$\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Equh\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equh\" name=\"EquationSource\"\u003e\n$$\\:\\begin{array}{c}RS{S}_{a}=\\sum\\:_{k=1}^{{K}_{a}}\\sum\\:_{i=1}^{3}{\\left({\\widehat{\\varvec{s}}}_{i,k}^{a}-{\\widehat{\\varvec{\\tau\\:}}}_{i,k}^{a}\\right)}^{2}\\:\\#\\left(8\\right)\\end{array}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003ewhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:K\\)\u003c/span\u003e\u003c/span\u003e is the total number of the focal mechanism data, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\widehat{\\varvec{s}}}_{i,k}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\widehat{\\varvec{\\tau\\:}}}_{i,k}\\)\u003c/span\u003e\u003c/span\u003e are \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:i\\)\u003c/span\u003e\u003c/span\u003e- component of the unit vectors \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{s}}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{\\tau\\:}}\\)\u003c/span\u003e\u003c/span\u003e for the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:k\\)\u003c/span\u003e\u003c/span\u003e-th focal mechanism, respectively, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:RSS\\)\u003c/span\u003e\u003c/span\u003e is the residual sum of squares between the unit vector \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{s}}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{\\tau\\:}}\\)\u003c/span\u003e\u003c/span\u003e and the subscript \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:b\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:a\\)\u003c/span\u003e\u003c/span\u003e correspond to before the largest foreshock (Period I) and after the mainshock (Period II and III), respectively. In this approach, a focal mechanism generally has two nodal planes. Therefore, to determine the appropriate fault plane, we calculated \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{\\tau\\:}}\\)\u003c/span\u003e\u003c/span\u003e at both planes and identified the one with the smaller \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:RSS\\)\u003c/span\u003e\u003c/span\u003e value as the fault plane.\u003c/p\u003e \u003cp\u003eThe grid search ranges are 1 kPa to 100 MPa and 0 to 1 for \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\xi\\:\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e, respectively. We determined the optimal pair of parameters by minimizing the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:RMS\\)\u003c/span\u003e\u003c/span\u003e. In addition, in Periods II and III, we used only events at least 1 km away from the mainshock fault plane.\u003c/p\u003e \u003c/div\u003e"},{"header":"4 Result","content":"\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003e4.1 Result of the grid search and the optimum stress field model\u003c/h2\u003e \u003cp\u003eAs a result of the grid search, the optimal pair of the parameters was determined to be \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:=0.71\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\xi\\:=5.62\\)\u003c/span\u003e\u003c/span\u003e MPa. This obtained \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e indicates that 29% of the co-seismic stress change was concentrated around the earthquake fault prior to the earthquake sequence. Furthermore, we estimated the confidence intervals for \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\xi\\:\\)\u003c/span\u003e\u003c/span\u003e using the bootstrap resampling. In this approach, focal mechanism solutions for each period were resampled 1000 times, allowing for duplication, and the 95% confidence interval was calculated. During this process, the regional stress field was also re-estimated for each resampled dataset (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb). As a result, the confidence intervals of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\xi\\:\\)\u003c/span\u003e\u003c/span\u003e were estimated to be 0.55\u0026ndash;1.00 and 3.16\u0026ndash;7.94 MPa, respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e). The confidence interval of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e suggests that a stress field with a stress concentration equivalent to 0\u0026ndash;45% of the co-seismic stress change before the earthquake sequence better explains the observed data.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWe calculated the deviatoric stress field in the target area using the previously determined optimal parameters (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:=0.71,\\:\\xi\\:=5.62\\)\u003c/span\u003e\u003c/span\u003e MPa) and Equations (1) and (2). Figure\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e shows the distribution of the differential stress at a depth of 12km before and after the mainshock, respectively. Before the earthquake sequence, there were two dominant high differential stress regions at the southern part of the Hinagu Fault and the central part of the Futagawa Fault at a depth of 12 km. In contrast, at the ends of the source fault plane (i.e. the southern edge of the Hinagu Fault, the northeastern of the Futagawa Fault and the junction of the Hinagu and Futagawa Faults), the differential stress was lower than that in the surrounding areas. On the other hand, after the mainshock, areas with high differential stress are located near the Hinagu and Futagawa Fault planes and the extension zones of both Faults, while differential stress is low in the central part of the Hinagu Fault and around the Futagawa Fault.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec9\" class=\"Section2\"\u003e \u003ch2\u003e4.2 The relationship between stress concentration and hypocenters of each earthquake\u003c/h2\u003e \u003cp\u003eTo verify the consistency of the stress concentration and observed seismic activity, we compared the distribution of stress concentration and hypocenters of earthquakes (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e). These hypocenters were determined using data manually picked by SEVO.\u003c/p\u003e \u003cp\u003eAs shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e, the low seismicity region corresponds to relatively low differential stress, while the high seismicity regions are associated with high differential stress during both the pre- and post-earthquake periods. In addition, by comparing the pre- (Period I) and post-earthquake sequence (Periods II and III), we found that stress built up, particularly at the southwestern end of the rapture zone of the Hinagu Fault, the western side of the Hinagu Fault, and the northern end of the Futagawa Fault. These areas were characterized by low seismicity in Period I and high seismicity in Periods II and III. Therefore, both spatially and temporally, there was a strong correspondence between the distribution of the differential stress and the seismic activity around the earthquake fault. However, the correspondence was not observed in some regions around the fault. For example, in Period I, although the differential stress around the central part of Futagawa Fault and the southern part of the Hinagu Fault is significantly higher than in the surrounding areas, seismic activity is not particularly high. Furthermore, to the north of the central part of the Futagawa Fault and the southeast part of the Hinagu Fault, an aftershock cluster with high activity exists, while the differential stress is low in Periods II and III. On the other hand, in the central part of the Futagawa Fault, a few aftershocks were observed, despite high differential stress in the stress field in Periods II and III. These inconsistencies between the distribution of hypocenters and differential stress is further addressed in the Discussion section.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec10\" class=\"Section2\"\u003e \u003ch2\u003e4.3 Misfits of the stress field model\u003c/h2\u003e \u003cp\u003eTo evaluate the optimal stress field model, we calculated the misfit angles between the unit slip vector \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{s}}\\)\u003c/span\u003e\u003c/span\u003e of focal mechanisms and the expected unit vector \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{\\tau\\:}}\\)\u003c/span\u003e\u003c/span\u003e oriented to the direction of the maximum shear stress on the fault plane from the stress tensor. We defined the misfit angle \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e as\u003cdiv id=\"Equi\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equi\" name=\"EquationSource\"\u003e\n$$\\:\\begin{array}{c}{\\theta\\:}_{mis}=\\text{arccos}\\left(\\widehat{\\varvec{s}}\\bullet\\:\\widehat{\\varvec{\\tau\\:}}\\right)\\:\\#\\left(9\\right)\\end{array}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eWe calculated \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{\\varvec{\\tau\\:}}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e at the two nodal planes of the focal mechanism and selected the plane with smaller \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e as the fault plane, in the same way as the previously defined RMS value. As a result of misfit angle \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e, some areas had large \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e values (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eIn the stress field prior to the earthquake sequence, the areas with large \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e were located to the north, away from the Futagawa Fault, and in the shallow region (0-10km depth) at the northeastern end of the Futagawa Fault. On the other hand, in the stress field in Periods II and III, some regions with particularly large \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e located in the southeast of the Hinagu Fault source area and near the junction of the Hinagu and Futagawa Faults. The potential causes of the large \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e value described above are also examined in detail in the Discussion section.\u003c/p\u003e \u003c/div\u003e"},{"header":"5 Discussion","content":"\u003cp\u003eAs described in Result section, there are some inconsistencies between the distribution of hypocenters and differential stress (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e). The inconsistency, in Period I, characterized by low seismic activity despite high differential stress in the central part of the Futagawa Fault and the southern part of the Hinagu Fault, may be attributed to higher crustal strength in these regions relative to its surrounding areas. This suggests that greater stress could be sustained without being released by earthquakes in these regions. In addition, other inconsistencies also observed in Periods II and III, namely, high seismic activity exists while the differential stress is low around the north of the central part of the Futagawa Fault and the southeast part of the Hinagu Fault and low seismic activity exist while the difference stress is particularly high in the central part of the Futagawa Fault. A possible explanation for this exception is the difference in geometry between the assumed and the actual fault planes. In this study, we adopted the flat fault plane model proposed by Asano and Iwata (\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). However, the actual fault plane may not be perfectly flat. In addition, the after-slip pointed out by some previous studies (e.g., Pollitz et al., \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Moore et al., \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Liu et al., \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2024\u003c/span\u003e) may have affected seismic activity by releasing stress. Alternatively, in this region, it is also possible that strain accumulated by the mainshock remains without releasing stress.\u003c/p\u003e \u003cp\u003eMoreover, as mentioned in Result section, there are some regions which have large \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e value in the target area. In Period I, large \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e values were located to the north away from the Futagawa Fault, and in the shallow region at the northeastern end of the Futagawa Fault (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003ea and \u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003eb). The former region is consisting of multiple clusters of earthquakes. Since these earthquakes are away from the mainshock fault, they may have been triggered by factors not considered in this study\u0026rsquo;s stress field model, such as local stress heterogeneity, resulting in large \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e. The large misfit in the latter region is located near Mt. Aso, which also could be attributed to the heterogeneity of the stress field due to the strong heterogeneous structure (e.g., Savage et al., \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). On the other hand, in Periods II and III, some regions with particularly large \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e located in the southeast of the Hinagu Fault source area and near the junction of the Hinagu and Futagawa Faults (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003ec and \u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003ed). These regions might be affected by stress heterogeneity and/or the fluid intrusion that existed before the earthquake sequence, as reported by Terakawa et al. (\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2025\u003c/span\u003e). They estimated the absolute background stress field in the focal area of the 2016 Kumamoto earthquake from the perspective of the effective friction coefficient and strain energy by combining the background stress field and co-seismic stress changes due to both the largest foreshock and the mainshock. They suggested that there were some earthquakes occurred due to factors not accounted for in their stress field model. Even though we took into account the co-seismic stress change due to the largest foreshock in addition to the mainshock, the average values of the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e around the Hinagu Fault decreased only slightly, and large \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e values remained in this area compared to the surrounding regions (Figure S4). Therefore, the large \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\theta\\:}_{mis}\\)\u003c/span\u003e\u003c/span\u003e in these regions might have been caused by factors other than the regional stress field or co-seismic stress changes during the largest foreshock and the mainshock.\u003c/p\u003e \u003cp\u003eSome previous studies estimated the stress field in both pre- and post-mainshock periods around the Hinagu and Futagawa Fault zones. Matsumoto et al. (\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2015b\u003c/span\u003e, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2018\u003c/span\u003e) reported the stress field around the Hinagu and Futagawa Fault zones and their surrounding area during the pre-mainshock period. They estimated the spatially heterogeneous stress field in the target area. In particular, Matsumoto et al. (\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2018\u003c/span\u003e) reported the presence of the lateral heterogeneous and depth dependent stress field around the source fault, which implies that stress concentration factor \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e is not completely 1. On the other hand, the spatial heterogeneity of the stress field around the focal area after the mainshock was also estimated in detail by some previous studies (e.g., Yoshida et al., \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2016\u003c/span\u003e; Mitsuoka et al., \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Therefore, the stress field during the post-earthquake period also cannot be expressed by only a uniform regional stress field, namely the model with \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:=0\\)\u003c/span\u003e\u003c/span\u003e is contradictory to these findings. Following them, the optimal model with \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:=0.71\\)\u003c/span\u003e\u003c/span\u003e, which has the stress concentration equivalent to 29% of the co-seismic stress change, is qualitatively consistent with the stress field estimated in previous studies.\u003c/p\u003e \u003cp\u003eHowever, there are differences in the absolute magnitude of the deviatoric stress tensor compared to previous studies (Figure S6). The magnitude of the regional deviatoric stress tensor representing the regional stress in the target area was slightly smaller than that estimated by Mitsuoka et al. (\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2020\u003c/span\u003e) for the pre-earthquake period (=\u0026thinsp;7.8 MPa) using focal mechanism data, but there was no significant difference exceeding the 95% confidence interval. However, Terakawa et al. (\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2025\u003c/span\u003e), based on the the effective friction coefficient and strain energy by combining the background stress field and co-seismic stress changes due to both the largest foreshock and the mainshock, estimated the absolute magnitudes of the deviatoric stress tensors are 37\u0026ndash;65 MPa and 39\u0026ndash;70 MPa around the Hinagu Fault and Futagawa Fault depth at 10 km, respectively. The difference in the differential stresses between this study and Terakawa et al. (\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2025\u003c/span\u003e) might be attributed to the evaluation of the stress field in the vicinity of the earthquake fault. We have estimated the absolute stress by excluding the focal mechanism data in this study located near the mainshock fault plane, within 1.0 km, during Periods II and III. When we use the events whose distances from the fault plane of the mainshock are 0.5 km or larger in the grid search analysis, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\xi\\:\\)\u003c/span\u003e\u003c/span\u003e increases and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\alpha\\:\\)\u003c/span\u003e\u003c/span\u003e decreases (Figure S7). This leads to a larger magnitude of the deviatoric stress tensor before the mainshock. However, in this case, the RMS residual and the confidence interval are larger than those in the original case. This implies the possibility that localized high stress existed close to the earthquake fault and radiated large seismic energy, which was not captured in this analysis.\u003c/p\u003e \u003cp\u003eAs mentioned in the introduction, one of the factors that can cause stress loading on the fault is the anelastic deformation in the lower crust beneath the fault zone. In Kyushu Island, the distribution of anelastic deformation rate in the lower crust prior to the Kumamoto earthquake sequence was estimated by Yuasa and Matsumoto (\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2023\u003c/span\u003e) based on GNSS observation data. They estimated a higher anelastic strain rate (\u0026gt;\u0026thinsp;0.4 \u0026micro;/yr) in the lower crust from the BSG to the focal area of the Kumamoto Earthquake than its surrounding areas. Because they used large blocks (20 km \u0026times; 20 km) to estimate the distribution of anelastic strain rate, it is difficult to examine the detailed distribution of anelastic strain rate around the earthquake fault. However, near the central part of the Futagawa Fault, where the stress concentration was found in this study before the mainshock, is located just above the block with the highest anelastic strain rate (~\u0026thinsp;0.6 \u0026micro;/yr) in the focal area. In addition, the stress change rate in the upper crust (at a depth of 7 km) was caused by the anelastic deformation in the lower crust, which was \u0026gt;\u0026thinsp;8 kPa/yr around the Futagawa Fault. Therefore, the stress concentration, especially large in the Futagawa Fault zone, during the pre-earthquake period might have been caused by the higher anelastic deformation in the lower crust beneath the source fault than it was in the rest area.\u003c/p\u003e \u003cp\u003eIn addition, the strength heterogeneity in the upper and lower crust is also the factor causing the localized stress concentration. Some previous studies have revealed the strength heterogeneity in the crust and proposed that it might have caused strain accumulation and generation of large earthquakes (Usui et al., \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2024\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eIn the focal area of the 2016 Kumamoto earthquake sequence, various crustal structure surveys, such as velocity structure and electrical resistivity structure, have been conducted, and the heterogeneity in crustal structure was pointed out (e.g., Shito et al., \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Aoyagi et al., \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Aizawa et al., \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). In particular, the electrical resistivity structure, sensitive to the presence of fluid, provides insights into the distribution of the mechanical strength in the crust. Aizawa et al. (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2021\u003c/span\u003e) revealed the low-resistivity zones in the deep part of the Hinagu Fault. They interpreted these low-resistivity zones as weak structures; the deep ones associated with high-temperature magmatic fluid and the shallow ones enriched with clay. In this study, the stress concentration regions prior to the mainshock were estimated to be in the southeastern part of the Hinagu and the central part of the Futagawa Faults, which are adjacent to the low-resistivity zone identified by Aizawa et al. (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). Therefore, around the fault zone, such fluid-rich weak zones in the crust might have contributed to the strain accumulation in these regions and stress concentration around them.\u003c/p\u003e \u003cp\u003eBased on the above, it can be inferred that the stress concentration prior to the mainshock around the focal area of the 2016 Kumamoto Earthquake was caused by strain accumulation in the weak zones of both the lower and upper crust (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e"},{"header":"6 Conclusion","content":"\u003cp\u003eWe estimated the stress field model in the focal area of the 2016 Kumamoto earthquake, including the stress concentration before the earthquake sequence using high quality focal mechanism data and the slip distribution data of the mainshock fault plane. We applied the method proposed by Matsumoto et al. (\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e2015a\u003c/span\u003e) to the target region and determined the magnitudes of the stress concentration and regional stress field by the grid search analysis. As a result, the optimal model had the stress concentration equivalent to approximately 29% of the co-seismic stress change around the fault zone before the mainshock. This stress concentration distributed around the weak strength structure in the crust indicated by previous studies. This result indicates that the deformation could occurs in the low-viscosity structures in the lower crust and low-elasticity structures in the upper crust, and the stress concentrated around there and loaded on the earthquake fault. In addition, there were also the stress concentration around the fault zone after mainshock, and it might have driven aftershock activity and post-seismic deformation at the extended part of the earthquake fault.\u003c/p\u003e \u003cp\u003eWe revealed the relationship between the distribution of the stress concentration and the strength heterogeneity in the crust around the fault zone of the 2016 Kumamoto earthquake sequence. Therefore, in the future studies, estimating the heterogeneity of the strength structure and the stress field around the fault zone with high resolution may provide the presence of the stress concentration and the potential of the future seismic activity.\u003c/p\u003e"},{"header":"Abbreviations","content":"\u003cp\u003eBSG: Beppu-Shimabara graben\u003c/p\u003e\n\u003cp\u003eMLT: the median tectonic kine\u003c/p\u003e\n\u003cp\u003eGNSS: Global Navigation Satellite System\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eEthics approval and consent to participate\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConsent for publication\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eList of abbreviations\u003c/strong\u003e\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eBSG: Beppu-Shimabara graben\u003c/p\u003e\n\u003cp\u003eMLT: the median tectonic kine\u003c/p\u003e\n\u003cp\u003eGNSS: Global Navigation Satellite System\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAvailability of data and materials\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eDataset of focal mechanism data is available upon reasonable request. If any readers might be interested in accessing the data, please contact Satoshi Matsumoto ([email protected]).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare that they have no competing interests.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, under its The Third Earthquake and Volcano Hazards Observation and Research Program (Earthquake and Volcano Hazard Reduction Research).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthors\u0026apos; contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eYN analyzed the data and designed the research with help from SM. All authors contributed to the discussion of the research and read the final manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgements\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe used seismic data from the JMA and Hi-net (NIED). This study was supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, under its The Third Earthquake and Volcano Hazards Observation and Research Program (Earthquake and Volcano Hazard Reduction Research).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthors\u0026apos; information\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eDepartment of Earth and Planetary Sciences, Graduate School of Science, Kyushu University, Fukuoka 819-0395, Japan\u003c/p\u003e\n\u003cp\u003eYushi Nagayama \u0026amp; Kono Taiki\u003c/p\u003e\n\u003cp\u003eInstitute of Seismology and Volcanology, Faculty of Science, Kyushu University, Kyushu University, Shimabara 855-0843, Japan\u003c/p\u003e\n\u003cp\u003eSatoshi Matsumoto, Takeshi Matsushima \u0026amp; Kentaro Emoto\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEndnotes\u003c/strong\u003e\u003c/p\u003e\n"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eAizawa K, Takakura S, Asaue H, Koike K, Yoshimura R, Yamazaki K, Komatsu S, Utsugi M, Inoue H, Tsukamoto K, Uyeshima M, Koyama T, Kanda W, Yoshinaga T, Matsushima N, Uchida K, Tukashima Y, Matsushima T, Ishigara H, Muramatsu D, Teguri Y, Shito A, Mtsumoto S, Shimizu H (2021) Electrical conductive fluid-rich zones and their influence on the earthquake initiation, growth, and arrest processes: observations from the 2016 Kumamoto earthquake sequence, Kyushu Island, Japan. Earth, Planets and Space 73:1-12. https://doi.org/10.1186/s40623-020-01340-w\u003c/li\u003e\n\u003cli\u003eAoyagi Y, Kimura H, Mizoguchi K (2020) Seismic velocity structure at the southern termination of the 2016 Kumamoto Earthquake rupture, Japan. Earth, Planets and Space 72:1-14. https://doi.org/10.1186/s40623-020-01276-1\u003c/li\u003e\n\u003cli\u003eAsano K, Iwata T (2016) Source rupture processes of the foreshock and mainshock in the 2016 Kumamoto earthquake sequence estimated from the kinematic waveform inversion of strong motion data. Earth Planets and Space 68:147. https://doi.org/10.1186/s40623-016-0519-9\u003c/li\u003e\n\u003cli\u003eAsano K, Iwata T (2021) Revisiting the source rupture process of the mainshock of the 2016 Kumamoto Earthquake and implications for the generation of near‐fault ground motions and forward‐directivity pulse. Bulletin of the Seismological Society of America 111.5: 2426-2440. https://doi.org/10.1785/0120210047\u003c/li\u003e\n\u003cli\u003eFukahata Y, Hashimoto M (2016) Simultaneous estimation of the dip angles and slip distribution on the faults of the 2016 Kumamoto earthquake through a weak nonlinear inversion of InSAR data. Earth, Planets and Space 68: 1-10. https://doi.org/10.1186/s40623-016-0580-4\u003c/li\u003e\n\u003cli\u003eHardebeck, J L, Shearer P M (2002) A new method for determining first-motion focal mechanisms. Bulletin of the Seismological Society of America 92.6: 2264-2276. https://doi.org/10.1785/0120010200\u003c/li\u003e\n\u003cli\u003eIio Y, Sagiya T, Kobayashi Y, Shiozaki I (2002) Water-weakened lower crust and its role in the concentrated deformation in the Japanese Islands. Earth and Planetary Science Letters 203.1:245-253. https://doi.org/10.1016/S0012-821X(02)00879-8\u003c/li\u003e\n\u003cli\u003eIio Y, Sagiya T, Kobayashi Y (2004) Origin of the concentrated deformation zone in the Japanese Islands and stress accumulation process of intraplate earthquakes. Earth, Planets and Space 56.8: 831-842. https://doi.org/10.1186/BF03353090\u003c/li\u003e\n\u003cli\u003eLiu Y, Cui X, Hu Y, Zhang J, Chen Y (2024) Integrated investigation on heterogeneous lower crust rheology in Kyushu and afterslip behavior following the 2016 Mw7. 1 Kumamoto earthquake. Geophysical Research Letters 51.7: e2023GL107606. https://doi.org/10.1029/2023GL107606\u003c/li\u003e\n\u003cli\u003eMatsumoto S, Katao H, Iio Y (2015a) Determining changes in the state of stress associated with an earthquake via combined focal mechanism and moment tensor analysis: Application to the 2013 Awaji Island earthquake, Japan. Tectonophysics 649: 58-67. https://doi.org/10.1016/j.tecto.2015.02.023\u003c/li\u003e\n\u003cli\u003eMatsumoto S, Nakao S, Ohkura T, Miyazaki M, Shimizu H, Abe Y, Inoue H, Nakamoto M, Yoshikawa S, Yamashita Y (2015b) Spatial heterogeneities in tectonic stress in Kyushu, Japan and their relation to a major shear zone. Earth, Planets and Space 67: 1-9. https://doi.org/10.1186/s40623-015-0342-8\u003c/li\u003e\n\u003cli\u003eMatsumoto S (2016) Method for estimating the stress field from seismic moment tensor data based on the flow rule in plasticity theory. Geophysical Research Letters 43.17: 8928-8935. https://doi.org/10.1002/2016GL070129\u003c/li\u003e\n\u003cli\u003eMatsumoto S, Nishimura T, Ohkura T (2016) Inelastic strain rate in the seismogenic layer of Kyushu Island, Japan. Earth, Planets and Space 68: 1-7. https://doi.org/10.1186/s40623-016-0584-0\u003c/li\u003e\n\u003cli\u003eMatsumoto S, Yamashita Y, Nakamoto M, Miyazaki M, Sakai S, Iio Y, Shimizu H, Goto K, Okada T, Ohzono M, Terakawa T, Kosuga M, Yoshimi M, Asano Y (2018) Prestate of stress and fault behavior during the 2016 Kumamoto earthquake (M7. 3). Geophysical Research Letters 45.2: 637-645. https://doi.org/10.1002/2017GL075725\u003c/li\u003e\n\u003cli\u003eMeneses‐Gutierrez A, Sagiya T, and Sekine S (2018) Crustal deformation process in the Mid‐Niigata region of the Niigata‐Kobe tectonic zone as observed by dense GPS network before, during, and after the Tohoku‐Oki earthquake. Journal of Geophysical Research: Solid Earth 123.7: 6072-6085. https://doi.org/10.1029/2018JB015567\u003c/li\u003e\n\u003cli\u003eMichael A J (1984) Determination of stress from slip data: faults and folds. Journal of Geophysical Research: Solid Earth 89.B13: 11517-11526. https://doi.org/10.1029/JB089iB13p11517\u003c/li\u003e\n\u003cli\u003eMitsuoka A, Shito A, Matsumoto S, Yamashita Y, Nakamoto M, Sakai S, Iio Y, Shimizu H, Goto K, Okada T (2020) Spatiotemporal change in the stress state around the hypocentral area of the 2016 Kumamoto earthquake sequence. Journal of Geophysical Research: Solid Earth 125.9: e2019JB018515. https://doi.org/10.1029/2019JB018515\u003c/li\u003e\n\u003cli\u003eMoore J DP, Yu H, Tahg C, Wang T, Barbot S, Peng D, Masuti S, Dauwels J, Hsu Y, Lambert V, Nanjundiah P, Wei S, Lindsey E, Feng L, Shibazaki B (2017) Imaging the distribution of transient viscosity after the 2016 M w 7.1 Kumamoto earthquake. Science 356.6334: 163-167. https://doi.org/10.1126/science.aal3422\u003c/li\u003e\n\u003cli\u003eNanjo K Z, Izutsu J, Orihara Y, Furuse N, Togo S, Nitta H, Okada T, Tanaka R, Kamogawa M, Nagao T (2016) Seismicity prior to the 2016 Kumamoto earthquakes. Earth, Planets and Space 68: 1-10. https://doi.org/10.1186/s40623-016-0558-2\u003c/li\u003e\n\u003cli\u003eOkada Y (1992) Internal deformation due to shear and tensile faults in a half-space. Bulletin of the seismological society of America 82.2: 1018-1040. https://doi.org/10.1785/BSSA0820021018\u003c/li\u003e\n\u003cli\u003ePollitz F F, Kobayashi T, Yarai H, Shibazaki B, Matsumoto T (2017) Viscoelastic lower crust and mantle relaxation following the 14\u0026ndash;16 April 2016 Kumamoto, Japan, earthquake sequence. Geophysical Research Letters 44.17: 8795-8803. https://doi.org/10.1002/2017GL074783\u003c/li\u003e\n\u003cli\u003eSavage M. K, Aoki Y, Unglert K, Ohkura T, Umakoshi K, Shimizu H, Iguchi M, Tameguri T, Ohminato T, Mori J (2016) Stress, strain rate and anisotropy in Kyushu, Japan. Earth and Planetary Science Letters 439: 129-142. https://doi.org/10.1016/j.epsl.2016.01.005\u003c/li\u003e\n\u003cli\u003eShito A, Matsumoto S, Shimizu H, Ohkura T, Takahashi H, Sakai S, Okada T, Miyazaki H, Kosuga M, Maeda Y, Yoshimi M, Asano Y, Okubo M (2017) Seismic velocity structure in the source region of the 2016 Kumamoto earthquake sequence, Japan. Geophysical Research Letters 44.15 (2017): 7766-7772. https://doi.org/10.1002/2017GL074593\u003c/li\u003e\n\u003cli\u003eShimizu H, Iio Y, Sakai S, Okada T, Takahashi H, Watanabe T, et al. (2016) Urgent joint seismic observation of the 2016 Kumamoto earthquake\u0026mdash;Seismic activities and their background. Paper presented at the Japan Geoscience Union 2016 Meeting, Japan. 22-26 May 2016 https://confit.atlas.jp/guide/event/jpgu2016/subject/MIS34-02/detail\u003c/li\u003e\n\u003cli\u003eTerakawa T, Asano K, Urata Y (2025) Constraint on the background stress in the source region of the 2016 Kumamoto earthquake sequence based on temporal changes in elastic strain energies and coseismic stress rotation. Geophysical Journal International 240.1: 174-188. https://doi.org/10.1093/gji/ggae376\u003c/li\u003e\n\u003cli\u003eUsui Y, Uyeshima M, Hase H, Ichihara H, Aizawa K, Koyama T, Sakanaka S, Ogawa T, Yamaya Y, Nishitani T, Asamori K, Ogawa Y, Yoshimura R, Takakura S, Mishina M Morita Y (2024) Three‐dimensional electrical resistivity structure beneath a strain concentration area in the back‐arc side of the northeastern Japan arc. Journal of Geophysical Research: Solid Earth 129.5: e2023JB028522. https://doi.org/10.1029/2023JB028522\u003c/li\u003e\n\u003cli\u003eYoshida K, Hasegawa A, Saito T, Asano Y, Tanaka S, Sawazaki K, Urata Y, Fukuyama E (2016) Stress rotations due to the M6. 5 foreshock and M7. 3 main shock in the 2016 Kumamoto, SW Japan, earthquake sequence. Geophysical Research Letters 43.19: 10-097. https://doi.org/10.1002/2016GL070581\u003c/li\u003e\n\u003cli\u003eYuasa Y, Matsumoto S (2023) Anelastic deformation in the lower crust and its implications for tectonic dynamics in Kyushu, Southwest Japan. Tectonophysics 846 : 229674. https://doi.org/10.1016/j.tecto.2022.229674\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"The 2016 Kumamoto Earthquake, Kyushu, stress field, stress concentration","lastPublishedDoi":"10.21203/rs.3.rs-6532782/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-6532782/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eTo evaluate future seismic activity around fault zones, understanding the stress state and its buildup is crucial. There are two major factors of earthquake occurrence, the stress concentration and the fault strength weakening, but it is difficult to distinguish between the two. Therefore, when estimating the stress field around the fault zone, it is important to consider the stress concentration prior to the earthquake. In this study, we modeled the stress field incorporating pre-seismic stress concentration at the focal area of the 2016 Kumamoto earthquake (Mj 7.3), Kyushu, Japan, using focal mechanism data and slip distribution on the mainshock fault plane. The estimated stress field shows the presence of stress concentration around the earthquake fault, corresponding to approximately 29% of the co-seismic stress change before the earthquake sequence. This stress concentration may reflect the strength heterogeneity and crustal deformation in the upper and lower crust, as suggested in previous studies. In addition, we also found the stress concentration after the mainshock, and its stress concentration may contribute to aftershock and post-seismic deformation.\u003c/p\u003e","manuscriptTitle":"Detecting the stress anomaly before occurrence the 2016 Kumamoto Earthquake, Kyushu Island, Japan","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-05-19 05:39:40","doi":"10.21203/rs.3.rs-6532782/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"32c055a7-1afb-4812-bec5-c970f97fb3a9","owner":[],"postedDate":"May 19th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2025-10-01T13:25:54+00:00","versionOfRecord":[],"versionCreatedAt":"2025-05-19 05:39:40","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-6532782","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-6532782","identity":"rs-6532782","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}