scispace - formally typeset
SciSpace - Your AI assistant to discover and understand research papers | Product Hunt

Journal ArticleDOI

3D controlled-source electromagnetic modeling in anisotropic medium using edge-based finite element method

01 Dec 2014-Computers & Geosciences (Pergamon Press, Inc.PUB1185Elmsford, NY, USA)-Vol. 73, pp 164-176

TL;DR: The method uses the edge-based vector basis functions, which automatically enforce the divergence free conditions for electric and magnetic fields, which is effective in modeling the seafloor bathymetry using hexahedral mesh.

AbstractThis paper presents a linear edge-based finite element method for numerical modeling of 3D controlled-source electromagnetic data in an anisotropic conductive medium. We use a nonuniform rectangular mesh in order to capture the rapid change of diffusive electromagnetic field within the regions of anomalous conductivity and close to the location of the source. In order to avoid the source singularity, we solve Maxwell's equation with respect to anomalous electric field. The nonuniform rectangular mesh can be transformed to hexahedral mesh in order to simulate the bathymetry effect. The sparse system of finite element equations is solved using a quasi-minimum residual method with a Jacobian preconditioner. We have applied the developed algorithm to compute a typical MCSEM response over a 3D model of a hydrocarbon reservoir located in both isotropic and anisotropic mediums. The modeling results are in a good agreement with the solutions obtained by the integral equation method. HighlightsThis paper develops a novel formulation of the edge-based finite element method for 3D modeling of marine CSEM data in anisotropic conductive medium.The method uses the edge-based vector basis functions, which automatically enforce the divergence free conditions for electric and magnetic fields.The developed method is effective in modeling the seafloor bathymetry using hexahedral mesh.

Summary (3 min read)

1. Introduction

  • For accurate interpretation of the subsurface structure using the MCSEM method, the bathymetry effect should be accurately simulated.
  • The edge-based finite element method uses vector basis functions defined on the edges of the corresponding elements.
  • More advanced preconditioners based on the approximated inverse of the stiffness matrix can be used to speed up the solvers.

4. Comparison with semi-analytical solution for a horizontally layered geoelectrical model

  • The frequency of the harmonic electric source is 0.5 Hz.
  • The grid is refined nearby the source, target layer and the surface of observation (see Fig. 3).
  • Fig. 4 shows a comparison of the anomalous electric field between the finite element solution and the 1D semi-analytical solution.
  • One can see that the finite element results are in a good agreement with the semi-analytical solution.
  • For this model with the specified source configuration, the y component of secondary electric field, the x and z components of secondary magnetic field are equal to 0 at y¼0.

5. Model of an off-shore hydrocarbon reservoir

  • Hz, which is a typical frequency for marine CSEM.
  • The EM receivers are located on the seafloor.
  • The sparsity pattern for the finite element stiffness matrix is shown in panel (a) of Fig.
  • From this figure, one can see that the matrix is very sparse, although the problem size is huge.
  • In the second model, the authors consider an anisotropic background conductivity and isotropic anomalous conductivity for the reservoir.

5.1. Model 1: isotropic background and isotropic reservoir

  • The numerical result obtained by the edge-based finite element method was compared with the integral equation solution.
  • One can clearly see an anomaly around x¼2 km where the offset is about 5 km.
  • Due to the page limits, the authors will only show the numerical modeling result for the frequency of 1 Hz in the following sections.
  • Due to the page limits, the authors only show a comparison of the convergences of QMR, GMRES and BiCGSTAB solvers for Model 1.
  • It took about 20 min to solve this model using the finite element method and about 3 min for the integral equation method on a PC with 2.6 GHz CPU.

5.2. Model 2: anisotropic background with isotropic reservoir

  • In a marine environment, the conductivity of sediments shows a strong transverse anisotropy due to the process of sedimentation with the longitudinal conductivity being larger than the transverse conductivity (Ellis et al., 2010; Ramananjaona et al., 2011).
  • The macroanisotropy is mainly caused by thin layering when bulk resistivity is measured so that the electric current prefers to travel parallel to the bedding planes (Ellis et al., 2010).
  • This transverse anisotropy could have a strong effect on the primary field and the anomalous field could also be distorted significantly.
  • The conductivity of seawater and air stays unchanged compared to the previous model.
  • As in the previous model, the authors compare the finite element result with the integral equation solution.

5.3. Model 3: isotropic background with anisotropic reservoir

  • In practical applications of the MCSEM method, not only the conductivity of the sea-bottom sediments, but also the reservoir conductivity can be anisotropic.
  • The reservoir anisotropy is usually weak in comparison to the background conductivity anisotropy (Brown et al., 2012).
  • The authors can see that the result obtained by the integral equation method is practically the same as the finite element solution for the anisotropic reservoir model.
  • Ex components of the background and total electric field with the frequency of 1 Hz.
  • The total memory requirement for solving this problem using finite element and integral equation methods were practically the same as for Model 1.

5.4. Model 4: anisotropic background with anisotropic reservoir

  • Finally, the authors study Model 4 having transverse anisotropy of both the background and reservoir conductivities.
  • Hz, with the anomaly observed around a point x¼3 km, which corresponds to the 6 km offset.
  • The total memory requirement of this problem for finite element and integral equation methods are the same as Model 1.
  • The effect of the anisotropy of background conductivity is manifested by shifting the anomaly to larger offset, while the increase of the anisotropy coefficient of the reservoir increases the amplitude of the anomaly without shifting the anomaly significantly.
  • Thus, their numerical modeling results confirm ones again that quantitative interpretation of the MCSEM data requires taking into account the effect of anisotropy on the observed data.

