# Controlled-source electromagnetic modelling using high order finite-
difference time-domain method on a nonuniform grid

Pengliang Yang1 and Rune Mittet2
1School of Mathematics, Harbin Institute of Technology, Harbin, China, 150001
E-mail<EMAIL_ADDRESS>
2Norwegian University of Science and Technology (NTNU), Norway
E-mail<EMAIL_ADDRESS>

###### Abstract

Simulation of 3D low-frequency electromagnetic fields propagating in the Earth
is computationally expensive. We present a fictitious wave domain high-order
finite-difference time-domain (FDTD) modelling method on nonuniform grids to
compute frequency-domain 3D controlled-source electromagnetic (CSEM) data. The
method overcomes the inconsistency issue widely present in the conventional
2nd order staggered grid finite difference scheme over nonuniform grid,
achieving high accuracy with arbitrarily high order scheme. The finite-
difference coefficients adaptive to the node spacings, can be accurately
computed by inverting a Vandermonde matrix system using efficient algorithm. A
generic stability condition applicable to nonuniform grids is established,
revealing the dependence of the time step and these finite-difference
coefficients. A recursion scheme using fixed point iterations is designed to
determine the stretching factor to generate the optimal nonuniform grid. The
grid stretching in our method reduces the number of grid points required in
the discretization, making it more efficient than the standard high-order FDTD
with a densely sampled uniform grid. Instead of stretching in both vertical
and horizontal directions, better accuracy of our method is observed when the
grid is stretched along the depth without horizontal stretching. The
efficiency and accuracy of our method are demonstrated by numerical examples.

##  1 Introduction

Marine controlled source electromagnetics (CSEM) provides valuable information
about subsurface resistivities and therefore potentially about pore fluids or
rocks. It is very useful to decipher subsurface properties to assist energy
exploration, in particular when combined with seismic data. The CSEM
technology relies on low-frequency electromagnetic field propagation to probe
the subsurface. The low-frequency electromagnetic (EM) field propagation does
not lend itself to an intuitive understanding in the same manner as seismic
field propagation does due to the diffusive nature of the EM field in
conductive media. Thus, three-dimensional modeling becomes an important tool
for the interpretation of CSEM data. Imaging of marine CSEM data is today
mainly done by inversion of the observed electric and/or magnetic fields.

The kernel of CSEM inversion is the numerical simulation of 3D electromagnetic
field propagation, which is computationally expensive. Reducing the simulation
time without compromising the accuracy is important. It can shorten the
turnaround time for an imaging project while reducing the investments in
computer hardware. The implementation of nonuniform grid schemes is a well
known strategy to reduce simulation time. To retain good accuracy we propose a
high-order finite-difference approach.

There are many studies on diffusive electromagnetic modelling using different
methods. Examples are the frequency-domain finite-difference method (Newman
and Alumbaugh,, 1995; Smith, 1996a, ; Mulder,, 2006; Streich,, 2009), the
frequency-domain finite-element method (Li and Key,, 2007; da Silva et al.,,
2012; Key,, 2016; Rochlitz et al.,, 2019), and the time-domain finite-
difference method (Oristaglio and Hohmann,, 1984; Wang and Hohmann,, 1993;
Taflove and Hagness,, 2005). A key fact in all numerical modelling methods is
that the computational cost and the memory requirement are connected and
cannot be splitted. A method can be very efficient if more computer memory is
available. The efficiency and accuracy of the modelling can be dramatically
hampered when the available computer resources are restricted.

Due to the diffusive nature of low-frequency CSEM fields, most of the 3D CSEM
modelling schemes resort to the frequency-domain solution of the Maxwell
equation to avoid the high computational cost dictated by the restrictive
stability condition for the direct solution in the time domain. Time-domain
methods are attractive options because they require less amount of memory than
frequency-domain modelling within a model of the same size. Another advantage
with time-domain solutions is that multiple frequencies can be extracted from
the same simulation. Both the frequency-domain finite-difference method and
the frequency-domain finite-element method formulate Maxwell equation as a
linear equation system, which may be solved using direct (Streich,, 2009) or
iterative (Smith, 1996b, ; Mulder,, 2006; Puzyrev et al.,, 2013) solvers. A
nice feature with a direct solver is that multiple right-hand sides are fast
to calculate after the system matrix has been factorized or inverted. However,
there are significant implementation challenges with this approach when
realistic size marine CSEM surveys are simulated. The memory requirements are
large even if the equation system is sparse.

The finite-difference time-domain (FDTD) modelling based on the staggered grid
proposed by Yee, (1966) has for several decades been a main workhorse for many
EM applications. The implementation of the numerical core is straight forward
and the computational efficiency is good for wave phenomena. The computational
efficiency for diffusive phenomena is rather poor in the time domain. The
system of partial differential equations can be considered stiff in this case
(Mittet,, 2010) and a very small time step is required to retain stability.
The computational efficiency can be improved significantly due to a
correspondence principle for wave and diffusion fields (Lee et al.,, 1989; de
Hoop,, 1996; Mittet,, 2010). Maaø, (2007) proposed a mixed wave and diffusion-
domain FDTD method to perform numerically efficient CSEM modelling. This
method allowed for large time steps compared to a purely diffusion-domain
solution. Mittet, (2010) proposed a high-order FDTD scheme by utilizing the
fictitious wave to diffusion-domain transformation. The simulation is
performed in the wave domain where the propagation velocity is proportional to
the square root of resistivity.

The Yee grid (staggered grid) FDTD scheme is often the method of choice due to
a good agreement with physics. It gives divergence free magnetic fields and
electric currents (Smith, 1996a, ). This standard scheme proposed by Yee,
(1966) is based on the second-order approximation of the first derivatives
assuming an equispaced mesh. On the uniform grid, moving from second-order
FDTD to high-order FDTD is straight forward (Mittet,, 2010), and gives
improved modelling accuracy due to the reduction of spatial dispersion errors.
To improve the modelling efficiency, the use of nonuniform grid is widespread
(Newman and Alumbaugh,, 1995; Mulder,, 2006), however it results in
inconsistencies for the grid staggering. As illustrated in Figure 1, this
inconsistency leads to only first-order local truncation error, even though
the global accuracy may be up to second order (Monk and Süli,, 1994). The
problem is persistent and has remained unresolved for high-order schemes.

Figure 1: The 2nd order staggered-grid finite-difference scheme on (a) uniform
grid and (b) nonuniform grid. The grid points on staggered grid are obtained
by taking the midpoints from nonuniform grid to ensure the 2nd order accuracy.
This accuracy is not guaranteed on nonuniform grid since the midpoints is
inconsistent. For example, $x_{3/2}=(x_{1}+x_{2})/2$,
$x_{5/2}=(x_{2}+x_{3})/2$, but $\tilde{x}_{2}=(x_{3/2}+x_{5/2})/2\neq x_{2}$
due to uniform grid spacing. This inconsistency leads to first-order local
truncation error using staggered-grid finite-difference scheme on the
nonuniform grid.

We propose an efficient 3D CSEM simulation method with high accuracy using
_high-order FDTD on a staggered, nonuniform grid_ , following the fictitious
wave domain approach (Mittet,, 2010). To resolve the inconsistency issue in
conventional 2nd order staggered grid approach, our key recognition is that
the order of local truncation error can be arbitrarily high also on a
nonuniform grid if the finite-difference operator coefficients are adapted
properly to the variable grid spacing. The derivative operator coefficients
are calculated by inverting a Vandermonde matrix system.

To gain good modeling efficiency, we transforms the diffusive Maxwell equation
into the fictitious wave domain. The efficiency of the method is restricted by
stability condition: the stepsize in time is proportional to the inverse of
the propagation velocity of the field, affected by the node spacing. The
gridding of the same physical domain leads to different number of gridpoints,
affecting the size of the linear system to be solved. The use of nonuniform
grid helps to reduce the number of gridpoint, thus reducing the computational
cost. Unfortunately, the stability condition for high order FDTD over non-
uniform grid is non-trivial. An important contribution of this paper is to
establish a new stability condition valid for arbitrarily high order FDTD
scheme on nonuniform grid. The stability condition shows the strong dependence
between time step and the finite difference coefficients computed by inverting
the Vandermonde matrix. This is not known in EM geophysics community, as far
as we know.

To generate the optimal grid using a power law, we design a recursion scheme
using fixed point iterations to find the optimal stretching factor. We prove
the recursion scheme is guaranteed to converge. The optimal factor found by
the recursive scheme allows accurate matching of the computational domain
using given number of mesh points. The high accuracy and efficiency of this
high-order FDTD method on a nonuniform grid will be exemplified by a number of
numerical tests using reference solutions.

##  2 Theory

We utilize the correspondence principle for electromagnetic wave and diffusion
fields (Lee et al.,, 1989; de Hoop,, 1996; Mittet,, 2010) to calculate the
CSEM response efficiently. The key to a high-order local truncation error on a
nonuniform grid relies on the solutions to a Vandermonde system giving
derivative operator coefficients that adapts to local grid properties. The
stability condition is then established for FDTD modelling on the nonuniform
grid.

The Maxwell equations in a quasi-static regime (i.e., with negligible effect
of displacement currents) are written in the time domain as

$\displaystyle\nabla\times\mathbf{E}+\mu\partial_{t}\mathbf{H}$
$\displaystyle=$ $\displaystyle-\mathbf{M},$ (1)
$\displaystyle-\nabla\times\mathbf{H}+\mathbf{\sigma}\mathbf{E}$
$\displaystyle=$ $\displaystyle-\mathbf{J},$

or in the frequency domain as

$\displaystyle\nabla\times\mathbf{E}-\mathrm{i}\omega\mu\mathbf{H}$
$\displaystyle=$ $\displaystyle-\mathbf{M},$ (2)
$\displaystyle-\nabla\times\mathbf{H}+\mathbf{\sigma}\mathbf{E}$
$\displaystyle=$ $\displaystyle-\mathbf{J},$

where $\mathbf{E}=(E_{x},E_{y},E_{z})^{T}$ and
$\mathbf{H}=(H_{x},H_{y},H_{z})^{T}$ are electric and magnetic fields. The
magnetic permeability is $\mu$. The conductivity is a symmetric $3\times 3$
tensor: $\sigma_{ij}=\sigma_{ji},\;i,j\in\\{x,y,z\\}$. An isotropic medium
means that only the diagonal elements of the conductivity tensor are non-zeros
and the same in all directions: $\sigma_{xx}=\sigma_{yy}=\sigma_{zz}$;
$\sigma_{ij}=0,i\neq j$. The vertical transverse isotropic (VTI) medium
implemented here still has only diagonal elements, but the vertical and the
horizontal conductivities may differ, i.e.,
$\sigma_{h}:=\sigma_{xx}=\sigma_{yy}$, $\sigma_{v}=\sigma_{zz}$. We use the
following Fourier transform convention,
$\partial_{t}\leftrightarrow-\mathrm{i}\omega$.

To speed up the FDTD modelling, we transform the above system from the
diffusion to the wave domain, following Mittet, (2010). The idea is to define
a fictitious dielectric permittivity in equation 2 as
$\sigma=2\omega_{0}\varepsilon$, yielding

$\displaystyle\nabla\times\mathbf{E}-\mathrm{i}\omega\mu\mathbf{H}$
$\displaystyle=$ $\displaystyle-\mathbf{M},$ (3)
$\displaystyle-\nabla\times\mathbf{H}+2\omega_{0}\varepsilon\mathbf{E}$
$\displaystyle=$ $\displaystyle-\mathbf{J},$

which gives the following relation after multiplying the second equation with
$\sqrt{-\mathrm{i}\omega/2\omega_{0}}$

$\displaystyle\nabla\times\mathbf{E}+\underbrace{\sqrt{-\mathrm{i}2\omega\omega_{0}}}_{-\mathrm{i}\omega^{\prime}}\mu\underbrace{\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}\mathbf{H}}_{\mathbf{H}^{\prime}}$
$\displaystyle=$ $\displaystyle-\mathbf{M},$ (4)
$\displaystyle-\nabla\times\underbrace{\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}\mathbf{H}}_{\mathbf{H}^{\prime}}+\underbrace{\sqrt{-\mathrm{i}2\omega\omega_{0}}}_{-\mathrm{i}\omega^{\prime}}\varepsilon\mathbf{E}$
$\displaystyle=$
$\displaystyle-\underbrace{\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}\mathbf{J},}_{\mathbf{J}^{\prime}}$

which translates into the wave and simulation domain as the time dependent
system

$\displaystyle\nabla\times\mathbf{E}^{\prime}+\mu\partial_{t}\mathbf{H}^{\prime}$
$\displaystyle=$ $\displaystyle-\mathbf{M}^{\prime},$ (5)
$\displaystyle-\nabla\times\mathbf{H}^{\prime}+\varepsilon\partial_{t}\mathbf{E}^{\prime}$
$\displaystyle=$ $\displaystyle-\mathbf{J}^{\prime}.$

We have introduced a prime to identify the fields in the wave domain.

From the electromagnetic fields in the wave domain, the frequency-domain
fields can be computed on the fly during modelling using the fictitious wave
transformation, exemplified by the electric field here,

$\displaystyle\mathbf{E}^{\prime}(\mathbf{x},\omega^{\prime})=\int_{0}^{T_{\max}}\mathbf{E}^{\prime}(\mathbf{x},t)e^{\mathrm{i}\omega^{\prime}t}\mathrm{d}t,$
(6)

where $T_{\max}$ is the final time until the field
$E^{\prime}(\mathbf{x},\omega^{\prime})$ reaches its steady state and where

$\displaystyle\omega^{\prime}=(1+\mathrm{i})\sqrt{\omega\omega_{0}}.$ (7)

In order to have results valid for the frequency domain we need to calculate
the Green’s functions. We need the following relation,

$\displaystyle\mathbf{E}^{\prime}$ $\displaystyle=$ $\displaystyle\mathbf{E},$
(8) $\displaystyle\mathbf{H}^{\prime}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}\mathbf{H},$
$\displaystyle\mathbf{M}^{\prime}$ $\displaystyle=$ $\displaystyle\mathbf{M},$
$\displaystyle\mathbf{J}^{\prime}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}\mathbf{J}.$

The Green’s functions are then obtained by normalizing the transformed
electric and magnetic fields with the source current,

$\displaystyle
G_{kj}^{E|J}(\mathbf{x},\omega|\mathbf{x}_{s})=\frac{E_{k}(\mathbf{x},\omega|\mathbf{x}_{s})}{J_{j}(\omega)}=\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}\frac{E_{k}^{\prime}(\mathbf{x},\omega|\mathbf{x}_{s})}{J_{j}^{\prime}(\omega)},$
(9) $\displaystyle
G_{kj}^{H|J}(\mathbf{x},\omega|\mathbf{x}_{s})=\frac{H_{k}(\mathbf{x},\omega|\mathbf{x}_{s})}{J_{j}(\omega)}=\frac{H_{k}^{\prime}(\mathbf{x},\omega|\mathbf{x}_{s})}{J_{j}^{\prime}(\omega)},$

where $G_{kj}^{E|J}(\mathbf{x},\omega|\mathbf{x}_{s})$ and
$G_{kj}^{H|J}(\mathbf{x},\omega|\mathbf{x}_{s})$ stand for the electrical and
magnetic Green’s function for angular frequency $\omega$ at spatial location
$\mathbf{x}$ with the source located at $\mathbf{x}_{s}$.

Equation 4 is a pure wave-domain equation and the time integration can easily
be discretized using the leap-frog method. We let the time be $t_{n}=n\Delta
t$ with $n$ the integer time variable and $\Delta t$ the time step. We also
introduce $N=n+\frac{1}{2}$ such that,

$\displaystyle{\mathbf{H}^{\prime}}^{N}$ $\displaystyle=$
$\displaystyle{\mathbf{H}^{\prime}}^{N-1}+\Delta
t\mu^{-1}(-\nabla\times{\mathbf{E}^{\prime}}^{n}-{\mathbf{M}^{\prime}}^{n}),$
(10) $\displaystyle{\mathbf{E}^{\prime}}^{n+1}$ $\displaystyle=$
$\displaystyle{\mathbf{E}^{\prime}}^{n}+\Delta
t\epsilon^{-1}(\nabla\times{\mathbf{H}^{\prime}}^{N}-{\mathbf{J}^{\prime}}^{N}).$

The time integration of these equations is second-order accurate. It is shown
in Mittet, (2010) that the calculation of the desired fields in the “real
world” diffusive domain is independent of the frequency content of the source
term used for calculating the fictitious fields. We exploit this fact and
achieve good accuracy for the time integration by transmitting a low-frequency
signal in the fictitious wave domain. Here we are concerned with the spatial
part of the simulation scheme so we turn to this topic next.

We use a similar notation for the space variables as for the time variables
where we write $x_{i}=x_{i-1}+\Delta x_{i}$ where $\Delta x_{i}$ is the node
separation between node $x_{i-1}$ and node $x_{i}$. Likewise, we assume a
forward staggered grid such that $x_{I}=x_{I-1}+\Delta x_{I}$ where $\Delta
x_{I}$ is the node separation between node $x_{I-1}$ and node $x_{I}$. For a
uniform staggered grid we have a constant node separation such that $\Delta
x_{I}=\Delta x_{i}=\Delta x$ and $I=i+\frac{1}{2}$. The $y$ and $z$ directions
can be described in the same way with lower case and upper case integer
arguments.

Calculation of the partial derivative of the field $f(x)$ can, in the
continuous case, be formulated as an integral operator by

$\displaystyle\partial_{x}f(x)=\int_{-\infty}^{\infty}\mathrm{d}x^{\prime}f(x+x^{\prime})\\{-\partial_{x^{\prime}}\delta(x^{\prime})\\}=\int_{-\infty}^{\infty}\mathrm{d}x^{\prime}f(x+x^{\prime})\alpha(x^{\prime}).$
(11)

The discrete formulation, with $f(x_{i}):=f(i)$, is,

$\displaystyle\partial_{x}f(x_{i})\approx
D_{x}f(i)=\sum_{l=-L}^{L}f(i+l)\alpha_{l}(i),$ (12)

where $\alpha_{l}(i)$ is a band-limited approximation to the operator
$\alpha(x^{\prime})$ in equation 11. The half length of the operator is $L$.
The argument $i$ is used to explicitly show that this operator will vary with
location for a nonuniform grid. For the staggered grid, we can then define
discretized forward, $D_{x}^{+}$, and backward, $D_{x}^{-}$, derivative
operators as,

$\displaystyle
D_{x}^{+}f(i)=\partial_{x}f(I)=\sum_{l=1}^{L}f(i+l)\alpha_{l}(i)-f(i-l+1)\alpha_{-l}(i),$
(13) $\displaystyle
D_{x}^{-}f(I)=\partial_{x}f(i)=\sum_{l=1}^{L}f(I+l-1)\alpha_{l}(I)-f(I-l)\alpha_{-l}(I),$

which is the form we need for nonuniform grids and which is investigated here.
The $f^{{}^{\prime}}$ implies spatial derivative in equation 13.

The operator simplifies for uniform grids where the operator becomes
independent of spatial location such that
$\alpha_{l}(i)=\alpha_{-l}(i)=\alpha_{l}$

$\displaystyle
D_{x}^{+}f(i)=\partial_{x}f(I)=\sum_{l=1}^{L}(f(i+l)-f(i-l+1))\alpha_{l},$
(14) $\displaystyle
D_{x}^{-}f(I)=\partial_{x}f(i)=\sum_{l=1}^{L}(f(I+l-1)-f(I-l))\alpha_{l}.$

If we use $L=1$ we have that $\alpha_{1}=1/\Delta x$ and equation 14
formulates the well known second-order accurate partial derivative operations,

$\displaystyle D_{x}^{+}f(i)=\partial_{x}f(I)=(f(i+1)-f(i))/\Delta x,$ (15)
$\displaystyle D_{x}^{-}f(I)=\partial_{x}f(i)=(f(I)-f(I-1))/\Delta x.$

Let us assume a non-magnetic subsurface so that $\mu$ has the same value as in
the vacuum, $\mathbf{M}=0$, while $\epsilon$ is a diagonal tensor
$\epsilon=\mbox{diag}(\epsilon_{ii}),i=x,y,z$. The staggering is as in Mittet,
(2010),

$\displaystyle{{H}^{\prime}}_{x}^{N}(i,J,K),\quad{{H}^{\prime}}_{y}^{N}(I,j,K),\quad{{H}^{\prime}}_{z}^{N}(I,J,k),$
(16)
$\displaystyle{{E}^{\prime}}_{x}^{n}(I,j,k),\quad{{E}^{\prime}}_{y}^{n}(i,J,k),\quad{{E}^{\prime}}_{z}^{n}(i,j,K),$
$\displaystyle{{J}^{\prime}}_{x}^{n}(I,j,k),\quad{{J}^{\prime}}_{y}^{n}(i,J,k),\quad{{J}^{\prime}}_{z}^{n}(i,j,K),$
$\displaystyle{\varepsilon}_{xx}(I,j,k),\quad{\varepsilon}_{yy}(i,J,k),\quad{\varepsilon}_{zz}(i,j,K),$

and the scheme implemented is,

$\displaystyle{{H}^{\prime}}_{x}^{N}$ $\displaystyle=$
$\displaystyle{{H}^{\prime}}_{x}^{N-1}-\frac{\Delta
t}{\mu}(D_{y}^{+}{{E}^{\prime}}_{z}^{n}-D_{z}^{+}{{E}^{\prime}}_{y}^{n}),$
(17) $\displaystyle{{H}^{\prime}}_{y}^{N}$ $\displaystyle=$
$\displaystyle{{H}^{\prime}}_{y}^{N-1}-\frac{\Delta
t}{\mu}(D_{z}^{+}{{E}^{\prime}}_{x}^{n}-D_{x}^{+}{{E}^{\prime}}_{z}^{n}),$
$\displaystyle{{H}^{\prime}}_{z}^{N}$ $\displaystyle=$
$\displaystyle{{H}^{\prime}}_{z}^{N-1}-\frac{\Delta
t}{\mu}(D_{x}^{+}{{E}^{\prime}}_{y}^{n}-D_{y}^{+}{{E}^{\prime}}_{x}^{n}),$
$\displaystyle{{E}^{\prime}}_{x}^{n+1}$ $\displaystyle=$
$\displaystyle{{E}^{\prime}}_{x}^{n}+\frac{\Delta
t}{\epsilon_{xx}}(D_{y}^{-}{{H}^{\prime}}_{z}^{N}-D_{z}^{-}{{H}^{\prime}}_{y}^{N}-{{J}^{\prime}}_{x}^{N}),$
$\displaystyle{{E}^{\prime}}_{y}^{n+1}$ $\displaystyle=$
$\displaystyle{{E}^{\prime}}_{y}^{n}+\frac{\Delta
t}{\epsilon_{yy}}(D_{z}^{-}{{H}^{\prime}}_{x}^{N}-D_{x}^{-}{{H}^{\prime}}_{z}^{N}-{{J}^{\prime}}_{y}^{N}),$
$\displaystyle{{E}^{\prime}}_{z}^{n+1}$ $\displaystyle=$
$\displaystyle{{E}^{\prime}}_{z}^{n}+\frac{\Delta
t}{\epsilon_{zz}}(D_{x}^{-}{{H}^{\prime}}_{y}^{N}-D_{y}^{-}{{H}^{\prime}}_{x}^{N}-{{J}^{\prime}}_{z}^{N}).$

The conventional second-order staggered grid FDTD scheme discretizes the
spatial derivatives as follows (Newman and Alumbaugh,, 1995; Mulder,, 2006):

$\displaystyle
D_{y}^{-}H^{\prime}_{z}=\frac{H^{\prime}_{z}(I,J,k)-H^{\prime}_{z}(I,J-1,k)}{\Delta
y_{J}},\quad
D_{z}^{-}H^{\prime}_{y}=\frac{H^{\prime}_{y}(I,j,K)-H^{\prime}_{y}(I,j,K-1)}{\Delta
z_{K}},$ $\displaystyle
D_{z}^{-}H^{\prime}_{x}=\frac{H^{\prime}_{x}(i,J,K)-H^{\prime}_{x}(i,J,K-1)}{\Delta
z_{K}},\quad
D_{x}^{-}H^{\prime}_{z}=\frac{H^{\prime}_{z}(I,J,k)-H^{\prime}_{z}(I-1,J,k)}{\Delta
x_{I}},$ $\displaystyle
D_{x}^{-}H^{\prime}_{y}=\frac{H^{\prime}_{y}(I,j,K)-H^{\prime}_{y}(I-1,j,K)}{\Delta
x_{I}},\quad
D_{y}^{-}H^{\prime}_{x}=\frac{H^{\prime}_{x}(i,J,K)-H^{\prime}_{x}(i,J-1,K)}{\Delta
y_{J}},$ $\displaystyle
D_{y}^{+}E^{\prime}_{z}=\frac{E^{\prime}_{z}(i,j+1,K)-E^{\prime}_{z}(i,j,K)}{\Delta
y_{j}},\quad
D_{z}^{+}E^{\prime}_{y}=\frac{E^{\prime}_{y}(i,J,k+1)-E^{\prime}_{y}(i,J,k)}{\Delta
z_{k}},$ $\displaystyle
D_{z}^{+}E^{\prime}_{x}=\frac{E^{\prime}_{x}(I,j,k+1)-E^{\prime}_{x}(I,j,k)}{\Delta
z_{k}},\quad
D_{x}^{+}E^{\prime}_{z}=\frac{E^{\prime}_{z}(i+1,j,K)-E^{\prime}_{z}(i,j,K)}{\Delta
x_{i}},$ $\displaystyle
D_{x}^{+}E^{\prime}_{y}=\frac{E^{\prime}_{y}(i+1,J,k)-E^{\prime}_{y}(i,J,k)}{\Delta
x_{i}},\quad
D_{y}^{+}E^{\prime}_{x}=\frac{E^{\prime}_{x}(I,j+1,k)-E^{\prime}_{x}(I,j,k)}{\Delta
y_{j}},$ (18)

where $\Delta x_{i}$, $\Delta y_{j}$ and $\Delta z_{k}$ are distances between
nodes on the reference grid, while $\Delta x_{I}$, $\Delta y_{J}$ and $\Delta
z_{K}$ are the distance between the grid points $(I,j,k)$ and $(I-1,j,k)$, the
distance between the grid points $(i,J,k)$ and $(i,J-1,k)$, and the distance
between the grid points $(i,j,K)$ and $(i,j,K-1)$. This scheme is second-order
accurate on a uniform grid with $\Delta x_{i}=\Delta x_{I}$ and likewise for
the other spatial directions. The scheme is consistent with using equation 14
with $L=1$ to approximate the derivatives.

A standard discretization method for nonuniform grids is to use the same
formulation as above, but where the node distance may vary along the same
spatial direction. It is well known that there are accuracy issues with this
implementation. To achieve second-order accuracy using equation 14 is not
possible. The midpoints between nodes do not align after going from the
reference grid to the staggered grid and back again. Since $\Delta
x_{i}\neq\Delta x_{I}$, we have a situation where the cell center does not
match on a nonuniform grid, as is illustrated in Figure 1b. Consequently, the
local truncation error of the resulting scheme can only reach first order.
Reduced accuracy will also be a problem if we implement a high-order FDTD
scheme ($L>1$) on a nonuniform grid, using derivative-operator coefficients
designed for a regular grid. This is unfortunate since the nonuniform grid is
potentially attractive for efficient modelling due to significant reduction of
the number of grid point. However, good accuracy can be restored if equation
13 is used instead of equation 14. The problem that remains is to calculate
the operator coefficients for equation 13.

###  2.1 Vandermonde matrix

The major difference between FDTD implementations on a uniform grid and on a
nonuniform grid lies in the design of the spatial-derivative operator
coefficients. The position of each field component on the nonuniform staggered
grid has been illustrated in Figure 2, which is similar to the staggered FDTD
on a uniform grid.

Figure 2: Electrical and magnetic fields on staggered grid (Yee,, 1966).

In order to compute the electromagnetic field as well as its derivatives with
arbitrary grid spacing, we have to do a polynomial interpolation using a
number of knots $x_{0},x_{1},\cdots,x_{n}$. According to the Taylor expansion,
we have

$\displaystyle f(x_{i})$ $\displaystyle=$ $\displaystyle
f(x)+f^{\prime}(x)(x_{i}-x)+\frac{1}{2}f^{\prime\prime}(x)(x_{i}-x)^{2}+\cdots+\frac{1}{n!}f^{(n)}(x)(x_{i}-x)^{n}+\cdots$
(19) $\displaystyle i$ $\displaystyle=$ $\displaystyle 0,1,\cdots,n.$

By defining $a_{i}(x):=f^{(i)}(x)/i!$, we end up with a polynomial of the
Newton form,

$\displaystyle f(x_{i})$ $\displaystyle=$ $\displaystyle
a_{0}(x)+a_{1}(x)(x_{i}-x)+a_{2}(x)(x_{i}-x)^{2}+\cdots+a_{n}(x)(x_{i}-x)^{n}+\cdots$
(20) $\displaystyle i$ $\displaystyle=$ $\displaystyle 0,1,\cdots,n.$

Let us consider $n+1$ distinct nodes $x_{0},x_{1},\cdots,x_{n}$ and drop the
terms $O((x_{i}-x)^{n+1})$. This builds a matrix system

$\displaystyle\underbrace{\begin{pmatrix}f(x_{0})\\\ f(x_{1})\\\ \vdots\\\
f(x_{n})\end{pmatrix}}_{\mathbf{f}}=\underbrace{\begin{pmatrix}1&x_{0}-x&(x_{0}-x)^{2}&\cdots&(x_{0}-x)^{n}\\\
1&x_{1}-x&(x_{1}-x)^{2}&\cdots&(x_{1}-x)^{n}\\\ \cdots\\\
1&x_{n}-x&(x_{n}-x)^{2}&\cdots&(x_{n}-x)^{n}\\\
\end{pmatrix}}_{\mathbf{V}^{T}(x_{0}-x,\cdots,x_{n}-x)}\underbrace{\begin{pmatrix}a_{0}(x)\\\
a_{1}(x)\\\ \vdots\\\ a_{n}(x)\end{pmatrix}}_{\mathbf{a}},$ (21)

where $\mathbf{V}^{T}(x_{0}-x,\cdots,x_{n}-x)$ is the transpose of a
Vandermonde matrix determined by $x_{0}-x,\cdots,x_{n}-x$. The above
expression implies that the function $f(x)$ and its derivatives up to the
$n$-th order at arbitrary location $x$ can be found by inverting the
Vandermonde matrix:
$(f(x),f^{\prime}(x),\cdots,f^{(n)}(x)/n!)^{T}=\mathbf{a}=[\mathbf{V}^{T}]^{-1}\mathbf{f}$.

It is well known that the Vandermonde matrix is highly ill-conditioned and
direct matrix inversion by Gaussian elimination should be avoided due to
numerical instabilities when the matrix size becomes large. Fortunately, there
exists an efficient algorithm based on the method of Björck and Pereyra,
(1970) to invert the Vandermonde matrix. In fact, the algorithm circumvents
the curse of severe ill-conditioning of the Vandermonde matrix to arrive at
arbitrarily high accuracy for the inversion (Demmel and Koev,, 2005). Compared
with Gauss-elimination of complexity $O(n^{3})$ , the Vandermonde matrix
inversion algorithm reduces the computational complexity to $O(n^{2})$.
Because the elements of the Vandermonde matrix are fully determined by the
interpolation nodes, there is no need to explicitly construct the matrix and
store it before inversion. The detailed implementation of this fast algorithm
is available in Golub, (1996, Algorithm 4.6.1).

Let the $i$-th row, $j$-th column of the inverse matrix
$[\mathbf{V}^{T}]^{-1}$ be $w_{ij}$, i.e.,
$([\mathbf{V}^{T}]^{-1})_{ij}=w_{ij},i,j=0,\cdots,n$. It also follows that

$\displaystyle\underbrace{\begin{pmatrix}a_{0}(x)\\\ a_{1}(x)\\\ \vdots\\\
a_{n}(x)\end{pmatrix}}_{\mathbf{a}}:=\underbrace{\begin{pmatrix}w_{00}&w_{01}&\cdots
w_{0n}\\\ w_{10}&w_{11}&\cdots w_{1n}\\\ \vdots\\\ w_{n0}&w_{n1}&\cdots
w_{nn}\\\
\end{pmatrix}}_{[\mathbf{V}^{T}]^{-1}}\underbrace{\begin{pmatrix}f(x_{0})\\\
f(x_{1})\\\ \vdots\\\ f(x_{n})\end{pmatrix}}_{\mathbf{f}}.$ (22)

The $i$-th row gives the explicit expression to find the $i$-th derivative

$\displaystyle\frac{1}{i!}f^{(i)}(x)=a_{i}(x)=\sum_{j=0}^{n}w_{ij}f(x_{j})=\langle
w_{i\cdot}|\mathbf{f}\rangle,\quad i=0,\cdots,n.$ (23)

The first row of the matrix $[\mathbf{V}^{T}]^{-1}$ (i.e., $w_{0j}$,
$j=0,\cdots,n$) is $f(x)$ which in analogy with equation 11 can be written in
discrete form as

$\displaystyle
f(x)=\int_{-\infty}^{\infty}\mathrm{d}x^{\prime}f(x+x^{\prime})\\{\delta(x^{\prime})\\}=\int_{-\infty}^{\infty}\mathrm{d}x^{\prime}f(x+x^{\prime})\beta(x^{\prime}).$
(24)

The discrete formulation is,

$\displaystyle f(x_{i})=\sum_{l=-L}^{L}f(i+l)\beta_{l}(i),$ (25)

where $\beta_{l}(i)$ are the coefficients for an interpolation operator
adapted to a nonuniform grid and identical to the $w_{0j}$ coefficients. The
second row of the matrix $[\mathbf{V}^{T}]^{-1}$ (i.e., $w_{1j}$,
$j=0,\cdots,n$) is $\partial_{x}f(x)$ and the continuous and discrete
representations are given in equations 11 and 12. The $\alpha_{l}(i)$
coefficients are for a derivative operator adapted to a nonuniform grid and
identical to the $w_{1j}$ coefficients. To be explicit, consider the staggered
finite-difference approximation of the first derivatives in $x$ direction
using $2L$ non-equidistant nodes. The finite-difference coefficients
$\alpha_{l}(x_{i})$ and $\alpha_{l}(x_{I})$, $l=-L+1,\cdots,L$ are the 2nd row
of the inverse of the matrices
$\mathbf{V}^{T}(x_{i+L}-x_{I},\cdots,x_{i-L+1}-x_{I})$ and
$\mathbf{V}^{T}(x_{I+L}-x_{i},\cdots,x_{I-L+1}-x_{i})$. Using the $2L$ nodes,
we achieve accuracy up to $2L$-th order in space.

In general we find that the operator coefficients (interpolation weights) for
$f^{(i)}(x)$ are $i!w_{ij}$. Given the points $x_{0},\cdots,x_{n}$ and $x$,
the Vandermonde matrix is determined and the operator coefficients can be
calculated. For the simulation we only need the derivative-operator
coefficients. It is noteworthy to mention that the Vandermonde matrix must be
non-singular to be inverted. For the derivative-operator coefficients we have
that the Vandermonde matrix is non-singular by construction of the staggered
grids. A finite-difference scheme on a staggered grid implies that the node
$x_{I}$, whose derivative is computed, will stay in the middle between the
selected nodes $x_{i}$, which ensures that the resulting Vandermonde matrix is
invertible.

The interpolation operators are useful for recording fields at arbitrary
locations and for the distribution of source contributions (Mittet,, 2017).
Just as for the derivative-operator coefficients, the interpolation
coefficients can be pre-calculated and then reused every time step. For the
interpolation operator we may find that some of the coordinates $x_{i}$
matches the interpolation points $x$, the inversion of Vandermonde matrix is
then not necessary. This might happen when the source or receiver positions
coincide with a finite-difference node. In this case, the interpolation
weights $w_{0i}$ to evaluate $f(x)$ should be exactly 1 at $x_{i}$ and $0$
elsewhere. In order to do 2D/3D simulation on nonuniform grid,
multidimensional interpolant is simply constructed by tensor products of many
1D interpolants.

The above procedure is significant as it allows us to use arbitrarily high-
order finite-difference scheme to accurately compute the electromagnetic
fields and their derivatives, typically with arbitrary grid spacing in the
rectilinear grid. This opens the door for CSEM modelling using high-order FDTD
on a nonuniform grid in a consistent framework. The computed finite-difference
coefficients may also be used to do high-order frequency-domain modelling on a
nonuniform grid, while the resulting sparse banded matrix has to be solved
accurately if sufficient computational resources are available.

###  2.2 Stability condition

Let us write down the FDTD scheme in equation 10 without source terms as
follows:

$\displaystyle\begin{cases}\mathbf{E}^{{}^{\prime}n+1}=\mathbf{E}^{{}^{\prime}n}+\Delta
t\epsilon^{-1}\nabla\times\mathbf{H}^{{}^{\prime}N}\\\
\mathbf{H}^{{}^{\prime}N+1}=\mathbf{H}^{{}^{\prime}N}-\Delta
t\mu^{-1}\nabla\times\mathbf{E}^{{}^{\prime}n+1}\end{cases},$ (26)

leading to

$\displaystyle\begin{split}\begin{bmatrix}\mathbf{E}^{{}^{\prime}n+1}\\\
\mathbf{H}^{{}^{\prime}N+1}\end{bmatrix}=\underbrace{\begin{bmatrix}\mathbf{I}&\Delta
t\epsilon^{-1}\nabla\times\\\ -\Delta t\mu^{-1}\nabla\times&\mathbf{I}-\Delta
t^{2}\epsilon^{-1}\mu^{-1}\nabla\times\nabla\times\end{bmatrix}}_{\mathbf{A}}\begin{bmatrix}\mathbf{E}^{{}^{\prime}n}\\\
\mathbf{H}^{{}^{\prime}N}\end{bmatrix}.\end{split}$ (27)

The numerical stability requires the eigenvalues of the amplification matrix
$\mathbf{A}$ to be less than or equal to 1. Assume the eigenvalue
decomposition for the amplification matrix is
$\mathbf{A}=\mathbf{\bar{V}}\mathbf{\Lambda}\mathbf{\bar{V}}^{T}$, where
$\mathbf{\bar{V}}=(\mathbf{V}_{E},\mathbf{V}_{H})^{T}$ is an unitary matrix
such that $\mathbf{\bar{V}}^{T}\mathbf{\bar{V}}=\mathbf{I}$. Then we have
$\mathbf{A}\mathbf{V}=\mathbf{V}\mathbf{\Lambda}$, yielding

$\displaystyle\begin{bmatrix}\mathbf{I}&\Delta t\epsilon^{-1}\nabla\times\\\
-\Delta t\mu^{-1}\nabla\times&\mathbf{I}-\Delta
t^{2}\epsilon^{-1}\mu^{-1}\nabla\times\nabla\times\end{bmatrix}\begin{bmatrix}\mathbf{V}_{E}\\\
\mathbf{V}_{H}\end{bmatrix}=\begin{bmatrix}\mathbf{V}_{E}\\\
\mathbf{V}_{H}\end{bmatrix}\Lambda.$ (28)

That is,

$\displaystyle\begin{cases}\Delta
t\epsilon^{-1}\mu^{-1}\nabla\times\mathbf{V}_{H}=\mathbf{V}_{E}(\mathbf{\Lambda}-\mathbf{I})\\\
-\Delta t\mu^{-1}\nabla\times\mathbf{V}_{E}-\Delta
t^{2}\epsilon^{-1}\mu^{-1}\nabla\times\nabla\times\mathbf{V}_{H}=\mathbf{V}_{H}(\mathbf{\Lambda}-\mathbf{I}).\end{cases}$
(29)

Multiplying the second sub equation $\mathbf{\Lambda}-\mathbf{I}$ from the
right and inserting the first sub equation gives

$-\Delta
t^{2}\mu^{-1}\epsilon^{-1}\nabla\times\nabla\times\mathbf{V}_{H}\mathbf{\Lambda}=\mathbf{V}_{H}(\mathbf{\Lambda}-\mathbf{I})^{2}.$
(30)

Denote $\mathbf{V}_{H,j}$ the $j$th column of $\mathbf{V}_{H}$ and
$\lambda_{j}$ the $j$th eigenvalue in $\mathbf{\Lambda}$. The above equation
reads

$\mu^{-1}\epsilon^{-1}\nabla\times\nabla\times\mathbf{V}_{H,j}=-\frac{(\lambda_{j}-1)^{2}}{\Delta
t^{2}\lambda_{j}}\mathbf{V}_{H,j},$ (31)

which shows that $-\frac{(\lambda_{i}-1)^{2}}{\Delta t^{2}\lambda_{i}}$ is the
eigenvalue of the matrix $(\mu\epsilon)^{-1}\nabla\times\nabla\times$
associated with the eigenvector $V_{H,j}$. This leads to

$(\lambda_{j}^{2}+(-2+\Delta
t^{2}c^{2}\nabla\times\nabla\times)\lambda_{j}+1)\mathbf{V}_{H,j}=(\lambda_{j}^{2}-(2+\Delta
t^{2}c^{2}\Delta)\lambda_{j}+1)\mathbf{V}_{H,j}=0,$ (32)

where we denote $c:=1/\sqrt{\mu\epsilon}$ and have applied
$\nabla\times\nabla\times\mathbf{F}=\nabla\nabla\cdot\mathbf{F}-\nabla\cdot\nabla\mathbf{F}=-\Delta\mathbf{F}$
due to Gauss theorem $\nabla\cdot\mathbf{F}=0$,
$\mathbf{F}=\mathbf{E},\mathbf{H}$ in the homogeneous, source free medium. The
roots of the above equation are

$\lambda_{j;1,2}=1+\frac{\Delta
t^{2}c^{2}\Delta}{2}\pm\frac{\mathrm{i}}{2}\sqrt{-\Delta
t^{2}c^{2}\Delta(4+\Delta t^{2}c^{2}\Delta)},$ (33)

which requires the following condition to be satisfied

$0\leq-\Delta t^{2}c^{2}\Delta\leq 4,$ (34)

in order to ensure $|\lambda_{j;1,2}|\leq 1$. Finally, we arrive at the same
stability condition as equation 41 of Mittet, (2010),

$\Delta
tc_{\max}\sqrt{(D_{x}^{\max})^{2}+(D_{y}^{\max})^{2}+(D_{z}^{\max})^{2}}\leq
2,$ (35)

where $D_{x}^{\max}$, $D_{y}^{\max}$ and $D_{z}^{\max}$ are the maximum value
of the the discretized first derivatives along $x$, $y$ and $z$ directions.
Let us emphasize this condition applies to both uniform and nonuniform grid.
The difference lies in the spatial derivative operator.

To proceed with the stability analysis, we represent the fields on the grid
via time harmonic plane waves

$\mathbf{E}^{\prime},\mathbf{H}^{\prime}\propto e^{-\mathrm{i}(\omega
t-k_{x}x-k_{y}y-k_{z}z)},$ (36)

where the amplitude has been omitted. Equation 13 becomes

$\begin{cases}D_{x}^{+}u(x_{I})=(\alpha_{L}(x_{I})e^{\mathrm{i}k_{x}(x_{i+L}-x_{I})}+\alpha_{L-1}(x_{I})e^{\mathrm{i}k_{x}(x_{i+L-1}-x_{I})}+\cdots+\alpha_{-L+1}(x_{I})e^{\mathrm{i}k_{x}(x_{i-L+1}-x_{I})})u(x_{I})\\\
D_{x}^{-}u(x_{i})=(\alpha_{L}(x_{i})e^{\mathrm{i}k_{x}(x_{I+L-1}-x_{i})}+\alpha_{L-1}(x_{i})e^{\mathrm{i}k_{x}(x_{I+L-1}-x_{i})}+\cdots+\alpha_{-L+1}(x_{i})e^{\mathrm{i}k_{x}(x_{I-L}-x_{i})})u(x_{i})\end{cases}.$
(37)

Hence, we end up with the maximum possible values for discrete first
derivative operators

$D_{x}^{\max}=\max\left(\sum_{l=-L+1}^{L}|\alpha_{l}(x_{I})|,\sum_{l=-L+1}^{L}|\alpha_{i}(x_{l})|\right)$
(38)

and similar estimations for $D_{y}^{\max}$ and $D_{z}^{\max}$. In case of a
uniform grid, the above expressions becomes much simpler

$\begin{cases}D_{x}^{+}u(x_{I})=(\alpha_{L}e^{\mathrm{i}k_{x}(L-1/2)\Delta
x}+\alpha_{L-1}e^{\mathrm{i}k_{x}(L-3/2)\Delta
x}+\cdots+\alpha_{-L+1}e^{\mathrm{i}k_{x}(-L+1/2)\Delta x}u(x_{I})\\\
D_{x}^{-}u(x_{i})=(\alpha_{L}e^{\mathrm{i}k_{x}(L-1/2)\Delta
x}+\alpha_{L-1}e^{\mathrm{i}k_{x}(L-3/2)\Delta
x}+\cdots+\alpha_{-L+1}e^{\mathrm{i}k_{x}(-L+1/2)\Delta x}u(x_{i})\\\
\end{cases},$

where $\Delta x$ stands for the uniform grid spacing in $x$ direction, while
the coefficients $\alpha_{l}$ (which can be computed by inverting a
Vandermonde matrix system according to Appendix AComputing uniform staggered-
grid finite difference coefficients via Vandermonde matrix inversion) are
independent of the location $x_{i}$.

###  2.3 Grid stretching

Our finite-difference modelling is carried out on a rectilinear mesh, which
can be generated from the tensor (outer) product of 1D non-equispaced meshes.
We use the geometrical progression to generate the 1D nonuniform grid,
following the work of Mulder, (2006, Appendix C). This is also often referred
to as power law grid stretching since the cell sizes stretch exponentially to
guarantee a smooth extension of the grid.

Assume we have the total grid length $L_{x}$ divided into $n$ intervals ($n+1$
nodes) with a common ratio $r>1$. Denote the smallest interval $\Delta
x=x_{1}-x_{0}$. Thus, the relation between $L_{x}$ and $\Delta x$ is

$L_{x}=(x_{1}-x_{0})+(x_{2}-x_{1})+\cdots+(x_{n}-x_{n-1})=\Delta
x(1+r+\cdots+\cdot r^{n-1})=\Delta x\frac{r^{n}-1}{r-1}.$ (39)

Given the total distance $L_{x}$, the smallest grid spacing $\Delta x$ and the
stretching factor $r$, we can compute an approximate value for the number of
nodes $n=\left\lceil\frac{\ln(1+\frac{L_{x}}{\Delta
x}(r-1))}{\ln(r)}\right\rceil$ following Mulder, (2006), where
$\left\lceil\cdot\right\rceil$ takes the ceiling integer value. This strategy
yields approximate solution as the value of $L_{x}$ is not exactly preserved.

Due to the stability requirement and the resulting computational cost in the
modelling, we are restricted to the smallest interval $\Delta x$ and a given
number of intervals $n$ to discretize over a certain distance $L_{x}$. The
question boils down to finding an optimal growth factor $r$. This problem is
more complicated since equation 39 does not yield an explicit expression for
the stretching factor $r$.

The relation in equation 39 is equivalent to

$r=\underbrace{\left(\frac{L_{x}}{\Delta
x}(r-1)+1\right)^{\frac{1}{n}}}_{g(r)},$ (40)

which inspires us to carry out a number of fixed point iterations until
convergence:

$r^{k+1}=g(r^{k}),\quad k=0,1,\cdots.$ (41)

Assume $r^{*}$ is the analytic solution such that $r^{*}=g(r^{*})$. Thanks to
Lagrange mean value theorem, the error estimation at ($k+1$)-th iteration is
linked with the error at $k$-th iteration via

$|e^{k+1}|=|r^{k+1}-r^{*}|=|g(r^{k})-g(r^{*})|=|g^{\prime}(\xi)(r^{k}-r^{*})|=|g^{\prime}(\xi)||e^{k}|,\quad\xi\mbox{
between }r^{k}\mbox{ and }r^{*}.$ (42)

It becomes evident that $e^{k}\rightarrow 0\;(k\rightarrow\infty)$ provided
that $|g^{\prime}(r)|<1$. Starting from any initial guess $r^{0}>1$, the fixed
point iteration scheme in equation 41 is guaranteed to converge since

$|g^{\prime}(r)|=\frac{L_{x}}{n\Delta x}\left(\frac{L_{x}}{\Delta
x}(r-1)+1\right)^{1/n-1}=\frac{1}{n}(r^{-(n-1)}+r^{-(n-2)}+\cdots+1)<1,$

thanks to the relations in equations 39 and 40.

We note that within areas of constant $r$, the derivative-operator
coefficients in equation 13 can be calculated from the $2L$ coefficients
$a_{-l}$ and $a_{l}$, such that

$\displaystyle\alpha_{-l}(i)$ $\displaystyle=$
$\displaystyle\frac{a_{-l}}{\Delta x\>r^{i}},$ (43)
$\displaystyle\alpha_{l}(i)$ $\displaystyle=$
$\displaystyle\frac{a_{l}}{\Delta x\>r^{i}}.$

Two solutions of the Vandermonde system are required in this case, one for
operators valid on the reference grid and one for operators valid on the
staggered grid.

###  2.4 Implementation

A number of techniques have been applied to achieve efficient and accurate 3D
CSEM simulation. According to equation 6, the frequency-domain CSEM response
is obtained by a time to frequency transform. The transform is using the
complex frequency given in equation 7. This gives exponential damping of late
arrivals as is discussed in more detail in Mittet, (2015). The lowest
frequency experiences the least damping in this transform and by that requires
the longest simulation time. This allows us to bound the number of times steps
to terminate the simulation when the lowest frequency component has converged.
The convergence means that the frequency-domain field obtained by time
integration has reached its steady state and later arrivals are damped to the
degree that they do not contribute to the time integral.

The source and receiver locations may be arbitrarily distributed over the
whole computational domain, not necessarily located at the nodes of the
finite-difference grid. In case they do, interpolation is not required and we
directly take the field from the grid; otherwise, we need interpolation
operators extracted from the first row of the relevant $[\mathbf{V}^{T}]^{-1}$
matrices.

The CFL condition to achieve stable FDTD modelling dictates the timestep for a
given spatial sampling. The air-wave travels at extremely high speed which
does not allow us to use the local finite-difference stencil to simulate it.
The air-water boundary condition is implemented in the Fourier-wavenumber
domain following the method proposed by Oristaglio and Hohmann, (1984) and
Wang and Hohmann, (1993) and using the extension to high-order schemes given
in Mittet, (2010).

To mimic wave propagation in unbounded domain, we use a truncated computation
mesh surrounded by an artificial absorbing boundary using convolutional
perfectly matched layer (PML) technique (Roden and Gedney,, 2000), except the
top air-water interface. To ease the implementation, we extend the domain with
equal spacing based on the last finite-difference cell size in the interior
domain. Our CPML implementation is very standard using all the parameter
setups given in Komatitsch and Martin, (2007).

##  3 Numerical examples

We now present three examples to demonstrate the merits of high-order finite
differences using nonuniform staggered grids. The first two examples are using
1D resistivity models under shallow water and deep water scenarios. In the 1D
case, the semi-analytic solution can be computed in the frequency-wavenumber
domain as a reference to benchmark our results. The third example takes into
account varying seafloor topography, in which a reference solution can be
computed using `emg3d` software (Werthmüller et al.,, 2019).

In all of our modelling, we use x-directed electrical dipole source. The
$E_{x}$ component at three commonly used frequencies - 0.25 Hz, 0.75 Hz and
1.25 Hz, are modelled. In total 12 layers of CPML are sufficient to achieve
nearly perfect absorbing effect. The computed electromagnetic fields are
normalized by the source current. The convergence check has been conducted
every 100 time steps to avoid modelling after the frequency-domain EM fields
stop evolving. We examine the amplitude error by inspecting the ratio of the
modelled field to the reference solution, $|E_{x}^{FD}|/|E_{x}^{ref}|$, which
should be close to unity if the modelling is precise. The phase difference is
computed by $\angle E_{x}^{FD}-\angle E_{x}^{ref}$ in degrees.

###  3.1 Moving to high-order schemes

The model shown in Figure 3 includes 5 layers: the top is the air, then a 325
m layer of sea water with a resistivity of 0.3 $\Omega$m, followed by three
layers of formation with increasing resistivity in the depth direction. The
whole model extends down to 5 km depth below the sea surface. A horizontal
electrical dipole source is deployed at 275 m water depth, its lateral
position being in the middle of the model. The EM fields are recorded using
201 receiver positions at the seafloor. The offsets range from -10 km to +10
km (receiver separation equals 100 m).

Figure 3: The resistivity model with air ($\rho=10^{12}$ $\Omega$m), shallow
column (325 m) of 0.3 $\Omega$m and 3 sediment layers ($\rho=1$ $\Omega$m in
[325, 1025] m; $\rho=2$ $\Omega$m in [1025, 1525] m; $\rho=4$ $\Omega$m
downwards).

We first compare the 3D CSEM response simulated by FDTD with the reference
solution calculated using `empymod` program (Werthmüller,, 2017). To limit the
factors affecting the modelling accuracy, this experiment has been done using
isotropic modeling and uniform grid spacing ($\Delta x=\Delta y=150$ m,
$\Delta z=50$ m). Panels (a) and (b) in Figure 4 shows the amplitude and the
phase of the modelled EM fields by FDTD of 2nd order, which is highly
consistent to the reference solution. At the offset beyond 1 km, the agreement
between finite-difference and the analytical solution is good. The kink in the
transition of the near-to-far offset in the amplitude response for all
frequencies is a manifestation of the strong air-wave effect. In the very near
field, where the receivers are located just below the transmitter, the phase
exhibits a $180^{o}$ jump due to the change in the direction of the electric
field immediately below the electric dipole source. Due to the extension of
finite-difference stencil, this phase rollover becomes smeared compared to the
reference solution. For the same reason, finite difference method introduces
significant errors in the near field in amplitude (which becomes more evident
for high-order schemes). Fortunately it is not an issue for practical 3D CSEM
applications since the CSEM inversion to deduce the subsurface resistivity is
mainly driven by far offset refractions. We thus focus on the analysis of the
amplitude and phase error beyond 1 km offset in the following.

Figure 4: Comparison between 2nd order FDTD (solid line) and reference
solution (dash line) for 3D CSEM simulation in the shallow water scenario. The
horizontal coordinates are offsets, while the vertical coordinates are (a)
Amplitude; (b) Phase.

The modelling results by FDTD of higher orders are not displayed in Figure 4,
since they are visually very similar to the reference solution. Instead, the
amplitude and phase error are computed using 1D reference solution to examine
the accuracy of our methods. In Figure 5a, c and e, we clearly see that the
2nd order FDTD gives the largest amplitude error for all frequencies; moving
from 2nd order to 4th order significantly reduces the amplitude error, while
moving to 6th order behaves even better, although the accuracy improvement
becomes less. The phase errors exhibit a similar behavior.

Figure 5: The amplitude and phase error of FDTD compared with 1D reference
solution for 0.25 Hz (a,b), 0.75 Hz (c,d) and 1.25 Hz (e,f).

In principle, the computational cost of the 4th order and the 6th order finite
differences will be double and triple of that for the 2nd order scheme. Table
1 lists the CPU time for these modelling exercises, running on a laptop
possessing Intel(R) Core(TM) i7-4710HQ CPU @ 2.50GHz. It shows the increase of
the computing cost by increasing the order of FDTD scheme. The computing time
of the 4th order scheme is less than two times the computing time of the 2nd
order scheme. Further increase of the FDTD order results in significant
increase of CPU time. Since moving to the 6th order scheme demands more
computation while the accuracy improvement is marginal, we stay with the 4th
order scheme from now on, as it gives the best compromise between increased
accuracy and computational load.

Table 1: Comparison of computing time using FDTD of orders 2, 4 and 6. FDTD | Order-2 | Order-4 | Order-6
---|---|---|---
Time (s) | 422.3 | 715.5 | 1090.9

###  3.2 The impact of nonuniform grid

The above example demonstrates the importance of higher order scheme to
achieve accurate CSEM modeling in the presence of strong air-wave. We now
examine the impact of grid non uniformity in achieving computational
efficiency for high-order staggered grid FDTD. To get rid of the impact of
air-wave, we consider a resistivity model in deep water scenario. As shown in
Figure 6a, the model has 1020 m water column of 0.3 $\Omega$m, followed by
formation of $1\Omega\mbox{m}$ down to 1900 m, and 120 m thickness of resistor
of 50 $\Omega$m. The background resistivity below the resistor is 2.5
$\Omega$m. To mimic a vertical transverse isotropic (VTI) Earth, all layers
below the seabed are assigned with an anisotropy ratio (defined as
$\lambda=\rho_{v}/\rho_{h}$ in this paper) $\lambda=1.5$. The source is placed
in the middle of the model, 40 m above the receivers sitting on the seabed.
The resistivity in the water column and the formation above 1200 m was
discretized with constant grid spacing: $\Delta x=\Delta y=150$ m, $\Delta
z=40$ m. From 1200 m down to the bottom of the resistivity model, the grid has
been stretched with different growing factors.

Figure 6: (a) The resistivity model in deep water: the 1020 m water column of
0.3 $\Omega$m followed by formation of $1\Omega\mbox{m}$ down to 1900 m, and
120 m thickness of resistor of 50 $\Omega$m, while the background resistivity
below the resistor is 2.5 $\Omega$m; (b) VTI anisotropy below seabed is 1.5.

Figure 7a displays the stretching factor $r$ of the nonuniform grid for
different number of grid points $n_{z}$. Note that the cell size grows very
fast with the factor $r$. For example, with $r=1.05$ after 40 cells we obtain
$r^{40}>7$ times the size of the first cell. Let us now increase $r$ (hence
decrease $n_{z}$ correspondingly) and analyze how it changes the total
modelling time for the resistivity model of the same physical length in $z$
direction. Figure 7b displays a significant decrease of the modelling time
associated with the reduction of the number of grid points $n_{z}$. We indeed
see that modelling using $n_{z}=66$ takes almost half of the simulation time
compared to modelling using $n_{z}=126$.

Figure 7: (a) The exponentially growing factor with the decreasing number of
grid points $n_{z}$; (b) With the decreasing number of grid points $n_{z}$,
the modelling time decreases dramatically.

Figure 8 overlays the FDTD modelled EM fields with the reference solution. We
see a very good agreement between the two in terms of both the amplitude and
phase. To gain an idea of the modelling accuracy reduction when using less
computing effort with nonuniform grid, we plotted the amplitude and phase
errors in Figure 9 corresponding to three different frequencies. These figures
clearly show that both the amplitude error and the phase error increase with
the increase of the frequencies. It is also interesting to note that
stretching the grid does not necessarily increase the amplitude error, but
does increase the phase discrepancy. This highlights the importance of
combination in examing both amplitude and phase.

Figure 8: Comparison of (a) amplitude and (b) phase between 4th order FDTD
(solid line) and reference solution (dash line) for 3D CSEM simulation in the
deep water scenario. Figure 9: The amplitude and phase error of FDTD compared
with 1D reference solution for 0.25 Hz (a,b), 0.75 Hz (c,d) and 1.25 Hz (e,f)
in deep water.

###  3.3 3D modelling with seafloor bathymetry

The above 1D example highlights the importance of nonuniform grid in
combination with high-order schemes to achieve efficient numerical modelling
with sufficient accuracy. Here we consider a more realistic 3D resistivity
model with bathymetry variations in horizontal directions, as shown in Figure
10a.

The model has seawater of 0.3 $\Omega$m, followed by 2 formation layers of
resistivity - 1 $\Omega$m and 2 $\Omega$m. A resistor of 50 $\Omega$m was
buried in the last formation layer to mimic a hydrocarbon bearing formation
located between 1800 m and 2000 m in depth, with the offset expanding from
-3000 m to 3000 m in both x and y directions.

Figure 10: The 3D resistivity model with seafloor bathymetry, followed by 2
formation layers of resistivity values 1 $\Omega$m and 2 $\Omega$m. The last
sediment layer includes a strong resistor of 100 $\Omega$m mimicking a
hydrocarbon bearing formation located within the depth range [1800, 2000] m
and the offset range [-3000, 3000] m in both x and y directions. The source
was placed at the center of the model while 101 receivers with 200 m spacing
are deployed along x direction.

In order to validate the modelling accuracy of the proposed method, we
simulated a reference solution by finite integration method in frequency
domain using `emg3d` (Werthmüller et al.,, 2019) The finite-integration
modelling extends the model tens of kilometers in each direction, to mimic
that the EM fields propagate to very far distance while avoiding possible
edge/boundary effects. In our finite-difference modelling, the PML boundary
condition attenuates the artificial reflections in the computational domain
within ten grid nodes to achieve the same behavior. A dipole source was placed
at 650 m depth in the middle of the model, while 101 receivers are sitting on
the curved seabed.

Figure 11 displays the comparison of the 3D CSEM modelling between our method
and the result from `emg3d`. From the Figures 11a and 11b, both the amplitude
and the phase from our method match very well with the reference solution. The
maximum amplitude discrepancy for all frequencies in Figure 11c is bounded
within 5% at most of the relevant offset. The maximum mismatch in phase is
less than 1 degree for 0.25 Hz and 2 degrees for 0.75 Hz and 1.25 Hz for the
EM field above the noise level (1e-15 V/m2), as can be seen in Figure 11d.
These demonstrate the good accuracy achieved by our method.

Figure 11: Comparison of the modelling results between the proposed method and
the reference solution.

##  4 Discussion

A natural idea to achieve higher modelling accuracy is to use denser sampling.
However, increasing the number of grid points in each dimension will lead to
exponentially growing computational overhead, i.e., double the sampling in x,
y and z coordinates results in eight times more nodes in the simulation. Over
the same domain the node separations are reduced with a factor of two. By
that, the stability criterion dictates a reduction in the time step by a
factor of two. The net result is an increase of computational cost with a
factor of sixteen. This is significantly more costly than considering high-
order FDTD. For a fixed error requirement, high-order FDTD has been
demonstrated to be much more efficient than a lower-order scheme with dense
sampling (Yefet and Petropoulos,, 2001). Assuming a linear scaling dependency
between computing time and the number of grid points, doubling the length of
the finite-difference operator in x, y and z directions will simply double the
computational cost.

Stretching along $z$-direction seems always to be beneficial in terms of
computational efficiency. It reduces the computational cost in two ways:
first, it leads to increased grid spacing, hence less grid points and larger
node spacing to discretize the resistivity model for simulation on the same
physical size. Meanwhile, larger grid spacing permits to use a larger temporal
step in terms of stability condition. This decreases the number of time steps
needed to reach the steady state of the frequency-domain EM fields.

It is natural to stretch the nonuniform grid in all three spatial directions,
with the motivation to decrease the computational cost further. In case there
is no a priori knowledge about the subsurface, a possible practice is to start
stretching from a given offset from the source location in the horizontal
directions. This approach is not followed up here due to the fact that grid
stretching complicates the calculation of the air-wave for time-domain codes.
Grid stretching is applicable also for frequency-domain finite-difference
codes. The airwave implementation is very different for frequency-domain
codes, where the air layer is part of the simulation domain. Frequency-domain
codes lend themselves easily to both vertical and horizontal grid stretching.

More research must be invested with respect to increasing the simulation
efficiency of time-domain finite-difference schemes focusing on airwave
implementation combined with horizontally nonuniform grids. The root of the
problem is that the most common airwave implementations require fast Fourier
transform (FFT) on regular grids. The fields are transformed to the wavenumber
domain, propagated into the air layer in this domain and then transformed back
to the space domain. A straight forward approach is to use interpolation
between the uniform grid and the nonuniform one, as illustrated in Figure 12.
The grid stretching in $x$\- and $y$\- directions certainly complicates the
implementation while introducing additional computational cost, which is
opposite to what we want to achieve. We have tested horizontal stretching for
the previous 1D model using the stretching factor 1.05 in both $x$\- and $y$\-
directions. The resulting amplitude and phase error in Figure 13 shows that
the numerical accuracy is highly degraded (maximum amplitude error for 0.25 Hz
is around 4%) compared with the result for a horizontally uniform grid
(maximum amplitude error is less than 1.5% for all calculated frequencies).
The running time is in fact longer than for the horizontally uniform grid
which has a higher number of nodes. Due to the nonuniform grid staggering,
several nodes on the uniform grid may reside in the same interval between two
neighboring nodes on the nonuniform grid. The error panel displays an
unsymmetrical pattern in Figure 13. Interpolating the fields from the uniform
grid back to nonuniform grid can thus produce less accurate solutions. No
stretching is therefore recommended along $x$\- and $y$\- for efficiency and
accuracy considerations until a more accurate solution to this problem is
developed.

Figure 12: Horizontal grid stretching requires interpolating between the
coarse nonuniform grid and a dense uniform grid, due to the equidistance
requirement of FFT in airwave manipulation. Figure 13: The amplitude and phase
error after horizontal stretching of the 1D model with a stretching factor
1.05 along both $x$\- and $y$\- directions.

It is noteworthy that fictitious wave domain is simply a mathematical tool to
compute correct frequency domain EM response efficiently by time stepping. The
fictitious time is different from real time. To compute the true EM time
series correctly, one needs to simulate a large number frequencies of the CSEM
fields and then perform inverse Fourier transform, as has been done in Mittet,
(2010) and Rochlitz et al., (2021).

##  5 Conclusion

We have presented a 3D CSEM modelling method using high-order FDTD on a
nonuniform grid. The key problem addressed in this work is the low accuracy
and inconsistency issue in standard 2nd order staggered FDTD scheme on
nonuniform grid. The strategy we propose is to adapt the interpolation weights
depending on the nodal distance. These finite difference coefficients can be
computed by inverting a Vandermonde matrix in an accurate and efficient
manner. This makes our approach different from the commonly used EM modelling
approaches. A new yet more generic stability condition has been established in
order to achieve stable FDTD modelling. In designing the nonuniform grid based
on geometrical progression, we develop a fixed point iteration to compute the
optimal growing factor which allows good match of the modelling domain in case
of grid stretching. The numerical examples demonstrate that there is a
significant improvement in accuracy by using a high-order FDTD scheme, while
combining it with a nonuniform grid reduces the computational cost without a
significant sacrifice of accuracy. We conclude that high-order finite
differences on nonuniform grid is a viable tool for full scale 3D CSEM
modelling applications. Since the key idea is to use high order finite-
difference coefficients adaptive to the node spacing, the method is expected
to be applicable also for finite-difference frequency-domain schemes.

## Acknowledgments

Pengliang Yang was supported by Chinese Fundamental Research Funds for the
Central Universities (AUGA5710010121) and National Natural Science Fundation
of China (42274156). Pengliang Yang thanks Dieter Werthermullter for the
assistance to produce the reference solution using `empymod` and `emg3d` to
validate the accuracy of the proposed method. The source code of this work can
be found in the github repository: https://github.com/yangpl/libEMM.

## Computing uniform staggered-grid finite difference coefficients via
Vandermonde matrix inversion

The method to invert Vandermonde matrices gives us a generic approach to
compute finite-difference coefficients with arbitrary grid spacing. A special
case is the regular grid spacing. In what follows, we show how the standard
staggered grid finite difference coefficients can be accurately computed also
within the same framework.

The Taylor series expansion of a function $f(x)$ can be written as

$\begin{cases}f(x+h)=f(x)+\frac{\partial f(x)}{\partial
x}h+\frac{1}{2!}\frac{\partial^{2}f(x)}{\partial
x^{2}}h^{2}+\frac{1}{3!}\frac{\partial^{3}f(x)}{\partial x^{3}}h^{3}+\ldots\\\
f(x-h)=f(x)-\frac{\partial f(x)}{\partial
x}h+\frac{1}{2!}\frac{\partial^{2}f(x)}{\partial
x^{2}}h^{2}-\frac{1}{3!}\frac{\partial^{3}f(x)}{\partial
x^{3}}h^{3}+\ldots\end{cases}.$ (44)

It leads to

$\begin{cases}\frac{f(x+h)+f(x-h)}{2}&=f(x)+\frac{1}{2!}\frac{\partial^{2}f(x)}{\partial
x^{2}}h^{2}+\frac{1}{4!}\frac{\partial^{4}f(x)}{\partial x^{4}}h^{4}+\ldots\\\
\frac{f(x+h)-f(x-h)}{2}&=\frac{\partial f(x)}{\partial
x}h+\frac{1}{3!}\frac{\partial^{3}f(x)}{\partial
x^{3}}h^{3}+\frac{1}{5!}\frac{\partial^{5}f(x)}{\partial
x^{5}}h^{5}+\ldots\end{cases}.$ (45)

Let $h=\Delta x/2$. This implies the 2nd order accuracy of centered finite
difference scheme using only two staggered nodes

$\begin{cases}\frac{\partial f(x)}{\partial x}=\frac{f(x+\Delta
x/2)-f(x-\Delta x/2)}{\Delta x}+O(\Delta x^{2})\\\ f(x)=\frac{f(x+\Delta
x/2)+f(x-\Delta x/2)}{2}+O(\Delta x^{2})\end{cases}.$ (46)

To approximate the 1st order derivatives as accurate as possible, we express
it using more consecutive nodes. Due to regular grid staggering, the
coefficients lying on symmetric positions should have the same coefficients.
This means the first-order derivative reads in the following form

$\begin{split}\frac{\partial f}{\partial x}=&b_{1}\frac{f(x+\Delta
x/2)-f(x-\Delta x/2)}{\Delta x}+\\\ &b_{2}\frac{f(x+3\Delta x/2)-f(x-3\Delta
x/2)}{3\Delta x}+\\\ &b_{3}\frac{f(x+5\Delta x/2)-f(x-5\Delta x/2)}{5\Delta
x}+\cdots.\end{split}$ (47)

Substituting the $f(x+h)$ and $f(x-h)$ with equation 44 for $h=\Delta
x/2,3\Delta x/2,\ldots$ results in

$\begin{split}\frac{\partial f}{\partial x}=&b_{1}\cdot 2\left(\frac{\Delta
x}{2}\frac{\partial f}{\partial x}+\frac{1}{3!}(\frac{\Delta
x}{2})^{3}\frac{\partial^{3}f}{\partial x^{3}}+\cdots\right)/{\Delta x}\\\
&+b_{2}\cdot 2\left(\frac{3\Delta x}{2}\frac{\partial f}{\partial
x}+\frac{1}{3!}(\frac{3\Delta x}{2})^{3}\frac{\partial^{3}f}{\partial
x^{3}}+\cdots\right)/{3\Delta x}\\\ &+b_{3}\cdot 2\left(\frac{5\Delta
x}{2}\frac{\partial f}{\partial x}+\frac{1}{3!}(\frac{5\Delta
x}{2})^{3}\frac{\partial^{3}f}{\partial x^{3}}+\cdots\right)/{5\Delta
x}+\ldots\\\ =&(b_{1}+b_{2}+b_{3}+b_{4}+\cdots)\frac{\partial f}{\partial
x}\\\ &+\frac{\Delta x^{2}}{3!\cdot
2^{2}}(b_{1}+3^{2}b_{2}+5^{2}b_{3}+7^{2}b_{4}+\cdots)\frac{\partial^{3}f}{\partial
x^{3}}\\\ &+\frac{\Delta x^{4}}{5!\cdot
2^{4}}(b_{1}+3^{4}b_{2}+5^{4}b_{3}+7^{4}b_{4}+\cdots)\frac{\partial^{5}f}{\partial
x^{5}}+\cdots.\end{split}$ (48)

Thus, taking first $L$ terms (corresponding to using $2L$ nodes) requires

$\begin{cases}b_{1}+b_{2}+b_{3}+\cdots+b_{L}&=1\\\
b_{1}+3^{2}b_{2}+5^{2}b_{3}+\cdots+(2L-1)^{2}b_{L}&=0\\\
b_{1}+3^{4}b_{2}+5^{4}b_{3}+\cdots+(2L-1)^{4}b_{L}&=0\\\ \cdots&\\\
b_{1}+3^{2L-2}b_{2}+5^{2L-2}b_{3}+\cdots+(2L-1)^{2L-2}b_{L}&=0\\\
\end{cases},$ (49)

which again builds up a Vandermonde-like system

$\displaystyle\underbrace{\begin{bmatrix}1&1&\ldots&1\\\
x_{1}&x_{2}&\ldots&x_{L}\\\ \vdots&&\ddots&\vdots\\\
x_{1}^{L-1}&x_{2}^{L-1}&\ldots&x_{L}^{L-1}\\\
\end{bmatrix}}_{\textbf{V}}\underbrace{\begin{bmatrix}b_{1}\\\ b_{2}\\\
\vdots\\\ b_{L}\\\ \end{bmatrix}}_{\textbf{b}}=\underbrace{\begin{bmatrix}1\\\
0\\\ \vdots\\\ 0\\\ \end{bmatrix}}_{\textbf{z}},$ (50)

in which $x_{i}=(2i-1)^{2}$, $i=1,\cdots,L$. The finite difference weights can
then be easily computed using Golub, (1996, Algorithm 4.6.2), see an
Octave/Matlab script for computing them in Yang, (2014, section 2.5). These
numerically computed weights may be cross-validated with the generic method by
Fornberg, (1988, Table 2). With regular grid spacing $\Delta x$, the weights
$b_{i}$ are dimensionless and can be connected to the coefficients
$\alpha_{i}$ in equation 13 involving a scaling factor $\Delta x$.

## References

  * Björck and Pereyra, (1970) Björck, A., and V. Pereyra, 1970, Solution of Vandermonde systems of equations: Mathematics of computation, 24, 893–903.
  * da Silva et al., (2012) da Silva, N. V., J. V. Morgan, L. MacGregor, and M. Warner, 2012, A finite element multifrontal method for 3D CSEM modeling in the frequency domain: Geophysics, 77, E101–E115.
  * de Hoop, (1996) de Hoop, A. T., 1996, A general correspondence principle for time-domain electromagnetic wave and diffusion fields: Geophysical Journal International, 127, 757–761.
  * Demmel and Koev, (2005) Demmel, J., and P. Koev, 2005, The accurate and efficient solution of a totally positive generalized vandermonde linear system: SIAM Journal on Matrix Analysis and Applications, 27, 142–152.
  * Fornberg, (1988) Fornberg, B., 1988, Generation of finite difference formulas on arbitrarily spaced grids: Mathematics of Computation, 51, 699–706.
  * Golub, (1996) Golub, G. H., 1996, Matrix computation, third edition: Johns Hopkins Studies in Mathematical Sciences.
  * Key, (2016) Key, K., 2016, MARE2DEM: a 2-D inversion code for controlled-source electromagnetic and magnetotelluric data: Geophysical Journal International, 207, 571–588.
  * Komatitsch and Martin, (2007) Komatitsch, D., and R. Martin, 2007, An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation: Geophysics, 72, SM155–SM167.
  * Lee et al., (1989) Lee, K. H., G. Liu, and H. Morrison, 1989, A new approach to modeling the electromagnetic response of conductive media: Geophysics, 54, 1180–1192.
  * Li and Key, (2007) Li, Y., and K. Key, 2007, 2D marine controlled-source electromagnetic modeling: Part 1—an adaptive finite-element algorithm: Geophysics, 72, WA51–WA62.
  * Maaø, (2007) Maaø, F., 2007, Fast finite-difference time-domain modeling for marine subsurface electromagnetic problems: Geophysics, 72, A19–A23.
  * Mittet, (2010) Mittet, R., 2010, High-order finite-difference simulations of marine CSEM surveys using a correspondence principle for wave and diffusion fields: Geophysics, 75, F33–F50.
  * Mittet, (2015) ——–, 2015, Seismic wave propagation concepts applied to the interpretation of marine controlled-source electromagnetics: Geophysics, 80, E63–E81.
  * Mittet, (2017) ——–, 2017, On the internal interfaces in finite-difference schemes: Geophysics, 82, T159–T182.
  * Monk and Süli, (1994) Monk, P., and E. Süli, 1994, A convergence analysis of Yee’s scheme on nonuniform grids: SIAM Journal on Numerical Analysis, 31, 393–412.
  * Mulder, (2006) Mulder, W., 2006, A multigrid solver for 3D electromagnetic diffusion: Geophysical prospecting, 54, 633–649.
  * Newman and Alumbaugh, (1995) Newman, G. A., and D. L. Alumbaugh, 1995, Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences: Geophysical Prospecting, 43, 1021–1042.
  * Oristaglio and Hohmann, (1984) Oristaglio, M. L., and G. W. Hohmann, 1984, Diffusion of electromagnetic fields into a two-dimensional earth: A finite-difference approach: Geophysics, 49, 870–894.
  * Puzyrev et al., (2013) Puzyrev, V., J. Koldan, J. de la Puente, G. Houzeaux, M. Vázquez, and J. M. Cela, 2013, A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling: Geophysical Journal International, 193, 678–693.
  * Rochlitz et al., (2021) Rochlitz, R., M. Seidel, and R.-U. Börner, 2021, Evaluation of three approaches for simulating 3-D time-domain electromagnetic data: Geophysical Journal International, 227, 1980–1995.
  * Rochlitz et al., (2019) Rochlitz, R., N. Skibbe, and T. Günther, 2019, custem: Customizable finite-element simulation of complex controlled-source electromagnetic data: Geophysics, 84, F17–F33.
  * Roden and Gedney, (2000) Roden, J. A., and S. D. Gedney, 2000, Convolution PML (CPML): An efficient FDTD implementation of the CFS–PML for arbitrary media: Microwave and optical technology letters, 27, 334–339.
  * (23) Smith, J. T., 1996a, Conservative modeling of 3-D electromagnetic fields, Part i: Properties and error analysis: Geophysics, 61, 1308–1318.
  * (24) ——–, 1996b, Conservative modeling of 3-D electromagnetic fields, Part II: Biconjugate gradient solution and an accelerator: Geophysics, 61, 1319–1324.
  * Streich, (2009) Streich, R., 2009, 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy: Geophysics, 74, F95–F105.
  * Taflove and Hagness, (2005) Taflove, A., and S. C. Hagness, 2005, Computational electrodynamics: The finite-difference time-domain method, 3rd ed.: Artech House.
  * Wang and Hohmann, (1993) Wang, T., and G. W. Hohmann, 1993, A finite-difference, time-domain solution for three-dimensional electromagnetic modeling: Geophysics, 58, 797–809.
  * Werthmüller, (2017) Werthmüller, D., 2017, An open-source full 3D electromagnetic modeler for 1D VTI media in Python: empymod: Geophysics, 82, WB9–WB19.
  * Werthmüller et al., (2019) Werthmüller, D., W. Mulder, and E. Slob, 2019, emg3d: A multigrid solver for 3d electromagnetic diffusion: Journal of Open Source Software, 4, 1463\.
  * Yang, (2014) Yang, P., 2014, A numerical tour of wave propagation: Technical report, Xi’an Jiaotong University.
  * Yee, (1966) Yee, K. S., 1966, Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media: IEEE Transactions on Antennas and Propagation, 14, 302–307.
  * Yefet and Petropoulos, (2001) Yefet, A., and P. G. Petropoulos, 2001, A staggered fourth-order accurate explicit finite difference scheme for the time-domain maxwell’s equations: Journal of Computational Physics, 168, 286–315.