Stress State Inferred from b Value and Focal Mechanism Distributions in the Aftershock Area of the 2005 West Off Fukuoka Prefecture Earthquake

The spatiotemporal stress states in the aftershock region of the 2005 West Off Fukuoka Prefecture Earthquake are examined via an analysis of the b values and focal mechanism solutions. The aftershocks are aligned roughly NW–SE, with the southeastern part of the aftershock region believed to correspond to the Kego Fault, which extends beneath the Fukuoka metropolitan area. This study reveals depth-dependent b values in the focal region, where the b values (b = 0.7–1.4) are generally higher above the mainshock depth (9.5 km) and lower (b = 0.5–1.0) at greater depths. The shallower region possesses a significant temporal increase in b values, whereas a lateral b value heterogeneity is observed in the deeper region. The b values (b ~ 1.0) near the mainshock are relatively high, whereas the northwestern and southeastern edges of the deep region have lower b values (b = 0.5–0.7). On the other hand, many of the focal mechanisms for the M≥3.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M\ge 3.5$$\end{document} events are located in the low b value area of the deep region. The stress tensor inversion results reveal a change in stress state from strike-slip to strike-slip/normal faulting. These findings imply that the stress state remains high and/or slightly decreased in the northwestern and southeastern parts of the deep region. These results and the findings of previous research on this earthquake sequence suggest that the likelihood of future large earthquakes along the southeastern part of the aftershock region should be considered relatively high.

