# 2D Multitransient Electromagnetic Response Modeling of South China Shale Gas Earth Model Using an Approximation of Finite Difference Time Domain with Uniaxial Perfectly Matched Layer

Authors: Olalekan Fayemi; Qingyun Di
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

In this study, we introduced multitransient electromagnetic (MTEM) method as an effective tool for shale gas exploration. We combined the uniaxial perfectly matched layer (UPML) equation with the first derivative diffusion equation to solve for a finite difference time domain (FDTD) UPML equation, which was discretized to form an algorithm for 3D modeling of earth impulse response and used in modeling MTEM response over 2D South China shale gas model. We started with stepwise demonstration of the UPML and the FDTD algorithm as an effective tool. Subsequently, quantitative study on the convergence of MTEM earth impulse response was performed using different grid setup over a uniform earth material. This illustrates that accurate results can be obtained for specified range of offset. Furthermore, synthetic responses were generated for a set of geological scenarios. Lastly, the FDTD algorithm was used to model the MTEM response over a 2D shale gas earth model from South China using a PRBS source. The obtained apparent resistivity section from the MTEM response showed a similar geological setup with the modeled 2D South China shale gas section. This study confirmed the competence of MTEM method as an effective tool for unconventional shale gas prospecting and exploitation.


---

## Body

## 1. Introduction

Shale gas as the natural gas that is found trapped within the shale formation has become an important source of natural gas in China and a major energy source in the rest of the world. China has the largest estimate of recoverable shale gas in the world, which is estimated at 31.57 trillion cubic meters [1–3]. The availability of these abundant resources as well as the determination to decrease the air pollution (CO2 emissions) created by coal burning, which is also poor in efficiency, has made the government focus on the developing technology in shale gas recovery. According to Nie et al. [4], the most favorable areas for shale gas accumulation and recovery are located in Sichuan Basin, Micangshan-Dabashan foreland and East Chongqing, West Hubei Province of Upper-Yangtze River, North Hubei Province of Middle-Yangtze River, and South Jiangsu Province of Lower-Yangtze River. Two of the three major shale gas reserve basins which account for the largest percentage of the countries reserve are found in South China. Also, research results show that the southern region has shale gas producing geological conditions [4, 5]. Hence, this is favorable area for China’s shale gas research and exploration. Therefore, in this paper, we studied the application of a relatively new EM geophysical technique used to reduce the risk involved in hydrocarbon prospecting. We determined its applicability in shale gas exploration (in unconventional petroleum systems), which is intended at mapping the unconventional source rocks and exploitation process (defined as the application of technology to increase the recovery of shale gas in existing reservoir) in South China.Although hydrocarbon prospecting techniques are more developed now than in past decades, extraction process is still a major challenge in the Chinese shale gas industry. Nonetheless, the development in oil recovery process such as horizontal drilling and hydraulic fracturing has played very important role in unconventional reservoir exploitation for abundant shale gas reserve for energy resources [6]. Production from low permeable shale plays can be enhanced by using the above mentioned technologies: (
1
) horizontal drilling that increases rock volume contact with wellbore thereby providing greater access to the gas trapped deep in the producing formation and (
2
) hydraulic fracture stimulation that increases permeability (by opening cracks (fractures) in the rock and allowing natural gas to flow from the shale into the well) and requires accurate understanding of the subsurface earth model. This can be achieved by analyzing the sensitivity of various rock parameters to geophysical methods via numerically generated synthetic geophysical data [7]. Synthetic responses can be generated for a set of geological setups using available rock physics/properties relationship and forward modeling. Most common geophysical data used to provide accurate 3D subsurface images are seismic data (which provides information about stress/fracture orientations [8]) and electromagnetic (EM) data (which provides information about subsurface geology and fluid saturation in relation to changes in resistivity [6, 9]). Hence, in this study, we introduced an approximation of finite difference time domain (FDTD) algorithm with uniaxial perfectly matched layer (UPML) for EM data modeling.FDTD method has been vastly used in electromagnetic (EM) modeling techniques [10–12]. This is because of its simplicity (i.e., straight forward solution) and its ability to model relatively complex structures. Both explicit solutions and implicit solution had been used to obtain numerical model of time domain transient electromagnetic response [10–14] within a specified region (computational region). However, the propagation of EM signal is indefinite irrespective of mesh termination. Therefore, it is essential to simulate open space region to accommodate the unlimited simulation signal, by initializing a boundary condition. The perfectly matched layers (PML) method proposed by Berenger [15] for the finite difference time domain (FDTD) method has been the most effective absorbing boundary condition (ABC), especially when a structure involving wave absorption of incident angles from 0 to near 90 degrees is being modeled [16, 17]. Over the years, several works have been done on the application of perfectly matched layer equation, which has led to improvements in this method and the development of more effective techniques such as stretched coordinate PML [18], uniaxial PML [19, 20], convolution PML [21], and nonsplit complex frequency shifted PML [22]. However, understanding the effect of diffusive propagation requires tools that are very different from those that are useful for understanding wave propagation. Hence, we harnessed the advantages of stable UPML (with physical interpretation) to compute a FDTD algorithm in which the model can be formulated and implemented without considering the matched layer. Using the UPML formulation and combining it with the 3D formulation for Maxwell’s first derivative (coupled) equation [11], explicit and stable FDTD approximate equation for EM diffusion equation was derived.The application of EM methods to complement seismic data and minimize exploration risk has proven to be very effective and is fast becoming a regular practice in the E&P industry. This is due to the limitations faced by conventional reflection seismic technique in some areas [3, 23, 24]. EM application in hydrocarbon exploration has grown rapidly in the past two decades from just MT study to the establishment of controlled source electromagnetic method (CSEM) as a major tool in offshore exploration [23, 25], long offset transient electromagnetic (LOTEM) method [26] as an onshore (land) EM technique, and multitransient electromagnetic (MTEM) method, for both land and marine studies [27]. MTEM is a relatively new time domain method which gives a high-resolution result [28], continuous 2D view of the subsurface resistivity distribution across a section [29], 3D volume resistivity image to supplement high-resolution seismic data in decision-making, and risk reduction [28] and has the capability to help maximize recoverable reserves and CO2 sequestration process, via 4D survey monitoring of water, gas, and/or chemicals injected into reservoirs during enhanced recovery programs [30, 31]. The comparison of both CSEM and MTEM methods by Anderson et al. [28] shows that the latter, which is used for both land and marine study (including deep offshore study), has several advantages over CSEM method. Some of the early EM methods such as induced polarization [32, 33], magnetotelluric (MT) methods [34], and etcetera (including DC method which is inform of well logging) represent precursor to the multitransient EM (MTEM) methodology operated by PGS, thus giving it a combined advantage over the other controlled source EM method.To start with, brief description of 3D FDTD formulation was given in the next section, which is followed by stepwise demonstration of the UPML and the FDTD algorithm as an effective tool for MTEM earth response modeling. The final exercise in this study established MTEM as an effective and efficient tool used in subsurface characterization for shale gas prospecting and exploitation via quick reconstruction of earth resistivity distribution.


## 2. Formulation of 3D Finite Difference Approximation

In 1994, Berenger [15] described an Absolute Boundary Condition (ABC) called perfectly matched layer (PML). This method involves the simulation of a region where the incoming wave is attenuated. The PML region serves as an absorbing layer for the incoming fields at every angle. Reflections entering the lossy regions are prevented because impedance is matched. Likewise, reflections from the grid boundary are prevented because the outgoing waves are absorbed [12]. However, the PML treatment described by Berenger relies on the split field formulation of Maxwell’s equations in the PML region and has been shown to be weakly stable [35]. This approach involves significant modification of standard Maxwell’s equations. Therefore, we adopted the uniaxial PML (UPML) which was derived based on the anisotropic unsplit perfectly matched layer (APML) described by Sacks et al. [36] and Gedney [20]. This does not require the splitting of Maxwell’s fields, hence, maintaining the strong well-posedness of original Maxwell’s system both in frequency and in time domain [36]. Maintaining the UPML advantages, a stable FDTD approximation was derived for EM diffusion equation.


### 2.1. Approximations within the UPML Region

The essence of a UPML is to achieve zero reflection at the boundary of interest. Therefore, four (4) assumptions were made to ensure zero reflection at both sides of the boundary. They are as follows:(1)
The loss within the matched layer is set to be greater than zero; that is,k
>
0.
(2)
The impedance across the boundary is set as 1.
(3)
[μ
r] = [ε
r] = [s
r], having the matrix form s
a
0
0
0
s
b
0
0
0
s
c where s is the fictitious electrical and magnetic property within the matched layer. Thus, θ
1 = θ
2 and n
1 = n
2.
(4)
Consider a uniaxial isotropic medium wheres
a = s
b = 1
/
s
c and s
b
s
c
=
1 for [s
c].Therefore, the UPML parameters can be written as(1)
s
x
=
s
x
-
1
0
0
0
s
x
0
0
0
s
x
,
s
y
=
s
y
0
0
0
s
y
-
1
0
0
0
s
y
,
s
z
=
s
z
0
0
0
s
z
0
0
0
s
z
-
1
,
s
r
=
s
x
·
s
y
·
s
z
=
s
y
s
z
s
x
0
0
0
s
x
s
z
s
y
0
0
0
s
x
s
y
s
z
,where s
x, s
y, and s
z are arbitrary parameters calculated using the mathematical expressions given,(2)
S
x
=
1
+
σ
x
j
ω
ε
0
,
σ
x
=
ε
0
2
Δ
t
x
L
x
p
,
S
y
=
1
+
σ
y
j
ω
ε
0
,
σ
y
=
ε
0
2
Δ
t
y
L
y
p
,
S
z
=
1
+
σ
z
j
ω
ε
0
,
σ
z
=
ε
0
2
Δ
t
z
L
z
p
,assuming the k
m
a
x (maximum linear absorption factor) equals 1. Δ
t is the time step to be used, r
=
x, y, or z, and p is either 3 or 4, while r
/
L
r
p is a function grid space and it ensures that the value of the arbitrary parameter increases outwardly (Figure 1). Figure 1 shows an image of a grid whose PML region is represented by the black stripes that highlight its cells and its main medium is colored blue.Figure 1
3D spatial grid with perfectly matched layer.


### 2.2. Introducing the UPML into the EM Diffusion Equation

