Response of a layered Earth to a plane wave

Purpose

Plane wave propagation through the earth is at the core of the Magnetotelluric problem. Here, we investigate the physics of plane wave propagation through a 1D layered earth and discuss the connection between the physical responses and what we observe in MT data

– +

Once Loop Reflect

Introduction

We present here a 1D modelisation of the Magnetotelluric waves in a layered Earth and the associated geophysical data. Our goal is to build better representation of the different physical phenomenons and better understanding of the resulting data. This work follows the derivation presented in [WH88] and is supported by interactive apps developed in a binder.

The magnetotelluric (MT) method is a widely used geophysical technique, in particular for imaging geothermal systems, that is sensitive to Earth structures as shallow as tens of meters to depths of hundreds of kilometers. It is a passive method that use plane waves generated mostly in the Earth’s Atmosphere. High frequency waves are mainly produced by lightning strikes all around the globe, traveling through the Earth’s Ionosphere that acts as a waveguide. Low frequency waves are produced through the interaction of the Earth’s Ionosphere with solar wind and Earth’s magnetic field.

In Magnetotelluric problems, the key diagnostic physical property is electrical conductivity \(\sigma\). In most cases we expect the contrasts of the others physical properties (magnetic permeability, dielectric permittivity) to be negligible compared to the electrical conductivity contrasts. As a result, at sufficiently low frequencies ( \(\prec 10^5 Hz\) ), the impact of the other physical properties contrasts on the EM response is expected to be negligible.

Setup

Fig. 127 1D layered Earth Model

In the movie above, we show an example of plane wave electromagnetic fields propagating in a 2-layered earth with:

• a layered Earth, each layer with its own physical properties \(\sigma_j, \varepsilon_j, \mu_j\)

• a plane wave traveling along the axis \(\mathbf{\hat{z}}\) coming from air, composed of an electric field \(\mathbf{E_x}\) and a orthogonal magnetic field \(\mathbf{H_y}\). These are the fields we will be measuring with our geophysical instruments at the surface to obtain information from the underground

We can see that several phenomenons are occuring. Just to mention few of them:

• the incoming wave (down component) is reflected at the surface

• once in the ground we observe a diffusive effect of the Earth on the wave

• the depth of investigation is regulate by the damping effect on the amplitudes, which is characterized by the skin depth. We notice the decay is more important in the second layer with a higher conductivity

• \(\mathbf{E_x}\) and \(\mathbf{H_y}\) are continuous. There is a phase shift between the two.

Governing Equations

The governing equation for Magnetotelluric problem can be obtained from Maxwell’s equations. We start with Faraday’s law and Ampere’s law:

(410)\[\nabla \times \mathbf{E_x} = - i \omega \mu \mathbf{H_y}\]

(411)\[\nabla \times \mathbf{H_y} = (\sigma + i \omega \varepsilon) \mathbf{E_x}\]

Knowing that that \(\mathbf{E}\) and \(\mathbf{H}\) are divergence free, according to Gauss’s Law for Electric Fields and Gauss’s Law for Magnetic Fields, we can combine the equations to write the Helmhotz (wave propagation) equation for both \(\mathbf{E}\) and \(\mathbf{H}\) fields:

(412)\[\nabla ^2 \mathbf{E_x} + k^2 \mathbf{E_x} = 0\]

(413)\[\nabla ^2 \mathbf{H_y} + k^2 \mathbf{H_y} = 0\]

with k the wavenumber:

(414)\[k = \sqrt{\omega ^2 \mu \varepsilon - i \omega \mu \sigma }\]

In the ground, we can generally assume that the displacement current is negligible, which means \(\sigma \gg \omega \varepsilon\). In this case

(415)\[k_{ground} \simeq (1-i) \sqrt{ \frac{\omega \mu \sigma}{2} }\]

In the air, the conductivity is almost 0

(416)\[k_{air} \simeq \omega \sqrt{ \mu_0 \varepsilon_0}\]

Taking the problem from the point of view of the electric field, we know the equation (412) has a solution in the form of:

(417)\[E_x (z) = U e^{i k z} + D e^{-i k z}\]

(418)\[H_y (z) = \frac{1}{- i \omega \mu} (\nabla \times \mathbf{E_x})_y = \frac{k}{ \omega \mu} (D e^{-i k z} - U e^{i k z}) = \frac{1}{Z} (D e^{-i k z} - U e^{i k z})\]

with \(\mathbf{E_x} = E_x \mathbf{\hat{x}}\) and U and D are the complex amplitudes of the Up and Down components of the field and \(Z = \frac{ \omega \mu}{k}\) the intrinsic impedance of the space.

Writing the solution for the j-th layer (See Fig. 127), we obtain:

(419)\[E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})}\]

(420)\[H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})})\]

which can be re-written in matrix form as:

(421)\[\begin{split}\left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) = \left(\begin{matrix} 1 & 1 \\ -\frac{1}{Z_j} & \frac{1}{Z_j} \end{matrix} \right) \left(\begin{matrix} U_j \\ D_j \end{matrix} \right) = P_j \left(\begin{matrix} U_j \\ D_j \end{matrix} \right)\end{split}\]

The transfert of the Up and Down components inside a layer can then be write as such

Fig. 128 Transfert of Up and Down components inside a layer, variables definition.

\[\begin{split}\left(\begin{matrix} U_j' \\ D_j' \end{matrix} \right) = \left(\begin{matrix} e^{i k h_j} & 0 \\ 0 & e^{-i k h_j} \end{matrix} \right) \left(\begin{matrix} U_j \\ D_j \end{matrix} \right) = T_j \left(\begin{matrix} U_j \\ D_j \end{matrix} \right)\end{split}\]

With the variables U, D, U’ and D’ defined as in (Fig. 128)

Using the continuity of the tangential \(\mathbf{E_x}\) and \(\mathbf{H_y}\) field at the interfaces, we find an iterative relation between the fields in consecutive layers:

\[\begin{split}\left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) = P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right)\end{split}\]

We are now only missing a Boundary Condition to be able to compute our MT forward modeling. A reasonable one is to set the Down Amplitude to 1 and the Up Amplitude to 0 in the last layer, as there is no reflection from an other interface below.

\[\begin{split}\left(\begin{matrix} U_n \\ D_n \end{matrix} \right) = \left(\begin{matrix} 0 \\ 1 \end{matrix} \right)\end{split}\]

We assume with this boundary condition that the last layer is a half-space. Knowing all the model parameters, the forward can now be solved by first use the matrix \(P_{n}\) to calculate the fields \(\mathbf{E_{x,n}}\) and \(\mathbf{H_{y,n}}\) and then propagate the field iteratively up to the top layer using the matrix \(P_j T_j P^{-1}_J\).

Building Intuition for MT problems

Skin Depth and Depth of investigation

Fig. 129 Depth of investigation in MT

Take the amplitude of the incident component of the electric wave, \(E_{x} (z) = D e^{Im(k) z}\).

The skin depth \(\delta\) is defined as the depth where the signal has decayed to a factor \(\frac{1}{e}(\simeq\) 36%).

\[e^{-i Im(k) \delta} = \frac{1}{e}\]

(422)\[\delta = \sqrt{ \frac{2}{\omega \mu \sigma}} \simeq \frac{500}{\sqrt{\sigma f}}\]

We see the skin depth is highly dependent on both the frequency of our signal and the conductivity of the Earth material. In air , the conductivity is almost 0, so we do not notice important decreased of the electromagnetic wave. In the ground, this is different.

• the more conductive, the faster the decay is. MT can see very deep in resistive environment.

• The lower the frequency, the slower the decay is. Lowest frequencies sample the deepest structures while high frequencies bring information on shallower structures.

In Fig. 129 and in the movie, we can see that even at very high frequency (20000 Hz), MT is still a deep exploration method in resistive environment (\(10^{-5} S/m\)) with a skin depth of about 1125m. Skin Depth is often use as an estimator for the depth of investigation of a survey.

Reflection and Transmission Coefficients

Fig. 130 Reflection at interface

Fig. 131 Notations for reflection

Let define at the j-th interface \((E^i ; H^i)\) as the incident waves, \((E^r ; H^r)\) as the refleted wave and \((E^t ; H^t)\) as the transmitted wave into the ground (Fig. 131)

Using the interface conditions for the tangential components of the electric, we can write:

(423)\[E^i + E^r = E^t\]

And same for the magnetic fields

(424)\[H^i + H^r = H^t\]

Using Faraday’s law, assuming variations in \(\mu\) are negligible, we also obtain from equation (424) :

(425)\[k_j E^i - k_j E^r = k_{j+1} E^t\]

Replacing the differents components of equation (425) with equation (423), we obtain the reflection coefficient R and the transmission coefficient T:

(426)\[R = \frac{E^r}{E^i} = \frac{k_j - k_{j+1}}{k_j + k_{j+1}}\]

(427)\[T = \frac{E^t}{E^i} = \frac{2 k_j}{k_j + k_{j+1}}\]

These coefficients tell us how much energy of the incoming has been reflected or transmitted.

Refraction angle

Fig. 132 Refraction and Reflection angles

In reality, the incident wave is coming from all the possible directions in the air. So how valid is our assumption of an incident vertical wave?

What is important is the refraction angle at the Air-Earth interface, the angle of the transmitted wave in the ground.

As any wave, electromagnetic waves follow Snell’s law, that we can derive from the Frequency Domain Maxwell’s equation.

Starting from an non orthogonal incident wave, modifyig the solution for (412) , we now get:

\[E^i(x,z) =||E^i|| e^{-i k_{iz} z} e^{-i k_{ix} x}\]

