Shear wave splitting and seismic velocity structure in the focal area of the earthquake swarm and their relation with earthquake swarm activity in the Noto Peninsula, central Japan

Seismic activity in the Noto region of Ishikawa Prefecture, central Japan, has increased since August 2020 and has continued as of August 2023. Stress changes due to subsurface sources and increases in fluid pressure have been discussed as the causes of the seismic activity increase. In this study, S-wave polarization anisotropy was investigated by S-wave splitting analysis using temporary and permanent stations located in the epicenter area. We also investigated the seismic wave velocity structure in the source region by analyzing seismic wave velocity tomography. The fast orientations of anisotropy (fast shear wave oscillation direction, FSOD) were generally NW–SE in the southern part of the focal area and east–west in the northern part. The NW–SE anisotropy generally coincides with the direction of the maximum horizontal compression axis, both near the surface and at earthquake depths. Therefore, stress-induced anisotropy can be the cause of the observed NW–SE anisotropy. On the other hand, faults with strike directions generally east–west have been identified, and structural anisotropy may be the cause of the observed east–west anisotropy. We examined the time variation of anisotropy at N.SUZH, one of the permanent stations. No significant time variation was observed in the FSOD. Larger anisotropy was observed, particularly for the activity in the western part of the focal area, from about June–September 2021 compared to the previous period. A high Vp/ Vs region was identified beneath the focal area, at a depth of 18 km. This high Vp/Vs region has slightly larger P-wave velocities than the surrounding area. Since Tertiary igneous rocks are distributed in the target area, the high Vp/Vs region may represent a Tertiary magma reservoir, suggesting that fluids released through the old magma reservoir are involved in this seismic swarm. This seismic activity started in the southern part of the area, where relatively immature fault structure exists, where stress-induced anisotropy is distributed, and where high Vp/Vs regions suggestive of fluid at depth are identified. Subsequently, seismicity became more active in the northern part, where structural anisotropy with well-developed fault structures is distributed.


Seismic velocity anisotropy and shear wave splitting
Seismic velocity anisotropy in the crust was first discovered by Anderson (1961) as a propagation azimuth anisotropy, using the principle of seismic wave propagation in an anisotropic medium.S-wave splitting, a form of polarization anisotropy, was first observed by Nur and Simmons (1969).Later, Ando et al. (1980) and Crampin and Booth (1985) examined it using earthquake sources in the crust, and it has been studied in various ways since then.Anisotropy in the crust that can cause S-wave splitting has been considered in multiple origins, including stress-induced (e.g., Crampin 1987;Crampin and Zatsepin 1997;Okada et al. 1994;Nur 1971; Zatsepin and Crampin 1997), fault-induced (e.g., Gledhill 1991;Savage et al. 1990;Boness and Zoback 2006;Li and Peng 2017), and structural (e.g., Zatsepin and Crampin 1997;Nur 1971;Boness and Zoback 2006;Li et al. 2014;Menke et al. 1994;Zinke and Zoback 2000), and mineralinduced due to preferential orientation of mineral crystal lattices by ductile deformation (e.g., Alford 1986;Aster and Shearer 1992;Sayers 1994).
In the earth's crust, microfractures are usually randomly oriented.When differential stresses are applied in the crust, microcracks oriented in directions other than the direction of maximum horizontal compressive stress (SHmax) are closed by compression of these cracks.This closure of the microcracks causes anisotropy in the direction of SHmax since only microcracks oriented in the direction of SHmax will remain open.Furthermore, only microcracks oriented toward SHmax can be newly formed, resulting in anisotropy in that direction.We call this anisotropy stress-induced anisotropy.
On the other hand, structure-induced anisotropy is caused by the orientation of microfractures parallel to the fault generated by fault movement (earthquakes) and layered minerals such as mica or cataclasite in the fault zone, which follow the strike of the fault, resulting in structure-induced anisotropy in the direction of the fault strike.By studying these anisotropies, we can obtain information on the crust and mantle's structure, composition, deformation, and stress state (Savage 1999).
In some studies, temporal change of shear wave splitting was reported.Mroczek et al. (2020) found that anisotropy changed in response to changes in pore pressure caused by water injection and production during geothermal power generation in the Taupo Volcanic Zone, North Island, New Zealand.Miller and Savage (2001), Gerst and Savage (2004), and Bianco et al. (2006) observed temporal changes in anisotropy related to volcanic eruptions in New Zealand and Italy.Illsley-Kemp et al. (2018) measured anisotropy before, during, and after volcanic eruptions in the Afar Basin, Ethiopia, as did Savage et al. (2015) in La Reunion, and they found that anisotropy changed in response to magma movement, suggesting that anisotropy can be used to monitor and predict magmatic activity at active volcanoes.Hiramatsu et al. (2005) found a temporal change in the delay time of shear wave splitting after the Mw 5.7 earthquake in and around the Philippine Sea plate subduction zone.They suggested that the crack opening or/and enlargement by the stress change due to the quake and the healing of the crack following it caused this temporal change in anisotropy.
Seismic anisotropy changes with time are often difficult to distinguish between changes in anisotropy with space combined with changes in earthquake locations, since often earthquakes move in response to stress changes, causing a change in the earthquake-station paths.Therefore, it is not always clear whether a spatial change in