Maxwell’s equation for EM fields in linear, isotropic, and source-free media, under the quasistatic approximation in the frequency domain is given as follows:(3)
∇
×
E
r
,
ω
=
-
j
ω
μ
H
r
,
ω
,
(4)
∇
×
H
r
,
ω
=
σ
r
E
r
,
ω
,where E and H are the electric and the magnetic fields, respectively, and μ and σ are the physical properties, namely, magnetic permeability and electric conductivity, respectively. “r” represents the 3D Cartesian coordinate system (x, y, and z), while ω is the angular frequency (
2
π
f
). Although the conventional EM diffusion equation to be solved is parabolic (i.e., the conventional time domain equation neglects the displacement current), however, in this study we solved an EM hyperbolic equation (coupled Maxwell’s equation which simulates the direct interaction between the magnetic and electric field) for the FDTD algorithm by incorporating a fictitious displacement current to ensure an algorithm that is nonerroneously explicit and unconditionally stable [11, 13] for both 2D and 3D computation. This gives a wave-like equation in a lossy medium but by choosing the appropriate fictitious value for electric permittivity, an explicit, always stable, time stepping solution is obtained. To ensure the stability and convergence of a diffused EM wave initialized by an electric dipole with current density J
s, an arbitrary value of electric permittivity given as(5)
γ
=
3
μ
m
i
n
Δ
t
Δ
m
i
n
2was used to replace the real permittivity [11]. Also, the time step was calculated from the equation given in the following: (6)
Δ
t
≤
μ
m
i
n
σ
a
1
/
2
min
⁡
Δ
x
,
Δ
y
,
Δ
z
2
,where a = 2, 4, and 6 for 1D, 2D, and 3D, respectively.This keeps the velocity of the fictitious wave field slower than that simulated by the finite difference scheme and prevents the fictitious displacement current from dominating the diffusive EM field by constraining the length of the time step.Therefore, (4) becomes(7)
∇
×
H
r
,
ω
=
σ
b
r
E
r
,
ω
+
j
ω
γ
E
r
,
ω
+
J
s
r
,
ω
.By incorporating the UPML parameters (from (1) and (2)) into (7), we obtain (8)
∇
×
H
r
,
ω
=
σ
b
r
E
r
,
ω
+
j
ω
γ
s
r
E
r
,
ω
+
J
s
r
,
ω
.Thus,(9)
∇
×
H
x
r
,
ω
=
σ
b
r
E
x
r
,
ω
+
j
ω
γ
1
+
σ
y
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
1
+
σ
x
j
ω
ε
0
-
1
E
x
r
,
ω
+
J
s
r
,
ω
,
for
E
x
,
(10)
∇
×
H
y
r
,
ω
=
σ
b
r
E
y
r
,
ω
+
j
ω
γ
1
+
σ
x
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
1
+
σ
y
j
ω
ε
0
-
1
E
y
r
,
ω
+
J
s
r
,
ω
,
for
E
y
,
(11)
∇
×
H
z
r
,
ω
=
σ
b
r
E
z
r
,
ω
+
j
ω
γ
1
+
σ
x
j
ω
ε
0
1
+
σ
y
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
-
1
E
z
r
,
ω
+
J
s
r
,
ω
,
for
E
z
.Computing for E
x, (9) was multiplied by the UPML parameter S
x and then divided by fictitious permittivity, γ, before expanding the operation/equation to obtain(12)
1
γ
∇
×
H
x
r
,
ω
+
σ
x
j
ω
ε
0
γ
∇
×
H
x
r
,
ω
=
σ
b
r
γ
E
x
r
,
ω
+
σ
b
r
σ
x
j
ω
ε
0
γ
E
x
r
,
ω
+
j
ω
E
x
r
,
ω
+
j
ω
σ
y
+
σ
z
j
ω
ε
0
E
x
r
,
ω
+
j
ω
σ
y
σ
z
j
ω
ε
0
2
E
x
r
,
ω
+
1
γ
1
+
σ
x
j
ω
ε
0
J
s
r
,
ω
.However, this study is aimed at using finite difference time domain method in the study of transient electromagnetic field; hence, for time domain modeling (aimed at using Yee grid for effective time stepping), (12) was converted from frequency domain to time domain, leaving the source out. Using the following operations,(13a)
∇
×
F
x
r
,
ω
⟷
∇
×
F
x
r
,
t
,
(13b)
1
j
ω
∇
×
F
x
r
,
ω
⟷
∫
τ
=
-
∞
t
∇
×
F
x
r
,
τ
d
τ
,
(13c)
j
ω
F
x
r
,
ω
⟷
∂
F
x
r
,
t
∂
t
,
(13d)
1
j
ω
F
x
r
,
ω
⟷
∫
τ
=
-
∞
t
F
x
r
,
τ
d
τ
,where F
x represents x component of the fields, (12) can be written as(14)
1
γ
∇
×
H
x
r
,
t
+
σ
x
ε
0
γ
∫
τ
=
-
∞
t
∇
×
H
x
r
,
τ
d
τ
=
σ
b
r
γ
E
x
r
,
t
+
σ
b
r
σ
x
ε
0
γ
∫
τ
=
-
∞
t
E
x
r
,
τ
d
τ
+
∂
E
x
r
,
t
∂
t
+
σ
y
+
σ
z
ε
0
E
x
r
,
t
+
σ
y
σ
z
ε
0
2
∫
τ
=
-
∞
t
E
x
r
,
τ
d
τ
.Other components can be obtained likewise.For the magnetic field equation,(15)
∇
×
E
r
,
ω
=
-
j
ω
μ
0
s
r
H
r
,
ω
.Thus,(16)
∇
×
E
x
r
,
ω
=
-
j
ω
μ
0
1
+
σ
y
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
1
+
σ
x
j
ω
ε
0
-
1
H
x
r
,
ω
,
for
H
x
,
(17)
∇
×
E
y
r
,
ω
=
-
j
ω
μ
0
1
+
σ
x
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
1
+
σ
y
j
ω
ε
0
-
1
H
y
r
,
ω
,
for
H
y
,
(18)
∇
×
E
z
r
,
ω
=
-
j
ω
μ
0
1
+
σ
x
j
ω
ε
0
1
+
σ
y
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
-
1
H
z
r
,
ω
,
for
H
z
.Computing for H
x, (16) was multiplied by the UPML parameter S
x and then divided by magnetic permeability, -
μ
0, before expanding the operation/equation to obtain(19)
-
1
μ
0
∇
×
E
x
r
,
ω
-
σ
x
j
ω
μ
0
ε
0
∇
×
E
x
r
,
ω
=
j
ω
H
x
r
,
ω
+
j
ω
σ
y
+
σ
z
j
ω
ε
0
H
x
r
,
ω
+
j
ω
σ
y
σ
z
j
ω
ε
0
2
H
x
r
,
ω
.For time domain modeling (aimed at using Yee grid for effective time stepping), (19) was converted from frequency domain to time domain, leaving the source out. Using the same operations given in (13a)–(13d), (19) can be written as(20)
-
1
μ
0
∇
×
E
x
r
,
t
-
σ
x
μ
0
ε
0
∫
τ
=
-
∞
t
∇
×
E
x
r
,
τ
d
τ
=
∂
H
x
r
,
t
∂
t
+
σ
y
+
σ
z
ε
0
H
x
r
,
t
+
σ
y
σ
z
ε
0
2
∫
τ
=
-
∞
t
H
x
r
,
τ
d
τ
.Equations (14) and (20) give the general format for both electric field and magnetic field, respectively. Therefore, to determine the fields along other axes, only the index of both the fictitious and the background conductivities is changed.


### 2.3. Finite Difference Approximation

However, considering Yee’s grid system whereby the electric fieldE
i at integer time indices, 0
,
1
,
2
,
3
,
…
,
n
-
1, n, along the time axis changes with a uniform or time varying time step D
t
i, in a leap frog manner over the magnetic field H
i with the same increment but at intermediate time indices, hence, the integral operation is updated such that (14) and (20) are written as (21)
1
γ
∇
×
H
x
r
,
t
+
σ
x
ε
0
γ
∫
τ
=
0
t
+
D
t
/
2
∇
×
H
x
r
,
τ
d
τ
=
σ
b
r
γ
E
x
r
,
t
+
σ
b
r
σ
x
D
t
2
ε
0
γ
∫
τ
=
0
t
E
x
r
,
τ
d
τ
+
E
x
t
+
D
t
/
2
r
,
τ
+
∂
E
x
r
,
t
∂
t
+
σ
y
+
σ
z
ε
0
E
x
r
,
t
+
σ
y
σ
z
D
t
2
ε
0
2
∫
τ
=
0
t
E
x
r
,
τ
d
τ
+
E
x
t
+
D
t
/
2
r
,
τ
,
-
1
μ
0
∇
×
E
x
r
,
t
-
σ
x
μ
0
ε
0
∫
τ
=
0
t
∇
×
E
x
r
,
τ
d
τ
=
∂
H
x
r
,
t
∂
t
+
σ
y
+
σ
z
ε
0
H
x
r
,
t
+
σ
y
σ
z
D
t
2
ε
0
2
∫
τ
=
0
t
-
D
t
/
2
H
x
r
,
τ
d
τ
+
H
x
t
r
,
τ
.Using the forward difference first-order derivative scheme, (21) can be written as(22)
1
γ
∇
×
H
x
r
,
t
+
σ
x
Δ
t
ε
0
γ
∫
t
=
0
t
+
Δ
t
/
2
∇
×
H
x
r
,
t
=
σ
b
r
2
γ
E
x
r
,
t
+
Δ
t
+
E
x
r
,
t
+
σ
b
r
σ
x
Δ
t
4
ε
0
γ
∫
t
=
0
t
E
x
r
,
t
+
E
x
r
,
t
+
Δ
t
+
E
x
r
,
t
+
E
x
r
,
t
+
Δ
t
-
E
x
r
,
t
Δ
t
+
σ
y
+
σ
z
2
ε
0
E
x
r
,
t
+
Δ
t
+
E
x
r
,
t
+
σ
y
σ
z
Δ
t
4
ε
0
2
∫
t
=
0
t
E
x
r
,
t
+
E
x
r
,
t
+
Δ
t
+
E
x
r
,
t
,
-
1
μ
0
∇
×
E
x
r
,
t
-
σ
x
Δ
t
μ
0
ε
0
∫
t
=
0
t
∇
×
E
x
r
,
t
=
H
x
r
,
t
+
Δ
t
/
2
-
H
x
r
,
t
-
Δ
t
/
2
Δ
t
+
σ
y
+
σ
z
2
ε
0
H
x
r
,
t
+
Δ
t
2
+
H
x
r
,
t
-
Δ
t
2
+
σ
y
σ
z
Δ
t
4
ε
0
2
∫
t
=
0
t
-
Δ
t
/
2
H
x
r
,
t
-
Δ
t
2
+
H
x
r
,
t
+
Δ
t
2
+
H
x
r
,
t
-
Δ
t
2
.Furthermore, (22) can be presented in the form(23)
E
x
r
,
t
+
1
=
C
1
e
E
x
r
,
t
+
C
2
e
∇
×
H
x
r
,
t
+
Δ
t
2
+
C
3
e
∫
∇
×
H
x
r
,
t
+
Δ
t
2
+
C
4
e
∫
E
x
r
,
t
,
H
x
r
,
t
+
Δ
t
2
=
C
1
h
H
x
r
,
t
-
Δ
t
2
+
C
2
h
∇
×
E
x
r
,
t
+
C
3
h
∫
∇
×
E
x
r
,
t
+
C
4
h
∫
H
x
r
,
t
-
Δ
t
2
,where C
1
e
=
C
a
/
C
0,(24)
C
a
=
1
Δ
t
-
ε
0
σ
b
r
+
σ
y
+
σ
z
γ
2
ε
0
γ
-
ε
0
σ
b
r
σ
x
Δ
t
+
σ
y
σ
z
Δ
t
γ
4
ε
0
2
γ
,
C
0
=
1
Δ
t
+
ε
0
σ
b
r
+
σ
y
+
σ
z
γ
2
ε
0
γ
+
ε
0
σ
b
r
σ
x
Δ
t
+
σ
y
σ
z
Δ
t
γ
4
ε
0
2
γ
,
C
2
e
=
1
γ
C
0
,
C
3
e
=
σ
x
ε
0
γ
C
0
,
C
4
e
=
ε
0
σ
b
r
σ
x
Δ
t
-
σ
y
σ
z
Δ
t
γ
4
ε
0
2
γ
C
0
,while C
1
h
=
C
c
/
C
d,(25)
C
c
=
1
Δ
t
-
σ
y
+
σ
z
2
ε
0
-
σ
y
σ
z
Δ
t
4
ε
0
2
,
C
d
=
1
Δ
t
+
σ
y
+
σ
z
2
ε
0
+
σ
y
σ
z
Δ
t
4
ε
0
2
,
C
2
h
=
-
1
μ
0
C
d
,
C
3
h
=
-
Δ
t
σ
x
ε
0
μ
0
C
d
,
C
4
h
=
-
σ
y
σ
z
Δ
t
4
ε
0
2
C
d
.Equations (22) are the main equations for 3D UPML diffusion equation. Other components can be expressed by following the above steps. These can further be reduced to 2D by equating one of the axes to zero, giving(26)
1
γ
∇
×
H
x
y
,
z
,
t
=
σ
b
y
,
z
2
γ
E
x
y
,
z
,
t
+
1
+
E
x
y
,
z
,
t
+
E
x
y
,
z
,
t
+
1
-
E
x
y
,
z
,
t
Δ
t
+
σ
y
+
σ
z
2
ε
0
E
x
y
,
z
,
t
+
1
+
E
x
y
,
z
,
t
+
σ
y
σ
z
Δ
t
4
ε
0
2
∫
t
=
0
t
E
x
y
,
z
,
t
+
E
x
y
,
z
,
t
+
1
+
E
x
y
,
z
,
t
.It can also be reduced to 2.5D by changing one of the axes from space domain (Cartesian coordinate) to wave number domain and calculating the fields along a plane, at the surface of the transformed axis. The Dufort-Frankel explicit method [37] is easy to implement using the above equation, since the curl operator is independent of the PML operation. The obtained 3D equation can be converted from the use of regular grid to irregular grid without having a complete overhaul of the equation. That is, the curl operator can be solved independently. For x component of the electric field, the curl operator for an irregular grid is given as(27)
∇
×
H
x
y
,
z
,
t
=
Δ
y
i
+
1
Δ
y
i
-
E
z
i
-
1
,
j
n
+
Δ
y
i
Δ
y
i
-
E
z
i
+
1
,
j
n
-
Δ
z
j
+
1
Δ
z
j
-
E
y
i
,
j
-
1
n
+
Δ
z
j
Δ
z
j
-
E
y
i
,
j
+
1
n
,where i and j are integer indices, 1
,
2
,
3
,
…
,
m
-
1, m, along y and z directions, respectively. While the main medium’s conductivity is calculated using the five-point Dufort-Frankel explicit difference scheme [38],(28)
σ
i
,
j
-
=
N
i
,
j
σ
i
,
j
+
N
i
+
1
,
j
σ
i
+
1
,
j
+
N
i
,
j
+
1
σ
i
,
j
+
1
+
N
i
+
1
,
j
+
1
σ
i
+
1
,
j
+
1
,where i and j are the indices for y and z directions, respectively. Δ
y
i
-
=
(
Δ
y
i
+
Δ
y
i
+
1
)
/
2, Δ
z
j
-
=
(
Δ
z
j
+
Δ
z
j
+
1
)
/
2, N
i
,
j
=
Δ
i
,
j
/
Δ, Δ
i
,
j
=
Δ
z
i
Δ
y
j, and others are obtained likewise.


