Identification of a nascent tectonic boundary in the San-in area, southwest Japan, using a 3D S-wave velocity structure obtained by ambient noise surface wave tomography

Abstract

We derived a three-dimensional S-wave velocity model for the San-in area of southwest Japan to examine heterogeneous structures such as tectonic faults. Many earthquakes occur in this area, but much of the activity has been relatively recent, so the fault distribution has yet to be fully clarified. Here, we used continuous ambient noise data from a dense seismic network, deployed from November 2009 to extract Rayleigh and Love wave dispersion data between station pairs, and then applied a direct surface wave inversion to the phase velocities of each station pair to determine a three-dimensional S-wave velocity model. In the resulting model, faults and a previously unrecognized tectonic boundary appeared as low-velocity anomalies or velocity boundaries, and the velocity anomalies were also associated with many past earthquake hypocenters. These results contribute to our understanding of heterogeneous structures caused by recent tectonic motion and of possible future tectonic activity, such as intraplate earthquakes. Surface wave tomography using ambient noise recorded in dense seismic networks could also be applied in other parts of the world to reveal new heterogeneous geological structures (i.e., unrevealed tectonic faults) and could contribute to disaster mitigation.

Introduction

In a convergent plate margin, subsurface structures are complex, and relatively new faults may become activated. However, few geophysical methods are available for identifying unrevealed faults or geologic boundaries over a wide area (i.e., on the scale of the Japanese Islands). For example, active-source seismic reflection and refraction surveys can reveal the distribution of faults at high resolution (e.g., Kamei et al. 2012), but these methods are less suitable for constructing a geological model covering a large region. By contrast, a seismic tomography approach using natural earthquake data can reveal geological boundaries or fault structures in a wide area (e.g., Yolsal-Çevikbilen et al. 2012), but the resolution is usually restricted by the distribution of natural earthquakes. To map previously unidentified faults in the crust and evaluate their seismic activity, a dense seismic network and a geophysical approach for exploring mesoscale heterogeneity (e.g., scale of intraplate faults) need to be developed.

In southwest Japan, the WNW-moving Philippine Sea plate is subducting beneath the Eurasia plate. As part of the oblique convergence of these two plates, several transcurrent fault systems (e.g., the Median Tectonic Line; Fitch 1972) have developed in western Japan (Tsuji et al. 2014). In the San-in area, which is located ~ 350 km north of the Nankai Trough (Fig. 1), strain partitioning has led to the development of tectonic structures in the southern margin of the continental plate that accommodate part of the relative motion of the Philippine Sea plate (Nishimura and Takada 2017). Tectonic structures in this area, such as the Northern Chugoku shear zone (Gutscher and Lallemand 1999) and the Southern Japan Sea fault zone (Itoh et al. 2002), indicate right-lateral strike–slip movement. Using Global Navigation Satellite System (GNSS) data, Nishimura and Takada (2017) identified a zone of concentrated deformation in this area that they called the San-in shear zone (Fig. 1). Although intraplate earthquakes such as the 2000 western Tottori earthquake and microseismic events have occurred along the Japan Sea coast (Gutscher 2001; Gutscher and Lallemand 1999; Kawanishi et al. 2009), no major active fault system has been identified in the shear zone (Headquarters for Earthquake Research Promotion (HERP) 2017). However, the right-lateral strike-slip movement and the presence of conjugate Riedel shears formed in relation to the WNW–ESE maximum compressional stress direction (Kawanishi et al. 2009) suggest that young active faults with strikes of ENE–WSW or NNW–SSE are developing in this shear zone (Fig. 1; Nishimura and Takada 2017; Okada 2002). Although fault features in this area have also been inferred from gravity anomalies (Honda et al. 2002), surface traces of these faults have not been identified (Okada 2002). Therefore, to reveal heterogeneities related to seismic activity, in this study, we estimated high-resolution subsurface structures in the San-in area.

Fig. 1

Topographic map of the San-in area and locations of seismic stations. Inverted red triangles indicate permanent stations (natural frequency 1 Hz), and inverted black triangles indicate temporary Manten project stations (natural frequency 2 Hz). The green line shows the Southern Japan Sea fault zone (SJSFZ), and the red line shows the Northern Chugoku shear zone (NCSZ). The San-in shear zone (SSZ) is shown by white shading. The white and green triangles show the locations of Mt. Daisen and Mt. Sanbe, respectively. Stations A and B are used to compute dispersion curve in Fig. 3

