Prediction of Conductive Anomalies Ahead of the Tunnel by the 3D-Resitivity Forward Modeling in the Whole Space

School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China Institute of Geophysics and Meteorology, University of Cologne, Cologne 50969, Germany CAS Center for Excellence in Comparative Planetology, University of Science and Technology of China, Hefei 230026, China National Geophysical Observatory at Mengcheng, Mengcheng 233500, China School of Resources and Environmental Engineering, Hefei University of Technology, Hefei 230009, China


Introduction
In recent years, with the development of high-speed railways and urban subways, it is often necessary to excavate tunnels in complex mountain areas. The geological survey is an indispensable task before the tunnel excavation. Due to the limited space of the tunnel, the prediction of geological disasters becomes very difficult and inaccurate, especially water inrush disaster, induced by those water-rich low resistivity structures, easily causes a lot of casualties and economic losses [1][2][3][4][5][6]. In addition, global mining coal resources have been exploited from shallow land to deep land, and the mining environment and construction have become more complex. The prediction of potentially disaster-producing water inrush structures ahead of the tunnel has become an urgent problem to be solved in the safe production of the coal mine [7].
Many geophysical methods in predicting geological disasters ahead of the tunnel have been developed in the last decade. The tunnel reflection tomography (TRT) is applied in geological advanced prediction ahead of the tunnel excavation [8]. The transient electromagnetic method (TEM) is also used for the mining industry [9,10]. The resistivity method for 3D anisotropic modeling is proposed to predict water-rich faults ahead of the tunnel [11]. Joint inversion of TEM and the resistivity method is investigated about advanced detection of the tunnel using the particle swarm optimization technique [12]. However, the DC resistivity method has an obvious advantage in that is logistically easy to apply in field prospecting compared to TRT and TEM.
The resistivity method is an effective way widely used in detecting water-bearing structures [7].
The resistivity method for forward-probing in tunnels has been developed rapidly in a variety of ways. Cheng et al. [13] proposed the basic principle and application of forward-probing in tunnel excavation by the resistivity method. Huang et al. [14] discussed the anomalous character of large polar distance in the tunnel using the resistivity method on the dipole-dipole array. Han et al. [15] presented the primary theory about the mine resistivity method to predict the sphere anomaly structure. Liu and Wu [16] investigated the anisotropic resistivity for advanced detection in tunnels. In this paper, we present the forward-probing prediction of a low-resistivity anomaly ahead of the tunnel by the DC-resistivity method.
In the DC-resistivity method, the structure and the geometry of the tunnel is the most important factor affecting the detection of the possible conductive anomaly ahead of the tunnel face, leading to the reduction of the signal-to-noise ratio [17]. To improve the detection, it is necessary to include the influence of the tunnel in the modeling because the tunnel is a high-resistance air cavity, which can distort the distribution of the stable DC-current field [18].
There are few prediction models to detect the conductive anomaly ahead of the tunnel face based on the tunnel effect in the existing literature about the application of the DC resistivity method in mines [16]. The tunnel effect was discussed [19]. However, it is based on a simple physical experiment and could not explain some phenomena of the field data. In this study, we propose the Monte Carlo method to test the quality of our modeling results. This paper is aimed at developing a 3D-DC model by considering the tunnel effect for forward-probing in tunnels to effectively forecast the possible water-rich structure ahead of the tunnel. A 3D finite element (FE) resistivity modeling for a whole-space tunnel model using the unstructured mesh is developed. We propose a predicted equation by simulating different tunnel models to locate the water-rich anomaly ahead of the tunnel face. Lastly, the quality of the predicted model is tested and evaluated by the Monte Carlo method in comparison with an available prediction model [20]. Considering the tunnel effect is of great significance to distinguish the actual geological anomalies and distortion anomalies [21].
The paper is structured as follows. Section 2 introduces the boundary and variational problems of the total potential, secondary potential approach, and the verification of the 3D DC-whole space algorithm. Section 3 firstly investigates the effect of different tunnel size and then proposes a predicted equation for a cube anomaly by a large number of simulations. Section 4 evaluates the quality of the predicted equation by comparing the results of our model and another prediction model using the Monte Carlo method. Section 5 is the conclusion.

