## 2 Forward Modeling Principle  
### Governing Equation  
Ignoring the effects of displacement current and magnetic permeability on the time-domain electromagnetic field, the Maxwell equations can be written in the following form [20; 27; 28]:  
\[\
abla\times\mathbf{E}(\mathbf{r},t)=-\frac{\partial\mathbf{B}(\mathbf{r},t)}{\partial t}, \tag{1a}\] \[\
abla\times\mu^{-1}\mathbf{B}(\mathbf{r},t)-\sigma(r)\mathbf{E}(\mathbf{r},t)=\mathbf{s}(\mathbf{r},t). \tag{1b}\]  
In these equations, \(\mathbf{E}(\mathbf{r},t)\) and \(\mathbf{B}(\mathbf{r},t)\) represent the electric field strength \(\left(\mathrm{V}\!\cdot\!\mathrm{m}^{-1}\right)\) and magnetic flux density \(\left(\mathrm{T}\right)\) at time \(t\) within the study area, respectively. \(\mu\) is the magnetic permeability of free space \(\left(\mathrm{H}\!\cdot\!\mathrm{m}^{-1}\right)\), \(\sigma(\mathbf{r})\) is the isotropic (or anisotropic) conductivity of the medium in the calculation domain \(\left(\mathrm{S}\!\cdot\!\mathrm{m}^{-1}\right)\), and \(\mathbf{s}(\mathbf{r},t)\) is the externally applied current density\(\left(\mathrm{A}\!\cdot\!\mathrm{m}^{-2}\right)\). Taking the curl of both sides of Equation (1a) and combining it with Equation (1b), we obtain the equation for the curl of the electric field [41; 42]:  
\[\
abla\times\
abla\times\mathbf{E}(\mathbf{r},t)+\mu\sigma(\mathbf{r})\frac{\partial\mathbf{E }(\mathbf{r},t)}{\partial t}+\mu\sigma(\mathbf{r})\frac{\partial\mathbf{s}(\mathbf{r},t)}{ \partial t}=0. \tag{2}\]  
This paper uses a tetrahedral mesh to discretize the study area and employs the vector finite element method to discretize Equation (2), yielding the governing equation for the electric field strength:  
\[\int\
u\ \left(\
abla\times\mathbf{N}_{i}\right)\cdot\left(\
abla\times\mathbf{N}_{j} \right)\mathrm{d}V\mathbf{E}+\mu\int\
u\ \sigma(\mathbf{r})\mathbf{N}_{i}\cdot\mathbf{N}_{j}\mathrm{d}V\frac{\partial\mathbf{E}(\mathbf{r},t )}{\partial t}=\mu\int\
u\ \ N_{i}\cdot N_{j}\mathrm{d}V\frac{\partial\mathbf{s}(\mathbf{r},t)}{ \partial t}, \tag{3}\]  
where \(N_{i}\) and \(N_{j}\) represent the interpolation basis function. Assume \(\mathbf{M}_{\mathrm{curl}}\) and \(\mathbf{M}_{\sigma}\) as follows [43]:  
\[\mathbf{M}_{\mathrm{curl}}=\int\
u\ \left(\
abla\times N_{i}\right)\cdot \left(\
abla\times N_{j}\right)dV \tag{4}\] \[\mathbf{M}_{\sigma}=\int\
u\ \mu_{0}\sigma N_{i}\cdot N_{j}dV\]  
In this paper, the electric field is defined on the edges of the tetrahedral mesh, and the first-kind boundary condition is applied, assuming that the electromagnetic field decays to zero at sufficiently distant boundaries. A step-off source is used as the excitation waveform. When calculating the secondary field after the current is turned off, the current term \(\mathbf{s}(\mathbf{r},t)\) is set to 0. Thus, Equation (3) can be written as follows [43]:  
\[\frac{\partial\mathbf{E}(\mathbf{r},t)}{\partial t}=-\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{ \mathbf{\mathbf{\mathbf{\  
The ODE, which depends solely on the time parameter \(t\), is represented by Equation (5). The approximation formula for solving \(\mathbf{E}(t)\) can be applied [31; 43; 44].  
\[\mathbf{E}(t)\approx\mathbf{E}_{m}(t)=\mathbf{V}_{m}e^{-\frac{t}{\gamma}(\mathbf{H}_{m,m}^{- 1}-\mathbf{I})}\parallel\mathbf{E}_{0}\parallel_{2}\mathbf{e}_{1}, \tag{6}\]  
where \(\gamma\) is the optimal shift and \(m\) is the minimal subspace order. \(\mathbf{e}_{1}=[1,0,\dots,0]_{m}^{T}\), and \(\mathbf{V}_{m}\), and \(\mathbf{H}_{m,m}\) are the orthonormal basis with \(m\) columns and the Hessenberg matrix, respectively. The Arnoldi technique is used in this paper to compute \(\mathbf{V}_{m}\) and \(\mathbf{H}_{m,m}\) in the above formula [45]. Additionally, \(e^{-\frac{t}{\gamma}(\mathbf{H}_{m,m}^{-1}-\mathbf{I})}\) has a small-dimensional order, which can be efficiently solved using the freely available code [46].