Abstract

The transient electromagnetic (TEM) method has long been applied in tunnel advanced prediction. However, it remains questionable to what extent a geologic anomaly body will influence the induced electromagnetic response in front of the heading face. The dilemma is partly because observed TEM data are frequently interpreted by empirical formulas or proportional relationships, and a quantitative measurement has not been established. In this paper, we strive to understand the TEM characteristics from a 3D finite-element time-domain (FETD) modeling aspect. The modeling algorithm is based on unstructured space meshing and unconditional stable time discretization, which ensures its accuracy and stability. The modeling algorithm is verified by a half-space model, in which the misfit of late-time channels that we are concerned with is generally below 1%. The algorithm has also been utilized to carry out the TEM response of tunnel models with different types of TEM devices. Through model studies, we find that both the traditional central-loop device and the recently developed weak-coupling opposing-coil device are feasible in tunnel advanced detection. Nevertheless, the latter type of device better distinguishes low-resistivity anomalies at 30 m ahead of the heading face with a relative difference (between models with and without the anomaly) of more than 1000% at certain time channels, compared with only a 10% difference of the central-loop device. Also, we conclude that the vertical electromagnetic field component should be recorded and interpreted together with the horizontal field to provide more convincing results.

1. Introduction

The tunnel advanced prediction method is aimed at determining blind faults, karsts, and other geologic bodies ahead of the working face [1]. Several types of geophysical prospecting methods, including seismic tomography [2], reverse time migration [3], direct current electric sounding [4], magnetic resonance sounding [5], and the transient electromagnetic (TEM, also called time-domain electromagnetic, TDEM) method [6], have been applied to tunnel advanced prediction. Among them, electromagnetic methods have the unique advantage that they are more sensitive to the composition of heading rocks rather than the distribution of geologic interfaces [7]. More explicitly, electromagnetic surveys can help find water-bearing faults and karsts with relatively high electric conductivity [8]. Thus, they play a crucial role in preventing water inrush disasters and protecting lives and property during the tunnel excavation process.

As one of the electromagnetic methods, the TEM method was first proposed by Keller in 1969 [9]. Since then, the TEM method has found an increasingly wide utilization, including the long-offset TEM method with dipole source [10], the land big-loop TEM method with square coil source [11], and the airborne TEM method with circular coil source [12]. The TEM method has shown its advantage at sites where space is highly limited since the transmission of pulse-type electromagnetic signals does not require a high-power device. The above aspects prompt applying the TEM method in tunnel advanced detection, and it has thrived over the past two decades [1316]. However, most of these interpretation results were carried out through empirical relationships between the buried depth of the geologic body and the arrival time of the electromagnetic wave. And thus, the geophysical society still holds the question of whether the response at a certain time point is caused by a specific anomaly body or, comparatively, to what extent the anomaly body will influence the recorded field.

3D modeling provides explicit theoretical support for geoprospecting methods on a physical basis. In the time domain, electromagnetic modeling is often carried out by the finite-difference time-domain (FDTD) method and the finite-element time-domain (FETD) method. With the FDTD method, the calculation domain is discretized with structured meshes, and a conditionally stable time-stepping technique is required [17]. In comparison, the FETD method allows unstructured meshing that can accurately fit the terrain and irregular geologic bodies [18]. And by introducing unconditionally stable stepping methods, the FETD method allows larger step size in late-time simulation, and thus, much fewer time steps are required [19]. These advantages of the FETD method ensure both its computational accuracy and efficiency and make it suitable for geoelectromagnetic modeling.

In the specific application of tunnel advanced prediction, an appropriate design of the TEM device is also crucial. For instance, it has been reported that TEM data at early time channels are sophisticated with primary-field disturbances, which may further result in a blind zone at shallow depths [20]. Meanwhile, the mutual inductance problem arises with the traditional central-loop device, as multiturn wires are twined tightly with each other [21]. To address these issues, weak-coupling TEM devices have recently been developed [22, 23]. And in this paper, TEM responses of both types of device are calculated to further confirm the feasibility of tunnel advanced prediction using the TEM method. In the following sections, the proposed 3D FETD algorithm is first illustrated, and its accuracy is verified by a half-space analytical solution. The feasibility study is then conducted through various model simulations.

2. Methods

With Faraday’s law and Ampère’s law we acquire the vector electromagnetic wave equation as

In these formulas, represents the external current source, in our case, the TEM source of the system. , , and are the magnetic permeability, the dielectric permittivity, and the electric conductivity of the medium, respectively. A rational approximation in low-frequency geoelectromagnetic simulation is to eliminate the second derivative term in Equation (1) [24], which gives

In the above equation, the boundary condition of a perfect electrical conductor on the domain boundary and the initial condition give the boundary-initial value problem to solve.

2.1. Discretization in Space and Time

