# The Finite Element Numerical Modelling of 3D Magnetotelluric

Authors: Ligang Cao; Yun Liu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The ideal numerical simulation of 3D magnetotelluric was restricted by the methodology complexity and the time-consuming calculation. Boundary values, the variation of weighted residual equation, and the hexahedral mesh generation method of finite element are three major causes. A finite element method for 3D magnetotelluric numerical modeling is presented in this paper as a solution for the problem mentioned above. In this algorithm, a hexahedral element coefficient matrix for magnetoelluric finite method is developed, which solves large-scale equations using preconditioned conjugate gradient of the first-type boundary conditions. This algorithm is verified using the homogeneous model, and the positive landform model, as well as the low resistance anomaly model.


---

## Body

## 1. Introduction

Both 1D and 2D MT numerical simulation technologies have brought benefits to economy with their maturity. However, they are both the approximate processing of a 3D numerical simulation, which is the closest to the real geological body model [1–9].With the development in computer science and computing technology, numerical simulation of 3D MT data has become the research focus. Shizhe [10] proposed the conjecture and theoretical basis of “Hexahedral mesh generation method for finite element” and “Boundary element method in 3D numerical simulation” [11]. However, there are only a few publications regarding successful scientific research on 3D numerical simulations [12–15]. The major difficulties lie in three aspects: (
1
) insufficient computing capability due to limited CPU memory; (
2
) immature parallel algorithm on multicore CPU’s; (
3
) limited boundary element 3D MT numerical simulation on pure terrain uniform earth. Shu and Chen [16] proposed a 3D forward method in horizontal electric dipole under the magnetotelluric response by using a vector finite element method. After combining scalar finite elements and vector finite elements together, Mitsuhata and Uchida [17] conducted forward computation and error analysis for the two 3D COMMEMI models using T
_
Ω formula, which verified the correctness of the program as well. Under the conditions of tetrahedral elements, Linping and Shikun [18] derived the analytical expression for the underground varying conductivity to calculate the three-dimensional magnetotelluric finite element equation. Nam et al. [19] proposed the 3D magnetotelluric response simulation based on hexahedral-vector finite element method and simulated the 3D terrain using isoparametric element method. As for the underground complex geological body, only a few papers are available from public sources [20–27].This paper deduces the hexahedral element coefficient matrix of the magnetotelluric finite element method based on Xu Shizhe’s work. It uses preconditioned conjugate gradient method to solve large-scale equations under the first-type boundary conditions. This algorithm can significantly improve the computation efficiency. Via computing analysis of homogeneous model, undulating terrain model, and 3D model, the 3D low resistivity abnormal morphology can be clearly reflected. As the frequency decreases, the errors between the calculation results and analytical solutions are found gradually decreasing.


## 2. The Boundary Problem of 3D MT

A 3D geology space model is shown in Figure1. The space is divided into two layers, air and underground, by an undulating (A
′
B
′
C
′
D
′). The surface of ABCD is the infinite boundary of upper air which is about 100 km height, and the surface of EFGH is infinite boundary of underground which is about 30 km beneath ground. In addition, there are four infinite vertical boundary surfaces: ABFE, BFGC, DCGH and ADHE, all of which are about 30 km away from the center area of measurement.Figure 1
3D space of the study.The source of natural alternating electromagnetic field is located on the surface of ABCD. When the electromagnetic wave passes through the atmosphere and then enters the underground, the whole 3D geology space becomes the study field of MT method due to continuous electric field on the tangential direction and magnetic fields on the normal direction.


### 2.1. The Boundary Value Problem of MT Method