Graphical Abstract
anisotropy is the cause of observed apparent temporal changes.(e.g., Aster et al. 1990).
The 2023 M6.5 earthquake and preceding swarm-like seismic activity in the Noto Peninsula, central Japan Seismic activity in the Noto region of Ishikawa Prefecture has increased since August 2020, has continued until August 2023 (Additional file 1: Figure S1), and still ongoing as of this writing of November 2023.The swarm consists of about four clusters.At first, it began at the southern cluster and extended to the western cluster in early 2021.The northwestern cluster started its activity in mid-2021 and spread to the northeastern cluster.The earthquake of M5.4 occurred on June 19, 2022, in the less active area between the northwestern and northeastern clusters.The largest earthquake of M6.5 occurred on May 5, 2023, near the M5.4 earthquake.This area is located at the southeastern margin of the Japan Sea, where extensional basin formation occurred in the Miocene and currently, compressional shorting with reverse faulting is ongoing (Additional file 1: Fig. S2; e.g., Sato 1994, Sibson 2009).No active volcano exists now.Nakajima (2022) used arrival time data at the routinely operated stations and performed seismic velocity tomography in and around the Noto Peninsula.He found a low-Vs area beneath the focal area of the swarm and suggested that fluid promoted the occurrence of this swarm activity.Nishimura et al. (2023) deployed the temporary GNSS stations and modeled the deformation source related to this swarm activity.One of the possible models is a tensile-shear crack in the hypocentral area, which represents an aseismic slip, promoting the earthquake swarm activity (e.g., Guglielmiet al. 2015).
These are some studies on the stress field in and around the Noto peninsula.Terakawa and Matsu'ura (2010) estimated the stress field by the inversion of moment tensor solutions from the National Institute for Earth Science and Disaster Mitigation.The reverse fault type with a SHmax oriented in the WNW-ESE direction is estimated in the northeastern Noto peninsula.
A geological map (Yoshikawa et al. 2002;Sawada et al. 2012;Fig. 1;Geological Survey of Japan, AIST, 2022) shows that most of the study area is characterized by Tertiary rocks.The sedimentary rocks are along the northern coastline, and the volcanic rocks are in the southern part of the area.Rhyolitic rocks are predominant in the volcanic rocks, but basaltic rocks are distributed in some limited areas.These rocks are distributed as elongated in the east-west direction.This distribution is because the east-west or ENE-WSW striking normal faults developed as the rift stage with north-south or NNW-SSE extension in the Tertiary.Some normal faults are activated as reverse faults called compressional inversion faults (e.g., Sato 1994;Okamura et al. 2019).

Temporary seismic observation and purpose of this study
There were the only two permanently operated stations (N.SUZH and SUZU) near the swarm focal area (Fig. 2).The average spatial distance between permanent stations in and around the focal area is about 20 km, so obtaining the structure on a scale smaller than 20 km is difficult.To get a detailed hypocenter distribution and crustal structure, Tohoku University and the Earthquake Research Institute of the University of Tokyo deployed four temporary stations (Sakai et al. 2022, Additional file 1: Table S1).Each station consists of a short-period three-component seismometer and a mobile-phone telemetry system.The seismographs were short-period seismographs (natural frequency: 1 Hz or 2 Hz).The sampling frequency was 100 Hz.This network began at the end of June and ended at the end of December 2022.
In the case of the Noto Peninsula earthquake swarm, stress changes due to subsurface sources and increases in fluid pressure are discussed as the causes of the seismic activity increase (e.g., Amezawa et al. 2023;Nishimura et al. 2023;Yoshida et al. 2023aYoshida et al. , 2023b)).Fluid pressure changes and stress changes could cause changes in shear wave splitting.In this study, S-wave polarization anisotropy was investigated by S-wave splitting analysis using temporary and permanent stations located in the epicenter area.Also, a fine structure of the anisotropic heterogeneity obtained from shear wave splitting would help us understand this earthquake swarm's process.
We also investigated the seismic wave velocity structure in the source region by analyzing seismic wave velocity tomography.There are some previous studies on seismic velocity tomography in the study area (e.g., Nakajima 2022).But the spatial resolution is about 20 km, which is the same as the average distance between stations, and it is difficult to know the detailed structure in the swarm focal area on a scale of about 10 km.For example, there is a low gravity (Bouguer) anomaly with a scale of 10 km in the focal region (Sawada et al. 2012).To understand the origin of this gravity anomaly, the seismic velocity structure on a comparable scale is necessary.In this study, we used travel time data from the temporary stations in the focal area to obtain a finer structure and its relation with the swarm earthquakes, gravity, and seismic anisotropy.

