Heterogeneous b-value Distributions Measured Over an Extensive Region from the Northern Okinawa Trough to Southern Kyushu Island, Japan

Spatiotemporal b-value maps are presented for an extensive region from the northern Okinawa trough to off southern Kyushu Island, Japan. This region shows high seismic activity associated with a tectonic setting characterized by seafloor spreading caused by the opening of a back-arc basin in the Ryukyu arc–trench system and the presence of a left-lateral shear zone on southern Kyushu Island. The obtained spatial and temporal distributions of b-value for the analyzed period are highly heterogeneous, reflecting the influence of the tectonic features and processes of the region. A comparison of the results with other geophysical observations suggests deep fluids associated with mantle upwelling below the seismogenic layer are the dominant control on the observed b-value variations. In addition, the hypocentral area of the largest earthquake (M7.1), which occurred in 2015, corresponds approximately to a region with low b-values (b = 0.5–0.7). Another region with low b-values (b = 0.5–0.7) occurs in the west of the analyzed region off southern Kyushu Island, where the occurrence of large-magnitude events is not clearly recognized. During the analysis period, there were few moderate to large earthquakes in the analyzed region. Furthermore, some active faults that have the potential to generate large earthquakes of ~ M7 are located in those parts of the analyzed region with low b-values. These findings suggest a high likelihood of nucleation of large earthquakes in the focal region.