Numerical Scheme
2.1. The Boundary and Variational Problems of the Total Potential. The total potential u for the basic equation and boundary value problems [22,23] can be expressed as ∇ is the Laplace operator. σ represents a conductivity, which is the reciprocal of the resistivity (σ = 1/ρ). r is the observation location, I is the current strength, δ is the Dirac delta function, r 0 is the current source location, Γ 0 denotes the Neumann boundary, Γ R is the mixed boundary, and n is the normal direction of the boundary.
By Galerkin's method [24], the boundary problems for Equation (1) can be described as ð2Þ Ω is the target. Equation (2) can then be separated as [25] where Ω e is the volume element [25]. Using the interpolation function u = U T e N [26] where N is the shape function, U e is the node potential of the element. Equation (3) is substituted as K 1 e is the volume element matrix, K 2 e is the face element matrix, and N e is the total number of element. To reveal their composition, K 1 e is a third-order matrix as an example. It reads [27] Assembling all local system matrices into a global system matrix [26], we get QU = P: ð6Þ Q = ∑ðK 1 e + K 2 e Þ represents the global system matrix, U represents the unknown potentials, and P is the source term.

Secondary Potential Approach.
To overcome the singularity of the source, we consider a secondary potential approach [28] by splitting the total potential u into primary potential u p and secondary potential u s [16].
The secondary potential for the basic equation and boundary value problems can be written [26] where σ s = σ p − σ, σ p is the background conductivity or primary conductivity, and σ s is the secondary conductivity or anomalous conductivity. Similar to Equation (1), the secondary potential for Equation (8) can be derived by Galerkin's method δF u s ð Þ = 0: By repeating the same procedure as in Section 2.1, a global system of linear equations with respect to u s can be derived: U p can be obtained by an analytical solution [20]. The conductivity σ in Equation (9) is replaced by σ s to calculate Q ′ . The matrix Q has a sparse property. Therefore, the conjugate gradient method [29] is a very efficient iterative method for solving Equation (10). Visual Studio 2015 and Parallel Studio XE 2018 provide the compilation environment of Fortran program for 3D-resitivity forward modeling in the whole space. The following numerical calculations were performed on a computer with Intel(R) Core (TM) i5-5350U CPU at 1.80 GHz and 8 GB RAM.

Algorithm Verification.
The analytical formulas of the potential u and apparent resistivity ρ s for a whole-space slab [20] using the pole-dipole configuration are given as where a is the thickness of the slab, d is the distance from the current source A to the slab, z is the distance between the current source and the observed point, ρ 1 denotes the resistivity of the background rock, and ρ 2 denotes the resistivity of the slab. K AMN = 4πAM•AN/MN is the geometrical factor for the pole-dipole array. ρ s is the apparent resistivity of pole-dipole array.Δu MN is the potential difference of observed points M and N. I is the current intensity.
To validate the accuracy of our algorithm, the wholespace slab model ( Figure 1) with an analytical formula [20] is used for comparison. The thickness and the horizontal distance of the slab anomaly in this model are considered as 2 m. The resistivities of the slab and the surrounding rocks are assumed as 10 Ωm and 100 Ωm, respectively. A poledipole array is used where the source electrode A with a current intensity of 1A is fixed in the origin and observed potential electrodes are distributed at intervals of 2 m on the left side of the current source. The distance of potential electrodes M and N moving along the opposite direction of the y-axis is 2 m. Apparent resistivities are calculated using our 3D DC-whole space algorithm for the model displayed in Figure 1. In addition, analytical apparent resistivities are calculated using Equation (11). Both curves as well as their relatively small errors are shown in Figure 2. A good agreement (maximum deviation of less than 0.3 percent) is achieved between analytical and calculated apparent resistivities indicating that our algorithm is feasible. In addition, there is a minimum value on the sounding curve as the electrode distance increases, which can be used for predicting the position of the conductive water-bearing anomaly ahead of the tunnel face, because a linear relationship between the actual distance of anomaly and the offset corresponding to the minimum on the sounding curve exists [16,20].

