Magnetotelluric method: theoretical deduction, analytical solutions and numerical modeling in a one-dimensional Earth

Método magnetotelúrico: dedução teórica, soluções analíticas e modelagem numérica em uma Terra 1D

Abstract

The magnetotelluric (MT) method is a passive geophysical technique that uses natural variations in electromagnetic fields to investigate the Earth’s subsurface resistivity structure. This article has a didactic objective: to guide students of physics and geophysics, at any level, through the theoretical foundations and mathematical derivations underlying the 1D MT method. Starting from Maxwell’s curl equations in the frequency domain, we derive wave-like differential equations for the electric and magnetic field components in a conductive medium. These equations lead to analytical solutions that describe the attenuation and phase shift of electromagnetic waves as they propagate through layered Earth models. From these results, we deduce expressions for the MT impedance and phase. The step-by-step presentation is complemented by simple computational implementations in Python for 1D models, enabling readers to connect the theoretical concepts with practical applications. By presenting the methodology in a clear and structured way, this work aims to support students and educators in developing a deeper understanding of the magnetotelluric method and its role in subsurface exploration.

Keywords:
Geophysics; Geophysical Method; Magnetotelluric Method; Maxwell Equations; Forward Modeling

Resumo

O método magnetotelúrico (MT) é uma técnica geofísica passiva que utiliza variações naturais em campos eletromagnéticos para investigar a estrutura de resistividade do subsolo da Terra. Este artigo tem um objetivo didático: guiar estudantes de física e geofísica, em qualquer nível, através dos fundamentos teóricos e derivações matemáticas que fundamentam o método de MT unidimensional. Partindo das equações rotacionais de Maxwell no domínio da frequência, derivamos equações diferenciais ondulatórias para as componentes dos campos elétrico e magnético em um meio condutor. Essas equações levam a soluções analíticas que descrevem a atenuação e o deslocamento de fase das ondas eletromagnéticas à medida que se propagam através de modelos terrestres em camadas. A partir desses resultados, deduzimos expressões para a impedância e a fase da MT. A apresentação passo a passo é complementada por implementações computacionais simples em Python para modelos unidimensionais, permitindo que os leitores conectem os conceitos teóricos com aplicações práticas. Ao apresentar a metodologia de forma clara e estruturada, este trabalho visa apoiar estudantes e educadores no desenvolvimento de uma compreensão mais aprofundada do método magnetotelúrico e seu papel na exploração do subsolo.

Palavras-chave:
Geofísica; Método Geofísico; Método Magnetotelúrico; Equações de Maxwell; Modelagem Direta

1. Introduction

Geophysics is an applied science that investigates subsurface structures by analyzing the physical properties of rocks, using a variety of instrumental techniques known as geophysical methods. These methods are classified according to the main physical property they measure, such as potential field methods (e.g., gravity and magnetics), seismic methods, and electromagnetic methods. The magnetotelluric (MT) method is part of a group of electromagnetic techniques that measure natural variations in the Earth’s electromagnetic fields.

The MT method was independently introduced by Tikhonov [¹] and Cagniard [²]. It is based on the simultaneous measurement of the horizontal components of the natural electric and magnetic field variations at the Earth’s surface to infer the subsurface geoelectric structure. Initially developed for one-dimensional models, where electrical resistivity varies solely with depth-the method has since undergone significant advancements, enabling its application to more complex two- and three-dimensional geological settings.

Although the MT problem is well established and has been extensively studied, this paper consolidates prior knowledge by presenting three different approaches to its resolution for 1D case. This paper provides a concise overview of the theoretical foundations and governing equations of the MT method. It also includes the derivation of analytical expressions for electromagnetic wave diffusion, leading to the formulation of apparent resistivity and phase equations. In addition, compute modeling equations are provides through simplified Python codes, based on theoretical formulations. All codes are made available for download, aiming to support students and researchers in understanding and applying the MT method.

2. Magnetotelluric Method

The MT method provides subsurface imaging with an intermediate resolution between reflection seismic and potential field methods [³]. Among geophysical techniques, MT stands out for its exceptional investigation depth, capable of probing from tens of meters to several hundred kilometers [⁴].

Eletromagnetic methods enable resource prospecting by measuring the effects of eletric current flowing in the surface. The MT method depends on naturally occurring varying electromagnetic fields present at the Earth’s surface. Understanding the sources of these electromagnetic signals is fundamental to interpreting MT data as they determine the frequency range and, consequently, the depth of investigation.

Electromagnetic signals used in MT span a wide range of frequencies, each associated with different physical processes as illustrated on Figure 1. High-frequency signals ($f>1\text{Hz}$) are generated in the region between the Earth’s surface and the ionosphere, primarily due to lightning discharges, electromagnetic phenomena commonly referred to as sferics. These signals dominate the audio-frequency range and are useful for shallow investigations.

Figure 1
Schematic illustration of natural electromagnetic sources used in magnetotellurics.

In contrast, low-frequency signals ($f<1\text{Hz}$) originate from the interaction between charged particles in the solar wind and the Earth’s magnetosphere and ionosphere, producing ultra-low frequency (ULF) variations that enable deep subsurface exploration.

A key aspect of MT is that it allows geophysicists to investigate multiple death layers by analyzing signals in different frequency bands. The ability of electromagnetic waves to penetrate the Earth depends on their frequency: lower frequencies reach deeper regions, while higher frequencies are limited to near-surface layers. This phenomenon is described by the concept of skin depth, which will be explored in the next sections. Consequently, the MT method is commonly subdivided into four principal frequency bands, Figure 2, each associated with specific signal sources and depth ranges.

Each frequency band is selected depending on the target depth and geological context. Modern MT surveys often combine measurements from multiple bands to construct a comprehensive image of the subsurface across several depth scales.

Long-Period MT (MTLP): Covers ultra-low frequencies typically below 1 Hz, down to ${10}^{-6}\text{Hz}$. These signals originate from global-scale interactions between the solar wind and the Earth’s magnetosphere and ionosphere. MTLP is particularly suited for crustal studies capturing both shallow and intermediate-depth structure.

Figure 2
Subdivisions of the magnetotelluric method according to frequency range: MTLP (long-period), MTBL (broadband), AMT (audio), and RMT (radio).

Broadband MT (MTBL): spanning approximately from ${10}^{-3}\text{Hz}$ to ${10}^3\text{Hz}$, broadband MT is the most widely used range in crustal studies. Capturing both shallow and intermediate-depth structures.