## 2.1. Approximations within the UPML Region

The essence of a UPML is to achieve zero reflection at the boundary of interest. Therefore, four (4) assumptions were made to ensure zero reflection at both sides of the boundary. They are as follows:(1)
The loss within the matched layer is set to be greater than zero; that is,k
>
0.
(2)
The impedance across the boundary is set as 1.
(3)
[μ
r] = [ε
r] = [s
r], having the matrix form s
a
0
0
0
s
b
0
0
0
s
c where s is the fictitious electrical and magnetic property within the matched layer. Thus, θ
1 = θ
2 and n
1 = n
2.
(4)
Consider a uniaxial isotropic medium wheres
a = s
b = 1
/
s
c and s
b
s
c
=
1 for [s
c].Therefore, the UPML parameters can be written as(1)
s
x
=
s
x
-
1
0
0
0
s
x
0
0
0
s
x
,
s
y
=
s
y
0
0
0
s
y
-
1
0
0
0
s
y
,
s
z
=
s
z
0
0
0
s
z
0
0
0
s
z
-
1
,
s
r
=
s
x
·
s
y
·
s
z
=
s
y
s
z
s
x
0
0
0
s
x
s
z
s
y
0
0
0
s
x
s
y
s
z
,where s
x, s
y, and s
z are arbitrary parameters calculated using the mathematical expressions given,(2)
S
x
=
1
+
σ
x
j
ω
ε
0
,
σ
x
=
ε
0
2
Δ
t
x
L
x
p
,
S
y
=
1
+
σ
y
j
ω
ε
0
,
σ
y
=
ε
0
2
Δ
t
y
L
y
p
,
S
z
=
1
+
σ
z
j
ω
ε
0
,
σ
z
=
ε
0
2
Δ
t
z
L
z
p
,assuming the k
m
a
x (maximum linear absorption factor) equals 1. Δ
t is the time step to be used, r
=
x, y, or z, and p is either 3 or 4, while r
/
L
r
p is a function grid space and it ensures that the value of the arbitrary parameter increases outwardly (Figure 1). Figure 1 shows an image of a grid whose PML region is represented by the black stripes that highlight its cells and its main medium is colored blue.Figure 1
3D spatial grid with perfectly matched layer.


## 2.2. Introducing the UPML into the EM Diffusion Equation

Maxwell’s equation for EM fields in linear, isotropic, and source-free media, under the quasistatic approximation in the frequency domain is given as follows:(3)
∇
×
E
r
,
ω
=
-
j
ω
μ
H
r
,
ω
,
(4)
∇
×
H
r
,
ω
=
σ
r
E
r
,
ω
,where E and H are the electric and the magnetic fields, respectively, and μ and σ are the physical properties, namely, magnetic permeability and electric conductivity, respectively. “r” represents the 3D Cartesian coordinate system (x, y, and z), while ω is the angular frequency (
2
π
f
). Although the conventional EM diffusion equation to be solved is parabolic (i.e., the conventional time domain equation neglects the displacement current), however, in this study we solved an EM hyperbolic equation (coupled Maxwell’s equation which simulates the direct interaction between the magnetic and electric field) for the FDTD algorithm by incorporating a fictitious displacement current to ensure an algorithm that is nonerroneously explicit and unconditionally stable [11, 13] for both 2D and 3D computation. This gives a wave-like equation in a lossy medium but by choosing the appropriate fictitious value for electric permittivity, an explicit, always stable, time stepping solution is obtained. To ensure the stability and convergence of a diffused EM wave initialized by an electric dipole with current density J
s, an arbitrary value of electric permittivity given as(5)
γ
=
3
μ
m
i
n
Δ
t
Δ
m
i
n
2was used to replace the real permittivity [11]. Also, the time step was calculated from the equation given in the following: (6)
Δ
t
≤
μ
m
i
n
σ
a
1
/
2
min
⁡
Δ
x
,
Δ
y
,
Δ
z
2
,where a = 2, 4, and 6 for 1D, 2D, and 3D, respectively.This keeps the velocity of the fictitious wave field slower than that simulated by the finite difference scheme and prevents the fictitious displacement current from dominating the diffusive EM field by constraining the length of the time step.Therefore, (4) becomes(7)
∇
×
H
r
,
ω
=
σ
b
r
E
r
,
ω
+
j
ω
γ
E
r
,
ω
+
J
s
r
,
ω
.By incorporating the UPML parameters (from (1) and (2)) into (7), we obtain (8)
∇
×
H
r
,
ω
=
σ
b
r
E
r
,
ω
+
j
ω
γ
s
r
E
r
,
ω
+
J
s
r
,
ω
.Thus,(9)
∇
×
H
x
r
,
ω
=
σ
b
r
E
x
r
,
ω
+
j
ω
γ
1
+
σ
y
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
1
+
σ
x
j
ω
ε
0
-
1
E
x
r
,
ω
+
J
s
r
,
ω
,
for
E
x
,
(10)
∇
×
H
y
r
,
ω
=
σ
b
r
E
y
r
,
ω
+
j
ω
γ
1
+
σ
x
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
1
+
σ
y
j
ω
ε
0
-
1
E
y
r
,
ω
+
J
s
r
,
ω
,
for
E
y
,
(11)
∇
×
H
z
r
,
ω
=
σ
b
r
E
z
r
,
ω
+
j
ω
γ
1
+
σ
x
j
ω
ε
0
1
+
σ
y
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
-
1
E
z
r
,
ω
+
J
s
r
,
ω
,
for
E
z
.Computing for E
x, (9) was multiplied by the UPML parameter S
x and then divided by fictitious permittivity, γ, before expanding the operation/equation to obtain(12)
1
γ
∇
×
H
x
r
,
ω
+
σ
x
j
ω
ε
0
γ
∇
×
H
x
r
,
ω
=
σ
b
r
γ
E
x
r
,
ω
+
σ
b
r
σ
x
j
ω
ε
0
γ
E
x
r
,
ω
+
j
ω
E
x
r
,
ω
+
j
ω
σ
y
+
σ
z
j
ω
ε
0
E
x
r
,
ω
+
j
ω
σ
y
σ
z
j
ω
ε
0
2
E
x
r
,
ω
+
1
γ
1
+
σ
x
j
ω
ε
0
J
s
r
,
ω
.However, this study is aimed at using finite difference time domain method in the study of transient electromagnetic field; hence, for time domain modeling (aimed at using Yee grid for effective time stepping), (12) was converted from frequency domain to time domain, leaving the source out. Using the following operations,(13a)
∇
×
F
x
r
,
ω
⟷
∇
×
F
x
r
,
t
,
(13b)
1
j
ω
∇
×
F
x
r
,
ω
⟷
∫
τ
=
-
∞
t
∇
×
F
x
r
,
τ
d
τ
,
(13c)
j
ω
F
x
r
,
ω
⟷
∂
F
x
r
,
t
∂
t
,
(13d)
1
j
ω
F
x
r
,
ω
⟷
∫
τ
=
-
∞
t
F
x
r
,
τ
d
τ
,where F
x represents x component of the fields, (12) can be written as(14)
1
γ
∇
×
H
x
r
,
t
+
σ
x
ε
0
γ
∫
τ
=
-
∞
t
∇
×
H
x
r
,
τ
d
τ
=
σ
b
r
γ
E
x
r
,
t
+
σ
b
r
σ
x
ε
0
γ
∫
τ
=
-
∞
t
E
x
r
,
τ
d
τ
+
∂
E
x
r
,
t
∂
t
+
σ
y
+
σ
z
ε
0
E
x
r
,
t
+
σ
y
σ
z
ε
0
2
∫
τ
=
-
∞
t
E
x
r
,
τ
d
τ
.Other components can be obtained likewise.For the magnetic field equation,(15)
∇
×
E
r
,
ω
=
-
j
ω
μ
0
s
r
H
r
,
ω
.Thus,(16)
∇
×
E
x
r
,
ω
=
-
j
ω
μ
0
1
+
σ
y
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
1
+
σ
x
j
ω
ε
0
-
1
H
x
r
,
ω
,
for
H
x
,
(17)
∇
×
E
y
r
,
ω
=
-
j
ω
μ
0
1
+
σ
x
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
1
+
σ
y
j
ω
ε
0
-
1
H
y
r
,
ω
,
for
H
y
,
(18)
∇
×
E
z
r
,
ω
=
-
j
ω
μ
0
1
+
σ
x
j
ω
ε
0
1
+
σ
y
j
ω
ε
0
1
+
σ
z
j
ω
ε
0
-
1
H
z
r
,
ω
,
for
H
z
.Computing for H
x, (16) was multiplied by the UPML parameter S
x and then divided by magnetic permeability, -
μ
0, before expanding the operation/equation to obtain(19)
-
1
μ
0
∇
×
E
x
r
,
ω
-
σ
x
j
ω
μ
0
ε
0
∇
×
E
x
r
,
ω
=
j
ω
H
x
r
,
ω
+
j
ω
σ
y
+
σ
z
j
ω
ε
0
H
x
r
,
ω
+
j
ω
σ
y
σ
z
j
ω
ε
0
2
H
x
r
,
ω
.For time domain modeling (aimed at using Yee grid for effective time stepping), (19) was converted from frequency domain to time domain, leaving the source out. Using the same operations given in (13a)–(13d), (19) can be written as(20)
-
1
μ
0
∇
×
E
x
r
,
t
-
σ
x
μ
0
ε
0
∫
τ
=
-
∞
t
∇
×
E
x
r
,
τ
d
τ
=
∂
H
x
r
,
t
∂
t
+
σ
y
+
σ
z
ε
0
H
x
r
,
t
+
σ
y
σ
z
ε
0
2
∫
τ
=
-
∞
t
H
x
r
,
τ
d
τ
.Equations (14) and (20) give the general format for both electric field and magnetic field, respectively. Therefore, to determine the fields along other axes, only the index of both the fictitious and the background conductivities is changed.


