Characteristics of seismicity in the southern Okinawa Trough and their relation to back-arc rifting processes

The southern part of the Okinawa Trough forms a narrow back-arc rift basin where evidence for submarine volcanoes and active hydrothermal venting is observed. The region is also known to cause large crustal earthquakes frequently which often accompany a rapid increase in seismicity rate. Although such swarm-like activities are common in active volcanic regions and are considered to be primarily induced by crustal fluid flows, potential interactions between tectonic, magmatic and hydrologic processes have been poorly examined in the southern Okinawa Trough despite these processes happening in the proximity. Here, I report the spatial and statistical characteristics of seismic activity in the southern Okinawa Trough and discuss their relation to other rifting-related phenomena. Most of the earthquakes with magnitude greater than 5 are localized around the rift axes (the Yaeyama Rift and the Yonaguni Rift) where seismic reflection data indicated the presence of solidified magmatic intrusions into the shallow sedimentary layers. I found the areas around the rift axes show low b values of < 0.8 and suggest that stress changes directly induced by dike intrusions beneath the rift axes control the occurrence of swarm activities. On the contrary, regions with high b values (> 1.2) are found around the Ishigaki Knoll and the Hatoma Knoll. These two areas are located between the rift axes and the Ryukyu Islands and correspond to potential submarine volcanoes proposed by seafloor bathymetry and seismic reflection images. This result may constrain the location of the volcanic front in the region.


Introduction
The Okinawa Trough in southwestern Japan is an active back-arc basin behind the Ryukyu subduction zone and provides an interesting research field to study crustal earthquakes associated with continental rifting (Fig. 1). The southern part of the trough, in particular, has been experiencing intensive back-arc rifting with a rifting rate of ~ 5 cm/year (Nishimura et al. 2004) and is actively producing thousands of shallow earthquakes each year (e.g., Nakamura and Katao 2003). Historically, a large earthquake with a magnitude of 7.2 occurred in 1938 within the rift basin and caused a tsunami with a height of 1.5 m (Hatori 1988). Recent active-source seismic studies obtained the detailed images of background fault structures generating these earthquakes (Arai et al. 2017;Nishizawa et al. 2019). However, characteristics of seismicity and physical properties of the seismogenic zone have been poorly understood, which implicitly hampers the evaluation of seismic and tsunami hazards in the region.
In addition to the seismic activity, the southern Okinawa Trough hosts active volcanism and hydrothermalism. Regional bathymetric and seismic reflection surveys suggest numerous submarine volcanoes in the southern Okinawa Trough. For example, Watanabe Arai Earth, Planets and Space (2021) 73:160 et al. (1995) revealed several knolls including the Ishigaki Knoll ~ 40 km north of the Ishigaki Island (Fig. 1b). Misawa et al. (2020) reported similar volcanic features at the Hatoma Knoll (Fig. 1b). In addition to these potential submarine volcanoes, seven hydrothermal vent sites have been confirmed in the rift basin to date, one of which is located on the Hatoma Knoll ( Fig. 1b; Beaulieu and Szafranski 2020). These volcanic and hydrothermal features are found around the rift axes and on the seamounts/knolls where a number of normal faults are developed. Although it is well known that in volcanic areas, stress changes and fluid circulation associated with magmatic activities can be a driving force of crustal earthquakes (e.g., Vidale and Shearer 2006;Yukutake et al. 2011), the relative roles of magmatism and fluid flows in the seismogenic processes has not been fully examined in the study area.
To characterize the seismogenic zones and better understand the underlying tectono-magmatic processes in the southern Okinawa Trough, I investigated spatial, temporal and statistical characteristics of crustal earthquakes that provide important information on the complicated interactions between aforementioned crustal activities. This paper especially focuses on b value (ratio of the number of larger earthquake relative to smaller ones) derived from the frequency-magnitude relationship of earthquakes (Gutenberg and Richter 1944). A number of observations demonstrate that b value is inversely correlated with the differential stress and is also indicative of the involvements of high-temperature fluids/magmas in the earthquake generation (Scholz 2015;Nanjo et al. 2018). In this study, I show that the statistical characteristics of seismicity are spatially heterogeneous and correlate well with fault distribution and volcanism. This consideration may also provide potential source areas of future large earthquakes in the southern Okinawa Trough. Figure 1c shows the distribution of earthquakes with a magnitude greater than 2 from the earthquake catalog provided by the Japan Meteorological Agency (JMA). The map clearly shows that the distribution is highly heterogeneous in space and, particularly, there are three locations where large earthquakes (M > 5) in the last two decades concentrate; (1) the Miyako Seamount, (2) the western end of the Yaeyama Rift and (3) the Yonaguni Rift. Around these areas, normal faults are densely distributed and, consistently, the large earthquakes mostly exhibit a normal-fault mechanism with a north-south tensional axis. Other areas with active seismicity are around the Ishigaki Knoll and the Hatoma Knoll. However, these two locations have not produced such large earthquakes in the last two decades (Fig. 1c).

