Magnetotelluric Modeling Research Articles

Three-dimensional edge-based finite element modeling of magnetotelluric data in anisotropic media with a divergence correction

When simulating a 3-D magnetotelluric (MT) response, the convergence rate slows at low frequencies. Previous works have shown that by implementing divergence correction, the convergence rate of the iterative solution can be accelerated, the computational accuracy can be enhanced. However, traditional correction schemes of MT are mostly based on isotropic modeling. Electrical anisotropy strongly affects the interpretation of MT data and anisotropic influence cannot be ignored. To our best knowledge, there has been no open literature to apply divergence correction in 3-D MT anisotropic modeling. In this study, based on the divergence-free property of the current density, a divergence correction technique is incorporated into the edge-based finite element (FE) method to speed up the modeling process for 3-D MT problems with axial anisotropic earth. By applying the divergence scheme, the computation time and iteration number for MT numerical simulation are largely reduced. The effectiveness of the anisotropic divergence correction technique on MT modeling is verified by numerical experiments.

Magnetotelluric Data Analysis Using Phase Tensor and Tipper Strike to Determine Geoelectrical Strike in “DKH” Geothermal Field

The magnetotelluric (MT) method is passive geophysical method that measures natural electromagnetic field to estimate the subsurface resistivity structure. Before conducting 2D MT data inversion, data analysis, dimensionality, including determining the direction of the geoelectrical strike must be performed. It is needed to fulfill the assumption of 2D MT modeling. 2D modeling of magnetotelluric data is an important part because the results can determine the quality of the subsurface interpretation. Therefore, this paper describes the determination of the dimensonality and geoelectrical strike using the phase tensor and tipper strike parameters. The invariant parameters of phase tensor including skew angle (β) and ellipticity are used to get the information about the dimensionality and geoelectrical strike. The results from the phase tensor are plotted in a rose diagram to determine the main direction of the geoelectrical strike. The phase tensor analysis contains ambiguity 90°, therefore tipper strike used to validate the analysis. From the analysis, the result shows the 1D or 2D structures of the MT data in high frequency, meanwhile the lower frequency, the majority of the data show 3D structures. The dominant geoelectrical strike is around N15°E and it is confirmed with the regional structure of the area. Determination of dimensionality and geoelectrical strike are applied to MT data from the geothermal field.

Three-dimensional magnetotelluric modeling using the finite element model reduction algorithm

This study presents a fast and efficient vector finite element method, based on the model reduction algorithm, to simulate the 3-D magnetotelluric response. Firstly, we derive a finite element governing equation based on electric field vector from Maxewell's equations Then we express the electric field by a function of frequency parameter to get a transfer function that can be approximated by the Krylov subspace model reduction algorithm. The transfer function after model reduction has far lower dimensions than the original one, and it can be directly solved by the Pardiso solver. Therefore, the computational efficiency is significantly improved. For validation, we compare the proposed algorithm with the 3-D direct finite element method (3DFEM) by a three-layer model, and the results show that the computation time of the former is less than 1/10 of the latter for the same accuracy. We also compare the electromagnetic response obtained by the proposed method with that obtained by two classical models, COMMEMI 3D-2A Model and Dublin Test Model 1. These results agree well, indicating that the proposed method is feasible and efficient for 3-D magnetotelluric forward modeling.

New insights into the High Agri Valley deep structure revealed by magnetotelluric imaging and seismic tomography (Southern Apennine, Italy)

