## 2 Forward Modeling Principle
### Large-Scale SAI-Krylov Algorithm
In Figure 10b, as the offset distance increases and the geological anomaly's influence grows, the early induced electromotive force (\(z\)-component) does not decay significantly in the double-logarithmic coordinate system. In fact, the electromagnetic signal slightly increases before it gradually decays. Additionally, due to the impact of the low-resistivity anomaly, the electromagnetic field at the receiver point (0 m, 1000 m, 0 m) undergoes a reversal in sign. The three-dimensional vertical contact band model further validates the effectiveness of the proposed method. Figure 10: The numerical results of 3D vertical contact zone: (**a**) loop source device [16; 28]; (**b**) grounded source device [14; 47]. Figure 9: Three-dimensional complex model with a vertical contact zone: (**a**) loop source device; (**b**) grounded source device. ### Large-Scale Undulating Stratigraphy Model

This paper constructs a two-dimensional large-scale model with five layers of complex interfaces. Figure 11 shows the resistivity distribution of the formations across the complex interface, which spans from \(-500\) m to \(500\) m along the \(x\)-axis. A long wire source is used to transmit and receive the TEM response \(\mathbf{E}_{x}\) and \(\partial\mathbf{B}_{z}/\partial t\). The long wire source has a length of 2 km, with its coordinates at (\(-1000\) m, \(100\) m, \(0\) m) and (\(1000\) m, \(100\) m, \(0\) m), and an emission current of 1 A. The receiving survey line, with a length of 1 km, is positioned along the \(x\)-axis, with endpoint coordinates at (\(-500\) m, \(0\) m, \(0\) m) and (\(500\) m, \(0\) m, \(0\) m). To achieve high-precision forward responses and properly model topographic fluctuations, the minimum mesh length used in this study is 5 m, and the grid division diagram is shown in Figure 11. The minimum grid size is 10 m, with the total number of unknowns reaching 7.4 million and a computation time of 14,984 s. On the other hand, the memory usage of this method is only 6.3 GB. Figure 12 shows cross-sectional views of the electric field \(\mathbf{E}_{x}\) at 0.01 ms, 0.25 ms, 1 ms, and 10 ms. At \(t=0.01\) ms, the electric field \(\mathbf{E}_{x}\) is mostly concentrated near the surface, and the interface topography of the first and second layers can be vaguely observed in Figure 12a. At \(t=0.25\) ms, the electric field \(\mathbf{E}_{x}\) gradually diffuses downwards, and the formation geometry of the second and third layers becomes clearer, although the deep-layer formation information is not fully reflected. As time progresses, at \(t=1\) ms, the geometry of the undulating formation aligns well with the contour distribution of the electric field \(\mathbf{E}_{x}\). The distribution of \(\mathbf{E}_{x}\) at \(t=10\) ms is similar to that at \(t=1\) ms, but the extreme value of the electric field \(\mathbf{E}_{x}\) decreases and diffuses deeper into the formation. Figure 13 shows cross-sectional views of the induced electromotive force \(\partial\mathbf{B}_{z}/\partial t\) at 0.01 ms, 0.25 ms, 1 ms, and 10 ms, respectively. Figure 13a is similar to Figure 12a, as both effectively reflect the shallow horizontal stratigraphic interfaces and show a slight response to the first undulating interface. However, the electromagnetic anomaly in Figure 13a is weaker than that in Figure 12a. Figure 13b,c capture the undulating strata and the underlying stratigraphic uplift, but compared to Figure 12b,c, the electromagnetic field anomalies show lower consistency with the undulating stratigraphic morphology. When the diffusion time reaches \(t=10\) ms, the electromagnetic anomaly in the shallow strata of Figure 13d disappears, making it difficult to accurately capture the structural morphology of the deeper strata. Figure 11: Two-dimensional undulating stratigraphy model: (**a**) model schematic diagram; (**b**) grid division diagram.