The electromagnetic waves from the magnetotelluric field source vertically incident into the underground as harmonic plane waves with a time harmonic factore
-
i
ω
t. Maxwell coupled magnetotelluric field equations in frequency domain as
(1)
∇
×
E
=
i
ω
μ
H
,
∇
×
H
=
(
σ
-
i
ω
ε
)
E
,
where E is the electric vector, H is the magnetic vector, ω is angular frequency, μ is the medium permeability when the permeability of air is 1.2567
×
1
0
-
6
H
/
m, σ is the medium conductivity or the reciprocal of medium resistivity, ε is the dielectric medium constant, and the dielectric constant of air is 8.85
×
1
0
-
12 F/m.Curling both sides of the electric equations (1), the differential electric equation can be obtained as follows:
(2)
∇
×
(
∇
×
E
)
-
k
2
E
=
0
,
∇
×
(
∇
×
H
)
-
k
2
H
=
0
,
where
(3)
k
2
=
i
ω
μ
σ
+
ω
2
ε
μ
.The equations of electric and magnetic have the same form that can be written in one unification equation. SubstitutingE and H for U yields
(4)
∇
×
(
∇
×
U
)
-
k
2
U
=
0
.The vectorU can be rewritten to the vector form with three components:
(5)
U
=
U
x
e
x
+
U
y
e
y
+
U
z
e
z
,
where, e
x, e
y, and e
z are the unit vectors in the directions of x, y, and z.


#### 2.1.1. When It Is Working in TE Model

(a)
In the surface of ABCD shown in Figure1, the field source is the electric vector E
x. As the ratio of electric and magnetic, U
x has no impact to the final result. Therefore, U
x can be either zero or nonzero value. Here U
x
=
1. U
y and U
z are all zero as well since there are no components of U in x and y directions.
(b)
In the four vertical boundaries shown in Figure1, the propagation direction of EM wave in magnetic field is vertical downward, being perpendicular to normal direction of the boundary, which is E
×
H
⊥
Γ. Therefore, ∂
U
/
∂
n equals zero.
(c)
In the surface of EFGH in Figure1, the underground can be assumed as a homogenous medium with a consideration that this surface is beyond the study area. The electromagnetic field attenuates in a normal way, which follows the index law during downward propagation. The vector U has no value in directions y and z. 													(6)
U
x
=
c
e
-
χ
z
,
U
y
=
0
,
U
z
=
0
,where c is a constant factor with any nonzero value, χ
=
-
i
ω
μ
σ, and σ is the homogenous conductivity below the infinite surface EFGH. Derivative of (6) is obtained as follows:
(7)
∂
U
x
∂
n
=
∂
U
x
∂
z
=
-
χ
U
x
.


#### 2.1.2. When It Is Working in TM Model

The field source is one component of the magnetic vectorH
x. When magnetic field passes through the atmosphere, attenuation barely happens. Therefore, the boundary conditions of ABCD can be replaced by A
′
B
′
C
′
D
′s in the TE model.


## 2.1. The Boundary Value Problem of MT Method

The electromagnetic waves from the magnetotelluric field source vertically incident into the underground as harmonic plane waves with a time harmonic factore
-
i
ω
t. Maxwell coupled magnetotelluric field equations in frequency domain as
(1)
∇
×
E
=
i
ω
μ
H
,
∇
×
H
=
(
σ
-
i
ω
ε
)
E
,
where E is the electric vector, H is the magnetic vector, ω is angular frequency, μ is the medium permeability when the permeability of air is 1.2567
×
1
0
-
6
H
/
m, σ is the medium conductivity or the reciprocal of medium resistivity, ε is the dielectric medium constant, and the dielectric constant of air is 8.85
×
1
0
-
12 F/m.Curling both sides of the electric equations (1), the differential electric equation can be obtained as follows:
(2)
∇
×
(
∇
×
E
)
-
k
2
E
=
0
,
∇
×
(
∇
×
H
)
-
k
2
H
=
0
,
where
(3)
k
2
=
i
ω
μ
σ
+
ω
2
ε
μ
.The equations of electric and magnetic have the same form that can be written in one unification equation. SubstitutingE and H for U yields
(4)
∇
×
(
∇
×
U
)
-
k
2
U
=
0
.The vectorU can be rewritten to the vector form with three components:
(5)
U
=
U
x
e
x
+
U
y
e
y
+
U
z
e
z
,
where, e
x, e
y, and e
z are the unit vectors in the directions of x, y, and z.


### 2.1.1. When It Is Working in TE Model