6. Modeling the effect of the bathymetry on the EM data

  • One of the advantages of the edge-based finite element method is its ability to model the bathymetry effect on the EM data.
  • Fig. 23 shows the hexahedral grid for the bathymetry model without a reservoir in the X–Z section at y¼0.
  • Fig. 26 shows a comparison of the anomalous field at y¼0 along the bathymetry, with the frequency of 1 Hz, computed using both the edge-based finite element and integral equation methods.
  • One can see that the results produced by these two methods show a good agreement.
  • Some minor difference may be related to the staircase approximation used in the integral equation method.

7. Conclusions

  • The authors have developed an edge-based finite element algorithm to solve the diffusive electromagnetic problem in the 3D anisotropic medium.
  • The authors use the edge-based vector basis functions, which automatically enforce the divergence free conditions for electric and magnetic fields.
  • The sparse finite element system is solved using the quasiminimum residual method with a Jacobian preconditioner.
  • The results of numerical study confirm the accuracy and the efficiency of a new code.

Did you find this useful? Give us your feedback

...read more

Content maybe subject to copyright    Report

3D controlled-source electromagnetic modeling in anisotropic medium
using edge-based nite element method
Hongzhu Cai
a
, Bin Xiong
b,
n
, Muran Han
a
, Michael Zhdanov
a,c,d
a
Consortium for Electromagnetic Modeling and Inversion (CEMI), University of Utah, Salt Lake City, UT 84112, USA
b
College of Earth Sciences, Guilin University of Technology, Guilin, Guangxi 541004, China
c
TechnoImaging, Salt Lake City, UT 84107, USA
d
Moscow Institute of Physics and Technology, Moscow 141700, Russia
article info
Available online 5 October 2014
Keywords:
Numerical solutions
Marine electromagnetics
Electromagnetic theory
Electrical anisotropy
Finite element
abstract
This paper presents a linear edge-based nite element method for numerical modeling of 3D controlled-
source electromagnetic data in an anisotropic conductive medium. We use a nonuniform rectangular
mesh in order to capture the rapid change of diffusive electromagnetic eld within the regions of
anomalous conductivity and close to the location of the source. In order to avoid the source singularity,
we solve Maxwell's equation with respect to anomalous electric eld. The nonuniform rectangular mesh
can be transformed to hexahedral mesh in order to simulate the bathymetry effect. The sparse system of
nite element equations is solved using a quasi-minimum residual method with a Jacobian precondi-
tioner. We have applied the developed algorithm to compute a typical MCSEM response over a 3D model
of a hydrocarbon reservoir located in both isotropic and anisotropic mediums. The modeling results are
in a good agreement with the solutions obtained by the integral equation method.
& 2014 Elsevier Ltd. All rights reserved.
1. Introduction
Controlled-source electromagnetic (CSEM) method has been
widely used in geophysical exploration on land for decades (Ward
and Hohmann, 1988; Zhdanov and Keller, 1994). The marine
controlled-source electromagnetic (MCSEM) method was also
applied for the off-shore hydrocarbon (HC) exploration (Srnka
et al., 2006; Constable and Srnka, 2007; Um and Alumbaugh,
2007; Andréis and MacGregor, 2008; Zhdanov, 2010). The subsur-
face conductivity structure could be very complex due to bathy-
metry and a lateral variation of the conductivity of the sea-bottom
sediments. In this case, a full 3D modeling of diffusive electro-
magnetic data is desirable to correctly interpret the eld MCSEM
data (Silva et al., 2012). The 3D electromagnetic modeling requires
solving the diffusive Maxwell's equations in a discretized form.
The most popular numerical techniques for EM forward modeling
are integral equation (IE), nite difference (FD), and nite element
(FE) methods.
Compared to the integral equation and nite difference meth-
ods, the nite element method is more suitable for modeling of
EM response in a complex geoelectrical structure. In a 3D scenario,
the subsurface can be discretized using either regular brick,
hexahedral, or tetrahedron elements. The electric and magnetic
elds within each element can be approximated by either linear or
higher order polynomial functions. Since the support of the nite
element basis function is small, the resulting stiffness matrix is
very sparse, which makes it easy to store. In the paper, we also
compared the sparsity pattern of the stiffness matrix created by
our nite element method and that from the nite difference
method. Although the nite element stiffness matrix is less sparse
than the nite difference stiffness matrix for the same model, the
bandwidth of nite element stiffness matrix is much narrower.
The node-based nite element method was applied in the past
to model EM data by solving the coupled equations for the vector
and scalar potentials and also by solving Maxwell's equations for
electric and magnetic elds (e.g., Zhdanov, 2009). However, for
accurate computations, the divergence free condition for the
electric and magnetic elds in the source free regions needs to
be addressed by an additional penalty term to alleviate possible
spurious solutions (e.g., Jin, 2002).
The advantage of the edge-based nite element method,
introduced by Nédélec (1980), is that the divergence free condi-
tions are satised automatically by an appropriate selection of the
basis functions. The basis function of the Nédélec element is a
vector function dened along the element edges and at the center
of each edge. The tangential continuity of either electric or
magnetic eld is imposed automatically on the element's
Contents lists available at ScienceDirect
journal homepage: www.elsevier.com/locate/cageo
Computers & Geosciences
http://dx.doi.org/10.1016/j.cageo.2014.09.008
0098-3004/& 2014 Elsevier Ltd. All rights reserved.
n
Corresponding author.
E-mail addresses: caihongzhu@hotmail.com (H. Cai),
xiongbin@msn.com (B. Xiong), michael.s.zhdanov@gmail.com (M. Zhdanov).
Computers & Geosciences 73 (2014) 164176

