Conjugate faulting and structural complexity on the young fault system associated with the 2000 Tottori earthquake

Young faults display unique complexity associated with their evolution, but how this relates to earthquake occurrence is unclear. Unravelling the fine-scale complexity in these systems could lead to a greater understanding of ongoing strain localization in young fault zones. Here we present high-spatial-resolution images of seismic sources and structural properties along a young fault zone that hosted the Tottori earthquake (Mw 6.8) in southwest Japan in 2000, based on data from a hyperdense network of ~1,000 seismic stations. Our precise micro-earthquake catalog reveals conjugate faulting over multiple length scales. These conjugate faults are well developed in zones of low seismic velocity. A vertically dipping seismic cluster of about 200 m length occurs within a width of about 10 m. Earthquake migrations in this cluster have a speed of about 30 m per day, which suggests that fluid diffusion plays a role. We suggest that fine structural complexities influence the pattern of seismicity in a developing fault system. Conjugate faulting and seismic velocity structure reveal the fine-scale complexity of fault growth on a young system that may be facilitated by diffusion of crustal fluids, according to a micro-earthquake catalogue from a hyperdense seismic network.

F ault systems develop with an increase in total fault displacement. During the initial stage of development, tectonic strain is accommodated by pervasive brittle deformation that can be associated with discontinuities such as bends, stepovers, and conjugate faults 1 . However, as the total displacement increases, slip gradually localizes into one major smooth fault 2 , eventually resulting to a release of tectonic strain along it. In a young fault system, fault geometries, and structural properties are likely to change in a complex manner. Complicate interactions between discontinuous faults and lateral structural properties may promote or inhibit fault nucleation and growth, and produce fault linkage and coalescence 1 . In addition, these interactions can control spatio-temporal pattern of seismic activity, slip distribution, and the location of large earthquake rupture areas, which are important in seismic hazard assessment 3,4 . Recent seismic and geodetic observations of fault networks have revealed that multiple ruptures along conjugate and subparallel faults can occur simultaneously or in quick succession [4][5][6][7][8] . Crustal fluids spreading through a fault network can also enhance fault-slip [9][10][11][12] . However, fault complexity along young fault systems is poorly understood owing to a lack of highspatial resolution images of seismic sources and structural properties at seismogenic depths.
To address this issue, we deployed a spatially dense seismic array of 1000 vertical short-period sensors in and around the source area of the 2000 Tottori earthquake (Japan Meteorological Agency (JMA) magnitude Mj 7.3; Mw 6.8), and obtained continuous waveform data from April 2017 to April 2018 (13 months) (Fig. 1). The Tottori earthquake took place along a left-lateral strike-slip fault within the San-in active shear zone, which is characterized by a high seismicity rate and a large shear strain rate of 1.0 × 10 −7 /year [13][14][15] . Along the San-in shear zone, a conjugate seismic lineation with WSW-ENE and NNW-SSE trends can be recognized from the JMA catalog (Fig. 1a). A Global Navigation Satellite System network has revealed that the shear deformation zone is as wide as 50 km along the seismic belt 15 . However, no major active faults have been identified at the surface in the San-in region (Fig. 1a). A detailed geomorphological study has clarified that conjugate active faults with WSW-ENE and NNW-SSE trends are distributed within the San-in shear zone, but that they are young growing faults with cumulative offsets of less than a few hundreds of meters 16 . Indeed, no obvious surface traces of active faults linked to the 2000 Tottori earthquake have been identified (Fig. 2a), indicating that a major fault structure has not developed in the study area. These observations suggest that the San-in shear zone is a young fault system with many faults that lie parallel and oblique to the shear zone 15 , but there is little constraint on the small-scale structures within the shear zone.
The average spatial interval of the hyperdense seismic network used in this study is~1 km, covering the entire aftershock area within a radius of~15 km (Fig. 1b) to be installed for more than 1 year, and above the aftershock area of a magnitude~7 earthquake 17,18 , to the best of our knowledge. Using the waveform data retrieved from the hyperdense network, we obtained arrival time data for 4033 micro-earthquakes with magnitudes greater than the completeness magnitude of M L −0.6 ( Fig. 1c). Dense, well-covered ray-paths from the many sources provide valuable opportunities to (1) precisely map the aftershock distribution following the 2000 Tottori earthquake (Fig. 1d), (2) delineate fine-scale P-wave velocity structures with a spatial resolution down to 500 m, (3) explore the spatial and temporal evolution of a small seismic cluster, and (4) understand the multiscale complexity of this young fault system at depth.