Data and method
We determined the seismic velocity structure using data from the six temporary and routinely operated stations by adopting the double-difference tomography method (Zhang andThurber 2003, 2006) to obtain the seismic velocity structure and earthquake locations simultaneously.The grid interval is 0.05 degrees horizontally and 6 km in depth.
We start the analysis with 1648 earthquakes from the Japan Meteorological Agency (JMA) with depths less than 40 km and distances between station and earthquake pair of less than about 1 degree from January 2018 to August 2022, when one of the temporary stations ended.This allows us to use all the stations in the focal area to obtain a stable seismic velocity structure (e.g., Koulakov 2009).Initial hypocenter locations, structure, and arrival time data are from the JMA catalog and their one-dimensional velocity model (Ueno et al. 2002), in addition to the arrival times at our temporary stations determined by the automated arrival time picking system of Horiuchi et al. (1992).Figure 2 shows the final hypocenter distribution and station map.The root-mean-square of the arrival time residual is decreased from 0.248 s to 0.061 s for the sum of both P and S waves.Next, using the seismic velocity structure obtained above, we relocated all the 11049 hypocenters between January 2018 and May 8, 2023.The root-mean-square of the arrival time residual of summing of both P-and S-waves is decreased from 0.225 s to 0.059 s.
We estimated the spatial resolution of the obtained model as one grid (6 km) in the area of the swarm's  Sawada et al. (2011).Red thick and thin lines denote the active and non-active faults, respectively.The rectangle denotes the spatial range of Fig. 5b.The ellipse with broken line shows the approximate location of the inherited caldera suggested by the gravity anomaly (see Fig. 3a).The moment tensor solution (NIED, 2023) of the M6.5 earthquake is shown on the right hypocenters to two grids (12 km) in the surrounding area by the checkerboard resolution test (CRT) (Additional file 1: Fig. S3).The hypocenter's location error, i.e. the difference between the true and obtained locations estimated through the CRT, is about 100 m horizontally and 280 m vertically.We also did a reconstruction resolution test and the image was well reconstructed in and around the hypocentral area except shallower than 6 km(Additional file 1: Fig. S4).
This seismic velocity tomography assumes isotropic velocity, which can be obtained when there is good recovery of the checkerboard and reconstructed pattern (c.f., Zhao et al. 2016).

Figure 3b
, c shows the horizontal slice of Vp at depths of 6 km and 18 km, respectively, for the final velocity model.Vp at 6 km has a low-velocity area of about 5.6 km in the central and southern part of the study area.This low-Vp region spatially corresponds to the low Bouguer gravity anomaly (Sawada et al. 2012) in Fig. 3a (inside the ellipse with the red dotted line).Thus this low Vp region is composed of low-density rocks.The surface geology of this area is Middle-Miocene to Late Pleistocene sedimentary rock (Fig. 1).Vp at 18 km has a high-velocity area just below the shallow low-Vp area at 6 km depth.The swarm of hypocenters are located along the northern boundary of this mid-crustal high-velocity area.