Seismic activity in the southern Okinawa Trough
The temporal characteristics are also intriguing. As shown in Fig. 2a, the southern Okinawa Trough experienced obvious rapid increases in seismicity rate several times since 2001. Among them, swarm events that occurred in 2002 and 2013 are located around the western end of the Yaeyama Rift and the Yonaguni Rift, respectively (Fig. 2b, d). The swarm activity in 2002 includes 26 earthquakes with magnitude greater than 5 within only 2 days from the onset of the swarm (Fig. 2c).
The events in 2013, on the other hand, sustained the high seismicity rate for several weeks (Fig. 2e). Similar swarm events around the Yonaguni Rift also occurred in 2001 and 2014. Although insufficient accuracy of the earthquake locations makes it difficult to pinpoint the source areas of the swarms, it is likely that these swarms are associated with back-arc rifting because, together with their spatial proximity to the rift axes, they occurred concurrently with the acceleration of the southward movement of the Ryukyu arc detected by the onshore geodetic observations (Tu and Heki 2017). Nakamura and Kinjo (2018) interpreted that the southward movement may be due to transient aseismic slips induced by dike intrusions in the Okinawa Trough. Earthquakes with magnitude greater than 2 and whose depth is shallower than 40 km are shown by gray dots. Earthquakes with a magnitude greater than 5 are highlighted by larger white circles. The black circle is the epicenter of a huge earthquake in 1938 (M7.2). The earthquake data are taken from the JMA catalog. Yellow squares are seismic stations used for the JMA catalog in the study area. White curves are faults, most of which are recognized as normal faults with east-west strike (The Headquarters for Earthquake Research Promotion 2016). d Seismic reflection image along the YA05 profile (Arai et al. 2017) showing an intrusive body beneath the Yaeyama Rift and a potential magmatic body (marked by a reflector with negative polarity) (See figure on next page.)

Spatial distribution of b value
To estimate b-value distribution, hypocentral and magnitude data from the JMA earthquake catalog were used. I derived the earthquakes which occurred in the area of 24-26°North and 122.5-126°East during January 2001 to August 2019. I then removed events whose location is deeper than 40 km and was determined by less than 5 stations because of their potential large location error. I used the events with a depth of down to 40 km in the analysis, but most of the events are likely to be crustal earthquakes because JMA events tend to be located deeper in offshore areas (Nakamura and Katao 2003). These criteria selected 52,070 events. The complete magnitude from this data set was estimated to be 2.0 using the goodness-of-fit test of 95% (Wiemer and Wyss 2000;Woessner and Wiemer 2005). This value is consistent with the previous study by Nakamura and Kinjo (2018). For sampling b value, I used the gridding technique by Schorlemmer et al. (2004). In this analysis, b value was calculated at each grid with 0.01 degree spacing with the maximum likelihood method (Utsu 1965;Aki 1965); where M and Mmin denote the mean magnitude and the minimum magnitude of the given samples and e indicates the base of natural logarithm, respectively. The standard deviation of b value was estimated using the following equation (Shi and Bolt 1982): where n is the total number of events of the given samples. To focus on the map-view distribution of b value, I used cylindrical volumes with a constant radius centered on the grid. In cases the number of events with magnitude equal to or larger than the complete magnitude in the volume was fewer than 50, I did not calculate b value.
To demonstrate the reliability of obtained distribution of b values, I performed a series of tests using different values for each parameter. First, I examined different values for the radius r of the cylindrical volumes (10, 12 and 15 km) and confirmed that an almost same pattern of b-value anomalies was reproduced (Additional file 1). However, samplings with radii smaller than 10 km significantly reduce the coverage. For these reasons, I decided to employ the b-value distribution with a radius of r = 12 km for interpretation and discussion. Tests using larger complete magnitude (2.3 and 2.5) and separate time windows confirmed that an almost same pattern of b-value distribution is recovered independently of these parameters (Additional file 1). I also examined the effects of uncertainty of earthquake locations (Additional file 1). Finally, the plot of the standard deviations show that the estimated errors of b values are smaller than 0.20 in the whole area where b value is calculated and supports that the obtained spatial variation in b value is significant (Additional file 1).