(a)
In the surface of ABCD shown in Figure1, the field source is the electric vector E
x. As the ratio of electric and magnetic, U
x has no impact to the final result. Therefore, U
x can be either zero or nonzero value. Here U
x
=
1. U
y and U
z are all zero as well since there are no components of U in x and y directions.
(b)
In the four vertical boundaries shown in Figure1, the propagation direction of EM wave in magnetic field is vertical downward, being perpendicular to normal direction of the boundary, which is E
×
H
⊥
Γ. Therefore, ∂
U
/
∂
n equals zero.
(c)
In the surface of EFGH in Figure1, the underground can be assumed as a homogenous medium with a consideration that this surface is beyond the study area. The electromagnetic field attenuates in a normal way, which follows the index law during downward propagation. The vector U has no value in directions y and z. 													(6)
U
x
=
c
e
-
χ
z
,
U
y
=
0
,
U
z
=
0
,where c is a constant factor with any nonzero value, χ
=
-
i
ω
μ
σ, and σ is the homogenous conductivity below the infinite surface EFGH. Derivative of (6) is obtained as follows:
(7)
∂
U
x
∂
n
=
∂
U
x
∂
z
=
-
χ
U
x
.


### 2.1.2. When It Is Working in TM Model

The field source is one component of the magnetic vectorH
x. When magnetic field passes through the atmosphere, attenuation barely happens. Therefore, the boundary conditions of ABCD can be replaced by A
′
B
′
C
′
D
′s in the TE model.


## 2.1.1. When It Is Working in TE Model

(a)
In the surface of ABCD shown in Figure1, the field source is the electric vector E
x. As the ratio of electric and magnetic, U
x has no impact to the final result. Therefore, U
x can be either zero or nonzero value. Here U
x
=
1. U
y and U
z are all zero as well since there are no components of U in x and y directions.
(b)
In the four vertical boundaries shown in Figure1, the propagation direction of EM wave in magnetic field is vertical downward, being perpendicular to normal direction of the boundary, which is E
×
H
⊥
Γ. Therefore, ∂
U
/
∂
n equals zero.
(c)
In the surface of EFGH in Figure1, the underground can be assumed as a homogenous medium with a consideration that this surface is beyond the study area. The electromagnetic field attenuates in a normal way, which follows the index law during downward propagation. The vector U has no value in directions y and z. 													(6)
U
x
=
c
e
-
χ
z
,
U
y
=
0
,
U
z
=
0
,where c is a constant factor with any nonzero value, χ
=
-
i
ω
μ
σ, and σ is the homogenous conductivity below the infinite surface EFGH. Derivative of (6) is obtained as follows:
(7)
∂
U
x
∂
n
=
∂
U
x
∂
z
=
-
χ
U
x
.


## 2.1.2. When It Is Working in TM Model

The field source is one component of the magnetic vectorH
x. When magnetic field passes through the atmosphere, attenuation barely happens. Therefore, the boundary conditions of ABCD can be replaced by A
′
B
′
C
′
D
′s in the TE model.


## 3. The Variation Problem of Weighted Residual Equation

Both sides of (4) are multiplied by δ
U and integrated on the whole region, showing the result as follows:
(8)
∫
Ω
[
∇
×
(
∇
×
U
)
-
k
2
U
]
·
δ
U
d
Ω
=
0
,
where
(9)
δ
U
=
δ
U
x
e
x
+
δ
U
y
e
y
+
δ
U
z
e
z
.According to the vector identities(10)
(
∇
×
A
)
·
B
=
∇
·
(
A
×
B
)
+
(
∇
×
B
)
·
A
.Equation (8) will be
(11)
∫
Ω
∇
×
(
∇
×
U
)
·
δ
U
d
Ω
=
∫
Ω
∇
·
[
(
∇
×
U
)
×
δ
U
]
d
Ω
+
∫
Ω
∇
×
δ
U
·
∇
×
U
d
Ω
.In which,(12)
∫
Ω
∇
·
[
(
∇
×
U
)
×
δ
U
]
d
Ω
=
∫
Γ
(
∇
×
U
)
×
δ
U
·
d
Γ
.According to the boundary condition(a)
in the surfaces ABCD andA
′
B
′
C
′
D
′, δ
U
=
0 due to fixed U in these surfaces;
(b)
in the four infinite vertical surfaces,(
∇
×
U
)
×
δ
U
·
d
Γ
=
0, because the reflection of the anomalous field can be ignored. Then the electric field U and δ
U have one component x only;
(c)
in the surface EFGH,(
∇
×
U
)
×
δ
U
=
-
(
∂
U
x
/
∂
z
)
δ
U
x
e
z.Substituting all the boundary conditions in (12) yields
(13)
∫
Ω
∇
×
δ
U
·
∇
×
U
d
Ω
-
∫
Ω
k
2
U
·
δ
U
d
Ω
+
∫
EFGH
χ
U
x
δ
U
x
d
Γ
=
0
.