## 2.3. Finite Difference Approximation

However, considering Yee’s grid system whereby the electric fieldE
i at integer time indices, 0
,
1
,
2
,
3
,
…
,
n
-
1, n, along the time axis changes with a uniform or time varying time step D
t
i, in a leap frog manner over the magnetic field H
i with the same increment but at intermediate time indices, hence, the integral operation is updated such that (14) and (20) are written as (21)
1
γ
∇
×
H
x
r
,
t
+
σ
x
ε
0
γ
∫
τ
=
0
t
+
D
t
/
2
∇
×
H
x
r
,
τ
d
τ
=
σ
b
r
γ
E
x
r
,
t
+
σ
b
r
σ
x
D
t
2
ε
0
γ
∫
τ
=
0
t
E
x
r
,
τ
d
τ
+
E
x
t
+
D
t
/
2
r
,
τ
+
∂
E
x
r
,
t
∂
t
+
σ
y
+
σ
z
ε
0
E
x
r
,
t
+
σ
y
σ
z
D
t
2
ε
0
2
∫
τ
=
0
t
E
x
r
,
τ
d
τ
+
E
x
t
+
D
t
/
2
r
,
τ
,
-
1
μ
0
∇
×
E
x
r
,
t
-
σ
x
μ
0
ε
0
∫
τ
=
0
t
∇
×
E
x
r
,
τ
d
τ
=
∂
H
x
r
,
t
∂
t
+
σ
y
+
σ
z
ε
0
H
x
r
,
t
+
σ
y
σ
z
D
t
2
ε
0
2
∫
τ
=
0
t
-
D
t
/
2
H
x
r
,
τ
d
τ
+
H
x
t
r
,
τ
.Using the forward difference first-order derivative scheme, (21) can be written as(22)
1
γ
∇
×
H
x
r
,
t
+
σ
x
Δ
t
ε
0
γ
∫
t
=
0
t
+
Δ
t
/
2
∇
×
H
x
r
,
t
=
σ
b
r
2
γ
E
x
r
,
t
+
Δ
t
+
E
x
r
,
t
+
σ
b
r
σ
x
Δ
t
4
ε
0
γ
∫
t
=
0
t
E
x
r
,
t
+
E
x
r
,
t
+
Δ
t
+
E
x
r
,
t
+
E
x
r
,
t
+
Δ
t
-
E
x
r
,
t
Δ
t
+
σ
y
+
σ
z
2
ε
0
E
x
r
,
t
+
Δ
t
+
E
x
r
,
t
+
σ
y
σ
z
Δ
t
4
ε
0
2
∫
t
=
0
t
E
x
r
,
t
+
E
x
r
,
t
+
Δ
t
+
E
x
r
,
t
,
-
1
μ
0
∇
×
E
x
r
,
t
-
σ
x
Δ
t
μ
0
ε
0
∫
t
=
0
t
∇
×
E
x
r
,
t
=
H
x
r
,
t
+
Δ
t
/
2
-
H
x
r
,
t
-
Δ
t
/
2
Δ
t
+
σ
y
+
σ
z
2
ε
0
H
x
r
,
t
+
Δ
t
2
+
H
x
r
,
t
-
Δ
t
2
+
σ
y
σ
z
Δ
t
4
ε
0
2
∫
t
=
0
t
-
Δ
t
/
2
H
x
r
,
t
-
Δ
t
2
+
H
x
r
,
t
+
Δ
t
2
+
H
x
r
,
t
-
Δ
t
2
.Furthermore, (22) can be presented in the form(23)
E
x
r
,
t
+
1
=
C
1
e
E
x
r
,
t
+
C
2
e
∇
×
H
x
r
,
t
+
Δ
t
2
+
C
3
e
∫
∇
×
H
x
r
,
t
+
Δ
t
2
+
C
4
e
∫
E
x
r
,
t
,
H
x
r
,
t
+
Δ
t
2
=
C
1
h
H
x
r
,
t
-
Δ
t
2
+
C
2
h
∇
×
E
x
r
,
t
+
C
3
h
∫
∇
×
E
x
r
,
t
+
C
4
h
∫
H
x
r
,
t
-
Δ
t
2
,where C
1
e
=
C
a
/
C
0,(24)
C
a
=
1
Δ
t
-
ε
0
σ
b
r
+
σ
y
+
σ
z
γ
2
ε
0
γ
-
ε
0
σ
b
r
σ
x
Δ
t
+
σ
y
σ
z
Δ
t
γ
4
ε
0
2
γ
,
C
0
=
1
Δ
t
+
ε
0
σ
b
r
+
σ
y
+
σ
z
γ
2
ε
0
γ
+
ε
0
σ
b
r
σ
x
Δ
t
+
σ
y
σ
z
Δ
t
γ
4
ε
0
2
γ
,
C
2
e
=
1
γ
C
0
,
C
3
e
=
σ
x
ε
0
γ
C
0
,
C
4
e
=
ε
0
σ
b
r
σ
x
Δ
t
-
σ
y
σ
z
Δ
t
γ
4
ε
0
2
γ
C
0
,while C
1
h
=
C
c
/
C
d,(25)
C
c
=
1
Δ
t
-
σ
y
+
σ
z
2
ε
0
-
σ
y
σ
z
Δ
t
4
ε
0
2
,
C
d
=
1
Δ
t
+
σ
y
+
σ
z
2
ε
0
+
σ
y
σ
z
Δ
t
4
ε
0
2
,
C
2
h
=
-
1
μ
0
C
d
,
C
3
h
=
-
Δ
t
σ
x
ε
0
μ
0
C
d
,
C
4
h
=
-
σ
y
σ
z
Δ
t
4
ε
0
2
C
d
.Equations (22) are the main equations for 3D UPML diffusion equation. Other components can be expressed by following the above steps. These can further be reduced to 2D by equating one of the axes to zero, giving(26)
1
γ
∇
×
H
x
y
,
z
,
t
=
σ
b
y
,
z
2
γ
E
x
y
,
z
,
t
+
1
+
E
x
y
,
z
,
t
+
E
x
y
,
z
,
t
+
1
-
E
x
y
,
z
,
t
Δ
t
+
σ
y
+
σ
z
2
ε
0
E
x
y
,
z
,
t
+
1
+
E
x
y
,
z
,
t
+
σ
y
σ
z
Δ
t
4
ε
0
2
∫
t
=
0
t
E
x
y
,
z
,
t
+
E
x
y
,
z
,
t
+
1
+
E
x
y
,
z
,
t
.It can also be reduced to 2.5D by changing one of the axes from space domain (Cartesian coordinate) to wave number domain and calculating the fields along a plane, at the surface of the transformed axis. The Dufort-Frankel explicit method [37] is easy to implement using the above equation, since the curl operator is independent of the PML operation. The obtained 3D equation can be converted from the use of regular grid to irregular grid without having a complete overhaul of the equation. That is, the curl operator can be solved independently. For x component of the electric field, the curl operator for an irregular grid is given as(27)
∇
×
H
x
y
,
z
,
t
=
Δ
y
i
+
1
Δ
y
i
-
E
z
i
-
1
,
j
n
+
Δ
y
i
Δ
y
i
-
E
z
i
+
1
,
j
n
-
Δ
z
j
+
1
Δ
z
j
-
E
y
i
,
j
-
1
n
+
Δ
z
j
Δ
z
j
-
E
y
i
,
j
+
1
n
,where i and j are integer indices, 1
,
2
,
3
,
…
,
m
-
1, m, along y and z directions, respectively. While the main medium’s conductivity is calculated using the five-point Dufort-Frankel explicit difference scheme [38],(28)
σ
i
,
j
-
=
N
i
,
j
σ
i
,
j
+
N
i
+
1
,
j
σ
i
+
1
,
j
+
N
i
,
j
+
1
σ
i
,
j
+
1
+
N
i
+
1
,
j
+
1
σ
i
+
1
,
j
+
1
,where i and j are the indices for y and z directions, respectively. Δ
y
i
-
=
(
Δ
y
i
+
Δ
y
i
+
1
)
/
2, Δ
z
j
-
=
(
Δ
z
j
+
Δ
z
j
+
1
)
/
2, N
i
,
j
=
Δ
i
,
j
/
Δ, Δ
i
,
j
=
Δ
z
i
Δ
y
j, and others are obtained likewise.