Data and methods
We used the MFAST automatic shear wave splitting calculation system of Savage et al. (2010).The code is based on Silver and Chan (1991,); MFAST performs a cluster analysis (Teanby et al. 2004) on the shear-wave splitting results as a function of time windows of various starting times and lengths.In MFAST, the splitting measurements are done for different time windows with the frequency filter that returns the best signal-to-noise ratio.Cluster analysis of different window measurements is used to select the final measurement and estimate its uncertainty automatically.Measurements are graded from A (best) to D (worst) depending on the consistency between the measurements in the different windows (Savage et al.We used MFAST on waveform data from permanent and temporary stations in the epicenter area from June 2022 to May 8, 2023.Waveforms with an incident angle greater than 35 degrees at the surface are not used in the analysis because of phase shifts of the waveform due to S-to-P conversion at the free surface.The incident angles were calculated with a velocity structure (Ueno et al. 2002) that had 2.84 km/s of S-wave velocity at the surface.We analyzed 24666 waveforms (event-station pairs) of 4742 events with grade A.

Results
Additional file 1: Figure S5 shows an example of a shear wave splitting measurement.Figure 5a shows a rose diagram for each station's fast shear wave oscillation direction (FSOD, φ).The dominant direction is east-west for northern stations, north-south in the southeast, WNW-ESE in the south, and NE-SW in the west.
We generated spatially averaged φ using TESSA (Johnson et al. 2011) (Fig. 5), weighting the φ values by the inverse of the distance (from the station to the grid cell) squared and assigning them to each grid block that each ray passes through.We plot the rose diagram of the φ values in grid and plot the mean direction (computed using circular statistics) from each grid block only if the standard deviation of the data is less than 30° and the standard error of the mean is less than 10°.These criteria allow us to exclude blocks that exhibit large scatter or multiple modes.In this case, we used a grid size of 1500 m.We did not use cells with fewer than 10 rays passing through.
From this spatially averaged map (Fig. 5), the anisotropy orientation is spatially heterogeneous, as expected from the individual station rose diagrams.In the northern part of the area, the anisotropy orientations are mostly oriented E-W.In the southern part of the area, the anisotropy orientations swing from NNE-SSW in the east to NW-SE in the center and NE-SW in the west.
This spatial change of the anisotropy orientation correlates with the geology.In the northern part of the area, i.e., the north coastal region, geologic features are elongated E-W (Fig. 1).This elongated geological pattern is related to the fault or the fold formation.These E-W The background grey boxes represent the shaded topography elongated fold and fault-related structures cause E-W fast seismic anisotropy in the macroscale if the size of the structure is enough smaller than the wavelength of a few kilometers to a few hundred meters.Thus the E-W-oriented anisotropy in the northern part of the area can be interpreted as structure-controlled anisotropy.
In the southern part of the area, the anisotropy orientations are mostly oriented NW-SE.This orientation is consistent with the SHmax orientation in this area (e.g., Terakawa and Matsu'ura 2010).Thus the NW-SE-oriented anisotropy in the southern part of the area can be interpreted as stress-induced anisotropy.
This difference between causes of anisotropy could be caused by variable deformation forming faults and folds.In the northern area, faults and folds were well-developed.In and around the caldera in the southern part, it might be difficult to form large-scale faults and folds.
The origin of the N-S trend of the first direction at SUZU in the southeastern part is unclear.One of the possible origins is the structural controlled anisotropy.In the study area, there are some short N-S striking faults along the southern coast of the Peninsula (see Fig. 1).If these N-S striking fault exist near SUZU, they might be a source of structural controlled anisotropy, which is oriented N-S.In future studies, estimating the anisotropy on a finer scale is necessary.
Previous studies discussed the temporal change of anisotropy related to seismic activity (e.g., Li and Peng 2017).Therefore, we examined the anisotropy as a function of time from 2020 to May 8, 2023, using data (1530 stationevent pair) from station N.SUZH, which is located in the swarm area.
Additional file 1: Figure S7 shows the temporal patterns of delay time, station-event distance, and φ.No significant temporal change in φ could be confirmed.We observed large delay times after June 2021 (2021.5).If the anisotropy is homogeneously distributed, the delay time may increase as the distance between the epicenter and the station increases.The relationship between time and the distance between the hypocenter and the station was shown as a reference.We found no change in the delay time proportional to the distance between the epicenter and the station (see also Additional file 1: Figure S9).
However, since the epicenters in this earthquake swarm are widely distributed in both azimuth and distance from the station, this temporal variation may be influenced by spatial sampling.The moving window median value of the delay time normalized by the length of ray path did not show a significant change (Additional file 1: Fig. S8).Therefore, we created Additional file 1: Figure S9, which plots the anisotropy parameters, i.e., the direction of anisotropy and the delay time, at the epicenter location for five periods.For all five time-windows, the relationship between the delay time and the hypocentral depth shows no distinct increase in delay time with depth, but the delay times are spatially varied.This suggests some possibilities that 1) the origin of anisotropy is distributed shallower than the focal depth; 2) the origin is localized at hypocenters, although they could not be distinguished due to the limited depth distribution of the hypocenters used in this study (see also Additional file 1: Fig. S10, in which delay times and fast directions as a function of depth for all are shown).Since seismic activity has been observed in the northern part of the focal area frequently after the latter half of 2021, some large delay times have been observed locally from that period.Thus, we considered that small regions of high anisotropy exist on the northern side of the epicenter area, and that time variation in anisotropy cannot be confirmed.

