Detecting fracture zones in granite with seismic and audio-magnetotelluric methods

Geological disposal is a feasible and safe method for dealing with the high-level radioactive waste problem at present. The Beishan area is the key area preselected for high-level radioactive waste geological disposal in Gansu Province, China. The Jijicao rock block is currently the most extensively studied area, where there are several fracture zones in borehole BS15. In 2011, a remarkable conductive anomaly near the BS15 was detected using a 4-km-long audio-magnetotelluric (AMT) profile with 50 m station space crossing the borehole. To study the anomaly, a short seismic survey was shot along the part of the AMT profile. This paper presents an example for detection of the fracture zones in granite using AMT and seismic methods at Beishan. The AMT data inversion model agrees with the borehole well-logging data and weak reflections are related to the fractured zones. The seismic detection may help in interpreting the cause of the conductive anomaly, and the conductive anomaly could determine whether there are fracture zones within the weak reflection areas. The simultaneous surveys can complement each other for detecting fracture zones. Besides this, such joint detection may be used to estimate whether there exist aquiferous fracture zones at depth.


Introduction
High-level radioactive waste (HLRW) mainly refers to highlevel liquid waste and solidifying objects generated by postprocessing of nuclear fuel. The geological disposal method is considered safe, reliable and technically feasible for disposing of HLRW (Wang et al. 2006). Aikin et al. (1977) gave four criteria that must be met before any crystalline rock body is used for a repository. In particular, the rock body must have very low in situ permeability because the fracture zones could be likely to result in the leakage of containment of radioactive wastes at depth. Geophysical exploration is the key step to check whether the preselected sites met such criteria for repository sites. Mair & Green (1981) identified fracture zones in the southern Canada Shield with a high-resolution seismic reflection survey profile across the central region of an Archaean granitic pluton, where it had previously been defined as relatively homogeneous by other data. On the island of Ävrö in south-eastern Sweden, three prominent dipping reflectors were observed and correlated with known fractures in boreholes using perpendicular reflection seismic profiles (Juhlin & Palm 1999). At Yucca Mountain (USA), Femina (2002) found a low electrical conductivity anomaly using  Chen et al. 2009). 1. Quaternary Holocene series alluvium sandy soil and gravel; 2. quaternary Pleistocene alluvium Gobi sand, gravel; 3. ancient Changchengian system biotite quartz schist; 4. Variscan medium-grained biotite monzonite granite; 5. Caledonian amphibolic diabase; 6. fault; 7. AMT survey line; 8. seismic survey line; 9. borehole and 10. study area. an electrical conductivity map, which correlated with active drainage ways and reflected lateral variations in allostratigraphy. At Sellafield in northwest England, which has a high level of cultural noise, Unsworth et al. (2000) detected a large zone of hypersaline groundwater extending 1 km inland toward the potential repository using CSAMT (controlledsource audio-magnetotelluric) data and indicated the effect of faults on the hydrogeology. Recently, some other geophysical methods have been developed to detect fracture zones in different fields, such as detecting overburden deformation with optical fiber sensing (Sun et al. 2018) and fracture-zone scope by microtremor signals (Du et al. 2019). However, a single geophysical method focuses on those applications.
Since the mid-1990s, many geophysical methods have been applied by the Beijing Research Institute of Uranium Geology (BRIUG) for determination of potential geological disposal sites of HLRW at Beishan (Gansu Province, China) and other places. The geophysical methods applied include AMT, MT, gravity, magnetic prospecting and so on. The AMT method has played an important role from the aspect of delineating the rock mass range, detecting fracture zones and so on (An et al. 2013;Xue et al. 2016). The preselected Beishan area is presently the most extensively studied in China. From 1990 to 1999, several favorable blocks were selected in this area by BRIUG. The Jijicao rock mass was mapped at a scale of 1:2000, mostly in 2006. In 2011, an AMT survey profile (AMT-northwest profile) was recorded through the BS15 borehole (figure 1) in the northwest direction and detected a remarkable low resistivity anomaly near the borehole. To study further the cause of the anomaly, a short seismic exploration profile (SEIS-northwest profile) along the same direction was shot through the BS15 borehole in 2012 (figure 1). This work provided an example of detection of such fracture zones in granite using audio-magnetotelluric (AMT) and seismic methods in China.

