Geothermal resource exploration by using Numerical simulation and comprehensive geophysical method

： 8 Geothermal energy is an important renewable clean energy resource with high development and usage potential. 9 Geothermal resources, on the other hand, are buried deep below, and mining hazards are significant. Geophysical 10 investigation is frequently required to determine the depth and location of geothermal resources. The Transient 11 Electromagnetic Method (TEM) and the Controlled Source Audio Frequency Magnetotellurics (CSAMT) have the 12 highest detection efficiency and accuracy of all electromagnetic exploration methods. This article initially explains 13 the algorithm theory of the finite difference technique before establishing a simplified geothermal system resistivity 14 model. Established on the simplified resistivity model, a simulation analysis of the ability of CSAMT and TEM to 15 distinguish target body faults at different resistivities and dip angles was performed, and the effectiveness and 16 difference of the two methods in detecting typical geothermal resource targets was verified. A complete exploratory 17 research of CSAMT and TEM was conducted in Huairen County, Shuozhou City, Shanxi Province, China, based on 18 theoretical analysis. Both approaches can reflect the geoelectric structure of the survey region, demonstrating the 19 efficacy of the two methods in detecting genuine geothermal resources. 20 21 structure. This study verifies the respective advantages and disadvantages of CSAMT and TEM in geothermal resource 328 detection.