## 3. Numerical Experiment

### 3.1. Varying the Number of Grids Used within the PML Region

In order to abate the reflection at the interface between two electrically contrasting media, the conductivitiesσ
x, σ
y, and σ
z increase gradually from a small value at the inner interface to a large value on the outer boundaries [15]. This means that the conductivity on the outer boundary is dependent on the grid spacing and the number of cells within the UPML region. Figure 2 shows the results of 2D modeling over a homogeneous medium (with resistivity of 100 Ωm), obtained by varying the number of cells within the PML region, thereby increasing the thickness of the absorbing area for the incident wave. The numbers of cells used were 5 × 5 and 15 × 15 of total cells 60 × 60 for both models. In practice, the condition of maintaining a lossy medium was ensured by assigning a background resistivity of 1 Ωm to the PML region. Hence, the conductivity of UPML region was kept within a reasonable range for absorption within the grid over a long period of time.Figure 2
2D model of a resistive medium with varying UPML size.
(a)
Electric field at 0.008 s using 5 cells for PML region
(b)
Electric field at 0.011 s using 5 cells for PML region
(c)
Electric field at 0.014 s using 5 cells for PML region
(d)
Electric field at 0.008 s using 15 cells for PML region
(e)
Electric field at 0.011 s using 15 cells for PML region
(f)
Electric field at 0.014 s using 15 cells for PML regionThe obtained results show that at early time between 0 s and 8.4 ms (<1500 iterations), the EM diffusion was uniform for both scenarios (Figures2(a) and 2(d)). However, at times equal to or greater than 11.2 ms (⩾2000 iterations), the effect of the matched layer was observed on the models (Figures 2(b) and 2(e)) as the inner boundaries act to absorb the incident signal with slight reduction in the maximum amplitude recorded on the 15 × 15 UPML cell’s model (Figure 2(e)). At larger times, 14 ms and above (⩾2500), continuous reduction in maximum amplitude was recorded on the 15 × 15 UPML cell’s model (Figure 2(f)) due to confinement of the main diffusion region (main medium). In general, the effectiveness of this method was independent on the number of cells within the UPML region (vis-à-vis the thickness or length of the PML space), as observed in Figures 2(a), 2(b), and 2(c), where there was no change in the diffusion pattern. Furthermore, the use of higher number of cells did not result in higher level of accuracy but a reduction in the number of cells within the main medium and the possibility of underestimating the amplitude at long iteration time. Therefore, the number of cells required within a UPML region for a model with 60 × 60 cells ranges between 5 × 5 cells for longer iteration period and 10 × 10 cells for an average number of iterations. Further discussion on the convergence of UPML is given in Section 3.3.2.


### 3.2. Comparing the PML Results with the Dirichlet Boundary Condition Result

Having established a base line for the UPML region in Section3.1, the Neumann boundary condition was used to test the accuracy of the UPML boundary treatment by comparing the FDTD results obtained using UPML boundary condition with the Neumann boundary condition’s results. The typical Neumann boundary condition used is that the directional derivative of the field normal to some boundary surface, termed the normal derivative, is zero. Figure 3 shows the result obtained from half-space of resistivity 100 Ωm, using the boundary conditions stated above. 10 cells were used for the UPML region (Figures 3(a)–3(c)).Figure 3
2D model of a resistive medium using UPML and Neumann boundary conditions.
(a)
Electric field at 0.005 s using UPML
(b)
Electric field at 0.007 s using UPML
(c)
Electric field at 0.012 s using UPML
(d)
Electric field at 0.005 s using Neumann boundary condition
(e)
Electric field at 0.007 s using Neumann boundary condition
(f)
Electric field at 0.012 s using Neumann boundary conditionIt presents the image of Ex field diffusion overx-z plane at 5 ms, 6.7 ms, and 12 ms. The EM diffusion over the given time using the Neumann BC shows a large disruption in the diffusion pattern, especially at the far edge of the grid due to fringe at the grids boundary (Figures 3(e) and 3(f)). The absence of such perturbation in the UPML sections (Figures 3(b) and 3(c)) shows the effectiveness of this UMPL method as an absorbing boundary condition, which is used in the FDTD algorithm. Hence, a 3D EM diffusion simulation was carried out over a homogeneous medium with resistivity 100 Ωm, using a 60 × 60 × 60 grid (Figure 4). 10 cells were used for the PML region on all sides and a dipole source in x direction was used with 900 iterations. Figure 4 shows the observed result at 5.5 ms, 8.2 ms, 11 ms, and 13.7 ms, respectively. This shows the sectional display of the EM field as it diffuses away from the dipole points (close to the zero origin across the x-section) in all directions. This is characterized by amplitude variation over a given location with time. The observed result is in agreement with the EM diffusion relation described by Liu et al. [38]. The EM field propagates outward as a function of time delay along the x-axis and diffuses in both y and z direction. From this figure, we can see that the transient field spreads upward, downward, and out from the dipole source with impedance along the zero section and few meters away.Figure 4
3D slice of electromagnetic field diffusion over homogeneous half-space at different time delays. (a) to (d) show the electric field at 5.5 ms, 8.2 ms, 11 ms, and 13.7 ms, respectively. The positions ofx-slice, y-slice, and z-slice used were 500 m, 30 m, and 30 m, respectively.
(a)
Electric field at 5.5 ms using 10 cells for PML region
(b)
Electric field at 8.2 ms using 10 cells for PML region
(c)
Electric field at 11 ms using 10 cells for PML region
(d)
Electric field at 13.7 ms using 10 cells for PML region


### 3.3. Modeling of MTEM Response over Homogeneous Medium

The relatively new method, multitransient electromagnetic method (MTEM), is a time domain controlled source electromagnetic method which uses transient current passed via grounded dipole electrodes to probe the subsurface geology. The returned time varying voltage response is measured simultaneously with multiple sets of dipole receivers, which are collinear with the current electrodes, and extends several kilometers (Figure5). The recorded response is a combination of the earth’s impulse response and the system response (combination of input current signal and receiver response), which can be measured by placing the same recording system used in impulse response measurement close to the input current [39]. However, in the absence of receiver response and uncorrelated noise, the recorded response can be expressed as(29)
V
r
,
t
=
Δ
x
s
Δ
x
r
i
r
,
t
∗
g
r
,
t
,where V
r
,
t is the measured voltage response, i
r
,
t is the measured input current, and g
(
r
,
t
) is the unknown earth impulse response.Figure 5
MTEM acquisition geometry.Simple processing methods such as the deconvolution of the measured voltage response for measured input current and air wave removal are used to retrieve the earth impulse response. The obtained earth impulse response is combined with the forward modeling response to obtain the true subsurface resistivity distribution over the volume of the earth via inversion process, with the focal point being midway between the source and the receiver.


#### 3.3.1.
2D Modeling

Using the FDTD algorithm derived above, on a 61 × 41 regular grid medium, with 100 m separation and 10 cells used for PML regions set at the far end of each axis, the MTEM earth response over a homogeneous subsurface material with resistivity of 10 Ωm was modeled. Figure6 shows the MTEM impulse response obtained at offsets ranging from 1000 m to 2200 m at an interval of 400 m using a Gaussian source. The FDTD responses (in red) were compared with the analytical results (in blue) obtained from Wilson’s [40] earth response equation. The obtained FDTD responses were in good agreement with the analytical responses; that is, the error margins observed were within the range of tolerance.Figure 6
MTEM earth impulse responses at offsets: (a) 1000 m; (b) 1400 m; (c) 1800 m; and (d) 2200 m.
(a)
(b)
(c)
(d)


#### 3.3.2. Numerical Study of UPML Convergence

In this section, we present numerical examples to illustrate the effectiveness of the above FDTD algorithm, vis-à-vis, the performance of the UPML method. FDTD simulation of MTEM earth impulse response was performed using different grid setup over a uniform earth material with resistivity of 10 Ωm at various offsets. The results were compared with the exact impulse responsew
0 for the respective offsets, obtained using the analytical equation [40], and the obtained relative error w
-
w
0
/
w
0 was plotted against time, while the root mean square error and root mean square relative error were plotted against receiver offset. This illustrates the general procedures for quantitative study of the convergence properties of the method. Consequently, the usefulness of the FDTD procedure for computing the MTEM earth impulse response was sufficiently verified by the comparison between the numerical result and the analytical data.Firstly, using a 76 × 56 regularly gridded medium, with 100 m separation and 10 cells used for UPML regions set at the far end of each axis, the MTEM earth response over a homogeneous subsurface material with resistivity of 10 Ωm was modeled. Figure7 shows the MTEM impulse response error obtained at offsets ranging from 1000 m to 4100 m at an interval of 300 m. The error margin varies between −0.18 and 0.16 at early time and late time, except over few centiseconds at the start of each response curve. The maximum error margin was recorded at offset 1000 m and it decreases as the offset increases. High level of convergence was recorded at time ranging between 0.05 s and 0.4 s at most of the offsets. This indicates that for different choices of offset, the numerical and analytical convergence is time dependent, more accurate over the peak time, and less accurate over few centiseconds at the start of each response.Figure 7
Convergence of the FDTD-computed MTEM earth impulse response error versus time for 2D homogeneous medium withρ = 10 Ωm.Furthermore, the performance of the UPML method as well as the effect of changes in the size of both the UPML region and the main medium on computed response was established by observing the numerical and analytical data convergence with offset. Using a 76 × 56 regularly gridded medium with 100 m separation, the MTEM earth response over a homogeneous subsurface material with resistivity of 10 Ωm was modeled. Two separate tests were conducted: (1) the number of cells used for the main medium was kept constant while the number of cells used for the UPML region was varied and (2) the number of cells used for the main medium was reduced as the number of cells used for the UPML region was increased. The root mean square (RMS) error and the root mean square relative (RMSR) error (which measures the absolute error as a ratio of absolute value rather than per period. As a result, it eliminates the problem of interpreting the measure of accuracy relative to the magnitude of the response values per offset) were plotted against the receiver offsets for each scenario to make comparison of the numerical response’s convergence over several offsets, ranging from 1000 m to 4000 m at an interval of 100 m.Figure8 shows the result obtained from test 1. 52 × 32 cells were used for the main medium while the number of cells used for the UPML region ranged from 5 × 5 to 25 × 25, using an increment of 5. The RMS error curves (Figure 8(a)) show no significant difference in accuracy of the numerical response with increase in the size of the UPML region. However, the exponential decay observed with offset depicts the change in error with respect to the magnitude of the response at each offset thereby showing the conventional decay in amplitude with offset. On the other hand, the RMSR error plot (Figure 8(b)) shows that at short offsets, the numerical response converges more as the UPML region increases (the error reduces with increase in PML size). However, this phenomenon changes over the large offsets where a maximum error of 0.22 was recorded when the UPML size is 25 × 25 (more than half the size of the main medium). Hence, the numerical results are expected to remain in good agreement with the exact results over a long offset provided the UPML region is less than half of the main medium.Figure 8
Convergence plot obtained by varying the size (L) of the UPML region at various offsets.
(a)
(b)Figure9 shows the result obtained from test 2, where the size of the main medium ranges from 70 × 50 to 40 × 20 as the number of cells used for the UPML region increases from 5 to 35 cells, respectively. The total number of cells used over the two sides is given as 76 × 56. The RMS error curves (Figure 9(a)) show a significant increase in value as the number of UPML cells increases from 30 to 35 and the main medium becomes smaller. This significant increase was also observed on the RMSR error plot (Figure 9(b)) where the value also increased with offset. This phenomenon (decrease in accuracy of the numerical response with increase in the size of the UPML region) was also observed when L equals 25 and 30. On the other hand, the RMSR error plot (Figure 9(b)) also shows that the numerical response converges more as the number of cells within the UPML region increases from 5 to 20, with the minimum error of 0.02 recorded at L = 20 (less than half the size of the main medium).Figure 9
Convergence plot obtained by varying the size (L) of the UPML region and the main medium at various offsets.
(a)
(b)Considering the results obtained above, the above FDTD algorithm has been found to be flexible and our numerical results are in good agreement with the exact results. For example, the numerical result in Figure9(b) shows that an accuracy of 0.03 with the exact result can be realized over a wide range of offsets. This illustrates that accurate results for a specified range of offset of interest can be guaranteed. Hence, it is accurate for the applications considered. Note that several factors affect the accuracy of numerical results including the total simulation time, overall grid size, and grid resolution. Hence, a quick study of the best UPML size for the medium in question is advised, preferably a value less than half the size of the main medium with few cells between them.