Geologic setting
The Jijicao rock mass covers an area of about 12 km 2 (figure 1) with the shape of long lenses in the northwestsoutheast direction, and which is 4-km long to the northwest and 3-km wide to the northeast. The land surface is between 1000 and 2000 m above sea level, where the geomorphological features are flat Gobi Desert and small hills. This rock block belongs to the Erdaojin-Hongqishan compound anticline of the Tianshan-Beishan folded belt in western China. The principal faulted structure in the block mass is a northeast-striking compression fault, with other compression and torsion faults. The outcrop in the Jijicao block is mainly composed of biotite monzonite granite in addition to a metamorphic residue body (the Proterozoic biotite schist) and Quaternary sediments (the Quaternary sediment includes valley type and low-lying land type) (Chen et al. 2009). The place where the AMT profile passed is mainly composed of monzonite granite (figure 1). Figure 2 shows the lithological column and well-logging curves of BS15. On the AMT profile at a distance of 1400 m (figure 1) and below the borehole BS15 (figure 2), there is mainly monzonite granite from the surface to a depth of 530 m. The Gamma Ray (GR) curve varies near 100 API except for some depths where there is higher GR and a spike anomaly is present due to radioactive enrichment by flowing water near the fracture zone. The low density, acoustic impedance and resistivity in the fractures at different depths demonstrate a decrease of the rock's elasticity. The resistivity curve presents a small electric contrast with thin serrated shapes that correspond approximately to the fractures at depths. In an S-shape, the amplitude of the resistivity curve starts falling gradually from the highest value at a depth of 100 m, then rises and falls again from 300 m deep to the bottom. There is metamorphic rock from 224.4 to 269.8 m deep in the borehole, the density of which is slightly higher than that of the surrounding granite. When the reflection coefficient curve exceeds 0.06, as calculated from interval transit time and density curves, the effective reflection can be recognised (Ludwig et al. 1970). Near depths (such as 208, 335, 370, 390, 418 and 515 m) with less density and speed than in the upper or lower regions, there is a significant contrast in acoustic impedance related to fractures or variation in the lithological composition.  figure 4; the discrete data are those observed; red for TM mode and blue for TE mode.)

AMT data acquisition and processing
AMT data were recorded by a Phoenix Geophysics V5 System 2000 TM with AMTC-30 induction coils. AMT data recorded over 25 min were sufficient with a tensor layout. There were 81 stations at intervals of 50 m along the profile.
The impedance of each station was calculated using a robust method with local magnetic field reference by SSMT 2000 software, based on over 100 cross-power spectra. In figure 3, the curves of apparent resistivity and phase were illustrated near the profile's beginning and borehole BS15. The apparent resistivity and phase data have high quality in the range of 8-1000 Hz while there are few data available above the frequency of 1000 Hz. This is because the natural source energy is very weak in that frequency range (Garcia & Jones 2002). This results in difficulties in the inversion process, as discussed later.
Because the nonlinear conjugate gradient (NLCG) algorithm (Rodi & Mackie 2001) has better computational speed and other advantages, the AMT data were inverted by this method using WinGLink TM software (Geosystem SRL 2008). The inversion model is composed of a series of rectangular cells: the shallow ones are of the same height, and below them the cell height increases with geometric progression.
The low quality of the impedance data above the frequency of 1000 Hz lowers the precision for delineating shallow abnormal bodies. In particular, the conductive anomaly close to the borehole will show a few differences when used with inversion models at different grid scales. The well-logging curves can help to determine the final geoelectrical model (figure 4), based on the characteristic of the resistivity curve, which decreases obviously at the bottom of the borehole. The geo-electric structures from transverse electric (TE) and transverse magnetic (TM) modes are generally similar, although with some differences. Although the TE mode is sensitive to good conductors, the resistivity values of the geo-electric section from the TE mode is lower than that of the TM mode and the spatial extent of the low resistivity from TM mode is smaller than that of the TE mode (the observed and model response pseudo-sections of TM and TE modes in figure S1and S2 in Appendix A). The geoelectric structure from joint modes resembles the real geoelectric structure more closely and weakens the prominent features of the single-mode TE or TM. The resistivity logging curve is merged in figure 4 (TE and TM mode). The varying shape of the resistivity curve (figure 4) is in accord with that of the section, especially from near 100 m depth (S with decreasing resistivity value) and at the bottom of the borehole (with the lowest value).