22
As a green, low-carbon clean energy, geothermal resources have huge reserves and are widely distributed in 23 China. However, currently utilized geothermal resources in China are less than 5‰ of the reserves, and further 24 development and utilization of geothermal resources have broad prospects (Shi B, 2005). The buried depth of 25 geothermal resources is generally large, most of which have exceeded 2000m, and the risks taken during mining are 26 great. In order to have a clearer understanding and understanding of geothermal resources, it is necessary to conduct 27 a geological survey of geothermal resources before mining to find out the depth, location and reserves of geothermal 28 resources. In the exploration of geothermal resources, geophysical methods are an effective method (Wang G et al.,29 2000; Yang L, 2017; Kang F, 2018). 30 Among the geophysical methods, Controlled Source Audio-frequency Magnetotellurics (CSAMT), Transient 31 Electromagnetic Method (TEM) are relatively effective methods for the exploration of geothermal resources, but 32 each method has its own limitations (Muraoka H et al., 1998;Spitzer K, 1995 However, most of the past geothermal resource exploration focused on engineering applications, and there were 37 few geophysical numerical simulation studies, that is why in this paper we focused to correlate the numerical 38 simulation with geophysical field data applicability. The ability of taking topography into account in a forward solver 39 greatly depends on the numerical method. The finite difference approach used to be a popular method in the past but 40 it is restricted only to be used in areas with flat topography (Spitzer K, 1995;Mufti IR, 1976). Alternatively, the 41 integral equation method is usually valid for simple structures (Dieter K et al., 1969;Okabe M, 1981). Including 42 complex topography and substructures in the modelling remains difficult with such methods. The finite element 43 approach was first introduced in the resistivity problem by Coggon (1971) and later further developed by many other 44 authors. It appears to be a flexible method to model both topography and possibly complicated geometries. Initially, 45 finite element modelling codes including topography were defined on distorted grids, using quadrangular or 46 hexahedral elements, and were thus limited to undulating topography .These approaches using block-oriented 47 meshes are also limited regarding the possibility of effective mesh refinements. Only recently, codes have been 48 developed on unstructured tetrahedral meshes; allowing for complex geometries and local mesh refinements, 49 especially around source locations. In parallel, surface integral methods are also available, with the advantage of only 50 discretizing the top interface 51 Independent of the numerical scheme, a resistivity forward solver consists of solving the electrical boundary 52 value problem for every source location. The solution of such a system of linear equations shows a singularity at the 53 source location. This leads to significant numerical errors in the vicinity of the electrode position. The singularity 54 removal technique introduced to splits the potential into a singular part, related to the source, and a residual part 55 taking into account a possible non-homogeneous conductivity model. The singular part is defined analytically for flat 56 topography. For non-flat topography, the singular component is usually computed numerically. Even with local mesh 57 refinement around the sources, a highly accurate solution is however difficult to be obtained as the objective is to 58 estimate a singular function (Günther T et al., 2006). With a thorough description and history of TEM and its 59 modeling, we discovered that TEM has never been utilized for geothermal research, and if it has been used, it has 60 never been coupled with numerical simulation. This research will demonstrate the relationship and significance of 61 numerical simulation and physical investigation in their respective fields of study. Furthermore, it is not required in 62 any investigation that numerical simulation produce equally exact and matched results with field data, but it does 63 provide a hint for future exploration in the field that is the subject of this study. We primarily simulate and evaluate 64 the anomalies induced by changes in the resistivity, inclination, and width of the low-blocking layer in the simplified 65 model of the geothermal system's resistivity using CSAMT and TEM. Evaluate the anomalies induced by changes in 66 Therefore, a uniform time grid is sampled, and the sampling interval is half a time step. When the grid form and time 86 sampling method are adopted, equation (3) Since the excitation source current is in the xoy plane, there is no z-direction component of the source current, 89 and only Jsx and Jsy components exist in the equation. 90 Furthermore, the difference scheme of the electric field FDTD iteration in the active medium. 91 The excitation source loading mode is only related to the x or y component of the electric field, so the iterative 94 formula of Ez is the same as that of the passive region. 95 The radiation conditions are the regional boundary conditions, and the converging boundary conditions are the 96 ground-air boundary and the stratum and anomalous body boundary. The connection criteria between the formation 97 and the aberrant body are automatically met in the Yee grid cell computation of FDTD, and no extra treatment is 98 required. Therefore, the ground-air boundary conditions are mainly considered here. 99 The truncated boundary of the underground space is easier to meet the conditions of the far field, but it is not 100 easy to meet the air area. In addition, there is a large difference between the aerial grid step length and the 101 underground grid step length. According to the principle of canceling the scale, a considerable part of the calculation 102 will be used in the non-key investigation area of the ground M-TEM. Therefore, a special solution is required. 103 In the quasi-steady state of the geophysical time-varying field, the air field satisfies the Laplace equation as: 119 The Fourier transform on the ground-air boundary uses the fast Fourier algorithm. In order to meet the 120 requirements of discrete Fourier uniform grids, different processing methods are adopted according to the 121 characteristics of grid division 122 1) Non-uniform grid. The uniform grid is interpolated by nonlinear interpolation algorithms such as cubic spline 123 and five-point cubic smoothing, and after fast Fourier transform and inverse transformation continuation, the 124 interpolation method is used to restore the original grid again. 125 2) Non-uniform grid. If the nonlinear interpolation has excessive fluctuations that affect the calculation accuracy, 126 you can change to the linear interpolation trial calculation. 127 3) Uniform grid. Fourier transform and inverse transform can be directly performed without interpolation, and a 128 grid that is more suitable for discrete Fourier transform can be obtained by linear interpolation according to 2n 129 interpolation. 130 In brief, when there is an issue with simulation accuracy, grid interpolation of the discrete Fourier transform, as 131 a component of simulation software, can be used to enhance accuracy and identify a factor that influences accuracy. 132 Through the iterative solution and the relationship between the electric field and the magnetic field, the response of 133 the electromagnetic field at each time can be obtained. 134

135
According to the relative relation of the resistivity of the basic elements of the geothermal system, a simplified 136 resistivity model of the geothermal system is designed for the common faults in the geothermal system that show low 137 resistance due to water conduction. 138 The numerical simulation program is the open source code SIMPEG for writing papers on the Geophysical 139 Tablet of the University of British Columbia, and the method is the finite difference method. addition, in order to simulate the measured data as much as possible, Gaussian noise with a mean value of 3% of the 153 response amplitude is added to the response obtained by the numerical simulation. In order to compare with the 154 subsequent transient electromagnetic simulation more intuitively, noise is added to the electromagnetic field 155 response obtained by the numerical simulation, and then the apparent resistivity is further calculated. Following are 156 the parameters upon which our modeling is based on. 157

A) CSAMT Response to the faults with different resistivities 158
In the simplified resistivity model of the geothermal system established, the resistivity of the fault is 10Ω•m. 159 usually in the stratum half-space, a fault can be regarded as an angularly inclined plate-like body extending 160 downward from the ground. There are a large number of cracks in the fault zone, which are filled by particles such as 161 breccia, and the resistivity of the filling is different, which determines Characteristics of fault resistivity. Therefore, 162 the resistivity of the fault is generally within a large range. We establish fault models with different resistivity values 163 (fault dip angle 30°, width 70m, and fault resistivity values 10, 50, 100Ω•m), and perform forward simulation to 164 obtain each response, then calculate its resistivity. All responses are shown in the Fig. 2 (a-c). 165

B) CSAMT Response to the faults with different dip angles 166
The dip of the fault is usually an important indicator of detection. Based on this, we establish fault models with 167 different dip angles (fault resistivity value is 10Ω·m, width is 70m, dip angles are 30°, 45°, 60°) and forward 168 modeling is performed, as shown in the Fig. 3(a-c). 169

