COMSOL 6.4 - Magnetotellurics

Magnetotellurics (MT) is a method for estimating the electric resistivity (or the reciprocal electric conductivity) profile of the Earth’s crust using the natural electromagnetic source provided by the ionosphere.

The resistivity of the Earth’s crust may vary over several orders of magnitude and gives a strong indication as to what kind of rocks are present. Metamorphic, igneous, or sedimentary rocks might have resistivities varying from 0.1 to 105 Ωm, depending on several petrophysical parameters.

The electromagnetic source in the MT method has its origin in the variations in the magnetic field of the Earth, caused by solar wind and electromagnetic noise. These variations induce electric currents in the Earth’s crust, proportional to the electric conductivity.

The electromagnetic source is random by nature, but nevertheless, the low frequency signals and the shape of the ionosphere generate plane-wave components that can be correlated over large areas. By simultaneous measurements of the natural source electromagnetic field as a function of time at different locations in an area of interest, one can statistically extract the local electromagnetic impedance as a function of frequency. This impedance can then be used to estimate the electric resistivity (or conductivity) as a function of depth.

This application is inspired by a study published in 1997 (Ref. 1). In that article, various scientific groups compared the software performance of a few models. This model, called COMMEMI-3D-2, has become one of the benchmarks for MT modeling.

Impedance Tensor

MT analysis is based on the impedance tensor Z defined as

here, H is the magnetic field and E denotes the electric field.

Normally, only the horizontal components of the fields are analyzed since the incident field is a plane-wave parallel to the Earth’s surface. When the z direction is the vertical axis, the relation for the horizontal components reads

The impedance tensor components are then analyzed to determine the electric properties of the subsurface. For example, if the subsurface can be approximated with a 1D model (resistivity variations as a function of depth), the following relations follow:

If the subsurface can be approximated with a 2D structure and the electric field is parallel to the strike (the axis in the horizontal plane along which the resistivity is constant), the mode is called transverse electric (TE). If the magnetic field is oriented along the strike, the mode is called transverse magnetic (TM). These two modes are uncoupled. Assuming that the x-axis is oriented along the strike, the impedance matrix component Zxy corresponds to the TE mode and Zyx describes the TM mode. In a 2D approximation, the impedance tensor’s components are

and

In the general 3D case, the diagonal components Zxx and Zyy are different and nonzero. These diagonal matrix elements are often analyzed to determine the dimensionality of the subsurface.

Traditionally, MT surveys are performed using electromagnetic sensors that are placed on the sea bottom or on land forming either lines or grids. The alignment of the sensors is often such that the line is perpendicular to the expected strike of the formation, such as a fault line.

In a right-handed coordinate system with the z-axis pointing downward from the surface of the model, the TE and TM modes relative to the coordinate system are then related to the impedance tensor, meaning that ZTE = Zxy and ZTM = Zyx. By positioning the sensors along a line perpendicular to the strike, the sensitivity is increased by measuring the biggest differences in induced currents across and parallel to the strike.

Apparent Resistivity

When modeling MT with a known source, and the applied magnetic field is aligned along the y-axis, the apparent resistivity components are calculated from

and the components in the impedance tensor are

Also, when the applied magnetic field is aligned along the x-axis, the apparent resistivity components are calculated from

and the components of the impedance tensor are

Therefore, in order to determine all the components of the impedance tensor Z, two simulations must be run with two different field polarizations.

Model Definition

The geometry and parameters for this application are taken from Ref. 1. Two rectangular inserts of high conductivity are placed in the upper layer of a three-layer earth. The main idea of the study is to estimate those conductivities by the MT method.

Figure 1: A 70 by 70 km piece of Earth’s crust, consisting of three layers of different conductivities. From top to bottom, the conductivities are 10, 100, and 0.1 Ωm. The upper layer has two rectangular inserts of conductivities 1 and 100 Ωm.

The application uses the Magnetic Fields interface in a 3D geometry. The magnetic field source is a plane wave that induces horizontal currents in the whole model. You impose suitable symmetry boundary conditions on the lateral sides to simulate an infinite domain. This assumption holds well if the domain is large enough around the area of interest. Make the computational domain larger if boundary effects are present. The imposed magnetic field is parallel or perpendicular to these boundaries depending on the chosen polarization.