## 4. Hexahedral Mesh Generation

Hexahedral element mesh is used to make Figure1 into a 3D geology space as shown in Figure 2. Obviously, the entire space of 3D geology is meshed into the amount of hexahedral elements.Figure 2
Hexahedral mesh.The vertices of each hexahedral element are sorted as shown in Figure3.Figure 3
The vertices of hexahedral mesh.


## 5. Unit Analysis

Using the finite element method to solve differential equations requires a mesh representation of the target area, such as the widely used hexahedral mesh. Figure4(a) is the parent element of the hexahedron and Figure 4(b) is the subelement. The coordinate correspondences between them are
(14)
x
=
x
0
+
a
2
ξ
,
y
=
y
0
+
b
2
η
,
z
=
z
0
+
c
2
ζ
,
where x
0, y
0, and z
0 is the center coordinate of sub element; a, b, c is the lengths of three edges. The differential relations between them are
(15)
d
x
=
a
2
d
ξ
,
d
y
=
b
2
d
η
,
d
z
=
c
2
d
ζ
.(a) The parent unit of hexahedron; (b) The subelement of hexahedron.
(a)
(b)The shape function is defined as(16)
N
i
=
1
8
(
1
+
ξ
i
ξ
)
(
1
+
η
i
η
)
(
1
+
ζ
i
ζ
)
,
where ξ
i, η
i, ζ
i are the coordinates of vertex i in the parent unit.The interpolation function is(17)
E
x
=
N
T
E
x
,
E
y
=
N
T
E
y
,
E
z
=
N
T
E
z
,
where N
T
=
(
N
1
,
…
,
N
8
), E
x
T
=
(
E
x
1
,
…
,
E
x
8
), E
y
T
=
(
E
y
1
,
…
,
E
y
8
), E
z
T
=
(
E
z
1
,
…
,
E
z
8
). Therefore
(18)
E
=
N
T
(
E
x
e
x
+
E
y
e
y
+
E
z
e
z
)
δ
E
=
N
T
(
δ
E
x
e
x
+
δ
E
y
e
y
+
δ
E
z
e
z
)
,
where δ
E
x, δ
E
y, and δ
E
z are the variations of E
x, E
y, and E
z.


## 6. The Form of Unit Coefficient Matrix