170
It is discussed above that, based on the numerical calculation of low-blocking layer models with different 171 resistivities and inclination angles; the ability of CSAMT to identify low-blocking layers was analyzed. As a 172 comparison, this section will carry out a numerical simulation analysis on the ability of the TEM method to identify 173 low-blocking layers. Similarly, starting from the resistivity and inclination of the fault, analyze its influence on the 174 TEM response. Because the apparent resistivity of the TEM method is divided into early and late stages, it is difficult 175 to obtain a stable apparent resistivity calculation formula covering all time channels. Therefore, the TEM analysis is 176 carried out for the derivative of the vertical magnetic field with respect to time (dB/dt). (2) CSAMT has stronger anti-interference ability than TEM, and comprehensive exploration is adopted in the 206 section with strong electromagnetic interference signal, which can improve the accuracy of data collection ( Fig. 6. 230 As shown in Fig. 7, four CSAMT survey lines are arranged in the exploration area. 231 In the CSAMT detection, the instrument used is the GDP32Ⅱ multifunctional electrical method workstation. 232 The transmitting pole distance AB=1500m, the transmitting current is 14-16A, the transmission distance is 6km; the 233 receiving point distance of line 59 and 64 is 50m, the receiving point distance of line 10 and line 60 is 100m, and the 234 signal frequency range is 0.125-8192Hz.In the TEM detection, the instrument used is also the GDP32Ⅱ 235 multifunctional electrical method workstation, and the center loop device is used for measurement. The transmitting 236 wire frame is a 600m×600m single-turn loop, powered by a generator, the fundamental frequency of the transmitting 237 source is 16Hz, and the transmitting current is 15A; the distance between the measuring points is 50m. After the 238 CSAMT measurement is completed, a TEM measurement line is arranged near the 60 line where the CSAMT 239 resistivity is more obvious, which is used to more accurately delineate the low resistance fracture zone, and combined 240 with the CSAMT and TEM results. 241

242
Different rock formations have different conductivity. Generally speaking, in most igneous rock regions, due to 243 the ancient diagenesis age, after long-term metamorphism, the formation is compact and complete, and the fractures 244 are not developed. The resistivity is relatively uniform in the horizontal direction, and gradually increases with the 245 increase of the depth in the vertical direction, which is not conducive to Looking for geothermal resources (Li H, 246 2010). 247 In some igneous rock areas, the pore water of the overlying loose layer, after receiving the vertical infiltration 248 replenishment of atmospheric precipitation and the intermittent leakage replenishment of surface water, can flow into 249 the deep igneous rock and the fissures of the igneous rock basement through the fault fracture zone and rock fissures. 250 Thermal storage aquifers form underground hot water and are stored in igneous rock formations. The resistivity of 251 eruptive igneous rocks (such as basalt, rhyolite, etc.) is generally low (20-500Ω•m), which is lower than that of 252 granite (greater than 500Ω•m), but the resistivity of igneous rocks is generally greater than that of sediments 253 Surrounding rock. Therefore, the difference in resistivity values between igneous rocks and sedimentary surrounding 254 rocks can be measured by geophysical prospecting methods to determine the location of igneous rocks and igneous 255 rocks, and then to search for geothermal resources (Zhang Q et al., 2016). 256 On the other hand, in areas where faults are developed, the resistivity value will change. The activity of the fault may 258 cause a fracture zone near the fault, and the fault will lead to the hydraulic connection between the upper and lower strata, thus 259 forming a low-resistivity anomaly area near the fault that is close to the dip angle of the fault. Such low-resistance anomalies 260 are in the resistivity. Generally, the cross-sectional view is relatively steep. By looking for such a steep anomaly with low 261 resistivity, it can reflect the existence of faults, and then search for geothermal resources (Li W et al., 2002;Garyet G et al., 262 Verified situation 264 In the processing of the measured data, we first use the data preprocessing software of GDP32Ⅱ to sort the collected data 265 and remove the dead pixels, and then use the CSAMT-2D software and TEM-1D software based on the OCCAM algorithm to 266 invert the data. The inversion results of CSAMT and TEM resistivity profiles are shown in Fig. 8(a-c). 267

268
The results show that the influence of the fault on the CSAMT response is mainly in the low frequency range, and the 269 influence on the TEM response is mainly in the late stage. Because CSAMT low frequency band can still obtain high 270 signal-to-noise ratio data, and the signal-to-noise ratio of the late TEM response is generally lower, so CSAMT has a stronger 271 ability to distinguish low-blocking layers at this depth than TEM, and CSAMT can also distinguish faults. Upper and lower 272 plate. On the other hand, the TEM response can clearly distinguish the abnormal response caused by the width of the fault, but 273 the CSAMT response cannot be distinguished, which shows that the TEM has a stronger ability to describe the details of the 274 target body in the fault with low resistivity than the CSAMT; In addition, the faults are mainly manifested as obvious 275 resistivity dislocations on the apparent resistivity section of CSAMT, and the resistivity characteristics of the fault fracture 276 zone cannot be directly judged through the apparent resistivity section. 277 Although both methods can identify the low-blocking layer of the model, the two methods show obvious 278 complementarity in the detection depth and detail description ability, and are suitable for comprehensive application in the 279 detection of geothermal resources. 280 In general, the morphology of the resistivity profile of CSAMT and TEM is basically the same, and both show that the 281 resistivity is medium-low resistance in the medium and shallow layers, and the resistivity gradually increases as the depth 282 increases, and the deep substrate shows high resistance. 283 Comparing the detection effects of CSAMT and TEM, the differences are: 284 (1) TEM is better than CSAMT in the detection effect at shallow depths on the surface. CSAMT is affected by topography, 285 the low-resistance shielding interference of the Quaternary muddy sand and clay layer is relatively large, which shows 286 abnormal medium resistance, and the resistivity curve presents characteristics such as distortion, which reduces the resolution 287 of shallow layers to a certain extent. Although there is a blind zone with a depth of about 100m in TEM, the low-resistance 288 shielding layer at the shallow surface is less interference, and the resistivity value shows obvious regularity at the shallow 289 surface, which is consistent with the actual geological characteristics. 290 (2) CSAMT is better than TEM in the detection effect at large depths. In the 4500-10000 point section, CSAMT has 291 better data stratification in the deep part, while the TEM resistivity curve in this section is slightly confused, and the TEM 292 shows a circle of the resistivity curve in the 9000-9500 point section. Closed high-resistance abnormal value. This 293 high-resistance abnormal value should be a false abnormality caused by interference. The secondary field potential of the late 294 TEM measurement track also fluctuates greatly in this section, lacking regularity, indicating that CSAMT has better 295 anti-interference ability TEM. 296 (3) TEM and CSAMT have their own advantages in the detection of water-conducting fault structures. 297 Since the pure TE field mode of TEM is especially sensitive to low-resistance targets (Xue et al. 2013), this is more 298 obvious in the reflection of F1 fault. CSAMT hardly reflects the F1 fault in the low-resistance area in the shallow part, but there 299 is obvious low-resistance anomaly on the TEM profile. It can also be seen from the TEM secondary field potential 300 multi-channel map that there is an obvious abnormal high value of the secondary field potential in this area. 301 The advantage of CSAMT for fault detection is mainly reflected in the detection of deep high-resistance basement 302 interruption layer. In the sections of 5000m-5700m and 9500m-10000m, it is inferred that the buried depth of the Quaternary 303 and Tertiary loose deposits is about 500m. The lower part is the basement of igneous rock, and the resistivity reflects the 304 characteristics of high resistance. In the CSAMT cross-section map, these two sections appear as steep gradient zones, which 305 are inferred to be F2 and F3 faults, which can more intuitively distinguish the upper and lower walls of the fault. 306 (4) CSAMT is better than TEM for the detection effect of igneous rock basement. In the 5800m-9300m section, CSAMT 307 has a very obvious high-resistance response to the igneous rock basement, while the high-resistance response of TEM is not 308 obvious at this position, and it does not highlight the igneous rock basement. 309 After the completion of the geophysical prospecting construction, the drilling verification was carried out at 3350 point. 310 The borehole encountered underground hot water at 1610m; then a pumping test was carried out. According to the results of 311 the pumping test, the unit output of underground hot water in this area was 233m3/d. The water temperature is 58°C. This with different resistivity values, it can be seen that when the fault resistivity is small (10, 50 Ω·m, less than the resistivity of the 418 first layer of medium), the apparent resistivity curve near the fault is relatively smooth; As the fault resistivity increases, the 419 apparent resistivity curve rises more and appears to be significantly higher. Generally speaking, the difference in resistivity in 420 the fault zone can have obvious characteristics in the apparent resistivity profile. 421 422 Fig. 3 The model of fault dip angle is 30°, 45°, 60° and the cross-section view of resistivity from the above analysis, it can be 425 seen that when the dip angle of the fault is 30°, 45°, 60°, "sags" will appear in the interrupted layer of the apparent resistivity 426 profile; The greater the dip of the fault, the greater the degree of "sag" in the apparent resistivity profile; the dip reflected by the 427 apparent resistivity in the apparent resistivity profile is greater than the dip of the actual fault. 428 429 430 431 Fig. 4 The model of fault resistivity is 10, 50, 100 Ω•m and the cross-section view of resistivity It can be seen from the above 432 analysis that when the resistivity of the fault zone is low, there is an obvious anomaly near the center of the fault, and the 433 maximum value of the anomaly is located at the midpoint of the fault. The lower the resistivity, the more obvious the abnormal 434 response. When the resistivity of the fault zone is 100Ω•m, the resistivity difference between it and the surrounding rock is 435 small, and no large anomaly can be generated on the resistivity profile. This shows that when the resistivity difference between 436 the fault and the surrounding rock is small, TEM is difficult to distinguish the existence of the fault. 437