Seismic velocity structure
We discuss the lithological interpretation of this high-Vp and high Vp/Vs anomalous area beneath the low-Vp, lowgravity anomaly, and southern subgroup of the earthquake swarm (Fig. 3).The values of Vp, and Vs in the upper to the middle crust at a depth of 6-18 km are in the range of about 5.4-6.4km/s, in a range of 3.2-3.4km/s, respectively (Figs. 3 and 4).These values are mostly consistent with granite-gneiss (Christenten, 1996).The surface geology suggests the existence of rhyolitic rock in a wider region within the focal area (Fig. 1, Sawada et al. 2012).Granite or granodiorite has a Vp of less than 6.3 km/s (Christenten 1996), so the average velocities are consistent with the surface rocks.However, it is difficult to explain the values of Vp of about 6.4-6.6 km/s in the anomalous area with granite.There are basaltic-andesite rocks in the western area of the focal area.Diabase has about Vp of about 6.78 km/s, Vs of about 3.76 km/s, Vp/Vs of about 1.80 at 600 MPa (Christenten, 1996).If overpressured fluid causes a decrease of Vp, Vs, and an increase of Vp/Vs, the obtained values of Vp, Vs, and Vp/ Vs (about Vp of 6.4-6.6 km/s, Vs of 3.4-3.8km/s, and Vp/Vs of 1.85) (Fig. 4) can be in a fluid-overpressured diabase rock.From Takei (2002), dlnVs/dlnVp, which is (1-Vs/Vs0)/(1-Vp/Vp0), where Vp0 and Vs0 are from dry diabase, is 1.71.From Fig. 5 of Takei (2002), this value of dlnVs/dlnVp corresponds to the pore aspect ratio of about 0.005 if we assume liquid in the pore as water, or 0.01 if we assume melt at a depth of 18 km.The volume fraction is estimated to be about 0.3% or 0.5% if the pore fluid is water or melt, respectively.
One of the possible other interpretations is that the anomaly (circle in Fig. 4) contains mafic granulite that was caused by a Miocene volcanic intrusion.From Christensen (1996), Vp, Vs and Vp/Vs of the mafic granulite at 600 MPa are 6.94, 3.82, and 1.82, respectively. From Takei (2002), dlnVs/dlnVp is 1.43, and this value of dlnVs/dlnVp corresponds to the pore aspect ratio of about 0.02 if we assume liquid in the pore as water, or 0.04 if we assume melt at a depth of 18 km.The volume fraction is estimated to be about 1% or 0.5% if the pore fluid is water or melt, respectively.
There are two possible causes of this high-Vp, high Vp/ Vs area (dashed-line circle in Fig. 4b-d): Tertiary magma, which caused the caldera activity in the focal region, or Quaternary magma at present.At present, no volcanic activity exists.Thus, the anomalous area suggests Tertiary magma, which formed the caldera activity.However, a current fluid path from a deep part may be possible.Umeda et al. (2009) showed a relatively large helium isotope ratio (He3/He4) value in the easternmost part of the Noto Peninsula.This large helium isotope ratio suggests that fluid is currently upwelling through the mantle.A regional-scale seismic velocity structure (e.g., Nakajima 2022) shows a continuous extent of the low Vs anomaly from a deep part to the focal area.Low Vp/Vs, where the hypocenters are located, has a Vp of about 6 km/s or smaller and Vs of about 4 km/s (marked by a broken-line square in Fig. 4b-d).These values correspond to quartzite.This low Vp/Vs area could be an area in which silicate was accumulated by fluid migration (e.g., Cox 1995;Sibson 1990).
In summary, these seismic velocity anomalies might have been formed as follows: (1) When the Sea of Japan opened, a caldera was formed in the area, and after it, the magma reservoir below the caldera slowly cooled and solidified.(2) As the magma cooled, water, in which SiO2 was dissolved, was released.The dissolved SiO2 was precipitated around the magma reservoir and formed quartz veins.The original magma became relatively quartz-poor and more mafic.(3) Currently, water was additionally supplied from deeper parts, causing earthquakes, and quartz veins further developed around them.(4) As a result, mafic rocks (high Vp/high Vp/Vs) remained in the caldera's center, and quartz concentrated in the surrounding areas, creating low Vp/ low Vp/Vs areas.

