# 3D fictitious wave domain CSEM inversion by adjoint source estimation

Pengliang Yang1
1 School of Mathematics, Harbin Institute of Technology, 150001, Harbin, China
E-mail<EMAIL_ADDRESS>
###### Abstract

Marine controlled-source electromagnetic (CSEM) method has proved its
potential in detecting highly resistive hydrocarbon bearing formations. A
novel frequency domain CSEM inversion approach using fictitious wave domain
time stepping modelling is presented. Using Lagrangian-based adjoint state
method, the inversion gradient with respect to resistivity can be computed by
the product between the forward and adjoint fields. Simulation of the adjoint
field using the same modelling engine is challenging as it requires time
domain adjoint source time functions while only a few discrete frequencies of
the data residual are available for the inversion. A regularized linear
inverse problem is formulated in order to estimate a long time series from
very few frequency samples. It can then be solved using linear optimization
technique, yielding a matrix-free implementation. Instead of computing adjoint
source time function one by one at each receiver location, a basis function
implementation has been developed such that the inverse problem can be solved
only once and reused every time to construct all time-domain adjoint sources.
The method allows computing all frequencies of the EM fields in one go without
heavy memory and computational overhead, making efficient 3D CSEM inversion
feasible. Numerical examples are employed to demonstrate the application of
our method.

##  1 Introduction

Controlled-source electromagnetic (CSEM) method is a well established
technology to do geophysical exploration (Chave and Cox, 1982; Constable et
al., 1986). It can also be configured for air-borne (Chang-Chun et al., 2015)
and cross-well (Alumbaugh and Newman, 1997) geometries. Land CSEM has been a
commonplace to find mineral deposit (Ward and Hohmann, 1988; Zhdanov and
Keller, 1994; Grayver et al., 2014). As a complement to seismic measurement,
marine CSEM have been successfully applied to detect hydrocarbon bearing
formations in oil and gas industry (Eidesmo et al., 2002; Ellingsrud et al.,
2002; Constable and Srnka, 2007; MacGregor and Tomlinson, 2014), using vessel-
towed dipole source over an array of receivers deployed on the seabed (hence
the name seabed logging). The key to these applications is the good
sensitivity of the electromagnetic signals to distinguish high contrast in
resistivity between saline-filled rocks and hydrocarbons (Abubakar et al.,
2008). This property makes CSEM an ideal tool for de-risking (MacGregor et
al., 2007) in reservoir exploration prior to drilling, as well as 4D
monitoring during the production (Shantsev et al., 2020). A topical review on
marine CSEM inversion has been given by Constable (2010).

CSEM inversion, also known as resistivity tomography, has been a standard
technique to translate CSEM data into electrical properties of the subsurface.
To determine the resistivity model, a nonlinear inverse problem is then
formulated to iteratively minimize the difference between the observed EM data
and the synthetic data derived from numerical modelling, using different
optimization schemes, such as nonlinear conjugate gradient method (Gribenko
and Zhdanov, 2007; Commer and Newman, 2008), quasi-Newton l-BFGS (Plessix and
Mulder, 2008; Schwarzbach and Haber, 2013), and Gauss-Newton method (Constable
et al., 1987; Abubakar et al., 2008; Zaslavsky et al., 2013).

The most computation intensive part of CSEM inversion is the numerical
simulation of 3D electromagnetic field. One can consider the time-domain
finite-difference method (Oristaglio and Hohmann, 1984; Wang and Hohmann,
1993; Taflove and Hagness, 2005), the frequency-domain finite-difference
method (Newman and Alumbaugh, 1995; Smith, 1996a; Mulder, 2006; Streich,
2009), and the frequency-domain finite-element method (Li and Key, 2007; da
Silva et al., 2012; Key, 2016; Rochlitz et al., 2019).

A plethora of CSEM modelling tools directly solve the linear system based on
the frequency domain solution of the discretized Maxwell equation. This avoids
the high computational cost by direct discretization of time-domain diffusive
Maxwell equation involving extremely large number of time steps dictated by
the restrictive stability condition. Both the frequency domain finite
difference method and finite element method formulate Maxwell equation as a
matrix-based linear system, which may be solved using direct (Streich, 2009)
or iterative (Smith, 1996b; Mulder, 2006; Puzyrev et al., 2013) solvers.
Modelling by direct solver is very attractive for multi-source CSEM problems
(Streich, 2009), as the matrix system can be factorized only once and reused
for all sources by forward-backward substitution. Since the inversion of the
large sparse matrix for 3D problems involves huge amount of memory resources,
direct methods may easily go beyond the memory capacity of a desktop computer.
Iterative solvers require much less memory storage but may be time-consuming
and difficult to converge due to ill-conditioning of the discretized Helmholtz
matrix when the conductivity/resistivity model is highly heterogeneous or the
modelling grid is severely stretched (Mulder, 2006). While the frequency
domain finite element method is flexible to address complex model geometries
thanks to the pre-computed meshes, the meshing in 3D geometries using
tetrahedral and hexahedral mesh itself is a challenging and time consuming
task.

In the time-domain, the fields are updated at each time step, allowing the
modelling over the same memory units. The frequency domain fields can be
integrated on the fly during time-stepping, such that multiple frequencies can
be extracted from the same simulation. This avoids repeating several times of
the modelling for different frequencies in frequency domain methods. It
motivates Maaø (2007) to propose a modified wave domain approach to
significantly speed up the computation of diffusive electromagnetic modelling.
The method has been adapted in Støren et al. (2008) to efficiently compute the
inversion gradient for 3D industrial scale applications. The modified wave
domain approach by Maaø (2007) has an attenuation/diffusive term. Inspired by
the fictitious wave domain approach initially proposed in Lee et al. (1989),
Mittet (2010) transformed the diffusion domain Maxwell equation into a pure
wave domain.

A straightforward time-domain discretization of the diffusive Maxwell equation
leads to the temporal sampling proportional to the square of grid spacing, due
to the requirement of the stability condition ($\Delta t\propto\Delta x^{2}$).
This means taking half the spatial sampling will lead to a quadratic increase
of the number of time steps for a simulation of the same duration.
Transforming the diffusive domain into wave domain using fictitious wave
domain approach allows a linear increase of time step with the use of finer
grid spacing ($\Delta t\propto\Delta x$). This yields a highly efficient
scheme to compute the same frequency domain EM field, as it significantly
reduces the required number of time steps. These advantages inspire us to
adopt the fictitious wave domain method for CSEM inversion.

The major contribution of this paper is to develop an efficient frequency
domain CSEM inversion scheme based on the fictitious wave domain modelling.
The gradient of the misfit functional for frequency domain CSEM inversion
requires both a forward and an adjoint field. Efficient computation of the
frequency domain adjoint field by time domain modelling manifests itself as a
major challenge in this development as it requires the time domain adjoint
source time functions at all receiver locations which are unavailable due to
the formulation of CSEM inversion in frequency domain. The key novelty of our
approach is to estimate the time-domain adjoint source based on only a few
discrete frequency samples of the data residual. To do so, we formulate a
regularized linear inverse problem, and then present an efficient, matrix-free
implementation using very few basis functions. The linear inverse problem can
be solved only once and the resulting solution can be reused to derive the
adjoint source time functions at all receiver locations. This novel
development delivers an efficient 3D CSEM inversion methodology in fictitious
wave domain. We finally apply our method to two examples for numerical
demonstration.

##  2 CSEM modelling in fictitious wave domain

Within a three dimensional space $X$, the diffusive Maxwell equations are
governed in the frequency domain

$\begin{cases}\nabla\times
E(\mathbf{x},\omega;\mathbf{x}_{s})-\mathrm{i}\omega\mu(\mathbf{x})H(\mathbf{x},\omega;\mathbf{x}_{s})=M(\mathbf{x}_{s},\omega),\quad\mathbf{x}\in
X,\omega\in\Omega\\\ \nabla\times
H(\mathbf{x},\omega;\mathbf{x}_{s})-\sigma(\mathbf{x})E(\mathbf{x},\omega;\mathbf{x}_{s})=J(\mathbf{x}_{s},\omega),\quad\mathbf{x}\in
X,\omega\in\Omega\end{cases},$ (1)

where the variables $\mathbf{x}$, $t$ and $\omega$ denote space, time and
frequency, respectively. The convention of Fourier transform
$\partial_{t}\leftrightarrow-\mathrm{i}\omega$ has been adopted above. The
electrical and magnetic fields ($E$ and $H$) are vectors consisting of 3
components in x, y and z directions; $J$ and $M$ stand for electrical and
magnetic sources at the location $\mathbf{x}_{s}$, respectively. The magnetic
permeability is $\mu$. The conductivity $\sigma$ is a symmetric $3\times 3$
tensor, i.e., $\sigma_{ij}=\sigma_{ji}$, $i,j\in\\{x,y,z\\}$. In the isotropic
medium, 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$. In the vertical transverse isotropic (VTI) medium,
the diagonal elements are different in horizontal and vertical directions:
$\sigma_{h}:=\sigma_{xx}=\sigma_{yy},\quad\sigma_{v}=\sigma_{zz}$, where
$\sigma_{h}$ and $\sigma_{v}$ stand for horizontal conductivity and vertical
conductivity, respectively. The resistivity is defined as the inverse of the
conductivity, i.e., $\rho_{ij}=1/\sigma_{ij}$.

The above diffusive Maxwell equation can be converted into wave domain
(Mittet, 2010), by defining a fictitious di-electrical permittivity in
equation (1) as $\sigma=2\omega_{0}\varepsilon$, while multiplying the 2nd
equation in (1) with $\sqrt{-\mathrm{i}\omega/2\omega_{0}}$:

$\begin{cases}\nabla\times
E^{\prime}(\mathbf{x},\omega;\mathbf{x}_{s})-\mathrm{i}\omega^{\prime}\mu
H^{\prime}(\mathbf{x},\omega;\mathbf{x}_{s})=M^{\prime}(\mathbf{x}_{s},\omega),\quad\mathbf{x}\in
X,\omega\in\Omega\\\ \nabla\times
H^{\prime}(\mathbf{x},\omega;\mathbf{x}_{s})+\mathrm{i}\omega^{\prime}\varepsilon
E^{\prime}(\mathbf{x},\omega;\mathbf{x}_{s})=J^{\prime}(\mathbf{x}_{s},\omega),\quad\mathbf{x}\in
X,\omega\in\Omega\end{cases}$ (2)