Discussion and conclusions
One of the most important findings in this study is that the rift axes in the southern Okinawa Trough (the Yaeyama Rift and the Yonaguni Rift) are characterized by low b values (Figs. 2, 3). The previous seismic surveys revealed on-axis intrusions beneath the rift axes (Arai et al. 2017;Figs. 1d, 4) and the coincidence of these low-b-value regions with earthquake swarms and intrusive bodies provides tectonic implications on the backarc rifting process. Figure 4 compares the variations in b values and the density of seismicity along the YA01 profile with the dike intrusion near the Yonaguni Rift. This figure demonstrates excellent spatial correlation of low b values and high seismicity rate with the intrusive body. Since the stress state in the seismogenic zone is one of the controlling factors of b values, and highly stressed areas, such as coseismic slip regions of a large earthquake, tend to show low b values (Scholz 2015), I suggest that the low b values around the rift axes can be explained by a high stress state caused by dike intrusions (Toda et al. 2002). This observation also supports the hypothesis by Nakamura and Kinjo (2018) that the swarm activities are induced by the dike intrusions. The integration of these results suggests that active intrusions of magma into seismogenic depths control the activation of crustal seismicity by enhancing stress concentration along the rift axes. One may have a concern that the low b values are an artifact due to the potential location difference between the large and small earthquakes arising from the different seismic stations likely used to locate earthquakes as a function of magnitude. This problem is often suspected in offshore areas far from the seismic network, such as outer rise regions of the Japan Trench (e.g., Marsan et al. 2013). Although this possibility as a cause for the low b values cannot be fully excluded, I can conclude that the event locations in the southern Okinawa Trough are better constrained for the following reasons. First, the western end of the Yaeyama Rift with low b values is only ~ 80 km away from the nearest station in the Ishigaki Island and multiple seismic stations are within ~ 100 km (Fig. 1c). These distances are even smaller than those in typical subduction zones (e.g., ~ 300 km in the Japan Trench). Second, the arched geometry of the islands (Yonaguni Island from west to Miyako Island to east) provides better azimuthal coverage of the seismic network (~ 120° as in Fig. 1c). Third, seismic velocity used to determine earthquake locations for the JMA catalog represents the average velocity structure of the onshore area of the Japanese Island (Ueno et al. 2002). This is similar to typical continental crust and thus suitable to the southern Okinawa Trough, which still holds continental crust (Arai et al. 2017) .  Fig. 3 a Map-view distribution of b values and b its relation to submarine volcanoes, hydrothermal vent sites and magmatic intrusions. A high-b-value zone between the rift axes and the island arc (red arrow) agrees well with the iso-depth contour of the subducting Philippine Sea plate of 80 km (Yamamoto et al. 2018) and may suggest a potential volcanic front. Contour lines for b values are drawn every 0.2. Other legends used here are same as Fig. 1 In contrast to the regions around the rift axes with low b values, high b anomalies (> 1.2), which is a typical feature for active volcanoes (Wiemer and McNutt 1997;Wyss et al. 1997;Nanjo et al. 2018), are detected around the Ishigaki Knoll and the Hatoma Knoll and south of the Yonaguni Rift. The Ishigaki Knoll and the Hatoma Knoll are thought to be submarine volcanoes based on seafloor bathymetry and seismic reflection images (Watanabe et al. 1995;Misawa et al. 2020). Arai et al. (2017) reported a potential magma chamber ~ 20 km northeast of the Ishigaki Knoll. Such a magma chamber is likely to contain more high-temperature fluids than the other areas and I suggest the volcanic fluids are a primary contributor for the high b values. The high b values in the south of the Yonaguni Rift are also consistent with the previous study by Lin et al. (2007). While previous studies suggest the existence of the volcanic front in the southern Okinawa Trough, its location has not been clearly defined yet (Lin et al. 2004;Nishizawa et al. 2019). Increasing evidence for linear distribution of volcanic features between the rift axes (the Yaeyama Rift and the Yonaguni Rift) and the Ryukyu Islands, together with high-b-value anomalies from this study, may support that there is a volcanic front along the Ishigaki Knoll and the Hatoma Knoll (Fig. 3). Interestingly, these regions with high b values agree well with the 80-km iso-depth contour line of the subducting Philippine Sea plate ( Fig. 3b; Yamamoto et al. 2018), which is slightly shallower than the global average of ~ 100 km (Syracuse and Abers 2006).
The difference in seismicity characteristics highlights the spatial variation in crustal deformation in the process of back-arc rifting. From the view point of structural seismology, the most important difference between the low-b-value regions (the Yaeyama Rift) and the high-b-value regions (the Ishigaki Knoll and the Hatoma Knoll) is the rifting structure. As discussed in Arai et al. (2017) in detail, the Yaeyama Rift is characterized as passive intrusion as it shows a clear topographic depression on the seafloor and the subsurface sedimentary layers dips toward the rift axis and its dipping angle becomes larger as it approaches the rift axis (Fig. 1d). This means that magma is not actively ascending beneath the Yaeyama Rift and alternatively the tectonic stretching primarily controls the rifting/faulting processes. So the on-axis magmatic intrusion is merely a consequence of crustal thinning, but is expected to increase the extensive stress in the surrounding areas once dike intrusions occur and thus contribute to producing earthquakes in its vicinity. Arai et al. (2017) also show that any low-velocity volume indicating a stable magmatic body is missing beneath the Yaeyama Rift. The lack is a typical feature of passive intrusion and implies that the crustal temperature is relatively low compared to the southerly volcanic regions. A hydrothermal vent is located in the central part of the Yaeyama Rift, but the earthquakes seem inactive here and more active at the both ends of the Yaeyama Rift (Fig. 1c). This may suggest that crustal faulting occurs more intensively at the eastern and western edges of the Yaeyama Rift. On the other hand, high-b-value regions correspond to the volcanic knolls (topographic highs) at the Ishigaki Knoll and the Hatoma Knoll. These bathymetric features suggest that active upwelling of magma occurs here, but extensive rifting does not play a role on the crustal deformation. As another piece of evidence, Arai et al. (2017) reported a  Figure S3a) and the magmatic intrusion along the YA01 profile. The seismic reflection image is from Arai et al. (2017) low-P-wave-velocity volume indicating a magma chamber ~ 10 km south of the Yaeyama Rift. The crust with a large amount of fluid and high temperature is unlikely to produce larger earthquakes, leading to high b values.
Finally, I point out an important implication on the seismic hazards in the southern Okinawa Trough. According to the JMA catalog, a huge earthquake (M7.2) occurred in 1938 around the Miyako seamount and the epicentral location agrees well with the area with low b values (Fig. 3). Although the earthquake location may not be accurate due to the limited onshore seismic network in the region, this coincidence implies that the b-value distribution can be used as potential source areas of future large earthquakes as proposed in other places (e.g., Nanjo et al. 2019). Such information on the likelihood of future earthquake is inevitably important to evaluate seismic and tsunami hazards. To fully understand the seismogenic process associated with the back-arc rifting processes and further refine the hazard assessment, it may be worth monitoring earthquakes using denser and close-in seafloor observation.
Additional file 1: Figure S1. Frequency-magnitude distributions of earthquakes in the southern Okinawa Trough. a The whole study area. b The Yaeyama Rift. c The Yonaguni Rift. d The Ishigaki Knoll and the Hatoma Knoll. Each location is indicated in the b-value map at the bottom. Red curves indicate cumulative numbers of events with a magnitude equal to or larger than each magnitude. Gray bars are event numbers for each magnitude bin. The complete magnitude was estimated to be 2.0 using the goodness-of-fit test of 95% (Wiemer and Wyss, 2000;Woessner and Wiemer 2005). Figure S2. Distribution of b value with different radii r of a 10 km and b 15 km. These two different parameters produce an almost same pattern of b value anomalies as that with r = 12 km, while samplings with r = 10 km significantly reduce the coverage. Contour lines are drawn every 0.2. Figure S3. Distribution of b value with different complete magnitude Mc of a 2.3 and b 2.5. These two different parameters also produce a similar pattern of b value anomalies as that with Mc = 2.0, while larger Mc significantly reduce the coverage. Contour lines are drawn every 0.2. Figure S4. Distribution of b value with different time period of a January 2001 to January 2010 and b February 2010 to August 2019. Although the number of seismic stations for JMA earthquake catalog increased from six to eight in February 2010 in the study area, any significant temporal change in b value distribution is not recognized. Contour lines are drawn every 0.2. Figure S5. Distribution of b value calculated by error-added locations of earthquakes. To examine the effects of location uncertainties on the b value distribution, I added a random location error of ± 5 km and ± 12 km to the longitude and the latitude of the earthquake locations from the JMA catalog and recalculated the b value distribution. The result successfully reproduces the general pattern of b value anomalies from the error-free locations (Fig. 3).