#### 3.3.3.
3D Modeling

With the above confirmation, together with results from several other source-receiver positions, 3D modeling of MTEM earth impulse response over a range of offsets and resistivity were carried out on a 61 × 61 × 41 grid (with a PRBS source) using the derived algorithm and the results were plotted as Figures10 and 11, respectively.Figure 10
MTEM earth response at different offsets.Figure 11
MTEM earth response over subsurface material with different resistivities.
(a)
Earth impulse response plotted separately on an arithmetic scale
(b)
Plot of earth impulse response on a log-log scaleFigure10 shows the decreasing nature of the earth response amplitude and an increase in peak time with offset over a 100 Ωm homogeneous medium. Therefore, at very large offset, the MTEM earth response becomes so small and negligible due to its rapid attenuation rate. The observed peak time ranges from 0.00122 s for an offset of 1000 m with amplitude of 5.6e
−6 V/A/m2/s to 0.005 s for an offset of 2000 m with amplitude of 1.33e
−7 V/A/m2/s (Figure 10). On the other hand, the earth response amplitude and peak time decrease and increase, respectively, with decrease in the subsurface resistivity (Figure 11(b)). The observed peak time ranges from 0.000121 s for resistivity of 1000 Ωm with amplitude of 4.9e
−4 V/A/m2/s to 0.012 s for resistivity of 10 Ωm with amplitude of 6.38e
−8 V/A/m2/s (Figure 11(a)).


## 3.1. Varying the Number of Grids Used within the PML Region

In order to abate the reflection at the interface between two electrically contrasting media, the conductivitiesσ
x, σ
y, and σ
z increase gradually from a small value at the inner interface to a large value on the outer boundaries [15]. This means that the conductivity on the outer boundary is dependent on the grid spacing and the number of cells within the UPML region. Figure 2 shows the results of 2D modeling over a homogeneous medium (with resistivity of 100 Ωm), obtained by varying the number of cells within the PML region, thereby increasing the thickness of the absorbing area for the incident wave. The numbers of cells used were 5 × 5 and 15 × 15 of total cells 60 × 60 for both models. In practice, the condition of maintaining a lossy medium was ensured by assigning a background resistivity of 1 Ωm to the PML region. Hence, the conductivity of UPML region was kept within a reasonable range for absorption within the grid over a long period of time.Figure 2
2D model of a resistive medium with varying UPML size.
(a)
Electric field at 0.008 s using 5 cells for PML region
(b)
Electric field at 0.011 s using 5 cells for PML region
(c)
Electric field at 0.014 s using 5 cells for PML region
(d)
Electric field at 0.008 s using 15 cells for PML region
(e)
Electric field at 0.011 s using 15 cells for PML region
(f)
Electric field at 0.014 s using 15 cells for PML regionThe obtained results show that at early time between 0 s and 8.4 ms (<1500 iterations), the EM diffusion was uniform for both scenarios (Figures2(a) and 2(d)). However, at times equal to or greater than 11.2 ms (⩾2000 iterations), the effect of the matched layer was observed on the models (Figures 2(b) and 2(e)) as the inner boundaries act to absorb the incident signal with slight reduction in the maximum amplitude recorded on the 15 × 15 UPML cell’s model (Figure 2(e)). At larger times, 14 ms and above (⩾2500), continuous reduction in maximum amplitude was recorded on the 15 × 15 UPML cell’s model (Figure 2(f)) due to confinement of the main diffusion region (main medium). In general, the effectiveness of this method was independent on the number of cells within the UPML region (vis-à-vis the thickness or length of the PML space), as observed in Figures 2(a), 2(b), and 2(c), where there was no change in the diffusion pattern. Furthermore, the use of higher number of cells did not result in higher level of accuracy but a reduction in the number of cells within the main medium and the possibility of underestimating the amplitude at long iteration time. Therefore, the number of cells required within a UPML region for a model with 60 × 60 cells ranges between 5 × 5 cells for longer iteration period and 10 × 10 cells for an average number of iterations. Further discussion on the convergence of UPML is given in Section 3.3.2.


## 3.2. Comparing the PML Results with the Dirichlet Boundary Condition Result

Having established a base line for the UPML region in Section3.1, the Neumann boundary condition was used to test the accuracy of the UPML boundary treatment by comparing the FDTD results obtained using UPML boundary condition with the Neumann boundary condition’s results. The typical Neumann boundary condition used is that the directional derivative of the field normal to some boundary surface, termed the normal derivative, is zero. Figure 3 shows the result obtained from half-space of resistivity 100 Ωm, using the boundary conditions stated above. 10 cells were used for the UPML region (Figures 3(a)–3(c)).Figure 3
2D model of a resistive medium using UPML and Neumann boundary conditions.
(a)
Electric field at 0.005 s using UPML
(b)
Electric field at 0.007 s using UPML
(c)
Electric field at 0.012 s using UPML
(d)
Electric field at 0.005 s using Neumann boundary condition
(e)
Electric field at 0.007 s using Neumann boundary condition
(f)
Electric field at 0.012 s using Neumann boundary conditionIt presents the image of Ex field diffusion overx-z plane at 5 ms, 6.7 ms, and 12 ms. The EM diffusion over the given time using the Neumann BC shows a large disruption in the diffusion pattern, especially at the far edge of the grid due to fringe at the grids boundary (Figures 3(e) and 3(f)). The absence of such perturbation in the UPML sections (Figures 3(b) and 3(c)) shows the effectiveness of this UMPL method as an absorbing boundary condition, which is used in the FDTD algorithm. Hence, a 3D EM diffusion simulation was carried out over a homogeneous medium with resistivity 100 Ωm, using a 60 × 60 × 60 grid (Figure 4). 10 cells were used for the PML region on all sides and a dipole source in x direction was used with 900 iterations. Figure 4 shows the observed result at 5.5 ms, 8.2 ms, 11 ms, and 13.7 ms, respectively. This shows the sectional display of the EM field as it diffuses away from the dipole points (close to the zero origin across the x-section) in all directions. This is characterized by amplitude variation over a given location with time. The observed result is in agreement with the EM diffusion relation described by Liu et al. [38]. The EM field propagates outward as a function of time delay along the x-axis and diffuses in both y and z direction. From this figure, we can see that the transient field spreads upward, downward, and out from the dipole source with impedance along the zero section and few meters away.Figure 4
3D slice of electromagnetic field diffusion over homogeneous half-space at different time delays. (a) to (d) show the electric field at 5.5 ms, 8.2 ms, 11 ms, and 13.7 ms, respectively. The positions ofx-slice, y-slice, and z-slice used were 500 m, 30 m, and 30 m, respectively.
(a)
Electric field at 5.5 ms using 10 cells for PML region
(b)
Electric field at 8.2 ms using 10 cells for PML region
(c)
Electric field at 11 ms using 10 cells for PML region
(d)
Electric field at 13.7 ms using 10 cells for PML region


## 3.3. Modeling of MTEM Response over Homogeneous Medium

The relatively new method, multitransient electromagnetic method (MTEM), is a time domain controlled source electromagnetic method which uses transient current passed via grounded dipole electrodes to probe the subsurface geology. The returned time varying voltage response is measured simultaneously with multiple sets of dipole receivers, which are collinear with the current electrodes, and extends several kilometers (Figure5). The recorded response is a combination of the earth’s impulse response and the system response (combination of input current signal and receiver response), which can be measured by placing the same recording system used in impulse response measurement close to the input current [39]. However, in the absence of receiver response and uncorrelated noise, the recorded response can be expressed as(29)
V
r
,
t
=
Δ
x
s
Δ
x
r
i
r
,
t
∗
g
r
,
t
,where V
r
,
t is the measured voltage response, i
r
,
t is the measured input current, and g
(
r
,
t
) is the unknown earth impulse response.Figure 5
MTEM acquisition geometry.Simple processing methods such as the deconvolution of the measured voltage response for measured input current and air wave removal are used to retrieve the earth impulse response. The obtained earth impulse response is combined with the forward modeling response to obtain the true subsurface resistivity distribution over the volume of the earth via inversion process, with the focal point being midway between the source and the receiver.


### 3.3.1.
2D Modeling

Using the FDTD algorithm derived above, on a 61 × 41 regular grid medium, with 100 m separation and 10 cells used for PML regions set at the far end of each axis, the MTEM earth response over a homogeneous subsurface material with resistivity of 10 Ωm was modeled. Figure6 shows the MTEM impulse response obtained at offsets ranging from 1000 m to 2200 m at an interval of 400 m using a Gaussian source. The FDTD responses (in red) were compared with the analytical results (in blue) obtained from Wilson’s [40] earth response equation. The obtained FDTD responses were in good agreement with the analytical responses; that is, the error margins observed were within the range of tolerance.Figure 6
MTEM earth impulse responses at offsets: (a) 1000 m; (b) 1400 m; (c) 1800 m; and (d) 2200 m.
(a)
(b)
(c)
(d)


### 3.3.2. Numerical Study of UPML Convergence