Results
Geometric complexity of faults and velocity structures. Nearly 20 years have passed since the 2000 Tottori earthquake, but the relocated hypocenters clearly reflect complex fault structures that extend along a dominantly NW-SE trend in the source area (Fig. 2a). The relocated seismicity shows that the seismogenic zone extends down to 15 km depth, gradually decreasing tõ 5 km depth toward the tips of the fault (Fig. 2b, c). This pattern of detected seismicity is similar to that of aftershocks that immediately followed the mainshock rupture 14,19 . This means that micro-seismicity (down to negative magnitudes) can be used to reconstruct geometric complexities produced by the mainshock rupture, even after~20 years. The JMA catalog shows a similar spatial pattern of seismicity ( Supplementary Fig. S1), but the geometric complexity (the cross-section numbers from 15 to 23 in Fig. 2a) is poorly defined due to the sparse coverage of the nationwide online seismic network.
Near the initiation point of the mainshock rupture, most events are well-aligned along a steeply dipping fault that strikes N35W-S35E (Fig. 2b). Near the mainshock epicenter, the major ruptured area is characterized as a high-velocity body ( Fig. 3 and Supplementary Figs. S2-S6). Within the major ruptured area, aftershocks are notably less common (Figs. 2c and 3). This is a common feature in the largeslip regions of many other earthquakes 20,21 .
In contrast, many relocated events on the northwest side of the source area have conjugate vertical alignments (blue arrows in Fig. 2a). The data show that the fault zone exhibits conjugate faulting over multiple length scales. The most prominent conjugate fault plane extends toward W20S from the main trend, with a horizontal length of~5 km along strike. This prominent aftershock lineation occurs close to the distinct boundary between a low-velocity body and a high-velocity body, on the northwestern side of the mainshock epicenter ( Fig. 3 and Supplementary Figs. S2-S6). The boundary between the lowand high-velocity bodies trends W20S-E20N, and the observed seismicity rate is higher along this boundary than in other areas. The magnitude of the velocity contrast between the NNW and SSE sides of this boundary increases from deep to shallow depths (Fig. 3). As discussed in previous studies 19,22 , the shallow remarkable low-velocity body might be composed primarily of non-alkali volcanic and pyroclastic rocks of the early to middle Miocene. However, the present model shows that the velocity contrast extends down to greater depths. The high-velocity body might be interpreted to be plutonic and metamorphic rocks from surface geological map.
Toward the southeastern edge of the fault zone, the relocated aftershocks begin to change dip slightly with depth. The aftershock orientations become listric, transitioning to a dip of 75°toward W35S at a depth of around 8 km (slices 04-10 in Fig. 2b). This dipping direction is opposite to the aftershock orientations dipping~75°toward E35N observed on the northwest side of the source area (slices 17-21 in Fig. 2b). We also  recognized a weak bifurcation geometry in the aftershocks at depths between 4 and 6 km (broken lines along the slice 02 in Fig. 2b). In terms of velocity structure, the southeast end of the aftershock area is located in a relatively low-velocity body for each of the studied depth intervals, except for the depths of 3, 4, and 10 km (Fig. 3 and Supplementary Figs. S2-S6).