The hypocenter was located at 130. 1616°E, 33.7434°N , and 9.5 km depth, with a best-fit focal mechanism solution showing left-lateral strike slip with a tension axis aligned N23 W-S23 E ). Many aftershocks occurred at 1-16 km depth, and were distributed along a 25-km-long NW-SE-trending linear feature Uehira et al. 2006). The largest aftershock (M J 5.8) occurred near the southeastern end of the aftershock region on 20 April 2005, followed by additional aftershock activity to the southeast (e.g., Uehira et al. 2006). Numerous research studies have also investigated the coseismic slip distribution of the mainshock, spatial distribution of static stress drops, postseismic deformation, and attenuation structure in the focal region of the 2005 West Off Fukuoka Prefecture Earthquake (e.g., Asano and Iwata 2006;Horikawa 2006;Iio et al. 2006;Matsumoto et al. 2009;Nakao et al. 2006;Nishimura et al. 2006). Iio et al. (2006) have provided important details for hazard assessments in the focal region via an investigation of the static stress drop distribution of the aftershocks using the waveforms from temporary seismic stations; they highlighted the possibility of a large stress concentration around the southeastern end of the aftershock area, which includes the largest aftershock.
Many active faults exist in northern Kyushu district (The Research Group for Active Faults of Japan 1991). The southeastern extension of the aftershock region is believed to correspond to the northwestern edge of the Kego Fault, which is known to be active fault (e.g., Okamura et al. 2009;Uehira et al. 2006). The Kego Fault runs beneath the Fukuoka metropolitan area with a population of about 1.5 million. Okamura et al. (2009) used acoustic exploration and sediment cores to investigate the paleoseismicity of the Kego Fault and estimated that the latest two major earthquakes occurred 4500-4000 yBP and 8500-6500 yBP, respectively. These short recurrence intervals are a stark contrast to the 15,500year recurrence interval estimated by Shimoyama et al. (2005). These geophysical and geological results suggest that the next major earthquake might impact segments of the Kego Fault that have not previously ruptured, which run through the central part of Fukuoka City. The National Institute of Advanced Industrial Science and Technology (AIST) (2005) calculated the Coulomb stress change induced by the 2005 mainshock, and found a positive stress change of up to 0.1-0.5 MPa around the southeastern and northwestern ends of the aftershock area. The probability of a M = 7.0 earthquake in the next 30 years on the Kego Fault (counted from 2005) increased about 7% due to the 2005 mainshock (AIST 2005). The metropolitan area would be devastated if such a major earthquake were to occur. However, details of the stress state have not been evaluated in and around the focal region, although the area remains seismically active.
Here the b values of the Gutenberg-Richter Law (Gutenberg and Richter 1944), which are indicators of the stress state (e.g., Scholz 1968Scholz , 2015, are investigated in the aftershock region of the 2005 West Off Fukuoka Prefecture Earthquake using a seismic catalog that was compiled from data recorded after the 2005 mainshock. Recently, the stress-dependent characteristics of the b values were used for real-time seismic hazard assessment. One of the main concerns after a large earthquake is whether or not a stronger subsequent event will occur. Gulia and Wiemer (2019) examined time series of b values in many aftershock sequences, and found that a larger subsequent event is more likely to occur if the b value decreases substantially. Real-time monitoring of the b value is thus considered useful to estimate whether a larger subsequent event is likely. Furthermore, the focal mechanisms are inverted to obtain the stress tensors and assess the present stress states in the study area. This inversion does not determine the magnitude of the deviatoric stress at seismogenic depths (e.g., Gephart and Forsyth 1984;Michael 1987Michael , 1984. However, the orientation of the maximum horizontal stress axis relative to the fault strike direction is useful to infer how the stress field evolves with time (e.g., Hardebeck and Hauksson 2001). The style of faulting, which depends on the directions of the principal stress axes, also indicates the absolute stress level from the Coulomb failure criterion (e.g., Jaeger and Cook 1979). This study aims to clarify the spatiotemporal stress states as well as resolve highly stressed areas with high likelihoods of nucleating large earthquakes in the focal region, through unified analysis of b values and stress tensor inversions.

b Values
The empirical relationship between the earthquake frequency and magnitude distribution, which is known as the Gutenberg-Richter Law, is a power law of the form: where NðMÞ is the cumulative number of earthquakes with magnitude equal to or greater than M, and a and b are constants (Gutenberg and Richter 1944;Ishimoto and Iida 1939). The constant b is termed the b value, and represents the ratio of small to large earthquakes. The b value depends on various factors, including the stress state, strength, material heterogeneity, and thermal gradient (e.g., Mogi 1962;Scholz 1968Scholz , 2015Urbancic et al. 1992;Warren and Latham 1970;Wyss 1973) and varies between 0.5 and 1.5 for different tectonic settings, with an average value of 1.0 (Frohlich and Davis 1993). The stressdependent characteristics of the b values are among the most widely studied earthquake parameters in the literature (e.g., Scholz 1968Scholz , 2015Schorlemmer and Wiemer 2005;Schorlemmer et al. 2004). The differential stress generally increases as the b value decreases, and vice versa. Many studies have examined the stress states in various tectonic settings using spatiotemporal b value distributions (e.g., Chiba 2019Chiba , 2020Ghosh et al. 2008;Nanjo and Yoshida 2018;Nanjo et al. 2016Nanjo et al. , 2019Schorlemmer and Wiemer 2005;Schorlemmer et al. 2004;Tormann et al. 2015). The highly stressed areas in any tectonic setting, such as regions with a large coseismic slip and large slip deficit rate, are usually characterized by low b values. For example, Schorlemmer and Wiemer (2005) examined the b value distributions in the focal area prior to the 2004 M w 6.0 Parkfield Earthquake on the San Andreas Fault (California, USA), and found a close spatial correspondence between low b values and large coseismic slip during the mainshock. Ghosh et al. (2008) found a negative correlation between the b value distribution and plate locking inferred from geodetic estimates of interface locking along the Cocos Plate near Nicoya Peninsula, Costa Rica. Nanjo et al. (2016) analyzed the b value distribution in the focal region of the 2016 Kumamoto Earthquake, Japan, prior to the mainshock, and discovered that a zone near the mainshock corresponds to a low b value region. These examples verify that spatiotemporal b value analysis is an important and useful tool for seismic hazard assessment in regions with a high likelihood of large earthquakes.

b Value Data and Computational Procedure
The unified earthquake catalog maintained by the Japan Meteorological Agency (JMA catalog) was used for the b value calculations in this study, which included 15,381 events that were recorded between 12:20 on 20 April 2005 (2005.3 in decimal year) and 23:59 on 31 January 2020 (JST) in the region bounded by 130.0°-130.5°E, 33.55°-33.85°N, and 0.0-20.0 km depth. Data within 1 month of the mainshock were excluded since the events recorded during this period were highly heterogeneous and the JMA catalog was incomplete. Figure 1 shows hypocenters of the M ! 0:5 events; this choice of magnitude threshold is explained below.
ZMAP (Wiemer 2001) was employed for the seismicity analysis of this dataset, including the b value computations. An estimation of the magnitude of completeness M c , which is the minimum magnitude in the earthquake catalog that satisfies the power law in Eq. (1), is crucial for the b value analysis. The detection limits of the JMA catalog have improved significantly since 2000; M c is typically 0.0-1.0 in the Kyushu district (Nanjo et al. 2010). The maximum curvature method (Wiemer and Wyss 2000) was employed to investigate temporal changes in M c in the study area using overlapping 500-event windows and a 20-event step. Generally, M c \0:5 throughout the analysis period (Fig. 2a). Gridding methods were further applied by calculating the spatial distribution of M c using the maximum curvature method (Wiemer and Wyss 2000). The M c crosssection for events within 5 km of line A-B in Fig. 1 was calculated using 0.5-km grid spacing and a constant sampling radius of 2.5 km, which represents the optimal radius value for b value computation of cross-sections, as explained below. M c varies across the study area, with M c ! 0:5 in the northwestern subregion (Fig. 2b). Two-step screenings were therefore employed to choose M c for the b value calculations. First, a total of 9324 M ! 0:5 events were used for further analysis, selected based on spatiotemporal M c distributions. Next, local M c values were recalculated at each grid point for M ! 0:5 events, using this catalog and the maximum curvature method (Wiemer and Wyss 2000); this was a precursory step to b value calculations due to observed spatial variations in M c . Events with magnitudes below each new M c value were discarded. In addition, the maximum curvature method is known to often underestimate M c by 0.2 on average (Woessner and Wiemer 2005); therefore, an M c correction of ?0.2 was applied. The analysis window was divided into two periods to investigate the temporal changes in b values: the approximately 8-month period after the mainshock (period 1: 4335 events, 12:20 JST on 20 April 2005-31 December 2005) and the 2006-2020 period (period 2: 4989 events, 1 January 2006-31 January 2020). However, it is insufficient to evaluate temporal changes in b values using only arbitrarily selected time windows. An additional analysis, which divided the analysis window into three periods, was performed to investigate the temporal changes in b values in detail: period a1 was 12:20 JST on 20 April 2005-23:59 on 31 August 2005, and had 3298 events; period a2 was 1 September 2005-31 December 2007, with 3017 events; period a3 was 1 January 2008-31 January 2020, with 3009 events. This analysis was only used for b value calculations because there are too few available focal mechanisms to constrain stress tensor inversions in three temporal subwindows.
Gridding methods were applied by calculating the maximum likelihood b values. Spatial b value Vol. 178, (2021) Stress State Inferred from b Value and Focal Mechanism Distributions distributions were constructed to produce plan-view maps and cross-sections using gridding intervals of 0.005°and 0.5 km for the plan-view maps and crosssections, respectively. Determining the optimal sampling radii for the b value calculations at each grid point is both necessary and nontrivial. Fewer grid points match the required minimum number of events to constrain the b values when the sampling radii are small. Conversely, the resolution of the regional b value heterogeneities is reduced when the sampling radii are too large. The optimal sampling radius is the largest radius that still resolves b value heterogeneities in detail (Schorlemmer et al. 2004). The optimal sampling radius was determined to be 2.5 km for the plan-view maps and cross-sections in this study after exploring a broad range of sampling radii.
The b value was then calculated at each grid point via the maximum likelihood method (Aki 1965 where M mean is the mean magnitude and M 0 ¼ M c À 0:05 for uniform 0.1 magnitude bins. The corresponding b values were not calculated if a grid node had \ 50 events within a 2.5-km radius. The heterogeneities in the b value distributions were quantitatively evaluated using the p test (Utsu 1992), which is defined as: À 2; N 1 and N 2 are the numbers of events within the volumes to be compared, and N ¼ N 1 þ N 2 . Smaller p values indicate more significant b value heterogeneities. Previous studies have found that b value heterogeneity is considered statistically significant if logp À 1:3ðp 0:0498Þ (Schorlemmer et al. 2004;Utsu 1992Utsu , 1999).

Stress Tensor Inversion Data and Computational Procedure
The focal mechanisms used for the stress tensor inversion were estimated from F-net data, which were recorded by the permanent seismic stations operated by the National Research Institute for Earth Science and Disaster Resilience (NIED). Twenty-three focal mechanisms (13 events during period 1 and 10 events during period 2), each with a variance reduction of [ 60% and trace data from [ 3 stations, were inverted to determine the best-fitting normalized stress tensor for each period. The stress tensor inversion was performed using the MSATSI software (Martinez-Garzon et al. 2014), a MATLAB package based on Hardebeck and Michael (2006). This method uses the linearized inversion scheme of Michael (1987), which is based on the Wallace-Bott hypothesis (Bott 1959;Wallace 1951), whereby an earthquake slip vector is parallel to the resolved shear stress on the fault plane, and a damped inversion method is used to avoid the apparent spatial variability resulting from dividing the analytical region into small subregions. The stress tensor inversion has four unknown parameters: the three principal stress axes, r 1 ; r 2 ; r 3 , and the stress ratio: which represents the relative stress magnitude. The uncertainties in the results were evaluated via 500 bootstrap resamplings (Michael 1987). The number of events used to estimate the uncertainty in each stress parameter was ! 20 in each trial. As mentioned in Sect. 2.2, stress tensor inversions were not performed for each of the three temporal subwindows because there were too few available focal mechanisms.

Results
The hypocenter distribution consists of a roughly NW-SE-oriented linear feature (Fig. 1). There is no clear migration of events between periods 1 and 2. The b value cross-sections are first calculated for the events within 5.0 km of transect A-B in Fig. 1 during both analysis periods (Fig. 3) based on the NW-SE alignment of the seismicity; this resolves depth-dependent characteristics with a transition at the mainshock depth (9.5 km). The shallower region has relatively high b values (b = 0.7-1.4) during both periods, whereas the b values (b = 0.5-1.0) are lower in the deeper region (Fig. 3). The reduction in b value (b = 0.5-0.7) is especially significant in the northwestern and southeastern parts of the deep region. Conversely, the characteristics of the b values near the mainshock hypocenter differ slightly from the in the areas near the mainshock (e.g., regions 5 and 7) during period 2 (Fig. 3). Namely, low b value areas at shallow depths near the mainshock are present during period 1 but disappear during period 2. The observed spatial b value heterogeneities for each period are statistically significant at the 99% confidence level based on p value tests (Utsu 1992), as shown in Fig. 4. Plan-view maps of the depth-dependent b value characteristics are shown in Fig. 5 for the shallow (0.0-9.5 km depth) and deep (9.5-20.0 km) regions during both periods, with the mainshock depth serving as the transition between these two regions. The basic features are the same as those in the crosssections: the shallow region possesses relatively high b values (Fig. 5a, c), whereas the b values in the deep region are low (Fig. 5b, d). This depth dependence is especially noticeable during period 2. The b values in the northwestern and southeastern parts of the aftershock region in the deep region are lower than those near the mainshock (Fig. 5b, d). The b values in the shallow region near the mainshock are lower than those in the surrounding shallow area during period 1 (Fig. 5a). Noticeable temporal increases in the b values are also observed in the shallow region (Fig. 5a,  c), whereas significant temporal variations are not found at greater depths (Fig. 5b, d).
The differences in b values between the two analyzed periods (Db ¼ b period2 À b period1 Þ were calculated separately for the shallow and deep regions (Fig. 6a, c) to qualitatively evaluate the temporal changes in b values; the corresponding p values indicate that these changes are statistically significant. The b values also exhibit statistically significant temporal increases ðp 0:0498Þ in many parts of the shallow region, including the area near the mainshock (Fig. 6b). Areas with significant temporal b value changes are not widely distributed at greater depths (Fig. 6d).
Results for the three temporal subwindows indicate no clear migration of events between them (Fig. S1). The general patterns in spatial and temporal distributions of b values are the same as when the analysis window is divided into two periods (Figs. S2-S4). However, using three subwindows captures some temporal evolution of the b values that was not clear when only two periods were used. The b value distributions in periods a1 and a3 are roughly similar to those of both subwindows in the two-period analysis (Figs. 3, 5, S2, S3); the b value distribution for period a2 appears to show a transition between periods 1 and 2. Most of the low b values in the shallow region near the mainshock disappear during period a2 (Figs. S2b, S3c). The observed low b value anomaly may be characteristic of an early stage during period 1. In addition, significant temporal increases in b values between periods a1 and a2 were found in many parts of the shallow region (Fig. S4a, . b values remained low at greater depths during all three periods; areas with significant temporal b value changes were not widely distributed (Fig. S4c, d, g, h). The relationship between the b values and focal mechanisms was then investigated (Figs. 5, S3). The focal mechanisms for the M ! 3:5 events, which are relatively large earthquakes for the analyzed region, were located in the low b value (b = 0.5-0.7) area of the deep region; this finding is consistent with the definition of the b value that moderate to large earthquakes are likely to occur in regions with low b value. A stress tensor inversion was performed for each period using these focal mechanisms. The stress state obtained for period 1 was dominated by strikeslip fault type (Fig. 7), with R ¼ 0:51 and an E-W-oriented maximum principal stress. A counterclockwise rotation of the maximum principal stress and changes in the stress ratio were determined for period 2, with a decrease in the stress ratio to R ¼ 0:29 and an ENE-WSW-oriented maximum principal stress. The rotation of the maximum principal stress was statistically significant at the 95% confidence level. This relatively low stress ratio represents a stress state that is indicative of strike-slip/normal faulting (Fig. 7). This observation implies that the stress state changed at greater depths, even though there were no significant temporal changes in b values.

Discussion
The spatiotemporal b value distributions are highly heterogeneous in the aftershock region of the 2005 West Off Fukuoka Prefecture Earthquake, with the stress tensor inversion results indicating that the stress state below 9.5 km depth changed from a strike-slip-dominated to strike-slip/normal faulting regime. These results imply that the stress state is heterogeneous across the study area. Possible causes of the observed stress heterogeneity in the analyzed region are discussed in this section by comparing the results obtained in this study with those from previous studies. The areas with a high likelihood of nucleating large earthquakes are also evaluated.
The key seismic characteristics found in this study include the overall depth dependence of the b value distributions. The shallow region possesses relatively high b values, whereas the b values in the deep region are low (Figs. 3, 5). Uehira et al. (2006) reported that the aftershock distribution in the central part of the study area was aligned on different fault planes that split at the mainshock depth, with a difference of * 10°between their respective strikes and dips. and meanings are the same as in Fig. 1 Furthermore, they noticed that the mainshock focal mechanism estimated from P-wave first-motion polarities was consistent with the deeper fault plane, whereas the moment tensor solution estimated from broadband seismic waveform data ) was consistent with the shallower fault plane. The large coseismic slip of the mainshock was mainly observed in the southeastern part of the aftershock region, near the mainshock, at shallower depths (e.g., Asano and Iwata 2006;Horikawa 2006). Uehira et al. (2006) suggested that the mainshock rupture initiated at the hypocenter and propagated downward, then ruptured the lower plane and propagated upward, with this last rupture segment releasing most of the associated seismic moment. These findings suggest that more accumulated stress was released at shallower depths. The high b values in the shallow region can therefore be interpreted as a low-stress state associated with significant stress release due to large coseismic slip. Nakao et al. (2006) found notable postseismic deformation on the coseismic fault at \ 3.0 km depth, with the postseismic slip on the fault accounting for up to * 13% of the coseismic slip. The b values in the shallow region are especially high in the uppermost layers of Fig. 3. These findings are all consistent with a low-stress state in the shallow region.
A comparison of the presented results with a b value analysis in the focal region of the 2016 Kumamoto Earthquake is useful for understanding the potential causes of the b value heterogeneities Their interpretation was that the heterogeneous b values in the focal region were due to postseismic deformation. The focal region was located at the boundary between an afterslip-dominated region and viscoelastic deformation-dominated region (Nanjo et al. 2019). They suggested that there was a local increase in the stresses at the boundary between regimes dominated by different postseismic deformation processes based on the aftershock decay law of Utsu (1961) and a seismicity rate model inferred from the stressing history (Dieterich 1994). Compared with the Kumamoto earthquake, the postseismic deformation associated with the Fukuoka earthquake was restricted to very shallow depths in the focal region (e.g., Nakao et al. 2006). The significant temporal changes in b values in the present study are mainly found in the shallow region, where afterslip occurred. Nakao et al. (2006) detected postseismic deformation using GPS time-series data from 21 March to 27 June 2005. Significant postseismic deformation was clearly seen around the aftershock area during the analysis period of the present study (after 2005-04-20). In addition, the Evidently, significant temporal increases in b values between periods a1 and a2 were observable in many parts of the shallow region (Fig. S4a, b). Therefore, the postseismic deformation following a mainshock may play an important role in the temporal evolution of the b values. Furthermore, the increase in b values near the mainshock in the shallow region, which approximately corresponds to regions 3 and 7 in Fig. 3, may also be ascribed to a decrease in the effective normal stress due to increased pore fluid pressure, which will be discussed below. These findings indicate that it is unlikely that the shallow region is the potential location of another large earthquake in the near future since a significant amount of stress was released by coseismic and postseismic slip. The heterogeneity of the b value distribution in the deeper region is then discussed. The southeastern extension of the aftershock region appears to correspond to the northwestern edge of the Kego Fault, which extends beneath the Fukuoka metropolitan area (e.g., Okamura et al. 2009;Uehira et al. 2006). An evaluation of the stress state in the deeper region is therefore considered crucial for the assessment of possible large earthquakes in the future. Matsumoto et al. (2009) investigated the seismic-wave attenuation (Q À1 ) structure in and around the aftershock region of the 2005 mainshock, and discovered a highly heterogeneous Q À1 structure in the focal region. Areas with large slip and high aftershock activity were located in low Q À1 regions, whereas a high Q À1 region was found along the southeastern edge of the aftershock region, near the largest aftershock. Matsumoto et al. (2009) suggested that the low and high Q À1 regions corresponded to relatively high and low fault strengths, respectively. The largest aftershock has the same strike angle as the Kego Fault . Matsumoto et al. (2009) posited a segment boundary between the mainshock fault and the Kego Fault, where the fault strength is weak, based on the difference between the strike angles of the mainshock and largest aftershock. Iio et al. (2006) investigated the spatial distribution of the static stress drops of the Fukuoka earthquake aftershocks from 23 March to 31 May 2005, and suggested the possibility of more stress being concentrated along the southeastern edge of the aftershock region. The b values along the southeastern edge of the deeper region are consistent with the results of Iio et al. (2006), and were low throughout the entire analysis period (Figs. 3, 5). These findings suggest that the southeastern extent of the aftershock region, including the largest aftershock, comprises a highly stressed area with a relatively low fault strength. This area may therefore have a high probability of experiencing a large earthquake in the future compared to northwestern part of the aftershock region, even though both subregions have low b values. The potential causes of the relatively high b values below the mainshock are considered based on the findings of other geophysical studies. Wang and Zhao (2006) investigated the 3D seismic velocity and Poisson's ratio structures in the epicentral area of the 2005 mainshock and found that the mainshock was located at the boundary between high-and lowvelocity regions in the upper crust, with the lower crust side corresponding to a low-velocity, high-Poisson's ratio region. They interpreted the low-velocity region beneath the mainshock to possess fluids associated with mantle upwelling due to the opening of the Okinawa Trough. This finding suggests that fluids in and around a mainshock source area may play an important role in earthquake generation (Wang and Zhao 2006). The relatively high b values in the region beneath the mainshock are therefore considered to be generated by a reduction in the effective normal stress due to abundant fluids in the focal region. The geophysical interpretation of Wang and Zhao (2006) is helpful for inferring the cause of the temporal increase in b values in the shallow region near the mainshock (regions 3 and 7 in Fig. 3). Fluid injection from the deeper region may continue if this increase in b values occurs due to an increase in pore pressure related to the fluids beneath the focal region.
Most of the focal mechanisms included in the stress tensor inversion were located in the Vol. 178, (2021) Stress State Inferred from b Value and Focal Mechanism Distributions northwestern and southeastern parts of the aftershock region, which possess low b value areas in the deeper region (Fig. 5). Unfortunately, the data are too sparse to invert for a stress tensor at each edge since only 13 and 10 events were available for periods 1 and 2, respectively. It should therefore be noted that the results of the stress tensor inversion show average stress states in both the deep northwestern and southeastern subregions during each period; despite this limitation, the inversion results revealed a change in the stress field that was characterized by counterclockwise rotation of the maximum principal stress axis and a change in R from strike-slip to strike-slip/ normal faulting. This change in R is useful to infer the differential stress near both edges of the aftershock region. The Coulomb failure criterion states that normal faulting events sustain the lowest differential stresses, strike-slip events are intermediate, and reverse faulting events have the highest stresses, assuming that the vertical stress is equal to the overburden pressure (e.g., Jaeger and Cook 1979). The observed change in R implies that the average differential stress decreased slightly in these subregions during the study period. Areas with the low b values appeared to decrease further during period 2 based on the b value cross-sections and plan-view maps (Figs. 3b, 6c, respectively), although the decreases in the low b value areas were not necessarily statistically significant at many nodes (Fig. 6d). The angle between the strike of the fault plane and maximum principal stress axis direction is used to infer the shear strength in the focal region (e.g., Hardebeck and Hauksson 2001). The maximum principal stress axis direction should be aligned $ 30 relative to the fault plane for a strong fault, which is controlled by the frictional force (Byerlee 1978). On the other hand, if the fault is weak, then the direction is predicted to be at a higher angle. The hypocenter distribution consists of a roughly NW-SE-oriented linear feature (Fig. 1). Therefore, the angle of r 1 relative to the mainshock fault strike was rotated, and increased, from period 1 to period 2 (Fig. 7a, c). This finding also supports the observation of a slight decrease in average differential stress. The potential cause of the slight decrease in the differential stress in the region with low b values is unknown because no significant postseismic deformation was observed in the deeper region ). However, a region with significant low b values does not necessarily lead to a future large earthquake. For example, Tormann et al. (2015) analyzed the b values along the subducting Pacific Plate off Japan over a period that included the 2003 M w 8.3 Tokachi-oki and 2011 M w 9.0 Tohoku-oki earthquakes, and found that the former occurred within an area of persistent low b values off Hokkaido, but not in a region of significant low b values. On the other hand, the latter occurred within a subregion of distinct low b values in the subsequent high-slip area of the mainshock (Tormann et al. 2015). The causes of these b value differences remain unknown. Whether a significant decrease in b value prior to a major earthquake is an ubiquitous phenomenon is a question that should be approached with caution, as some previous studies of major earthquakes found precursory increases in b values (e.g., Smith 1981). This point will be addressed in future research, where it can be considered in the context of other geophysical studies. It is noteworthy, however, that El-Isa and Eaton (2014) reviewed previous research on spatial and temporal variations in b values and found that the vast majority of previous papers support a decrease in b value before a major earthquake. A significant and/or moderate decrease in b value is therefore useful information when drawing inferences about stress accumulation. Here the deeper region along the southeastern edge of the aftershock region is considered to have a high likelihood of triggering future large earthquakes due to the low b values; however, there is a possibility that the differential stress decreased slightly in the focal region. Temporal changes in the b values that were related to the presence of fluids were also found near the mainshock hypocenter. An increase in b values related to a reduction in effective stress, which is due to the presence of abundant fluids, may also be an important precursor to the next large earthquake. Wiemer et al. (1998) reported that increased pore pressure lowered the effective stress and increased the b value, a phenomenon associated with the onset of an earthquake swarm at Long Valley Caldera, California. The possibility that seismic activity increases due to abundant fluids in the focal region cannot be ruled out.

Conclusions
This study examined the spatiotemporal b value distributions in the aftershock region of the 2005 West Off Fukuoka Prefecture Earthquake. The b value distributions possess distinct depth-dependent characteristics, with the mainshock depth forming the transition zone. Relatively high b values (b = 0.7-1.4) are distributed throughout the shallow region, whereas the b values (b = 0.5-1.0) in the deep region are low. The high b values in the shallow region are ascribed to significant stress release due to large coseismic and postseismic slip. Although the b values are generally low in the deep region, the area below the mainshock possesses relatively high b values (b * 1.0). This anomaly may reflect a decrease in the effective stress due to an increase in pore fluid pressure that is caused by abundant fluids in the focal region. The b values (b = 0.5-0.7) below 9.5 km depth in the northwestern and southeastern parts of the aftershock region remained low throughout the analysis period, although a stress tensor inversion using focal mechanisms revealed a slight decrease in the average differential stress in these areas. It is possible that the southeastern part of the aftershock region possesses a high likelihood of hosting future large earthquakes.