\[E^r(x,z) = ||E^r|| e^{i k_{rz} z} e^{-i k_{rx} x}\]

\[E^t(x,z) = ||E^t|| e^{-i k_{tz} z} e^{-i k_{tx} x}\]

The equation (423) is still valid, for all x. This is possible if and only if \(||k_{ix}||=||k_{rx}||=||k_{tx}||\)

For the reflected wave

\[||k_{ix}||=||k_{rx}||\]

\[||k_{air}|| *sin (\theta_i) = ||k_{air}||*sin (\theta_r)\]

\[\theta_i = \theta_r\]

We find the intuitive result that the wave is reflected at the same angle than the incident wave

For the transmitted wave

\[||k_{ix}||=||k_{tx}||\]

\[||k_{air}|| *sin \theta_i = ||k_{earth}||*sin \theta_t\]

\[\theta_t = sin^{-1} (\frac{||k_{air}||}{||k_{earth}||} *sin (\theta_i)) \simeq 0\]

As \(\frac{||k_{air}||}{||k_{earth}||}\) is a really small number as the conductivity of the earth is usually several order of magnitude higher than the one of the air, \(\theta_t \simeq 0\). Any wave that hits the Earth gets refracted vertically because of the extreme contrast in conductivity, regardless of the angle of incidence.

Field Acquisition

In MT, the source is unknown but we are avoiding the problem by measuring the ratio of the fields, which cancel the amplitude of the source. The data are acquired usually at the surface. We define an apparent impedance:

(428)\[\hat{Z}_{xy} = \frac{E_x}{H_y}\]

Notice this is a complex number, with a norm and an angle.

Impendance matrix

We saw that in 1D, the horizontal orthogonal components of the electric and magnetic fields \(\mathbf{E_x}\) and \(\mathbf{H_y}\) are linked through the Faraday’s law and Ampere’s law. We can then write the same types of relationship for \(\mathbf{E_y}\) and \(\mathbf{H_x}\) and write the system in a matrix form:

\[\begin{split}\left(\begin{matrix} E_{x} \\ E_{y} \end{matrix} \right) = \left(\begin{matrix} 0 & \hat{Z}_{xy} \\ -\hat{Z}_{xy} & 0 \end{matrix} \right) \left(\begin{matrix} H_x \\ H_y \end{matrix} \right)\end{split}\]

which can be generalised:

\[\begin{split}\left(\begin{matrix} E_{x} \\ E_{y} \end{matrix} \right) = \left(\begin{matrix} \hat{Z}_{xx} & \hat{Z}_{xy} \\ \hat{Z}_{yx} & \hat{Z}_{yy} \end{matrix} \right) \left(\begin{matrix} H_x \\ H_y \end{matrix} \right)\end{split}\]

The matrix linking the component of \(\mathbf{E}\) and \(\mathbf{H}\) is called the impedance matrix.

On field, we do not know a priori the orientation of the source wave. This orientation can also changes over times if the source wave is polarised. We usally record both horizontals components of each field. If the Earth is purely 1D, a simple rotation of the matrix would allow to find the antisymetric matrix and thus obtain the apparent impedance \(\hat{Z}_{xy}\).

Note: for a pure 2D Earth, the impedance matrix is also purely off-diagonal (with the right rotation if needed) but is not anymore antisymetric. In 3D the impedance matrix is a full matrix.

Data

Fig. 133 MT data for a 2 layers Earth

Apparent Resistivity

The apparent resistivity is obtained through the amplitude of the apparent Impedance \(\hat{Z_{xy}}\).

(429)\[\rho_{app} = \frac{1}{\mu_0 \omega} |\hat{Z_{xy}}|^2\]

For a half-space, \(\rho_{app} = \rho_{earth}\) :

\[\hat{Z}_{xy} = \frac{\omega \mu}{k_{earth}} = (1+i) \sqrt{\frac{\omega \mu}{2 \sigma_{earth}}}\]

\[\rho_{app} = \frac{1}{\mu_0 \omega} |1+i|^2 \frac{\omega \mu}{2 \sigma_{earth}} = \rho_{earth}\]

For a nonhomogeneous earth, \(\rho_{app}\) at a particular frequency is an average of the conductivity of the earth on about a sphere with a radius equal to the skin depth.

Phase

The phase is obtained through the angle of the apparent Impedance \(\hat{Z}_{xy}\).

(430)\[\Theta =tan^{-1} \frac{Im(\hat{Z}_{xy})}{Re(\hat{Z}_{xy})}\]

for a half-space,

\[\Theta = tan^{-1} \frac{Im({Z_{xy}})}{Re({Z_{xy}})} = tan^{-1} 1 = \frac{\pi}{4}\]

If \(\sigma\) increases at depth, then \(\Theta\) increases before returning to 45°

If \(\sigma\) decreases at depth, then \(\Theta\) decreases before returning to 45°

Survey Design

Interpretation

Pratical Consideration