Efficient and accurate magnetotelluric forward modeling with explicit divergence projection using second-order finite elements on tetrahedral meshes

Efficient and accurate magnetotelluric (MT) forward modeling is essential for three-dimensional MT inversion and data interpretation. While existing MT forward solvers on unstructured tetrahedral meshes predominantly employ direct solvers, the development of efficient iterative methods for high-order elements remains significantly underexplored compared to first-order elements, despite their superior accuracy. To address this, we present a unified iterative framework for both first- and second-order finite element (FE) formulations. The complex-valued FE system is transformed into an equivalent real-valued form and solved using the flexible generalized minimal residual (Flexible GMRES, FGMRES) method with an optimal block-diagonal preconditioner. This preconditioning is performed by the auxiliary-space Maxwell solver (AMS) based on the Hiptmair–Xu decomposition. To eliminate spurious divergence and accelerate convergence, particularly for models with complex topography, an explicit divergence projection (EDP) scheme is introduced to filter out non-physical components from the intermediate preconditioned vectors during FGMRES iterations. This is achieved by solving a Poisson equation for a scalar potential using the preconditioned conjugate gradient method with an algebraic multigrid preconditioner, utilizing the divergence of the conductivity-scaled vector as the source term. Numerical experiments validate our algorithm against established open-source finite difference and FE solvers, confirming that both FE approaches achieve the required accuracy, with the second-order scheme generally yielding smoother and more precise solutions. Results indicate that at lower frequencies, AMS-based preconditioning becomes sensitive to air conductivity, where excessively low values can trigger numerical instabilities that compromise solver convergence. Consequently, an appropriate choice of air conductivity is critical for sustaining numerical stability and superior solver performance. Furthermore, while complex topography typically degrades solver efficiency, the EDP technique mitigates this performance penalty by suppressing the accumulation of spurious divergence. This stabilization of iteration counts ensures robust and efficient convergence across a broad frequency range, even in the presence of challenging terrain.