The magnetotellurics analysis is based on the decomposition of the incoming plane waves into two waves with perpendicular polarizations. Hence, you solve these two polarizations in two independent physics interfaces. Then you solve the two physics interfaces in a single solver sequence to get all results in a single dataset. The solutions from the two polarizations can thus be obtained by a parametric sweep over several frequencies.

Results and Discussion

The components of the apparent resistivity can now be extracted from the top layer of the model and plotted as surface plots for each frequency. This is where the measurements are made when real data is collected in magnetotelluric surveys. As discussed earlier, the apparent resistivity at a certain position is equal to the resistivity below it in the case of a uniform half-space (1D model).

Because the skin depth is small at high frequencies, the values of the apparent resistivity should be equal to the material resistivity in areas where a half-space approximation is valid. This can be seen in Figure 2 and Figure 3. In the areas with uniform resistivities, far from the “fault lines,” the apparent resistivity is almost equal to the resistivity of the material right underneath.

For a magnetotelluric survey, the electromagnetic sensors are usually deployed along a line perpendicular to the expected fault line. Then, the apparent resistivity at a given frequency is plotted as a function of position along the line. See Figure 4 where data is extracted along a line crossing the three fault lines in the model. The results in this plot are comparable with the results published in the reference article (Ref. 1).

At lower frequencies, the discrepancy between the apparent resistivity and the material resistivity values increases.

Some magnetotellurics surveys also measure the z-component of the magnetic field. The plot of Hz is called the tipper plot. A rule of thumb is that a nonzero value of the tipper indicates a large change in resistivity from one side of a fault line to another. When moving along the direction of the field vector at a given phase value, the sign of the tipper indicates if one is moving from a region of higher resistivity to a low resistivity region (tipper is positive) or the opposite (tipper is negative). See Figure 6 where the tipper has been plotted on the domain, for the frequency of 0.01 Hz.

Figure 2: Logarithm of apparent resistivity ρxy, top view of the 3D domain.

Figure 3: Logarithm of apparent resistivity ρyx, top view of the 3D domain.

Figure 4: Apparent resistivity across strike (ρxy and ρyx), for the two frequencies 0.1 and 0.01 Hz.

Figure 5: Skin depth (m) at the center of the 3D model for the frequency of 0.01 Hz.

Figure 6: Tipper plot, Hz component of magnetic field (A/m) for the frequency of 0.01 Hz.

It is worth noting that magnetotelluric data analysis is typically done at frequencies ranging from 0.1 mHz to 10 Hz in an evenly spaced logarithmic scale. In this example, only two frequencies are modeled for simplicity. To calculate the response at other frequencies, the model should be adjusted by adapting the model size and the mesh accordingly. This is left as an exercise to the user, but here are a few things to consider. For very high frequencies, the deep part of the model can be ignored. For low frequencies, the lateral dimensions should be larger. Also, evaluate the skin depth to be sure that the model is deep enough for the very low frequencies.

Reference

1. M. Zhdanov, I.M. Varentsov, J.T. Weaver, N.G. Golubev, and V.A. Krylov, “Methods for Modelling Electromagnetic Fields. Results from COMMEMI — The International Project on the Comparison of Modelling Methods for Electromagnetic Induction,” J. Appl. Geophys., vol. 37, pp. 133–271, 1997.

ACDC_Module/Devices,_Resistive/magnetotellurics

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

The parameters Lx and Ly can be used to control the size of the domain to be studied. For very low frequencies, these parameters may have to take values up to 100 km. For high frequencies, the domain can be much smaller.

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
Lx	70[km]	70000 m	Domain size in x direction
Ly	70[km]	70000 m	Domain size in y direction
Lh	20[km]	20000 m	Height of bottom layers
h_box	10[km]	10000 m	Box height
w_box	20[km]	20000 m	Box width
d_box	40[km]	40000 m	Box depth

Geometry 1

Block 1 (blk1)

In the toolbar, click