Previous studies, using seismic waves generated by local earthquakes, have applied P- and S-wave tomography to construct seismic velocity models in the San-in area. Zhao et al. (2004) suggested that heterogeneous crustal structures in the 2000 western Tottori earthquake area are related to fluids and arc magma derived from nearby Quaternary volcanoes. Zhao et al. (2018) also suggested that melting of the Philippine Sea slab leads to upward migration of magmatic fluids, affecting both earthquakes and Quaternary volcanic activity in this area. Shibutani et al. 2005, who used dense seismic observations to estimate P- and S-wave velocity structures at high resolution in the area of the 2000 western Tottori earthquake and its aftershocks, identified small heterogeneities near the hypocenters. These heterogeneous structures might have controlled the rupture processes of both the mainshock and aftershocks. However, to identify additional small-scale crustal structures such as young faults in the whole San-in area, it is crucial to estimate seismic velocity structures at higher resolution in a wide area.

In the western Tottori area, many seismic stations have been established to elucidate the distribution of seismic activity and stress state. To detect unidentified heterogeneous structures, therefore, we constructed a high-resolution shallow S-wave velocity model by ambient noise surface wave tomography using data acquired by several seismic arrays, including the dense Manten project array (Iio et al. 2018). Ambient noise surface wave tomography can be used to estimate subsurface heterogeneity even where available earthquakes for body-wave tomography are limited.

Data and methods

We used continuous seismic data recorded by three-component seismometers at 49 permanent stations operated separately by the National Research Institute for Earth Science and Disaster Resilience (NIED), the Japan Meteorological Agency (JMA), and the Disaster Prevention Research Institute of Kyoto University, and by 127 temporary Manten Project stations (Iio et al. 2018) (Fig. 1). Specifically, we used ambient noise data recorded during the month of July 2015, because we found that the maximum number of seismic stations was available in this period. For the 1 month data, we applied the zero-crossing method (Ekström et al. 2009) to estimate Rayleigh and Love wave phase velocities between station pairs. Then, we used the direct surface wave inversion method (Fang et al. 2015) to obtain the three-dimensional (3D) S-wave velocity structure. Each processing step is described in detail below.

Preprocessing and cross-correlation

We divided continuous ambient noise data recorded during each day during July 2015 into 10-min segments with a 50% overlap and corrected the data for the instrumental response. We then computed normalized cross-spectra for each segment as follows:

$$\rho \left( f \right) = \frac{{u_{i} \left( f \right)u_{j}^{*} \left( f \right)}}{{\left| {u_{i} \left( f \right)} \right|\left| {u_{j} \left( f \right)} \right|}} ,$$

(1)

where $\rho$ is the cross-spectrum, $f$ is frequency, $u$ is the seismic record in the frequency domain from station $i$ or $j$, and the asterisk denotes the complex conjugate. To estimate Rayleigh and Love wave dispersion curves, we used vertical cross-spectra from vertical component records and transverse cross-spectra from transverse component records. The transverse cross-spectra were calculated by rotation of the records of the two orthogonal horizontal components. Then, the normalized cross-spectra for data collected during the entire month were stacked. Surface wave propagation between station pairs could be clearly observed in the stacked cross-correlation functions (Fig. 2).

Fig. 2

Cross-correlation functions of a vertical (Z–Z) and b transverse (T–T) components in the 0.2–1.1 Hz frequency band showing empirical Green’s function between station pairs. Each panel shows traces randomly chosen from among those of all station pairs; thus, not all station pairs are included

Surface wave phase velocity measurement

We applied the zero-crossing method (Ekström et al. 2009), which is based on the spatial autocorrelation (SPAC) method (Aki 1957), to estimate phase velocities between station pairs in the frequency domain. With time-domain methods, the interstation distance must generally be larger than three wavelengths to satisfy the high-frequency approximation (Bensen et al. 2007), but with frequency-domain methods the interstation distance can be one wavelength or less (Ekström et al. 2009).

The SPAC method, which assumes that surface waves are dominant in the ambient noise and that ambient noise sources are isotropically distributed, models the real part of vertical cross-spectrum by using the Bessel function of the first kind and zeroth order as follows:

$${\text{Real}}\left( {\rho_{Z} \left( {f,\Delta x} \right)} \right) = J_{0} \left( {\frac{2\pi f\Delta x}{{C_{\text{R}} \left( f \right)}}} \right) ,$$

(2)

where $\rho_{Z}$ is a vertical cross-spectrum, $\Delta x$ is the interstation distance, $J_{0}$ is the Bessel function of the first kind and zeroth order, and $C_{\text{R}} \left( f \right)$ is the Rayleigh wave phase velocity. The Love wave contribution to the real part of transverse cross-spectrum can be expressed as follows (Ekström 2017):

$${\text{Real}}\left( {\rho_{\text{T}} \left( {f,\Delta x} \right)} \right) = \frac{1}{2}J_{0} \left( {\frac{2\pi f\Delta x}{{C_{\text{L}} \left( f \right)}}} \right) - \frac{1}{2}J_{2} \left( {\frac{2\pi f\Delta x}{{C_{\text{L}} \left( f \right)}}} \right) ,$$

(3)

where $\rho_{\text{T}}$ is a transverse cross-spectrum, $J_{2}$ is the Bessel function of the first kind and second order, and $C_{\text{L}} \left( f \right)$ is the Love wave phase velocity. In the real part of transverse cross-spectrum, there is a contribution from Rayleigh waves (Aki 1957) but it would be negligible because we used the station pairs whose interstation distances are longer than two wavelengths, as we describe below. The theoretical background of the zero-crossing method is same as that of the SPAC method, but the zero-crossing method estimates the phase velocity by fitting zeros of the observed cross-spectrum to zeros of the theoretical cross-spectrum described by Eqs. (2) and (3) as follows:

$$C_{m} \left( {f_{n} } \right) = \frac{{2\pi f_{n} \Delta x}}{{Z_{n + 2m} }} ,$$

(4)

where $f_{n}$ is the frequency of the $n$th observed zero on a cross-spectrum, $Z_{n + 2m}$ is the $n$th zero of the theoretical cross-spectrum, and $m \left( { = 0, \pm 1, \pm 2, \ldots } \right)$ is the number of missed or extra zero points. The zeros of the cross-spectrum are used to decrease the effects of spectral attenuation, nonlinear filtering, and background noise (Ekström et al. 2009). To obtain stable zero points on cross-spectra, we applied a moving average, depending on interstation distance, as a smoothing filter (Kästle et al. 2016) (Fig. 3a).

Fig. 3

a Observed vertical cross-spectrum between station A and B (see Fig. 1), which are separated by a distance of 18.986 km. Green dots indicate zero-crossing frequencies. b Phase velocity estimated from the observed cross-spectrum in a by the zero-crossing method. Blue dots show all possible velocities, and red dots show the picked velocities. The black line shows the reference velocity derived by the extended SPAC method, and the green line is the reference dispersion curve obtained from the phase velocity map

Because several phase velocities are estimated at each frequency (blue dots in Fig. 3b), phase velocity is not uniquely determined by only the zero-crossing method. To overcome this problem, we determined a unique phase velocity at the zero-crossing in the low-frequency range closest to a reference velocity. It is easier to determine a unique phase velocity in the low-frequency range because the interval between possible velocities is larger at certain low frequencies (e.g., Fig. 3b) than the interval in the higher frequency range. We computed a reference velocity without lateral variation at 0.25–0.45 Hz (black line in Fig. 3a) by the extended SPAC method (Ling and Okada 1993) using the cross-spectra of all station pairs. Because of lateral velocity variation (i.e., heterogeneity) in the study area, however, it is difficult to use a reference velocity to determine a unique phase velocity. Since phase velocities between station pairs with a long interstation distance (e.g., more than 80 km) are not much affected by localized velocity variations, a unique phase velocity can be picked for those pairs from a reference velocity without lateral variation. Thus, after estimating a few continuous phase velocities from the first picked phase velocity, we constructed a two-dimensional (2D) reference phase velocity map in the 0.25–0.45 Hz frequency range by applying surface wave tomography to phase velocities between station pairs with an interstation distance of more than 80 km.

