Detection of urban hidden faults using group-velocity ambient noise tomography beneath Zhenjiang area, China

The Mufushan-Jiaoshan fault (MJF) is a hidden active fault located on the north side of the Ningzhen Mountain Range and developed along the Yangtze River in Zhenjiang area, China. In this paper, the structure of MJF is detected and studied using group-velocity ambient noise tomography. In the study area (18 km × 25 km), 47 short-period seismic stations were deployed with the average station spacing of about 3 km and 24 days (from 27 February to 22 March 2019) of continuous ambient-noise recordings were collected. And 510 group velocity dispersion curves in the period band 0.5–5 s were extracted using the vertical component data. And then the three-dimensional shear-wave velocity structure was inverted using group dispersion data by the direct surface-wave tomographic method. Our results are consistent with the geological background of the study area, showing that in the depth range of 0.6–1.5 km, the north side of MJF presents a relatively high velocity, and the south side presents a distribution pattern of high and low velocity. While in the depth range of 1.5–2.0 km, the shear-wave velocity (Vs) model is relatively simple with relatively low velocity on the north side and relatively high velocity on the south side. And the gradient zone of Vs may be the location of the main fracture surface of MJF. The good correspondence between the Vs model and the fault structure indicates that the ambient noise tomography method can be used as an effective method for detecting hidden faults in urban environments.