based on the following correspondence relation

$E^{\prime}=E,\;H^{\prime}=\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}H,\;M^{\prime}=M,\;J^{\prime}=\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}J,\;\omega^{\prime}=(1+\mathrm{i})\sqrt{\omega\omega_{0}}.$
(3)

The time domain counterpart of equation (2) reads

$\begin{cases}\nabla\times
E^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})+\mu(\mathbf{x})\partial_{t^{\prime}}H^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})=M^{\prime}(\mathbf{x}_{s},t^{\prime}),\quad\mathbf{x}\in
X,t^{\prime}\in T\\\ \nabla\times
H^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})-\varepsilon(\mathbf{x})\partial_{t^{\prime}}E^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})=J^{\prime}(\mathbf{x}_{s},t^{\prime}),\quad\mathbf{x}\in
X,t^{\prime}\in T,\end{cases}$ (4)

which allows us to do efficient modelling using leap-frog finite-difference
time-domain (FDTD) method over the staggered grid. From the electromagnetic
fields in the time domain, the frequency domain fields can be integrated
during modelling using the transformation

$u(\mathbf{x},\omega;\mathbf{x}_{s})=\int_{0}^{T_{\max}}u(\mathbf{x},t^{\prime};\mathbf{x}_{s})e^{\mathrm{i}\omega^{\prime}t^{\prime}}\mathrm{d}t^{\prime}=\int_{0}^{T_{\max}}u(\mathbf{x},t^{\prime};\mathbf{x}_{s})e^{-\sqrt{\omega\omega_{0}}t^{\prime}}e^{\mathrm{i}\sqrt{\omega\omega_{0}}t^{\prime}}\mathrm{d}t^{\prime},$
(5)

where
$u(\mathbf{x},t^{\prime};\mathbf{x}_{s})\in\\{E^{\prime}_{x},E^{\prime}_{y},E^{\prime}_{z},H^{\prime}_{x},H^{\prime}_{y},H^{\prime}_{z}\\}$
denotes the electric/magnetic field at the spatial location $\mathbf{x}\in X$
and the fictitious time $t^{\prime}\in[0,T_{\max}]$ excited by a source from
the spatial location $\mathbf{x}_{s}$; $T_{\max}$ is the final time that the
frequency domain field $u(\mathbf{x},\omega;\mathbf{x}_{s})$ reaches its
steady state. Equation (5) shows that the complex-valued frequency
$\omega^{\prime}=(1+\mathrm{i})\sqrt{\omega\omega_{0}}$ leads to an
exponentially decay factor in the time integration, which is the key to
realize the attenuation effect of diffusive EM field during time evolution.

To model the CSEM response as precise as possible, we resort to our newly
developed high-order finite-difference time-domain method over non-uniform
grid (Yang and Mittet, 2023; Yang, 2023). To satisfy homogeneous boundary
condition approximately, perfectly matched layers (PML) (Komatitsch and
Martin, 2007) are padded surrounding the domain of interest to absorb the
reflections in the truncated domain. The implementation of (5) can be carried
out on the fly during timestepping modelling thanks to the discrete time
Fourier transform (DTFT):

$u(\mathbf{x},\omega_{k};\mathbf{x}_{s})=\sum_{n=0}^{N_{t}-1}u(\mathbf{x},t^{\prime}_{n};\mathbf{x}_{s})\exp(-\sqrt{\omega_{k}\omega_{0}}t^{\prime}_{n})\exp(\mathrm{i}\sqrt{\omega_{k}\omega_{0}}t^{\prime}_{n}),$
(6)

where the discretized fictitious time is $t^{\prime}_{n}=n\Delta t$ with time
step $\Delta t$ and time index $n=0,\cdots,N_{t}$, while the interested
frequencies are $\omega_{k}$, $k=1,\cdots,N_{\omega}$. Note that the total
number of time steps $N_{t}$ is much large than the total number of discrete
frequencies $N_{\omega}$, i.e., $N_{t}\gg N_{\omega}$. The modelling
terminates when the electric field
$E^{\prime}(\mathbf{x},\omega;\mathbf{x}_{s})$ and the magnetic field
$H^{\prime}(\mathbf{x},\omega;\mathbf{x}_{s})$ evolves to the stead state.
This implies that we may check the convergence of the field regularly during
the modelling to avoid additional timesteppings which have negligible
contributions to the time integral. It is noteworthy that multiple frequencies
can be integrated on the fly during the same timestepping procedure. These
advantages make the method attractive for forward modelling.

##  3 Fictitious wave domain CSEM inversion

The CSEM inversion is an iterative nonlinear optimization procedure. At each
iteration, the gradient of the data misfit with respect to conductivity can be
computed via the product between forward and adjoint fields. To use fictitious
wave domain modelling engine for adjoint simulation, we present a novel
approach to estimating the adjoint source time functions by formulating a
regularized linear inverse problem. An in-depth analysis shows that this
linear inverse problem needs to be solved only once to result in a limited
number of basis functions, which can then be used to construct the time-domain
adjoint source at different receiver locations. This key development enables
efficient and practical 3D CSEM inversion by fictitous wave domain modelling
engine.

###  3.1 CSEM inverse problem

The Maxwell equations forms a linear system as

$\underbrace{\begin{bmatrix}-\sigma&\nabla\times\\\
\nabla\times&-\mathrm{i}\omega\mu\end{bmatrix}}_{\mathbf{A}(m)}\underbrace{\begin{bmatrix}E\\\
H\end{bmatrix}}_{\mathbf{u}}=\underbrace{\begin{bmatrix}J\\\
M\end{bmatrix}}_{\mathbf{f}},$ (7)

where the electromagnetic propagator $\mathbf{A}(m)$ is a linear operator
applied to the vector field $\mathbf{u}$, which gathers electrical and
magnetic fields (each component has been denoted by
$u:=u(\mathbf{x},\omega;\mathbf{x}_{s})$ in the previous section). The source
vector with nonzero excitation at the source location $\mathbf{x}_{s}$ is
prescribed by $\mathbf{f}:=\mathbf{f}(\mathbf{x},\omega;\mathbf{x}_{s})$. Note
that all field variables are functions of both frequency and space
($u(\mathbf{x},\omega;\mathbf{x}_{s}),\mathbf{x}\in X,\omega\in\Omega$), while
the medium property is only a function of space ($m(\mathbf{x}),\mathbf{x}\in
X$), which can be the conductivity $\sigma_{ij}$ or the permeability $\mu$.

To find the resistivity of the subsurface, the data misfit is defined in
least-squares sense

$\phi_{d}(m)=\frac{1}{2}\Big{\|}\underbrace{\begin{bmatrix}\mathbf{W}_{1}&0\\\
0&\mathbf{W}_{2}\end{bmatrix}}_{\mathbf{W}}\Big{(}\underbrace{\begin{bmatrix}E^{obs}(\mathbf{x}_{r},\omega;\mathbf{x}_{s})\\\
H^{obs}(\mathbf{x}_{r},\omega;\mathbf{x}_{s})\end{bmatrix}}_{\mathbf{d}(\mathbf{x}_{r},\omega;\mathbf{x}_{s})}-\underbrace{\begin{bmatrix}E(\mathbf{x}_{r},\omega;\mathbf{x}_{s})\\\
H(\mathbf{x}_{r},\omega;\mathbf{x}_{s})\end{bmatrix}}_{\mathbf{d}_{syn}(\mathbf{x}_{r},\omega;\mathbf{x}_{s})[m]}\Big{)}\Big{\|}^{2}=\frac{1}{2}\|\mathbf{W}(\mathbf{d}-\mathbf{R}\mathbf{u})\|^{2},$
(8)

where $\mathbf{d}(\mathbf{x}_{r},\omega;\mathbf{x}_{s})$ denotes the observed
EM data at the receiver location $\mathbf{x}_{r}$ due to the source at the
location $\mathbf{x}_{s}$, while the synthetic data
$\mathbf{d}_{syn}(\mathbf{x}_{r},\omega;\mathbf{x}_{s})[m]:=\mathbf{u}(\mathbf{x}_{r},\omega;\mathbf{x}_{s})=\mathbf{R}\mathbf{u}(\mathbf{x},\omega;\mathbf{x}_{s})$
(simulated with the model parameter $m$) are extracted by the restriction
operator $\mathbf{R}$ from the modelled wavefield at the receiver location

$u(\mathbf{x}_{r},\omega;\mathbf{x}_{s})=\int_{X}u(\mathbf{x},\omega;\mathbf{x}_{s})\delta(\mathbf{x}-\mathbf{x}_{r})\mathrm{d}\mathbf{x},\quad\mathbf{x}\in
X,\omega\in\Omega$ (9)

where $u(\mathbf{x}_{r},\omega;\mathbf{x}_{s})$ is one component of the vector
field $\mathbf{u}(\mathbf{x}_{r},\omega;\mathbf{x}_{s})$. In the remainder of
the paper, we shall drop the dependence of
$(\mathbf{x}_{r},\omega;\mathbf{x}_{s})$ and the summation over sources and
receivers without loss of clarity unless clearly stated if necessary. Two
weighting matrices, $\mathbf{W}_{1}$ and $\mathbf{W}_{2}$ together forming the
diagonal weighting matrix $\mathbf{W}$, are employed to weight the electric
and magnetic fields respectively. These weighting matrices may be specified
according to the uncertainty model in terms of the real acquisition and the
equipment (Mittet and Morten, 2012). If the magnetic data are not considered
in the inversion, we simply set $\mathbf{W}_{2}=0$.

Taking the first derivative of equation (7) with respect to model parameter
$m$ gives

$\frac{\partial\mathbf{u}}{\partial
m}=-\mathbf{A}^{-1}(m)\frac{\partial\mathbf{A}(m)}{\partial m}\mathbf{u}.$
(10)

The gradient of the data misfit with respect to the model parameter $m$ is

$\begin{split}\frac{\partial\phi_{d}(m)}{\partial
m}=&\Re\langle\mathbf{W}\mathbf{R}\frac{\partial\mathbf{u}}{\partial
m},\mathbf{W}(\mathbf{R}\mathbf{u}-\mathbf{d})\rangle\\\
=&-\Re\langle\mathbf{W}\mathbf{R}\mathbf{A}^{-1}(m)\frac{\partial\mathbf{A}(m)}{\partial
m}\mathbf{u},\mathbf{W}(\mathbf{R}\mathbf{u}-\mathbf{d})\rangle\\\
=&\Re\langle\frac{\partial\mathbf{A}(m)}{\partial
m}\mathbf{u},\underbrace{(\mathbf{A}^{\dagger})^{-1}(m)\mathbf{R}^{\dagger}\mathbf{W}^{\dagger}\mathbf{W}(\mathbf{d}-\mathbf{R}\mathbf{u}}_{\mathbf{v}})\rangle\\\
=&\Re\langle\frac{\partial\mathbf{A}(m)}{\partial
m}\mathbf{u},\mathbf{v}\rangle\\\
=&\Re\sum_{s}\sum_{\omega}\bar{\mathbf{v}}^{\mathrm{T}}\frac{\partial\mathbf{A}(m)}{\partial
m}\mathbf{u},\end{split}$ (11)

where $\Re$ takes the real part of a complex number, $\dagger$ is the complex
conjugate transpose. It is noteworthy that the newly introduced variable
$\mathbf{v}$, coined adjoint variable, or co-state variable in optimal control
theory, must satisfy

$\mathbf{A}^{\dagger}(m)\mathbf{v}=\mathbf{R}^{\mathrm{T}}\mathbf{W}^{\mathrm{T}}\mathbf{W}(\mathbf{d}-\mathbf{R}\mathbf{u}).$
(12)

This shows that $\mathbf{v}$ satisfies another Maxwell equation based on
adjoint Maxwell operator $\mathbf{A}^{\dagger}(m)$. The right hand side of the
adjoint equation acts as the virtual source to emanate the adjoint field.

Equation (11) is a generic expression applicable in fully anisotropic medium
for the model parameters $m\in\\{\sigma_{ij},\mu\\}$. Let us point out that
the above gradient expression may be scaled with the local cell volume $\Delta
V(\mathbf{x})$ due to the discretization of the spatial integral
$\int_{X}\mathrm{d}\mathbf{x}$. In the CSEM settings, the magnetic
permeability $\mu$ is considered as the same constant as in the vacuum. We are
interested in retrieving the conductivity $\sigma$ (or the resistivity
$\rho=1/\sigma$) which may exhibit anisotropy. Equation (7) implies

$\frac{\partial\mathbf{A}(m)}{\partial\sigma}=\frac{\partial}{\partial\sigma}\begin{bmatrix}-\sigma&\nabla\times\\\
\nabla\times&-\mathrm{i}\omega\mu\end{bmatrix}=\begin{bmatrix}-1&0\\\
0&0\end{bmatrix}.$ (13)

Denote the adjoint field
$\mathbf{v}=(\underline{E},\underline{H})^{\mathrm{T}}$, where the underline
is utilized to distinguish it from the forward field. The gradient in equation
(11) becomes

$\frac{\partial\phi_{d}}{\partial\sigma_{ij}}=-\Re\sum_{s}\sum_{\omega}\underline{\bar{E}}_{j}\cdot
E_{i},\qquad i,j\in\\{x,y,z\\}.$ (14)

For commonly used VTI medium, the gradients of the misfit with respect to
horizontal conductivity $\sigma_{h}$ and vertical conductivity $\sigma_{v}$
are

$\frac{\partial\phi_{d}}{\partial\sigma_{h}}=-\Re\sum_{s}\sum_{\omega}(\underline{\bar{E}}_{x}\cdot
E_{x}+\underline{\bar{E}}_{y}\cdot
E_{y}),\quad\frac{\partial\phi_{d}}{\partial\sigma_{v}}=-\Re\sum_{s}\sum_{\omega}\underline{\bar{E}}_{z}\cdot
E_{z}.$ (15)

In isotropic medium ($\sigma_{ii}=\sigma$, $i=x,y,z$), we have the simple
gradient expression

$\frac{\partial\phi_{d}}{\partial\sigma}=-\Re\sum_{s}\sum_{\omega}\underline{E}^{\dagger}\cdot
E=-\Re\sum_{s}\sum_{\omega}(\underline{\bar{E}}_{x}\cdot
E_{x}+\underline{\bar{E}}_{y}\cdot E_{y}+\underline{\bar{E}}_{z}\cdot E_{z}).$
(16)

###  3.2 Adjoint modelling in fictitious wave domain

To facilitate the computation of the adjoint field and the inversion gradient,
one may take the conjugate of equation (12):

$\mathbf{A}^{\mathrm{T}}(m)\bar{\mathbf{v}}=\overline{\mathbf{R}^{\mathrm{T}}\mathbf{W}^{\mathrm{T}}\mathbf{W}(\mathbf{d}-\mathbf{R}\mathbf{u})}.$
(17)

Since $\mathbf{A}^{\mathrm{T}}(m)=\mathbf{A}(m)$, the above operation allows
us to use the same modelling engine to compute the conjugate of the adjoint
field, as long as proper adjoint source can be provided for adjoint
simulation. Unfortunately, for CSEM inversion in the frequency domain, only
few discrete frequencies are available, leading to the adjoint source (i.e.
the data residual) specified only at few discrete frequencies. A time domain
adjoint source is therefore not directly accessible.

To still use fictitious wave domain timestepping modelling for adjoint
simulation, let us now split the weighted data residual into electrical and
magnetic components
$\overline{\mathbf{R}^{\mathrm{T}}\mathbf{W}^{\mathrm{T}}\mathbf{W}(\mathbf{d}-\mathbf{R}\mathbf{u})}=(\overline{\delta
d_{E}},\overline{\delta d_{H}})^{\mathrm{T}}$ and repeat the substitution and
multiplication operations as in Section 2:

$\begin{cases}\nabla\times\underbrace{\overline{\underline{E}}}_{\underline{\overline{E}}^{\prime}}+\mu\underbrace{\sqrt{-\mathrm{i}2\omega\omega_{0}}}_{-\mathrm{i}\omega^{\prime}}\underbrace{\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}\underline{\overline{H}}}_{\underline{\overline{H}}^{\prime}}=\overline{\delta
d_{E}},\\\
\nabla\times\underbrace{\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}\underline{\overline{H}}}_{\underline{\overline{H}}^{\prime}}-\underbrace{\sqrt{-\mathrm{i}2\omega\omega_{0}}}_{-\mathrm{i}\omega^{\prime}}\varepsilon\underline{\overline{E}}=\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}\overline{\delta
d_{H}}.\end{cases}$ (18)

We need to consider the time domain counterpart of the above system for
efficient adjoint modelling by fictitious timestepping modelling. Switching
from diffusive frequency domain back to fictitious time domain with wave
should be understood as a linear transformation rather than the true inverse
Fourier transform. This is because the mapping of fictitious transform is not
bijective, as the number of time steps $N_{t}$ for numerical simulation is
significantly larger than the number of discrete frequencies $N_{\omega}$ used
for CSEM investigation ($N_{\omega}\ll N_{t}$). There exists an infinite
number of time series which may match the few discrete frequencies. Searching
for the inverse is therefore an under-determined problem.

Let us denote the right hand side of equation (18) ($\overline{\delta d_{E}}$
or $\sqrt{\frac{-\mathrm{i}\omega}{2\omega_{0}}}\overline{\delta d_{H}}$) at a
specific receiver location as $s(\omega)$. To retrieve a long time series for
adjoint modelling, a linear inverse problem is formulated to convert the
adjoint source $s(\omega)$ from frequency to time domain based on DTFT in
equation (6) by minimizing the following misfit functional

$\begin{split}\Psi&=\sum_{k=1}^{N_{\omega}}|s(\omega_{k})-\sum_{n=0}^{N_{t}-1}e^{-\sqrt{\omega_{k}\omega_{0}}t^{\prime}_{n}}e^{\mathrm{i}\sqrt{\omega_{k}\omega_{0}}t^{\prime}_{n}}s(t^{\prime}_{n})|^{2}+\gamma\sum_{n=0}^{N_{t}-1}|s(t^{\prime}_{n})|^{2}\\\
&=\sum_{k=1}^{N_{\omega}}\Big{(}|\Re\\{s(\omega_{k})\\}-\sum_{n=0}^{N_{t}-1}e^{-\sqrt{\omega_{k}\omega_{0}}t^{\prime}_{n}}\cos(\sqrt{\omega_{k}\omega_{0}}t^{\prime}_{n})s(t^{\prime}_{n})|^{2}\\\
&\qquad+|\Im\\{s(\omega_{k})\\}-\sum_{n=0}^{N_{t}-1}e^{-\sqrt{\omega_{k}\omega_{0}}t^{\prime}_{n}}\sin(\sqrt{\omega_{k}\omega_{0}}t^{\prime}_{n})s(t^{\prime}_{n})|^{2}\Big{)}+\gamma\sum_{n=0}^{N_{t}-1}|s(t^{\prime}_{n})|^{2}\end{split}$

where $\Im$ takes the imaginary part of the complex variable. If we discretize
the time as $t_{n}=n\Delta t$ and define

$\mathbf{s}(\omega):=\begin{bmatrix}\Re\\{s(\omega_{1})\\}\\\ \vdots\\\
\Re\\{s(\omega_{N_{\omega}})\\}\\\ \Im\\{s(\omega_{1})\\}\\\ \vdots\\\
\Im\\{s(\omega_{N_{\omega}})\\}\end{bmatrix},\mathbf{s}(t^{\prime}):=\begin{bmatrix}s(t^{\prime}_{0})\\\
\vdots\\\ s(t^{\prime}_{N_{t}-1})\end{bmatrix}$ (19)

and

$\mathbf{B}:=\begin{bmatrix}e^{-\sqrt{\omega_{1}\omega_{0}}t^{\prime}_{0}}\cos(\sqrt{\omega_{1}\omega_{0}}t^{\prime}_{0})&\cdots&e^{-\sqrt{\omega_{1}\omega_{0}}t^{\prime}_{N_{t}-1}}\cos(\sqrt{\omega_{1}\omega_{0}}t^{\prime}_{N_{t}-1})\\\
\vdots&\ddots&\vdots\\\
e^{-\sqrt{\omega_{N_{\omega}}\omega_{0}}t^{\prime}_{0}}\cos(\sqrt{\omega_{N_{\omega}}\omega_{0}}t^{\prime}_{0})&\cdots&e^{-\sqrt{\omega_{N_{\omega}}\omega_{0}}t^{\prime}_{N_{t}-1}}\cos(\sqrt{\omega_{N_{\omega}}\omega_{0}}t^{\prime}_{N_{t}-1})\\\
e^{-\sqrt{\omega_{1}\omega_{0}}t^{\prime}_{0}}\sin(\sqrt{\omega_{1}\omega_{0}}t^{\prime}_{0})&\cdots&e^{-\sqrt{\omega_{1}\omega_{0}}t^{\prime}_{N_{t}-1}}\sin(\sqrt{\omega_{1}\omega_{0}}t^{\prime}_{N_{t}-1})\\\
\vdots&\ddots&\vdots\\\
e^{-\sqrt{\omega_{N_{\omega}}\omega_{0}}t^{\prime}_{0}}\sin(\sqrt{\omega_{N_{\omega}}\omega_{0}}t^{\prime}_{0})&\cdots&e^{-\sqrt{\omega_{N_{\omega}}\omega_{0}}t^{\prime}_{N_{t}-1}}\sin(\sqrt{\omega_{N_{\omega}}\omega_{0}}t^{\prime}_{N_{t}-1})\\\
\end{bmatrix},$ (20)

the objective of the above linear inverse problem can be compactly written in
matrix form

$\Psi=\|\mathbf{s}(\omega)-\mathbf{B}\mathbf{s}(t^{\prime})\|^{2}+\gamma\|\mathbf{s}(t^{\prime})\|^{2},$
(21)

where a regularization term taking into account the minimum energy of the
solution has been penalized by a parameter $\gamma$ to stabilize the linear
inversion and to determine a unique solution. Note that the complex factors
have been split into real and imaginary part. To solve this ill-posed problem,
regularization has been introduced to recover the well-posedness during the
inversion of the matrix $\mathbf{B}\in\mathbb{R}^{2N_{\omega}\times N_{t}}$.
The solution $\mathbf{s}(t)$ can easily be found using linear optimization
algorithms, e.g. LSQR (Paige and Saunders, 1982) and the conjugate gradient
method for the resulting normal equation (CGNR) (Saad, 2003). Based on the
frequency and time index, the each element of $\mathbf{B}$ can be formed on
the fly when computing matrix vector product, implying a matrix free
implementation.

It should be noted that fictitious wave domain modelling is simply a
mathematical tool to efficiently compute frequency domain EM fields.
Disguising a diffusive phenomenon as a wave event will not change the
diffusive nature of the underlying physics. Indeed, Mittet (2010) proposed
another method to calculate the adjoint source time function by superposition
of delayed causal wavelets. However, it yields a nonlinear inverse problem
which is very difficult to solve. Clearly, embedding a linear inverse problem
in a nonlinear inverse problem is a much better option than embedding double
nonlinear inverse problems. The truncated singular value decomposition (SVD)
have been applied (Støren et al., 2008) to invert for an adjoint source,
following the modified wave domain formulation of Maaø (2007). Mathematically,
our approach using damped least-squares should give the same solution as the
one from truncated SVD, since SVD is equivalent to finding the Moore-Penrose
pseudo-inverse for $\mathbf{B}$ (Björck, 1996, theorem 1.2.10 in p. 15, and
section 2.7.2, p. 101). However, our empirical experience shows that the
approximate inverse given by SVD suffers from numerical instability while
consuming too much memory. The above linear inversion method gets rid of these
issues, incorporating all frequencies in one run with better memory and
computational efficiency.

###  3.3 A matrix-free basis function implementation

Let us denote
$\mathbf{B}^{+}=(\mathbf{B}^{\mathrm{T}}\mathbf{B}+\gamma\mathbf{I})^{-1}\mathbf{B}^{\mathrm{T}}\in\mathbb{R}^{N_{t}\times
2N_{\omega}}$ the Moore-Penrose pseudo-inverse of
$\mathbf{B}\in\mathbb{R}^{2N_{\omega}\times N_{t}}$ such that

$\mathbf{s}(t)=\mathbf{B}^{+}\mathbf{s}(\omega).$ (22)

The matrix $\mathbf{B}^{+}$ may be written down using the $2N_{\omega}$
columns:

$\mathbf{B}^{+}=\begin{bmatrix}\mathbf{b}_{1},\cdots,\mathbf{b}_{N_{\omega}},\mathbf{b}_{N_{\omega}+1},\cdots,\mathbf{b}_{2N_{\omega}}\end{bmatrix}$
(23)

in which $\mathbf{b}_{k}\in\mathbb{R}^{N_{t}}$ is the $i$-th column of
$\mathbf{B}^{+}$. As a result, equation (22) translates into

$\mathbf{s}(t)=\sum_{k=1}^{N_{\omega}}\Re\\{s(\omega_{k})\\}\mathbf{b}_{k}+\Im\\{s(\omega_{k})\\}\mathbf{b}_{k+N_{\omega}}.$
(24)

It then becomes clear that the columns of $\mathbf{B}^{+}$ are the very few
number of basis functions to construct the adjoint source time function. Given
$\Delta t$ and $N_{t}$, these basis functions are uniquely determined, hence
independent of the frequency spectrum $s(\omega)$ (which are simply the
coefficients for the linear combination of these basis functions). This means
the iterative solution procedure needs to be performed only once to find all
$\mathbf{b}_{k}$, which can then be used for each receiver locations to infer
the corresponding adjoint source time function. Equation (24) also gives a
recipe to efficiently compute these functions:

  * •

If we set $\Re\\{s(\omega_{k})\\}=1$, $\Im\\{s(\omega_{k})\\}=0$ and
$\Re\\{s(\omega_{j})\\}=\Im\\{s(\omega_{j})\\}=0$ ($j\neq k$), we obtain
$\mathbf{s}(t)=\mathbf{b}_{k}$.

  * •

If we set $\Im\\{s(\omega_{k})\\}=1$, $\Re\\{s(\omega_{k})\\}=0$ and
$\Re\\{s(\omega_{j})\\}=\Im\\{s(\omega_{j})\\}=0$ ($j\neq k$), we obtain
$\mathbf{s}(t)=\mathbf{b}_{k+N_{\omega}}$.

By repeatedly feeding the iterative optimization algorithm $2N_{\omega}$ times
using different Dirac delta input
$\mathbf{s}(\omega)=(0,\cdots,0,1,0,\cdots,0)^{\mathrm{T}}$, we obtain the
$2N_{\omega}$ columns of $\mathbf{B}^{+}$.

Figure 1 gives an example of these basis functions using the parameters
$\Delta t=0.002$, $N_{t}=1000$ and $f_{0}=1$ and $N_{\omega}=3$. To solve this
linear optimization problem, Claerbout’s conjugate gradient algorithm
(Claerbout and Fomel, 2008, chapter 2.3.6) (which is essentially a CGNR
algorithm according to Saad (2003, chapter 8.3, pp. 266-268)) has been
applied. As can be seen from Figure 2, within 15 CG iterations, the proposed
method leads to highly accurate estimation for these basis functions, with the
error less than $10^{-8}$. In this test, the regularization parameter $\gamma$
was chosen to be $10^{-3}$. It is found the shape of the basis functions do
not change much using different choices of $\gamma$, indicating that the
method is robust.

For different receiver locations, the data residuals in the frequency domain
may be dramatically different, but these basis functions are the same. We
therefore only solve the above regularized linear inverse problem once to
construct all adjoint source time functions through their linear combination.
The frequency domain data residuals are the weighting coefficients to modulate
these basis.

Figure 1: The basis functions using the parameters $\Delta t=0.004$ s,
$N_{t}=1000$ and $N_{\omega}=3$. The amplitude of these basis functions decays
to zero with the increase of the time, which is quite suitable for frequency
domain EM fields converging to steady state after simulation using sufficient
number of time steps. Figure 2: The convergence of CG iterations for the
solution of basis functions

One thing to remark is that the adjoint modelling in the time domain should
have reverse time order compared with the forward modelling. Taking the
conjugate reverse the time again. This is the reason why a decaying trend of
the estimated basis function is observed in Figure 1.

###  3.4 Connections and distinctions with existing methods

Note the incident field $E_{i}$ is the product between the Green’s function
with the source current

$E_{i}(\mathbf{x},\omega;\mathbf{x}_{s})=\sum_{k}G_{ik}^{EE}(\mathbf{x},\omega;\mathbf{x}_{s})J_{k}(\mathbf{x}_{s},\omega),$
(25)

while the conjugate of adjoint field can be specified as

$\overline{\underline{E}}_{j}(\mathbf{x},\omega;\mathbf{x}_{s})=\int_{X}(\sum_{k}G_{jk}^{EE}(\mathbf{x},\omega;\mathbf{x}_{r})\overline{\delta
d^{E}_{k}(\mathbf{x}_{r},\omega;\mathbf{x}_{s})}+\sum_{i}G_{ji}^{EH}(\mathbf{x},\omega;\mathbf{x}_{r})\overline{\delta
d^{H}_{i}(\mathbf{x}_{r},\omega;\mathbf{x}_{s})})\mathrm{d}\mathbf{x}_{r},$
(26)

where $i,j,k\in\\{x,y,z\\}$ indicate different components of the EM field,
while $G_{ij}^{EE}(\mathbf{x},\omega;\mathbf{x}_{s})$ and
$G_{ik}^{EH}(\mathbf{x},\omega;\mathbf{x}_{s})$ stand for the $i$th component
of the electrical Green’s function for angular frequency $\omega$ at the
spatial location $\mathbf{x}$ due to the electrical and magnetic sources
directed in $j$\- and $k$\- th directions. This implies that we can still use
fictitious wave modelling for adjoint simulation: one can, for each frequency
and each receiver, first compute the Green’s function (which is independent of
the source time function) in frequency domain, and then form the adjoint field
by linear combination of them using the data residual. This method, ensured by
the superposition principle thanks to the linearity of the Maxwell equation,
enables the approach still working out. Such a computing scheme repeats many
times of simulation depending on the number of receivers, which is very
inefficient when the number of receivers becomes extremely large.

Equation (17) shows that the adjoint field can be computed using only one
modelling, provided that a time-domain adjoint source time function is
available and then injected at once to do time-domain simulation. The spirit
of our method is to perform frequency domain inversion using time domain
modelling engine, similar to Sirgue et al. (2008, 2010) for seismic full
waveform inversion (FWI). However, there are a number of differences.
Switching between time and frequency domain is straightforward in seismic FWI
using the definition of discrete inverse Fourier transform. One can convert
frequency domain data residual using the definition of inverse discrete
Fourier transform

$s(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}s(\omega)e^{-i\omega
t}\mathrm{d}\omega\approx\sum_{k=1}^{N_{\omega}}s(\omega_{k})e^{-i\omega_{k}t}.$
(27)

Of course, the same issue presented in Sirgue et al. (2010) approach: only
limited number of frequencies are available to reconstruct a long time series.
The last equality in (27) assumes all the absent frequencies are zeros. Since
Fourier basis is orthonormal, the computation of adjoint source time function
under this assumption leads to a minimum energy solution equivalent to solving
the following least-squares minimization problem

$\min_{s(t)}\sum_{k=1}^{N_{\omega}}\|s(\omega_{k})-\sum_{j=1}^{N_{t}}s(t_{j})e^{\mathrm{i}\omega_{k}t_{j}}\|^{2}.$
(28)

For fictitious wave domain approach, plugging the complex-valued frequency
$\omega^{\prime}=(1+\mathrm{i})\sqrt{\omega\omega_{0}}$ into exponential
factor yields an exponentially decay factor

$e^{\mathrm{i}\omega^{\prime}t^{\prime}}=e^{\mathrm{i}(1+\mathrm{i})\sqrt{\omega\omega_{0}}t^{\prime}}=e^{-\sqrt{\omega\omega_{0}}t^{\prime}}e^{\mathrm{i}\sqrt{\omega\omega_{0}}t^{\prime}},$
(29)

leaving the converted EM field possessing strong attenuation behavior, while
making the basis no more orthogonal. Precisely speaking, the transformation in
(5) is a Laplace transform rather than normal Fourier transform. As a result,
the inverse fictitious transformation is not well defined, that is,

$u(\mathbf{x},t^{\prime};\mathbf{x}_{s})\neq\frac{1}{2\pi}\int_{\Omega}\mathrm{d}\omega
u(\mathbf{x},\omega;\mathbf{x}_{s})e^{-\mathrm{i}\omega^{\prime}t^{\prime}}.$
(30)

Since direct use of inverse Fourier transform as done by Sirgue et al. (2010)
does not apply, a matrix-free iterative solution of adjoint source time
function in this work is evidently desirable.

###  3.5 The final algorithm

The total misfit functional for CSEM inversion consists of both data fitting
and model regularization:

$\phi(m)=\phi_{d}(m)+\beta\phi_{m}(m),$ (31)

where $\phi_{m}(m)$, penalized by the parameter $\beta$, is the model misfit
term to enforce smoothness on the inverted model. A popular choice is to use
Tikhonov regularization to minimize the roughness of the model compared with a
reference model $m_{ref}$, i.e.,
$\phi_{m}(m)=\frac{1}{2}\|\nabla^{\alpha}(m-m_{ref})\|^{2}$, where the
anisotropic first order difference operator
$\nabla^{\alpha}=(\alpha_{x}\partial_{x},\alpha_{y}\partial_{y},\alpha_{z}\partial_{z})^{\mathrm{T}}$
with coefficients $\alpha_{x}$, $\alpha_{y}$ and $\alpha_{z}$ can be tuned by
the user. This leads to the total gradient of the misfit functional prescribed
by $\frac{\partial\phi}{\partial m}=\frac{\partial\phi_{d}}{\partial
m}+\frac{\partial\phi_{m}}{\partial m}$. Based on the gradient information
computed above, a descent direction $\delta m^{k}$ can be constructed to
update the model parameters iteratively

$m^{k+1}=m^{k}+\alpha\delta m^{k},$ (32)

where $\alpha$ is the step length estimated by line search method. In this
paper, the descent direction $\delta m^{k}$ at the $k$-th iteration is
computed using l-BFGS algorithm (Nocedal and Wright, 2006) by storing the
gradients in the previous iterations.

It is worth noting that the inversion gradient computed using the proposed
approach is approximate rather than exact. Consequently, the descent direction
estimated using l-BFGS after several iterations may have difficulty to succeed
in the line search procedure. We therefore propose to restart l-BFGS using the
steepest descent direction after the first failure of the line search. The
algorithm will eventually be terminated if the maximum number of iterations is
reached.

Using conductivity as model parameter creates some kinds of ill-conditioning
due to its value varying in a large range (the seawater is around 0.3
$\Omega\cdot$m while a thin layer of hydrocarbon bearing sediment is more than
100 $\Omega\cdot$m). To capture the high contrast of the parameter variations,
the inversion is re-parameterized using the logarithmic transformation
$m:=\ln\rho$, which is a dimensionless parameter with comparable magnitude
varying in a much smaller dynamic range. Indeed, there exists other types of
parameter scaling, see another example in Abubakar et al. (2008). This re-
parametrization gives a relation $\rho=e^{m}$. Switching the parametrization
is trivial thanks to the chain rule:

$\frac{\partial\phi}{\partial m}=\frac{\partial\rho}{\partial
m}\cdot\frac{\partial\sigma}{\partial\rho}\cdot\frac{\partial\phi}{\partial\sigma}=-\frac{1}{\rho}\frac{\partial\phi}{\partial\sigma}.$
(33)

During the iterative inversion, the model parameter may be bounded within an
upper bound $m_{\max}$ and a lower bound $m_{\min}$ according to a priori
knowledge. These bounds help to stabilize the inversion in a more physically
sensible manner.

##  4 Application examples

In this section we present two numerical examples to demonstrate the
application of our method. The two examples are designed to validate the
gradient expressions for the isotropic (the first example) and the VTI
anisotropic (the second example) cases. In both examples, I use in-line and
azimuth data with $E_{x}$ and $E_{y}$ components. This motivates us to utilize
the uncertainty model proposed in Morten et al. (2009) to compute the diagonal
weight matrix $\mathbf{W}$ in order to capture the varying sensitivity of
broadside EM data. The observed data are poluted by 3% of Gaussian white noise
before inversion. The regularization parameters are configured with
$\beta=0.01$ using a cooling factor (between 0.7 and 1) through iterations. We
set $\alpha_{x}=\alpha_{y}=1$ and $\alpha_{z}=0.1$ in both tests. We choose
the initial model as the reference model for both inverse exercises. To
mitigate the weak sensitivity of the EM field to the deeper part of the model,
a depth preconditioning has been applied following the work of Plessix and
Mulder (2008). The code is parallelized over the sources to achieve the best
scaling performance since the sources are independent of each other.

###  4.1 A two-block model

Our first example is a two-block land model similar to the one presented in
Grayver et al. (2013). The model has a background resistivity of 5
$\Omega\cdot$m, including two anomalies: a low resistivity inclusion of 1
$\Omega\cdot$m and a high resistivity inclusion of 100 $\Omega\cdot$m, as
shown in Figure 3a. The physical dimension of the model spans over a 3D domain
$X=X_{1}\times X_{2}\times X_{3}$, where $X_{1}=[-2000,2000]$ m,
$X_{2}=[-2000,2000]$ m and $X_{3}=[0,1500]$ m. There are 16 receivers and 256
transmitters with equal distance on the surface of this model. Figure 4 gives
the survey layout sheet of the acquisition geometry. The model has been
discretized with the regular grid using $\Delta x=\Delta y=60$ m and $\Delta
z=25$ m. The observed data are generated at two frequencies: 0.25 Hz and 1 Hz.

Figure 3: (a) The true land model (5 $\Omega\cdot$m background with two
anomalous inclusions of 1 $\Omega\cdot$m (left) and 100 $\Omega\cdot$m (right)
); (b) The homogeneous initial model of 5 $\Omega\cdot$m. The models are
clipped at 25 $\Omega\cdot$m for display purpose. Figure 4: Survey layout
sheet for the two-block model. In total, 16 transmitters (marked by dots) and
256 receivers (labeled with triangles) are deployed.

This 3D CSEM inversion runs for 30 iterations. In order to check how the
synthetic data fits the observations, we plot both the amplitude and the phase
of the synthetic data at all receivers together with the observed data using
the initial model and the retrieved resistivity model after inversion. In
Figure 5, the synthetic data created from the homogeneous initial model cannot
match the observed data very well. In Figure 6, the synthetic data generated
from inverted resistivity volume indeed are better aligned with the observed
data. Note that the horizontal axis in Figures 5 and 6 are receiver indices
rather than offset. This is due to the broadside survey configuration. To
better visualize the data matching connected to acquisition geometry, we
design a scatter plot with hot colors indicating the data error at each
receiver location. Figures 7 shows that at the beginning of the inversion, the
significant misfit (defined as
$\|\mathbf{W}(\mathbf{d}_{obs}-\mathbf{d}_{syn})\|$ at each receiver location)
is very large. Since we mute the data within 600 m offset, the significant
misfit surrounding transmitter Tx-10 are zeros. The plot in Figure 7 also
highlights the importance of azimuthal data has much larger significant misfit
based than inline directions. After 30 inversion iterations, the significant
misfit becomes much lower, see Figure 8.

Figure 5: Comparison between observed data and synthetic data from initial
model according to (a) amplitude and (b) phase for transmitter Tx-10. The
horizontal axis is the index of the receivers rather than offset, since we
consider broadside configuration. It is clear that both the amplitude and
phase does not match very well. Note also that the level of the data matching
has a strong correlation with the receiver indices. Figure 6: Comparison
between observed data and synthetic data from inverted model according to (a)
amplitude and (b) phase for transmitter Tx-10. The synthetic data are now
better aligned with the observed data after overlapping display. Figure 7:
Scatter plot of the significant misfit for transmitter Tx-10 at iteration 1.
Figure 8: Scatter plot of the significant misfit for transmitter Tx-10 at
iteration 30. The significant misfit becomes much lower after inversion.

The above data comparison gives us a confidence of our inversion scheme
working properly, according to a specific source (Tx-10). In Figure 9, we plot
the normalized data misfit which is a global measure of the inversion. We see
that after 30 iterations, it arrives at a relatively low normalized misfit.
According to Figure 10a, the frequency of large significant misfit at receiver
locations are distributed in a large range from 0 to 15. After the inversion,
it has been compressed a lot (most of them are clustered within 3), as shown
in Figure 10b. The final inverted model is capable to retrieve both the low
resistivity inclusion (Figure 11a) and the high resistivity inclusion (Figure
11b).

Figure 9: The convergence history in terms of the normalized misfit for the
inversion of the two-block model. Figure 10: The histogram of the significant
misfit at (a) iteration 1 and (b) iteration 30. Figure 11: Inversion result
for (a) low resistivity anomaly and (b) high resistivity anomaly.

###  4.2 A marine CSEM example

In practical CSEM survey for reservoir exploration, a number of towlines using
very low-frequency (0.1-10 Hz), high-energy electrical source will be deployed
with hundreds of receivers. The source-receiver offset extends up to 10 km is
very standard. Here, a 3D CSEM inversion is carried out for 30 iterations
based on a model of size $X=X_{1}\times X_{2}\times X_{3}$, where
$X_{1}=X_{2}=[-9000,9000]$ m and $X_{3}=[0,3500]$ m. Figure 12 illustrates the
survey layout sheet for the test: 10 towlines (5 in x direction and 5 in y
direction) are deployed in order to cover the region of interest. Each towline
has 81 source locations with a separation of 200 m. The receivers (25 in
total) are deployed at the crossings of these towlines. We consider the
reciprocity to switch the source and receiver, to achieve efficiency in
inversion.

As shown in Figure 13a, this synthetic model includes the seawater of 0.3
$\Omega\cdot$m. The resistivity of the background sediment is varying from 1.5
to 2.5 $\Omega\cdot$m along the depth. The most striking feature is a resistor
of 10 $\Omega\cdot$m at shallow part and a disk-shaped resistor of 100
$\Omega\cdot$m sitting at the depth between 2200 m and 2350 m, which mimics
the canonical reservoir. The starting model for the inversion in Figure 13b
takes a homogeneous value 1.5 $\Omega\cdot$m to mimic the situation that no
accurate a priori information is known for inversion. A seafloor bathymetry
with varying depth has been embedded in the model. The model is densely
gridded around bathymetry and above, with the grid spacing $\Delta z=25$ m,
see Figure 14. The grid has been stretched with a constant growing factor at a
certain distance below the seabed.

Three frequencies (0.25 Hz, 0.75 Hz and 2.25 Hz) have been used to perform
this inversion. Figure 15 shows that this CSEM inversion was converging well.
The comparison between Figure 16 and Figure 17 shows that the significant
misfit has been largely reduced for the receiver at the location (0, 0) m.
Note that in Figure 17c, the data misfit for the third frequency at the near
offset is still large, indicating the presence of shallow resistor, which has
not been well recovered in the inverted model (cf. Figure 19). To recover also
the shallow resistor, we should use more near offset data rather than drop off
them since 1000 m. A frequency dependent weighting strategy can be applied to
boost the importance of high frequencies, as they have shallow penetrating
depth according to skin depth. Since the scatter plot of the signficant misfit
does not reflect the global misfit, the significant misfit at each receiver
location has been shown in Figure 18: at the beginning of the inversion, they
are far from 1; after the inversion, the significant misfit at all receiver
locations are close to 1, which is the target value after adding noise into
the observed data.

After 30 l-BFGS iterations, the inversion was successful to recover the
resevoir in both vertical resistivity ($\rho_{v}=1/\sigma_{h}$) and the
horizontal resistivity ($\rho_{h}=1/\sigma_{v}$), but the shallow resistor is
difficult to obtain. The imprint of the shallow resistor can be found just
above the resevoir in our reconstructed model: the horizontal location is
correct but the depth is obviously incorrect. It can be seen that the
recovered anomaly in $\rho_{v}$ is of 100-150 m shallower than the true
location. This highlights the low resolution of CSEM inversion and the
ambuiguity of the depth. The inverted vertical resistivity in Figure 19a is
much better resolved than the horizontal resistivity in Figure 19b. The value
of the retrieved $\rho_{v}$ is higher than the retrieved $\rho_{h}$, while
both are smaller than 100 $\Omega\cdot$m and thicker than the truth. This is
due to the effect of anomalous transverse resistance (ATR) (Mittet and Morten,
2012, equation 11): because the CSEM response is proportional to the product
of resistivity difference times the thickness, lower resistivity with larger
thickness can create similar amplitude response compared with higher
resistivity with thin depth distribution. There is also weak increase of the
resistivity in $\rho_{h}$ in Figure 19b, suggesting that the initial model of
homogeneous background lower than the true model, as can be seen by comparing
with Figure 19c.

Figure 12: Survey layout sheet for marine CSEM inversion. Figure 13: (a) The
true resistivity model (note that the resistivity of the sediment below the
seabed is varying from 1.5 to 2.5 $\Omega\cdot$m); (b) The initial resistivity
model (the resistivity of the sediment takes a homogeneous value 1.5, to mimic
the situation that no accurate a priori information is known). The models are
clipped at 5 $\Omega\cdot$m for display purpose. Figure 14: The vertical
section of the nonuniform grid to mesh the 3D marine CSEM resistivity Figure
15: The convergence history of the normalized misfit for marine CSEM
inversion. Figure 16: Scatter plot of the significant misfit for receiver at
(0, 0) m in iteration 1. Figure 17: Scatter plot of the significant misfit for
receiver at (0, 0) m in iteration 30. The significant misfit becomes much
lower after inversion. Figure 18: The significant misfit at each receiver
location before and after inversion Figure 19: The result of 3D marine CSEM
inversion: (a) vertical resistivity $\rho_{v}$, (b) horizontal resistivity
$\rho_{h}$ and (c) true resistivity model. The display is in logorithmic
scale. Note that in the inverted vertical resistivity, there is a small uphill
variation just above the deep resistor which has a good correspondence to the
shallow resistor in horizontal direction, where the highest resistivity has
been retrieved at the center.

##  5 Discussion

In our method, the modelling jobs for each source can be parallelized
independently with perfect scalability. Since it is based on time stepping,
the EM fields at the next time step will overwrite the memory unit at the
previous step. This makes the memory usage of the method very economically.
Our computation was performed on a Intel(R) Xeon(R) Gold 6258R CPU @ 2.70GHz,
possessing 56 CPU cores (each core has 2 threads) and 128 GB memory. For the
two-block model of 3D grid size 139*139*99 (after padding with PML boundaries
on each side), combination of forward and adjoint modelling for building the
gradient associated with one source only takes the memory of approximately 402
MB to compute all frequencies. (The memory requirement using frequency domain
finite difference direct solver for such a large model size could not even fit
128 GB memory.) This implies that we can simulateneously launch more than 100
modelling jobs (corresponding to 100 independent sources) to fully explore the
small memory footprint of fictitous wave domain modelling in performing large
scale 3D inversion. For 16 sources, we can use 4 threads to parallelize the
computation by OpenMP using shared memory. The modelling takes less time at
the beginning of the inversion, as the model is conductive and has less
inhomogeneties. At the later iterations, each gradient building takes more
computing time than before, because the resistive anomalies gradually come
into the model, the time step $\Delta t$ becomes smaller and the number of
timesteps becomes larger due to stability condition. The inversion with 30
iterations can be completed within 10 hours. On average, building one gradient
(equivalent to one forward plus one adjoint solve) takes roughly 20 minutes.

It should be noted that the time series generated by our method only allows
the matching of the field at given frequencies. The above showcases simply
demonstrate that the proposed method works well, making fictitious wave domain
an attractive modelling engine for 3D practical inversion. However, one has to
be aware that the choice of regularization parameters $\beta$, $\alpha_{x}$,
$\alpha_{y}$, $\alpha_{z}$ are of paramount importance to the final inversion
result. Lelièvre and Farquharson (2013) have made a thorough investigation and
concluded that regularization of noisy data usually requires a lot of apriori
knowledge about the geologic structure, such that the related parameters can
be specified properly. Our choices may be far from perfect, and better choices
should be possible to significantly improve the inversion result, while
avoiding the false positives in imaging. These are, however, irrelevant to the
core contribution of this work.

The derivation in the appendix A develops an alternative adjoint formulation
to avoid solving the regularized linear inverse problem. Equation (40) gives
an explicit expression to build CSEM inversion gradient directly by cross-
correlation of fictitious time domain forward and adjoint fields. However, 3D
implementation of such a strategy is not trivial as it is mandatory to
simultaneously access fictitious time domain forward and adjoint field at each
time step. This is challenging due to the opposite directions in time between
forward and adjoint fields. In doing so, one may think of storing all 3D
wavefield snapshots during thunsands of times and accessing them one by one in
reverse order, which is extremely inefficient due to slow memory acess speed
and heavy IO traffic. A feasible approach is to use the wavefield
reconstruction method by storing the boundary values (Yang et al., 2016b), or
resort to optimal checkpointing strategy (Symes, 2007; Yang et al., 2016a).
Wavefield reconstruction by boundary storage implies each gradient computation
involves three times of simulation, two for forward fields and one for adjoint
fields. Optimal checkpointing requires much higher recomputation ratio and
thus even more computationally expensive. In summary, these strategies demand
either significant memory consumption for storing boundaries of forward fields
at every time step, or sophisticated operations to achieve optimal
combinatorial recomputation of the forward fields, thus left as the future
work. The scheme used in this paper is a hybrid combination of time domain
modelling and frequency domain inversion. To build the gradient, only two
times of modelling are needed, one for forward fields and one for adjoint
fields. Consequently, the gradient computed by product between forward and
adjoint fields in frequency domain becomes the most efficient option among all
equivalent formulations.

##  6 Conclusion

A nonlinear CSEM inversion based on fictitious wave domain modelling using
timestepping simulation has been presented. To construct the gradient using
forward and adjoint fields, the adjoint source time functions are estimated by
solving a regularized linear inverse problem. A matrix-free implementation
using the basis functions is proposed such that the linear inverse problem is
iteratively solved only once and then the resulting basis functions can be
reused to build up all adjoint source time functions at every receiver
location. The method enables efficient 3D CSEM inversion which are
demonstrated by numerical examples. The method will become more attractive for
land CSEM which normally requires much more number of frequencies than marine
CSEM. The proposed method can equally be applied to magnetotelluric (MT)
inversion to decipher the subsurface structure in even larger scale using much
lower frequencies of the EM recordings.

## Computer Code Availability

Name of the code: `libEMMI`

Functionality: Library for 3D controlled-source electromagnetic modelling and
inversion

Contact: Pengliang Yang, Harbin Institute of Technology, China

E-mail<EMAIL_ADDRESS>
System requirements: Linux OS

Programming language: C, Fortran

Result visualization: Python3, Madagascar (Scons) and Gnuplot

Compilation requirement: gcc compiler, mpicc compiler and make

Software required: fftw3 (http://fftw.org/)

The code has been released publically via github repository
https://github.com/yangpl/libEMMI.

## Acknowledgements

The author acknowledges the support from Chinese Fundamental Research Funds
for the Central Universities (AUGA5710010121) and National Natural Science
Fundation of China (42274156). The author is indebted to Rune Mittet for the
inspiring discussions and the encouragement of this work. The alternative
adjoint formulation in the appendix A was due to René-Édourt Plessix.

##  Appendix A An alternative adjoint formulation

The above approach uses diffusive Maxwell equation in frequency domain as the
state equation, which holds no matter what kind of modelling engine is used as
the solver. The misfit in (8) can also be expanded as

$\phi_{d}(m)=\frac{1}{2}\sum_{s}\sum_{k}|\mathbf{W}_{1}(E^{obs}(\mathbf{x}_{r},\omega_{k};\mathbf{x}_{s})-RE(\mathbf{x},\omega_{k};\mathbf{x}_{s}))|^{2}+|\mathbf{W}_{2}(H^{obs}(\mathbf{x}_{r},\omega_{k};\mathbf{x}_{s})-RH(\mathbf{x},\omega_{k};\mathbf{x}_{s}))|^{2}.$

Based on equation (4) and the correspondence between wave and diffusive EM
fields, the Lagrangian functional translates a constrained minimization into
an unconstrained optimization

$\begin{split}\mathcal{L}=&\frac{1}{2}\sum_{s}\sum_{k}|\mathbf{W}_{1}(E^{obs}(\mathbf{x}_{r},\omega_{k};\mathbf{x}_{s})-RE(\mathbf{x},\omega_{k};\mathbf{x}_{s}))|^{2}+|\mathbf{W}_{2}(H^{obs}(\mathbf{x}_{r},\omega_{k};\mathbf{x}_{s})-RH(\mathbf{x},\omega_{k};\mathbf{x}_{s}))|^{2}\\\
+&\Re\sum_{s}\sum_{k}\int_{X}\mathrm{d}\mathbf{x}\overline{\lambda_{1}(\mathbf{x},\omega_{k};\mathbf{x}_{s})}\bigg{(}E(\mathbf{x},\omega_{k};\mathbf{x}_{s})-\int_{0}^{T_{\max}}dt^{\prime}E^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})e^{-\sqrt{\omega_{k}\omega_{0}}t^{\prime}}e^{\mathrm{i}\sqrt{\omega_{k}\omega_{0}}t^{\prime}}\bigg{)}\\\
+&\Re\sum_{s}\sum_{k}\int_{X}\mathrm{d}\mathbf{x}\overline{\lambda_{2}(\mathbf{x},\omega_{k};\mathbf{x}_{s})}\bigg{(}H(\mathbf{x},\omega_{k};\mathbf{x}_{s})-\sqrt{-\frac{2\omega_{0}}{\mathrm{i}\omega_{k}}}\int_{0}^{T_{\max}}dt^{\prime}H^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})e^{-\sqrt{\omega_{k}\omega_{0}}t^{\prime}}e^{\mathrm{i}\sqrt{\omega_{k}\omega_{0}}t^{\prime}}\bigg{)}\\\
+&\sum_{s}\int_{0}^{T_{\max}}dt^{\prime}\int_{X}\mathrm{d}\mathbf{x}\underline{H}^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})\cdot\bigg{(}\nabla\times
E^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})+\mu^{\prime}\partial_{t}H^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})-M^{\prime}(\mathbf{x}_{s},t^{\prime})\bigg{)}\\\
+&\sum_{s}\int_{0}^{T_{\max}}dt^{\prime}\int_{X}\mathrm{d}\mathbf{x}\underline{E}^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})\cdot\bigg{(}\nabla\times
H^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})-\epsilon\partial_{t}E^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})-J^{\prime}(\mathbf{x}_{s},t^{\prime})\bigg{)},\end{split}$
(34)

where $\lambda_{1}(\mathbf{x},\omega;\mathbf{x}_{s})$ and
$\lambda_{2}(\mathbf{x},\omega;\mathbf{x}_{s})$ are Lagrangian multipliers in
the frequency domain, while
$\underline{E}^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})$ and
$\underline{H}^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})$ are Lagrangian
multipliers in the fictitious time domain.

Setting the Lagrangian with respect to state variables to zero
($\partial\mathcal{L}/\partial\lambda_{1}=\partial\mathcal{L}/\partial\lambda_{2}=\partial\mathcal{L}/\partial\underline{E}^{\prime}=\partial\mathcal{L}/\underline{H}^{\prime}=0$)
gives exactly the state equations. The adjoint equations are obtained by
setting the Lagrangian with respect to state variables, i.e.,

$\displaystyle\frac{\partial\mathcal{L}}{\partial
E(\mathbf{x},\omega_{k};\mathbf{x}_{s})}=0$
$\displaystyle\Leftrightarrow\lambda_{1}(\mathbf{x},\omega_{k};\mathbf{x}_{s})=R^{\dagger}\mathbf{W}_{1}^{\dagger}\mathbf{W}_{1}(E^{obs}(\mathbf{x}_{r},\omega_{k})-RE(\mathbf{x},\omega_{k};\mathbf{x}_{s})),$
(35a) $\displaystyle\frac{\partial\mathcal{L}}{\partial
H(\mathbf{x},\omega_{k};\mathbf{x}_{s})}=0$
$\displaystyle\Leftrightarrow\lambda_{2}(\mathbf{x},\omega_{k};\mathbf{x}_{s})=R^{\dagger}\mathbf{W}_{2}^{\dagger}\mathbf{W}_{2}(H^{obs}(\mathbf{x}_{r},\omega_{k};\mathbf{x}_{s})-RH(\mathbf{x},\omega_{k};\mathbf{x}_{s})),$
(35b) $\displaystyle\frac{\partial\mathcal{L}}{\partial
E^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})}=0$
$\displaystyle\Leftrightarrow\nabla\times\underline{H}^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})+\epsilon\partial_{t}\underline{E}^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})=\Re\sum_{k}\overline{\lambda_{1}(\mathbf{x},\omega_{k};\mathbf{x}_{s})}e^{-\sqrt{\omega_{k}\omega_{0}}t^{\prime}}e^{\mathrm{i}\sqrt{\omega_{k}\omega_{0}}t^{\prime}},$
(35c) $\displaystyle\frac{\partial\mathcal{L}}{\partial
H^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})}=0$
$\displaystyle\Leftrightarrow-\mu\partial_{t}\underline{H}^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})+\nabla\times\underline{E}^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})=\Re\sum_{k}\sqrt{-\frac{2\omega_{0}}{\mathrm{i}\omega_{k}}}\overline{\lambda_{2}(\mathbf{x},\omega_{k};\mathbf{x}_{s})}e^{-\sqrt{\omega_{k}\omega_{0}}t^{\prime}}e^{\mathrm{i}\sqrt{\omega_{k}\omega_{0}}t^{\prime}}.$
(35d)

Equations (35a) and (35b) reveals that $\lambda_{1}$ and $\lambda_{2}$ are
nothing more than the electric and magnetic components of the weighted data
residual, while equations (35c) and (35d) form a new set of Maxwell system. It
is important to remark that in deriving the above equations, we have tacitly
applied integration by parts in time and space, assuming zero initial
condition of the forward EM fields

$E^{\prime}(\mathbf{x},t^{\prime}=0)=H^{\prime}(\mathbf{x},t^{\prime}=0)=0,\quad\mathbf{x}\in
X,$ (36)

zero final conditions of the adjoint fields

$\underline{E}^{\prime}(\mathbf{x},t^{\prime}=T_{\max})=\underline{H}^{\prime}(\mathbf{x},t^{\prime}=T_{\max})=0,\quad\mathbf{x}\in
X,$ (37)

and homogeneous boundary condition in space for both forward and adjoint
fields

$E^{\prime}(\mathbf{x},\cdot)=H^{\prime}(\mathbf{x},\cdot)=\underline{E}^{\prime}(\mathbf{x},\cdot)=\underline{H}^{\prime}(\mathbf{x},\cdot)=0,\quad\mathbf{x}\in\partial
X.$ (38)

At the saddle point, we again obtain

$\frac{\partial\phi_{d}}{\partial\epsilon}=\frac{\partial\mathcal{L}}{\partial\epsilon}=-\sum_{s}\int_{0}^{T_{\max}}\mathrm{d}t^{\prime}\underline{E}^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s})\cdot\partial_{t}E^{\prime}(\mathbf{x},t^{\prime};\mathbf{x}_{s}).$
(39)

Since $\sigma_{ij}=2\omega_{0}\epsilon_{ij}$, the application of the chain
rule yields

$\frac{\partial\phi_{d}}{\partial\sigma_{ij}}=\frac{\partial\phi_{d}}{\partial\epsilon_{ij}}\frac{\partial\epsilon_{ij}}{\partial\sigma_{ij}}=-\sum_{s}\frac{1}{2\omega_{0}}\int_{0}^{T_{\max}}\mathrm{d}t^{\prime}\underline{E}^{\prime}_{j}(\mathbf{x},t^{\prime};\mathbf{x}_{s})\cdot\partial_{t}E^{\prime}_{i}(\mathbf{x},t^{\prime};\mathbf{x}_{s}).$
(40)

With the gradient at hand, one can then construct the descent direction based
on nonlinear minimization methods to solve the inverse problem.

## References

  * Abubakar et al. (2008) Abubakar, A., Habashy, T., Druskin, V., Knizhnerman, L., Alumbaugh, D., 2008.  2.5 D forward and inverse modeling for interpreting low-frequency electromagnetic measurements.  Geophysics 73, F165–F177.
  * Alumbaugh and Newman (1997) Alumbaugh, D., Newman, G., 1997\.  Three-dimensional massively parallel electromagnetic inversion-II. Analysis of a crosswell electromagnetic experiment.  Geophysical Journal International 128, 355–363.
  * Björck (1996) Björck, Å., 1996.  Numerical methods for least squares problems.  SIAM, Society for Industrial and Applied Mathematics, Philadelphia.
  * Chang-Chun et al. (2015) Chang-Chun, Y., Xiu-Yan, R., Yun-He, L., Yan-Fu, Q., Chang-Kai, Q., Jing, C., 2015\.  Review on airborne electromagnetic inverse theory and applications.  Geophysics 80, W17–W31.
  * Chave and Cox (1982) Chave, A.D., Cox, C.S., 1982\.  Controlled Electromagnetic Sources for Measuring Electrical Conductivity Beneath the Oceans 1. Forward Problem and Model Study.  Journal of Geophysical Research 87, 5327–5338.
  * Claerbout and Fomel (2008) Claerbout, J.F., Fomel, S., 2008\.  Image estimation by example: geophysical soundings image construction: multidimensional autoregression.  Stanford University.
  * Commer and Newman (2008) Commer, M., Newman, G.A., 2008\.  New advances in three-dimensional controlled-source electromagnetic inversion.  Geophysical Journal International 172, 513–535.
  * Constable (2010) Constable, S., 2010.  Ten years of marine CSEM for hydrocarbon exploration.  Geophysics 75, 75A67–75A81.
  * Constable and Srnka (2007) Constable, S., Srnka, L.J., 2007\.  An introduction to marine controlled-source electromagnetic methods for hydrocarbon exploration.  Geophysics 72, WA3–WA12.
  * Constable et al. (1986) Constable, S.C., Cox, C.S., Chave, A.D., 1986.  Offshore electro-magnetic surveying techniques, in: 56th Annual International Meeting, SEG, Expanded Abstracts, Society of Exploration Geophysicists. pp. 81–82.
  * Constable et al. (1987) Constable, S.C., Parker, R.L., Constable, C.G., 1987.  Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data.  Geophysics 52, 289–300.
  * Eidesmo et al. (2002) Eidesmo, T., Ellingsrud, S., MacGregor, L., Constable, S., Sinha, M., Johansen, S., Kong, F., Westerdahl, H., 2002\.  Sea bed logging (SBL), a new method for remote and direct identification of hydrocarbon filled layers in deepwater areas.  First break 20.
  * Ellingsrud et al. (2002) Ellingsrud, S., Eidesmo, T., Johansen, S., Sinha, M., MacGregor, L., Constable, S., 2002\.  Remote sensing of hydrocarbon layers by seabed logging (SBL): Results from a cruise offshore Angola.  The Leading Edge 21, 972–982.
  * Grayver et al. (2013) Grayver, A.V., Streich, R., Ritter, O., 2013.  Three-dimensional parallel distributed inversion of CSEM data using a direct forward solver.  Geophysical Journal International 193, 1432–1446.
  * Grayver et al. (2014) Grayver, A.V., Streich, R., Ritter, O., 2014.  3D inversion and resolution analysis of land-based CSEM data from the Ketzin CO2 storage formation.  Geophysics 79, E101–E114.
  * Gribenko and Zhdanov (2007) Gribenko, A., Zhdanov, M., 2007\.  Rigorous 3D inversion of marine CSEM data based on the integral equation method.  Geophysics 72, WA73–WA84.
  * 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., Martin, R., 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., Liu, G., Morrison, H., 1989.  A new approach to modeling the electromagnetic response of conductive media.  Geophysics 54, 1180–1192.
  * Lelièvre and Farquharson (2013) Lelièvre, P.G., Farquharson, C.G., 2013\.  Gradient and smoothness regularization operators for geophysical inversion on unstructured meshes.  Geophysical Journal International 195, 330–341.  URL: https://doi.org/10.1093/gji/ggt255, doi:10.1093/gji/ggt255.
  * Li and Key (2007) Li, Y., Key, K., 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.
  * MacGregor et al. (2007) MacGregor, L., Barker, N., Overton, A., Moody, S., Bodecott, D., 2007.  Derisking exploration prospects using integrated seismic and electromagnetic data-A Falkland Islands case study.  The Leading Edge 26, 356–359.
  * MacGregor and Tomlinson (2014) MacGregor, L., Tomlinson, J., 2014\.  Marine controlled-source electromagnetic methods in the hydrocarbon industry: A tutorial on method and practice.  Interpretation 2, SH13–SH32.
  * 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 and Morten (2012) Mittet, R., Morten, J.P., 2012\.  Detection and imaging sensitivity of the marine CSEM method.  Geophysics 77, E411–E425.
  * Morten et al. (2009) Morten, J.P., Bjørke, A.K., Støren, T., 2009.  CSEM data uncertainty analysis for 3D inversion, in: SEG Technical Program Expanded Abstracts 2009, Society of Exploration Geophysicists. pp. 724–728.
  * Mulder (2006) Mulder, W., 2006.  A multigrid solver for 3D electromagnetic diffusion.  Geophysical prospecting 54, 633–649.
  * Newman and Alumbaugh (1995) Newman, G.A., Alumbaugh, D.L., 1995\.  Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences.  Geophysical Prospecting 43, 1021–1042.
  * Nocedal and Wright (2006) Nocedal, J., Wright, S.J., 2006\.  Numerical Optimization.  2nd ed., Springer.
  * Oristaglio and Hohmann (1984) Oristaglio, M.L., Hohmann, G.W., 1984\.  Diffusion of electromagnetic fields into a two-dimensional earth: A finite-difference approach.  Geophysics 49, 870–894.
  * Paige and Saunders (1982) Paige, C.C., Saunders, M.A., 1982\.  ALGORITHM 583 LSQR : Sparse linear equations and least squares problems.  ACM Transactions on Mathematical Software 8, 195–209.
  * Plessix and Mulder (2008) Plessix, R.E., Mulder, W.A., 2008\.  Resistivity imaging with controlled-source electromagnetic data: depth and data weighting.  Inverse Problems 24, 034012\.
  * Puzyrev et al. (2013) Puzyrev, V., Koldan, J., de la Puente, J., Houzeaux, G., Vázquez, M., Cela, J.M., 2013\.  A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling.  Geophysical Journal International 193, 678–693.
  * Rochlitz et al. (2019) Rochlitz, R., Skibbe, N., Günther, T., 2019.  custEM: Customizable finite-element simulation of complex controlled-source electromagnetic data.  Geophysics 84, F17–F33.
  * Saad (2003) Saad, Y., 2003.  Iterative Methods for Sparse Linear Systems.  SIAM, Philadelphia.
  * Schwarzbach and Haber (2013) Schwarzbach, C., Haber, E., 2013\.  Finite element based inversion for time-harmonic electromagnetic problems.  Geophysical Journal International 193, 615–634.
  * Shantsev et al. (2020) Shantsev, D.V., Nerland, E.A., Gelius, L.J., 2020.  Time-lapse CSEM: how important is survey repeatability?  Geophysical Journal International 223, 2133–2147.
  * da Silva et al. (2012) da Silva, N.V., Morgan, J.V., MacGregor, L., Warner, M., 2012\.  A finite element multifrontal method for 3D CSEM modeling in the frequency domain.  Geophysics 77, E101–E115.
  * Sirgue et al. (2008) Sirgue, L., Etgen, J.T., Albertin, U., 2008.  3D Frequency Domain Waveform Inversion using Time Domain Finite Difference Methods, in: Proceedings 70th EAGE, Conference and Exhibition, Roma, Italy, p. F022.
  * Sirgue et al. (2010) Sirgue, L., Etgen, J.T., Albertin, U., Brandsberg-Dahl, S., 2010\.  System and method for 3D frequency domain waveform inversion based on 3D time-domain forward modeling.  US Patent 7,725,266.
  * Smith (1996a) Smith, J.T., 1996a.  Conservative modeling of 3-D electromagnetic fields, Part i: Properties and error analysis.  Geophysics 61, 1308–1318.
  * Smith (1996b) Smith, J.T., 1996b.  Conservative modeling of 3-D electromagnetic fields, Part II: Biconjugate gradient solution and an accelerator.  Geophysics 61, 1319–1324.
  * Støren et al. (2008) Støren, T., Zach, J.J., Maaø, F.A., 2008.  Gradient calculations for 3D inversion of CSEM data using a fast finite-difference time-domain modelling code, in: 70th EAGE Conference and Exhibition incorporating SPE EUROPEC 2008.
  * 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.
  * Symes (2007) Symes, W.W., 2007.  Reverse time migration with optimal checkpointing.  Geophysics 72, SM213–SM221.  URL: http://link.aip.org/link/?GPY/72/SM213/1, doi:10.1190/1.2742686.
  * Taflove and Hagness (2005) Taflove, A., Hagness, S.C., 2005\.  Computational Electrodynamics: The Finite-Difference Time-Domain Method.  3rd ed., Artech House.
  * Wang and Hohmann (1993) Wang, T., Hohmann, G.W., 1993\.  A finite-difference, time-domain solution for three-dimensional electromagnetic modeling.  Geophysics 58, 797–809.
  * Ward and Hohmann (1988) Ward, S.H., Hohmann, G.W., 1988\.  Electromagnetic theory for geophysical applications, in: Electromagnetic Methods in Applied Geophysics: Volume 1, Theory. Society of Exploration Geophysicists, pp. 130–311.
  * Yang (2023) Yang, P., 2023.  libEMM: A fictious wave domain 3D CSEM modelling library bridging sequential and parallel GPU implementation.  Computer Physics Communications 288, 108745.  doi:https://doi.org/10.1016/j.cpc.2023.108745.
  * Yang et al. (2016a) Yang, P., Brossier, R., Métivier, L., Virieux, J., 2016a.  Wavefield reconstruction in attenuating media: A checkpointing-assisted reverse-forward simulation method.  Geophysics 81, R349–R362.  doi:10.1190/geo2016-0082.1.
  * Yang et al. (2016b) Yang, P., Brossier, R., Virieux, J., 2016b.  Wavefield reconstruction from significantly decimated boundaries.  Geophysics 80, T197–T209.  doi:10.1190/GEO2015-0711.1.
  * Yang and Mittet (2023) Yang, P., Mittet, R., 2023.  Controlled-source electromagnetics modelling using high order finite-difference time-domain method on a nonuniform grid.  Geophysics 88, E53–E67.  doi:10.1190/geo2022-0134.1.
  * Zaslavsky et al. (2013) Zaslavsky, M., Druskin, V., Abubakar, A., Habashy, T., Simoncini, V., 2013.  Large-scale gauss-newton inversion of transient csem data using the model order reduction framework.  Geophysics 78, E161–E171.
  * Zhdanov and Keller (1994) Zhdanov, M.S., Keller, G.V., 1994\.  The geoelectrical methods in geophysical exploration. volume 31.