Feasibility of Tunnel TEM Advanced Prediction: A 3D Forward Modeling Study

,


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 highpower device.The above aspects prompt applying the TEM method in tunnel advanced detection, and it has thrived over the past two decades [13][14][15][16].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.

Methods
With Faraday's law and Ampère's law we acquire the vector electromagnetic wave equation as In these formulas, J 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, E 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 n edge is the total edge number of the discretized domain.And global integral matrices A, C, and S have the form of • N e j dΩ, 9 where the superscript e represents a local element and Ω represents the whole simulation domain.

Geofluids
To discretize the ODE system (Equation ( 7)) in the time domain, we apply the unconditionally stable Newmark-β method with β = 1/4 [27] and obtain a time-stepping equation where E n represents the field vector at the n th step.The stepping matrices P, Q and R have the form of where Δt 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 Δt 0 = t p /20 in which t p represents the pulse width of the TEM source, and we double Δt every 25 steps in this research to ensure stability.It is also worth noticing that within these steps at which Δt does not change, the sparse matrices P, Q, and R also keep unchanged.Therefore, the symbolic factorization [29] of P, 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 n is derived by If needed, the total magnetic field at a specific time point t c 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 I denotes the global edges that the source wire lays on.n is the direction of the current, and N I • nI values 1 when N I and nI are in the same direction and values -1 when they are in the opposite direction [30].
In the above equation, the relationship holds, and p 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 erf 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 n always has the same direction within a single coil, yet N points from the node with a lower identifier to the node with a higher identifier according to a common definition [30].Thus, N • n 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.

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 R = 0 5 m is placed at the origin with its axis in the z-direction, and the vertical induced electromotive force ε is calculated at the origin.The transmitting current I 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 σ = 10 −2 S m −1 , μ = 4π × 10 −7 Hm −1 half-space.Then, the FETD solution is verified by the analytical solution [38].where erf represents the error function and parameter θ = μσ/4t.We notice from Figure 2 that the FETD solution is well coincident with the analytical solution after 2 × 10 −5 s.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 t p = 5 × 10 −7 s 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.

Synthetic Tunnel Model Studies
In this section, the characteristics of tunnel TEM response are studied by synthetic model studies.In this model, a σ = 1Sm −1 brick abnormal body of 20 m wide and 5 m thick, is placed within the σ = 10 −2 Sm −1 background.As shown in Figure 3, the dimension of the tunnel is 2 m × 2 m, and it extends from the leftmost of the model to x = 0 m.Both the abnormal body and the tunnel are symmetric with respect to the x-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 lowresistivity bodies [39][40][41].Here, the single-TX transmitting 5 Geofluids device is first studied.The coil has a radius of 0.5 m, and it is located at x = −1 m with its axis in the x-direction.The transmitting current of 1A keeps unchanged.Figure 4(a) gives the ε x responses at the receiving point (−1 m, 0 m, and 0 m) with the low-resistivity body located at x = 10 m, 20 m, and 30 m, respectively.And percentage relative differences of ε x 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 t p must be chosen (in this example 10 −6 s) 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 ε x field at 2 5 × 10 −5 s is given in Figure 5.We can observe that the transmitted electromagnetic field concentrates if a low-resistivity body exists.Contrarily, the ε x field flattens out over the entire domain.This further suggests that the ε x anomaly recorded by our TEM device is merely a repercussion of field enhancement in the distance.And if the ε x 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 ε x is calculated similarly to that of Figure 4(b).In the model without a tunnel, its air conductivity 10 −8 Sm −1 of the tunnel is replaced by the background conductivity 1 0 −2 Sm −1 .And in both models, the same low-resistivity body is located at x = 20 m. 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 con-tribution 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.

Different Field Components.
It is also noteworthy to understand the characteristics of the vertical ε z response.Using the same single-TX model, the ε z 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 x = 20 m 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  7 Geofluids ε z field recorded by the single-TX device can be significantly affected by the anomaly at x = 10 m, while it is hard to be affected in the other two cases.

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 ε x 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.

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 lowresistivity bodies.Nevertheless, the latter device performs  8 Geofluids better in the relative difference of the axial ε x field at the cost of a lower field amplitude.The presence of the highresistivity tunnel has little influence on TEM advanced prediction.Meanwhile, it is also found that the horizontal field (the ε x field in this paper) should be interpreted along with the vertical field (the ε z 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 weakcoupling devices are promising for future fieldwork.

Figure 1 :Figure 2 :
Figure 1: Sketch of (a) traditional small-loop TEM device with one transmitting coil and (b) small-loop TEM device with two opposing coils.Red arrows represent the direction N of global edges, green arrows represent the direction n of transmitting current, and numbers are global node identifiers.

Figure 3 :Figure 4 :
Figure 3: (a) Sketch of the synthetic tunnel TEM model and (b) cross-section of its tetrahedral meshing result.

Figure 5 :Figure 6 :
Figure 5: Spatial ε x response (a) with a low-resistivity body located at x = 20 m and (b) without any low-resistivity body in the xOz plane of the single-TX device at 2 5 × 10 −5 s.

Figure 7 :Figure 8 :
Figure 7: (a) ε z response of the single-TX device.(b) Percentage relative difference of the ε z response between models with and without the low-resistivity anomaly at x = 10 m, 20 m.

Figure 9 :Figure 10 :
Figure 9: TEM response of the opposing-TX device.(a) ε x response of the synthetic tunnel model without and with the low-resistivity anomaly at different positions.(b) Percentage relative difference of the ε x response between models with and without the low-resistivity anomaly.