Following the classical finite element method [25], we have the system of equations as

Here, represents the global electric field vector under the vector finite element basis [26]. With this set of basis functions, the electric field to be measured can be approximated by where is the total edge number of the discretized domain. And global integral matrices , , and have the form of where the superscript represents a local element and represents the whole simulation domain.

To discretize the ODE system (Equation (7)) in the time domain, we apply the unconditionally stable Newmark- method with [27] and obtain a time-stepping equation where represents the field vector at the step. The stepping matrices , and have the form of where is the current step size.

In TEM modeling, specifically, Um et al. [28] pointed out that fixed step size cannot cover the wide time magnitudes from early-time transient response to late-time diffusive field, and thus a step-changing scheme should be imposed. Based on this, the initial step is defined as in which represents the pulse width of the TEM source, and we double every 25 steps in this research to ensure stability. It is also worth noticing that within these steps at which does not change, the sparse matrices , , and also keep unchanged. Therefore, the symbolic factorization [29] of , which is the most time-consuming part of 3D TEM modeling, is only performed a limited (typically less than 20) times throughout the simulation. Finally, the electric field is calculated by Equation (8), and the induced electromotive force in a certain direction is derived by

If needed, the total magnetic field at a specific time point can then be deduced by the time integral.

2.2. Numerical Implementation

The construction of the source term (Equation (11)) is challenging. Since if the wire passes through an element arbitrarily, a complicated integral has to be carried out. One simplification is to directly subdivide the transmitting wire, by which the element integral only has to be calculated on the edges. With this consideration, Equation (11) could be rewritten as in which the subscript denotes the global edges that the source wire lays on. is the direction of the current, and values 1 when and are in the same direction and values -1 when they are in the opposite direction [30]. In the above equation, the relationship holds, and represents the electric moment of the source. As a typical pulse type used in TEM signal transmission, the step pulse is simulated in this research. We approximately formulate the unit step pulse with [31]. where represents the error function.

In this paper, 3D unstructured tetrahedral meshing is applied, and the open-source software Gmsh is used to generate the consequent Delaunay triangulation [32, 33]. With this said, the two types of TEM devices of our concern are illustrated in Figure 1. It can be seen from Figure 1(a) that always has the same direction within a single coil, yet points from the node with a lower identifier to the node with a higher identifier according to a common definition [30]. Thus, equals -1 on the edge that links points 1 and 26, and it equals 1 elsewhere with regard to this example. Hereinafter, the traditional small-loop TEM device with one transmitting coil shown in Figure 1(a) is referred to as the single-TX device, and the device with two opposing coils shown in Figure 1(b) is referred to as the opposing-TX device.

In a typical 3D geoelectromagnetic modeling problem, the total edge number can reach up to several million. And specific to the time-domain variation, an iteration over several hundreds of time steps has to be calculated. Therefore, the efficient solving of the huge sparse system occupies a pivotal place in the simulation. A sparse system is traditionally solved by iterative solvers, and various types of algorithms have been developed [34]. However, the iterative method has several drawbacks when solving a system with an extremely ill-conditioned left-hand side matrix or multiple right-hand side vectors. With a direct solver, comparatively, the sparse matrix is first handled by numerical factorization, and solutions of different right-hand sides are then derived by forward and backward substitution [35]. This separate design of the direct solver services our needs as stated in the above-proposed algorithm. And PARDISO [36], a parallel realization of direct sparse solver, is applied in this research. The performance of the sparse linear algebra library is also crucial to the overall computational time of 3D TEM modeling. In our program, the construction of sparse system matrices, the addition of sparse matrices, and the sparse matrix-vector multiplication are realized by a high-performance sparse BLAS library [37]. The following simulations are carried out with an 8-core Intel I7-11700 CPU.

3. Verification of Algorithm Accuracy

To test the reliability of the proposed algorithm, the TEM response of a half-space model is first calculated. In this model, a single-TX device of radius is placed at the origin with its axis in the -direction, and the vertical induced electromotive force is calculated at the origin. The transmitting current is 1A. The model is discretized into 162719 tetrahedrons with 190779 edges in total, and it consists of an above air half-space and a below , half-space. Then, the FETD solution is verified by the analytical solution [38]. where represents the error function and parameter .

We notice from Figure 2 that the FETD solution is well coincident with the analytical solution after . However, the numerical solution slightly deviates from the ideal case at early time channels. This is because both the primary transmitting field and the secondary inductive field contribute to the FETD solution, while only the secondary field contributes to the analytical solution. This is the FETD simulation of pulse width in this verification, which is more of an indication of the real world rather than an ideal situation. With regard to computational speed, the total simulation ends within 120 seconds, in which only 11 numerical factorizations and 254 backward substitutions are required.

4. Synthetic Tunnel Model Studies