We present an electrical resistivity model obtained from a 2D Magnetotelluric survey across a large sector of the Southern Apennine in the High Agri Valley (HAV), a NW-SE trending intra-mountain basin, with a high seismogenic potential. The intensive hydrocarbon exploitation (Val d'Agri oilfield) makes this area also affected by induced seismicity. In this HAV sector, the injection of salt-water in an unproductive disposal well (Costa Molina 2) causes localized swarms of microearthquakes; a second cluster of continuous induced seismicity is also observed SW of the Pertusillo Lake and it is associated to the seasonal fluctuations of the reservoir's water level.The major insight inferred from this study concerns a better understanding of the geological and tectonic framework in the HAV. The electrical resistivity model images the subsurface as conductive sedimentary sequences (Allochthonous Units) upon the carbonate Apulian Platform Unit characterized by higher resistivity values. Both these units appear composed of thrust-and-fold system deepening with larger wavelength anticlines N-E toward. Most of the structures identified in the magnetotelluric model are rather superficial and confined within the Allochthonous Units. A sudden break of the Apulian platform under the central part of the MT profile defines a conductive zone possibly associated to a major SW-dipping reverse fault or to several branches, as closely spaced thrust-sheets cutting eastern flanks of the Agri Valley.Additional information on the HAV deep structures comes from the joint interpretation of the resistivity model and a 3D seismic tomographic model obtained from the inversion of passive seismic data collected in the period 2002–2018. The availability of this elastic representation of the subsurface allowed us to perform a cluster analysis on the electrical resistivity and seismic P-wave velocity distribution within the subsoil. This joint quantitative interpretation unveiled new insights, otherwise hidden by individual models, on the subsurface structure distinguishing some rheological zones in terms of barriers and asperities.

Melting of subducted sediments reconciles geophysical images of subduction zones

Sediments play a key role in subduction. They help control the chemistry of arc volcanoes and the location of seismic hazards. Here, we present a new model describing the fate of subducted sediments that explains magnetotelluric models of subduction zones, which commonly show an enigmatic conductive anomaly at the trenchward side of volcanic arcs. In many subduction zones, sediments will melt trenchward of the source region for arc melts. High-pressure experiments show that these sediment melts will react with the overlying mantle wedge to produce electrically conductive phlogopite pyroxenites. Modelling of the Cascadia and Kyushu subduction zones shows that the products of sediment melting closely reproduce the magnetotelluric observations. Melting of subducted sediments can also explain K-rich volcanic rocks that are produced when the phlogopite pyroxenites melt during slab roll-back events. This process may also help constrain models for subduction zone seismicity. Since melts and phlogopite both have low frictional strength, damaging thrust earthquakes are unlikely to occur in the vicinity of the melting sediments, while increased fluid pressures may promote the occurrence of small magnitude earthquakes and episodic tremor and slip.

Synthesis of hydrothermal alteration, rock mechanics and geophysical mapping to constrain failure and debris avalanche hazards at Mt. Ruapehu (New Zealand)

ABSTRACT Composite volcanoes can progressively weaken through hydrothermal alteration, which may lead to volcano collapse, forming far-reaching debris avalanches. This work reviews and synthesises the type and extent of hydrothermal alteration on Mt. Ruapehu, New Zealand, from rock mechanical, mineralogical, hyperspectral imaging and aero-magnetic studies for input into the first debris avalanche models of Mt. Ruapehu. Mt. Ruapehu shows surface weathering processes, forming goethite, hematite and phyllosilicate assemblages, and hydrothermally altered rock, with phyllosilicate, Fe- oxides and sulphates, pyrite, jarosite, alunite, anhydrite and sulphur. The surface alteration mapped by hyperspectral remote sensing is linked to the deep-seated (≤500 m) alteration by aero-magnetic inversion modelling. Textural analysis and mapping show that alteration weathering and porosity distribution are a result of the both the primary lithology, age, and location relative to active or fossil hydrothermal systems. The synthesis of the aero-magnetic data, previous magnetotelluric models, hyperspectral and hydrothermal mapping, and rock mechanics studies enabled the construction of a schematic model to provide inputs for modelling potential flank failure zones. These models suggest that the southwest and west flanks of Mt. Ruapehu have the highest probability of collapsing and generating debris avalanches.

Subsurface geoelectrical structure from 3-d inversion of magnetotelluric data of Gisenyi geothermal field, western part of Rwanda

In this paper, a set of Magnetotelluric (MT) data collected from 69 stations at the Gisenyi geothermal field in the western part of Rwanda were used to derive a 3D electrical resistivity distribution model to investigate deep geologic structures related to the geothermal resources. The MT method is a geophysical technique commonly used in geothermal exploration projects. The aim of this study is to investigate subsurface electrical resistivity distribution in the Gisenyi geothermal area to gain an understanding of the process related to geothermal activities and to model deep crustal structures. MT data analysis using phase tensors revealed the nature of the subsurface structures of the geothermal area prior to performing a 3D MT inversion. The dimensionality analysis results revealed a 3D conductivity structure in the deeper portion. The ModEM 3D MT inversion code was used to infer the subsurface electrical distribution from the MT responses. The recovered 3D resistivity distribution model showed a conductive region (<10 Ω.m) at deeper depths beneath a thick high-resistivity zone with a resistivity more than 100 Ω.m. The 3D model indicates patches of high and low resistivity at shallow depths related the presence of recent volcanic lavas and basement rocks from the surface to depths of a few hundred meters. The MT model revels important subsurface structures that can be related to the geothermal resource in the area. The deep-seated conductive body is presumably associated to the heat source of the Gisenyi geothermal system.

Deep Structure of the Lithosphere in the Central Tien Shan along the Son-Kul Magnetotelluric Sounding Profile

New results of a detailed study of the deep structure of the Central Tien Shan along the Son-Kul magnetotelluric (MT) profile crossing the Son-Kul Lake are reported. Based on the results of magnetotelluric data modeling, the regional and local geoelectric anomalies in the lithosphere are studied and their quantitative characteristics are given. Geological interpretation of the geoelectric cross-section was carried out, which supported the existing ideas about the block–hierarchical structure of the upper part of the Earth’s crust. This corresponds to the tectonophysical concepts of the sequential inserted subordination of large and smaller elements of the zone–block structure consisting of stable blocks and limiting mobile zones, which are distinguished by the high dislocation of the geological substrate. The integral pattern of the distribution and morphology of zones of high electrical conductivity in this segment of the Central Tien Shan crust may reflect discretely localized palm tree–type structures associated with the evolution of transgressive suture zones of localized deformation during the Hercynian and Alpine tectogenesis.

Three-dimensional magnetotelluric modeling in general anisotropic media using nodal-based unstructured finite element method

Magnetotelluric sounding is an efficient and economical method for detecting Earth's deep structure because of its large surveying depth. We use a finite element approach to solve magnetotelluric fields in three-dimensional general anisotropic media based on the vector–scalar potentials and unstructured meshes. The implementation of a three-dimensional unstructured mesh generator TetGen with high-quality local mesh refinement can easily handle complicated models with high performance. The electromagnetic field equations are expressed in terms of Coulomb-gauged vector–scalar potentials, and the nodal-based finite element method can be used to compute the electromagnetic fields with high efficiency. The symmetric successive over-relaxation preconditioned conjugate gradient solver is used to solve the resulting linear equation system. The magnetotelluric responses in one-, two-, and three-dimensional medias are computed through our proposed approach. These results are then compared with the analytical solutions, two- and three-dimensional adaptive finite element solutions. Numerical examples demonstrate that these results are in good agreement. Because an unstructured mesh is incorporated in this approach, it has great potential for simulating the magnetotelluric responses in geometrically complicated situations.

Two-Dimensional magnetotelluric modelling based on the numerical manifold method

As a numerical method based on the Galerkin variation, the numerical manifold method (NMM) adopts a dual covering system, which is convenient for constructing high-order manifold elements and performing adaptive analysis. In this paper, the NMM is used to simulate the electromagnetic response based on the structure of the undulating terrain existing in actual geophysical exploration. The relevant system equations are derived, and a penalty function method is proposed for dealing with boundary conditions. The simulation results validate the effectiveness and accuracy of this method. Through the simulation of ground barriers and ridge models, the effect of undulating terrain on the apparent resistivity and phase is investigated, which is potentially useful for examining the complex geoelectric model characteristics as well as for terrain correction of the detection results.

Lithospheric resistivity structure of the 2001 Bhuj earthquake aftershock zone

SUMMARY The Bhuj area, in the Kutch region of western India, is a unique intraplate seismic zone in the world where aftershock activity associated with a large magnitude earthquake (7.7 Mw Bhuj earthquake on 26 January 2001) has persisted over a decade and up till today. We studied the lithospheric resistivity structure of the Bhuj earthquake aftershock zone to gain more insight into the structure and processes influencing the generation of intraplate seismicity in broad and, in particular, to detect the deep origin and upward migration channels of fluids linked to the crustal seismicity in the area. A lithospheric resistivity model deduced from 2-D and 3-D inversions of long-period magnetotelluric (MT) data shows low resistive lithospheric mantle, which can be best explained by a combination of a small amount of interconnected melts and aqueous fluid in the upper mantle. The MT model also shows a subvertical modestly conductive channel, spatially coinciding with the Kutch Mainland Fault, which we interpret to transport fluids from the deep lithosphere to shallow crust. We infer that pore pressure buildup aids to achieve the critical stress conditions for rock failure in the weak zones, which are pre-stressed by the compressive stress regime generated by ongoing India–Eurasia collision. The fluidized zone in the upper mantle beneath the area perhaps provides continuous fluid supply, which is required to maintain the critical stress conditions within the seismogenic crust for continued seismicity.

Integrating Magnetotelluric and Seismic Images of Silicic Magma Systems: A Case Study From the Laguna del Maule Volcanic Field, Central Chile

AbstractImaging silicic systems using geophysics is challenging because many interrelated factors (e.g., temperature, melt fraction, melt composition, geometry) can contribute to the measured geophysical anomaly. Joint interpretation of models from multiple geophysical methods can better constrain interpretations of the subsurface structure. Previously published resistivity and shear wave velocity (Vs) models, derived separately from magnetotelluric (MT) and surface wave seismic data, respectively, have been used to model the restless Laguna del Maule Volcanic Field, central Chile. The Vs model contains a 450 km3 low‐velocity zone (LVZ) interpreted as a region with an average melt fraction of 5–6%. The resistivity model contains a conductor (C3) interpreted as a region with a melt fraction &gt;35%. The spatial extents of the LVZ and C3 overlap, but the geometries and interpretations of these features are different. To resolve these discrepancies, this study investigates the resolution of the MT data using hypothesis testing and constrained MT inversions. It is shown that the MT data are best fit with discrete conductors embedded within the larger LVZ. The differences between the MT and seismic models reflect resolution differences between the two data sets as well as varying sensitivities to physical properties. The MT data are sensitive to smaller volumes of extractable mush that contain well‐connected crystal‐poor melt (C3). The seismic data have lower spatial resolution but image the full extent of the poorly connected crystal‐rich magma storage system. The combined images suggest that the LdMVF magma plumbing system is thermally heterogeneous with coexisting zones of warm and cold storage.

Validity of the dispersion relations in magnetotellurics. Part II: synthetic and field data

The impedance tensor and tipper are shown to be non-causal in some classes of 2-D and 3-D magnetotelluric (MT) models, including those widely encountered off shore. At the same time, examination of a large database of onshore MT records yields only a handful of non-causal tensors, being non-causal due to the influence of nearby artificial conductors rather than some exotic geological conditions. This result implies that the actual chance to encounter a non-causal impedance tensor or tipper in terrestrial MT exploration is vanishingly small, thus securing the universal applicability of the dispersion relation technique for quality assessment of the mentioned transfer functions.

Modular finite volume approach for 3D magnetotelluric modeling of the Earth medium with general anisotropy

Electrical property of the Earth should be best depicted by anisotropic 3D models, while many numerical issues remain to be solved before it becomes practical, such as implementation complexity and computational efficiency. Magnetotellurics (MT) is the primary method to explore electrical property of the deep Earth, and MT forward modeling with 3D anisotropy has been investigated primarily by finite element methods (FEM) and finite difference methods (FDM). Here, we present a new implementation of finite volume method (FVM) accounting for general anisotropy. Theoretically, FVM can combine the advantages of FEM in mesh flexibility and superior accuracy, and the advantages of FDM in straightforward implementation and economic computation load. For demonstration, a modular framework to approximate Maxwell's equations on a structured staggered-grid (mesh flexibility ignored), which solves electric fields tangential to cell edges, is implemented in an object-oriented paradigm. Synthetic examples are presented to verify its accuracy and feasibility, which includes (quasi-) analytic solution of simplified models and numerical solution of a dyke model. A variant of COMMEMI 3D-2 model with general anisotropy is simulated to probe the characteristic behaviors of EM fields under such circumstance. This approach makes the dependence of model operators on anisotropic conductivity parameters explicit and can simplify further sensitivity derivation, and takes a step forward for the feasible development of MT inversion algorithms concerning 3D anisotropy.

MATE: An Analysis Tool for the Interpretation of Magnetotelluric Models of the Mantle

AbstractInterpretation of electrical conductivity anomalies observed in magnetotelluric models provides an important opportunity to understand the nature of the lithospheric mantle and its dynamics. Over the course of the last two decades, a great number of experimental petrology studies have been carried out which can be utilized to construct electrical conductivity distribution models for a given composition and geotherm. We have developed an open‐source software (MATE, Mantle Analysis Tool for Electromagnetics) with an easy‐to‐use graphical interface that creates such theoretical models. The program is developed in such a way that additional effects and models can be added very easily. To investigate the conductivity distribution of the cratonic mantle, a series of experiments was made. Results indicate that it is of utmost importance to analyze the magnetotelluric models using accurate compositions, water distributions, and geometric models. Hence, using only olivine conductivity models can lead to erroneous interpretations of both conductivity and estimated water content. Analysis of the potential causes for conductive anomalies shows that the upper and lower lithospheric mantle can be interpreted separately with the transition between them at 75–125 km. Conductive anomalies in the upper lithospheric mantle (&lt;1,000 Ωm) are likely to be explained by the existence of well‐connected minor phases associated with metasomatic fluids, whereas in the lower lithospheric mantle, hydration and/or well‐connected minor phases (e.g., phlogopite or amphibole) can explain conductive anomalies.

Three-dimensional magnetotelluric numerical simulation of realistic geologic models

We have developed a workflow for constructing realistic mesh-based magnetotelluric (MT) models from 3D geologic models. The routine is developed for unstructured meshes that adapt to the complex shapes of geologic bodies including 3D surfaces and volumes in realistic modeling scenarios. The methodology is applied to the complexly altered Lalor volcanogenic massive sulfide deposit in Manitoba, Canada. The host rock envelope of the Lalor deposit is compartmentalized into lithostratigraphic units leading to a watertight model. This model then is meshed into unstructured tetrahedral meshes suitable for synthetic geophysical modeling of the MT method. Subsequently, two 3D resistivity models are generated from wireline logs: (1) a host rock background model in which each tetrahedral cell is attributed with the average resistivity of each lithostratigraphic unit and (2) a heterogeneous background-ore model in which the resistivity values of the cells are resampled from a 3D curvilinear grid model, generated by computing sequential Gaussian simulations from the resistivity data for each unit of a 3D lithofacies model produced by categorical kriging. To calculate the synthetic response of this model for MT, a numerical-modeling code is developed based on solving the vector-scalar potential formulation of the electromagnetic diffusion equation using the finite-element method on unstructured meshes. After validating the numerical method for the Commemi test model, the MT response of the Lalor model is investigated. A reasonable agreement is observed between the survey field data and the data synthesized from our constructed heterogeneous model. Using an investigation of the inductive and galvanic parts, we conclude with the ideal frequency range for detecting the ore deposit. We also conclude with and visualize the importance of regional-scale alteration zones around the ore deposits and model inhomogeneities in boosting the detectability of the ore formations through feeding electrical currents as a result of galvanic field dominance at depth.

Shallow structure of Los Humeros (LH) caldera and geothermal reservoir from magnetotellurics and potential field data

SUMMARYThis study focuses in the analysis of the internal structure of the upper 3 km of Los Humeros (LH) caldera and the relation of electrical and hydrothermal anomalies. For this purpose, we measured, processed and interpreted 78 broad-band magnetotelluric (MT) soundings. We performed a 3-D inversion of the data set (ModEM) using all MT soundings, although only half of the available frequencies per sounding due to limited computed power. We also carried out the 2-D inversions (NLCG) of the invariant determinant along two orthogonal profiles (EW and NS) crossing the caldera structure; their comparison yields similar resistivity and structural models results. The resistivity modelling is complemented with the results of a joint 3-D inversion of an accurate gravity database of 720 stations, and total field aeromagnetic data (SGM) from the caldera crater. The combined results provide novel details about the structure of the shallow geothermal reservoir of the resurgence caldera complex hosting the active hydrothermal system. Density and resistivity models show the existence of a composed crater basin structure separated by an EW high-density structure; the northern basin is associated to the LH crater, whereas the southern basin associates to the emergent Los Potreros (LP) caldera basin. The magnetization model indicates that there is a common source for the magnetic volcanic products observed at the caldera surface, and that the LP fault is the more magnetized fault of the geothermal system. The propylic zoning under the geothermal field, which according to the MT model results has resistivities above ∼100 Ω-m, was extrapolated using this and additional criteria to obtain the distribution of other hypothetical propylitic zones of hydrothermal potential.

An efficient parallel algorithm for 3D magnetotelluric modeling with edge-based finite element

Three-dimensional magnetotelluric modeling algorithm of high accuracy and high efficiency is required for data interpretation and inversion. In this paper, edge-based finite element method with unstructured mesh is used to solve 3D magnetotelluric problem. Two boundary conditions—Dirichlet boundary condition and Neumann boundary condition—are set for cross-validation and comparison. We propose an efficient parallel algorithm to speed up computation and improve efficiency. The algorithm is based on distributed matrix storage and has three levels of parallelism. The first two are process level parallelization for frequencies and matrix solving, and the last is thread-level parallelization for loop unrolling. The algorithm is validated by several model studies. Scalability tests have been performed on two distributed-memory HPC platforms, one consists of Intel Xeon E5-2660 microprocessors and the other consists of Phytium FT2000 Plus microprocessors. On Intel platform, computation time of our algorithm solving Dublin Test Model-1 with 3,756,373 edges at 21 frequencies is 365 s on 2520 cores. The speedup and efficiency are 1609 and 60% compared to 100 cores. On Phytium platform, scalability test shows that the speedup from 256 cores to 86,016 cores has been increased to 11,255.

Neoproterozoic amalgamation and Phanerozoic reactivation of Central/Western Hoggar (Southern Algeria, Tuareg Shield) lithosphere imaging using Magnetotelluric data

Hoggar, as part of the Tuareg shield, is a key area in the assembly of Western Gondwana during the Neoproterozoic. To better understand the tectonic process, 61 broadband magnetotelluric (MT) stations were acquired in Central/Western Hoggar to image the underlying lithosphere. We performed 2-D and 3-D inversions of the impedance tensor and the vertical magnetic field transfer function using data from three profiles. The resistive crust is interrupted by several conductors related to the known shear zones of the area. Some conductors can be associated with Cretaceous-Cenozoic fluid flow combined with mineralization (graphite, magnetite or more locally gold). MT models suggest three distinct blocks. (1) An oriental block composed of Tefedest, Laouni, Aouilène and Silet juvenile terrane characterized by a moderately conductive lithospheric mantle (LM). (2) A central block formed by the In Tedëini and Tin Zaouatene terranes, previously considered as Neoproterozoic juvenile, exhibit resistive characteristics typical of an old continental LM. Finally, (3) an occidental block corresponding to the In Ouzzal terrane with a highly conductive LM. The oriental and occidental blocks were strongly remobilized during the Neoproterozoic, and perhaps even locally during the Cretaceous and Cenozoic. The main shear zones of the study area are (1) the Iskel shear zone, a major suture zone separating the oriental and central blocks. The MT model suggests that the geometry of this lithospheric shear zone is characterized by an eastward dipping angle, according to the central block a passive margin character where Silet arc was accreted. Consequently, east from Silet a collision with the oriental block may have been in an active margin environment. (2) The Adrar fault and the West Ouzzalian shear zone correspond to terrane boundaries of the In Ouzzal. They are associated with deep-reaching high conductivity, consistent with large scale suture zones. The resistivity models provide new geometric constraints for the terranes and highlight new Pan-African and Phanerozoic tectonic features.

A Chebyshev pseudo-spectral approach for simulating magnetotelluric TM-mode responses on 2D structures

We introduce a Chebyshev pseudo-spectral (CPS) approach for 2D TM-mode magnetotelluric modelling, which was the first applied on the magnetotellric data in 2D model. In this method, 2D partial differential equations are solved using the Chebyshev differentiation matrix over the Gauss-Chebysev-Lobatto points. In order to test computational efficiency and its accuracy, we compare the simulatation of this method to the analytical solution and finite element (FE) and finite difference (FD) methods. The comparison of CPS numerical solution and corresponding analytical solution on a homogeneous half-space model and a layered model show that the CPS approach can provide an accurate modelling result. The discrete nodes for modelling were also analyzed on two models, which suggest that a few nodes can produce modelling results with a high accuracy. Especially for 2D and 3D modelling problems, the reduction in the number of nodes can improve the computational efficiency significantly. However, the CPS approach will lead to an unsymmetrical and non-sparse matrix, which is not favored by the iterative solver. We further compare the numerical result of the CPS method to the FD method and the averaged numerical solutions of the COMMEMI on the COMMEMI 2D-1 model. We find the numerical result of CPS method is closer to the averaged numerical solutions of the COMMEMI than that of observed by the FD. Finally, a realistic smooth model is considered to compare the modelling results using the CPS method to the published data. All these comparison results suggest that the CPS method can not only give a high accuracy for modelling results and also provide more detail information. Especially, a few nodes can be required in this method relative to the FD and FE method, which can decrease the relative errors. We then deduce that CPS method might be an alternative approach to simulate the magnetotelluric responses in 2D structures.