Audio-Magnetotellurics (AMT): Operating in the range of roughly ${10}^1$ to ${10}^4\text{Hz}$, AMT is dominated by high-frequency electromagnetic signals caused by lightning discharges (sferics). This band is well-suited for resolving shallow subsurface features, such as sedimentary basins, aquifers, or near-surface faults.

Radio-Magnetotellurics (RMT): Extending from approximately $10$ to ${10}^6\text{Hz}$, RMT uses man-made radio transmitters as the main source of signal. It provides high-resolution imaging of very shallow structures, typically within the first tens of meters below the surface.

Although RMT operates in the frequency range from approximately $10$ to ${10}^6\text{Hz}$, using man-made radio transmitters as the primary source of signal, this article focuses specifically on natural-source methods, such as those driven by ionospheric and magnetospheric variations.

Unlike other MT subdivisions, which relies on naturally occurring electromagnetic fields, RMT introduces controlled artificial sources, which may require distinct mathematical formulations, especially in terms of boundary conditions and far-field approximations.

Similarly, methods like Controlled-Source Audio Magnetotellurics (CSAMT) also fall outside the scope of this work, as they depend on transmitter-receiver setups with designed signal characteristics. These techniques are valuable for high-resolution imaging of shallow subsurface layers, but their reliance on non-natural sources introduces different theoretical assumptions.

For these reasons, RMT and CSAMT are not addressed here, and the theoretical developments, derivations, and modeling presented in this article are exclusively based on the natural-source MT method.

2.1. Boundary conditions

Nabighian [⁵] affirms that the resultant field of electromagnetic problems is the sum of the primary and secondary fields. Each of the fields must satisfy Maxwell’s equations, or equations derived therefrom, plus appropriate conditions to be applied at boundaries between the homogeneous regions involved in the problem, e.g., at the air–earth interface.

Boundary conditions can be divided into two main categories: those derived from the normal components and those from the tangential components of the electromagnetic fields. Each category leads to specific constraints on the behavior of the fields $\overrightarrow{B}$, $\overrightarrow{D}$, $\overrightarrow{E}$, $\overrightarrow{H}$, and $\overrightarrow{J}$, as well as on the scalar potentials $V$ and $U$, and the vector potentials $\overrightarrow{F}$ and $\overrightarrow{A}$. The main boundary conditions are summarized schematically in Figure 3. For detailed information, check Griffiths [⁶] section 7.3.6.

Figure 3
Flowchart of boundary condition categories used in magnetotelluric theory.

2.2. Normal $\overrightarrow{B}$

The normal component of the magnetic field does not present discontinuity when passing from one medium to another. This is due to the absence of magnetic monopoles and follows directly from Maxwell’s equation (7).

2.3. Normal $\overrightarrow{D}$

The normal component of the electric displacement field $\overrightarrow{D}$ is generally discontinuous across an interface due to the accumulation of surface charge density ${\rho }_s$. This condition is expressed as $D_{n2}-D_{n1}={\rho }_s$ as presented on Figure 4.

Figure 4
Electric displacement field in two mediums.

In MT method, it is common to assume the absence of free surface charge at the interfaces between geological layers, leading to the simplified boundary condition. However, this condition becomes particularly relevant in the presence of dielectric interfaces, such as the transition between air and soil, or between rock and a fluid with different electrical conductivity. In such cases, discontinuities in the normal component of the electric displacement field $\overrightarrow{D}$ may arise due to surface charge accumulation.

2.4. Tangencial $\overrightarrow{E}$

The tangential condition for the electric field is expressed as ${\overrightarrow{E}}_{t1}={\overrightarrow{E}}_{t2}$. That is, the tangential electric field is continuous across an interface, provided there is no abrupt variation in potential or surface current sources.

2.5. Tangencial $\overrightarrow{H}$

The tangential condition for the magnetic field is given by ${\overrightarrow{H}}_{t1}={\overrightarrow{H}}_{t2}$, this means that the tangential component of the magnetic field is continuous across an interface, as long as there is no surface current density present.

In geophysical practice, particularly in the MT method, it is commonly assumed that the magnetic permeability is equal to that of free space, i.e., $\mu ={\mu }_0$, since most geological materials are not strongly magnetic. This leads to ${\mu }_1={\mu }_2$ and $H_{n1}=H_{n2}$. Thus, in MT modeling, this boundary condition is often simplified to $H_n$ continuous.

However, this simplification may not hold in cases involving materials with significant magnetization, such as iron-rich rocks or models with strong contrasts in magnetic permeability $\mu$.

2.6. Current density $\overrightarrow{J}$

The boundary condition for the current density vector $\overrightarrow{J}$ is expressed as $J_{n1}=J_{n2}$. That is, the normal component of the current density is continuous across an interface. Strictly speaking, this result is valid for direct currents. However, it is generally acceptable for Earth materials up to frequencies of approximately ${10}^5$ Hz where the displacement current can be neglected [⁵],

(1) $$\left|i\omega ϵ\overrightarrow{E}\right|\ll \left|\sigma \overrightarrow{E}\right|.$$

2.7. Scalar potentials

For static fields, the eletrical scalar potential and magnetic one can be matematically described by,

(2) $$\overrightarrow{E}=-\nabla V,$$

(3) $$\overrightarrow{H}=-\nabla U.$$

The scalar electric potential $V$ and the scalar magnetic potential $U$ must be continuous across an interface. This ensures that the work done in moving a small test charge or magnetic dipole across the boundary is path independent, where $V_1=V_2$ and $U_1=U_2$.

In geophysical applications, this condition is relevant in low-frequency or static approximations, where fields are derived from scalar potentials rather than from time-varying vector fields.

2.8. Vector potentials

The Schekkunoff vector potentials are more complex consideration of fields that facilitate the boundary condition applications by separating the variables of the system in symmetry problems. There are two types of vector potentials: the magnetic and electric.

Chart 1 Types of sources for Schelkunoff vector potentials.

Magnetic Vector Potential $\overrightarrow{\mathbf{F}}$	Electric Vector Potential $\overrightarrow{\mathbf{A}}$
Magnetic dipole perpendicular to a layered Earth Infinite cylinder in a uniform inducing magnetic field	Electric dipole perpendicular to a layered Earth Cylinder in a uniform electric field
Sphere in a uniform inducing magnetic field Horizontal magnetic dipole	Sphere in a uniform electric field Horizontal electric dipole over a layered Earth