. Geologic structure, faults, river systems, earthquake location, and station location in the study area and its adjacent regions. The inset figure shows the locations of the study area and its adjacent regions (red box) and the NCB, LYP represent North China Block, and Lower Yangtze Plate, respectively. (The software used to create the map is MAPGIS 6.7, http://www.mapgi s.com.cn/). In order to get more reasonable results, the standard ambient noise processing procedure 25 were followed to process our data. Firstly, we re-sampled the vertical component data, which was cut into hourly segments, at a sampling rate of 10 Hz to reduce the cost of computational time. Secondly, we removed the mean and trend of the data, and also band-pass filtered the data in the frequency band from 0.5 to 5 s. Thirdly, after these steps, to suppress the effects of earthquakes, the spectral whitening and temporal normalization were performed 26 . Finally, we conducted cross correlation to obtain the time domain cross-correlation functions (CFs) between two stations for each hourly data with the lag time from − 40 to 40 s. Using normalized linear stacking method, we stacked all the hourly CFs from each station pair and obtained the time-domain empirical Green's functions from the time derivative of the CFs 8,9 . Figure 3a show the CFs with signal-to-noise ratio (SNR) greater than five in the period band 0.5-5 s after stacking and the surface-wave signals of about 2.5 km/s can be clearly observed. Here, the SNR is defined as the ratio of the maximum amplitude of the signal window of CFs (− 20 to ~ 20 s) and the average absolute amplitude of the noise window (− 40 to ~ 20 s and 20 to ~ 40 s). In addition, the additional arrivals with higher velocities (5 km/s) might be higher-mode surface waves and the phases at − 30 s and 30 s might be related to some kinds of interference signals, which cause little impacts on the cross-correlation calculation results 27 . In the period band of 0.5-5 s used in this study, the ambient noise is usually regarded to be related to the manmade sources (e.g., traffic, cultural, and human activities), the scattering from teleseismic earthquakes, and the "leaking" from some second microseismic peak 28,29 . To check the distribution of the ambient noise sources in the study area, the normalized amplitude of the positive-and negative-time CFs for all station pairs were computed (Fig. 3b). The ambient noise sources from the east are slightly stronger than those from other directions (Fig. 3b), which may be related to the East China Sea. But from the overall perspective, most values of normalized amplitudes are between 0.5 to 1.0, which means that there are small amplitude differences between positive-and negative-time CFs and that the noise source distribution will not cause large bias in dispersion measurements 30 .

Scientific Reports
After stacking the positive-time and negative-time parts of the time-domain empirical Green's functions linearly, we used the quick tracing method EGFAnalysisTimeFreq based on frequency-time analysis and image analysis technique to first construct the velocity-period spectrogram for each TDEGF and then extract the fundamental mode Rayleigh wave group velocity dispersion curves of each station pair by picking the peak of the envelope function of the narrow-band filtered signal 31,32 . In order to satisfy the far field approximation of surface wave propagation, the interstation distance is required at least 2 times the wavelength 31 . All the 510 group velocity dispersion curves in the period band 0.5-5 s for the fundamental mode Rayleigh waves are shown in Fig. 4. It can be seen that, the quantity of dispersion data decreases with increasing period and there is no measurement can be made above the 5 s period due to relatively short interstation distance. In addition, as inferred from the   T01   T02  T03   T04   T05   T06   T07   T08   T09   T10   T11   T12   T13   T14   T15   T16   T17   T18   T19   T20   T21   T22   T23   T24   T25   T26  T27   T28   T29   T30   T31   T32   T33   T34   T35   T36   T37   T38   T39   T40   T41   T42   T43   T44   T45   T46  www.nature.com/scientificreports/ dispersion curve, the velocity variation in the study region appear to be very large; For example, at the period of 0.5 s, group velocity varies from 1.3 to 3.7 km/s.

Method
In this study, we used the direct surface-wave tomography method, which is based on frequency-dependent ray tracing and a wavelet-based sparsity-constrained tomography inversion, to invert group velocity dispersion data for the 3D shear velocity structure 33 . In comparison with the traditional two-step surface wave inversion method, this method avoids the intermediate step of inversion for group velocity maps and considers ray-bending effects of surface waves due to complex velocity structure based on the fast-marching method 34 . The goal of this method is to find a model m that minimizes the differences δt i (ω) between the observed times t obs i (ω) and the model prediction t i (ω) for all frequencies ω and the travel-time difference at angular frequency ω with respect to a reference model for path i is given by Fang, et al. 33 www.nature.com/scientificreports/ where t i (ω) is the calculated travel time from a reference model which can be updated in the inversion, and v ik are the bilinear interpolation coefficients along the ray path associated with the ith travel-time data, the group velocity C k (ω) and its perturbation δC k (ω) at the kth two-dimensional (2D) surface grid point at angular frequency ω . The group velocity perturbation can be written as Based on the empirical relationships proposed by Brocher 35 , compressional wave velocity and density are related to shear wave velocity, thus Eq. (1) can then be expressed as where θ k represents the 1D reference model at the kth surface grid point on the surface and ∝ k z j , β k z j , and ρ k z j are, respectively, the compressional velocity, shear velocity, and mass density at the jth node in the depth direction, and M = KJ represents the total number of grid points of the 3D model 33 . And R ′ ∝ z j and R ′ ρ z j are the scaling factors derived from the empirical relationships 35 .
The direct surface-wave tomography Eq. (3) can be further expressed in matrix form as where d represents the surface-wave travel-time residual vector for all ray paths and periods, and G and M represent the data sensitivity matrix and the model parameter vector, respectively.

Inversion details and model resolution tests
We here used the group dispersion data and an average 1D model, which was determined using the average group velocity dispersion curves, as the initial model to invert the near surface shear velocity structures. For the inversion grid, we parameterized the study area into 18 by 18 grid points on the horizontal plane and 9 grid points along the depth direction. The grid intervals along the latitude and longitude are 0.014° and 0.025°, respectively. At the depth direction, 9 grid points were set from 600 m to 2 km. To evaluate the model resolution, the path coverage test and checkerboard resolution test were performed. Figure 5 presents the path coverage of group velocity measurements at four selected periods (0.5, 1.0, 2.0, 3.0) based on the final 3D velocity model. It can be clearly seen that the path number decreases with the increase of www.nature.com/scientificreports/ the period (from 484 at 0.5 s to 59 at 3.0 s). And generally speaking, the ray coverage at these periods is relatively good except for some marginal area, indicating that in the most of the study area, this dataset has the capability to resolve the structure. To further test the model resolution at different depths, the checkerboard resolution test was performed. The synthetic checkerboard model was set with − 40% and + 40% anomalies and each anomaly has a size of 3 km and 5 km in the latitude and at longitude direction, respectively. Figure 6 shows the initial checkerboard model and its recovery after inversion at five different depths. We can notice that, the resolution of checkerboard models decreases as the inversion depth increases and for each inversion depth, checkerboard models are recovered better in the central part of the study area compared with the marginal area, which is closely related to the path coverage. The results of path coverage test and checkerboard test suggest that, in the depth range of 0.6-1.6 km, the structures of most regions can be resolved well, while that of some marginal area is not resolved so well (less well fitted) due to insufficient ray paths. And then we inverted the real group-velocity dispersion data using the direct surface-wave tomography method mentioned above. After eight iterations, the root-mean-square RMS residual value decreases from about 1.1 s to about 0.75 s and remains stable, indicating that the inversion converges to a stable condition. And compared with that before iterations, the residual distribution after iterations is obviously closer to zero and follows the Gaussian distribution, which also means our inversion process has good convergence and the data fitting is relatively good (Fig. 7). Figure 8 shows the V s model at different depths in the study area. It can be seen that (1) The crust medium in the study area presents strong lateral velocity heterogeneity. Taking MJF as the boundary, the V s on the north side is distributed in the near north-south direction with relatively low velocity, which may be related to the Yizheng depression, while on the south side, it is distributed in the near east-west direction with relatively high velocity, which could be related to the NMR; (2) Starting from the depth of 1.0 km, two high-velocity anomalies with a velocity of about 3.3 km/s began to appear in the northwest and southeast of the study area, and extended to the center forming a whole at a depth of 1.6 km. Since Mesozoic, the Lower Yangtze fault block has developed strong intermediate-acid magmatic activity and the MJF has invaded a lot of magmatic rocks near Shiye Town 22 . We think that the high-velocity anomaly in the northwest may be related to the magmatic rocks intruded in this area. While the high-velocity anomaly in the southeast probably be relevant to Silurian sandstone (S 1 g, S 2f. ) and Devonian shale (D 3 w) which widely developed in this area. And this lithologic distribution possibly be related to the intersection of DJF and RSF in this area (Fig. 1); (3) The MJF is developed in the V s anomaly gradient zone. Specifically, the west section of MJF is located on the transition zone between two relatively low-velocity areas on the north and south sides; the middle section is located on the transition zone between the relatively www.nature.com/scientificreports/ high-velocity area on the north side and the relatively low-velocity area on the south side; and the east section is located on the edge of the low-velocity anomaly on the north side.

Results and discussion
In order to study the relationship between MJF and V s in depth direction, we drew three cross-sections of V s model along AA' (15 km), BB' (15 km), and CC' (13.5 km) along NW-SE direction (Fig. 9), and compared the AA' section with bedrock geological section DD' (Fig. 10). The profile DD' is located about 2 km east of AA' section, with a total length of 18.9 km, and the direction is basically the same as that of section AA' . See Fig. 1 for the location of profiles AA' , BB' , CC' , and DD' .
We can obviously see from Fig. 10 that in the depth range of 0.6-1.5 km, the V s model along profile AA' has a good correspondence with bedrock geological profile DD' . According to the distribution of velocity anomaly www.nature.com/scientificreports/ along profile AA' , we divide it into five parts with dotted lines, which correspond to the Pukou Formation sandstone of Mesozoic (K 2 p), quartz diorite of late Yanshan (δo 5 3 (1) ), sandstone of middle Silurian (S 2 ), sandstone of lower Silurian (S 1 ), and sandstone of middle Silurian (S 2 ) in bedrock geological profile DD' from north to south, respectively. In addition, it is worth noting that the tendency of velocity interface in profile AA' is also consistent with that of lithology interface in profile DD' (such as the interface between K 2 p and δo 5 3(1) , δo 5 3 (1) and S 2 , S 1 and S 2 ). The good correspondence between profiles AA' and DD' shows that the V s model can reflect the change of regional geology to a certain extent. And the V s models along profiles BB' and CC' are similar to those of profile AA' , and also have a good corresponding relationship with bedrock lithology. Taking the MJF as the boundary: on the north side of profile BB' is the Pukou Formation sandstone of Mesozoic (K 2 p), showing a relatively high velocity, on the south side is the quartz diorite of late Yanshan (δo 5 3 (1) ) and the Huangmaqing formation sandstone of Middle Triassic (T 2 h), both exhibiting relatively low velocity, and the latter is slower than the former; on the north side of CC' section is the Pukou Formation sandstone of Mesozoic (K 2 p), showing a relatively high velocity, and on the south side is the quartz diorite of late Yanshan (δo 5 3 (1) ) showing a relatively low velocity and porphyritic quartz diorite (ηοπ 5 3(1) ) showing a relatively high velocity. While in the depth range of 1.5-2.0 km, the velocity distribution characteristics of the three sections are similar, showing relatively low velocity in the north and relatively high velocity in the south. According to the characteristics of V s distribution, it can be determined that the lower regions under station T09 (section AA'), T20 (section BB'), and the area between T43 and T44 (section CC') are the development zones of MJF.
In general, the V s model in study area is obviously non-uniform in both horizontal and vertical directions. In the depth range of 0.6-1.5 km, the medium on the north side of MJF fault presents relatively high velocity, while on the south side, it presents a pattern of alternating low and high velocity distribution, which is consistent with the multi-stage geological tectonic activities during the formation of MJF fault and NMR; in the depth range of 1.5-2.0 km, the velocity structure is relatively simple, with a relatively low velocity in the north and a relatively high velocity in the south, indicating that the lithology tends to be unified with the increase of depth, which is consistent with the structural background of Yizheng depression in the north of the fault and NMR in the south.

Conclusion
In this study, we used group-velocity ambient noise tomography to investigate the structure of MJF in Zhenjiang area, China. 510 group velocity dispersion curves in the period band 0.5-5 s were used to invert 3D V s structure by the direct surface-wave tomographic method. Our results reveal that, in the study area, the distribution of the V s in different depth regions has different characteristics, which is consistent with the geological background of the area. And there are obvious differences in V s distribution between the north and south sides of MJF. According to the three cross-sections of V s model, the location of the development zone of MJF is also determined. The good agreement between the V s model and the distribution of MJF as well as the geological background in the study indicate that, the ambient noise tomography method can be used as an effective method for detecting hidden faults in urban environments.