The marine magnetotelluric (MMT) method is a crucial tool in electrical exploration, capable of inferring the subsurface electrical structure of the Earth by measuring naturally varying electromagnetic (EM) signals at the seafloor or ocean surface [1,2]. This technique has been widely employed in offshore research, including studies of the Earth’s tectonics hydrocarbon exploration, and so forth [3,4,5].

Conductivity and magnetic permeability are both essential physical properties of the Earth. Previous studies have shown that the susceptibility of certain geological materials significantly exceeds that of free space [6]. Common minerals exhibit a wide range of magnetic susceptibility variations, with magnetite minerals displaying exceptionally high values (~19.2 SI) [7]. For instance, Figure 1 [8] illustrates that the magnetic susceptibility of magnetite rock ranges from 0.6 SI to 1.0 SI in Huayang, Shaanxi Province, China. Although geological targets often possess both conductivity and permeability characteristics [9], the magnetic permeability in practical and numerical studies is usually assumed to be the free-space permeability ${\mu }_0$ ($4\pi \times {10}^{-7} \mathrm{H}/\mathrm{m}$). Many studies have shown that the influence of magnetic permeability in magnetite-rich regions cannot be ignored. Zhang and Oldenburg [10] demonstrated that assuming a susceptibility value of 0 for strong magnetized examples may result in an incorrect electrical model. Thomas [11] revealed that susceptibilities greater than 0.005 can significantly affect the EM responses in layered Earth models. Sasaki et al. [12] pointed out that ignoring susceptibility may lead to misinterpretation of models with high susceptibility values. Zhang et al. [13] successfully inverted two-dimensional data while simultaneously considering both conductivity and susceptibility, highlighting the significant impact of local magnetic anomalies on MT responses.

In recent years, there has been increased attention on magnetotelluric (MT) studies investigating the influence of conductive anisotropy [14,15,16], particularly in the development of three-dimensional (3D) forward modeling algorithms. Weidelt et al. [17] presented a staggered-grid finite difference algorithm for MT modeling in conductive anisotropic media; however, this algorithm is not convenient for dealing with irregular anomalies. Xu [18] proposed a node-based finite element (FE) method for MT modeling, but its accuracy is insufficient as it fails to satisfy the required condition of discontinuous normal electrical fields at electrical interfaces. To overcome this disadvantage of the node-based FE method, Xiao et al. [19] developed an edge-based FE method with rectangular meshes; however, this technique also faces challenges in dealing with irregular anomalies and terrain. Subsequently, Liu et al. [20] introduced a 3D adaptive FE method employing tetrahedron elements, enabling improved simulation of topography and irregular anisotropic anomalies. Additionally, Xiao et al. [21,22] contributed A-$\phi$ and T-$\Omega$ algorithms based on EM potentials. Zhou et al. [23] further advanced the field by presenting a divergence scheme rooted in the edge-based FE method, while Bai et al. [24] developed a mixed algorithm that combines nodal and edge-based approaches. However, there has been limited focus on incorporating conductive anisotropy and non-zero magnetic susceptibility. Xiao et al. [25] proposed a 3D MT modeling algorithm based on the edge-based FE method, utilizing the total field formula to obtain the EM field. Yu et al. [26] provided another 3D MT modeling scheme based on the A-$\phi$ formula, where the EM field is obtained by differentiating potential. Notably, neither of these algorithms have been parallelized, which limits their computational efficiency and scalability.

To address the aforementioned problems, we developed a robust and scalable multi-level parallelized 3D MMT modeling scheme that is capable of handling conductive and magnetic anisotropic anomalies. Based on the second field formula, the edge-based FE method is adopted to discretize the electric governing equation. To ensure accuracy, the algorithm enables the directly computation of secondary fields by solving equations with the SuperLU_dist direct solver [27]. To enhance efficiency, a multi-level parallel algorithm has been developed. This paper is organized as follows: Firstly, we provided a concise overview of the methodology for MMT modeling in conductive and magnetic anisotropic media. Subsequently, we presented crucial components necessary for achieving accurate and efficient implementation based on the secondary field approach. Following this, we conducted tests on different typical models, and we analyzed the impact of conductive anisotropy and magnetic susceptibility.