Dipoles perpendicular to a layered Earth allow vertical stratification and 1D modeling; The cylindrical and spherical geometries are used to derive analytical solutions in 2D and 3D cases, serving as benchmarks for validating numerical codes. Also the horizontal dipoles better represent field conditions in near-surface geophysics, such as in magnetotelluric or controlled-source EM methods.

Understanding the physical origin of each potential and its associated source helps clarify why different potentials are preferred in different modeling contexts. This conceptual distinction also underlies many numerical implementations, such as finite-difference or finite-element schemes, which discretize these potentials rather than the fields directly.

3. Maxwell’s Equations

All electromagnetic phenomena are governed by empirical Maxwell’s equations. These are first-order linear differential equations that can be combined by the empirical constitutive relations with the aid of algebric manipulation, which reduce the number of basic vector functions from five to two [²]. Where Faraday’s Law is described by,

(4) $$\nabla \times \overrightarrow{E}=-\frac{\partial \overrightarrow{B}}{\partial t},$$

stating that a time-varying magnetic flux induces an opposing electric field. Ampère’s Law is formulated as:

(5) $$\nabla \times \overrightarrow{H}=\frac{\partial \overrightarrow{D}}{\partial t}+\overrightarrow{J},$$

describing how magnetic fields arise from electric currents and time-varying electric flux. Mathematically, Gauss’s Law takes the form:

(6) $$\nabla \cdot \overrightarrow{D}=\rho ,$$

that relates electric fields to the presence of electric charges. Gauss’s Law for magnetism is represented as:

(7) $$\nabla \cdot \overrightarrow{B}=0,$$

which formalizes the nonexistence of magnetic monopoles.

Where $\overrightarrow{E}$ is the electric field intensity vector (in V/m), $\overrightarrow{B}$ is the magnetic flux density vector (in Wb/m² or T), $\overrightarrow{D}$ is the electric displacement vector (in C/m²), $\overrightarrow{H}$ is the magnetic field intensity vector (in A/m), $\rho$ is the electric charge density (in C/m³) and $\nabla$ is the rotational operator:

(8) $$\nabla =({\displaystyle \frac{\partial }{\partial x}},{\displaystyle \frac{\partial }{\partial y}},{\displaystyle \frac{\partial }{\partial z}}).$$

4. Constitutive Relationships

Constitutive relationships or material equations describe the electrical behavior of materials under the influence of electromagnetic fields, supplementing Maxwell’s equations to determine unambiguously vectorial fields.

(9) $$\overrightarrow{D}=\epsilon \overrightarrow{E},$$

(10) $$\overrightarrow{B}=\mu \overrightarrow{H},$$

(11) $$\overrightarrow{J}=\sigma \overrightarrow{E},$$

where $ϵ$ is called dielectric permittivity and is a tensor attributed to the motion of electrons, nuclei, and polar molecules due to the application of an electric field and its tendency for electrical balance [⁵]. $\mu$ is defined as magnetic permeability, which describes a material’s facility to allow a magnetic field formation in its interior. $\sigma$ is known as conductivity, characterizing the capacity of a material to transport an electrical charge when a potential difference exists.

The Equation (11) is known as Ohm’s Law. To simplify the electromagnetic problem, induced magnetization is neglected, and the magnetic permeability is assumed to be that of free space: $\mu ={\mu }_0=4\pi \times {10}^{-7}N/A^2$.

5. Diffusion Equations in Frequency Domain

Considering the Eqs. (4) and (5) in the frequency domain, they can be respectively rewritten as:

(12) $$\nabla \times \overrightarrow{E}=-i{\mu }_0\omega \overrightarrow{H},$$

(13) $$\nabla \times \overrightarrow{H}=\sigma \overrightarrow{E}+i\omega \epsilon \overrightarrow{E},$$

being $\partial /\partial t=i\omega$ and $\omega =2\pi f$. Applying curls in both Eqs.,

(14) $$\nabla \times \nabla \times \overrightarrow{E}=\nabla \times \left(-i{\mu }_0\omega \overrightarrow{H}\right),$$

(15) $$\nabla \times \nabla \times \overrightarrow{H}=\nabla \times \left(\sigma \overrightarrow{E}+i\omega ϵ\overrightarrow{E}\right).$$

Using the vector identity,

(16) $$\begin{array}{l}\nabla \times \left(\nabla \times \overrightarrow{A}\right)=\nabla \left(\nabla \cdot \overrightarrow{A}\right)-\nabla \cdot \nabla \overrightarrow{A},\\ =\nabla \left(\nabla \cdot \overrightarrow{A}\right)-\nabla \cdot \nabla \overrightarrow{A},\end{array}$$

we can rewrite Eqs. (14) and (15) as,

(17) $$\nabla \left(\nabla \cdot \overrightarrow{E}\right)-{\nabla }^2\overrightarrow{E}=-i{\mu }_0\omega \nabla \times \overrightarrow{H},$$

(18) $$\nabla \left(\nabla \cdot \overrightarrow{H}\right)-{\nabla }^2\overrightarrow{H}=\sigma \nabla \times \overrightarrow{E}+i\omega ϵ\nabla \times \overrightarrow{E}.$$

Knowing that there is no accumulation of electric charges, we have $\nabla \cdot E=0$, and $\nabla \cdot H=0$ according to Gauss’s Law.

(19) $${\nabla }^2\overrightarrow{E}=i{\mu }_0\omega \nabla \times \overrightarrow{H},$$

(20) $$-{\nabla }^2\overrightarrow{H}=\sigma \nabla \times \overrightarrow{E}+i\omega ϵ\nabla \times \overrightarrow{E}.$$

Replacing Ampere’s and Faraday’s Laws Eqs. (12) and (13) on Eqs. (19) and (20). To Eq. (19),

(21) $${\nabla }^2\overrightarrow{E}=i{\mu }_0\omega \left(\sigma \overrightarrow{E}+i\omega ϵ\overrightarrow{E}\right),$$

(22) $$\Rightarrow {\nabla }^2\overrightarrow{E}=\underset{\text{diffusion term}}{\underbrace{i{\mu }_0\omega \sigma \overrightarrow{E}}}+\underset{\text{wave term}}{\underbrace{{\left(i\omega \right)}^2{\mu }_0ϵ\overrightarrow{E}}},$$