interfaces while the normal components are still can be discontin-
uous (Jin, 2002). In this paper, we present the formulation of
Maxwell's equation for the electric eld directly using edge-based
nite element approach and the continuity of tangential electric
eld can be imposed directly. We can also formulate Maxwell's
equation for magnetic eld in the same way and the continuity of
tangential magnetic eld will be imposed directly in the formula-
tion. In spite of the fact that the edge-based nite element method
was widely used in electrical engineering for over 30 years, it
started gaining the interest from the geophysical community
recently only. Mitsuhata and Uchida (2004) implemented an
edge-based nite element modeling algorithm for solving 3D
magnetotelluric problem. Schwarzbach et al. (2011) applied linear
and higher order edge element for the modeling of marine CSEM
data using tetrahedron discretization to better simulate the sea-
oor bathymetry. Silva et al. (2012) proposed a nite element
multifrontal method which is very efcient for 3D CSEM modeling
in the frequency domain. One needs to note that all these
formulations of the 3D CSEM problem assume the subsurface
conductivity to be isotropic.
In the marine environment, the subsurface conductivity is
usually characterized by strong anisotropy due to sedimentation.
Generally the subsurface is more conductive in the horizontal
direction than in the vertical direction for a horizontally stratied
medium. The anisotropy of conductive sediment can affect the
response of the electromagnetic eld in a marine CSEM survey and
this effect has already been well studied (Ramananjaona et al.,
2011; Ellis et al., 2010; Brown et al., 2012; Newman et al., 2010).
Obviously, to accurately interpret the marine CSEM data, the
conductivity anisotropy needs to be considered in the forward
modeling. There are already series of papers published on 2.5D
and 3D modeling of marine CSEM data in the anisotropic medium
using node-based nite element and nite difference methods
(Kong et al., 2008; Weiss and Newman, 2002).
Meanwhile, another challenge arising in the interpretation of
MCSEM data is strong distortion of the data by the effect of
seaoor bathymetry (e.g., Sasaki, 2011). For accurate interpretation
of the subsurface structure using the MCSEM method, the bathy-
metry effect should be accurately simulated. The nite element
method is very well suited to solve this problem.
In this paper, we formulate the 3D marine CSEM problem using
the linear edge-element method in the anisotropic medium. In a
general case, we assume that the model has a triaxial conductivity
anisotropy. In order to compare the EM response with integral
equation solution (Zhdanov et al., 2006), we also consider trans-
verse anisotropy in our model study. To avoid the source singu-
larity, we solve Maxwell's equations with respect to anomalous
electric
eld. The background EM eld for the layered background
model is computed using Hankel transforms ( Anderson, 1989;
Guptasarma and Singh, 1997). For simplicity, we use a rectangular
element for the at seaoor model. In order to simulate the
bathymetry effect, the rectangular element is transformed into
hexahedral one by shifting the vertical coordinate. The sparse
nite element system of equations is solved using a quasi-mini-
mum residual method (QMR) with a Jacobian preconditioner.
To validate our code, we rst test it for a 1D model with an
analytical solution. For a full 3D anisotropic problem, we compare
the numerical results from our method and integral equation
solutions.
2. Formulation of the EM eld equations with respect to
anomalous eld in anisotropic medium
The low frequency electromagnetic eld, considered in geo-
physical application, satises the following Maxwell's equations
(Zhdanov, 2009):
ωμ∇× =iEH
(1)
0
σ∇× = +
¯
HJ E
(2)
s
where we adopt the harmonic time dependence
ω
e
i
t
,
ω
is the
angular frequency,
μ
0
is the free space magnetic permeability,
J
s
is
the distribution of source current, and the term
σ
¯
E
is the induced
current in the conductive earth,
σ
¯
is the conductivity tensor which
is dened as follows:
σ
σ
σ
σ
¯
=
00
00
00
(3)
x
y
z
In (3),
σσσ,,
xyz
are principle conductivities. Actually, our formula-
tion works for a general anisotropy case where the tensor has six
independent components. For simplicity, we consider that our
coordinate axes coincide with the principal axes of the conductiv-
ity tensor. In marine environment, we consider a transverse
anisotropy that
σσ σ=≠.
(4)
xyz
In a case of a total eld formulation of numerical modeling using
the nite element method, the grid needs to be rened in order to
capture the rapid change of the primary current. To overcome this
difculty, anomalous eld formulation is desirable. In the anom-
alous eld formulation of diffusive EM eld problem, the total eld
is decomposed into background and anomalous elds (Zhdanov,
2009)
=+EE E,
(5)
ba
σσσ Δ
¯
=
¯
+
¯
.
(6)
b
Based on this decomposition, one can derive the following
equation for the anomalous electric eld:
σ
ωμσ ωμ
Δ
∇×∇×
¯
=
¯
iiEE E.
(7)
aa b
From (7), we can see that the source term for this equation is the
primary electric eld, which is much smoother than the source
current. The normal electric eld can be computed analytically for
a full-space and half-space background conductivity. For a
general layered earth model, the normal eld can be computed
semi-analytically by using a digital lter to calculate Hankel
transforms.
The differential equation for anomalous electric eld can be
solved by using integral equation, nite difference or nite
element method. Once the anomalous electric eld is solved
numerically, the anomalous magnetic eld can be obtained by
using Faraday's law (Silva et al., 2012)
ωμ=∇×
iHE()
(8)
a a
1
3. Edge-based nite element analysis
The edge-based nite element method uses vector basis
functions dened on the edges of the corresponding elements.
Similar to the conventional node-based nite element method,
the modeling domain can be discretized using rectangular,
tetrahedron, hexahedron or other complex elements (Jin,
2002).
For simplicity, we will discuss the rectangular grid rst. Fig. 1 is
an illustration of the bricks grid that we used with node and edge
H. Cai et al. / Computers & Geosciences 73 (2014) 164176 165

