Fast 3D simulation of magnetotelluric data in anisotropic media using a rational Krylov method

This study introduces a computational method aimed at accelerating multifrequency 3D magnetotelluric (MT) forward modeling in 3D conductivity structures with general anisotropy. The approach combines an edge-based finite-element (FE) method with a rational Krylov subspace method. In this method, the MT source term is expressed as a planar current source, and the frequency-dependent electric field response is the product of a transfer function and a constant vector of the current source. Consequently, the rapid approximate electric field for multifrequency calculations can be obtained by constructing the orthogonal basis of the rational Krylov subspace. By applying this technique, the large sparse matrix resulting from FE discretization can be projected onto a much lower-order matrix, significantly reducing its size to tens or hundreds. Incorporating the direct solver PARDISO to construct the orthogonal basis enables a considerably faster solution compared with the previous FE method. Moreover, the algorithm is implemented using Julia, a high-level programming language known for its readability, maintainability, and extensibility. The validity of the algorithm is demonstrated through a 1D anisotropic model, and two additional 3D models are designed to further evaluate the applicability and effectiveness of our algorithm. Numerical experiments illustrate a notable increase in modeling speed compared with the traditional edge-based FE method.