Geo-electric features
In figure 1, there are several faults, most of which are related to geo-electric models: TE, TM and TE and TM modes. Fault F2 intersects the AMT profile at a distance of 500 m. Along with the F3 fault direction toward the southwest, this fault may intersect the AMT-northwest profile at a distance of 900 m. Faults F4 and F6 intersect the AMT-northwest profile at distances near 3500 and 3900 m, respectively. At four such locations, conductive anomalies are present. The anomaly near F2 extends down over 200 m. Between F4 and F6, the model resistivity declines gradually from the surface with depth. Near the distance of 1400 m of the profile, there is a similarly conductive anomaly in the different models. Those anomalies of the TM and TE and TM modes extend below the depth of 500 m and are more conductive than the bilateral parts. Although the reason why the conductive anomaly is caused is related to fractures known from the BS15 borehole data, it is still difficult to understand the genesis; that is, what really causes that anomaly in the host rock of monzonite granite, especially its shape.

Seismic exploration data acquisition
The seismic profile was located between the distances of 733 and 1850 m of the AMT profile (figure 1). The shots were recorded by a Strata Visor Nzxp 24 of the US Geometrics Company at 1024 ms in length and 1-ms sampling interval. The correlation time length was 8 s with the pre-amplifier in the 'all pass open' mode. The first shot was located at a distance of 1900 m, and the first geophone was at a distance of 1720 m. Each trace used one immobile geophone and the trace interval was 5 m. The recording layout had 144 tracks that recorded synchronously. The shot spacing was 30 m with a total of 38 shots.
The car-borne Minivib T-15 000 (US IVI Company) was used as an energy source for excitation with scanning bandwidth from 10 to 400 Hz and 8 s scanning length. Such excitation was repeated 3-4 times with a maximum ground force of 2718 kg.