Shear wave splitting
The spatial variation of the stress field controls the geometry of the preferably opened cracks.For example, the open crack aligns horizontally under the reversefault type stress regime.Under the strike-slip fault type stress regime, the open crack aligns vertically.Thus if the crack density is the same, the delay time or anisotropy strength for the near-vertical ray path will be larger under the strike-slip fault type stress regime than under the reverse-fault type stress regime.
To check whether this variation occurs, we have calculated focal mechanisms and inverted for stress tensors (Additional file 1: Figure S11).We determined focal mechanisms from the P-wave initial motion by the HASH (Hardebeck and Shearer) code of Hardebeck and Shearer (2002).The left figure shows the focal mechanism at each epicenter, which is reverse-or strike-slip type.P-axis orientations are predominately SE-NW.
Then, we determined the stress tensor by the inversion method of Michael (1987) and Hardebeck and Michael (2006).We also did the inversion by dividing the area into four sub-areas.The solutions of the stress tensors are a reverse type with Shmax orientation of NW-SE.This orientation is similar to the previous study (e.g., Terekawa and Matsu'ura 2010).The stress ratio (= (ph2-ph3)/(ph1-ph3)) is almost 0.5.Thus, while the anisotropy in the southwest could be caused by the stress regime, the spatial variation of delay time or anisotropy strength is not caused by the stress-regime variation at depth.
Surface topography as a spatial density variation causes inhomogeneity of the stress.Some studies suggest shallow crustal anisotropy is controlled by the gravity-induced stress caused by the surface topography (Mt.Fuji, Central Japan, Araragi et al. 2015; southern Hikurangi subduction zone, New Zealand, Evanzia et al. 2017; Kaikoura earthquake, New Zealand; Graham et al. 2020).We calculated the Shmax direction and magnitude of the gravity-induced stress field with Hirschberg et al. (2019) based on Flesch et al. 's (2001) method.Using the topography of the Noto Peninsula and the thickness of the crust (mean density 2.67 g/cm3), we created a model of the depth of the Moho by synthesizing Matsubara et al. (2017) with CRUST1.0 (Laske et a. 2013).
Additional file 1: Figure S12 shows, from left to right, the topography, the gravity-induced stress field averaged from the surface to a depth of 25 km, and the Moho depth.The Noto Peninsula area has a rising topography with a long axis approximately northeastsouthwest when viewed from the seafloor.However, the topography on the peninsula is fairly flat compared to the mountains to the southeast on the mainland.The Moho discontinuity gradually becomes shallower from the central mountainous region in the southeast direction to the Japan Sea in the northwest direction.The direction of the maximum horizontal compression axis of the gravity-induced stress field near the Peninsula is in the northwest-southeast direction.This direction of the maximum horizontal compression axis may be due to the orthogonal direction of the topographic rise's long axis and the Moho surface's inclined direction.The magnitude of stress (second invariant) is less than 10 MPa.The direction of the maximum horizontal compressive axis of the gravity-induced stress field is almost consistent with the direction of the maximum horizontal compressive axis obtained from the stress tensor inversion of the earthquake.This direction is almost compatible with the direction of the maximum horizontal compressive strain axis obtained from GNSS observations (e.g., Sagiya et al. 2000).Thus, the stress field is generally considered tectonically generated, but there may be a small contribution from the gravity-induced stress field.The east-west anisotropy along the northern coast could be structural since it can match the east-west direction of the fault-fold structure along the north coast.
The regional difference in anisotropy seems to be related to the seismic velocity structure (Additional file 1: Fig. S13).The NW-SE oriented stress-induced anisotropy is distributed where high Vp/Vs at a depth of 18 km, and the structure-controlled anisotropy is distributed in the surrounding area where low Vp/Vs at a depth of 18 km.The regional difference in anisotropy also seems to be related to the degree of swarm activity.The activity began in the southern part of the study area, but the of activity was not so high.From late 2021, the activity expanded to the northern part and the swarm became more active.If the anisotropy existed in the depth of the earthquake, the preexisting structure would affect the degree of activity.It was less active where there was no distinct fault structure (small stress-induced anisotropy) in the southern part, but more active in the northern part where fault structure (large structure-controlled anisotropy) is distinct.Note that inferring from the aftershock distribution (Fig. 4) and the moment tensor solution (Fig. 1), the fault plane of the M6.5 earthquake in the northern part corresponds the southwestward dipping nodal plane with a high dip angle (e.g., ~ 56°, National Research Institute for Earth Science and Disaster Resilience 2023) and can be interpreted as a compressional inversion fault on which overpressured fluid is expected to be necessary for its reactivation (e.g., Sibson and Ghisetti 2018).This high-seismicity northern area including the M6.5 earthquake can be considered a long-term heavily fractured zone including quartz vein, with fluid, as suggested by the seismic velocity structure and other previous studies.The high pore pressure of the fluid reduces the effective normal stress and promotes the occurrence of earthquakes.This would mean the earthquakes occurred more actively on the pre-existing fault/structure in the northern part.

Conclusions
Using data from temporarily deployed seismic stations, we have performed seismic velocity tomography and shear wave splitting analysis for the Noto Peninsula earthquake swarm.Details of the obtained results are as follows.
On the seismic velocity structure, we identified a high Vp/Vs region beneath the focal area at a depth of 18 km.This high Vp/Vs region has slightly larger P-wave velocities than the surrounding area.Since Tertiary igneous rocks are distributed in the target area, the high Vp/Vs region may represent a tertiary magma reservoir, suggesting that liquids released through the magma reservoir are involved in this seismic swarm.
On the spatial distribution of S-wave polarization anisotropy, stress-induced anisotropy can cause the observed NW-SE anisotropy in the southern part.Structural east-west and north-south anisotropy may be another cause in the northern, eastern, and western regions.
No significant time variation was observed in the orientation of anisotropy.Larger anisotropy was locally observed, in particular for the activity in the northern part of the focal area, which was mainly active after October 2021.
The seismic activity started in the southern part of the area where anisotropy is consistent with the stress direction, i.e., stress anisotropy is distributed, and high Vp/Vs regions suggestive of fluid, which could trigger the seismicity, at depth are identified.The activity migrated from south to west and north (Additional file 1: Fig. S1).Subsequently, earthquakes became more active in the northern part, where structural anisotropy with well-developed fault structures is distributed.This observation suggests heterogeneous anisotropic seismogenic crust would control the gradual swarm activity.
station the temporary stations, respectively.Fig. S2.Active tectonics in northern Honshu: seismic tectonic sketch map showing epicenters of recent reverse fault ruptures and basin boundary faults, exposed areas of pre-Neogene basement, actively growing regional anticlines, Quaternary volcanoes, and volcanic front, relationship to subducting plate boundaries in the Japan Trench, and estimated maximum horizontal stress trajectories (modified from Sibson, 2009, original mainly following Sato, 1994).Red thin square denotes the area of Fig. S1.Fig. S3.Result of the checkerboard resolution test.In the checkerboard resolution test, we adopted the three-dimensional seismic velocity model of the checkerboard pattern with +-5% for P-wave and -+5% (negative) for S-wave, 1.57 or 1.91 for Vp/Vs, respectively.We calculated the travel time for each ray path of the data and added the random noise of 0.18s for P-wave and 0.20s for S-wave, respectively.Thus we did the same procedure as the real data for the synthetic travel time data to obtain the result of the checkerboard resolution test.Hypocenter location error was defined as the difference between the true (applied) location and the relocated location by the checkerboard resolution test.The result of the checkerboard resolution is shown instead of the actual Vp image in (d) and (e) in Fig. 2. Fig. S4.Result of the reconstruction resolution test.In the reconstruction resolution test, we adopted the obtained three-dimensional seismic velocity structure as the model.We calculated the travel time for each ray path of the data and added the random noise of 0.18s for P-wave and 0.20s for S-wave, respectively.Thus we did the same procedure as the real data for the synthetic travel time data to obtain the result of the reconstruction resolution test.The result of the reconstruction resolution is shown instead of the checkerboard resolution e in (d) and (e) in Fig. S3.S1.List of the temporary stations.

Fig. 1
Fig.1Geological map modified fromSawada et al. (2011).Red thick and thin lines denote the active and non-active faults, respectively.The rectangle denotes the spatial range of Fig.5b.The ellipse with broken line shows the approximate location of the inherited caldera suggested by the gravity anomaly (see Fig.3a).The moment tensor solution (NIED, 2023) of the M6.5 earthquake is shown on the right

Figure 4
Figure 4 shows the Vp and Vp/Vs structure in cross section compared to the Vp structure at 18 km depth in map view.In the map view a, there is a high Vp/Vs of over 1.8 beneath the southern part of the focal area.The crosssections b and c clearly show a south-eastward-dipping alignment of aftershocks of the M6.5 earthquake.There is a high Vp/Vs and high Vp area (marked by a brokenline ellipse) at depths > 15 km in the southern part of the focal area.The high-Vp and high Vp/Vs anomalous area beneath the earthquake swarm has Vp of about 6.4-6.6 km/s, Vs of 3.4-3.8km/s, and Vp/Vs of 1.85 at maximum.

Fig. 2 a
Fig. 2 a Hypocenter distribution.Hypocenters are colored circles in map view.The circle's color and size denote the earthquake's depth and magnitude.b Station distribution with tomography grid.Stations are plotted as the square in the map view.One broadband station is shown as bold square.The color denotes the number of the double-difference of S-wave at each station.Cross denotes the grid used for the seismic tomography analysis.c A close-up view of the station distribution with their code

Fig. 3
Fig. 3 Result of seismic velocity tomography: horizontal slice.a Bouguer gravity anomaly, modified from Sawada et al. (2012).The ellipse with broken line shows the approximate location of the inherited caldera suggested by the low gravity.b The middle and c right figures show the horizontal slice of Vp at depths of 6 km and 18 km, respectively.Earthquakes within a ± 3 km depth range are plotted.a The white and black stars in c denote the M5.4 earthquake on June 19, 2022, and the M6.5 earthquake on May 5, 2023, respectively.Grey and black dots denote the earthquake from 2019 before the M6.5 earthquake and after the M6.5 earthquake, respectively.Red small squares in b and c denote the seismic station used in this study Examples of high-quality A-grade measurements of regional events recorded at station TU.KKMS.(a) Filtered waveforms of the eastern (e), northern (n), and vertical (z) components.The vertical solid line is the S arrival time.The dashed lines are the minimum start time (1) and maximum end time (4) of the window used for processing.(b) For the window shown in gray, the original filtered waveform (top two curves) and the waveform corrected with dt determined by SC91 (bottom two curves) are shown rotated in the early S-wave polarization direction (p), and its vertical direction (p⊥) determined by SC91.The vertical red solid line is the S arrival time.The two dashed lines on either side of the vertical line indicate the range of start windows (1 and 2) and end windows (3 and 4) allowed in the SC91 measurements.(c) Values of φ and dt determined for each measurement window as a function of the window number shown on the abscissa.(d) For all clusters of five or more measurements, the large cross is the cluster finally chosen as the optimal value.(e) The (top) waveform and (bottom) particle motion of the (left) original and (right) corrected waveforms obtained by the finally selected window.(f ) Contours of the smallest eigenvalue of the covariance matrix of the finally selected measurement with SC91.Gray boxes in a, b, and e represent the time window for the final measurement.Fig. S6.Result of the time window and filter automatedly determined by MFAST.An example at the station N.SUZH is shown.(a) Histogram of the length of the time window.(b) Histogram of the starting time of the window from S-wave arrival.(c) The lower cut-off frequency of the filter.(d) The higher cut-off frequency of the filter.Fig. S7.Shear wave splitting: temporal change at station N.SUZH with decimal years on the horizontal axis.Delay time between two split shear waves divided by the distance (b) Distance between the hypocenter of the earthquake and the observed station, (c) the orientation of anisotropy, (d) delay time.Fig. S8.Median values of the delay time in 7-day moving windows at station N.SUZH.(a) Cumulative value of median of delay time normalized by the path length.The horizontal axis is the number of the window.(b) The median of delay time is normalized by the length of the path.The horizontal axis is the number of windows.(c) The median of delay time is normalized by the length of the path.The horizontal axis is the year.(d) The delay time, which is not normalized by the length of the path.The horizontal axis is the year.Fig. S9.Shear wave splitting: temporal change: map plots for five time windows.The observation at station N.SUZH is shown.The direction of anisotropy is the direction of the bar.The delay time is shown as the color of a small circle at the center of the bar.They are plotted at the epicenter location.Delay time versus hypocenter depth is also shown.The white star in (d) denotes the M5.4 earthquake on June 19, 2022.The white star in (e) denotes the M6.5 earthquake on May 5, 2023.Fig. S10.(a) Delay time and (b) fast shear wave oscillation direction (the orientation of anisotropy) versus hypocenter depth at station N.SUZH.(c) The direction of anisotropy is the direction of the bar.The delay time is shown as the color of a small circle at the center of the bar.They are plotted at the epicenter location.Fig. S11.Focal mechanisms and stress tensors.They are shown in the lower hemisphere projection.The left figure (a) shows the focal mechanism at each epicenter.We obtained 93 mechanisms from May 2021 to May 8, 2023.The right figure (b) shows the stress tensors in sub-areas with a 0.0625 x 0.0625 degrees size.Grey and white dots denote the P-and T-axis, respectively.Fig. S12.Gravitational-induced stress.This figure shows, from left to right, the topography, the gravityinduced stress field averaged from the surface to a depth of 25 km, and the depth of the Moho surface.In the gravitational stress map, red and blue arrows denote the maximum and minimum compressional stress axis, respectively.Fig. S13.The spatial average of the polarization of a fast shear wave from shear wave splitting is calculated by weighting inversely proportional to the square of the distance from the station.The red rose diagram is normalized, and the yellow bars are the average polarization.Rose diagrams are plotted at the center of each block; blocks with fewer data than the threshold are not plotted.Map view of Vp/Vs at a depth of 18km is presented as background.The ellipse with broken line shows the approximate location of the inherited caldera suggested by the gravity anomaly in Fig 3a.Table