The regional integral of (13) can be decomposed into the integral of each unit:
(19)
∑
Ω
∫
e
∇
×
δ
E
·
∇
×
E
d
Ω
-
∑
Ω
∫
Ω
k
2
E
·
δ
E
d
Ω
+
∑
EFGH
∫
1234
χ
E
x
δ
E
x
d
Γ
=
0
.In one unite
(20)
∇
×
E
=
E
e
T
(
α
e
x
+
β
e
y
+
γ
e
z
)
,
∇
×
δ
E
=
δ
E
e
T
(
α
e
x
+
β
e
y
+
γ
e
z
)
,
where
(21)
E
e
T
=
(
E
x
T
,
E
y
T
,
E
z
T
)
,
δ
E
e
T
=
(
δ
E
x
T
,
δ
E
y
T
,
δ
E
z
T
)
α
T
=
[
0
,
-
(
∂
N
∂
z
)
T
,
(
∂
N
∂
y
)
T
]
,
β
T
=
[
-
(
∂
N
∂
z
)
T
,
0
,
(
∂
N
∂
x
)
T
]
,
γ
T
=
[
-
(
∂
N
∂
y
)
T
,
(
∂
N
∂
x
)
T
,
0
]
,
(
∂
N
∂
x
)
T
=
(
∂
N
1
∂
x
,
…
,
∂
N
8
∂
x
)
∂
N
i
∂
x
=
∂
N
i
∂
ξ
2
a
,
(
∂
N
∂
y
)
T
=
(
∂
N
1
∂
y
,
…
,
∂
N
8
∂
y
)
∂
N
i
∂
y
=
∂
N
i
∂
η
2
b
,
(
∂
N
∂
z
)
T
=
(
∂
N
1
∂
z
,
…
,
∂
N
8
∂
z
)
∂
N
i
∂
z
=
∂
N
i
∂
ζ
2
c
,
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
i
=
1
,
…
,
8
.
Substituting (20) into the first integral of (19)
(22)
∫
e
∇
×
δ
E
·
∇
×
E
d
Ω
=
δ
E
e
T
K
1
e
E
e
,
where
(23)
K
1
e
=
∫
e
(
α
α
T
+
β
β
T
+
γ
γ
T
)
d
Ω
.
Its matrix form is
(24)
K
1
e
=
[
K
11
K
12
K
13
K
21
K
22
K
23
K
31
K
32
K
33
]
,
K
11
=
∫
e
[
(
∂
N
∂
y
)
(
∂
N
∂
y
)
T
+
(
∂
N
∂
z
)
(
∂
N
∂
z
)
T
]
d
Ω
,
K
22
=
∫
e
[
(
∂
N
∂
x
)
(
∂
N
∂
x
)
T
+
(
∂
N
∂
z
)
(
∂
N
∂
z
)
T
]
d
Ω
,
K
33
=
∫
e
[
(
∂
N
∂
x
)
(
∂
N
∂
x
)
T
+
(
∂
N
∂
y
)
(
∂
N
∂
y
)
T
]
d
Ω
,
K
21
=
-
∫
e
(
∂
N
∂
x
)
(
∂
N
∂
y
)
T
d
Ω
,
K
12
=
(
K
21
)
T
,
K
31
=
-
∫
e
(
∂
N
∂
x
)
(
∂
N
∂
z
)
T
d
Ω
,
K
13
=
(
K
31
)
T
,
K
32
=
-
∫
e
(
∂
N
∂
y
)
(
∂
N
∂
z
)
T
d
Ω
,
K
23
=
(
K
32
)
T
.The expressions of the elements in each matrix are listed as(25)
K
1
1
i
j
=
∫
e
(
∂
N
i
∂
η
∂
N
j
∂
η
4
b
2
+
∂
N
i
∂
ζ
∂
N
j
∂
ζ
4
c
2
)
a
b
c
8
d
ξ
η
ζ
=
1
36
(
a
c
b
β
i
j
+
a
b
c
v
i
j
)
,
K
2
2
i
j
=
∫
e
(
∂
N
i
∂
ξ
∂
N
j
∂
ξ
4
a
2
+
∂
N
i
∂
ζ
∂
N
j
∂
ζ
4
c
2
)
a
b
c
8
d
ξ
η
ζ
=
1
36
(
b
c
a
α
i
j
+
a
b
c
v
i
j
)
,
K
3
3
i
j
=
∫
e
(
∂
N
i
∂
ξ
∂
N
j
∂
ξ
4
a
2
+
∂
N
i
∂
η
∂
N
j
∂
η
4
b
2
)
a
b
c
8
d
ξ
η
ζ
=
1
36
(
b
c
a
α
i
j
+
a
c
b
β
i
j
)
,
K
2
1
i
j
=
-
∫
e
∂
N
i
∂
ξ
∂
N
j
∂
η
c
2
d
ξ
η
ζ
=
-
c
24
m
i
j
,
K
3
1
i
j
=
-
∫
e
∂
N
i
∂
ξ
∂
N
j
∂
ζ
b
2
d
ξ
η
ζ
=
-
b
24
n
i
j
,
K
3
2
i
j
=
-
∫
e
∂
N
i
∂
η
∂
N
j
∂
ζ
a
2
d
ξ
η
ζ
=
-
a
24
l
i
j
,
where
(26)
α
i
j
=
|
4
2
1
2
-
4
-
2
-
1
-
2
4
2
1
-
2
-
4
-
2
-
1
4
2
-
1
-
2
-
4
-
2
4
-
2
-
1
-
2
-
4
4
2
1
2
4
2
1
4
2
4
|
,
β
i
j
=
|
4
2
-
2
-
4
2
1
-
1
-
2
4
-
4
-
2
1
2
-
2
-
1
4
2
-
1
-
2
2
1
4
-
2
-
1
1
2
4
2
-
2
-
4
4
-
4
-
2
4
2
4
|
,
γ
i
j
=
|
4
-
4
-
2
2
2
-
2
-
1
1
4
-
4
-
2
2
2
-
2
-
1
4
-
4
-
2
2
2
-
2
4
-
4
-
2
2
2
4
-
4
-
2
2
4
2
-
2
4
-
4
4
|
,
m
i
j
=
|
2
1
-
1
-
2
2
1
-
1
-
2
1
2
-
2
-
1
1
2
-
2
-
1
1
2
-
2
-
1
1
2
-
2
-
1
2
1
-
1
-
2
2
1
-
1
-
2
-
2
-
1
1
2
-
2
-
1
1
2
-
1
-
2
2
1
-
1
-
2
2
1
-
1
-
2
2
1
-
1
-
2
2
1
-
2
-
1
1
2
-
2
-
1
1
2
|
,
n
i
j
=
|
2
-
2
-
1
1
2
-
2
-
1
1
2
-
2
-
1
1
2
-
2
-
1
1
1
-
1
-
2
2
1
-
1
-
2
2
1
-
1
-
2
2
1
-
1
-
2
2
-
2
2
1
-
1
-
2
2
1
-
1
-
2
2
1
-
1
-
2
2
1
-
1
-
1
1
2
-
2
-
1
1
2
-
2
-
1
1
2
-
2
-
1
1
2
-
2
|
,
l
i
j
=
|
2
-
2
-
2
2
1
-
1
-
1
1
2
-
2
-
2
2
1
-
1
-
1
1
-
2
2
2
-
2
-
1
1
1
-
1
-
2
2
2
-
2
-
1
1
1
-
1
1
-
1
-
1
1
2
-
2
-
2
2
1
-
1
-
1
1
2
-
2
-
2
2
-
1
1
1
-
1
-
2
2
2
-
2
-
1
1
1
-
1
-
2
2
2
-
2
|
.Substituting (18) into the second integral of (19):
(27)
∫
e
k
2
E
·
δ
E
d
Ω
=
δ
E
e
T
K
2
e
E
e
,
where
(28)
K
2
e
=
[
K
11
0
0
0
K
22
0
0
0
K
33
]
,
K
11
=
K
22
=
K
33
=
∫
e
k
2
N
N
T
d
Ω
,
K
1
1
i
j
=
∫
e
k
2
N
i
N
j
a
b
c
8
d
ξ
d
η
d
ζ
,
i
,
j
=
1
,
…
,
8
=
1
144
a
b
c
k
2
o
i
j
,
o
i
j
=
|
6
3
3
6
2
1
1
2
6
6
3
1
2
2
1
6
3
1
2
2
1
6
2
1
1
2
6
3
3
6
6
6
3
6
3
6
|
.When the boundary surface unit 1234 is in the boundary surface EFGH, the boundary integral is(29)
∫
1234
χ
E
x
δ
E
z
d
Γ
=
δ
E
e
T
K
3
e
E
e
,
where
(30)
K
3
e
=
[
K
11
0
0
0
0
0
0
0
0
]
,
K
11
=
∫
1234
χ
N
N
T
d
Γ
,
K
1
1
i
j
=
∫
1234
χ
N
i
N
j
a
b
4
d
ξ
d
η
,
i
,
j
=
1
,
…
,
4
=
1
36
a
b
χ
p
i
j
,
p
i
j
=
|
4
2
1
2
4
2
1
4
2
4
|
,
where N
i, N
j are the shape functions of quadrilateral plane element rather than the shape functions of hexahedral plane element. Its form is shown in the following:
(31)
N
i
=
1
4
(
1
+
ξ
i
ξ
)
(
1
+
η
i
η
)
.Substituting (14) into (22) and (27), (31) into (29), respectively, each element of the matrix can then be calculated.