In this section, we present numerical examples to illustrate the effectiveness of the above FDTD algorithm, vis-à-vis, the performance of the UPML method. FDTD simulation of MTEM earth impulse response was performed using different grid setup over a uniform earth material with resistivity of 10 Ωm at various offsets. The results were compared with the exact impulse responsew
0 for the respective offsets, obtained using the analytical equation [40], and the obtained relative error w
-
w
0
/
w
0 was plotted against time, while the root mean square error and root mean square relative error were plotted against receiver offset. This illustrates the general procedures for quantitative study of the convergence properties of the method. Consequently, the usefulness of the FDTD procedure for computing the MTEM earth impulse response was sufficiently verified by the comparison between the numerical result and the analytical data.Firstly, using a 76 × 56 regularly gridded medium, with 100 m separation and 10 cells used for UPML regions set at the far end of each axis, the MTEM earth response over a homogeneous subsurface material with resistivity of 10 Ωm was modeled. Figure7 shows the MTEM impulse response error obtained at offsets ranging from 1000 m to 4100 m at an interval of 300 m. The error margin varies between −0.18 and 0.16 at early time and late time, except over few centiseconds at the start of each response curve. The maximum error margin was recorded at offset 1000 m and it decreases as the offset increases. High level of convergence was recorded at time ranging between 0.05 s and 0.4 s at most of the offsets. This indicates that for different choices of offset, the numerical and analytical convergence is time dependent, more accurate over the peak time, and less accurate over few centiseconds at the start of each response.Figure 7
Convergence of the FDTD-computed MTEM earth impulse response error versus time for 2D homogeneous medium withρ = 10 Ωm.Furthermore, the performance of the UPML method as well as the effect of changes in the size of both the UPML region and the main medium on computed response was established by observing the numerical and analytical data convergence with offset. Using a 76 × 56 regularly gridded medium with 100 m separation, the MTEM earth response over a homogeneous subsurface material with resistivity of 10 Ωm was modeled. Two separate tests were conducted: (1) the number of cells used for the main medium was kept constant while the number of cells used for the UPML region was varied and (2) the number of cells used for the main medium was reduced as the number of cells used for the UPML region was increased. The root mean square (RMS) error and the root mean square relative (RMSR) error (which measures the absolute error as a ratio of absolute value rather than per period. As a result, it eliminates the problem of interpreting the measure of accuracy relative to the magnitude of the response values per offset) were plotted against the receiver offsets for each scenario to make comparison of the numerical response’s convergence over several offsets, ranging from 1000 m to 4000 m at an interval of 100 m.Figure8 shows the result obtained from test 1. 52 × 32 cells were used for the main medium while the number of cells used for the UPML region ranged from 5 × 5 to 25 × 25, using an increment of 5. The RMS error curves (Figure 8(a)) show no significant difference in accuracy of the numerical response with increase in the size of the UPML region. However, the exponential decay observed with offset depicts the change in error with respect to the magnitude of the response at each offset thereby showing the conventional decay in amplitude with offset. On the other hand, the RMSR error plot (Figure 8(b)) shows that at short offsets, the numerical response converges more as the UPML region increases (the error reduces with increase in PML size). However, this phenomenon changes over the large offsets where a maximum error of 0.22 was recorded when the UPML size is 25 × 25 (more than half the size of the main medium). Hence, the numerical results are expected to remain in good agreement with the exact results over a long offset provided the UPML region is less than half of the main medium.Figure 8
Convergence plot obtained by varying the size (L) of the UPML region at various offsets.
(a)
(b)Figure9 shows the result obtained from test 2, where the size of the main medium ranges from 70 × 50 to 40 × 20 as the number of cells used for the UPML region increases from 5 to 35 cells, respectively. The total number of cells used over the two sides is given as 76 × 56. The RMS error curves (Figure 9(a)) show a significant increase in value as the number of UPML cells increases from 30 to 35 and the main medium becomes smaller. This significant increase was also observed on the RMSR error plot (Figure 9(b)) where the value also increased with offset. This phenomenon (decrease in accuracy of the numerical response with increase in the size of the UPML region) was also observed when L equals 25 and 30. On the other hand, the RMSR error plot (Figure 9(b)) also shows that the numerical response converges more as the number of cells within the UPML region increases from 5 to 20, with the minimum error of 0.02 recorded at L = 20 (less than half the size of the main medium).Figure 9
Convergence plot obtained by varying the size (L) of the UPML region and the main medium at various offsets.
(a)
(b)Considering the results obtained above, the above FDTD algorithm has been found to be flexible and our numerical results are in good agreement with the exact results. For example, the numerical result in Figure9(b) shows that an accuracy of 0.03 with the exact result can be realized over a wide range of offsets. This illustrates that accurate results for a specified range of offset of interest can be guaranteed. Hence, it is accurate for the applications considered. Note that several factors affect the accuracy of numerical results including the total simulation time, overall grid size, and grid resolution. Hence, a quick study of the best UPML size for the medium in question is advised, preferably a value less than half the size of the main medium with few cells between them.


### 3.3.3.
3D Modeling

With the above confirmation, together with results from several other source-receiver positions, 3D modeling of MTEM earth impulse response over a range of offsets and resistivity were carried out on a 61 × 61 × 41 grid (with a PRBS source) using the derived algorithm and the results were plotted as Figures10 and 11, respectively.Figure 10
MTEM earth response at different offsets.Figure 11
MTEM earth response over subsurface material with different resistivities.
(a)
Earth impulse response plotted separately on an arithmetic scale
(b)
Plot of earth impulse response on a log-log scaleFigure10 shows the decreasing nature of the earth response amplitude and an increase in peak time with offset over a 100 Ωm homogeneous medium. Therefore, at very large offset, the MTEM earth response becomes so small and negligible due to its rapid attenuation rate. The observed peak time ranges from 0.00122 s for an offset of 1000 m with amplitude of 5.6e
−6 V/A/m2/s to 0.005 s for an offset of 2000 m with amplitude of 1.33e
−7 V/A/m2/s (Figure 10). On the other hand, the earth response amplitude and peak time decrease and increase, respectively, with decrease in the subsurface resistivity (Figure 11(b)). The observed peak time ranges from 0.000121 s for resistivity of 1000 Ωm with amplitude of 4.9e
−4 V/A/m2/s to 0.012 s for resistivity of 10 Ωm with amplitude of 6.38e
−8 V/A/m2/s (Figure 11(a)).


## 3.3.1.
2D Modeling

Using the FDTD algorithm derived above, on a 61 × 41 regular grid medium, with 100 m separation and 10 cells used for PML regions set at the far end of each axis, the MTEM earth response over a homogeneous subsurface material with resistivity of 10 Ωm was modeled. Figure6 shows the MTEM impulse response obtained at offsets ranging from 1000 m to 2200 m at an interval of 400 m using a Gaussian source. The FDTD responses (in red) were compared with the analytical results (in blue) obtained from Wilson’s [40] earth response equation. The obtained FDTD responses were in good agreement with the analytical responses; that is, the error margins observed were within the range of tolerance.Figure 6
MTEM earth impulse responses at offsets: (a) 1000 m; (b) 1400 m; (c) 1800 m; and (d) 2200 m.
(a)
(b)
(c)
(d)


## 3.3.2. Numerical Study of UPML Convergence

In this section, we present numerical examples to illustrate the effectiveness of the above FDTD algorithm, vis-à-vis, the performance of the UPML method. FDTD simulation of MTEM earth impulse response was performed using different grid setup over a uniform earth material with resistivity of 10 Ωm at various offsets. The results were compared with the exact impulse responsew
0 for the respective offsets, obtained using the analytical equation [40], and the obtained relative error w
-
w
0
/
w
0 was plotted against time, while the root mean square error and root mean square relative error were plotted against receiver offset. This illustrates the general procedures for quantitative study of the convergence properties of the method. Consequently, the usefulness of the FDTD procedure for computing the MTEM earth impulse response was sufficiently verified by the comparison between the numerical result and the analytical data.Firstly, using a 76 × 56 regularly gridded medium, with 100 m separation and 10 cells used for UPML regions set at the far end of each axis, the MTEM earth response over a homogeneous subsurface material with resistivity of 10 Ωm was modeled. Figure7 shows the MTEM impulse response error obtained at offsets ranging from 1000 m to 4100 m at an interval of 300 m. The error margin varies between −0.18 and 0.16 at early time and late time, except over few centiseconds at the start of each response curve. The maximum error margin was recorded at offset 1000 m and it decreases as the offset increases. High level of convergence was recorded at time ranging between 0.05 s and 0.4 s at most of the offsets. This indicates that for different choices of offset, the numerical and analytical convergence is time dependent, more accurate over the peak time, and less accurate over few centiseconds at the start of each response.Figure 7
Convergence of the FDTD-computed MTEM earth impulse response error versus time for 2D homogeneous medium withρ = 10 Ωm.Furthermore, the performance of the UPML method as well as the effect of changes in the size of both the UPML region and the main medium on computed response was established by observing the numerical and analytical data convergence with offset. Using a 76 × 56 regularly gridded medium with 100 m separation, the MTEM earth response over a homogeneous subsurface material with resistivity of 10 Ωm was modeled. Two separate tests were conducted: (1) the number of cells used for the main medium was kept constant while the number of cells used for the UPML region was varied and (2) the number of cells used for the main medium was reduced as the number of cells used for the UPML region was increased. The root mean square (RMS) error and the root mean square relative (RMSR) error (which measures the absolute error as a ratio of absolute value rather than per period. As a result, it eliminates the problem of interpreting the measure of accuracy relative to the magnitude of the response values per offset) were plotted against the receiver offsets for each scenario to make comparison of the numerical response’s convergence over several offsets, ranging from 1000 m to 4000 m at an interval of 100 m.Figure8 shows the result obtained from test 1. 52 × 32 cells were used for the main medium while the number of cells used for the UPML region ranged from 5 × 5 to 25 × 25, using an increment of 5. The RMS error curves (Figure 8(a)) show no significant difference in accuracy of the numerical response with increase in the size of the UPML region. However, the exponential decay observed with offset depicts the change in error with respect to the magnitude of the response at each offset thereby showing the conventional decay in amplitude with offset. On the other hand, the RMSR error plot (Figure 8(b)) shows that at short offsets, the numerical response converges more as the UPML region increases (the error reduces with increase in PML size). However, this phenomenon changes over the large offsets where a maximum error of 0.22 was recorded when the UPML size is 25 × 25 (more than half the size of the main medium). Hence, the numerical results are expected to remain in good agreement with the exact results over a long offset provided the UPML region is less than half of the main medium.Figure 8
Convergence plot obtained by varying the size (L) of the UPML region at various offsets.
(a)
(b)Figure9 shows the result obtained from test 2, where the size of the main medium ranges from 70 × 50 to 40 × 20 as the number of cells used for the UPML region increases from 5 to 35 cells, respectively. The total number of cells used over the two sides is given as 76 × 56. The RMS error curves (Figure 9(a)) show a significant increase in value as the number of UPML cells increases from 30 to 35 and the main medium becomes smaller. This significant increase was also observed on the RMSR error plot (Figure 9(b)) where the value also increased with offset. This phenomenon (decrease in accuracy of the numerical response with increase in the size of the UPML region) was also observed when L equals 25 and 30. On the other hand, the RMSR error plot (Figure 9(b)) also shows that the numerical response converges more as the number of cells within the UPML region increases from 5 to 20, with the minimum error of 0.02 recorded at L = 20 (less than half the size of the main medium).Figure 9
Convergence plot obtained by varying the size (L) of the UPML region and the main medium at various offsets.
(a)
(b)Considering the results obtained above, the above FDTD algorithm has been found to be flexible and our numerical results are in good agreement with the exact results. For example, the numerical result in Figure9(b) shows that an accuracy of 0.03 with the exact result can be realized over a wide range of offsets. This illustrates that accurate results for a specified range of offset of interest can be guaranteed. Hence, it is accurate for the applications considered. Note that several factors affect the accuracy of numerical results including the total simulation time, overall grid size, and grid resolution. Hence, a quick study of the best UPML size for the medium in question is advised, preferably a value less than half the size of the main medium with few cells between them.