In this section, the characteristics of tunnel TEM response are studied by synthetic model studies. In this model, a brick abnormal body of 20 m wide and 5 m thick, is placed within the background. As shown in Figure 3, the dimension of the tunnel is , and it extends from the leftmost of the model to . Both the abnormal body and the tunnel are symmetric with respect to the -axis.

4.1. Different Positions of Anomalies

One of the major concerns about the feasibility of the tunnel TEM prediction method is how well it could identify abnormal low-resistivity bodies [3941]. Here, the single-TX transmitting device is first studied. The coil has a radius of 0.5 m, and it is located at with its axis in the -direction. The transmitting current of 1A keeps unchanged. Figure 4(a) gives the responses at the receiving point (−1 m, 0 m, and 0 m) with the low-resistivity body located at , 20 m, and 30 m, respectively. And percentage relative differences of between the three models with anomaly and the model without anomaly are given in Figure 4(b). It could be noticed that with the increase of buried depth, the relative difference between model responses with and without the low-resistivity body drops dramatically. Also, a small enough must be chosen (in this example ) to avoid the interference of the primary TEM field and to ensure the max relative difference occurs after the transmitting current turns off.

For a deeper understanding of the TEM response within a whole-space tunnel environment, the spatial distribution of the field at s is given in Figure 5. We can observe that the transmitted electromagnetic field concentrates if a low-resistivity body exists. Contrarily, the field flattens out over the entire domain. This further suggests that the anomaly recorded by our TEM device is merely a repercussion of field enhancement in the distance. And if the data in some predrilling holes are collected, and more precise advanced detection results will be given.

4.2. The Presence of the Tunnel

Another concern about the feasibility of tunnel TEM detection is the presence of the empty tunnel. To evaluate this, model responses with and without the tunnel are compared in Figure 6, and the relative difference of is calculated similarly to that of Figure 4(b). In the model without a tunnel, its air conductivity of the tunnel is replaced by the background conductivity . And in both models, the same low-resistivity body is located at . Figure 6 shows that the existence of the tunnel has a negligible effect on the simulation result. This is because the tunnel can be regarded as an anomaly of high electric resistivity, which typically has minimal contribution to the total electromagnetic field. Therefore, we conclude that the presence of an empty tunnel does not impact the feasibility of tunnel TEM advanced detection.

4.3. Different Field Components

It is also noteworthy to understand the characteristics of the vertical response. Using the same single-TX model, the field is again calculated as shown in Figures 7 and 8. The most obvious feature is that the response curve goes through several sign reversals (Figure 7(a)) because of the nature of small-loop TEM devices. We also notice that the response patterns of the two models with the anomaly at and 30 m are quite similar to that of the model without the low-resistivity anomaly. This can be further explained in Figure 8, as the field recorded by the single-TX device can be significantly affected by the anomaly at , while it is hard to be affected in the other two cases.

4.4. Different Types of Transmitting Coils

Here, the feasibility of the novel opposing-TX device is further studied. To achieve this, the same tunnel TEM model as the last simulation is utilized, with the only difference being that we replace the 0.5 m-radius coil with two coils carrying opposite directions of transmitting current. The two coils are spaced 0.2 m apart. And the center of the device, which is also the receiving point, keeps unchanged at (−1 m, 0 m, and 0 m). The response of the same abovementioned models is calculated as shown in Figures 9 and 10. We conclude from the comparison between Figures 4 and 9 that the opposing-TX device behaves with a much larger percentage difference between models with and without the low-resistivity body over the entire observation period. However, its field amplitude is lower than that of the traditional single-TX device. From this forward modeling point of view, the opposing-TX device, which records pure secondary TEM field, is also suitable for tunnel advanced prediction.

5. Conclusions

This paper applies the FETD method to simulate the TEM response in a tunnel environment, offering a theoretical assessment of the effectiveness of the TEM method for advanced prediction in tunnels. Through tunnel model studies, we find that both the traditional single-TX device and the opposing-TX device are capable of revealing low-resistivity bodies. Nevertheless, the latter device performs better in the relative difference of the axial field at the cost of a lower field amplitude. The presence of the high-resistivity tunnel has little influence on TEM advanced prediction. Meanwhile, it is also found that the horizontal field (the field in this paper) should be interpreted along with the vertical field (the field in this paper). Based on these 3D simulations, we conclude that the TEM method is crucial and feasible for tunnel advanced prediction, and that weak-coupling devices are promising for future fieldwork.

Data Availability

We confirm that no additional data is available.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.

Acknowledgments

This research is funded by the National Natural Science Foundation of China (No. 42304080), the Joint Open Fund of Mengcheng National Geophysical Observatory (No. MENGO-202301), the Anhui Provincial Natural Science Foundation (No. 2308085QD126), and the Fundamental Research Funds for the Central Universities, China (No. WK2080000173).