Introduction
High levels of seismic activity occur in an extensive region encompassing the northern Okinawa trough to off southern Kyushu Island, Japan, and events of moderate to large size have been frequently observed in this area [e.g., Earthquake Research Committee (ERC) webpage; Kakuta & Goto, 2002]. The northeastern part of the Okinawa trough is characterized by seafloor spreading caused by the opening of a back-arc basin in the Ryukyu arc-trench system (e.g., Letouzey & Kimura, 1986;Nakahigashi et al., 2004). This spreading caused the invasion of mantle fluids into the crust owing to the upwelling of mantle in the back-arc of Kyushu (e.g., Nakada & Kamata, 1988;Seno, 1999;Yanagi & Maeda, 1998). The three-dimensional (3D) velocity structure for the extensive area from western Japan to eastern China has also been derived, with this structure depending on partial melting associated with the upwelling of hot mantle around the northeastern edge of the extending Okinawa trough (Sadeghi et al., 2000). A large earthquake of magnitude M7.1, which is the largest event recorded in the focal region since the installation of dense seismic observations (ERC, 2015), occurred on 14 November 2015 to the west off Satsuma Peninsula on southern Kyushu Island. Nishizawa et al. (2019) demonstrated that the hypocentral area is undergoing tectonic rifting by using seismic reflection and refraction surveys. The southern part of Kyushu Island is located adjacent to the northern part of the Ryukyu arc and includes a local left-lateral shear zone cutting across about 32°N (Nishimura & Hashimoto, 2006;Nishimura et al., 2018;Wallace et al., 2009). Moderate to large doublet earthquakes, 5 km and 48 days apart, occurred in the active shear zone in Kagoshima Prefecture in southern Kyushu Island on 26 March (M6.5) and 13 May (M6.3) 1997, followed by intense aftershock activity (ERC 1997). Umeda et al. (2014) conducted 3D magnetotelluric modeling above 40 km depth around the source region of the 1997 Kagoshima Supplementary Information The online version contains supplementary material available at https://doi.org/10.1007/s00024-022-02958-5. doublet earthquakes and proposed that the invasion of mantle fluids into the crust may have led to the occurrence of these events. These findings suggest that seismicity in the focal region is influenced by several tectonic processes.
The estimation of local seismicity is important for understanding the stress state within the focal region with respect to the abovementioned tectonic context. Seismicity analysis is also useful for seismic hazard assessment. However, analyses of seismicity in the marginal parts of Kyushu Island relative to other parts of the island are lacking. Among various seismicity parameters, the Gutenberg-Richter relation (G-R law), which generates b-values, has been investigated by many studies. The G-R law is a power law that describes the distribution of earthquake frequencymagnitude: where N is the cumulative number of earthquakes with a magnitude of ! M, and a and b are constants (Gutenberg & Richter, 1944;Ishimoto & Iida, 1939). The constant b is the relative occurrence of large to small events and is commonly termed the ''b-value''. The proportion of small earthquakes increases as the b-value increases, and vice versa. b-values vary depending on tectonic setting, but the average is b * 1.0 (e.g., Frohlich & Davis, 1993). b-values are known to be inversely dependent on differential stress (e.g., Scholz, 1968Scholz, , 2015. Highly stressed areas, such as regions with a large slip deficit rate and large coseismic slip, commonly show low b-values (e.g., Ghosh et al., 2008;Schorlemmer & Wiemer, 2005;Tormann et al., 2015). Furthermore, decreases in bvalues have been reported for the period leading up to the occurrence of moderate to large earthquakes (e.g., Enescu & Ito, 2001;Nanjo et al., 2012Nanjo et al., , 2016Nuannin et al., 2005), and these decreases have been interpreted as reflecting a form of stress increase prior to the seismic events (Scholz, 1968;Tormann et al., 2015;Wyss, 1973). In contrast, high b-values have frequently been observed in volcanic and geothermal settings owing to a reduction in effective normal stress resulting from high-pore-pressure regimes (e.g., Bridges & Gao, 2006;Chiba & Shimizu, 2018;Farrell et al., 2009). Spatiotemporal mapping of bvalues would therefore be useful to infer the stress state under various tectonic contexts. For example, Enescu and Ito (2001) reported that a noticeable decrease in b-values was observed approximately 2 years before the occurrence of the M j 7.2 Hyogo-ken Nanbu earthquake (Kobe earthquake) in western Japan on 17 January 1995. Tormann et al. (2015) comprehensively analyzed the spatiotemporal distribution of b-values along the subducting Pacific Plate off northeast Japan. They found that parts of the eventual rupture area of the 2011 M j 9.0 Tohoku earthquake exhibited stress increase (b-value decrease) up to the occurrence of the 2011 event, almost complete stress release (b-value increase) during the event, and a rapid stress recovery after the event (b-value decrease). Furthermore, those authors identified a relatively homogeneous band of high bvalues at depths of \ 100 km, where dehydration and partial melting of the subducting slab release material that ascends and then eventually forms the magmatic material (e.g., Wyss et al., 2001).
Previous studies have conducted b-value analyses in parts of Kyushu Island. Nanjo et al. (2016) investigated the distribution of b-values in the focal region of the 2016 M7.3 Kumamoto Earthquake, central Kyushu, on 16 April 2016, prior to the mainshock, and discovered that a zone near the mainshock corresponds to a low-b-value region. Nanjo et al. (2019) subsequently investigated the distribution of b-values in the focal region of the 2016 Kumamoto Earthquake sequence after the mainshock and identified a highly stressed area further to the south of the mainshock hypocentral area. The area of low-b-values after the mainshock is considered to reflect the history of stressing owing to postseismic deformation (Nanjo et al., 2019). Chiba (2021) investigated stress states in the aftershock area of the 2005 M7.0 West Off Fukuoka Prefecture earthquake, northern Kyushu, on 20 March 2005, through unified analysis of b-values and stress tensor inversion using focal mechanisms. This study found that the stress state decreases slightly in the focal region during the period after the mainshock via an increase in b-values and a change in the stress field. These findings suggest that b-value analysis is useful for seismic hazard assessment based on the characteristics of stress dependence. Here, the spatial distribution of b-values is presented for the extensive region from the northern Okinawa trough to off southern Kyushu Island. Spatiotemporal stress states, as well as highly stressed areas with high likelihoods of nucleating large earthquakes in the focal region, are evaluated from relationships between b-values and other geophysical observations.

Data and Methods
The unified earthquake catalog maintained by the Japan Meteorological Agency (JMA catalog) was used for the calculation of b-values in this study. Station coverage is sparse in the west of the analyzed region. Precise hypocenter data are required for the measurement of accurate b-value distributions. Therefore, a total of 129,847 regular earthquake events, which consist of relatively high-to moderateprecision events with flags of K, k, and A according to the ''hypocenter determination flag classification'' in the JMA data file (https://www.data.jma.go.jp/svd/ eqev/data/bulletin/data/format/hypfmt_e.html) were used in this study. The quality flags of K and k represent high-and middle-precision hypocenters closely examined via manual inspection, respectively. The quality flag of A represents automatically processed hypocenters with middle precision. The analysis period is 1 January 2000 to 30 April 2021, and the spatial extent of the region is bounded by 30.3°-32.5°N, 127.8°-131.0°E, and 0.0-30.0 km depth. The detection limits of the JMA catalog have improved substantially since 2000, with the magnitude of completeness M c , which is the minimum magnitude in the earthquake catalog that satisfies the power law in Eq. (1), typically being 0.0-1.0 for Kyushu Island (Nanjo et al., 2010). Data since 2000 were thus used for this study. Figure 1a shows hypocenters of events with M ! 1:3; this choice of magnitude threshold is explained below. Location errors are 0.63 km in the horizontal dimension and 1.81 km in the vertical dimension for the analyzed region.
The program ZMAP (Wiemer, 2001) was employed for analysis of the selected seismic dataset. Gridding methods were applied by calculating spatial changes in the magnitude of completeness, M c , and the maximum-likelihood b-values, as follows. Values of M c were first determined. As the M7.1 earthquake occurred west off Kyushu Island on 14 November 2015 during the analysis period, the analysis window was divided into two periods (first period: 1 January 2000 to immediately before the M7.1 event; second period: 1 January 2016 to 30 April 2021) to investigate temporal changes in b-values. Values of M c were calculated for each period. Events that occurred during the period immediately after the M7.1 earthquake (about 6 weeks) were discarded because the catalog is highly heterogeneous for that interval. The maximum curvature method (Wiemer & Wyss, 2000) was employed to investigate spatial changes in M c in the study area using events of all magnitudes, with a gridding interval of 0.05°and a constant sampling radius of 25 km for the plan-view map, which is the same as the optimal radius value for b-value calculation, as explained below. M c has values of \ 1.0 in the eastern part of the analyzed region, but has values of [ 2.0 in many areas of the western part of the region for both analyzed periods, as the station coverage is sparse in the western part (Fig. 2). If a larger M c is used, the power law relation will hold for all areas of the analyzed region. However, the larger M c leads to a decrease in the amount of earthquake data and therefore increases the size of the uncertainty of the b-values. Two-step screening was thus employed for calculation of the b-values to prevent large uncertainties: (1) A total of 20,438 M ! 1:3 events (first period: 12,624 events; second period: 7814 events) in the analyzed region were firstly selected for further analysis; and (2) M c was recalculated at each grid point via the maximum curvature method (Wiemer & Wyss, 2000) prior to the b-value analysis owing to the observed spatial variation in M c , and events with magnitudes below each new M c value were discarded. It is known that the maximum curvature method often underestimates M c by 0.2 on average (Woessner & Wiemer, 2005), so a magnitude correction of ?0.2 was applied to the recalculated M c values to improve the precision of the b-value distributions.
Using a grid with a 0.05°spacing, b-value distributions were then calculated for the selected events. The appropriate radius must be determined to extract the events around each grid node for the bvalue calculations. Fewer grid points match the Vol. 179, (2022) Heterogeneous required minimum number of events for constraining the b-values when the sampling radii are small. Conversely, the resolution of regional b-value heterogeneities is reduced when the sampling radii are large. The optimal sampling radius is the largest radius that gives the most detailed resolution of the heterogeneity in b-value (Schorlemmer et al., 2004). The optimal sampling radius was determined to be 25 km for the present study after exploring a broad range of sampling radii for both analysis periods. The b-value was then estimated at each grid point using the maximum-likelihood method (Aki, 1965), as follows: where M mean is the mean magnitude, and M 0 ¼ M c À 0:05 for a 0.1 magnitude interval. The corresponding b-values were not calculated if a grid node contained fewer than 50 events within a sampling radius of 25 km.
The heterogeneities of the b-value distributions were evaluated quantitatively using the p-test (Utsu, 1992), which is defined as: , the numbers) of events in each group compared, and N ¼ N 1 þ N 2 . More significant b-value heterogeneities are indicated by smaller pvalues and are considered statistically significant if logp À 1:3ðp 0:0498Þ (Schorlemmer et al., 2004;Utsu, 1992Utsu, , 1999. The relationship between focal mechanisms and b-values was also investigated in this study. A total of 264 focal mechanisms with M C 3.3 (114 and 150 events during the first and second periods, respectively) from the National Research Institute for Earth Science and Disaster Resilience (NIED) were used, each with a variance reduction of [ 60% and waveform data from more than three stations, and were plotted by classification with respect to the plunges of P-, B-, and T-axes, following Frohlich (1992). In the classification scheme, strike-slip faults have a B-axis plunge greater than 60°, normal faults a P-axis Vol. 179, (2022) Heterogeneous b-value Distributions Measured greater than 60°, and reverse faults a T-axis greater than 50°.

Results
Figure 3a, b shows the spatial distribution of bvalue for the first and second periods, respectively. Both the spatial and temporal distributions of b-value are highly heterogeneous in the analyzed region. Volcanic areas such as Kirishima and Sakurajima show very high b-value anomalies (b = 1.2-1.6) during both periods. It should be noted that b-value computation was not performed in the vicinity of the 2015 M7.1 earthquake prior to the occurrence of the mainshock, owing to the low number of events there (Fig. 3a). The observed spatial heterogeneities in bvalue at selected points during each period are statistically significant at the 98% confidence level, as computed using the p-test (Eq. 3; Utsu, 1992) (Figs. 3a, b, and 4). However, the b-value generally decreased in the western and central parts of the analyzed region from the first to second periods. c Figure 3 Spatial distributions of b-values estimated using the maximumlikelihood method (Aki, 1965) for a the first period (1 January 2000 to immediately before the M7.1 event of 14 November 2015) and b the second period (1 January 2016 to 30 April 2021). Focal mechanisms from F-net data are superimposed on the b-value maps. The sizes of the focal mechanisms are proportional to their magnitudes. Green, blue, orange, and gray focal mechanisms represent strike, normal, reverse, and other types, respectively, based on the criterion of Frohlich (1992). Small black squares represent the regions where the p-test was performed for each period. Solid black ellipses represent characteristic regions with low b-values (labeled regions A, B, C, D, and a). Moment release distributions were calculated for the focal mechanisms within the characteristic regions during the second period ( Table 1). The dotted black ellipse represents the hypocentral area of the 1997 Kagoshima doublet earthquakes (labeled region b). The other symbols and meanings are the same as in Fig. 1a. c P-wave velocity perturbations at a depth of 40 km in the back-arc area of Kyushu ( modified from the original figure by Sadeghi et al., 2000). Differences in b-value between the two analyzed periods were calculated to qualitatively evaluate the temporal changes in b-value (Fig. 5a). The corresponding p-values for the difference map indicate that these changes are statistically significant in the region bounded by 31:4 to31:8 N and 129:8 to130:5 E (Figs. 5b, S1), corresponding to the eastern part of region D, as explained below.
In Fig. 3a, b, focal mechanisms from F-net data overseen by the NIED are superimposed on the bvalue maps for comparison with the results of the present study. Many events comprise focal mechanisms with strike-slip and normal faulting, consistent with previous studies (e.g., Kakuta & Goto, 2002). The focal mechanisms are generally located in the moderate-to low-b-value regions. This finding is consistent with the definition of the b-value whereby moderate to large events are likely to occur in regions with low b-values. Here, four characteristic regions with low b-values (b = 0.5-0.7) during the second period (labeled regions A, B, C, and D in Fig. 3b) were selected for more detailed investigation. Similarly, a characteristic region (labeled region a) with low b-values, which is located slightly to the south of region D, and the epicentral area of the 1997 Kagoshima doublet earthquakes with moderate bvalues (labeled region b) during the first period were also considered for further investigation (Fig. 3a). Regions A and B are regions where seismicity has been activated since the 2015 M7.1 event. The seismicity in region C continued through both periods. Earthquake density maps for each period demonstrate these heterogeneous seismic activities (Fig. S2). bvalues for the western part of region D are low for the first period, although those for the eastern part are not necessarily low (Figs. 3a, b, 5). However, the eastern part shows much lower b-values in the second period compared with the first (Figs. 3a, b, 5, S1), but the number of moderate to large focal mechanisms is small (Fig. 3b). Here, the energy release within the characteristic regions with low b-values using these focal mechanisms during the second period was calculated via the equation logE ¼ 4:8 þ 1:5M, where M is magnitude and E is energy in N Á m (Hanks & Kanamori, 1979;Kanamori, 1977). The moment release from the F-net focal mechanisms in region D is in the order of 10 15 Nm, and is substantially smaller (by 1-3 orders of magnitude) than those of other regions (Table 1). Meanwhile, the epicentral area of the 1997 Kagoshima doublet earthquakes (region b) shows moderate b-values (b * 1.0). As discussed in the following section, seismicity in the Kagoshima epicentral area may be strongly affected by the activity of fluids driven by upwelling asthenosphere from the Okinawa trough.

Discussion
The b-values in and around the volcanic areas of the study region, including Kirishima and Sakurajima, are very high, consistent with the general finding that b-values are typically higher for volcanic areas than for other tectonic settings (e.g., Frohlich & Davis, 1993;McNutt, 2005;Wyss et al., 1997). Chiba and Shimizu (2018) investigated spatial and temporal variations in b-values in and around Shinmoe-dake volcano, one of the most active volcanoes of the Kirishima volcano group, using data from the JMA volcanic earthquake catalog, and identified detailed spatiotemporal b-value variations associated with the volcanic activity. Nanjo et al. (2018) investigated proposed that hot fluids beneath the hypocentral area played an important role in the generation of the swarm activity. Thus, the relatively high b-values in and around the volcanic area calculated in the present study are also considered to reflect highly heterogeneous crustal structure and stress states related to volcanic activity. Detailed perturbations in generally high-b-value anomalies associated with volcanic activity in and around Kyushu Island were also discussed in the abovementioned studies.
Here, the potential causes of the observed heterogeneous spatial and temporal distributions of b-value are discussed with reference to the stress dependence of b-values and by comparison with other geophysical studies. Some previous studies have found the invasion of mantle fluid into the crust caused by rifting of the Okinawa trough and the existence of a local shear zone developed in the analyzed region (e.g., Nakada & Kamata, 1988;Nishimura & Hashimoto, 2006;Nishimura et al., 2018;Sadeghi et al., 2000;Umeda et al., 2014;Wallace et al., 2009;Yanagi & Maeda, 1998). Sadeghi et al. (2000) determined the 3D velocity structure to the west off Kyushu Island thorough tomography analysis. Their valuable tomographic results for the back-arc region of Kyushu were derived from depths of C 40 km, and plan-view maps of P-wave velocity perturbations at three depth slices of 40, 65, and 90 km were constructed (Fig. 3 of Sadeghi et al., 2000). The depth range of events used in the present study was 0 to 30 km. It is thus useful to compare the plan-view map of b-value distribution of this study and the seismic velocity perturbation at 40 km depth of Sadeghi et al. (2000) ( Fig. 3c, d). The comparison shows that low-b-value anomalies are located in areas that lack low-velocity anomalies at 40 km depth. In particular, the locations of regions A, B, C, D, and a, which are regions characterized by very low b-values during the first and second analysis periods, correspond to areas that lack such anomalies below the seismogenic layer (Fig. 3c, d). Furthermore, regions where b-values decreased significantly between the first and second periods (Figs. 5 and S1) also seem to be located within areas that lack low-velocity anomalies at 40 km depth. Sadeghi et al. (2000) interpreted the formation of the low-velocity anomalies in terms of partial melting associated with the upwelling of hot mantle around the northeastern edge of the extending Okinawa trough. However, many of the events used in the b-value analysis are located at depths of 20 km (Fig. 1b). There is a pronounced gap in depth between the observed b-values in maps and the velocity structure at 40 km depth. Here, b-values as a function of depth are calculated for characteristic regions with low b-values (regions a, A, B, C, and D) and the epicentral area of the 1997 Kagoshima doublet earthquakes (region b) in Fig. 3 to check whether the shallow seismicity was affected by deep fluids. The b-value distributions at depth are calculated using an event sliding window method, with overlapping 100-event windows shifted by five events per spatial step. Results are shown in Fig. S3. The bvalues in region b are higher for all depth ranges (Fig. S3b) compared with other regions. Region b corresponds to an area with low-velocity anomalies below the seismogenic layer (Fig. 3c). Umeda et al. (2014) performed 3D inversion of magnetotelluric sounding data in the source region of the 1997 Kagoshima doublet earthquakes and identified a nearvertical conductive zone with a width of 20 km  extending to the base of the crust and perhaps into the upper mantle toward the Okinawa trough. This vertical conductive zone extends to a depth of 10-40 km (Fig. 5 of Umeda et al., 2014). Furthermore, elevated 3 He/ 4 He ratios in groundwater samples, which are indicators of the emission of mantle-derived helium from the seismic source region, have also been observed in the focal region (Umeda et al., 2014). These high b-values in region b are thus considered to reflect the existence of fluids associated with the invasion of mantle fluids into the crust. Furthermore, a marine seismic structure survey conducted in the northernmost Okinawa trough also supports the existence of igneous activity in this area, associated with back-arc rifting (Nakahigashi et al., 2004). These findings suggest that the observed heterogeneity in b-values is controlled mainly by the upwelling of deep fluids. Regions A and B, where the b-value decreased significantly in many parts from the first to second periods (Figs. 5, S1), are located near the hypocenter of the 2015 M7.1 event. It is possible that stress changes due to the 2015 event affected the observed decrease in b-value, as seismicity in these regions has been activated since the M7.1 event (Fig. S2). Possible mechanisms of the decreased b-values in the eastern part of region D (Figs. 5, S1) remain unknown, but factors such as fluid depletion beneath the seismogenic layer and local stress concentration might have influenced the decreased b-values.
Many larger events (focal mechanisms with M C 3.3 from F-net data) occur in and around the hypocentral area of the 1997 Kagoshima doublet earthquakes during the first period, although the bvalues in this area are not necessarily low. Many focal mechanisms in this area are characterized by left-lateral strike-slip (Fig. 3a). This hypocentral area corresponds spatially to the left-lateral shear zone identified by Wallace et al. (2009), the position of which corresponds approximately to a block boundary between the North Ryukyu block and the outer arc of southwestern Japan determined by Nishimura and Hashimoto (2006). Nishimura and Hashimoto (2006) showed 2-4-mm/year left-lateral slip deficits along this block boundary and explained them in terms of an E-W-trending discontinuity in GPS velocities across southern Kyushu. Nishimura et al. (2018) also investigated strain partitioning and interplate coupling in and around the plate margin in southwestern Japan, including the Nankai Trough, using block modeling and onshore and offshore geodetic data. Their block geometry-based on active fault traces, shallow seismicity, and past block models-is more detailed than that of Nishimura and Hashimoto (2006). The study of Nishimura et al. (2018) clarified the high slip deficit rate (C 10 mm/ year) on the block boundary in southern Kyushu established by Nishimura and Hashimoto (2006) and corresponding to the left-lateral shear zone identified by Wallace et al. (2009). These findings suggest an increase in the local shear stress in the focal region. As mentioned above, Umeda et al. (2014) identified a near-vertical conductive zone in the source region of the 1997 Kagoshima doublet earthquakes using 3D inversion of magnetotelluric sounding data. Those authors pointed out that the possible existence of aqueous fluids in and below the seismogenic layer could change the effective stress in the focal region, and change the local stress state, leading to the occurrence of the 1997 Kagoshima doublet earthquakes. Furthermore, the Kagoshima earthquakes seem to be located in a region of relatively low velocity identified by Sadeghi et al. (2000), which those authors interpreted as being associated with partial melting associated with the upwelling of hot mantle around the northeastern edge of the extending Okinawa trough (Fig. 3c, d). The relatively high bvalues at depth in region b also support the existence of deep fluids (Fig. S3b). These tectonic processes involving local stress accumulation and/or a decrease in effective normal stress due to fluid supply from the low-velocity zone below the seismogenic layer have therefore been proposed as potential causes of the moderate b-values in the focal region. In northern Italy in September 1997, the Umbria-March sequence, with one with M w 5.7 and 9 h later another with M w 6.0, occurred (e.g., Ekström et al., 1998). Miller et al. (2004) found that seismic activities associated with this doublet earthquake sequence are not well explained by static stress transfer. They suggested that this sequence may be generated by fluid pressure changes from a deep source. Using a nonlinear diffusion model, Miller et al. (2004) demonstrated that the triggering amplitude of fluid 908 K. Chiba et al. Pure Appl. Geophys. pressure was in the order of 10-20 MPa in the case of this sequence, which overwhelms the typical value (0.1-0.2 MPa) in the static stress transfer model. Similarly, the 1997 Kagoshima earthquakes comprised a doublet sequence, and the largest asperity of the second event was located in a zone of stress shadow related to the first event (Horikawa, 2001). These findings suggest that the reduction in effective normal stress that results from the supply of fluid from the low-velocity zone below the seismogenic layer exerts a stronger influence on variation in bvalues in the focal region compared with a local stress increase caused by a local shear zone. The bvalues in the hypocentral area of the 1997 Kagoshima doublet earthquakes for the second period are the same as and/or slightly lower than those for the first period, but there are few focal mechanisms in the focal region (Fig. 3b). The 2016 Kumamoto earthquake sequence occurred in the area adjacent to the analyzed region during the second period. Local stress changes involved in that sequence might also have affected the observed heterogeneity in b-value in the studied region.
Regarding the low-b-value anomalies to the west off Kyushu Island, the hypocenter of the 2015 M7.1 earthquake is located along the margin of a low-velocity zone identified by Sadeghi et al. (2000) (Fig. 3c, d). The aftershock area of the 2015 M7.1 earthquake with low b-values, such as regions A and B, does not seem to expand into the low-velocity zones identified by Sadeghi et al. (2000) (Fig. 3d). These findings imply that a rupture has been initiated near the outer edge of the low-velocity zones and that the extent of the rupture is controlled by the velocity structure below the seismogenic layer. Aizawa et al. (2021) investigated the electrical resistivity structure in and around the area of the 2016 Kumamoto earthquake sequence and found that the ruptures that were initiated along the outer edge of the low-resistivity fluid-rich zones tended to cause relatively large earthquakes compared with those that initiated either away from or within the low-resistivity zones. This interpretation could also be applied to the present study, as it considers deep magmatic fluid zones to be one of the potential causes of low-resistivity zones, similar to the potential cause of the low velocity below the seismogenic layer. The key concept of this interpretation is that the spatial difference in the prefailure pressure/temperature (PT) gradient of the pore fluids contributes to the propagation and subsequent cessation of the earthquake rupture (Aizawa et al., 2021). When the rupture is initiated at one of the cracks near the edge of a low-resistivity zone, the rupture can propagate along other cracks, and these rupturing cracks can migrate from a high-PT zone to a low-PT zone because of the PT gradient between the two zones. This rupture propagation is considered to cease during attempted migration to a high-PT zone. The depth range of the events used in this study is not the same as that used by Sadeghi et al. (2000) for establishing their velocity structure. However, seismicity in the analyzed region is considered to be affected by the aqueous structure below the seismogenic layer on the basis of the findings of previous studies (e.g., Nakahigashi et al., 2004;Umeda et al., 2014) and the b-value distributions in the depth dimension (Fig. S3). The PT gradient of pore fluid pressure at 40 km depth in the uppermost mantle associated with the upwelling of hot mantle would affect the extent of rupture of events that are located at depths 30km for the present study. Consequently, it is considered that ruptures involved in moderate to large events do not extend across areas of moderate to high b-values (i.e., the low-velocity area below the seismogenic layer) in the analyzed region.
The above interpretation may be useful for seismic hazard assessment in region D, which shows persistent low b-values in its western part and a significant decrease in b-values in its eastern part from the first to second periods (Figs. 3a,b,5). This local low-b-value area seems to be located in a pocket surrounded by a low-velocity structure (Fig. 3d), which suggests that the invasion of mantle fluids into the crust is less pronounced in that area, resulting in high effective normal stress. Indeed, the b-values in region D are low and essentially invariant with depth (Fig. S3f). The dimensions of region D are about 100 km 9 50 km in the longitudinal and latitudinal directions, respectively. If a rupture occurs on the outer edge of the low-velocity zone, the occurrence of a large event of * M7 is possible since earthquake magnitude is related to surface rupture length, subsurface rupture length, and rupture area via log-linear regressions (Wells and Coppersmith, 1994). Vol. 179, (2022) Heterogeneous b-value Distributions Measured Furthermore, some active faults, such as the Koshiki and Ichiki faults, are located between the northwestern part of Satsuma Peninsula to off western Kyushu Island (ERC 2013a, b;Figs. 1, 3, 5). These faults are expected to have a potential earthquake occurrence of maximum magnitude * M7.5 based on geological and geophysical survey results, but the occurrence of such large events has not been recorded for these faults (ERC 2013a, b). The present study showed a lack of large-to moderate-size earthquakes in the focal region ( Fig. 3b and Table 1). This lack of such earthquakes suggests a large amount of stress accumulation in region D compared with other low-bvalue regions and might indicate a high likelihood of nucleation of large earthquakes in the focal region. The 1997 Kagoshima doublet earthquakes (M6.5 and M6.3 events) are also located on the outer edge of a very low-velocity region, but the velocity around the hypocenter locations is also quite low (Fig. 3c, d). It is therefore expected that the PT gradient of pore fluid pressure in the focal region is probably not high. However, the gradient in P-wave velocity in and around the occurrence region of the 2015 M7.1 earthquake and region D of the present study is large. It is quite likely, therefore, that region D may be a location of a future large earthquake.

Conclusions
In this study, maps of the spatial variation in bvalues were presented for the extensive region encompassing the northern part of the Okinawa trough to off southern Kyushu Island, Japan, which is characterized by tectonic settings consisting of seafloor spreading owing to a back-arc basin in the Ryukyu arc-trench system and a local left-lateral shear zone cutting across southern Kyushu Island. Both spatial and temporal distributions of b-value are highly heterogeneous in the analyzed region. A comparison of the presented results with other geophysical observations suggests that deep fluids associated with mantle upwelling below the seismogenic layer are the dominant control on generation of b-value variations. Results also show the existence of a region with low b-values (b = 0.5-0.7) to the west off southern Kyushu Island. There are few moderate to large earthquakes in this region, although the low b-values inherently express a high likelihood of moderate to large events. Furthermore, some active faults in the region are expected to have a potential occurrence of large earthquakes of * M7. These findings suggest a high likelihood of nucleation of large earthquakes in the focal region. Real-time monitoring of b-values may be a useful tool for seismic hazard assessment in the analyzed region.

Acknowledgements
Data for the study were obtained from the unified earthquake catalog of the Japan Meteorological Agency (JMA) (https://www.data.jma.go.jp/svd/ eqev/data/bulletin/index.html). The analysis of bvalues was performed using the MATLAB software package ZMAP (Wiemer, 2001) (http://www.seismo. ethz.ch/static/stat_2010_website/stat-website-pre2010/ www.earthquake.ethz.ch/software/zmap.html). Figures were prepared using GMT (Wessel and Smith, 2013). This study was partially funded by the Tokio Marine Kagami Memorial Foundation, Japan (EKAF320500). The author thanks the editor, Yangfan Deng, and two anonymous reviewers for helping to improve the manuscript.

Declarations
Conflict of interest The author declares no known competing financial interests or personal relationships that are directly or indirectly related to the work submitted for publication.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.