## 7. The Integration of the Coefficient Matrix

Integrating all the hexahedral elements obtained above gives the integration of the coefficient matrix which is shown below:(32)
∑
∂
U
e
T
(
K
1
e
-
K
2
e
+
K
3
e
)
U
e
=
∑
∂
U
e
T
K
e
U
e
=
∂
U
T
∑
K
e
U
=
∂
U
T
K
U
=
0
.The arbitrary∂
U
T is not always zero, which gives the following large-scale equations:
(33)
K
U
=
0
,
where K is the coefficient matrix and U is the unknown variables.


## 8. The Computation of Apparent Resistivity

The expansions of (2) in the Cartesian coordinate system are presented below:
(34)
H
x
=
1
i
ω
μ
(
∂
E
z
∂
y
-
∂
E
y
∂
z
)
,
H
y
=
1
i
ω
μ
(
∂
E
z
∂
z
-
∂
E
z
∂
x
)
,
H
z
=
1
i
ω
μ
(
∂
E
y
∂
x
-
∂
E
z
∂
y
)
.WhenE
x, E
y, and E
z of each vertex have been obtained, H
x, H
y, and H
z can be computed from the approximate values of each vertex’s differential quotient using the difference method. Therefore, the apparent resistivity equation is derived as:
(35)
ρ
x
y
=
i
ω
μ
(
E
z
H
y
)
2
.