For the surface wave tomography, we conducted a least-squares inversion that included a small amount of damping to minimize the roughness of the slowness variation by assuming straight ray paths of surface wave propagation between two stations (Ekström 2014). The reference phase velocity between each station pair was then calculated by computing the path-averaged slowness from the reference map. The resulting reference phase velocity (green line in Fig. 3b) differs for each station pair. We used the reference phase velocity of each station pair to choose the phase velocity dispersion curve between two stations (red dots in Fig. 3b).

At high frequencies, cross-spectra can be modeled by a cosine function. Therefore, we used the following zero-crossing frequency step criterion in picking phase velocities (Kästle et al. 2016). To choose the next phase velocity at frequency $f_{k + 1}$, the zero-crossing frequency step (${\text{df}}$) was calculated as follows:

$${\text{df}} = f_{k + 1} - f_{k} = \frac{{C_{\text{ave}} }}{2\Delta x} ,$$

(5)

where $C_{\text{ave}}$ is the average reference phase velocity computed by averaging reference velocities with and without lateral variation. If the frequency step was much smaller or larger than ${\text{df}}$, the phase velocity $C\left( {f_{k + 1} } \right)$ was skipped and we chose the next phase velocity $C\left( {f_{k + 2} } \right)$ in our algorithms. When the unstable frequency step began to repeat, we stopped picking phase velocities. After measurement of all station pairs, we checked each dispersion curve manually and removed the dispersion curves of station pairs that had not been picked accurately. The phase velocity dispersion curves of Rayleigh and Love waves from 0.2 to 1.0 Hz are shown in Fig. 4a, b, respectively.

Fig. 4

Surface wave dispersion curves: a Rayleigh wave phase velocities. b Love wave phase velocities. In each panel, black curves show the phase velocities at each frequency, and the red curve shows the number of measurements at each frequency

Direct inversion of the surface wave dispersion curve

To determine the 3D S-wave velocity structure, we directly inverted the surface wave phase velocity using DSurfTomo (Fang et al. 2015; Li et al. 2016). This approach avoids the intermediate process of constructing a 2D phase velocity map by directly inverting the surface wave travel-time data. It uses the fast-marching method (Rawlinson and Sambridge 2004) to take account of the ray-path-bending effects of surface waves caused by a complex velocity structure.

Following Fang et al. (2015), the travel-time difference $\delta t_{i}$ at frequency $f$ along path $i$ is given by

$$\delta t_{i} \left( f \right) = t_{i}^{\text{obs}} \left( f \right) - t_{i} \left( f \right) \approx - \mathop \sum \limits_{k = 1}^{K} \nu_{ik} \frac{{\delta C_{k} \left( f \right)}}{{C_{k}^{2} \left( f \right)}} ,$$

(6)

where $t_{i}^{\text{obs}} \left( f \right)$ is the observed surface wave travel time, $t_{i} \left( f \right)$ is the travel time calculated from a reference model that can be updated in the inversion, $K$ is the number of grid points in the horizontal direction, $\nu_{ik}$ is bilinear interpolation coefficients along the ray path associated with the $i$th travel-time data, $C_{k} \left( f \right)$ is the phase velocity, and $\delta C_{k} \left( f \right)$ is its perturbation at the $k$th 2D surface grid point at frequency $f$. By using one-dimensional (1D) depth kernels of Rayleigh and Love wave phase velocity data with respect to the P-wave velocity $\alpha$, the S-wave velocity $\beta$, and density $\rho$ at each surface grid node, Eq. (6) can be transformed as follows:

$$\begin{aligned} \delta t_{i} \left( f \right) = \mathop \sum \limits_{k = 1}^{K} \left( { - \frac{{\nu_{ik} }}{{C_{k}^{2} }}} \right)\mathop \sum \limits_{j = 1}^{J} \left. {\left[ {R^{\prime}_{\alpha } \left( {z_{j} } \right)\frac{{\partial C_{k} }}{{\partial \alpha_{k} \left( {z_{j} } \right)}} + R^{\prime}_{\rho } \left( {z_{j} } \right)\frac{{\partial C_{k} }}{{\partial \rho_{k} \left( {z_{j} } \right)}} + \frac{{\partial C_{k} }}{{\partial \beta_{k} \left( {z_{j} } \right)}}} \right]} \right|_{{\theta_{k} }} \delta \beta_{k} \left( {z_{j} } \right) = \mathop \sum \limits_{l = 1}^{M} G_{il} m_{l} , \\ \end{aligned}$$