Seismic data analysis and processing
Seismic data analysis is the primary step of seismic data processing, which affects imaging quality, resolution and reliability of results. There are many kinds of disturbed wave in the shot seismic records due to the changes of near surface conditions and the influence of various complex factors in the acquisition process. A shot record (figure 5) shows the direct waves are not obvious due to the influence of the quaternary weathering layer on the surface within the small offset. The energy of surface wave and linear interference wave generated by near surface shallow medium is strong (from the first to . The reflected wave group is continuous from the 20th to 64th trace with high SNR. However, it is interrupted near the 40th trace, and the following reflected waves and surface waves are mixed. Between the 65th and 100th traces, the surface wave, random disturbance waves and effective waves interfere with each other with low SNR (yellow shaded area in figure 5). Between the 101st and 144th traces (figure 5, pink shaded area), the surface wave energy is weakened while the energy of high-frequency random noise wave is strengthened. In figure 6, the surface waves are concentrated in frequency components in the range of 20-40 Hz, the reflected waves in the range of 40-70 Hz and the high-frequency distance waves in the range of 150-200 Hz.
After general pre-processing of the raw data, common scattering-point (CSP) gather is gained with the scattering imaging method (Wang et al. 2007), and then the prestacking time profile is gained after velocity analysis (figures S3-S7 in Appendix B). The scattering image is a Kirchhoff pre-stack time migration for subsurface geologic bodies based on CSP gathers of equivalent offset (Bancroft et al. 1998). The method is composed of two steps: the first step is to form a CSP gathers, and weight amplitude with spherical diffusion correction factor, tilt correction factor and phase correction factor; the second step is to return the energy of the hyperbolic distribution to the scattering point through Kirchhoff principle. The time-depth transform is conducted after calibration with well-logging data 6 Figure 7. Scattering imaging time section, the green dot shows the BS15 well position. and lithological information (figure S8 in Appendix C). Because of the small topographic variation and bedrock outcropping, the elevation static correction amount can be ignored.

Seismic exploration results
When the incident wave comes along the normal incident direction and meets low-velocity and low-density media, it presents a negative reflection coefficient and represents a negative event.
In the time migration section (figure 7), the events near the BS15 borehole at the middle of the profile (below 100 ms, from 1150 to 1650 m) discontinue and significantly differ from the upper and laterals with low SNR and weak reflection energy. There is no outstanding reflection phase although the stack times are the highest.
There is a fault in figure 7 from (1000 m, 0 ms) to (1200 m, 150 ms) (see black line); as well as from (1300 m, 0 ms) to (1250 m, 150 ms). The bulge of the events near 100, 190 and 330 ms (green rectangle) in the middle part of the profile represent deformation because of underground 7 Figure 8. Joint exploration results. Left, depth section obtained by logging data and velocity analysis and right, 2D resistivity inversion section of the AMT data (part). rock compression. Due to adoption of the 5-m trace space to acquire high-resolution seismic data and to use of the welllogging data to calibrate a reflecting interface, the horizontal and vertical resolutions and the depth profile (left of figure 8) were obtained by time-depth conversion.

Discussion
According to the time-depth conversion results and geoelectrical models from the resistivity inversion, the following understanding was obtained: (1) Faults revealed. On the left side of figure 8, there is a fault from the surface at the distance of 1050 m toward the southeast to the depth of 500 m (magenta line). There is a spatial relationship between this fault and F3 in figure 1, and it is inferred that this results from extension of the F3 fault toward the southwest. This fault in the deeper area (below the depth of 500 m) is clearer, which is consistent with the well-logging data of the BS15 borehole. Moreover, F3 in figure figure 8 where the events represent continuous convex or concave shapes that are interrupted by weak reflections. The location where the conductive anomaly is outlined, is related to F3 (figure 1). The bent events usually are bound up with deformations. Accordingly, it may be inferred that this conductive anomaly may be the deformation that forms the fractured zones.
(3) The relationship between the geo-electric variance and the events. The purple rectangles mark where the geoelectric contours vary rapidly, and are projecting a timedepth profile as in figure 8. Where the resistivity varies more, the events always represent interrupted, weak reflection or sudden change. (4) The joint detection results may infer an aquifer at depth. The fractured granite may be aquiferous. These two pieces of evidence, the low resistivity anomaly detected by the AMT method and clear event variation detected by seismic exploration, may be conducive to inferring an aquiferous fracture.

Conclusions
It is one of the key steps for the selection of a HLRW repository to investigate the fine geological structure underground. From the joint geophysical surveying example near the BS15 borehole in the Jijicao block, it is clear that the conductive anomaly detected by AMT results from sub-horizontal compression deformation based on seismic surveying; moreover, it is demonstrated that joint surveys could complement each other significantly. The weak reflections imaged by the seismic method can be interpreted as the fractured zone, unitary lithology or an energy problem, but it is hard to identify which one causes the weak reflection using only the single seismic method. However, the conductive anomaly detected by AMT is always related to the fractured zones. Therefore, with the joint survey, it is possible to determine whether there is a fractured zone underground.
It is difficult to determine what causes the conductive anomaly using only the AMT method or other resistivity methods, although its shape and extension can be outlined roughly through the geo-electrical section. While the seismic method can image whether the terrain is deformed, the conductive anomaly may be related to a fractured zone if the terrain was deformed in the seismic image.
For instance, between the distance of 2400 to 3300 m in the AMT profile, there exists another remarkable conductivity feature of which the dimensions are much larger than that near the borehole. However, based only on the geo-electrical model, it is difficult to make a reliable judgment about which, choosing between alternation and compressed fracture, is associated with this anomaly. of HLRW in China' funded by SACTIND (the State Administration of Science, Technology and Industry for National Defence, China), NSFC (grant no. 41790443 and grant no. 41641040). Professor Yuanxin Jin and Professor Weiming Chen of CNNC key laboratory on geological disposal of HLRW provided geological data and help. With help from engineers Yunfeng Li and Ming Zhang, the field work was more convenient. Last, special thanks are offered to Professor Handong Tan of CUGB, who gave critical suggestions and improved the paper. Specially, the authors thank the Editors and three anonymous referees for their constructive reviews that have greatly improved the manuscript.

Appendix A. The pseudo-section of the observed and response data of all AMT sites
The pseudo-sections of all sites at the AMT-NW profile are shown in figure S1 (apparent resistivity) and S2 (phase), which also present the response pseudo-sections. The observed data near 1000 Hz present poor quality for the weak energy of the natural source, which increases the difficulties in the inversion process. With a different inversion model grid, the final models are gained with a rms of 15.8 after a total of 200 iterations following three restarting.

Appendix B. Common scattering point (CSP) gathers and velocity spectrum
In this paper, the CSP gathers are used to extract migration velocity information. The purpose of forming CSP gathers is to increase the number of coverage and the signal-to-noise ratio. Figures S3-S6 show the velocity analysis for CSP gathers in No. 51, 101, 181 and 321. Figure S7shows the velocity field of scattering imaging.

Appendix C. Calibration
Seismic horizon calibration is a bridge connecting seismic and logging data, and plays an important role in reservoir parameter predictions such as seismic inversion. When the wavelet extraction and horizon calibration are more accurate, a prediction result with high accuracy is obtained. The result of calibration directly affects the corresponding relationship between seismic reflection layer and geological horizon. Data is used from from sonic logging and density logging for a synthetic seismic record to calibrate the horizon and establish the corresponding relationship between the seismic events and the layers, or lithology. The calibration is shown in figure  S8. The reflection of metamorphic rock and fracture zone is in good agreement with logging.