## 3 Three-Dimensional Finite Element Method Forward Modeling Theory  
### Three-Dimensional Finite Element Method Algorithm  
Boreholes in roadways are typically not arranged vertically but inclined at a certain angle relative to the working face. This configuration places higher demands on the forward modeling algorithms. The finite difference (FD) and integral equation (IE), which rely on structured grids, have limitations in accurately representing this inclined borehole model. As a result, the finite-element (FE) method using unstructured tetrahedral grids [41] is employed for modeling. In this approach, the computational domain \(\Omega\) is initially divided into a finite number of tetrahedral meshes, as illustrated in Figure 2. The forward modeling of the cross-hole grounded-wire-source TEM is based on full-space theory, which simplifies the background model as it does not need to consider the influence of air.  
Figure 1: Schematic diagram of cross-hole grounded-wire-source transient electromagnetic method.  
Figure 2: Grid subdivision for full-space EM model.  
Starting from [PERSON]'s equations and using the Galerkin weighted residual method [42], the residual \(R\) of the electric field diffusion equation can be expressed as  
\[R^{e}=\] \[\iint_{\Omega^{\prime}}\mathbf{w}\cdot\left(\frac{1}{\mu}\
abla \times\
abla\times\mathbf{E}(t)+\sigma\frac{\partial\mathbf{E}(t)}{\partial t }+\frac{\partial\mathbf{j}(t)}{\partial t}\right)dv\, \tag{1}\]  
where \(\mu\) and \(\sigma\) are the magnetic and electric permeability, respectively, \(E_{(t)}\) and \(j_{(t)}\) are the electric field and current density at the time \(t\), \(\Omega^{e}\) is the volume of the \(e\)-th element, and \(\mathbf{w}\) is the weight function. Furthermore, the interpolated basis function \(\mathbf{N}\) is the same as the weight function; Equation (1) can then be rewritten as  
\[\begin{split}& R_{i}^{e}=\frac{1}{\mu}\sum_{j=1}^{6}E_{i}^{e}(t) \iint_{\Omega^{\prime}}\
abla\times\mathbf{N}_{i}\cdot\
abla\times\mathbf{N}_{ j}dv\\ &+\sum_{j=1}^{6}\frac{\partial E_{j}^{e}(t)}{\partial t}\iint_{ \Omega^{\prime}}\mathbf{N}_{i}\sigma\mathbf{N}_{j}dv+\iint_{\Omega^{\prime}} N_{i}\cdot\frac{\partial\mathbf{j}(t)}{\partial t}dv\end{split}\, \tag{2}\]  
The global system is constructed by assembling the elements in the computational domain with the weighted residuals \(R=0\), so that Equation (2) can be expressed as a linear equation:  
\[\mathbf{M}\frac{d\mathbf{E}(t)}{dt}+\mathbf{S}\mathbf{E}(t)+\mathbf{J}=0, \tag{3}\]  
where \(\mathbf{M}\) is the mass matrix, \(\mathbf{S}\) is the stiffness matrix, and \(\mathbf{J}\) is the source term matrix. The time derivative in Equation (3) is discretized by second-order backward Euler equation to obtain the following form:  
\[(3\mathbf{M}+2\Delta t\mathbf{S})\mathbf{E}^{i+2}=\mathbf{M}\Big{(}4\mathbf{ E}^{i+1}-\mathbf{E}^{i}\Big{)}-2\Delta t\mathbf{J}^{i+2}, \tag{4}\]  
where \(\mathbf{E}^{i+2}\), \(\mathbf{E}^{i+1}\), and \(\mathbf{E}^{i}\) denote the electric field at the \((i+2)\)-th, \((i+1)\)-th, and \(i\)-th moments, and \(\mathbf{J}^{i+2}\) is the source term at the \((i+2)\)-th moment. \(\Delta t\) denotes the time step, which is taken as a step growth in the calculations in this paper [43]. In this way, only Equation (4) needs to be solved to obtain the values of the electric field at all edges of the computational region and, thus, we can obtain the electromagnetic field parameters at the locations of the receivers. This leads to the following large linear equation:  
\[\mathbf{K}\mathbf{E}=\mathbf{b}, \tag{5}\]  
The matrix form can be expressed as  
\[\begin{pmatrix}3M+2\Delta t_{1}S\\ -4M&3M+2\Delta t_{2}S\\ M&-4M&3M+2\Delta t_{3}S\\ M&-4M&3M+2\Delta t_{4}S\\ &&\ddots&\ddots&\ddots\\localized mesh refinement conducted at the receivers, sources, and boundary surfaces of the anomalous bodies.