## 2 Forward Modeling Principle
### Large-Scale SAI-Krylov Algorithm  
As seen in step 3 of Algorithm 1, linear equations must be solved when constructing the matrices \(\mathbf{V}_{m}\) and \(\mathbf{H}_{m,m}\). The direct method becomes challenging to apply when the simulated geological model is large-scale or consists of a high number of meshes. Therefore, an iterative strategy is employed in this paper to solve step 3 [17]. ```
Input:\(\mathbf{C}\in\mathbb{R}^{n\times n}\), \(\mathbf{E}_{0}\in\mathbb{R}^{n}\), \(m\)
```  
1. set \(\mathbf{v}_{1}:=\mathbf{E}_{0}/\|\mathbf{E}_{0}\|_{2}\)
2. Do \(j=1\), \(2\),..., \(m\)
3. compute \(\mathbf{w}_{j}:=(\mathbf{I}+\gamma\mathbf{C})\mathbf{v}_{j}\)
4. Do \(i=1\), \(2\),..., \(j\)
5. \(h_{i,j}:=\left(\mathbf{w}_{j},\mathbf{v}_{i}\right)\)
6. \(\mathbf{w}_{j}:=\mathbf{w}_{j}-h_{i,j}\mathbf{v}_{i}\)
7. End do
8. \(h_{i+1,j}:=\|\mathbf{w}_{j}\|_{2}\). If\(h_{j+1,j}=0\) Stop
9. \(\mathbf{v}_{j+1}:=\mathbf{w}_{j}/h_{j+1,j}\)
10. End do
11. \(\mathbf{V}_{m}:=\left[\mathbf{v}_{1},\mathbf{v}_{2},\dots,\mathbf{v}_{m}\right]\), \(\mathbf{H}_{m+1,m}:=\left[\mathbf{h}_{1},\mathbf{h}_{2},\dots,\mathbf{h}_{m}\right]\) ```
Output:\(\mathbf{V}_{m}\), \(\mathbf{H}_{m,m}\)
```

**Algorithm 1** Constructing reduced-order basis and projection matrix

By rewriting step 3 of Algorithm 1, we obtain the following:

\[(\mathbf{M}_{\sigma}+\gamma\mathbf{M}_{\text{curl}})\mathbf{\omega}_{j}=\mathbf{ M}_{\sigma}\mathbf{v}_{j}, \tag{7}\]

Here, the PCG technique (Algorithm 2) is employed to solve Equation (7) in this study. Rewriting Equation (7) as \(\mathbf{A}\mathbf{x}=\mathbf{b}\), where \(\mathbf{A}=\mathbf{M}_{\sigma}+\gamma\mathbf{M}_{\text{curl}}\) and \(\mathbf{b}=\mathbf{M}_{\sigma}\mathbf{v}_{j}\), the expression for the preprocessing matrix \(\mathbf{M}\) is as follows:

\[\mathbf{M}=(2-\omega)^{-1}(\mathbf{D}/\omega+\mathbf{L})(\mathbf{D}/\omega)^{ -1}(\mathbf{D}/\omega+\mathbf{L})^{T}, \tag{8}\]

Here, \(\mathbf{D}\) represents the diagonal matrix, \(\mathbf{L}\) denotes the lower triangular matrix, and \(\omega\) is the relaxation factor, where \(0<\omega<2\). \[ \begin{split} W=\mathbf{D}/\omega+\mathbf{L},& \mathbf{V}=(2-\omega)\mathbf{D}/\omega\\ y=\mathbf{W}^{-1}g,z=\mathbf{W}^{-T}d\end{split} . \tag{9}\]The modeling simulations in this paper were all conducted on a workstation equipped with 128 GB of RAM and a 2.4 GHz processor. ### The Selection of the Initial Value and Tolerance for the PCG

This section describes the process for selecting the initial values for the iterative algorithm. The time range considered in this study is \(\big{[}10^{-5},10^{-2}\big{]}\)s. The conductivity of air and earth is \(\sigma_{air}=10^{-4}\) S/m and \(\sigma_{earth}=10^{-2}\) S/m, respectively. The excitation device used in this study employs a grounded source, with a length of 100 m. According to [PERSON] et al., the optimal shift for this study is \(\gamma=2.6\times 10^{-7}\), and the Krylov subspace order is \(m=320\)[26]. To ensure numerical accuracy, the tolerance \(\epsilon\) is chosen to be sufficiently small, with \(\epsilon=10^{-21}\). Figure 1a shows the relative error curves for different initial values, with the relative error varying significantly depending on the initial value. As seen in Figure 1a, the relative error decreases as the initial value decreases and eventually stabilizes. When the initial value is 100, the relative error changes dramatically, surpassing 100% around \(t=1\) ms. The relative error curves for initial values of \(10^{-1}\) and \(10^{-2}\) show similar characteristics to the curve with an initial value of \(10^{0}\). However, the relative error trend becomes much gentler for initial values of \(10^{-3}\) and \(10^{-4}\), though it still does not meet the precision criteria for TEM forward simulations. When the initial value is further reduced, the relative error continues to decrease, and the curves become more stable. Figure 1b shows the relative error curves for different tolerance values \(\epsilon\).