(7)

where $\theta_{k}$ represents the 1D reference model at the $k$ th surface grid point, and $\alpha_{k} \left( {z_{j} } \right)$, $\beta_{k} \left( {z_{j} } \right)$, and $\rho_{k} \left( {z_{j} } \right)$ are the P-wave velocity, S-wave velocity, and mass density, respectively. $J$ is the number of grid points in the depth direction, and the total number of grid points in the 3D model is $M = KJ$. This direct inversion uses an empirical relationship based on Brocher (2005) to relate perturbations of P-wave velocity and density to perturbations of S-wave velocity with scaling factors $R^{\prime}_{\alpha }$ and $R^{\prime}_{\rho }$ (see Fang et al. 2015 for details). Equation (7) can also be formulated in matrix form as

$${\mathbf{d}} = {\mathbf{Gm}},$$

(8)

where ${\mathbf{d}}$ is the residual vector of the surface wave travel time for all ray paths at different frequencies, ${\mathbf{G}}$ is the data sensitivity matrix, and ${\mathbf{m}}$ is the model parameter vector. An estimate of the solution ${\mathbf{m}}$ can be obtained in the same way as for a classical inversion function [Eq. (8)] by using the LSQR (least squares with QR-factorization) algorithm (Paige and Saunders 1982) and regularized by applying both damping and first-order smoothing. We selected the values of both parameters by trial and error, by considering the patterns of the heterogeneities in inverted S-wave velocity models.

To obtain the 3D S-wave velocity structure by direct inversion, we constructed an initial S-wave velocity model (Fig. 5) that was calculated by multiplying the measured Rayleigh wave phase velocity at a depth of one-third the wavelength by 1.1 (i.e., a one-third wavelength transformation) (Fang et al. 2015). We defined the initial S-wave velocities by averaging converted phase velocity within ± 200 m at the depth of each grid (e.g., between the red lines in Fig. 5). The initial S-wave velocities for depths shallower or deeper than the investigation depth derived from the one-third wavelength transformation were obtained by extrapolating the linear trend close to the boundaries. We selected the surface wave travel-time data in each period according to distance and wavelength limits such that $\Delta x_{ \text{min} } = 2\lambda$ and $\Delta x_{ \text{max} } = 5\lambda$, where $\Delta x_{ \text{min} }$ is the minimum interstation distance, $\Delta x_{ \text{max} }$ is the maximum interstation distance, and $\lambda$ is the wavelength, because measurements of short paths have large relative errors and uncertainties in the measurements of long paths during short periods are also large (Ekström 2014). The grid size of our model was 0.025° by 0.025° laterally, with vertical levels at 0, 0.1, 0.2, 0.4, 0.8, 1.2, 1.8, 2.4, 3.2, 4.0, 5.0, 6.0, and 7.0 km depth. To consider topographic effects on our inverted velocity model, we subtracted altitude from depth at each grid point and show the resulting S-wave velocity structure below sea level.

Fig. 5

Initial S-wave velocity model. Black dots show all transformed phase velocity dispersion curves on a depth–S-wave velocity approximation profile. The blue line and dots show the average S-wave velocity model for this profile, which was used as the initial S-wave velocity model in the inversion. The red lines indicate, as an example, the depth range used to compute the average S-wave velocity at 2.4 km depth

Results

We built a 3D S-wave velocity model by direct inversion of Rayleigh and Love wave dispersion curves separately because of the possibility of radial anisotropy (Fig. 6). Figure 7 shows the path coverage of both Rayleigh and Love waves based on the final inverted S-wave velocity models (Fig. 6). Although ray-path coverage decreases with an increase in frequency in the whole study area, the spatial coverage of ray paths remains good in the area with dense station coverage even at the higher frequency. The major velocity anomalies that we discuss in the following can be observed by changing the smoothing and damping parameters (Additional file 1: Figures S2 and S3). To confirm the resolution of our inverted S-wave velocity model, we conducted a checkerboard resolution test using anomaly sizes of ~ 0.2° (Fig. 8a–d) and ~ 0.1° (Fig. 8e–h). The amplitude of the velocity anomaly was ~ 10%. Both Rayleigh and Love wave results clearly resolved the larger velocity anomalies (Fig. 8a–d). However, the smaller velocity anomalies were retrieved only in the northern part of the study area, where the density of seismic stations was high (Fig. 8e–h). At shallower depth, Love wave results have higher accuracy than Rayleigh wave results, while Rayleigh wave results have higher accuracy at deeper depth (Fig. 8). This trend can be explained by the different depth sensitivities of Rayleigh and Love waves (Fig. 9).