GVp(%)
In the central fault section there is a thin, low-velocity band that lies subparallel to the alignment of the majority of the aftershocks. This low-velocity band is offset by~3 km to the west, as indicated by the white arrows in Fig. 3. A similar low-velocity band can still be recognized at a spatial grid spacing of 1 km ( Supplementary Fig. S3), albeit less clearly than the band in Fig. 3. The horizontal and vertical dimensions of this low-velocity band are~5 km and~3 km, respectively. No aftershocks have been detected along the low-velocity band, indicating that it may be located within a stress shadow zone produced by the 2000 Tottori earthquake 23 . The spatial resolution of the tomographic image is limited to 500 m, so it is difficult to accurately image the width of the buried low-velocity band, but the width of the low-velocity band appears to be comparable to or less than 500 m.
Spatial and temporal evolution of a tiny seismic cluster. The hyperdense seismic network of this study makes it possible to relocate micro-earthquakes accurately and to depict numerous small seismic clusters (Figs. 2 and 3). Here, we focus on a tiny seismic cluster (boosted by a M L 2.2 event) near the edge of the prominent conjugate fault on the NNW side of the fault zone (a small red rectangle around the slice 19 in Fig. 2a). It has a leftlateral strike-slip focal mechanism and displays a clear lineation (Fig. 4a). Most of the relocated events at~5.65 km depth are aligned along a vertically dipping plane that is subparallel to the mainshock fault plane (Fig. 4b). The tiny delineated fault is 200 m long in the direction of strike and 150 m long in the direction of dip (Fig. 4c). The lineation localizes to a width of around 10 m (Inset in Fig. 4a), which is larger than the horizontal uncertainty (6.5 m: 95% confidence limit) of the relative hypocenter relocation (see "Methods").
The tiny seismic lineation consists of several small sub-clusters that are present in the central to northern parts of the structure (gray arrows in Fig. 4a). The fault-length of each sub-cluster is 30 m and the sub-clusters are separated by~15 m in the horizontal direction. High-accuracy focal mechanism solutions show that one of the nodal planes (NNW-SSE trend) is parallel to the direction of the major alignment ( Supplementary Fig. S7). The earthquake alignment and focal mechanisms clearly indicate that this seismic sequence was macroscopically produced by leftlateral strike-slip movement.
The spatio-temporal evolution of this tiny fault zone can be used to define two sequences (Fig. 4d-f  (mainshock), and another characterized by swarm-like seismicity that initiated several days after the first mainshock-aftershock sequence. Based on the distribution of the immediate aftershocks, the M L 2.2 event is inferred to have initiated near the center (Fig. 4d) and down-dip (Fig. 4e) of the rupture area. Interestingly, the second swarm-like seismicity sequence migrated along the strike of the fault toward its southeast edge, as well as in the down-dip direction. The migration distance was around 100 m over 4 days. Thus, the migration speed of the second sequence was~30 m/d.

Discussion
The ultra-high density of the seismic network shows that conjugate vertical lineations of relocated seismic events could be defined using very small-magnitude earthquakes (down to M L −0.6). Sub-vertical alignments can be identified in a crosssection along the fault-strike of the 2000 Tottori earthquake (Fig. 2c). These are interpreted to be conjugate faults related to the mainshock fault plane. Conjugate vertical faults were also observed during the 2019 Ridgecrest earthquake sequence in California 8 and the 2012 Indian Ocean earthquake 24 , in strikeslip stress domains. A mixed rupture of conjugate strike-slip and thrust faults was depicted during the 2018 northern Osaka earthquake, Japan 25 .
Most studies on the source process of the 2000 Tottori earthquake 26,27 have suggested that the largest coseismic-slip occurred around the central section (near the rupture initiation point), with a small amount of coseismic slip at both ends of the aftershock area (Fig. 2c). The central fault section at depths shallower than 6 km corresponds to the high-velocity body, with no obvious low-velocity fault zone visible along the aftershock lineation at 500 m spatial resolution (−3 < Y < 6 km in Fig. 3), showing relatively weak aftershock activity (Fig. 2c). Thus, the 2000 Tottori mainshock must have ruptured an immature fault zone (close to intact rock), causing the largest coseismic slip to occur in the high-velocity body at depths shallower than 6 km. Previous tomography studies have also investigated structural heterogeneities on fault planes, revealing that high-velocity features with low aftershock activity correspond to areas of large coseismic slip 20,28,29 . Laboratory measurements of an exhumed fault zone 30 reveal that microfractures in the fluid-alteration zone are sealed leading to relatively high P-wave velocities compared with the surrounding wall rock. Sealing of microfractures in fault zone after the dynamic rupture might contribute to postseismic recovery of seismic wave velocity, resulting to vanishing of a low-velocity damaged zone associated with the 2000 Tottori earthquake.
Even after~20 years since the mainshock, the aftershock activity in the large-slip area has been persistently low (Fig. 2c), suggesting the potential use of micro-seismic activity as an indicator of large-slip area of historical earthquakes.
High-resolution seismic tomography has revealed a thin, lowvelocity band (the length is~5 km and the down-dip height is 3 km) that is subparallel to the mainshock fault and is offset bỹ 3 km (white arrows in Fig. 3). Low velocities can be caused by fractures filled with crustal fluids. Thus, the thin, low-velocity band is interpreted to be a blind fault zone within the brittle crust, where intense and pervasive fracturing has accumulated from previous events. This localized low-velocity feature (<0.5 km wide) matches with fault zone property imaged by fault-trapped waves along several active faults, including those associated with the 1992 Landers earthquake and the 1999 Hector Mine earthquake within the eastern California shear zone in USA 31,32 . In contrast, the Calico fault located midway between the Landers and Hector Mine ruptures has a ∼1.5-km-wide low-velocity zone 33,34 . Differences in the inferred widths of low-velocity zones may be due to intrinsic variations in the width of damaged zones around different faults and/or different spatial resolutions of geophysical explorations 33 .
The distinct low-velocity anomaly is present at the northwest side of the aftershock area. The aftershock distributions are also more complex than in the central section ( Fig. 3 and Supplementary Fig. S2). Multiple conjugate faults are developed on the northwest side of the Tottori mainshock rupture area, where intense aftershocks followed the mainshock event 35 . The lowvelocity zone could represent developing fracture networks at the fault edge that are filled with crustal fluids. The mainshock rupture initiated in the high-velocity body at the central section, propagated to the northwest side 27 , and caused shear deformation in the soft volcanic rocks (low-velocity zone), resulting to the development of the conjugate faults.
We revealed the spatio-temporal migration of earthquakes in the tiny cluster both along fault-strike and down dip (Fig. 4). The migration occurred at a speed of~30 m/d, down to the deep part of the M L 2.2 source area. During the migration, each seismic subcluster was activated after the propagation front had passed. Earthquake migration is believed to be due to fluid diffusion, aseismic slow-slip transients, or a combination of the two [36][37][38][39] . The earthquake migration speed within the tiny seismic cluster is much slower than the typical migration speed of slow-slip events along subduction zones 40 . Therefore, the earthquake migration could reflect infiltration of crustal fluid into the down-dip extension of the coseismic rupture area of the M L 2.2 event.
Brittle deformation associated with the M L 2.2 event could have induced a local stress concentration in the deep parts of the fault, facilitating the development of fractures. Crustal fluid might have infiltrated these fractures and triggered the migrating seismicity 41 . In addition, cascading ruptures for which small earthquakes trigger each other by stress transfer could play a role to induce the earthquake migration 42 . However, cascading ruptures may be less important as far as we assume a plausible range of static stress drop, as each migrating earthquake is not located at the edge of the rupture zone of the previous event ( Supplementary Fig. S8).
Diffusional earthquake migration has been observed in fluidinduced earthquake experiments 36 , during tectonic swarm activity 12 , and in active volcanic areas 43 . Upwelling of deep fluids into the San-in shear zone has also been implied from seismic tomography and geo-chemical analysis 44 .
From the center of the M L 2.2 fault to the north, three sets of sub-clusters can be recognized. Each cluster is~30 m long and rotates slightly counterclockwise from the main trend (gray arrows in Fig. 4a). The sub-clusters are almost parallel to each other and have a spacing of~15 m. These stepped, subparallel faults are similar to en echelon cracks, which can be observed at the surface following a strike-slip earthquake rupture 45 . The orientation of the sub-cluster relative to the main fault suggests that it could be Riedel shears formed in association with leftlateral strike-slip faulting. Riedel shear structures often form at low levels of cumulative displacement. With a progressive increase in cumulative displacement, Riedel structures will organize into a dense elongate network, eventually coalescing to form a smooth mature fault zone 46,47 .
The conjugate faults, potential blind fault, heterogeneous velocity structure, and Riedel shears indicate that the 2000 Tottori earthquake occurred as a part of a fault growth process in a young fault system. It appears that the fault system will develop into a mature, large-scale tectonic boundary over long time. The spatiotemporal pattern of seismic activity and the fault-slip distribution appear to be controlled by complexities within the young fault system. The smallest scale of crustal deformation that is resolved by our current seismic network (~10 m) is approaching that of a geological survey. Future hyperdense seismic networks could unravel deformation across multiple length scales, with meter or sub-meter spatial resolution. This could lead to a greater understanding of ongoing strain localization in young fault systems possessing unique complexity.

Methods
Hyperdense seismic observation. We deployed 1000 temporary seismic stations in and around the source region of the 2000 Tottori earthquake from April 2017 to April 2018. The average spatial interval of the hyperdense seismic network was~1 km, covering the entire aftershock area (Fig. 1). Each seismic station was equipped with a short-period vertical-component sensor with a natural frequency of 2.0 or 4.5 Hz. The waveform data were recorded continuously at a sampling rate of 100 Hz and was stored on an SD card in each data logger. The internal clock accuracy was maintained to within 1 ms by calibrating at 6-h intervals, using the GPS receiver on each data logger. Approximately 800 stations transferred continuous data every 6 h via mobile internet, using a SIM card installed in each data logger. The data loggers were operated by dry batteries for over 1 year. This low-powerconsumption data logger has been newly developed in cooperation with Kinkei System Corporation for the purpose of hyperdense seismic observation.
Seismic tomography analysis. To detect micro-earthquakes, we applied an automatic processing workflow (composed of phase picking and phase association) to the continuous waveform data retrieved from the hyperdense seismic network. All of the sensors were vertical components, so we focused exclusively on P-wave arrival time data. The initial hypocenter was calculated using P-wave travel times, assuming the one-dimensional velocity structure routinely used by JMA.
To relocate the detected earthquakes and estimate velocity structure, we applied Double-difference travel-time tomography code 48 to the arrival time data. The earthquakes with a sufficient number of arrival times (i.e., the number of P-arrivals must be greater than 23; the median value of P-wave arrivals for each event is 209) and a low average residual of P-wave travel times were then selected for the tomographic step. The total number of selected events was 4033 (Fig. 1b), which is about nine times the number of events recorded in the JMA catalog for the same period (455 events) ( Supplementary Fig. S1). This is a significant improvement in the completeness of magnitude down to M L −0.6 (Fig. 1c). The number of absolute arrival times of P-waves used for tomography analysis was 1.13 million, and the number of differential travel-time data for the catalog-picked P-waves exceeded 3.64 million. The JMA velocity model was used as the starting model for tomographic analysis (Supplementary Table S1). To obtain more precise differential P travel-time data, we applied the waveform cross-correlation technique for earthquake pairs with a spatial separation of less than 2 km. We applied a two-way 4-10 Hz Butterworth filter to waveform data and used a 1.2 s window length beginning 0.6 s before the automatically picked arrival time. We then computed the differential travel-time for each earthquake pair to subsample precision using time domain cross-correlation. After selecting the differential P travel-time with cross-correlation exceeding 0.85, we obtained a dataset of more accurate differential travel-time data that contained 8.47 million P-wave observations.
The grid nodes for the tomography were located at −300.  (Fig. 1d). Near the mainshock failure area, the grid spacing was 0.5 km (Fig. 3j), which is much finer than that used in previous studies 19 .
We also used a velocity model estimated by ref. 19 as the starting model for the tomographic analysis. However, the final velocity model does not show any significant change from the model that was inverted using the JMA velocity model as the initial model (Fig. 3). That is because the present study has used more than 20 times seismic stations and huge number of arrival time data, compared with ref. 19 .
Spatial resolution of the P-wave velocity model. We tested different spatial smoothing weights of 1, 2, 3, 5, 7, 10, 15, 20, 30 45, and 60 (see ref. 48 for details) to verify the robustness of the velocity model. Despite changing the spatial weights, the main features of the tomograms remained similar in areas with good ray coverage ( Supplementary Figs. S5 and S6). Considering the trade-off between roughness and stabilization of the model (Supplementary Fig. S9), we chose the velocity model with a weighting of 10 for our preferred final model (Fig. 3). The root mean square of the absolute travel-time residual was reduced from 0.057 s to 0.022 s after 28 iterations.
To evaluate the spatial resolution of the final velocity model, we conducted a checkerboard resolution test. We computed synthetic travel times for a threedimensional velocity model using perturbations of ±3% amplitude of the initial velocity model at each inversion grid, adding random noises uniformly distributed from −0.025 s to 0.025 s for P-waves. We then inverted the synthetic data using the same distribution of hypocenters and stations as the real data, starting from the initial velocity structure, and verified how well the test model is recovered (Supplementary Figs. S2 and S4). We demonstrated that P-wave velocity anomalies are well recovered in the areas where there are sufficient seismic events ( Supplementary Figs. S2 and S4) and the derivative weight sum (DWS) 49 values are sufficiently large. Note that the entire aftershock area is well resolved, with horizontal and vertical spatial grid intervals of 0.5 km and 1.0 km, respectively, owing to the high density of seismic stations.
In order to assess the reliability of the model resolutions for the velocity anomalies described in the present study (Fig. 3), we also conducted a sensitivity (spike) test using a validation model shown in Supplementary Fig. S10a. Lowvelocity anomalies were assigned to grid nodes at the northwestern side of the mainshock epicenter (−4% at depths from 2 to 6 km, and −2% at depths from 7 to 10 km). Furthermore, a thin low-velocity band of −3% are given to grid nodes at depths from 5 to 8 km, and a high-velocity zone with +4% are assigned to grid nodes along the fault-strike nearby the mainshock epicenter (at depths from 2 to 10 km) (Supplementary Fig. S10a). The calculation of the synthetic travel-time and the inversion method were the same as those employed for the checkerboard model resolution test. The results shown in Supplementary Fig. S10b indicate that the assumed velocity anomalies of P-wave were well recovered, although the recovery of the greater depths was weaker than that of the shallow depths. Based on these results, we consider that the presence of the velocity anomalies discussed in the present study is reliable.
Relocation of a seismic cluster associated with M L 2.2. We relocated a tiny seismic cluster, boosted by an event with M L 2.2. We used differential arrival times from the nearby seismic stations, with an epicentral distance of 8 km from the cluster centroid (189 stations among 1000). This seismic cluster contains 66 events that occurred from 24 December 2017 to 30 January 2018 (Fig. 4). We applied a double-difference relocation algorithm 50 to differential travel times constructed by automatic picking (46 thousand P-wave data) and waveform cross-correlation technique (84 thousand P-waves with cross-correlation coefficient ≥ 0.85). We used the one-dimensional velocity structure routinely employed in the JMA location procedures. To estimate the relative location uncertainties (95% confidence limit) in the horizontal and vertical directions, we conducted a jackknife test. After generating 300 sets of differential travel times with 80% of randomly selected catalog and cross-correlation data, we applied the same relocation algorithm to all of the resampled datasets. After removing events with relatively large location uncertainties, we finally obtained the highly accurate hypocenters of 56 events. The relative location errors (95% confidence limit) are estimated to be 3 m and 12 m in longitudinal and latitudinal directions, respectively, and 19 m vertically.
Focal mechanisms. We estimated focal mechanisms of several earthquakes with M L greater than 1.3 included in the tiny relocated seismic cluster, boosted by the M L 2.2 event (Fig. 4a), applying the HASH code 51 to the first-motion polarities picked by automatic processing. The takeoff angles were computed using the onedimensional velocity structure routinely used in the JMA location procedure. The focal mechanisms have the highest grade (quality A) and are well constrained based on large number of polarity data (>640). Ref. 52 has discussed a non-double couple component of several focal mechanisms using waveform data retrieved by the hyperdense seismic network consisting of~1000 seismic stations.

Data availability
All seismic data needed to evaluate the conclusions in the paper are available in the Zenodo repository, https://zenodo.org/deposit/4059423. Additional data related to this paper are available in the paper and/or the Supplementary Materials and can be requested from the authors.

Code availability
The code of the "HASH" used to determine focal mechanism can be found at https:// www.usgs.gov/software/hash-12. The code of the "hypoDD" used to relocate hypocenters can be found at https://www.usgs.gov/software/hypodd. The code of the "tomoDD" used to conduct seismic tomography with simultaneous relocation of hypocenters is available from H.Z. (email: zhang11@ustc.edu.cn) on request. We plotted figures using GMT5 (http://gmt.soest.hawaii.edu/projects/gmt). The other custom codes used in this study are available from A.K. on request.