## 9. Example Analysis

### 9.1. Homogeneous Model

Figure5 shows a homogeneous model: the resistivity is 100
Ω
·
m; the grid size is 9000
×
9000m2; the dot pitch is 1000 m; the line spacing is 1000 m. The measuring frequency is 1000 Hz, 500 Hz, 100 Hz, and 50 Hz, and the observation points are the five points of the midline (2000 m, 3000 m, 4000 m, 5000 m, and 6000 m).Figure 5
Homogeneous model.Table1 lists the numerical computation results. According to the results, the average of the relative error is 2.81% in 1000 Hz and 0.16% in 100 Hz. Gradually decreasing relative error is observed as the frequency decreases.Table 1
The result of numerical computation.
X/Km
Fre/Hz
Rho/Ω·m
Phase/(°)
RE/(%)
2000
1000
97.17
44.67
2.83
2000
500
98.26
44.66
1.74
2000
100
100.17
44.82
0.17
2000
50
99.36
45.11
0.64
3000
1000
97.16
44.68
2.84
3000
500
98.24
44.66
1.76
3000
100
100.14
44.83
0.14
3000
50
99.29
45.13
0.71
4000
1000
97.17
44.68
2.83
4000
500
98.24
44.66
1.76
4000
100
100.13
44.83
0.13
4000
50
99.23
45.14
0.77
5000
1000
97.19
44.68
2.81
5000
500
98.26
44.66
1.74
5000
100
100.15
44.84
0.15
5000
50
99.26
45.12
0.74
6000
1000
97.25
44.67
2.75
6000
500
98.34
44.65
1.66
6000
100
100.21
44.85
0.21
6000
50
99.4
45.09
0.6


### 9.2. Positive Landform Model

In Figure6, a positive landform model has a resistivity of 100
Ω
·
m and the topographic fall is 60 m; the grid size is 9000
×
9000m2; the dot pitch is 1000 m; and the line spacing is 1000 m. The measuring frequencies are 1000 Hz and 100 Hz, and the observation points are the nine points on the midline.Figure 6
Positive landform model.Results of the numerical computation are shown in Figure7. The highest value is the pseudoabnormality in the positive landform. The distortion error is 39.5% at 1000 Hz and 12.8% at 100 Hz. With a decreasing frequency, the effect of topographic distortion decreases gradually.Figure 7
The result of numerical computation.


### 9.3. The Model with Low Resistance Body