Fig. 6

Absolute S-wave velocity values and their perturbation from the velocity average in lateral S-wave velocity slices at 0.5 km and 3 km depth. a, c, e, and g Rayleigh wave results. b, d, f, and h Love wave results. In each panel, black lines show prefecture borders, and depth is shown above the upper left corner. a White and red triangles show the locations of Mt. Daisen and Mt. Sanbe, respectively, the red line shows the Kamakurayama south fault (KSF), the black bold line shows the 2000 western Tottori earthquake fault, and blue lines show other previously recognized faults in this area. c The locations of sections X–X′ and Y–Y′, shown in Fig. 10, are indicated. f The dashed black rectangle indicates a lineament through Shimane to eastern Hiroshima. The results without topography correction are shown in Additional file 1: Figure S1

Fig. 7

Ray-path coverage for two selected periods of Rayleigh waves (a and c) and Love waves (b and d). In each panel, frequency is shown above the upper left corner. Inverted black triangles indicate seismic stations, and blue lines indicate ray paths

Fig. 8

Checkerboard resolution test for Rayleigh waves (a, c, e, and g) and Love waves (b, d, f, and h). Depth is shown above the upper left corner of each panel. The anomaly size was ~ 0.2° (a–d) and ~ 0.1° (e–h), and the velocity anomaly amplitude was ~ 10%

Fig. 9

Depth sensitivity kernels for phase velocity of a Rayleigh waves and b Love waves at three different frequencies, computed using the code from Herrmann (2013)

Low-velocity anomalies can be observed in results derived from both Rayleigh and Love waves (Fig. 6). The low-velocity zone observed in the northern part of the study area along the coast (Shimane peninsula) was also observed in the S-wave velocity structure estimated by S-wave tomography using data of natural earthquakes by Zhao et al. (2004), as well as in the S-wave velocity structure constructed from ambient noise by Nishida et al. (2008). The lateral resolution of our velocity model is higher than that of Nishida et al. (2008), however, because we used shorter wavelength surface waves estimated by using data from a denser seismic array in our tomographic analysis. Another characteristic structure in our model is a low-velocity lineament with an ENE–WSW strike at relatively shallow depth in the eastern part of the San-in shear zone (around the red line in Fig. 6a). Furthermore, in the S-wave velocity model derived from Love waves, we identified a continuous low-velocity zone oriented NW–SE and extending from the Sea of Japan to eastern Hiroshima (black dashed rectangle labeled A in Fig. 6f).

The 3D S-wave velocity distributions estimated by using Rayleigh waves (Fig. 6a, c, e, g) and Love waves (Fig. 6b, d, f, h) are similar. The slight differences may be caused in part by the different depth resolutions of Rayleigh and Love waves, which are determined by their depth sensitivities (Fig. 9), and by the radial anisotropy of S-waves (e.g., Dreiling et al. 2018; Luo et al. 2013; Ojo et al. 2017).

Interpretation

The low-velocity zone along the coast in the northern San-in area (i.e., around the Shimane Peninsula and Mt. Daisen) (Fig. 6) is attributable to the geologically young rocks in that region; Mt. Daisen is a Quaternary volcano, and the Shimane Peninsula is composed of Neogene and Quaternary rocks (Geological Survey of Japan, AIST 2015). The velocity boundary on the south side of this low-velocity area might be related to a geological boundary generated by the opening of the Sea of Japan (Geological Survey of Japan, AIST 2015) or to the tectonic boundary proposed by Gutscher and Lallemand (1999) and Nishimura et al. (2018). This velocity boundary has a northward dip, and beneath the Shimane Peninsula the low-velocity anomalies appear to have a synclinal structure (Fig. 10). Sawada et al. (2001) observed this dip and synclinal structure in seismic and borehole data obtained near our vertical X–X′ section (Fig. 10a, b). In the eastern Shimane Peninsula, Yamauchi and Iwata (1998) also delineated a dipping velocity boundary and a synclinal structure of geologically young rocks, based on geological and geophysical data obtained in the area around our Y–Y′ vertical section (Fig. 10c, d). Thus, our S-wave velocity model revealed the 3D geometry of the geology and a tectonic boundary in this area.

Fig. 10

Vertical slices through the S-wave velocity structure beneath the Shimane Peninsula. a, c Rayleigh wave results. b, d Love wave results. The locations of these slices are shown in Fig. 6c

Although velocity tended to be higher in the southern part of the study area, we observed some low-velocity anomalies, including lineaments. One of the dominant anomalies is observed around the Kamakurayama south fault (KSF; red line in Fig. 6a). The fault was identified by a low S-wave velocity anomaly at relatively shallow depth (~ 0.5 km; Fig. 6a, b, e, f), maybe because it has been recently active (e.g., Nishimura and Takada 2017). The low-velocity anomaly clearly identified in our results, however, is longer than the previously identified fault (Research Group for Active Faults of Japan 1991), suggesting a continuation of the fault. The KSF might have been generated by right-lateral strike-slip movement with an ENE–WSW strike. Fault systems such as this one revealed by our results are common features of the San-in shear zone (Nishimura and Takada 2017). The detection of the KSF as a low-velocity anomaly could be partly attributable to the dense seismic station network around the fault (Fig. 1).

The wide low-velocity area observed from Shimane to eastern Hiroshima (black dashed rectangle A in Fig. 6f) corresponds to sedimentary rocks in the Miyoshi Basin (Geological Survey of Japan, AIST 2015; Fig. 6b). However, we showed that a clear low-velocity anomaly extends toward the NW from eastern Hiroshima to the Sea of Japan coast. This low-velocity lineament includes an active volcano (Mt. Sanbe; red triangle in Fig. 6a) (Geological Survey of Japan, AIST 2015). Furthermore, the structure is more visible in the tomographic image at 3 km depth, compared to that at 0.5 km depth (Fig. 6). To evaluate seismic activity associated with the distribution of heterogeneities identified here, we plotted hypocenters (the JMA catalog) on the seismic velocity perturbation maps, and the result showed that the low-velocity anomalies are well consistent with hypocenter distributions (Fig. 11). In particular, many small earthquakes have occurred along the low-velocity zone from eastern Hiroshima to Mt. Sanbe (Fig. 11a, b). Especially earthquakes intensively occurred at ~ 30 km southeast of Mt. Sanbe where the low-velocity lineament is discontinuous.

Fig. 11

Hypocenters of earthquakes shallower than 5 km depth (white dots) superimposed on the S-wave velocity structure. a, c Rayleigh wave results. b, d Love wave results. Depth is shown above the upper left corner. Hypocenters of earthquakes with moment magnitude > 1 occurring during 1983–2016 are shown. Green triangle shows the Mt. Sanbe

Because the maximum compressional stress direction in the San-in shear zone is WNW–ESE (white-shaded area in Fig. 1; Nishimura and Takada 2017; Okada 2002), conjugate Riedel shears in this area may be oriented NW–SE, consistent with the orientation of this low-velocity lineament. Therefore, this zone may be related to a nascent microplate boundary associated with conjugate Riedel shears in the San-in shear zone (Nishimura and Takada 2017; Watanabe 2004). Because we identified the nascent tectonic boundary using ambient noise surface wave tomography, this approach could be used to reveal intraplate fault systems over wide areas and possibly contribute to disaster mitigation.

Conclusion

To clarify small-scale heterogeneous structures (~ 10 km) in the San-in area (> 100 km), we constructed a high-resolution 3D S-wave velocity structure by surface wave tomography of dense seismic network data. A fault striking ENE–WSW, a microplate boundary in the western San-in area, and a geological or tectonic boundary on the northern side of the area appear in our S-wave velocity model as low-velocity anomalies or velocity boundaries. These features were related to right-lateral strike-slip movement and conjugate Riedel shears in the shear concentration zone, and they agree with hypocenter distributions. Thus, our high-resolution S-wave velocity model contributes to our understanding of tectonic and geological features in the San-in area.

The resolution of the derived S-wave velocity model is suitable for identifying geological boundaries in a wide area. Thus, our approach could also be applied to ambient noise recorded by dense seismic networks in other areas to identify faults and geological and tectonic boundaries.