indexing. Following the work by Jin (2002), we denote the edge
length in
xy z,,
directions as
l
ll,,
x
e
y
e
z
e
and the center of the
element as
xyz
(
,,
)
c
e
c
e
c
e
. The tangential component of the electric
eld is assigned to the center of each edge, the
xyz,,
compo-
nents of the electric eld inside the rectangular prism can be
expressed as follows:
∑∑
===
== =
ENEENEENE,,,
(9)
x
e
i
xi
e
xi
e
y
e
i
yi
e
yi
e
z
e
i
zi
e
zi
e
1
4
1
4
1
4
where the edge basis functions are dened by the following
expressions (Jin, 2002):
=++N
ll
y
l
yz
l
z
1
22
,
(10)
x
e
y
e
z
e
c
e
y
e
c
e
z
e
1
=−++N
ll
yy
l
z
l
z
1
22
,
(11)
x
e
y
e
z
e
c
e
y
e
c
e
z
e
2
=++N
ll
y
l
yz z
l
1
22
,
(12)
x
e
y
e
z
e
c
e
y
e
c
e
z
e
3
=−++N
ll
yy
l
zz
l
1
22
,
(13)
x
e
y
e
z
e
c
e
y
e
c
e
z
e
4
=++N
ll
z
l
zx
l
x
1
22
,
(14)
y
e
z
e
x
e
c
e
z
e
c
e
x
e
1
=−++N
ll
zz
l
x
l
x
1
22
,
(15)
y
e
z
e
x
e
c
e
z
e
c
e
x
e
2
=++N
ll
z
l
zx x
l
1
22
,
(16)
y
e
z
e
x
e
c
e
z
e
c
e
x
e
3
=−++N
ll
zz
l
xx
l
1
22
,
(17)
y
e
z
e
x
e
c
e
z
e
c
e
x
e
4
=++N
ll
x
l
xy
l
y
1
22
,
(18)
z
e
x
e
y
e
c
e
x
e
c
e
y
e
1
=−++N
ll
xx
l
y
l
y
1
22
,
(19)
z
e
x
e
y
e
c
e
x
e
c
e
y
e
2
=++N
ll
x
l
xz y
l
1
22
,
(20)
z
e
x
e
y
e
c
e
x
e
c
e
y
e
3
=−++N
ll
xx
l
yy
l
1
22
.
(21)
z
e
x
e
y
e
c
e
x
e
c
e
y
e
4
Eq. (9) can be written in a more compact form as
=
=
EEN,
(22)
e
i
i
e
i
e
1
12
where
=
^
=
^
=
^
++
Nx Ny NzNN N,,
(23)
i
e
xi
e
i
e
yi
e
i
e
zi
e
48
for
=i 1, 2, 3, 4
.
It is easy to verify that the vector edge basis functions are
divergence free but not curl free
∇· = × NN00, .
(24)
i
e
i
e
The vector basis functions are also continuous at the element
boundaries. As a result, the divergence free condition of the electric
eld in the source free region and the continuity conditions are
imposed dir ectly in the edge-based nite element formulation.
By substituting (9) and (22) into (7), and using Galerkin's
method, one can nd the weak form of the original differential
equation as follows:
ωμσ ωμ σ= · ∇×∇×
¯
−Δ
¯
Ω
RiidvNEEE[],
(25)
ii s s p
where
ω
is the modeling domain.
After applying the rst vector Green's theorem, we can nd the
discretized form of (25) for each element
ωμ σ ωμ σ=−
¯
−Δ
¯
=
RKEiMEiME[],
(26)
i
e
e
N
e
s
e
e
e
s
e
e
e
p
e
1
e
where
K
e
and
M
e
are the local stiffness matrices dened as follows:
=∇×·×
Ω
KdvNN()(),
(27)
ij
e
i
e
j
e
e
Ω
MdvNN ,
(28)
ij
e
i
e
j
e
e
and
Ω
e
indicates the domain for one element. The integrals in (27)
and (28) can be calculated analytically for the rectangular ele-
ments (Jin, 2002).
After assembling the local element matrices in (26) into a
global system, one can obtain a sparse linear system of equations
as follows:
=Ae b.
(29)
In order to get a unique solution for this equation, proper
boundary conditions need to be added. Following the work of Jin
(2002) and Silva et al. (2012), we consider the homogeneous
Dirichlet boundary conditions in the edge element formulation
|=
Ω
e0
(30)
which holds approximately for the anomalous electric eld
at a distance from the domain with the anomalous conductivity.
Fig. 1. An illustration of rectangular brick element. The number with circle
indicates the node index and the number without circle is the index of edges.
H. Cai et al. / Computers & Geosciences 73 (2014) 164176166