Tunnel Model.
To study the influence of the tunnel on the apparent resistivity curve, we introduce a tunnel factor α to show the degree of the tunnel impact [14,21]. It reads α denotes the tunnel factor of tunnel effect, ρ tunnel denotes the observed apparent resistivity in the FE simulation for the tunnel as an only high resistivity anomaly in background rock, and ρ l denotes the background rock resistivity. To remove impacts of the tunnel, we use the tunnel factor α to obtain a modified apparent resistivity ρ c . Here, ρ s is the calculated apparent resistivity in the FE simulation under the condition of the water-bearing anomaly and the tunnel interference. ρ c denotes a modified apparent resistivity. A 3D resistivity model in whole space on unstructured tetrahedral meshes using finite element method (FEM) for different tunnel sizes is simulated in Figure 3. The tunnel face areas are 2m × 2m, 3m × 3m, and 4m × 4m. Tunnel resistivity is 10 8 Ωm and the background surrounding rocks resistivity is 500 Ωm. The depth of the tunnel is 500 m. 90 potential electrodes in pole-dipole array are distributed with an interval of 2 m, and the distance of potential electrodes M and N is also 2 m. The current source electrode A is fixed. An open-source software Gmsh can be used to produce unstructured meshes of the model [30]. Unstructured meshes have been proven effective with high accuracy and computation efficiency [26]. In addition, the local refinement strategy of the measurement electrodes in the tunnel has an advantage for achieving a balance between calculation time and solution accuracy [31,32].
The numerical results by the pole-dipole array in Figure 4 show that the effect of the tunnel on apparent resistivity curves is related to the geometric size of the tunnel. The effect is the most serious near the tunnel face, and the greater the length of the cross-section is, the more serious the effect of the tunnel becomes. With the increase of the electrode distance, the curve began to become flat at approximately 10 m away from the tunnel face and gradually up to be a stable value (α = 1). At this electrode space, the effect of the tunnel can be ignored, and we can think that the tunnel does not affect the DC-current field distribution behind a very large electrode distance. In future prospecting applications, the data of large electrode offsets are generally accurate without the tunnel effect while the data of small electrode offsets should not be used as much as possible except for the nearby anomalies. If the water-bearing anomaly is very close to the tunnel face with a small distance, we have proposed a prediction model in the following study to forecast accurately the position of the anomaly ahead of the tunnel in the FE simulation by considering the tunnel effect.

Geofluids
probing in tunnels, we design a cube anomaly model ahead of the tunnel in Figure 5, which is the most realistic case in the field surveys. The tunnel face area is 2 m × 2 m, 2 m × 3 m, and 3 m × 3 m. Tunnel resistivity is 10 8 Ωm and the background surrounding rocks resistivity is 500 Ωm. A 10 m × 10 m × 10 m cube anomaly with a low resistivity of 20 Ωm is considered at a distance of d = 10 m. The center of the anomaly and the tunnel with a depth of 500 m is at the same level. The current source electrode A is fixed in the origin of the y-axis and 90 potential electrodes with MN = 2 m are distributed at intervals of 2 m behind the tunnel face. Because of limited tunnel dimensions, a pole-dipole array is generally used to measure the electric potential and to derive the corresponding apparent resistivity curve. The minimum value of the apparent resistivity curve indicates that there is an anomaly ahead of the tunnel and the electrode distance of minimum point has a linear relationship with the actual distance of the anomaly, which can predict the position of the low-resistivity anomaly ahead of the tunnel face.
The apparent resistivity soundings for the tunnel area s = 4 m 2 , 6 m 2 , and 9 m 2 at d = 10 m are shown in Figure 6. The corresponding modeling results for d = 6 m, 8 m, and 14 m showed a similar pattern; therefore, they are not shown here. Firstly, in the FE simulations without the water-bearing anomaly ahead of the tunnel, models only including highresistivity tunnels with different dimensions were simulated to get different tunnel factor curves. Then, models by adding the water-bearing anomaly in front of the tunnel in the FE simulation were used to calculate the observed apparent resistivity curves. Lastly, the corresponding tunnel factors were introduced to correct the observed apparent resistivities to derive the modified apparent resistivities, which can highlight the low resistivity anomaly ahead of the tunnel face more accurately. From Figure 6 we can see that the electrode distance (x min ) to the apparent resistivity minimum point in curves becomes larger with the increase of the tunnel area s for the same actual distance d. The x min is highly correlated with the tunnel area and the actual distance, which can be used for the prediction of the anomaly location. It is noteworthy that this prediction method using the minimum value on apparent resistivity curves is only valid if there is no conductive interference structure beneath the tunnel.
Generally, the shape of the tunnel is circular [33][34][35]. But the square tunnel is also reasonable, which was often simulated for the tunnel prediction in many journals [14,20,[36][37][38]. To investigate the influence of the tunnel shape on the prediction of anomalies ahead of the tunnel, we used a circular tunnel model in Figure 7 to compare it with a square tunnel model in the simulation. For the same tunnel area of s = 12:56 m 2 , apparent resistivity curves of the circular tunnel and square tunnel are very consistent in Figure 8(a). Their maximum relative error is less than 0.25% and the relative error trends zero as the electrode distance increases in Figure 8(b). We can think that the tunnel shape could not affect the results when the tunnel area keeps unchanged (e.g., the square tunnel equals approximately the circular tunnel). For greater convenience in geometric modeling, we prefer the square tunnel.
To investigate the possible influence of the background rock resistivity on the prediction of anomalies ahead of the tunnel, we used the square tunnel model (s = 4 m 2 ) in Figure 5 with different resistivity of background rock (e.g., 200 Ωm, 300 Ωm, 400 Ωm, and 500 Ωm). As shown in Table 1, the electrode distances of the minimum value in apparent resistivity curves are the same regardless of the difference of background resistivity. The value of background resistivity only affects the minimum value, which is independent of the electrode distance. Figure 7 shows the simulated results of the above anomaly model at the tunnel area s = 6 m 2 , 9 m 2 , and 12 m 2 and at d = 6 m, 8 m, 10 m, and 14 m: The black geometric points represent the numerical results of the electrode distance x min to the minimum point in curves and the lines represent the corresponding fitted curves. The different data of the electrode distance x min to the minimum point, the tunnel area, and the real distance of the anomaly are calculated by 3D-resitivity forward modeling in the whole space, and then they are analyzed and fitted by the multiple linear regression method using an 6 Geofluids open-source software Origin2019b. A formula is proposed to predict the distance of the low-resistivity anomaly ahead of the tunnel as

A Proposed Equation for Prediction.
where d pre represents the prediction distance of the waterbearing anomaly, x min is the electrode distance of the minimum apparent resistivity, and jsj only denotes the number value of the tunnel cross-sectional area s. The regression coefficients are c 1 = 1:07, c 2 = −1:54, and c 3 = 6:01. The R -square is a wide measure of fit goodness. The R-square of Equation (13) is 0.95, which shows the fit is reliable. The dotted line represents the predicted formula proposed by another researcher on numerical simulations using the water-bearing model displayed in Figure 5 [20]. They suggested a predicted formula with coefficients of c 1 = 0:25 and c 2 = c 3 = 0:00; from these coefficients, we can see that Huang et al. (2006) ignored the effect of the tunnel in the simulation so that their coefficient c 1 is much smaller than the proposed value in this study, because the coefficient c 1 is directly related to the electrode distance of corresponding minimum apparent resistivity. If c 1 is too small, then the predicted distance will be also very small. However, for dif-ferent cases of the tunnel face area, the predicted results are very different. As shown in Figure 9, when the simulated electrode distance x min is 9 m for the water-bearing anomaly at d = 6 m and for the tunnel area s = 6 m 2 , our predicted distance is d pre = 6:4 m, which is very consistent with the actual distance (d = 6 m) of the anomaly with only 6.7% of the error. However, the predicted results from Huang et al. (2006) are d pre = 2:25 m, which is much less than the actual distance (d = 6 m) of the anomaly resulting in large errors, which may lead to serious water inrush disaster and huge economic loss.

Evaluation on the Monte Carlo Method
The Monte Carlo method dealing with random problems has been widely used in geophysics [39,40]. It is very popular for time-consuming simulations on high-speed computers [41,42]. The Monte Carlo method is used for random problems in this study. Due to the connectivity of the water, the resistivity distribution of the anomaly ahead of the tunnel may be very complicated. We can use the Monte Carlo method to simulate the randomly distributed resistivity anomaly by a large number of random models. We divide the water-bearing anomaly into 1000 independent units with the dimension of 1 m × 1 m × 1 m. Each element of the anomaly has a random resistivity value less than the surrounding resistivity of 500 Ω m. In this way, we can simulate the characteristics of low resistivity water-bearing anomaly with random distribution. If the number of the simulations is large enough (e.g., the simulated number is 10,000), the probability distribution of the predicted distances for the water-bearing anomaly will be close to the actual situation, and thus the feasibility of the prediction equation in this study can be evaluated quantitatively. Based on the same pole-dipole array, different resistivity distributions of the anomaly ahead of the tunnel may lead to different predicted results. The function of the prediction distance for i th random model is defined as where d pre i stands for the i th prediction distance of waterbearing anomaly and f are the coefficients of our proposed prediction equation and the coefficients proposed by Huang et al. (2006); x i min means the i th offset (electrode distance) of the minimum apparent resistivity; ρ i 1 , ρ i 2 , ⋯, ρ i 1000 denote the i th random resistivity distribution of the anomaly; and i stands total frequency of samples varying from 1 to 10,000. Figure 10 shows the distribution of randomly generated resistivity in three different planes to simulate the characteristics of a low resistivity water-bearing anomaly ahead of the tunnel in the actual formation.

Geofluids
We calculated two groups of 10,000 random resistivity models using Equation (13) and using the equation suggested by Huang et al. (2006) by searching the offset of the minimum apparent resistivity in each FE simulation, to count all simulated prediction distances of the anomaly in the histogram. The key parameters of the Monte Carlo method are mean μ, standard deviation σ, and probability density distributions function f ðxÞ. They are determined by the following expressions: Here, n = 10, 000. The normal distribution curve is symmetrical, which can show the important information clearly. Figure 11 shows the probability of a random variable in different intervals. Figure 12 shows the results using Equation (13) and using the equation given by Huang et al. (2006) for 10,000 random anomaly models for the actual distance d = 20 m and for the tunnel face area s = 4 m 2 . The abscissa denotes the predicted distance from the anomaly to the tunnel face by two different predicted formulas, and the ordinate stands the frequency (the number of random models calculations). The predicted distances by the reference [20] range from 11 m to 16 m in Figure 12(a), and they are much less than the actual distance (d = 20 m) of the anomaly with a maximum error of 9 m, which can likely cause serious water inrush disaster if this prediction equation is used during the actual tunnel excavation. How-ever, the predicted distances in this study range from 14 m to 25 m in Figure 12(b), most of which are very close to the actual distance (d = 20 m) of the anomaly with a relatively narrow error range and accord with the normal distribution of probability perfectly (red line), which proves the predicted formula in this study is more suitable to derive the distance of the anomaly than the former (e.g. Huang et al., 2006). Figure 13 shows the detailed results of the prediction distance errors for 10,000 random low-resistivity models in this study. The majority of the prediction distances match well with the actual distance of the anomaly (d = 20 m). 56.52% of the simulated results are approximately predicted within 5% of the error for random models using the predicted Equation (13), which means that these predicted distances ranged from 19 m to 21 m. 81.55% of the results are counted within 10% of the error changing from 18 m to 22 m. 94.30% of the results are forecasted within the error of less than 15%. These results show that our predicted model is reliable and sufficiently accurate in predicting the distance of the waterrich anomaly ahead of the tunnel. Meanwhile, the random models calculated by the Monte Carlo method also make   Figure 13: Pie chart of predicted error distribution in this study.
10 Geofluids our predicted conditions closer to the actual situation in tunneling.

Conclusion
In actual forward-probing tunnel prospecting, water inrush is mostly popularly induced by the anomaly ahead of the tunnel face. Because of the water connectivity, the resistivity distribution in water-rich anomalies can be arbitrary and random. Moreover, the existence of the tunnel can distort the distribution of the potential field, which causes the error of predicted distance in the water-bearing anomaly. The tunnel effect has been discussed in this study. A predicted 3D DC-resistivity model considering the tunnel effect for forward-probing in tunnels using unstructured FEM is developed firstly in this paper. On this basis, the quality of our predicted formula is evaluated using the Monte Carlo method by simulations of 10,000 random models. Statistics show 94.30% of the results are predicted within 15% of the error, indicating that our predicted results are much better than those of the existing predicted model. The Monte Carlo method proves to be a quantitative description of the reliability for forward-probing in tunnels.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest
The authors declare no conflicts of interest.