(23) $$\therefore {\nabla }^2\overrightarrow{E}=\left[i{\mu }_0\omega \sigma +{\left(i\omega \right)}^2{\mu }_0ϵ\right]\overrightarrow{E}.$$

To Eq. (20),

(24) $${\nabla }^2\overrightarrow{H}=-\sigma \left(-i{\mu }_0\omega \overrightarrow{H}\right)-i\omega ϵ\left(-i{\mu }_0\omega \overrightarrow{H}\right),$$

(25) $$\Rightarrow {\nabla }^2\overrightarrow{H}=\underset{\text{diffusion term}}{\underbrace{i{\mu }_0\omega \sigma \overrightarrow{H}}}+\underset{\text{wave term}}{\underbrace{{\left(i\omega \right)}^2ϵ{\mu }_0\overrightarrow{H}}},$$

(26) $$\therefore {\nabla }^2\overrightarrow{H}=\left[i{\mu }_0\omega \sigma +{\left(i\omega \right)}^2ϵ{\mu }_0\right]\overrightarrow{H}.$$

By reordering the equations and replacing $k^2=\left[{\omega }^2ϵ{\mu }_0-i{\mu }_0\omega \sigma \right]$ we obtain the Helmholtz equations:

(27) $${\nabla }^2\overrightarrow{E}+k^2\overrightarrow{E}=0,$$

(28) $${\nabla }^2\overrightarrow{H}+k^2\overrightarrow{H}=0.$$

Using the quasi-static approximation, where the displacement currents are negligible due their small magnitudes compared to conductivity currents $\left(\omega ϵ<<\sigma \right)$, the wave term is neglected and $k^2=i{\mu }_0\omega \sigma$. By this consideration, the diffusion equations are described by,

(29) $${\nabla }^2\overrightarrow{E}-i{\mu }_0\omega \sigma \overrightarrow{E}=0,$$

(30) $${\nabla }^2\overrightarrow{H}-i{\mu }_0\omega \sigma \overrightarrow{H}=0.$$

For a linearly polarized field $E_x$ or $H_x$ in a homogeneous medium without sources, the scalar Helmholtz equations are:

(31) $${\nabla }^2\overrightarrow{E}=\frac{{\partial }^2{E}_x}{\partial x^2}+\frac{{\partial }^2{E}_x}{\partial y^2}+\frac{{\partial }^2{E}_x}{\partial z^2},$$

(32) $${\nabla }^2\overrightarrow{H}=\frac{{\partial }^2{H}_x}{\partial x^2}+\frac{{\partial }^2{H}_x}{\partial y^2}+\frac{{\partial }^2{H}_x}{\partial z^2}.$$

6. Skin Depth

Knowing that $k^2=\left[{\omega }^2ϵ{\mu }_0-i{\mu }_0\omega \sigma \right]$ and that k is complex, it can be rewritten as $k^2={\left(\alpha +i\beta \right)}^2$, we have,

(33) $$\left[{\omega }^2ϵ{\mu }_0-i{\mu }_0\omega \sigma \right]={\left(\alpha +i\beta \right)}^2,$$

(34) $$\Rightarrow {\omega }^2ϵ{\mu }_0-i{\mu }_0\omega \sigma ={\alpha }^2+2i\alpha \beta -{\beta }^2.$$

Creating the polynomial system,