For the numerical modeling, the distance, where conditions (30)
hold, can be determined based on the skin depth of the eld. One
can refer to the work of Silva et al. (2012) and Puzyrev et al. (2013)
for more details.
We use the quasi-minimum residual method with the Jacobian
preconditioner to solve the linear system of (29). In order to
capture the rapid change of electromagnetic eld close to the
source region and target area and to minimize the computational
cost, we use a non-uniform rectangular grid.
As a matter of fact, there exist many other choices of iterative
solvers for the large linear system of equations. The commonly
used iterative solvers include GMRES and BiCGSTAB methods
besides QMR. GMRES is an Arnoldi-based method which only
requires one matrixvector multiplication for the every iteration.
However, this method requires large memory, because it needs all
the previously generated Arnoldi vectors to be saved (Saad, 2003;
Puzyrev et al., 2013). BiCGStab (Van der Vorst, 1992) and QMR
(Freund and Nachtigal, 1991) are both Lanczos-based approaches.
Comparing with GMRES, these two methods require two matrix
vector multiplications in the every iteration, but the memory
requirements for these two methods are much less comparing to
the GMRES method (Puzyrev et al., 2013). In our application, the
stiffness matrix is complex symmetric, and in this case the QMR
method requires just one matrixvector multiplication per itera-
tion (Puzyrev et al., 2013; Weiss and Newman, 2002). Therefore,
the QMR method is suitable for our application in terms of both
computation time and memory requirements.
The convergence behavior of Krylov subspace based iterative
solvers strongly depends on the conditioner number of the matrix.
The computation time for solving the linear system of equations
can be reduced by applying preconditioner to improve the condi-
tioner number of the matrix (Van der Vorst, 2003). There are also
many choices of preconditioners. Among these preconditioners,
the Jacobian preconditioner is the simplest one which does not
require extra computation (Axelsson, 1994). This type of precondi-
tioners is demonstrated to be effective for a general case and
should be used when there are no other available preconditioners
(Axelsson, 1994). More advanced preconditioners based on the
approximated inverse of the stiffness matrix can be used to speed
up the solvers. In this paper, we have adopted the Jacobian
preconditioner for simplicity and because it provided an adequate
result for demonstration of our algorithm. In future, we will
consider a more complex choice of the preconditioner.
We should note that the rectangular elements are not quite
suitable for bathymetry modeling. It is more appropriate to use the
hexahedral elements or transformed rectangular prism elements
in order to simulate the bathymetry. In order to compute the
stiffness matrix for the hexahedral element, we transform this
element into a cubic element centered at the origin. Fig. 2a shows
an arbitrary hexahedral element in the xyz-coordinate system,
while Fig. 2b shows the transformed cubic element in the
ξηζ
-
coordinate system.
The transformation can be described by the following formulas
(Jin, 2002):
ξηζ=
=
xN x(, , ) ,
(31)
i
i
e
i
e
1
8
ξηζ=
=
yN y(, , ) ,
(32)
i
i
e
i
e
1
8
ξηζ=
=
zN z(, , ) ,
(33)
i
i
e
i
e
1
8
where the scalar node-based shape function,
N
i
,isdened as
follows:
ξηζ ξξ ηη ζζ=+ + +N ( , , ) (1 )(1 )(1 ),
(34)
i
e
i
i
i
1
8
and i is a local node index of the element.
The vector basis functions for the edge-based elements can be
dened accordingly as follows:
ηη ζζ ξ
ξξ ζζ η
ξξ ηη ζ
=+ +
=+ +
=+ +
l
i
l
i
l
i
N
N
N
8
(1 )(1 ) , 1 4,
8
(1 )(1 ) , 5 8,
8
(1 )(1 ) , 9 12,
(35)
i
e
i
e
i
i
i
e
i
e
ii
i
e
i
e
i
i
where
l
i
e
is the length of the
i
th edge of the element.
It was shown above that for vector nite element analysis one
needs to evaluate two volume integrals (27) and (28). For a
hexahedral element, these two volume integrals can be trans-
formed to the integral in the
ξηζ
-coordinate system by using the
Jacobian transform
∫∫∫
ξηζ ξ η ζ∇∇·×||
−−−
KJdddNN()()(,,) ,
(36)
ij
e
i
e
j
e
1
1
1
1
1
1
∫∫∫
ξηζ ξ η ζ||
−−−
MJdddNN (, , ) ,
(37)
ij
e
i
e
j
e
1
1
1
1
1
1
where
ξηζ
J
(, , )
is the Jacobian matrix and
ξηζ||
J
(, , )
is its determi-
nant
ξηζ
ξξξ
ηηη
ζζζ
=
J
xyz
xyz
xyz
(, , ) .
(38)
The components of Jacobian matrix can be found easily by taking
derivatives of (31)(through 33) with respect to
ξ
,
η
and
ζ
coordinates
ξηζ
ξηηζζ
ηξξζζ
ξξξηη
=
++
++
++
=
J
xyz
xyz
xyz
(, , )
1
8
(1 )(1 )
,,
(1 )(1 )
,,
(1 )(1 )
,,
.
(39)
k
k
k
k
k
e
k
e
k
e
k
kk
k
e
k
e
k
e
kk
k
k
e
k
e
k
e
1
8
Fig. 2. (a) Hexahedral element in the xyz-coordinate system. (b) The same element
transformed into a cubic element in the ξηζ-coordinate system.
H. Cai et al. / Computers & Geosciences 73 (2014) 164176 167

We can apply curl operator to (35). After doing some algebra, we
can nd the expression for the curl of the vector nite basis
functions
η η ζ ζ ηζη ζ ηζζ η ξ∇∇×= + + + ×
≤≤
l
i
N
8
[],
14,
(40)
i
e
i
e
i
i
i
i
i
i
ξ ξ ζ ζ ξζξ ζ ξζζ ξ η∇∇ ∇× = + + + ×
≤≤
l
i
N
8
[],
58,
(41)
i
e
i
e
i i ii ii
ξ ξ η η ξηξ η ξηη ξ ζ∇∇ ∇× = + + + ×
≤≤
l
i
N
8
[],
912.
(42)
i
e
i
e
i
i
i
i
i
i
Once the vector basis functions and their curl are obtained, the
stiffness matrix in (36) and (37) can be evaluated numerically by
using the three-dimensional Gaussian integral (Jin, 2002)
∫∫∫
∑∑
ξηζ ξ η ζ
ξη ζ=
−−−
===
fddd
WW W f
(, , )
(, , ),
(43)
i
n
j
n
k
n
ijk i
j
k
1
1
1
1
1
1
111
x 23
where
ξ
η,
i
j
and
ζ
k
are Gaussian integral points;
n
1
,
n
2
and
n
3
are
the number of Gaussian integral points along
ξ
,
η
and
ζ
axes; and
W
i
,
W
j
and
W
k
are weighting factors. In our application, the
polynomial order of the integrand is less than 3, therefore, it is
sufcient to choose
===nnn3
123
for higher accuracy. As a
result, 27 Gaussian points are selected in each element to numeri-
cally compute the local stiffness matrix.
4. Comparison with semi-analytical solution for a horizontally
layered geoelectrical model
In this section, we validate our algorithm by considering an
isotropic horizontally layered geoelectrical model as shown in
Fig. 3. The background is a seawater-sediment model with air
earth interface at z¼ 0 and the depth of seawater is 1000 m. The
conductivities of air, seawater and sediments are
−−
1
0Sm
6
1
,
3
.3 Sm
1
and
1
Sm
1
, respectively. An innite horizontal resistive
layer with a conductivity of
0
.01 Sm
1
is embedded in the sedi-
ments from the depth of 1400 m to 1500 m. The electromagnetic
eld is excited by an horizontal electric dipole oriented in the x
direction with the moment of 100 Am and located in the seawater
with the coordinates (0, 0, 950) m which is 50 m above seaoor.
The frequency of the harmonic electric source is 0.5 Hz. The
computation domain is selected based on the skin depth criteria
as
Ω =− ≤≤xy z{ 3 km , 3 km; 0.5 km 3 km}
. We use a
nonuniform rectangular grid to discretize this domain. The grid
is rened nearby the source, target layer and the surface of
observation (see Fig. 3).
For such an isotropic model, the anomalous eld caused by the
target layer can be computed semi-analytically using the Hankel
transform (Ward and Hohmann, 1988; Anderson, 1989; Zhdanov
and Keller, 1994; Guptasarma and Singh, 1997). Fig. 4 shows a
comparison of the anomalous electric eld between the nite
element solution and the 1D semi-analytical solution. Fig. 5 shows
a comparison of the anomalous magnetic eld between the nite
element solution and the 1D semi-analytical solution. One can see
that the nite element results are in a good agreement with the
semi-analytical solution. For this model with the specied source
conguration, the y component of secondary electric eld, the x
and z components of secondary magnetic eld are equal to 0 at
y¼ 0. As such, the nite element solution of these three compo-
nents are compared with the 1D semi-analytical solution at
=−y 50 m
where their values are not 0.
Fig. 3. Rectangular mesh for a horizontally layered geoelectrical model. The red
resistive layer is embedded in sediments below seawater. The red star in this gure
indicates the excitation dipole source. (For interpretation of the references to color
in this gure caption, the reader is referred to the web version of this paper.)
Fig. 4. A comparison between nite element result and the 1D semi-analytical
solution for the secondary electric eld with a frequency of 0.5 Hz on the seaoor.
The upper panel shows a comparison for the x component of secondary electric
eld at y¼ 0; the middle panel shows a comparison for the y component of
secondary electric eld at
=−y 50
m
; the lower panel shows a comparison for the
z component of secondary electric eld at y ¼ 0.
H. Cai et al. / Computers & Geosciences 73 (2014) 164176168

Figures (29)
Citations
More filters

Journal ArticleDOI
TL;DR: An edge-based finite element method for 3D CSEM modeling which is effective in modeling complex geometry such as bathymetry and capable of dealing with anisotropic conductivity is developed.
Abstract: We solve the 3D controlled-source electromagnetic (CSEM) problem using the edge-based finite element method. The modeling domain is discretized using unstructured tetrahedral mesh. We adopt the total field formulation for the quasi-static variant of Maxwell's equation and the computation cost to calculate the primary field can be saved. We adopt a new boundary condition which approximate the total field on the boundary by the primary field corresponding to the layered earth approximation of the complicated conductivity model. The primary field on the modeling boundary is calculated using fast Hankel transform. By using this new type of boundary condition, the computation cost can be reduced significantly and the modeling accuracy can be improved. We consider that the conductivity can be anisotropic. We solve the finite element system of equations using a parallelized multifrontal solver which works efficiently for multiple source and large scale electromagnetic modeling. HighlightsThis paper develops an edge-based finite element method for 3D CSEM modeling.The algorithm is capable of dealing with anisotropic conductivity.We adopt a total field formulation and propose a new advanced boundary condition.We use a parallelized multifrontal method to solve the system of equations.The developed method is effective in modeling complex geometry such as bathymetry.

45 citations


Cites background or methods from "3D controlled-source electromagneti..."

  • ...The electric field inside each element can be represented as: ∑ EE N= .e i i e i e =1 6 (7) The vector basis functions are continuous on element boundaries and the continuity conditions are imposed directly (Jin, 2002, 2014; Silva et al., 2012; Cai et al., 2014, 2015)....

    [...]

  • ...The iterative solvers were widely used for solving 3D EM forward modeling problem for less memory requirement (Axelsson, 1994; Badea et al., 2001; Cai et al., 2014; Freund and Nachtigal, 1991; Puzyrev et al., 2013; Saad, 2003)....

    [...]

  • ...For simplicity, we consider our coordinate axes coincide with the principal axes of the conductivity tensor (Cai et al., 2014)....

    [...]

  • ...…online 22 November 2016 crossmark have observed that the edge-based FE starts to gain more interests from geophysical community (Mitsuhata and Uchida, 2004; Nam et al., 2007; Mukherjee and Everett, 2011; Schwarzbach et al., 2011; Silva et al., 2012; Kordy et al., 2016; Cai et al., 2014, 2015)....

    [...]


Journal ArticleDOI
Abstract: Unstructured tetrahedral grids with local refinement facilitate the use of total-field solution approaches to geophysical electromagnetic (EM) forward problems. These approaches, when combined with the vector finite-element (FE) method and with refinement near transmitters and receivers, can give accurate solutions and can easily handle realistic models with complex geometry and topography. We have applied this approach to 3D forward modeling for fixed- and moving-loop configurations. MUMPS, a direct solver, was used to solve the linear system of equations generated by FE analysis. A direct solver is particularly suited to the moving-loop configuration for which the right side is different for every transmitter loop, but for which the coefficient matrix is unchanged. Therefore, the coefficient matrix need only be factorized once, and then the system can be solved efficiently for all different right sides. We compared our results with several typical scenarios from the literature: a conductive bric...

36 citations


Journal ArticleDOI
Abstract: Time-domain airborne EM data is currently interpreted based on an isotropic model. Sometimes, it can be problematic when working in the region with distinct dipping stratifications. In this paper, we simulate the 3D time-domain airborne EM responses over an arbitrarily anisotropic earth with topography by edge-based finite-element method. Tetrahedral meshes are used to describe the abnormal bodies with complicated shapes. We further adopt the Backward Euler scheme to discretize the time-domain diffusion equation for electric field, obtaining an unconditionally stable linear equations system. We verify the accuracy of our 3D algorithm by comparing with 1D solutions for an anisotropic half-space. Then, we switch attentions to effects of anisotropic media on the strengths and the diffusion patterns of time-domain airborne EM responses. For numerical experiments, we adopt three typical anisotropic models: 1) an anisotropic anomalous body embedded in an isotropic half-space; 2) an isotropic anomalous body embedded in an anisotropic half-space; 3) an anisotropic half-space with topography. The modeling results show that the electric anisotropy of the subsurface media has big effects on both the strengths and the distribution patterns of time-domain airborne EM responses; this effect needs to be taken into account when interpreting ATEM data in areas with distinct anisotropy.

26 citations


Journal ArticleDOI
TL;DR: An edge-based finite-element time-domain (FETD) modeling method to simulate the electromagnetic fields in 3D dispersive medium and considers the Cole-Cole model in order to take into account the frequency-dependent conductivity dispersion.
Abstract: The induced polarization (IP) method has been widely used in geophysical exploration to identify the chargeable targets such as mineral deposits. The inversion of the IP data requires modeling the IP response of 3D dispersive conductive structures. We have developed an edge-based finite-element time-domain (FETD) modeling method to simulate the electromagnetic (EM) fields in 3D dispersive medium. We solve the vector Helmholtz equation for total electric field using the edge-based finite-element method with an unstructured tetrahedral mesh. We adopt the backward propagation Euler method, which is unconditionally stable, with semi-adaptive time stepping for the time domain discretization. We use the direct solver based on a sparse LU decomposition to solve the system of equations. We consider the Cole-Cole model in order to take into account the frequency-dependent conductivity dispersion. The Cole-Cole conductivity model in frequency domain is expanded using a truncated Pade series with adaptive selection of the center frequency of the series for early and late time. This approach can significantly increase the accuracy of FETD modeling.

23 citations


Journal ArticleDOI
Abstract: We have developed a goal-oriented adaptive unstructured finite-element method based on the scattered field for 3D frequency-domain airborne electromagnetic (AEM) modeling. To guarantee the EM field divergence free within each element and the continuity conditions at electrical material interfaces, we have used the edge-based shape functions to approximate the electrical field. The posterior error for finite-element adaptive meshing procedure is estimated from the continuity of the normal component of the current density, whereas the influence functions are estimated by solving a dual forward problem. Because the imaginary part of the scattered current is discontinuous and the real part is continuous, we use the latter to estimate the posterior error. For the multisources and multifrequencies problem in AEM, we calculate the weighted posterior error for each element by considering only those transmitter-receiver pairs that do not adhere to our convergence criteria. Finally, we add a minimum volume ...

21 citations


References
More filters

Book
01 Apr 2003
TL;DR: This chapter discusses methods related to the normal equations of linear algebra, and some of the techniques used in this chapter were derived from previous chapters of this book.
Abstract: Preface 1. Background in linear algebra 2. Discretization of partial differential equations 3. Sparse matrices 4. Basic iterative methods 5. Projection methods 6. Krylov subspace methods Part I 7. Krylov subspace methods Part II 8. Methods related to the normal equations 9. Preconditioned iterations 10. Preconditioning techniques 11. Parallel implementations 12. Parallel preconditioners 13. Multigrid methods 14. Domain decomposition methods Bibliography Index.

12,575 citations


"3D controlled-source electromagneti..." refers methods in this paper

  • ...However, this method requires large memory, because it needs all the previously generated Arnoldi vectors to be saved (Saad, 2003; Puzyrev et al., 2013)....

    [...]


Journal ArticleDOI
TL;DR: Numerical experiments indicate that the new variant of Bi-CG, named Bi- CGSTAB, is often much more efficient than CG-S, so that in some cases rounding errors can even result in severe cancellation effects in the solution.
Abstract: Recently the Conjugate Gradients-Squared (CG-S) method has been proposed as an attractive variant of the Bi-Conjugate Gradients (Bi-CG) method. However, it has been observed that CG-S may lead to a rather irregular convergence behaviour, so that in some cases rounding errors can even result in severe cancellation effects in the solution. In this paper, another variant of Bi-CG is proposed which does not seem to suffer from these negative effects. Numerical experiments indicate also that the new variant, named Bi-CGSTAB, is often much more efficient than CG-S.

4,439 citations


"3D controlled-source electromagneti..." refers methods in this paper

  • ...BiCGStab (Van der Vorst, 1992) and QMR (Freund and Nachtigal, 1991) are both Lanczos-based approaches....

    [...]


Book
01 Mar 1993
Abstract: A new edition of the leading textbook on the finite element method, incorporating major advancements and further applications in the field of electromagneticsThe finite element method (FEM) is a powerful simulation technique used to solve boundary-value problems in a variety of engineering circumstances. It has been widely used for analysis of electromagnetic fields in antennas, radar scattering, RF and microwave engineering, high-speed/high-frequency circuits, wireless communication, electromagnetic compatibility, photonics, remote sensing, biomedical engineering, and space exploration.The Finite Element Method in Electromagnetics, Third Edition explains the methods processes and techniques in careful, meticulous prose and covers not only essential finite element method theory, but also its latest developments and applicationsgiving engineers a methodical way to quickly master this very powerful numerical technique for solving practical, often complicated, electromagnetic problems.Featuring over thirty percent new material, the third edition of this essential and comprehensive text now includes:A wider range of applications, including antennas, phased arrays, electric machines, high-frequency circuits, and crystal photonicsThe finite element analysis of wave propagation, scattering, and radiation in periodic structuresThe time-domain finite element method for analysis of wideband antennas and transient electromagnetic phenomenaNovel domain decomposition techniques for parallel computation and efficient simulation of large-scale problems, such as phased-array antennas and photonic crystalsAlong with a great many examples, The Finite Element Method in Electromagnetics is an ideal book for engineering students as well as for professionals in the field.

3,699 citations


"3D controlled-source electromagneti..." refers background or methods in this paper

  • ...…E N E E N E, , , (9) x e i xi e xi e y e i yi e yi e z e i zi e zi e 1 4 1 4 1 4 where the edge basis functions are defined by the following expressions (Jin, 2002): = + − + −N l l y l y z l z 1 2 2 , (10) x e y e z e c e y e c e z e 1 ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ = − + + −N l l y y l z l z 1 2 2 ,…...

    [...]

  • ...The integrals in (27) and (28) can be calculated analytically for the rectangular elements (Jin, 2002)....

    [...]

  • ...Following the work of Jin (2002) and Silva et al. (2012), we consider the homogeneous Dirichlet boundary conditions in the edge element formulation | =Ω∂e 0 (30) which holds approximately for the anomalous electric field at a distance from the domain with the anomalous conductivity....

    [...]

  • ...Similar to the conventional node-based finite element method, the modeling domain can be discretized using rectangular, tetrahedron, hexahedron or other complex elements (Jin, 2002)....

    [...]

  • ...The tangential continuity of either electric or magnetic field is imposed automatically on the element's interfaces while the normal components are still can be discontinuous (Jin, 2002)....

    [...]


Journal ArticleDOI
Abstract: We present here some new families of non conforming finite elements in ?3. These two families of finite elements, built on tetrahedrons or on cubes are respectively conforming in the spacesH(curl) andH(div). We give some applications of these elements for the approximation of Maxwell's equations and equations of elasticity.

2,774 citations


"3D controlled-source electromagneti..." refers methods in this paper

  • ...The advantage of the edge-based finite element method, introduced by Nédélec (1980), is that the divergence free conditions are satisfied automatically by an appropriate selection of the basis functions....

    [...]


01 Jan 2000
Abstract: This Key Note presents a summary of the development of the Finite Element Method in the field of Electromagnet ic Engineering, together with a description of several contributions of the authors to the Finite Element Method and its application to the solution of electromagnetic problems. First, a self-adaptive mesh scheme is presented in the context of the quasi-static and full-wave analysis of general anisotropic multiconductor arbitrary shaped waveguiding structures. A comparison between two a posteriori error estimates is done. The first one is based on the complete residual of the differential equations defining the problem. The second one is based on a recovery or smoothing technique of the electromagnetic field. Next, an implementation of the first family of Nedelec's curl-conforming elements done by the authors is outlined. Its features are highlighted and compared with other curl-conforming elements. A presentation of an iterative procedure using a numerically exact radiation condition for the analysis of open (scattering and radiation) problems follows. Other contributions of the authors, like the use of wavelet like basis functions and an implementation of a Time Domain Finite Element Method in the context of two-dimensional scattering problems are only mentioned due to the lack of space.

2,182 citations


"3D controlled-source electromagneti..." refers background or methods in this paper

  • ...…E N E E N E, , , (9) x e i xi e xi e y e i yi e yi e z e i zi e zi e 1 4 1 4 1 4 where the edge basis functions are defined by the following expressions (Jin, 2002): = + − + −N l l y l y z l z 1 2 2 , (10) x e y e z e c e y e c e z e 1 ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ = − + + −N l l y y l z l z 1 2 2 ,…...

    [...]

  • ...The integrals in (27) and (28) can be calculated analytically for the rectangular elements (Jin, 2002)....

    [...]

  • ...The transformation can be described by the following formulas (Jin, 2002):...

    [...]

  • ...Following the work of Jin (2002) and Silva et al. (2012), we consider the homogeneous Dirichlet boundary conditions in the edge element formulation | =Ω∂e 0 (30) which holds approximately for the anomalous electric field at a distance from the domain with the anomalous conductivity....

    [...]

  • ...Similar to the conventional node-based finite element method, the modeling domain can be discretized using rectangular, tetrahedron, hexahedron or other complex elements (Jin, 2002)....

    [...]


Frequently Asked Questions (2)
Q1. What have the authors contributed in "3d controlled-source electromagnetic modeling in anisotropic medium using edge-based finite element method" ?

This paper presents a linear edge-based finite element method for numerical modeling of 3D controlledsource electromagnetic data in an anisotropic conductive medium. The authors use a nonuniform rectangular mesh in order to capture the rapid change of diffusive electromagnetic field within the regions of anomalous conductivity and close to the location of the source. 

Future work will be aimed at the implementation of the high order finite elements and at the use of the unstructured tetrahedral and hexahedron meshes to include seafloor bathymetry and complex geoelectrical structures in the modeling of the MCSEM data.