## 3.3.3.
3D Modeling

With the above confirmation, together with results from several other source-receiver positions, 3D modeling of MTEM earth impulse response over a range of offsets and resistivity were carried out on a 61 × 61 × 41 grid (with a PRBS source) using the derived algorithm and the results were plotted as Figures10 and 11, respectively.Figure 10
MTEM earth response at different offsets.Figure 11
MTEM earth response over subsurface material with different resistivities.
(a)
Earth impulse response plotted separately on an arithmetic scale
(b)
Plot of earth impulse response on a log-log scaleFigure10 shows the decreasing nature of the earth response amplitude and an increase in peak time with offset over a 100 Ωm homogeneous medium. Therefore, at very large offset, the MTEM earth response becomes so small and negligible due to its rapid attenuation rate. The observed peak time ranges from 0.00122 s for an offset of 1000 m with amplitude of 5.6e
−6 V/A/m2/s to 0.005 s for an offset of 2000 m with amplitude of 1.33e
−7 V/A/m2/s (Figure 10). On the other hand, the earth response amplitude and peak time decrease and increase, respectively, with decrease in the subsurface resistivity (Figure 11(b)). The observed peak time ranges from 0.000121 s for resistivity of 1000 Ωm with amplitude of 4.9e
−4 V/A/m2/s to 0.012 s for resistivity of 10 Ωm with amplitude of 6.38e
−8 V/A/m2/s (Figure 11(a)).


## 4.
2D MTEM Response Modeling of South China Shale Gas Earth Model

Prior to any field observation using this technique, we decided to model the MTEM response (using the above algorithm) over an existing shale gas model from South China and observe the correlation between the obtained apparent resistivity and the model. The shale gas model was built on a 81 × 42 grid space (Figure12), using a grid spacing of 300 m along the x direction and 150 m along the z direction, 2500 time steps, and an offset varying with respect to the receiver position. 10 cells were used for both the far left and lower boundary PML region. Also, 2 extra cells were used to separate the PML regions from the main grid space. In addition, the upward continuation technique by Wang and Hohmann [11] was applied at the boundary between the earth surface and the air. A highly conductive layer with resistivity of 30 ohm m and maximum thickness of 350 m, confined by a layer with resistivity of 300 ohm m at the bottom and a thick, highly resistive layer with resistivity of 2000 ohm m, was used to represent a confined shale gas reserve. The highly resistive layer with an average thickness of 900 m runs from the surface to the top of the shale gas reserve. The shale layer’s structure shows an evident faulting process within the area, which led to the trapping of a section of the shale gas reserve within the very high resistive layer between the distances 7 km and 18 km. The pseudorandom binary waveform with order n
=
10, having three cycles, was used to generate the short dipole current source used for the modeling. Considering the large variation in resistivity distribution over the section, the time step was calculated from the equation given in the following:(30)
Δ
t
=
c
∗
μ
m
i
n
σ
a
1
/
2
min
⁡
Δ
x
,
Δ
z
2
,where c
=
n
∗
(
σ
b
/
σ
m
i
n
) and n ranges between 0.15 and 0.5, to keep the velocity of the EM field gradational with change in Δ
t at the near surface (to reduce the effect of the high resistivity) and across the resistivity distribution, as against the conventional finite difference scheme. The resistivity ratio ensures the record of wide range of response time.Figure 12
2D earth model for South China shale gas.Figure13 shows the unprocessed output data at offsets of 600 m, 2400 m, 4200 m, and 5400 m, with the midpoint of the current dipole (Tx) at 4200 m. The MTEM responses were presented as offset-time profile plots (Figure 14), 2D common offset sections (Figure 15), and 2D apparent resistivity section (Figure 16).Figure 13
Unprocessed earth response plot against time at offsets 600 m, 2400 m, 4200 m, and 5400 m with Tx at 4200 m. The apparent resistivity values were calculated from the peak time obtained from each trace [40]. The obtained impulse response over a range of offsets at Tx 3000 m and 9000 m is shown in Figure 14.Figure 14
Processed earth impulse response plot at offsets ranging from 600 m to 9000 m at an interval of 600 m, with Tx at (a) 3000 m and (b) 9000 m.
(a)
(b)Figure 15
2D MTEM response displayed as common offset sections of the earth impulse response at offsets (a) 600, (b) 1800, (c) 2400, (d) 7800, and (e) 9000 m. The observed changes in amplitude reflect a corresponding change in subsurface resistivity.
(a)
(b)
(c)
(d)
(e)Figure 16
2D apparent resistivity section for South China shale gas model.The two offset-time plots shown in Figure14 represent two types of major offset-time plot observed, showing the change in earth impulse response with depth. Evident from the peak times on the earth impulse response plot for Tx at 3000 m (Figure 14(a)) are 3 major segments representing the change in subsurface resistivity across the vertical section from the source point in a diagonal direction. The top layer is marked by a gentle slope with early peak time which increases with offset (offset: 600 m–1800 m). Considering the range of peak time values, this segment was inferred to be the highly resistive region close to the surface. Below this are the profiles from offsets 3000 m and 3600 m, which are marked by very steep slope with the highest peak time range. This segment was inferred to be the low resistivity region that represents the shale gas layer. The third segment is the profiles from offsets between 4800 m and 9000 m, whose peak times formed a steep slope. These profiles, with respect to the early profiles, have moderately high peak times and are assumed to fall within the moderately low resistivity region within the geological section. Likewise, 3 major segments were observed on Figure 14(b). The first segment (offset: 600 m–1200 m) is characterized by a gentle slope, with very low peak times, and was inferred to be a highly resistive layer. Directly below it is the steep sloped curved structure (2400 m–3000 m), which was followed by a steep sloped peak times with high values (4200 m–9000 m), representing the low resistivity medium. Finally, between each of these sections is a profile, which represents a transition from one layer to another.Figures15(a)–15(e) show a progressive change in amplitude of the earth impulse response with distance, over a range of offsets from 600 m to 9000 m. Figure 15(a) shows lateral variation in amplitude, with high value recorded across the section except between distance 5500 m and 7000 m where the value was remarkably low due to the presence of a low resistive body close to the surface. This amplitude distribution pattern was observed over a range of offsets from 600 m to 1800 m. However at offset 1800 m (Figure 15(b)), a curved relatively low amplitude region, which represents the top of a dipping low resistivity body, was observed. The vertical extension of the low amplitude region between distances 12000 m and 14000 m was apportioned to be the effect of an underlying low resistive material. As the offset increases, the amplitude distribution changed in response to the earth model, with reduction in high amplitude region (Figure 15(c)). The observed amplitude distribution at offset 2400 m, corresponds to the resistivity structure between the depth 600 m to 900 m with little variation.At an offset of 7800 m (Figure15(d)), only a fraction of the subsurface material gave a relatively high amplitude response, which was between distances 9200 m and 11,000 m and 14,000 m to 15,800 m (with a conductive body in the middle). Although this amplitude distribution pattern continued throughout the remaining sections (Figures 15(d) and 15(e)), both the width of the high amplitude regions and their distance from the origin increased proportionately with offset. This pattern represents the structurally induced gently dipping conductive material (shale gas reservoir), which is confined by the highly resistive geological body (Figure 11). In general, the amplitude distribution with offset gave a good representation of the subsurface resistivity distribution.Figure16 shows a 2D image of the apparent resistivity plot on a common midpoint-offset section. The apparent resistivity distribution on the section indicates the presence of 3 distinct geoelectric mediums. The first medium is represented by a yellowish red to dark red coloration, which reflects a resistive to highly resistive medium, respectively. This runs majorly from the surface to about a quarter of the section in vertical direction, with downward extensions between horizontal distances of about 14.5 km to 17.5 km (dipping at an angle less than 90°). The second medium is the conductive (deep blue color) shale gas formation. This runs from the surface between horizontal distances of about 5.5 km to 7 km and runs parallel with the gently dipping resistive material between horizontal distances of about 11.5 km to 15.5 km. It also runs across the section below the highly resistive layer with short discontinuity between horizontal distances 19 km and 21 km. The last medium, which has intermediate resistivity value between the mediums mentioned above, represents the base material. These 3 distinct geoelectric mediums also characterizes the inversion result obtained by Li Diquan (not yet published), using the wide-field electromagnetic method. Although the obtained apparent resistivity section was similar to the wide-field electromagnetic method’s inversion result, however, the MTEM section gave a more distinct representation of the dipping shale gas medium but at a higher angle. On the whole, the obtained apparent resistivity section was in good agreement with the South China shale gas model, thus, showing the capability of the relatively new method as an effective geophysical tool to complement seismic data in shale gas prospecting. Although few discrepancies were observed (especially the inadequate mapping of the dipping resistive material), however, the 2D modeling of MTEM earth response has been successful in the demonstration of good knowledge of the subsurface resistivity cross-section. However, some aspects of this modeling need to be improved, for example, the complex deconvolution process for earth impulse and step response recovery to accommodate noise suppression and anisotropy effect.


## 5. Conclusions

In this paper, we solved the finite difference time stepping UPML solution to the electromagnetic diffusion equation. This method when tried over a half-space gave a good solution in modeling electromagnetic diffusion within a homogeneous earth material. Furthermore, the convergence properties computed by the FDTD method as judged by comparison to analytical results shows that the method is good and our numerical results are in good agreement with the analytical results. Hence, the use of UPML equation in MTEM modeling proved productive in terms of convergence and accuracy. However, the number of cells within each PML region must be determined or chosen carefully to achieve a minimum interference with the main medium. The South China shale gas case study given above shows the high resolving power of MTEM method. The effects of the subsurface targets were noticed on the apparent resistivity section at offsets of about three times its depth. Hence, the apparent resistivity section above can be developed rapidly to locate and characterize the shale gas reserve in terms of qualitative and quantitative study and reduce the risk involved in hydrocarbon prospecting. Also, the obtained apparent resistivity either combined with inferred depth estimation from prior information or obtained via skin-depth equation can be used as a starting model for full waveform inversion, which is recommended for accurate quantitative interpretation. In conclusion, this study has confirmed the capability and efficiency of the MTEM method as an onshore unconventional shale gas reservoir mapping technique and a potential geophysical tool for energy resource prospecting in general. Thus, time lapse MTEM study will be an effective tool for shale gas exploitation via quick reconstruction of earth resistivity distribution. Also, this modeling method will be very useful in determining the best field array setup for geophysical survey over shale gas reservoir under different structural and stratigraphic setup.



---

*Source: 6863810-2016-02-29.xml*