Figure8 shows a model containing a low resistance body with 10
Ω
·
m resistivity; resistivity of the rest of the model is 100
Ω
·
m. The horizontal geometry of the anomalous body is 3000
m
×
3000
m; upper surface is 250 m away from the ground and its height is 750 m; the grid size is 9000
×
9000m2; the dot pitch is 1000 m; the line spacing is 1000 m (Figure 9). The measuring frequencies are 1000 Hz and 100 Hz, and the observation points are the nine points of the mid-line.Figure 8
The model with low resistance body.(a) Horizontal top view. (b) Positive side view.
(a)
(b)The result of numerical computation is shown in Figure10 with the curves showing low resistivity abnormality at both 1000 Hz and 100 Hz.Figure 10
The result of numerical computation.


## 9.1. Homogeneous Model

Figure5 shows a homogeneous model: the resistivity is 100
Ω
·
m; the grid size is 9000
×
9000m2; the dot pitch is 1000 m; the line spacing is 1000 m. The measuring frequency is 1000 Hz, 500 Hz, 100 Hz, and 50 Hz, and the observation points are the five points of the midline (2000 m, 3000 m, 4000 m, 5000 m, and 6000 m).Figure 5
Homogeneous model.Table1 lists the numerical computation results. According to the results, the average of the relative error is 2.81% in 1000 Hz and 0.16% in 100 Hz. Gradually decreasing relative error is observed as the frequency decreases.Table 1
The result of numerical computation.
X/Km
Fre/Hz
Rho/Ω·m
Phase/(°)
RE/(%)
2000
1000
97.17
44.67
2.83
2000
500
98.26
44.66
1.74
2000
100
100.17
44.82
0.17
2000
50
99.36
45.11
0.64
3000
1000
97.16
44.68
2.84
3000
500
98.24
44.66
1.76
3000
100
100.14
44.83
0.14
3000
50
99.29
45.13
0.71
4000
1000
97.17
44.68
2.83
4000
500
98.24
44.66
1.76
4000
100
100.13
44.83
0.13
4000
50
99.23
45.14
0.77
5000
1000
97.19
44.68
2.81
5000
500
98.26
44.66
1.74
5000
100
100.15
44.84
0.15
5000
50
99.26
45.12
0.74
6000
1000
97.25
44.67
2.75
6000
500
98.34
44.65
1.66
6000
100
100.21
44.85
0.21
6000
50
99.4
45.09
0.6


## 9.2. Positive Landform Model

In Figure6, a positive landform model has a resistivity of 100
Ω
·
m and the topographic fall is 60 m; the grid size is 9000
×
9000m2; the dot pitch is 1000 m; and the line spacing is 1000 m. The measuring frequencies are 1000 Hz and 100 Hz, and the observation points are the nine points on the midline.Figure 6
Positive landform model.Results of the numerical computation are shown in Figure7. The highest value is the pseudoabnormality in the positive landform. The distortion error is 39.5% at 1000 Hz and 12.8% at 100 Hz. With a decreasing frequency, the effect of topographic distortion decreases gradually.Figure 7
The result of numerical computation.


## 9.3. The Model with Low Resistance Body

Figure8 shows a model containing a low resistance body with 10
Ω
·
m resistivity; resistivity of the rest of the model is 100
Ω
·
m. The horizontal geometry of the anomalous body is 3000
m
×
3000
m; upper surface is 250 m away from the ground and its height is 750 m; the grid size is 9000
×
9000m2; the dot pitch is 1000 m; the line spacing is 1000 m (Figure 9). The measuring frequencies are 1000 Hz and 100 Hz, and the observation points are the nine points of the mid-line.Figure 8
The model with low resistance body.(a) Horizontal top view. (b) Positive side view.
(a)
(b)The result of numerical computation is shown in Figure10 with the curves showing low resistivity abnormality at both 1000 Hz and 100 Hz.Figure 10
The result of numerical computation.


## 10. Conclusions and Recommendations

During the research of magnetotelluric finite element numerical simulation algorithm,(
1
) the boundary values of 3D magnetotelluric numerical simulations are analyzed; (
2
) the preconditioned conjugate gradient is adopted for solving large-scale equations under the first-type boundary conditions; (
3
) the stability and the iterative convergence speed are ensured in 3D magnetotelluric forward calculations. This paper provides an effective tool for theoretical studies on distribution and characteristics of magnetotelluric field response of a 3D geology structure. Meanwhile, it also has established critical foundation for future 3D inversion algorithm research.



---

*Source: 890516-2014-08-06.xml*