(35) $$\{{\alpha }^2-{\beta }^2={\omega }^2{\mu }_0ϵ,$$

(36) $$2\alpha \beta =-{\mu }_0\omega \sigma .$$

Squaring both equations,

(37) $$\{({\alpha }^2-{\beta }^2)=({\omega }^2{\mu }_0∍{)}^2,$$

(38) $${\left(2\alpha \beta \right)}^2={\left(-{\mu }_0\omega \sigma \right)}^2.$$

Developing Eqs. (37) and (38),

(39) $$\{{\alpha }^4-2{\alpha }^2{\beta }^2+{\beta }^4={\omega }^4{\mu }_0^2{ϵ}^2,$$

(40) $$4{\alpha }^2{\beta }^2={\mu }_0^2{\omega }^2{\sigma }^2.$$

Adding the Eqs. (39) and (40),

(41) $${\alpha }^4-2{\alpha }^2{\beta }^2+{\beta }^4+4{\alpha }^2{\beta }^2={\omega }^4{\mu }_0^2{ϵ}^2+{\mu }_0^2{\omega }^2{\sigma }^2,$$

(42) $$\Rightarrow {\alpha }^4+2{\alpha }^2{\beta }^2+{\beta }^4={\omega }^4{\mu }_0^2{ϵ}^2\left(1+\frac{{\sigma }^2}{{\omega }^2{ϵ}^2}\right),$$

(43) $$\therefore {\left({\alpha }^2+{\beta }^2\right)}^2={\omega }^4{\mu }_0^2{ϵ}^2\left(1+\frac{{\sigma }^2}{{\omega }^2{ϵ}^2}\right).$$

Taking the square root,

(44) $$\left({\alpha }^2+{\beta }^2\right)={\omega }^2{\mu }_0ϵ\sqrt{1+\frac{{\sigma }^2}{{\omega }^2{ϵ}^2}}.$$

Adding the Eq. (35) to (44),

(45) $$2{\alpha }^2={\omega }^2{\mu }_0ϵ\sqrt{1+\frac{{\sigma }^2}{{\omega }^2{ϵ}^2}}+{\omega }^2{\mu }_0ϵ,$$

(46) $$\Rightarrow {\alpha }^2=\frac{{\omega }^2{\mu }_0ϵ}{2}\left[\sqrt{1+\frac{{\sigma }^2}{{\omega }^2{ϵ}^2}}+1\right],$$

(47) $$\therefore \alpha =\omega \sqrt{\frac{{\mu }_0ϵ}{2}\left[\sqrt{1+\frac{{\sigma }^2}{{\omega }^2{ϵ}^2}}+1\right]}.$$

Subtracting the Eq. (35) to (44),

(48) $$2{\beta }^2={\omega }^2{\mu }_0ϵ\sqrt{1+\frac{{\sigma }^2}{{\omega }^2{ϵ}^2}}-{\omega }^2{\mu }_0ϵ,$$

(49) $$\Rightarrow {\beta }^2=\frac{{\omega }^2{\mu }_0ϵ}{2}\left[\sqrt{1+\frac{{\sigma }^2}{{\omega }^2{ϵ}^2}}-1\right],$$

(50) $$\therefore \beta =\omega \sqrt{\frac{{\mu }_0ϵ}{2}\left[\sqrt{1+\frac{{\sigma }^2}{{\omega }^2{ϵ}^2}}-1\right]}.$$

where $\alpha$ is the attenuation factor and $\beta$ is the phase factor. The units of $\alpha ,\beta ,k$ are respectively, $Np/m$, $rad/m$ and $m^{-1}$.

Nabighian [⁵] affirms that when conductive currents dominate over displacement currents, as is customary in electrical prospecting, $\alpha$ and $\beta$ are identical real quantities defined by

(51) $$\alpha =\beta =\sqrt{\frac{\omega {\mu }_0\sigma }{2}}.$$

The skin depth $\delta =1/\beta$ is the distance over which amplitude of EM field is reduced by a factor $1/e$, in quasi-static approximation,

(52) $$\delta =\sqrt{\frac{2}{\omega {\mu }_0\sigma }},$$

replacing ${\mu }_0$ and $\omega =2\pi f$, we have:

(53) $$\delta \approx 503\sqrt{\frac{\rho }{f}}.$$

$\delta$ has the dimension of meter, $\rho$$\left(\Omega \cdot m\right)$ is the resistivity of the medium, and $f\left(Hz\right)$ is the frequency.

7. Solution of the Wave Equations

In a uniform Earth, with homogeneous signal sources located far enough to be considered at infinity, the incident electromagnetic waves are assumed to propagate vertically, i.e., parallel to the $z-axis$. Thus, the solutions to the diffusion equation for the electric and magnetic fields can be written as:

(54) $$\overrightarrow{E}(z,t)={\overrightarrow{E}}_0{e}^{-i\left(kz-\omega t\right)}={\overrightarrow{E}}_0{e}^{-i\left[\left(\alpha +i\beta \right)z-\omega t\right]},$$

(55) $$\overrightarrow{H}(z,t)={\overrightarrow{H}}_0{e}^{-i\left(kz-\omega t\right)},$$

where ${\overrightarrow{H}}_0$ and ${\overrightarrow{E}}_0$ are initial amplitudes. Substituting $k$ into Equations (54) to (55),

(56) $$\overrightarrow{E}(z,t)={\overrightarrow{E}}_0{e}^{-i\left[\left(\alpha +i\beta \right)z-\omega t\right]},$$

(57) $$\therefore \overrightarrow{E}(z,t)={\overrightarrow{E}}_0{e}^{-i\alpha z}{e}^{-\beta z}{e}^{i\omega t},$$

(58) $$\overrightarrow{H}(z,t)={\overrightarrow{H}}_0{e}^{-i\left[\left(\alpha +i\beta \right)z-\omega t\right]},$$

(59) $$\therefore \overrightarrow{H}(z,t)={\overrightarrow{H}}_0{e}^{-i\alpha z}{e}^{-\beta z}{e}^{i\omega t}.$$

The Eqs. (57) and (59) highlights three important components of wave behavior in conductive media [1]: since $\beta$ is real, the $e^{-\beta z}$ gets smaller and $z$ gets bigger. It represents the attenuation, responsible for the exponential decay of the field with depth; the therm $e^{-i\alpha z}$ states that the wave varies sinusoidally with $z$; the term $e^{i\omega t}$ states that the wave varies sinusoidally with $t$.

8. Impedance Tensor

The electromagnetic impedance tensor $\mathbf{Z}$ is defined as the ratio between the horizontal complex components of the electric fields and the magnetic fields, measured at the same location and for a given angular frequency $\left(\omega \right)$. Mathematically, we have,

(60) $$\mathbf{Z}=\frac{\overrightarrow{E}(z,t)}{\overrightarrow{H}(z,t)}=\frac{{\overrightarrow{E}}_0{e}^{-i\alpha z}{e}^{-\beta z}{e}^{i\omega t}}{{\overrightarrow{H}}_0{e}^{-i\alpha z}{e}^{-\beta z}{e}^{i\omega t}},$$

(61) $$\therefore \mathbf{Z}=\frac{{\overrightarrow{E}}_0}{{\overrightarrow{H}}_0}.$$

However, the notation in the equations (60) and (61) is not standard since the ratio of two vectors is not strictly defined. For a valid impedance estimation, the electric and magnetic fields must be orthogonal and originate from the same plane wave.

Consequently, the impedance is defined only for specific field components. In that case, the relationship becomes the full impedance tensor:

(62) $$\left[\begin{array}{c}{E}_x\\ E_y\end{array}\right]=\left[\begin{array}{cc}{Z}_{xx}& Z_{xy}\\ Z_{yx}& Z_{yy}\end{array}\right]\left[\begin{array}{c}{H}_x\\ H_y\end{array}\right].$$

9. Cases of Impedance: TE and TM Modes

In the context of magnetotellurics, the propagation of electromagnetic plane waves in the Earth can be described by two fundamental modes, depending on the orientation of the electric and magnetic fields. The transversal electric (TE) and transversal magnetic (TM) modes, presented in Figure 5.

The TE mode occurs for the electric field $y$-component, that is, $\overrightarrow{E}=(0,E_y,0)$, while the magnetic field has components in the $x$ and $z$ directions, that is, $\overrightarrow{H}=(H_x,0,H_z)$. The TM mode arises when the electric field has the $x$ and $z$-component, that is, $\overrightarrow{E}=(E_x,0,E_z)$, while the magnetic field has components in the $y$ direction, that is, $\overrightarrow{H}=(0,H_y,0)$.

Figure 5
Two types of wave propagation in a two-dimensional Earth model. TE mode (left) and TM mode (right).

10. One-Dimensional Earth

In a 1D model, the conductivity depends only on the depth. The incidence of electromagnetic waves is assumed to be vertical (plane waves), and the fields decouple in the TE and TM modes. In the 1D model, the components of $Z$ obey the following relationships, where $Z_{xy}=-Z_{yx}$.

(63) $$\left[\begin{array}{c}{E}_x\\ E_y\end{array}\right]=\left[\begin{array}{cc}0& Z_{xy}\\ Z_{yx}& 0\end{array}\right]\left[\begin{array}{c}{H}_x\\ H_y\end{array}\right],$$

10.1. Equations to TE mode

From Eq. (12) and considering the vectors $\overrightarrow{E}=(0,E_y,0)$ and $\overrightarrow{H}=(H_x,0,H_z)$, then

(64) $$\left|\begin{array}{ccc}\widehat{i}& \widehat{j}& \widehat{k}\\ \frac{\partial }{\partial x}& \frac{\partial }{\partial y}& \frac{\partial }{\partial z}\\ 0& E_y& 0\end{array}\right|=\frac{\partial E_y}{\partial x}\widehat{k}-\frac{\partial E_y}{\partial z}\widehat{i}.$$

Then,

(65) $$\nabla \times \overrightarrow{E}=\frac{\partial E_y}{\partial x}\widehat{k}-\frac{\partial E_y}{\partial z}\widehat{i}.$$

According to Faraday’s law, Eq. (12) is described by:

(66) $$\left[\begin{array}{c}-\frac{\partial E_y}{\partial z}\\ 0\\ \frac{\partial E_y}{\partial x}\end{array}\right]=-i\omega {\mu }_0\left[\begin{array}{c}{H}_x\\ 0\\ H_z\end{array}\right],$$

and how electric field only depends from depth, $\partial E_y/\partial x=0$, so we have

(67) $$\frac{\partial E_y}{\partial z}=i\omega {\mu }_0{H}_x.$$

From Eq. (13) and considering a quasi-static condition (assumes that at low frequencies the displacement current term is negligible compared to the conduction current), where $\nabla \times \overrightarrow{H}=\sigma \overrightarrow{E}$ and how magnetic field only depends from depth, we obtain:

(68) $$\left|\begin{array}{ccc}\widehat{i}& \widehat{j}& \widehat{k}\\ \frac{\partial }{\partial x}& \frac{\partial }{\partial y}& \frac{\partial }{\partial z}\\ H_x& 0& H_z\end{array}\right|=\frac{\partial H_z}{\partial y}\widehat{i}+\frac{\partial H_x}{\partial z}\widehat{j}-\frac{\partial H_x}{\partial y}\widehat{k}-\frac{\partial H_z}{\partial x}\widehat{j}.$$

(69) $$\Rightarrow \frac{\partial H_x}{\partial z}=\sigma E_y.$$

Differentiating Eq (67) with respect to $z$, we obtain:

(70) $$\frac{{\partial }^2{E}_y}{\partial z^2}=i\omega {\mu }_0\frac{\partial H_x}{\partial z}.$$

Substituting into the expression above:

(71) $$\frac{{\partial }^2{E}_y}{\partial z^2}=i\omega {\mu }_0\sigma E_y.$$

This yields a second-order differential equation, in which the complex wave number is defined as $k^2=-i{\mu }_0\omega \sigma$. The equation becomes:

(72) $$\frac{{\partial }^2{E}_y}{\partial z^2}=-k^2{E}_y\Rightarrow \frac{{\partial }^2{E}_y}{\partial z^2}+k^2{E}_y=0.$$

This is a homogeneous linear second-order differential equation with constant coefficients. The general form of such an equation is:

(73) $$y^{\prime \prime }+a{y}^{\prime }+by=0.$$

The general solution is:

(74) $$y\left(x\right)=A{e}^{{\lambda }_{1}x}+B{e}^{{\lambda }_{2}x}.$$

Solving the characteristic equation:

(75) $${\lambda }^2+a\lambda +b=0\Rightarrow {\lambda }^2-k^2=0\Rightarrow \lambda =\pm k.$$

Thus, the solution for the electric field is:

(76) $$E_y\left(z\right)=A{e}^{kz}+B{e}^{-kz}.$$

To determine $H_x\left(z\right)$, we use the relation:

(77) $$\frac{\partial E_y}{\partial z}=i\omega {\mu }_0{H}_x\Rightarrow H_x\left(z\right)=\frac{1}{i\omega {\mu }_0}\frac{\partial E_y}{\partial z}.$$

Differentiating Eq. (76) with respect to $z$:

(78) $$\frac{\partial E_y}{\partial z}=Ak{e}^{kz}-Bk{e}^{-kz}.$$

Substituting into the expression for $H_x$:

(79) $$H_x\left(z\right)=\frac{1}{i\omega {\mu }_0}\left(Ak{e}^{kz}-Bk{e}^{-kz}\right).$$

The magnetotelluric impedance is defined as:

(80) $$Z_{yx}=\frac{E_y}{H_x}=\frac{A{e}^{kz}+B{e}^{-kz}}{{\displaystyle \frac{-k}{i\omega {\mu }_0}}\left(-A{e}^{kz}+B{e}^{-kz}\right)}.$$

Since electromagnetic waves attenuate with depth in a conductive medium, the exponentially increasing solution $e^{kz}$ must be discarded to avoid unphysical divergence as $z\to \infty$. Therefore:

(81) $$Z_{yx}=\frac{B{e}^{-kz}}{{\displaystyle \frac{-k}{i\omega {\mu }_0}}B{e}^{-kz}}=\frac{-i\omega {\mu }_0}{k}.$$

Substituting the definition of $k$, we obtain:

(82) $$Z_{yx}=\frac{-i\omega {\mu }_0}{\sqrt{-i\omega {\mu }_0\sigma }}\Rightarrow Z_{yx}=\sqrt{\frac{-i\omega {\mu }_0}{\sigma }}.$$

The apparent resistivity is then given by:

(83) $$Z_{yx}^2={\left(\frac{-i\omega {\mu }_0}{k}\right)}^2=\frac{i\omega {\mu }_0}{\sigma }\Rightarrow \frac{1}{\sigma }=\frac{Z_{yx}^2}{i\omega {\mu }_0}.$$

Thus, the expression for the apparent resistivity becomes:

(84) $${\rho }_a=\frac{1}{\omega {\mu }_0}{\left|Z_{yx}\right|}^2.$$

Finally, the phase of the impedance, which quantifies the phase shift between the electric and magnetic fields, is defined by:

(85) $$\varphi =\arctan \left(\frac{\mathrm{Im}\left(Z\right)}{\mathrm{Re}\left(Z\right)}\right)$$

10.2. Equations to TM mode

From Eq. (12) and considering the vectors $\overrightarrow{E}=(E_x,0,E_z)$ and $\overrightarrow{H}=(0,H_y,0)$ then,

(86) $$\left|\begin{array}{ccc}\widehat{i}& \widehat{j}& \widehat{k}\\ \frac{\partial }{\partial x}& \frac{\partial }{\partial y}& \frac{\partial }{\partial z}\\ E_x& 0& E_z\end{array}\right|=\frac{\partial E_x}{\partial z}\widehat{j}+\frac{\partial E_z}{\partial y}\widehat{i}-\frac{\partial E_x}{\partial y}\widehat{k}.$$

Then,

(87) $$\nabla \times \overrightarrow{E}=\frac{\partial E_x}{\partial z}\widehat{j}+\frac{\partial E_z}{\partial y}\widehat{i}-\frac{\partial E_x}{\partial y}\widehat{k},$$

how electric field only depends on depth, $\partial E_z/\partial y=\partial E_x/\partial y=0$, so we have:

(88) $$\frac{\partial E_x}{\partial z}=-i\omega {\mu }_0{H}_y.$$

From Eq. (13) and considering a quasi-static condition where $\omega ϵ\ll \sigma$,then $\nabla \times \overrightarrow{H}=\sigma \overrightarrow{E}$, we get:

(89) $$\left|\begin{array}{ccc}\widehat{i}& \widehat{j}& \widehat{k}\\ \frac{\partial }{\partial x}& \frac{\partial }{\partial y}& \frac{\partial }{\partial z}\\ 0& H_y& 0\end{array}\right|=\frac{\partial H_y}{\partial x}\widehat{k}-\frac{\partial H_y}{\partial z}\widehat{i}.$$

According to Ampere’s law, Eq. (13),

(90) $$\left[\begin{array}{c}-\frac{\partial H_y}{\partial z}\\ 0\\ \frac{\partial H_y}{\partial x}\end{array}\right]=\sigma \left[\begin{array}{c}{E}_x\\ 0\\ E_z\end{array}\right].$$

As the magnetic field only depends on depth, $\partial H_y/\partial x=0$ we have,

(91) $$\frac{\partial H_y}{\partial z}=-\sigma E_x.$$

Differentiating Eq. (88) with respect to $z$, we get

(92) $$\frac{{\partial }^2{E}_x}{\partial z^2}=-i\omega {\mu }_0\frac{\partial H_y}{\partial z}.$$

Replacing Eq. (91) in (92),

(93) $$\frac{{\partial }^2{E}_x}{\partial z^2}=i\omega {\mu }_0\sigma E_x.$$

Replacing $k^2=-i{\mu }_0\omega \sigma$ and reordering the equation,

(94) $$\frac{{\partial }^2{E}_x}{\partial z^2}=-k^2{E}_x\Rightarrow \frac{{\partial }^2{E}_x}{\partial z^2}+k^2{E}_x=0.$$

Again, we have a differential equation with constant coefficients of the second order with $\lambda =\pm k$. Then,

(95) $$E_x\left(z\right)=A\cdot e^{kz}+B\cdot e^{-kz}.$$

To find $H_y\left(z\right)$, we can use Eq. (88):

(96) $$\frac{\partial E_x}{\partial z}=-i\omega {\mu }_0{H}_y\Rightarrow H_y\left(z\right)=\frac{-1}{i\omega {\mu }_0}\frac{\partial E_x}{\partial z}.$$

Differentiating Eq. (95) with respect to $z$,

(97) $$\frac{\partial E_x}{\partial z}=A\cdot k\cdot e^{kz}-B\cdot k\cdot e^{-kz}.$$

Applying Eq. (95) in (96),

(98) $$H_y\left(z\right)=\frac{-1}{i\omega {\mu }_0}\left(A\cdot k\cdot e^{kz}-B\cdot k\cdot e^{-kz}\right).$$

The impedance is defined as:

(99) $$Z_{xy}=\frac{E_x}{H_y}=\frac{A\cdot e^{kz}+B\cdot e^{-kz}}{{\displaystyle \frac{-k}{i\omega {\mu }_0}}\left(A\cdot e^{kz}-B\cdot e^{-kz}\right)}.$$

Considering $e^{kz}$ is nullified, then,

(100) $$Z_{xy}=\frac{E_x}{H_y}=\frac{B\cdot e^{-kz}}{{\displaystyle \frac{-k}{i\omega {\mu }_0}}\left(-B\cdot k\cdot e^{-kz}\right)}=\frac{i\omega {\mu }_0}{k}.$$

Thus, we have,

(101) $$Z_{xy}=\frac{i\omega {\mu }_0}{\sqrt{i\omega {\mu }_0\sigma }}\Rightarrow Z_{xy}=\sqrt{\frac{i\omega {\mu }_0}{\sigma }}.$$

The apparent resistivity follows the same form:

(102) $${\rho }_a=\frac{1}{\omega {\mu }_0}{\left|Z_{xy}\right|}^2.$$

The phase of the impedance is identical to Eq. (85).

11. Recurrence Formulas for Impedance in 1D Earth Models

A variety of equations have been developed to solve 1D magnetotelluric models. These equations primarily compute the transfer function (impedance) at the top of the N-th layer. The solution is obtained iteratively, starting from the bottom layer, which is assumed to be a homogeneous half-space. The solution is then recursively propagated upward through the overlying layers until it reaches the surface.

This article presents three of these equations. Additionally, all of them are accompanied by Python code, which is available on GitHub, to facilitate a better understanding of the underlying physical modeling. The program can be downloaded from link: here.

11.1. Wait recursion (1954)

Wait [⁷] criticizes the validity of Cagniard [²] analysis of the behavior of telluric currents, refining the conditions under which the concept of field impedance is valid. Wait [⁷] affirms that the harmonic components of the electric field and the magnetic field tangential to the ground are only proportional to one another if the fields are sufficiently slowly (not changing abruptly from one point to another) varying over the surface of the ground. The visualization of the medium, the distribution of physical properties and the model parameters is shown on Figure 6.

Figure 6
Wait’s interpretation and notation of a layered ground.

Wait results were extended to include the effects of a layered ground with both conductivity and susceptibility variations. Finally, the corresponding transient problem for a two-layer horizontally stratified earth is solved, and a general equation for N-layers could be derived.

We have two types of impedance $\left(\eta \right)$, the intrinsic ${\eta }_j$ which is the impedance from an homogeneous medium which only depends from the properties of the layer and the entry ${\eta }^{\left(j\right)}$ that represents the impedance from the top of the ${\eta }^{\left(j\right)}$-layer. Other properties from equations are: the conductivity ${\sigma }_j$, the permittivity ${\mu }_j$, the subsurface propagation constant, here namely ${\gamma }_j$ and the angular frequency ${\omega }_j$, described by Eq. (103), from which we use frequency $f_i$ so,

(103) $${\omega }_i=2\pi f_i,{\eta }_j=\sqrt{\frac{i{\omega }_i{\mu }_0}{{\sigma }_j}}.$$

The other one is the entry impedance ${\eta }^{\left(j\right)}$

(104) $${\gamma }_j=\sqrt{i{\omega }_i{\mu }_0{\sigma }_j},$$

(105) $${\eta }^{\left(j\right)}={\eta }_j\cdot \frac{{\eta }^{\left(j+1\right)}+{\eta }_j\tanh \left({\gamma }_j{h}_j\right)}{{\eta }_j+{\eta }^{\left(j+1\right)}\tanh \left({\gamma }_j{h}_j\right)}.$$

Wait’s recursion compute the impedance at the top of the N-th layer. The equation is solved iteratively, starting from the bottom layer, which is assumed to be a homogeneous half-space, and recursively propagating the solution upward through the overlying layers until the surface.

11.2. Constable recursion (1987)

Constable [⁸] introduced an alternative recursive formulation for calculating the MT impedance, in which the quantity $c$, somewhat informally referred to as the “complex C value”, represents the layer-to-layer impedance response. This approach provides a natural and direct means of expressing MT data in terms of impedance. The recursion is defined as:

(106)$$\begin{array}{c}c=\frac{1}{k_1}\coth [k_1{t}_1+{\coth }^{-1}(\frac{k_1}{k_2}\coth \{k_2{t}_2\cdots \\ +\coth (k_{N-1}{t}_{N-1}+{\coth }^{-1}\left(\frac{k_{N-1}}{k_N}\right)\cdots \left\}\right],\end{array}$$

where $k_j=\sqrt{{\displaystyle \frac{i\omega {\mu }_0}{{\rho }_j}}}$ is the complex wave number for layer $j$, with ${\rho }_j$ denoting the electrical resistivity of layer $j$ (in $\Omega \cdot \text{m}$) and $t_j$ representing its thickness (in meters).

11.3. Grandis recursion (1999)

Grandis [⁹] proposed an alternative algorithm for one-dimensional magnetotelluric response calculation. Instead of using the classical recursive formula that directly relates the impedance at the surface of two successive layers, the proposed method employs a recursive formulation based on the electromagnetic fields at the top of two consecutive layers. This approach leads to a matrix multiplication scheme that offers improved numerical stability and flexibility in certain computational scenarios.

The surface impedance $Z_1$, representing the magnetotelluric response at the Earth’s surface, is calculated from the elements of the total transfer matrix $S$, which encapsulates the cumulative effect of all subsurface layers. This total transfer matrix $S$ is obtained by sequentially multiplying the individual transfer matrices $T_j$ associated with each of the $N-1$ finite layers in the model:

(107) $$S=T_1\times T_2\times \cdots \times T_{N-1},$$

where each $T_j$ relates the electromagnetic fields across the $j$-th layer. The final expression for the surface impe-dance is given by:

(108) $$Z_1=\frac{s_{11}{Z}_{I,N}+s_{12}}{s_{21}{Z}_{I,N}+s_{22}},$$

where $s_{ij}$ are the elements of the total transfer matrix $S$, and $Z_{I,N}$ is the intrinsic impedance of the $N$-th layer, modeled as a homogeneous half-space extending to infinity.

The original implementation of the algorithm was developed in Fortran 77. For the purpose of study and analysis, it was translated into Python 3.8.

12. Sintetic Application

To demonstrate the application of the modeling codes and analyze their responses, we applied them to a synthetic model illustrated in Figure 7. The corresponding results are presented in Figure 8. All modeling approaches yielded identical response curves, reinforcing the consistency and reliability of the implemented methods.

Figure 7
Synthetic 1D model of physical properties (resistivity) used in magnetotelluric simulations.

Figure 8
Apparent resistivity and phase curves obtained from the synthetic model shown in Figure 9.

13. Conclusions

The basic theory of the magnetotelluric (MT) method was presented, covering fundamental concepts such as signal sources, frequency range, and boundary conditions. The development of the method’s core equations began with Maxwell’s curl equations in the frequency domain and proceeded step by step, incorporating key concepts such as constitutive relations, diffusion equations, skin depth, and the specific characteristics of the impedance tensor.

The one-dimensional (1D) case we also discussed leading to analytical expressions for the impedance, where $Z_{xy}=-Z_{yx}$, as well as for the apparent resistivity and phase shift. Furthermore, three classical formulations for 1D modeling were presented: the recursive approaches of Wait, Constable, and Grandis. In support of teaching and learning, these methods were implemented in Python, providing practical tools to help students and educators bridge the gap between theoretical derivations and computational applications.

Data Availability

All the data supporting the results of this study were published in the article itself.

Referências

• ^[1] A.N. Tikhonov, Doklady 73, 295 (1950).
• ^[2] L. Cagniard, Geophysics 18, 605 (1953).
• ^[3] S. Constable, Geophysics 75, 75A67 (2010).
• ^[4] K. Vozoff, Geophysics 37, 98 (1972).
• ^[5] M.N. Nabighian (ed.), Electromagnetic methods in applied geophysics: Volume 1, Theory (Society of Exploration Geophysicists, Tulsa, 1988).
• ^[6] D.J. Griffiths, Eletrodinâmica (Pearson Addison-Wesley, São Paulo, 2011).
• ^[7] J.R. Wait, Geophysics 19, 281 (1954).
• ^[8] S.C. Constable, R.L. Parker and C.G. Constable, Geophysics 52, 289 (1987).
• ^[9] H. Grandis, Computers & Geosciences 25, 119 (1999).

Edited by

• Editor-in-Chief:
  Marcello Ferreirahttps://orcid.org/0000-0003-4945-3169

Publication Dates

• Publication in this collection
  09 Jan 2026
• Date of issue
  2025

History

• Received
  12 July 2025
• Reviewed
  28 Sept 2025
• Accepted
  28 Oct 2025