In the window for , locate the section.

In the text field, type Lx.

In the text field, type Ly.

In the text field, type Lh.

Locate the section. From the list, choose .

In the text field, type -2*Lh.

Right-click and choose .

Block 2 (blk2)

In the toolbar, click

In the window for , locate the section.

In the text field, type Lx.

In the text field, type Ly.

In the text field, type Lh.

Locate the section. From the list, choose .

In the text field, type -Lh.

Right-click and choose .

Click the

button in the toolbar.

Block 3 (blk3)

In the toolbar, click

In the window for , locate the section.

In the text field, type Lx.

In the text field, type Ly.

In the text field, type h_box.

Locate the section. From the list, choose .

In the text field, type -h_box/2.

Right-click and choose .

Block 4 (blk4)

In the toolbar, click

In the window for , locate the section.

In the text field, type w_box.

In the text field, type d_box.

In the text field, type h_box.

Locate the section. From the list, choose .

In the text field, type -w_box/2.

In the text field, type -h_box/2.

Right-click and choose .

Block 5 (blk5)

In the toolbar, click

In the window for , locate the section.

In the text field, type w_box.

In the text field, type d_box.

In the text field, type h_box.

Locate the section. From the list, choose .

In the text field, type w_box/2.

In the text field, type -h_box/2.

Right-click and choose .

Click the

button in the toolbar.

Form Union (fin)

In the window, right-click and choose .

Import the definition of the components of the resistivity tensor from a text file.

Definitions

Variables 1

In the toolbar, click

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file magnetotellurics_variables.txt.

x Boundaries

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundaries 1, 4, 7, and 25–27 only.

In the text field, type x Boundaries.

y Boundaries

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundaries 2, 5, 8, and 11–13 only.

In the text field, type y Boundaries.

Top

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundaries 10, 17, and 22 only.

In the text field, type Top.

Magnetic Fields (mf)

Ampère’s Law in Solids 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

The magnetic field source is a plane wave that induces horizontal currents in the whole model. Use a combination of the default Magnetic Insulation boundary condition and the Symmetry Plane boundary condition and to make the domain look like an infinite domain. This assumption should hold well if the domain is large enough around the area of interest with 3D effects. Make Lx and Ly larger if the boundary effects are present. The imposed magnetic field should be parallel to the Magnetic Insulation boundaries and perpendicular to the antisymmetry boundaries.

The magnetotelluric analysis is based on decomposition of the incoming plane waves into two waves with perpendicular polarizations. Solve these two polarizations in two independent physics interfaces. Then solve the two physics interfaces in a single solver sequence to get all results in a single dataset. The solutions from the two polarizations can thus be obtained for several frequencies.

For the other magnetic field polarization, the boundary conditions should be swapped.

Symmetry Plane 1

In the toolbar, click

and choose .

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

The magnetic field is a plane wave with arbitrary amplitude and polarization along the x-axis or y-axis. For numerical stability, it makes sense to impose a magnetic field with a large amplitude. Because the magnetotelluric analysis is based on normalization between the various field components, the amplitude is not significant.

Tangential Magnetic Field 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Specify the H0 vector as


0	x
1000	y
0	z

Magnetic Fields 2 (mf2)

In the window, under click .

Ampère’s Law in Solids 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Symmetry Plane 1

In the toolbar, click

and choose .

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Tangential Magnetic Field 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Specify the H0 vector as


1000	x
0	y
0	z

Materials

Now define the four kinds of rocks in this model. The relative permeability and permittivity are set to one. The resistivities, entered as conductivities, are 100, 10, 1 and 0.1 Ωm.

Rock 100ohmm

In the window, under right-click and choose .

Select Domains 2 and 5 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	0.01	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic

In the text field, type Rock 100ohmm.

Rock 10ohmm

Right-click and choose .

Select Domain 3 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	0.1	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic

In the text field, type Rock 10ohmm.

Rock 1ohmm

Right-click and choose .

Select Domain 4 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	1	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic

In the text field, type Rock 1ohmm.

Rock 0.1ohmm

Right-click and choose .

Select Domain 1 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	10	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic