Influence of crustal lithology and the thermal state on microseismicity in the Wakayama region, southern Honshu, Japan

Here we investigate the influence of the lithology and thermal state of the upper crust on earthquake distributions beneath the Wakayama region, southern Honshu, Japan, to better understand the influence of crustal conditions on regional seismogenesis. The earthquakes are concentrated in the deeper sections of mafic belts and shallower sections of pelitic belts, based on a comparison of relocated hypocenters and estimated subsurface geological structures. We compare the frictional properties of pelitic rocks and basalt, as obtained from petrological experiments, with the hypocenter depth distributions in pelitic and mafic belts to assess the control of crustal lithology on the depth extent of regional seismicity. The earthquake distributions are consistent with the temperature ranges over which the respective rock types are expected to exhibit a velocity-weakening behavior, based on the petrological experiments. The results suggest that the occurrence of shallow intraplate earthquakes is controlled by the temperature- and lithology-dependent friction of the upper crust.


Introduction
The upper crust is generally considered a zone of seismic heterogeneity. This heterogeneity has been linked to physical conditions within the upper crust, including the heterogeneous distributions of active faults and subsurface lithologies, and spatial heterogeneities in the stress and temperature fields. Previous studies have reported a correlation between the distributions of inland earthquakes and the physical conditions (e.g., temperature, structure, and lithology) of lithotectonic belts (e.g., Magistrale and Zhou 1996;Albaric et al. 2009;Hauksson and Meier 2019). These studies have suggested that the lithological properties of the crust controls the depth extent and magnitude of earthquakes.
The strength profile and rock frictional properties of the crust are closely related to the observed seismicity (e.g., Scholz 1988). The relationship between the lower limit of the seismogenic zone and the depth to the brittle-ductile transition in the crust has been extensively investigated in seismically active lithotectonic belts (e.g., Albaric et al. 2009;Hauksson and Meier 2019). However, these previous studies did not discuss whether these lithotectonic belts extend downward through the seismogenic layer. Sibson (1977) originally proposed a crustal model based on a crustal strength profile derived from geological observations and experimental studies. Improved models have since been developed by comparing this model with natural earthquake Open Access *Correspondence: maeda.sumire@bosai.go.jp 1 National Research Institute for Earth Science and Disaster Resilience, 3-1 Tennodai Tsukuba, Ibaraki 305-0006, Japan Full list of author information is available at the end of the article observations and additional experimental studies (e.g., Sibson 1982;Scholz 1988;Verberne et al. 2020). The crustal strength profile varies with the predominant minerals in the crust (e.g., Scholz 1988;Wintsch and Yeh 2013), which indicates that the crustal strength profile is dependent on crustal lithology.
Another fundamental property of the crust is the frictional nature of crustal rocks (e.g., Dieterich 1979a, b;Ruina 1983;Scholz 1988Scholz , 1998. Scholz (1988) proposed a crustal model of rock frictional properties based on experimental studies, with the seismogenic layer corresponding to the region, where the friction rate parameter, a-b (Dieterich 1979a, b;Ruina 1983), is negative. Numerous experimental results have provided information on rock frictional properties for a variety of crustal lithologies (e.g., Blanpied et al. 1991;He et al. 2007;Hartog et al. 2012;Zhang et al. 2017;Verberne et al. 2020), and it is becoming possible to discuss the relationship between the seismogenic layer and rock frictional properties in a variety of crustal lithology.
It is generally thought that the lower limit of the background seismicity is controlled by the frictional properties of the rock and the limit is located shallower than the brittle-ductile transition zone (Scholz 1988). Therefore, the background seismicity should be sensitive to the change in the frictional properties better than the brittle-ductile transition although the rupture of large earthquakes can reach down to the transition zone (Scholz 1988). Moreover, the upper limit of the background seismicity cannot be explained by the vertical strength profile expected from a simple brittle model, while a rate-and state-dependent friction model suggests the existence of the upper limit as indicated by Scholz (1988).
Here we attempt to quantitatively interpret the variations in the depth extent of shallow inland seismicity by focusing on differences in lithotectonic belts and temperature-dependent rock friction to better understand the relationship between the depth extent of seismogenesis and the lithological properties of the crust at depth. We first describe the geological background of the study area, the Wakayama region, southern Honshu, Japan (Fig. 1). We then compare the relocated hypocenters, which are constrained via the doubledifference relocation technique (Waldhauser and Ellsworth 2000), with the estimated subsurface geological structure, as derived from Bouguer gravity anomalies, to clarify the detailed relationship between the hypocenters and lithotectonic belts. We conclude by discussing the possibility of this relationship being controlled by the temperature-dependent rock friction of lithotectonic belts.

Geological background
The area to the south of the Median Tectonic Line (MTL) in southwestern Japan is known as the "Outer Zone of Southwest Japan. " This outer zone consists mainly of an accretionary prism (sediments, volcanic rocks, and equivalent metamorphic rocks) above a descending slab. The Sambagawa, Chichibu, and Shimanto belts strike E-W in the Wakayama region, and are separated by the MTL, Aridagawa Tectonic Line (ATL), and Butsuzo Tectonic Line (BTL) from north to south, respectively (Fig. 1). The Sambagawa belt consists mainly of pelitic and mafic schists (e.g., Takagi 1986), and is traditionally divided into three zones (spotted, non-spotted, and Mikabu zones) based on lithology, geological structure, and metamorphic grade (e.g., Kurimoto 1993Kurimoto , 2013Nakayama 1983). The spotted zone is composed mainly of relatively high-grade mafic schist (metabasalt) that is are especially common around Mt. Ryumon (Fig. 1). The non-spotted zone is composed mainly of relatively low-grade pelitic schist. The Mikabu zone is composed mainly of relatively low-grade mafic rocks (basalt) and pelitic rocks. The mafic rocks typically occur north of the ATL in the Mikabu zone, whereas the pelitic rocks are mainly distributed south of the non-spotted zone. The weakly metamorphosed Chichibu belt consists mainly of alternating marine sandstone and mudstone (e.g., Wallis and Okudaira 2016). The non-metamorphosed Shimanto belt consists mainly of sedimentary and mixed rocks (e.g., Taira et al. 1988).

Hypocenter relocation
We used the double-difference hypocenter location technique (hypoDD; Waldhauser and Ellsworth 2000) to relocate the earthquakes, because this approach provides better constraints on the detailed structures of active faults in a high-seismicity area than absolute location methods (e.g., Bouchaala et al. 2013). We relocated 8761 events (M ≥ 0.0) that occurred during the period from 1 January 2015 to 31 December 2016 using the JMA2001 one-dimensional seismic velocity model (Ueno et al. 2002) and the compressional-and shear-wave arrivals from the JMA unified catalog (left panel of Fig. 2). Because the seismicity in the Wakayama region is almost stable in the long period as shown in Additional file 1: Fig. S1, we regarded this 2-year seismicity as the representative seismicity in this region. The inset in Additional file 1: Figure S2a shows the locations of the seismograph stations (https:// www. data. jma. go. jp/ svd/ eqev/ data/ bulle tin/ deck. html) used to relocate the hypocenters in this study. For further details on the hypoDD relocation method, see Waldhauser and Ellsworth (2000).

Estimation of subsurface lithological structure
Gravity methods provide information on the subsurface geology in a given region. The Bouguer gravity anomalies across the Wakayama region, which are shown in the left panel of Fig. 3 (Geological Survey of Japan 2013), coincide with the key lithological structures ( Fig. 1). Higher Bouguer anomalies are present, where mafic rocks are the predominant lithology in the spotted and Mikabu zones. The Bouguer correction removes the gravitational effect of the topographic mass between the geoid and the ground surface, yielding a positive correlation between topographic height and positive gravity anomalies, as shown in Fig. 3. This indicates that the assumed density of 2.67 g cm -3 is either too low or the material beneath the topographically high areas possesses a higher density. Higher Bouguer anomalies are also observed across the southern Kii Peninsula (Fig. 3), associated with the subducted Philippine Sea Plate.
We imaged the two-dimensional (2D) density structures along five survey lines (Fig. 3) across the Wakayama region using the Talwani method (Talwani et al. 1959;Cady 1980) and gravity data from the Japan Gravity Database (Geological Survey of Japan 2013). The purpose of the density structure analysis was to evaluate the subsurface structure of the two E-W-trending high-gravityanomaly zones in the Wakayama region. These zones, which evidently correspond to mafic rocks, have higher densities than those of zones of sedimentary rocks (e.g., Iwaya and Kano 2005). We, therefore, constructed models by assuming that the high-density mafic rocks at the surface extend to seismogenic depths.  We first investigated the influence on gravity anomalies of the Philippine Sea Plate. The Bouguer anomaly across the study area extends downward with the subducting Philippine Sea Plate (e.g., Hagiwara 1986). We focused on the gravity variations within the spotted and Mikabu zones after removing the linear trend related to the subducting Philippine Sea Plate (Additional file 1: Fig. S3). The linear trend was estimated via a least-squares analysis of the gravity data along the five survey lines shown in Fig. 3. We used an assumed density of 2.67 g cm -3 for both the Bouguer correction in Fig. 3 and the density of the pelitic rocks around the mafic rocks. We used a density of 2.83 g cm -3 for the mafic rocks in the study region, based on the calculated mean density of exposed igneous rocks in accretionary complexes across Japan (Iwaya and Kano 2005). The calculated density values can be found in the PROCK database (https:// gbank. gsj. jp/ prock/ welco me. html). The subsurface structures of the two zones with a high gravity anomaly were determined by trial-anderror until a reasonable correspondence was obtained between the detrended data and calculated gravity anomalies.

Results
The JMA catalog hypocenters and relocated hypocenters ( Fig. 2; See also Additional file 1: Fig. S2) show a clear tightening of the relocated hypocenters relative to the original JMA hypocenters. We divided the region into 0.04° × 0.04° bins, and estimated the upper limit, lower limit, and thickness of the seismogenic layer based on the relocated earthquakes within each bin. We defined the D10, D90, and D90-D10 indices as the upper limit, lower limit, and thickness of the seismogenic layer, respectively, with D10 and D90 representing the depths above which 10% and 90% of the relocated crustal earthquakes occurred, respectively. We then shifted the bins horizontally by 0.01°, estimated the upper limit, lower limit, and thickness of the seismogenic layer for each bin via the above-mentioned procedure, and superimposed the results to investigate their detailed, smoothed distribution. The distributions of the estimated D10, D90, and D90-D10 indices for the bins containing > 50 earthquakes are shown in Fig. 4.
The average depths to the upper and lower limits of the seismogenic layer are 5.05 km and 7.16 km, respectively, and the average thickness of the seismogenic layer  Fig. S5 for lines a-a' , c-c' , d-d' , and e-e'). The upper panel compares the detrended observed anomalies (red cross symbols) and calculated anomalies (cyan plus symbols), as synthesized from the 2D density structure model in the lower panel. The assumed densities are 2.67 g cm -3 for the background and 2.83 g cm -3 for the regions in green is 2.12 km in this region. It is not easy to evaluate the uncertainty in the absolute depth of the relocated hypocenters, because we used the double-difference method which basically provided us relative location errors only. However, the errors in the average depths to the upper and lower limits of the seismogenic layer are expected to be less than 1 km in the absolute sense, because the standard errors of the hypocenters from the JMA catalog used in this study are 0.50 km in the horizontal and 0.86 km in the vertical direction on average.
The northwestern part of this region is characterized by shallow upper and lower limits of the seismogenic layer (Fig. 4a, b). This area is considered to be strongly affected by medium-sized earthquakes that occur at shallow depths, as Mj4-class earthquakes have occurred in a linear NE-SW belt across the region at ~ 5 km depth (Additional file 1: Fig. S4). The Chichibu belt is characterized by a thick seismogenic layer (Fig. 4c) due to the presence of planar hypocenter clusters between 4 and 10 km depth. Note that this study does not consider the coastal areas, where M4-5 earthquakes have occurred and medium-sized faults (i.e., planar clusters of hypocenters) exist, as the focus here is the influence of lithological properties.
The right panels in Fig. 3 show cross sections of Bouguer anomalies and our best-fit structural model along survey line b-b' (See Additional file 1: Fig. S5 for survey lines a-a' , c-c' , d-d' , and e-e'). The maximum thicknesses of the two high-density zones, which correspond to mafic rocks, are estimated to be ~ 7 km. Ohno et al. (1989) measured the densities of greenschist samples in the Sambagawa belt and obtained an average value of 2.87 g cm -3 , which is slightly higher than the average density of 2.83 g cm -3 reported by Iwaya and Kano (2005). This slightly higher value yields a smaller maximum thickness for the two high-density zones of ~ 4.5 km. However, mafic and pelitic rocks are commonly interlayered in the area of mafic rocks (e.g., Ito et al. 1996). It may, therefore, be appropriate to use a density that is lower than the average density of mafic rocks when estimating the maximum thicknesses of the two high-density zones. Here we assume that the two high-density zones consist of a 60:40 ratio of alternating mafic and pelitic rocks, based on Ito et al. (1996). We determine an average density of 2.76 g cm -3 based on the densities of the mafic and black pelitic rocks obtained by Ohno et al. (1989), which yields an estimated maximum thickness of 20 km. Therefore, the thickness of mafic rocks in the spotted and Mikabu zones is estimated to be between 4.5 and 20 km, based on the possible rock densities and their associated gravity anomalies. We, therefore, consider that the surface rock bodies in this area likely extend to the seismogenic layer. Figures 2 and 4 show that the microearthquake hypocenters occur mainly in the Sambagawa and Chichibu belts. This feature is consistent with the one noted by Maeda et al. (2018), who analyzed the seismicity rate distribution in this region using the JMA catalog for the period from 2005 to 2013. We also observe a low-seismicity zone, where mafic schist and mafic rocks are present in the spotted and Mikabu zones, excluding the area around Mt. Ryumon. These areas of predominantly mafic schist and mafic rock bodies coincide with high Bouguer anomalies (left panel in Fig. 3). These characteristics are also evident in the hypocenter distribution, even before relocation (left panels in Fig. 2). The mafic rocks in the spotted and Mikabu zones, which define the two zones of high gravity anomaly in the Wakayama region, are ~ 7 km thick. We infer that the surface mafic rocks in these zones extend to the seismogenic layer, resulting in the observed low seismicity. The one exception is the Mt. Ryumon area, where the hypocenters are concentrated at depth (Figs. 2 and 4). Furthermore, region b in Fig. 4a, which corresponds to the area of mafic rock, is characterized by a low-seismicity layer with deep upper and lower limits. These results show that the earthquakes are mainly concentrated in the shallow zone of the pelitic rock belts, whereas the seismicity rate is low in the shallow zone but high in the deeper zone of the mafic rock belts in the Wakayama region.

Discussion
Our study results suggest that the seismicity in the Wakayama region is controlled by the lithological properties of the crust. Here we compare the frictional properties of the pelitic and mafic rocks, as obtained from petrological experiments, with the hypocenter depths in the pelitic and mafic lithotectonic belts in the study area to investigate the influence of crustal lithology on the depth extent of seismicity. Petrological experiments have shown that pelitic rocks at 250-400 °C and an effective normal stress of 170 MPa exhibit a velocity-weakening behavior (Hartog et al. 2012). However, basalt, which is a mafic rock, exhibits a velocity-weakening behavior at 300-600 °C, an effective normal stress of 50 MPa, and a pore fluid pressure of 100 MPa (Zhang et al. 2017). These results indicate that pelitic and mafic rocks exhibit a velocity-weakening behavior for different temperature ranges and effective normal stresses. The temperature ranges of 250-400 °C and 300-600 °C correspond to depth ranges of 4.68-7.74 km and 5.70-11.90 km, respectively, when using the heat flow of 129 mW m -2 measured at NKMH station ( Fig. 4a; Matsumoto 2007) and the lithospheric thermal structure estimated by Tanaka (2009). These depths are very close to the depth ranges of the earthquakes analyzed in this study. Figure 5 shows a depth-frequency diagram of the representative hypocenters within each rock zone that are close to NKMH station, thereby comparing the seismicity in the pelitic and mafic rocks with the thermal structure of the crust. The seismogenic layer in each rock unit corresponds to the respective velocity-weakening temperature range obtained from the petrological experiments. Furthermore, the depths of high-seismicity zones coincide with the 300 °C and 400 °C isotherms for the pelitic and mafic rocks, respectively, which is consistent with the experimental results. Abundant hot fluid is expected just beneath the hypocenters in the Wakayama region (e.g., Matsumoto 2007;Yoshida et al. 2011;Omuralieva et al. 2012;Kato et al. 2014;Maeda et al. 2018), and low effective normal stresses that are similar to experimental values are easily achieved at these depths. Therefore, the pelitic rocks at depth in the Wakayama region are expected to experience the temperature and pressure conditions under which velocity-weakening behavior is dominant. The upper and lower limits of the seismogenic layer in pelitic rocks and the upper limit of the seismogenic layer in mafic rocks are consistent with the petrological results. However, the lower limit of the seismogenic layer in mafic rocks does not correspond with the petrological results.
The model of the crustal strength profile, which is another fundamental property of the upper lithosphere, can also explain the dependence of the crustal lithology on the depth extent of the observed seismicity. The crustal strength increases linearly with confining pressure, or depth, in the domain of brittle failure. Furthermore, the strength varies with the minerals in the crust. There are  Fig. 4). The histograms are created using the relocated hypocenters within the rectangles labelled 'a ' and 'b' in Fig. 4a. The light blue and green bars on the right side of the figure indicate the temperature ranges, where the pelitic and mafic rocks show velocity-weakening properties, respectively, as inferred from petrological experiments. The black bars in the background indicate the temperature ranges of the experiments. The light blue and green triangles indicate the temperatures that cause the greatest degree of velocity-weakening (i.e., when the friction rate parameter, a-b, is at a minimum) three main minerals in the rocks of the upper crust (feldspar, quartz, and mica), with feldspar being the strongest and mica being the weakest (e.g., Fig. 1 of Wintsch and Yeh 2013). Albaric et al. (2009) showed that both earthquake magnitude and depth to the brittle-ductile transition are related to the composition of crustal rock, with mafic rocks (diabase and mafic granulite) that consist mainly of mafic minerals tending to possess higher crustal strengths and deeper brittle-ductile transitions. The mafic rocks in the spotted and Mikabu zones contain more mafic minerals than the other zones, which consist mainly of metapelitic rocks that contain more mica and quartz. Therefore, we conclude that the lower limit of the seismogenic layer within mafic rocks is governed by the depth to the brittle-ductile transition.
Our results indicate that earthquakes occur at shallower depths in pelitic rocks and greater depths in mafic rocks as a result of differences in lithological properties (mineral composition, rheology, and rock friction) and the thermal state of the upper crust. The friction parameters made it possible for us to discuss the upper limit as well as the lower limit of the background seismicity, which varies depending on the lithological properties. We demonstrated that it is important to use the friction parameters to understand the influence of crustal lithology and the thermal state on microseismicity.

Conclusion
We investigated the spatial relationship between seismicity and lithological structures in the Wakayama region, southern Honshu, Japan, to elucidate the lithological controls on the occurrence of shallow intraplate earthquakes. We found that the seismicity was high in both shallow pelitic rocks and deep mafic rocks. The observed variations in hypocenter distributions can be explained by the lithofacies and thermal state of these two seismogenic zones.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1186/ s40623-021-01503-3.  Figure S2. (a) Original hypocenter distribution from the JMA unified catalog. (b) Relocated hypocenter distribution using the double-difference algorithm of Waldhauser and Ellsworth (2000) with the JMA2001 seismic velocity model (Ueno et al. 2002). The distribution of seismograph stations used in the relocation is shown in the inset of (a). The two dotted red lines delineate the boundaries between the spotted, non-spotted, and Mikabu zones from north to south. Figure  S3. Bouguer anomalies across the study region. The inset map shows the Bouguer anomalies across the entire study region and the locations of survey lines; the colors and contours are as in Fig. 3. The observed Bouguer anomalies (red circles) and detrended data (yellow circles) are shown for survey lines a-a' to e-e' . Figure S4. Depth distribution of the F-net (NIED) focal mechanism solutions (Mj ≥ 4.0). The presented focal mechanism solutions span the period from 1 January 2000 to 15 March 2021. Figure  S5. Best-fit 2-D gravity models along survey lines a-a' to e-e' . Two panels are shown for each survey line: the upper panels provide a comparison of the detrended observed anomalies (red cross symbols) and calculated anomalies (cyan plus symbols), and the lower panels show the best-fit density structure models that were used to obtain the calculated Bouguer anomalies. Assumed densities of 2.67 g cm -3 and 2.83 g cm -3 are used for the background and green regions, respectively. The inset map shows the observed Bouguer anomaly map for the entire study region and the locations of survey lines.