PDE-based geophysical modelling using finite elements: examples from 3D resistivity and 2D magnetotellurics

Numerical methods are routinely employed to calculate the geophysical response of models with known physical properties. Particularly 2D and 3D models in general require numerical solvers, except for simple geometries which exhibit a certain symmetry such as spheres and cylinders. Commonly tailor-made numerical modeling software is used within each specific geophysical discipline and consequently there exists a large collection of proprietary and open-source geophysical modeling software. In this paper we present a general finite-element solver, called escript, which is capable of solving a wide range of different geophysical modeling problems if they can be written in terms of partial differential equations (PDEs) with suitable boundary and initial conditions. Specifically we employ escript for modeling problems in applied geophysics, here demonstrated with examples of 3D direct current (DC) resistivity and 2D magnetotelluric forward modeling. Note that escript is also capable of PDE-constrained inversion for geophysical applications, however not demonstrated here. This paper was motivated by the work on multi-geophysical coupling in a hydrogeophysical context involving DC Resistivity and magnetotelluric measurements, and is the first to show in mathematical detail how escript is used to solve PDE-based geophysical problems. The paper serves as an introduction into PDE-based modeling with escript for geophysical applications. Based on the descriptions here, other more complicated multiphysics problems can be relatively easily solved.

The escript package is an extension of Python, and can be imported just like any other package in Python. By using an object oriented programming style it provides a flexible and easy-to-use programming environment for numerical simulations, whilst allowing for fast solutions of large models by performing time-intensive calculations in C++ and C. The code has been efficiently parallelised using MPI (mpi-forum.org), OpenMP (openmp.org) and hybrid modes; scripts can run on desktop computers as well as highly parallel supercomputers without changes. PDEs are solved on structured or unstructured 2D or 3D meshes, which facilitates modelling of complex and geologically realistic structures if required. Escript's abstract interface allows geoscientists to focus on solving the actual problem without being experts in numerical modelling; escript is actively developed at the Centre for Geoscience Computing, University of Queensland and is released under the Open Software Licence v3.0. Up-to-date information, documentation and support are found on https://launchpad.net/escript-finley.

General-purpose finite element solvers have found wide use especially in engineering fields using well-known commercial software packages such as Comsol, Abaqus or Ansys among others; these packages have also been applied for geophysical modelling (e.g. Dyksterhuis et al 2005, Zhang 2010, Butler and Sinha 2012, Ikard et al 2012). Open-source projects which are used in the geophysical community include FEniCS, OpenFoam, Code_Aster and Moose, among others. FEniCS started out as a project to automate the solution of mathematical models and takes a modular approach for solving PDEs. In a geophysical context, examples of its use are in turbulent flow and pumping test modelling (e.g. Mortensen et al 2011, Bakhos et al 2014). More examples can be found in Logg et al (2012). OpenFoam offers a wide range of modules in the context of fluid flows involving chemical reactions (e.g. Berselli et al 2015). Code_Aster focuses on structural mechanics problems (e.g. Rohmer 2014), and was also used for coupled hydro-mechanical modelling (Seyedi et al 2009). Moose is a relatively new finite element framework developed by the Idaho National Laboratory. It started as a project to tackle nuclear physics problems but has found applications in geophysical areas as well (e.g. Guo et al 2013). Most of the mentioned packages are GUI driven, which offers ease-of-access, however repeatability is arguably better served via a programming interface such as FEniCS or escript.

The escript solver has been used in a geophysical context for joint inversion of magnetic and gravity data (Gross and Kemp 2013), for erosion in porous media (Muhlhaus et al 2014), for reactive transport modelling in porous media (Poulet et al 2012), for shear band formation (Finzi et al 2013), and the growth of lava domes (Gross et al 2007). In this paper we demonstrate the application of escript in applied geophysics to solve the 3D DC resistivity forward problem and the 2D magnetotelluric forward problem. Both methods evaluate the geoelectric subsurface resistivity structure however with different approaches. In DC resistivity, current is injected into the ground via electrodes, which sets up an electric potential as a function of the subsurface resistivity distribution. The MT method utilises the naturally occurring fluctuations of the geo-electromagnetic fields as an inductive source to image the subsurface resistivity.

Given a known 3D resistivity distribution, the DC resistivity forward modelling solves a Poisson-type PDE for the electric potential. This involves Dirac delta functions as source terms and mixed boundary conditions. One of the main difficulties in calculating the DC resistivity response (i.e. solving for the electric potential) arises from the source term containing Dirac delta functions, which leads to singularities in the solution. Several different numerical approaches have been followed to solve the resistivity equation. For example Ma (2002) employed a boundary element method, a special case of the integral equation method, for DC resistivity modelling of 3D inhomogeneous bodies buried in a layered earth. The boundary element method was also employed by Blome et al (2009) for singularity removal to calculate the primary potential. Most numerical schemes employ finite-difference (FDM) or finite-element (FEM) methods, which are suitable for modelling complex 3D geometries (Li and Spitzer 2002). Dey and Morrison (1979) developed the first 3D finite-difference algorithm based on self-adjoint difference equations, which were obtained by integrating the continuity equation over elemental volumes. They furthermore incorporated a mixed boundary condition based on the asymptotic behaviour of the potential field in a homogeneous medium, which had the advantage of reducing the spatial extent of the domain and hence computational costs.

The finite-difference scheme from Lowry et al (1989) mathematically removed the Dirac delta singularities from the modelling process and reintroduced them later. This approach is also known as the secondary potential approach. In this approach the sought-after total potential is split into a primary potential, containing the Dirac sources, and a non-singular secondary potential caused by resistivity contrasts. The primary potential can be solved via analytical solutions, applicable for flat topography geometries.

Pain et al (2002) developed parallelised 2D/3D finite-element code, where the source-terms at nodes are set to rescaled basis functions or are approximated by Gaussian functions if the sources are not coinciding with mesh-nodes. Li and Spitzer (2002) used the singularity removal technique for a finite-difference and finite-element scheme. They used a horizontally layered earth or a vertical contact for the primary potential. The results they obtained were very similar for both schemes. Rücker et al (2006) employed quadratic-shape finite elements together with the secondary potential approach using a homogeneous half-space. The resistivity was set to the resistivity values at the location of the source-electrodes. Incorporation of topography into the primary potential was numerically evaluated via mesh-refinement at the location of the source-electrodes.

Penz et al (2013) employed a finite-difference scheme together with a modified secondary potential approach. They also used an analytic solution for the primary potential for models with non-flat topography. To account for the topography they modified the boundary conditions for the computation of the secondary potential. Ren and Tang (2014) used this technique for an adaptive finite-element solver, which automatically refined the mesh around the current electrodes based on error-estimates.

We implemented the secondary field approach in escript using unstructured meshes. Both, the primary and secondary potentials are calculated numerically for flat and non-flat topography. In order to calculate the primary potential with adequate accuracy, the unstructured mesh was refined around the locations of the source electrodes.

The second example solves the 2D magnetotelluric (MT) forward problem. The electromagnetic source fields originate in high altitudes and are therefore commonly considered to be plane wave fields with sounding periods ranging from  ∼10⁻³ seconds to  ∼10⁵ seconds, where longer periods translate to an increasing depth of investigation. The MT forward problem corresponds to a diffusion equation and is usually solved in the frequency domain. The plane wave source is incorporated via boundary conditions. This differs markedly from the DC resistivity problem. 2D models are characterised by an electric strike direction along which the electrical resistivity is constant. This leads to a decomposition of the electromagnetic fields into two independent modes: electric field intensity E parallel to strike, denoted as transverse electric (TE) mode or as E-polarisation; and magnetic induction B parallel to strike, which is denoted transverse magnetic (TM) mode, or B-polarisation (Chave and Jones 2012, pp.37). Because the magnetotelluric problem is solved in the frequency domain, it also involves complex arithmetic. In escript, complex quantities have to be decomposed into real and imaginary parts.

One of the earliest finite-element programs for 2D MT is presented in Kisak and Silvester (1975), which describes the procedures for the total field parallel to strike for E- and B-polarisation. Wannamaker et al (1987) describe a 2D MT modelling algorithm based on a secondary field approach, similar to the resistivity procedure described earlier. The primary or background field is typically formulated as the response of a 1D homogeneous or layered half-space. Various different numerical implementations of the magnetotelluric method, such as the integral equation and finite-difference methods, are reviewed by Zhdanov et al (1997). Franke-Börner (2003) implemented a 2D finite-element program in Matlab using the total field procedure for E- and B-polarisation. The paper by Key and Weiss (2006) describes an adaptive finite-element method using unstructured grids. In the adaptive finite-element method, a solution is calculated iteratively using successively refined grids until converging to the desired tolerance. Rung-Arunwan and Siripunvaraporn (2010) describe a parallel direct solver algorithm to solve the 2D MT problem based on a domain-decomposition approach. In Key and Ovall (2011) a parallel goal-oriented adaptive finite-element method is presented for 2.5D controlled-source electromagnetics and 2D magnetotelluric.

Depending on whether the total field approach or the secondary field approach is implemented, different boundary conditions apply. The secondary field approach generally employs homogeneous Dirichlet conditions, i.e. the values for the fields at the boundaries are set to zero; this requires that the domain boundaries are far away from the anomalous region. The secondary field approach explicitly involves source terms in the PDE, which are usually based on 1D solutions and the resistivity difference between host and anomalous region (Chave and Jones 2012, p.305). Boundary conditions for the total field formulation can also be based on 1D layered earth solutions and are implemented as non-homogeneous Dirichlet conditions; alternatively homogeneous Neumann conditions can be implemented if the sides of the domain are vertical (see Chave and Jones 2012, p.307). We implemented the total field approach in escript for TE and TM mode employing inhomogeneous Dirichlet boundary conditions, based on 1D layered earth solutions.

A linear second order boundary-value-problem (BVP) for the unknown scalar function u is processed in escript in the form of a templated PDE with coefficients A, B, C, D, X and Y. The coefficients comprise the physical properties and are generally dependent on their locations and on time/frequency. In order to solve a PDE the user has to write a script in Python which defines the values for the PDE coefficients typically as constants or as location dependent functions within the PDE domain. It is also possible to define the PDE coefficients depending on the solution of another PDE. This way non-linear PDEs, time dependent PDEs and coupled PDEs can be solved. Using tensorial notation and the Einstein summation convention, the templated PDE for a single-valued solution u is of the following form (Gross et al 2007):

$$\begin{eqnarray}&&-{{\left({{A}_{jl}}{{u}_{,l}}\right)}_{,\,j}}-{{\left({{B}_{j}}u\right)}_{,\,j}}+{{C}_{l}}{{u}_{,l}}+Du=-{{X}_{j,\,j}}+Y+\underset{\text{s}}{\sum}\,y_{\text{s}}^{\text{d}}\,{{\delta}_{{{\mathbf{x}}_{\text{s}}}}}\quad \text{in} ~\Omega \end{eqnarray} \tag{ 1 }$$

The notation u_{, l} denotes the partial derivative with respect to the lth spatial coordinate; ${{\mathbf{x}}_{\text{s}}}$ are locations of a point source (or sink) within the domain with production rates $y_{\text{s}}^{\text{d}}\,$; ${{\delta}_{{{\mathbf{x}}_{\text{s}}}}}$ denotes the Dirac delta function with center ${{\mathbf{x}}_{\text{s}}}$. Summation is performed over the spatial indices j and l. For example the first term in above template conforms to the following classical notation, with $\mathbf{A}\,\hat{=}\,\left({{A}_{jl}}\right)$ a second-rank tensor:

$$\begin{eqnarray}&&{{\left({{A}_{jl}}{{u}_{,l}}\right)}_{,j}}\,\Leftrightarrow \,\underset{j,l}{\sum}\,\frac{\partial}{\partial {{x}_{j}}}\left({{A}_{jl}}\frac{\partial}{\partial {{x}_{l}}}u\right)\,\Leftrightarrow \,\nabla \centerdot \left(\mathbf{A}\centerdot \nabla u\right)\end{eqnarray}$$

Boundary (and initial) conditions are required in order to obtain a solution in some domain $\Omega$. The boundary of $\Omega$ is denoted $\Gamma={{\Gamma}^{\text{D}}}{\cup}^{}{{\Gamma}^{\text{N}}}$, and ${{\Gamma}^{\text{D}}}{\cap}^{}{{\Gamma}^{\text{N}}}=\varnothing$, where ${{\Gamma}^{\text{D}}}$ specifies the parts of the boundary where essential (Dirichlet) boundary conditions apply and along which u has a prescribed value; ${{\Gamma}^{\text{N}}}$ denotes the boundary where natural (Neumann) boundary conditions apply and along which values of the derivative of u are specified. The following natural boundary conditions are considered on the domain boundary ${{\Gamma}^{\text{N}}}$:

$$\begin{eqnarray}&&{{n}_{j}}\left({{A}_{jl}}{{u}_{,l}}+{{B}_{j}}u\right)+du={{n}_{j}}{{X}_{j}}+y~\text{on}~{{\Gamma}^{\text{N}}}\end{eqnarray} \tag{ 2 }$$

where n_j denotes a component of the outer unit normal vector on ${{\Gamma}^{\text{N}}}$ and the coefficients A, B and X are the same as defined in the PDE, equation (1). The scalar coefficients d and y need to be defined on the entire boundary $\Gamma$ of the domain but are used on ${{\Gamma}^{\text{N}}}$ only. On ${{\Gamma}^{\text{D}}}$ essential (Dirichlet) boundary conditions prescribe the value of the solution of the form

$$\begin{eqnarray}&&u=r~\text{on}~{{\Gamma}^{\text{D}}}\end{eqnarray} \tag{ 3 }$$

where r is a function defining the values at these locations. Conventionally, the essential boundary conditions on any location of the domain boundary overwrite any natural boundary conditions that may have been set on this location via equation (2). Note that both types of boundary conditions can be combined and can be homogeneous or non-homogeneous. It is pointed out that escript also supports a template for systems of PDEs. This formulation is not presented here, but is used in the example for solving the magnetotelluric PDE; the escript users guide provides more details (Gross et al 2015).

In order to solve a PDE in escript for the unknown scalar dependent variable u, the governing equations have to be formulated in weak form. As opposed to the classical strong formulation of the BVP, which states the conditions at every point over a domain that a solution must satisfy, the weak form states the conditions that the solution must satisfy in a weighted average sense. The weak form of the BVP is defined as an inner product of u with a set of suitable weighting function, w, which vanish on ${{\Gamma}_{\text{D}}}$, i.e. w  =  0 on ${{\Gamma}_{\text{D}}}$. Carrying out the integration of the escript template over the domain $\Omega$, integrating by parts and applying Gauss' theorem, gives,

$$\begin{eqnarray}&&{{{\int}^{}}_{\Omega }}\text{d}\Omega \;{{w}_{,j}}\,\left({{A}_{jl}}{{u}_{,l}}+{{B}_{j}}u\right)\,+{{{\int}^{}}_{\Omega }}\text{d}\Omega \;w\,\left({{C}_{l}}{{u}_{,l}}+Du\right)\\ &&\quad \,-{{{\int}^{}}_{\Gamma}}\text{d} \Gamma\;w\,{{n}_{j}}\left({{A}_{jl}}{{u}_{,l}}+{{B}_{j}}u\right)={{{\int}^{}}_{\Omega }}\text{d}\Omega \;\left({{w}_{,j}}\,{{X}_{j}}+w\,Y\right)\\ &&\quad +\underset{\text{s}}{\sum}\,y_{\text{s}}^{\text{dirac}}w\left({{\mathbf{x}}_{\text{s}}}\right)\,-{{{\int}^{}}_{\Gamma}}\text{d} \Gamma\;w\,{{n}_{j}}\,{{X}_{j}}\end{eqnarray} \tag{ 4 }$$

Applying the natural boundary conditions from equation (2) and using the fact that w  =  0 on ${{\Gamma}_{\text{D}}}$ gives the weak form of the BVP, subject to essential constraints, equation (3):

$$\begin{eqnarray}&&{{{\int}^{}}_{\Omega }}\text{d}\Omega \;\,{{w}_{,j}}\left({{A}_{jl}}{{u}_{,l}}+{{B}_{j}}u\right)\,+{{{\int}^{}}_{\Omega }}\text{d}\Omega \;w\,\left({{C}_{l}}{{u}_{,l}}+Du\right)\,+{{{\int}^{}}_{\Gamma}}\text{d} \Gamma\;w\,\text{d}u\\ &&\quad ={{{\int}^{}}_{\Omega }}\text{d}\Omega \;\left({{w}_{,j}}\,{{X}_{j}}+w\,Y\right)+\underset{\text{s}}{\sum}\,y_{\text{s}}^{\text{dirac}}w\left({{\mathbf{x}}_{\text{s}}}\right)+{{{\int}^{}}_{\Gamma}}\text{d} \Gamma\;w\,y\end{eqnarray} \tag{ 5 }$$

The weak form PDE, equation (5), is numerically solved using the Ritz–Galerkin method, which converts the continuous weak form PDE into a discrete formulation. The solution u is approximated by a function ${{\tilde{u}}^{h}}$ which is a linear combination in some basis functions $\phi _{n}^{h}$ ($n=1,\ldots,N$). The construction of the basis functions is based on a subdivision of the domain into smaller pieces of maximum size h, so-called elements, sharing vertices $\mathbf{x}_{n}^{h}$ ($n=1,\ldots,N$), called nodes. The structure describing the location and the connectivity of the nodes to form the elements is called a mesh. For the sake of simplicity it is assumed that all locations ${{\mathbf{x}}_{\text{s}}}$ of point sources are also mesh nodes. The basis function $\phi _{n}^{h}$ corresponding to node $\mathbf{x}_{n}^{h}$ is constructed as a linear polynomial for each element, which is continuous across element boundaries and has value zero at all other nodes, i.e. $\phi _{m}^{h}\left(\mathbf{x}_{n}^{h}\right)={{\delta}_{mn}}$ for $m,n=1,\cdots N$, with Kronecker delta ${{\delta}_{mn}}$, see for instance (Zienkiewicz et al 2013) for details. With this construction the solution approximation is set to

$$\begin{eqnarray}&&{{\tilde{u}}^{h}}=\underset{m=1}{\overset{N}{\sum}}\,u_{m}^{h}\phi _{m}^{h}\end{eqnarray} \tag{ 6 }$$

where by construction $\tilde{u}_{m}^{h}={{u}^{h}}\left(\mathbf{x}_{m}^{h}\right)$ is the value of ${{\tilde{u}}^{h}}$ at node $\mathbf{x}_{m}^{h}$. Substituting equation (6) into the weak form, equation (5), and choosing the basis functions $\phi _{n}^{h}$ for the function w leads to a system of linear equations

$$\begin{eqnarray}&&{{\mathbf{K}}^{h}}\;{{\mathbf{u}}^{h}}={{\mathbf{F}}^{h}}\end{eqnarray} \tag{ 7 }$$

where ${{\mathbf{u}}^{h}}={{\left(u_{m}^{h}\right)}_{m=1\ldots N}}$ is the vector of basis function coefficients, ${{\mathbf{K}}^{h}}={{\left(K_{nm}^{h}\right)}_{m,n=1\ldots N}}$ is the stiffness matrix and ${{\mathbf{F}}^{h}}={{\left(F_{n}^{h}\right)}_{n=1\ldots N}}$ is the right-hand-side vector. For a node $\mathbf{x}_{n}^{h}$ which is not located within ${{\Gamma}_{\text{D}}}$ one obtains for the mth element ${{\Omega }_{m}}$:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} K_{nm}^{h}={{{\int}^{}}_{{{\Omega }_{m}}}}\text{d}\Omega ~{{A}_{jl}}\,\phi _{n,j}^{h}\,\phi _{m,l}^{h}\,+{{{\int}^{}}_{{{\Omega }_{m}}}}\text{d}\Omega ~{{B}_{j}}\,\phi _{n,j}^{h}\,\phi _{m}^{h}\; \\ \qquad \quad +{{{\int}^{}}_{{{\Omega }_{m}}}}\text{d}\Omega ~{{C}_{l}}\,\phi _{n}^{h}\,\phi _{m,l}^{h}\,+{{{\int}^{}}_{{{\Omega }_{m}}}}\text{d}\Omega ~D\,{{\phi}_{n}}\,{{\phi}_{m}}\,+{{{\int}^{}}_{{{\Gamma}_{m}}}}\text{d} \Gamma\;\text{d}\phi _{n}^{h}\,{{\phi}_{m}} \end{array}\end{eqnarray} \tag{ 8 }$$

for the entries in the stiffness matrix and

$$\begin{eqnarray}&&F_{n}^{h}={{{\int}^{}}_{{{\Omega }_{m}}}}\text{d}\Omega ~{{X}_{j}}\,\phi _{n,j}^{h}+{{{\int}^{}}_{{{\Omega }_{m}}}}\text{d}\Omega ~Y\,\phi _{n}^{h}+y_{n}^{\text{dirac}}+{{{\int}^{}}_{{{\Gamma}_{m}}}}\text{d} \Gamma\;y\,\phi _{n}^{h}\end{eqnarray} \tag{ 9 }$$

for the right-hand-side vector entries. To consider the essential boundary conditions, equation (3), one sets

$$\begin{eqnarray}&&K_{nm}^{h}={{\delta}_{nm}}\quad \text{and}\quad F_{n}^{h}=r\left(\mathbf{x}_{n}^{h}\right)\end{eqnarray} \tag{ 10 }$$

for those rows n where $\mathbf{x}_{n}^{h}$ are located on ${{\Gamma}_{\text{D}}}$. By construction of the basis function the stiffness matrix is sparse, i.e. the number of non-zero entries in the matrix is very small compared to the total size of the stiffness matrix. For these type of matrices iterative methods are most suitable (Weiss 1996) but escript also provides the option to use direct solvers to solve the linear system, equation (7), e.g. with umfpack (Davis 2004). In the discretisation approach discussed here it is assumed that values for B and C are sufficiently small so no additional stabilisation of the finite element matrix is required. For these cases escript provides a special solver using algebraic flux control, see Gross et al (2013).

For execution on a parallel compute cluster the rows of the stiffness matrix ${{\mathbf{K}}^{h}}$ are distributed across the processors. The distribution of the rows correspond to a distribution of mesh nodes and the attached elements can be seen as a subdivision of the domain into smaller subdomains that are assigned to processors for parallel calculation. Efficiency of this domain decomposition approach is controlled by the costs for exchanging nodal FEM data across sub-domain boundaries and an equally distributed computational load across processors (Saad 1996). For structured meshes this can easily be implemented but for unstructured meshes a compromise of these two objectives can be achieved only. In escript the ParMetis library is used for mesh decomposition (Schloegel et al 2002). It is pointed out that scripts for escript do not refer to the actual mesh or its distribution across a set of processors. As a consequence scripts can be run with different mesh types, e.g. structured or unstructured, and across different processor configurations without modifications to the script.

3.1. DC resistivity in 3D

The electric potential over a given three-dimensional domain $\Omega$ is obtained as the solution of a poisson-type PDE, which is derived in weak form and numerically solved using escript. Following (Pridmore et al 1981), the partial differental equation derives from the steady-state equation of continuity of current density, $\nabla \centerdot {{\mathbf{J}}^{c}}=0$, where ${{\mathbf{J}}^{c}}$ is the conduction current density in the subsurface. Ohm's law can be substituted for the current density, ${{\mathbf{J}}^{c}}=\sigma \mathbf{E}$, with σ the electrical conductivity and $\mathbf{E}$ the electric field intensity. Considering artificial sources of current density, ${{\mathbf{J}}^{\text{s}}}$, and writing the electric field as the gradient of a scalar potential, u, leads to the Poisson equation for the electric potential:

$$\begin{eqnarray}&&\nabla \centerdot (\sigma \nabla u)=\nabla \centerdot {{\mathbf{J}}^{\text{s}}}\end{eqnarray} \tag{ 11 }$$

For a single point source at location ${{\mathbf{x}}_{0}}$ and a measurement location $\mathbf{x}$ the non-stationary source current density, ${{\mathbf{J}}^{\text{s}}}$, follows the continuity equation:

$$\begin{eqnarray}&&\nabla \centerdot {{\mathbf{J}}^{\text{s}}}=-\frac{\partial {{q}_{\text{e}}}}{\partial t}\delta \left(\mathbf{x}-{{\mathbf{x}}_{0}}\right)=-I\delta \left(\mathbf{x}-{{\mathbf{x}}_{0}}\right)\end{eqnarray} \tag{ 12 }$$

where ${{q}_{\text{e}}}$ is the electrical charge density and I the source current strength. For N point sources at locations ${{\mathbf{x}}_{i}}$ with current strengths I_i, the potential field, $u\left(\mathbf{x}\right)$, is the solution of:

$$\begin{eqnarray}&&-\nabla \centerdot (\sigma \nabla u)=\underset{s=1}{\overset{{{N}_{\text{s}}}}{\sum}}\,{{I}_{i}}{{\delta}_{{{\mathbf{x}}_{\text{s}}}}}\equiv Q\end{eqnarray} \tag{ 13 }$$

Note that for the case of a pair of charging electrodes at locations $\mathbf{x}_{i}^{+}$ and $\mathbf{x}_{i}^{-}$ with current strengths $I_{i}^{+}=-I_{i}^{-}$ of opposite polarity, one sets for the source term Q

$$\begin{eqnarray}&&Q=\underset{s=1}{\overset{{{N}_{\text{s}}}}{\sum}}\,\left(I_{i}^{+}{{\delta}_{\mathbf{x}_{\text{s}}^{+}}}+I_{i}^{-}{{\delta}_{\mathbf{x}_{\text{s}}^{-}}}\right)\end{eqnarray} \tag{ 14 }$$

The Dirac delta functions in the source term lead to singularities and can be tackled by refining the mesh near the current source points and by splitting up the potential into a primary potential, ${{u}^{\text{p}}}$, and a secondary potential, ${{u}^{\text{s}}}$, for a background or primary conductivity ${{\sigma}^{\text{p}}}$ and an anomalous or secondary conductivity ${{\sigma}^{\text{s}}}=\sigma -{{\sigma}^{\text{p}}}$, respectively (Lowry et al 1989):

$$\begin{eqnarray}&&u={{u}^{\text{p}}}+{{u}^{\text{s}}}\end{eqnarray} \tag{ 15 }$$

The primary potential can be calculated analytically for flat-earth models as the response for a uniform half-space with conductivity ${{\sigma}^{\text{p}}}$. For topographic models ${{u}^{\text{p}}}$ is calculated numerically, applying a local mesh refinement at the source locations (see Rücker et al 2006). Here, we calculate the primary potential for flat and non-flat topographies numerically. The conductivity ${{\sigma}^{\text{p}}}$ is chosen as the conductivity near the source points (Zhao and Yedlin 1996). The primary potential ${{u}^{\text{p}}}$ satisfies the PDE equation (13):

$$\begin{eqnarray}&&-\nabla \centerdot \left({{\sigma}^{\text{p}}}\nabla {{u}^{\text{p}}}\right)=Q\end{eqnarray} \tag{ 16 }$$

The PDE for the secondary potential follows from substituting equation (15) into equation (13) and use equation (16), yielding:

$$\begin{eqnarray}&&\nabla \centerdot \left(\sigma \nabla {{u}^{\text{s}}}\right)=\nabla \centerdot \left(\left({{\sigma}^{\text{p}}}-\sigma \right)\,\nabla {{u}^{\text{p}}}\right)\end{eqnarray} \tag{ 17 }$$

To complete the problem, boundary conditions have to be applied for the partial differential equations of the primary and secondary potential. The boundary conditions are introduced below in combination with the weak form. The weak form is derived from the weighted residual method and is simply obtained by multiplying the PDEs, equations (16) and (17), with suitable weight functions, w, and integrating over the domain $\Omega$; subsequent use of integration by parts and Gauss' theorem provides the weak form PDE, see also figure 2. The weak form PDE for the primary potential evaluates to:

$$\begin{eqnarray}&&-{{{\int}^{}}_{\Gamma}}\text{d} \Gamma\;{{\sigma}^{\text{p}}}\,\omega \,\mathbf{n}\centerdot \nabla {{u}^{\text{p}}}+{{{\int}^{}}_{\Omega }}\text{d}\Omega \;{{\sigma}^{\text{p}}}\,\nabla \omega \centerdot \nabla {{u}^{\text{p}}}={{{\int}^{}}_{\Omega }}\text{d}\Omega \;\omega \,Q\end{eqnarray} \tag{ 18 }$$

and for the secondary potential one gets:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} -{{{\int}^{}}_{\Gamma}}\text{d} \Gamma\;\sigma \,\omega \,\mathbf{n}\centerdot \nabla {{u}^{\text{s}}}+{{{\int}^{}}_{\Omega }}\text{d}\Omega \;\sigma \,\nabla \omega \centerdot \nabla {{u}^{\text{s}}} \\ \quad ={{{\int}^{}}_{\Gamma}}\text{d} \Gamma\;\left({{\sigma}^{\text{p}}}-\sigma \right)\,\omega \,\mathbf{n}\centerdot \nabla {{u}^{\text{p}}}+{{{\int}^{}}_{\Omega }}\text{d}\Omega \;\left({{\sigma}^{\text{p}}}-\sigma \right)\,\nabla \omega \centerdot \nabla {{u}^{\text{p}}} \end{array}\end{eqnarray} \tag{ 19 }$$

Homogeneous natural and essential boundary conditions for the primary and secondary potentials can now be applied. At the air-earth boundary, ${{\Gamma}^{\text{N}}}$, the electric flux $\mathbf{n}\centerdot \nabla u$ is zero, since air is virtually non-conductive, i.e. Lowry et al (1989):

$$\begin{eqnarray}\left. \begin{array}{*{35}{r}} \sigma \,\mathbf{n}\centerdot \nabla {{u}^{\text{s}}} & =0 \\ {{\sigma}^{\text{p}}}\,\mathbf{n}\centerdot \nabla {{u}^{\text{p}}} & =0 \\ \left({{\sigma}^{\text{p}}}-\sigma \right)\,\mathbf{n}\centerdot \nabla {{u}^{\text{p}}} & =0 \end{array}\right\}~\text{on} ~{{\Gamma}^{\text{N}}}\end{eqnarray} \tag{ 20 }$$

For all other boundaries, ${{\Gamma}^{\text{D}}}$, it is assumed that the domain is adequately large for the potentials to have dissipated to zero, i.e.:

$$\begin{eqnarray}&&u=0~\text{on}~{{\Gamma}^{\text{D}}}\end{eqnarray} \tag{ 21 }$$

Therefore the integrals to be solved reduce to the integrals over the domain $\Omega$; for the primary and the secondary potential, in escript's formulation, one gets:

$$\begin{eqnarray}\begin{array}{*{35}{l}} \left(\text{primary}\right)\quad {{{\int}^{}}_{\Omega }}\text{d}\Omega \;{{\sigma}^{\text{p}}}\,{{\omega}_{,j}}\,u_{,j}^{\text{p}} & ={{{\int}^{}}_{\Omega }}\text{d}\Omega \;\omega \,Q \end{array}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \left(\text{secondary}\right)\quad {{{\int}^{}}_{\Omega }}\text{d}\Omega \;\sigma \,{{\omega}_{,j}}\,u_{,j}^{\text{s}} & ={{{\int}^{}}_{\Omega }}\text{d}\Omega \;\left({{\sigma}^{\text{p}}}-\sigma \right)\,{{\omega}_{,j}}\,u_{,j}^{\text{p}} \end{array}\end{eqnarray} \tag{ 23 }$$

The integrals are solved within the escript framework for all admissible (i.e. adequately smooth) weight functions ω. The next step is therefore to match the PDE coefficients, see equation (5), of the escript PDE template to the corresponding quantities as defined in the weak form formulations for the primary and secondary potentials. The boundary integrals in the template formulation evaluate to zero due to the homogeneous boundary conditions on ${{\Gamma}^{\text{N}}}$, i.e. the escript PDE coefficients d and y are identical zero. Via comparison, the non-zero PDE coefficients for the primary potential are:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} A_{jl}^{\text{p}}={{\sigma}^{\text{p}}}{{\delta}_{jl}},\quad y_{\text{s}}^{\text{d}}\,{{\delta}_{{{\mathbf{x}}_{\text{s}}}}}=Q\quad \left(\text{for}~s~\text{electrodes}\right) \end{array}\end{eqnarray} \tag{ 24 }$$

where $y_{\text{s}}^{\text{d}}\,{{\delta}_{{{\mathbf{x}}_{\text{s}}}}}$ are the Dirac point source terms defined in equation (1) and where ${{\delta}_{jl}}$ is Kronecker's delta function, so that the general tensor expression ${{A}_{jl}}=\sigma$ for all j  =  l. This can be seen by writing out the terms in equations (22), viz.,

$$\begin{eqnarray}&&{{\sigma}^{\text{p}}}\,{{\omega}_{,j}}\,u_{,j}^{\text{p}}={{\sigma}^{\text{p}}}{{\partial}_{x}}\,\omega \,{{\partial}_{x}}\,{{u}^{\text{p}}}+{{\sigma}^{\text{p}}}{{\partial}_{y}}\,\omega \,{{\partial}_{y}}\,{{u}^{\text{p}}}+{{\sigma}^{\text{p}}}{{\partial}_{z}}\,\omega \,{{\partial}_{z}}\,{{u}^{\text{p}}}\\ &&\qquad \quad \\ &&\;={{w}_{,j}}\,{{\sigma}^{\text{p}}}\,{{\delta}_{jl}}\,u_{,l}^{\text{p}}\end{eqnarray} \tag{ 25 }$$

For the secondary potential one obtains:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} A_{jl}^{\text{s}}=\sigma {{\delta}_{jl}},\quad {{X}_{j}}=\left({{\sigma}^{\text{p}}}-\sigma \right)\,u_{,j}^{\text{p}} \end{array}\end{eqnarray} \tag{ 26 }$$

the source term for the secondary potential, X_j, involves the solution for the primary potential, thus a coupled PDE system is solved together with essential boundary conditions enforced as per equation (21). The escript formulation of the DC resistivity problem is thus:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\int}^{}}_{\Omega }}\text{d}\Omega \;\,{{w}_{,j}}A_{jl}^{\text{p}}\,u_{,l}^{\text{p}}\, & ={{{\int}^{}}_{\Omega }}\text{d}\Omega \;w\,y_{\text{s}}^{\text{d}}\,{{\delta}_{{{\mathbf{x}}_{\text{s}}}}} \end{array}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\int}^{}}_{\Omega }}\text{d}\Omega \;\,{{w}_{,j}}A_{jl}^{\text{s}}\,u_{,l}^{\text{s}}\, & ={{{\int}^{}}_{\Omega }}\text{d}\Omega \;\left({{\sigma}^{\text{p}}}-\sigma \right)\,{{\omega}_{,j}}\,u_{,j}^{\text{p}} \end{array}\end{eqnarray} \tag{ 28 }$$

The secondary potential approach was implemented in escript for a layered half-space, a half-sphere model and a prism model. A scripted example for the prism model can be found in appendix.

The layered half-space response is shown in figure 1 and depicts a Schlumberger sounding over a highly resistive basement overlain by a conductive middle layer and a resistive top. Model parameters are described in the figure's caption. Note that both, the primary and secondary potentials were computed numerically with escript, where the primary potential was calculated with the resistivity of the top layer. Deviations between the analytical and numerical response are noticeable near the 10 m mark, coinciding with the transition of resistivities from top to middle layer.

Zoom In Zoom Out Reset image size

The second example is the half-sphere model from Rücker et al (2006); model parameters and response for a pole–pole spread traversing the half-sphere are displayed in figure 2. See the caption for a description of the model parameters. The analytical response is calculated via the formulae for the electric potential of a sphere in a full space by Zhdanov and Keller (1994, p.216); the potential for a half-sphere buried in a half-space is then obtained by multiplying the full-space potential by two. The equations for the full space potential are repeated here:

$$\begin{eqnarray}U(r,\theta )=\frac{I{{\rho}_{\text{p}}}}{4\pi R}+\frac{I{{\rho}_{\text{p}}}}{4\pi}\left\{\begin{array}{*{35}{l}} \underset{n=1}{\overset{\infty}{\sum}}\,\,{{K}_{n}}\,\frac{{{a}^{2n+1}}}{{{(rd)}^{n+1}}}\,{{P}_{n}}(\cos \theta ), & r>a \\ \underset{n=1}{\overset{\infty}{\sum}}\,\,{{K}_{n}}\,\frac{1}{d}\,{{\left(\frac{r}{d}\right)}^{n}}\,{{P}_{n}}(\cos \theta ), & r\leqslant a \end{array}\right. \end{eqnarray}$$

where I is the current strength, ${{\rho}_{\text{p}}}$ is the resistivity of the host, P_n are Legendre polynomials with argument $\cos \theta$; a is the sphere radius and R is the distance from the source to the observation point, d the distance from the sphere centre to the source point and r the distance from the sphere centre to the observation point; θ is the angle between d and r. The reflection coefficient K_n contains the sphere resistivity, ${{\rho}_{\text{s}}}$, and is defined as

$$\begin{eqnarray}&&{{K}_{n}}=\frac{n\left({{\rho}_{\text{s}}}-{{\rho}_{1}}\right)}{(n+1){{\rho}_{\text{s}}}+n{{\rho}_{1}}}\end{eqnarray}$$

Zoom In Zoom Out Reset image size

The electrical potential, based on above formulae, was calculated for 25 terms and multiplied with 2 to get the half-sphere solution, see Zhdanov and Keller (1994, p.219). The half-sphere solution matches the analytical response better on the left half of the profile which contains the source electrode; small deviations are observable on the right half.

Because the primary field was numerically evaluated, the 3D resistivity computations required a large number of nodes to account for the singularities at the electrode locations. In order to calculate the response relatively accurate for the half-sphere model a total of $50\,642$ were employed, which is similar to the number of nodes ($49\,341$) employed by Rücker et al (2006) for the total field computation of the half-sphere model. Incorporating analytical solvers for the primary field can markedly reduce the required number of nodes in case of flat-earth topographies.

The third example compares numerically evaluated apparent resistivities of a buried cube in a half-space. The escript values are compared to the integral equation (IE) values by Pridmore et al (1981, figure 5). Figure 3 shows the finite element mesh used to evaluate the primary and secondary electrical potentials with escript. The mesh has a total of $50\,435$ nodes with the same dimensions as for the half-sphere model $\left(2~\text{km}\times 2~\text{km}\times 1~\text{km}\right)$; depth-to-top of the ${{\left(2~\text{m}\right)}^{3}}$ cube is 0.5 m. The apparent resistivities are shown in figure 4 and were calculated at 21 electrodes for a dipole–dipole spread from  −10 m to 10 m along a centre line over the cube with a step length of 1 m. Shown are the apparent resistivity curves for three different electrode-pair separations n. Overall, the apparent resistivity values of both schemes are close.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. 2D magnetotelluric

Whereas the DC resistivity technique is driven by galvanic ground-coupling, the magnetotelluric method is driven by inductive coupling of naturally occurring electromagnetic plane-wave fields with the subsurface resistivity structures. The behaviour of magnetotelluric fields is governed by the low-frequency approximation of Maxwell's equations in a source free medium, yielding the following diffusion equations (e.g. Chave and Jones 2012, p.303):

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \nabla \times \frac{1}{\sigma}\,\nabla \times \mathbf{H}+\text{i}\omega \mu \mathbf{H}=0 \end{array}\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \nabla \times \frac{1}{\mu}\,\nabla \times \mathbf{E}+\text{i}\omega \sigma \mathbf{E}=0 \end{array}\end{eqnarray} \tag{ 30 }$$

where $\mathbf{H}$ is the magnetic field, $\mathbf{E}$ the electric field, σ is the electrical conductivity and $\mu ={{\mu}_{0}}$ the magnetic permeability assumed to have the value of free-space, ${{\mu}_{0}}=4\pi {{10}^{-7}}$; ω is the angular frequency and i the imaginary unit.

For the 2D model in figure 5 there is no variation along-strike in the x-direction, implying that all derivatives with respect to x are zero. The plane-wave fields are mutually orthogonal and decouple into two independent modes with respect to the strike direction: the TE-mode (transverse electric) and the TM-mode (transverse magnetic), which are both characterised by a single field-component parallel to the strike direction, E_x and H_x, respectively. These modes are also known as E-polarisation and B-polarisation.

Derivation of the partial differential equations for TE- and TM-mode follows from the quasi-static diffusion equations (29) and (30), where $\nabla ={{\left(0,{{\partial}_{y}},{{\partial}_{z}}\right)}^{\intercal}}$ and the relevant parallel electric or magnetic field is denoted as $\mathbf{u}=\left({{u}_{x}},0,0\right)$. The curl–curl expressions in the diffusion equation for both 2D modes evaluate then to a single term:

$$\begin{eqnarray}&&{{\left(\nabla \times f\,\nabla \times \mathbf{u}\right)}_{x}}=-{{\partial}_{y\,}}f\,{{\partial}_{y}}{{u}_{x}}-{{\partial}_{z\,}}f\,{{\partial}_{z}}{{u}_{x}}\end{eqnarray} \tag{ 31 }$$

where $f=f(\,y,z)$ is the property function which is either $1/\mu$ or $1/\sigma$ for the TE or TM mode respectively. The PDE for 2D magnetotellurics therefore simplifies to a scalar partial differential equation of the form

$$\begin{eqnarray}&&-\nabla \centerdot \left(\,f\;\nabla u\right)+g\,u=0\end{eqnarray} \tag{ 32 }$$

where the coefficients and the field are defined below:

$$\begin{eqnarray}\begin{array}{*{35}{l}} \text{TE:} & \quad f=\frac{1}{\mu}\quad g=\text{i}\omega \sigma \quad u={{E}_{x}} \end{array}\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \text{TM:} & \quad f=\frac{1}{\sigma}\quad g=\text{i}\omega \mu \quad u={{H}_{x}} \end{array}\end{eqnarray} \tag{ 34 }$$

Equation (32) defines the modelled field as the solution of the PDE. The orthogonal auxiliary fields are subsequently calculated from the PDE solution via Ampere's law (derivatives of E_x) and Faraday's law (derivatives of H_x):

$$\begin{eqnarray}\begin{array}{*{35}{l}} \text{TE:}\quad {{H}_{y}} & =-\frac{1}{\text{i}\omega \mu}\,{{\partial}_{z}}{{E}_{x}}\quad {{H}_{z}} & =\frac{1}{\text{i}\omega \mu}\,{{\partial}_{y}}{{E}_{x}} \end{array}\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \text{TM:}\quad {{E}_{y}} & =\frac{1}{\sigma}\,{{\partial}_{z}}{{H}_{x}}\quad {{E}_{z}} & =-\frac{1}{\sigma}\,{{\partial}_{y}}{{H}_{x}} \end{array}\end{eqnarray} \tag{ 36 }$$

Apparent resistivity and phase are defined for both modes via the magnetotelluric impedances, i.e. the ratio of the transverse electric and magnetic fields,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{Z}_{\text{TE}}}=-\frac{{{E}_{x}}}{{{H}_{y}}}\;\text{and}\;{{Z}_{\text{TM}}}=\,\frac{{{E}_{y}}}{{{H}_{x}}} \end{array}\end{eqnarray} \tag{ 37 }$$

and apparent resistivity and phase are given as:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\rho}_{m}}=\frac{1}{\omega \mu}|{{Z}_{m}}{{|}^{2}}\;\text{and}\;{{\phi}_{m}}={{\tan}^{-1}}\,\left(\frac{\Im \mathfrak{m}\,{{Z}_{m}}}{\Re \mathfrak{e}\,{{Z}_{m}}}\right) \end{array}\end{eqnarray} \tag{ 38 }$$

where the subscript m denotes either the TE- or TM-mode.

In 2D MT the modelled field quantities, E_x and H_x are both continuous across the entire domain; also all other components of the TE mode are continuous under the assumption that $\mu ={{\mu}_{0}}$. Only the TM-mode electrical field component E_y, equation (36), is discontinuous across the vertical contact. Telluric currents flow across structures in the TM mode and along structures in the TE mode. The TM mode is therefore characterised by anomalies of galvanic nature due to charge build-up, whereas TE mode anomalies are of inductive nature (Berdichevsky et al 1998). In consequence TM mode data is more sensitive to near-surface resistive structures and TE mode data more sensitive to deeper conductive structures. The modelling domain, $\Omega$, is set up differently for TE- and TM-mode. Due to the virtually infinite impedance at the air-earth interface, there are no induced electric field components in the air-layer, as such the air-layer can be ignored for TM modelling. In case of the magnetic field components of the TE mode, an air-layer must be included to account for magnetic field perturbations from longitudinal current gradients.

For all boundaries $\Gamma$ of the modelling domain, inhomogeneous boundary conditions of the Dirichlet type apply, i.e. essential boundary conditions are enforced (figure 6). On the left and right sides of the domain, the assumption is made that the conductivity distribution is a function of depth only and the electromagnetic fields are calculated analytically via 1D solutions on each side (see Wait 1953). The conditions on the top and bottom, respectively, are calculated via interpolation of the values at the left and right sides. For the TM-mode, the boundary condition is satisfied with a homogeneous Dirichlet condition at the air-earth interface with H_x  =  1 A m⁻¹. For the TE-mode, the electric field at the top of the air-layer is calculated via the free-space value of the electric field.

Zoom In Zoom Out Reset image size

The equations all contain complex terms. For escript these have to be decomposed into real and imaginary parts, $u={{u}^{r}}+\text{i}\,{{u}^{i}}$, where u is any of the complex field components and i is the imaginary unit. Substituting $u={{u}^{r}}+\text{i}\,{{u}^{i}}$ into the PDE, equation (32), yields,

$$\begin{eqnarray}&&\nabla \centerdot \left(\,f\,\nabla {{u}^{r}}+\text{i}\,f\,\nabla {{u}^{i}}\right)+g\,{{u}^{r}}+\text{i}\,g\,{{u}^{i}}=0\end{eqnarray} \tag{ 39 }$$

and substituting the terms for f and g provides the solution for TE- and TM-mode, which are written as two PDEs for real and imaginary part. For the TE-mode one obtains:

$$\begin{eqnarray}\left. \begin{array}{*{35}{l}} \nabla \centerdot \left({{\mu}^{-1}}\,\nabla E_{x}^{r}\right)+\omega \sigma E_{x}^{i} & =0 \\ \nabla \centerdot \left({{\mu}^{-1}}\,\nabla E_{x}^{i}\right)-\omega \sigma E_{x}^{r} & =0 \end{array}\;\right\}\end{eqnarray} \tag{ 40 }$$

and for the TM-mode:

$$\begin{eqnarray}\left. \begin{array}{*{35}{l}} \nabla \centerdot \left({{\sigma}^{-1}}\,\nabla H_{x}^{r}\right)+\omega \mu H_{x}^{i} & =0 \\ \nabla \centerdot \left({{\sigma}^{-1}}\,\nabla H_{x}^{i}\right)-\omega \mu H_{x}^{r} & =0 \end{array}\right\}\end{eqnarray} \tag{ 41 }$$

Thus in escript a coupled PDE must be solved for both modes which simultaneously solves for the real and imaginary part of the electric or magnetic field. To decompose apparent resistivity and phase into real and imaginary parts, the ratio of electric to magnetic fields are specified as (see equation (37)):

$$\begin{eqnarray}&&\frac{E}{H}=\frac{{{E}^{r}}+\text{i}{{E}^{i}}}{{{H}^{r}}+\text{i}{{H}^{i}}}=\frac{{{E}^{r}}{{H}^{r}}+{{E}^{i}}{{H}^{i}}}{{{\left({{H}^{r}}\right)}^{2}}+{{\left({{H}^{i}}\right)}^{2}}}+\text{i}\frac{{{E}^{i}}{{H}^{r}}-{{E}^{r}}{{H}^{i}}}{{{\left({{H}^{r}}\right)}^{2}}+{{\left({{H}^{i}}\right)}^{2}}}\end{eqnarray} \tag{ 42 }$$

where the subscripts indicating TE- or TM-mode have been omitted. Using the relationship $|z{{|}^{2}}={{x}^{2}}+{{y}^{2}}$ for a complex number z, the apparent resistivity can now be defined via the decomposed real and imaginary terms,

$$\begin{eqnarray}&&\rho =\frac{1}{\mu \omega}|\frac{E}{H}{{|}^{2}}=\frac{1}{\mu \omega}\left\{\frac{{{\left({{E}^{r}}\right)}^{2}}+{{\left({{E}^{i}}\right)}^{2}}}{{{\left({{H}^{r}}\right)}^{2}}+{{\left({{H}^{i}}\right)}^{2}}}\right\}\end{eqnarray} \tag{ 43 }$$

The phase evaluates to,

$$\begin{eqnarray}&&\phi ={{\tan}^{-1}}\left\{\frac{{{E}^{i}}{{H}^{r}}-{{E}^{r}}{{H}^{i}}}{{{E}^{r}}{{H}^{r}}+{{E}^{i}}{{H}^{i}}}\right\}\end{eqnarray} \tag{ 44 }$$

Note that the appropriate components have to be substituted for TE and TM mode respectively.

For implementation in escript, the weak form is derived and matched against the escript template. The weak form of equations (40) and (41) are combined into a vector-valued PDE system and evaluates to (see equation (18)):

$$\begin{eqnarray}&&{{{\int}^{}}_{\Gamma}}\text{d} \Gamma\;f\,\varphi \,\mathbf{n}\centerdot \nabla u+{{{\int}^{}}_{\Omega }}\text{d}\Omega \;f\,\nabla \varphi \centerdot \nabla u+{{{\int}^{}}_{\Omega }}\text{d}\Omega \;\varphi \,g\,u=0\end{eqnarray} \tag{ 45 }$$

where φ denotes the weighting function (to avoid confusion with the angular frequency ω), $f={{\mu}^{-1}}$ and $g=\omega \,\sigma$ for TE-mode, and $f={{\sigma}^{-1}}$ and $g=\omega \,\mu$ for TM-mode. The solution of the PDE, u, is a vector-valued quantity and yields ${{\left(E_{x}^{r},E_{x}^{i}\right)}^{\intercal}}$ and ${{\left(H_{x}^{r},H_{x}^{i}\right)}^{\intercal}}$, according to the TE- or TM-mode equations. In this context the functions f and g are tensors and the weighting function φ is vector-valued as well. At all boundaries the values of the fields are determined via analytical 1D solutions and the boundary integral term over $\Gamma$, which denotes the natural boundary condition, is zero.

Because equation (45) is a PDE system of two equations, it is cast to a more general weak form for PDE-systems in the templated escript formulation where the coefficient A is now defined as a rank-4 tensor and D a rank-2 tensor (see equation (5)),

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\int}^{}}_{\Omega }}\text{d}\Omega \;{{\varphi}_{i,\,j}}\,{{A}_{ijkl}}\,{{u}_{k,l}}-{{{\int}^{}}_{\Omega }}\text{d}\Omega \;{{\varphi}_{i}}\,{{D}_{ik}}{{u}_{k}}=0 \end{array}\end{eqnarray} \tag{ 46 }$$

only four tensor elements of PDE coefficients A_ijkl are non-zero and two elements of D_ik. These terms can be identified by writing out the equations (40) and (41); for example writing out equation (41) yields:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{y}}\left({{\sigma}^{-1}}\,{{\partial}_{y}}\,H_{x}^{r}\right)+{{\partial}_{z}}\left({{\sigma}^{-1}}\,{{\partial}_{z}}H_{x}^{r}\right)+\omega \mu H_{x}^{i} & =0 \end{array}\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{y}}\left({{\sigma}^{-1}}\,{{\partial}_{y\,}}H_{x}^{i}\right)+{{\partial}_{z}}\left({{\sigma}^{-1}}\,{{\partial}_{z}}H_{x}^{i}\right)-\omega \mu H_{x}^{r} & =0 \end{array}\end{eqnarray} \tag{ 48 }$$

By inspection, the non-zero values of A_ijkl and the two non-zero elements of D_ik evaluate to:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{A}_{ijkl}} & ={{\sigma}^{-1}}\quad \text{for}\;{{A}_{0000}},{{A}_{0101}},{{A}_{1010}},{{A}_{1111}} \end{array}\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{D}_{01}} & =\omega \mu \quad \text{and}\quad {{D}_{10}}=-\omega \mu \end{array}\end{eqnarray} \tag{ 50 }$$

The 2D MT problem was implemented in escript based on above description. In the following the escript implementation is tested via comparison with an analytical solution for a layered half-space and with the numerical results of two published forward models.

Figure 7 shows the response of a three-layered half-space, with escript computations carried out for TM-mode excitation. Model parameters are described in the figure's caption. The analytical solution is based on Wait's recursion formula for a layered earth (Wait 1953). Small deviations occur at short periods, whereas at all other periods the curves virtually overlap.

Zoom In Zoom Out Reset image size

The following two examples compare the escript apparent resistivity responses with published numerical solutions by the commemi project (Zhdanov et al 1997). The compared models are the commemi-1 and commemi-4 2D models, where the apparent resistivity values for comparison are taken from those by Weaver and Brewitt–Taylor; identified as '3 Canada' in tables B.8 and B.33 in the paper by Zhdanov et al (1997, p.194 and p.219) The phase values for the commemi-1 model are calculated from the field values in tables B.9 and B.10 of the same paper. Phase values for the commemi-4 model are compared with phase values as calculated from field values as reported by Franke-Börner (2003, tables 5.7 and 5.8).

The commemi-1 test model is shown in figure 8 together with the employed model parameters. The model simulates the response of a single conductive prismatic body in a more resistive uniform host. This model tests the ability to deal with highly varying values since the resistivity contrast is high, and the prism is large and near the surface, which therefore causes a large anomalous response. The escript responses for TE- and TM-mode are displayed in figures 9(a) and (b), respectively. The apparent resistivities as well as phases are compared with the published results and are in good agreement for both modes. Numerical values are given in table 1 for comparison. The phase values for comparison were calculated via equation (44) from the reported field values in Zhdanov et al (1997).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. commemi-1 escript apparent resistivity and phase at selected coordinates for te- and tm-mode compared with results by Weaver.

	x (km)	20.0	20.5	21.0	22.0	24.0	28.0	36.0
${{\rho}_{a}}$	escript (te)	8.08	14.06	49.91	95.88	104.10	100.31	100.12
	Weaver (te)	8.07	14.10	52.50	95.70	104.00	100.00	100.00
${{\rho}_{a}}$	escript (tm)	9.66	45.06	94.57	98.31	99.58	99.89	99.87
	Weaver (tm)	9.86	46.40	94.80	98.30	99.70	100.00	100.00
ϕ	escript (te)	76.00	71.49	65.90	53.60	46.11	44.99	45.03
	Weaver (te)	75.83	71.38	65.77	53.24	46.07	44.94	44.94
ϕ	escript (tm)	71.33	49.96	44.58	44.79	45.01	44.96	44.96
	Weaver (tm)	71.24	49.80	44.65	45.17	45.06	45.00	45.0

The next example employs the commemi-4 test model, which exhibits a more complex geometry as shown in figure 10; model parameters are given in the figure. The model consists of a total of five distinct regions, including an inclined boundary and a thin surficial layer, with different electrical resistivities varying between 2.5 $\Omega$m and 1000 $\Omega$m. The stratification on the left and right differs so that essential (Dirichlet) boundary conditions are calculated separately for both boundaries. The computed escript responses are compared in figures 11(a) and (b), and compared with Weaver's results in table 2 at selected coordinates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. commemi-4 escript apparent resistivity and phase at selected coordinates for te- and tm-mode as compared with the results by Weaver.

	x (km)	20.0	23.0	24.0	25.0	32.0	35.0
${{\rho}_{a}}$	escript (te)	12.70	12.00	8.84	6.87	6.70	6.27
	Weaver (te)	12.70	12.00	8.80	6.84	6.67	6.25
${{\rho}_{a}}$	escript (tm)	11.30	11.51	8.87	6.73	6.75	5.68
	Weaver (tm)	11.40	11.50	9.03	6.78	6.80	5.71
ϕ	escript (te)	45.92	53.54	56.84	59.71	62.50	54.69
	Weaver (te)	45.88	53.50	56.84	59.69	62.39	54.63
ϕ	escript (tm)	44.71	45.11	52.60	62.47	62.15	50.36
	Weaver (tm)	44.94	45.40	52.57	62.61	62.43	49.87

In order to calculate to the desired accuracy, a total of $58\,044$ nodes were used (te-mode). The number of nodes can be reduced when using the secondary field approach with modified boundary conditions, instead of the total field approach employed here. The apparent resistivity and phase for both modes in figure 11 reflect the layered structure at the left side and represents the graben-structure with a steep decline in apparent resistivity near the 24 km mark. The embedded graben-like structure tapers off as a layer to the right model boundary and causes strong current-channeling effects within the 2.5 $\Omega$m layer in both, TE- and TM-mode, as shown in Franke et al (2007). The apparent resistivities for both modes reflect the channelling effect above the 2.5 $\Omega$m layer, but more so for the TM-mode, where a significant decrease in apparent resistivity is noticeable near the 32 km mark. The escript values at selected coordinates are nearly identical to those of Weaver as shown in table 2 (note that the phase values for comparison are here based on the field values as reported in Franke-Börner (2003)).

General-purpose finite element solvers find increasing application in the geophysical disciplines as these offer a single interface to tackle different geophysical problems. These solvers are useful for data interpretation and for research, but can also be a useful tool in educational settings. We developed the parallel FEM solver escript as an extension of Python with a focus on solving boundary value problems in geophysics. A key concept of escript is the abstraction from the underlying spatial discretisation method (i.e. finite elements), providing a convenient programming environment. In principle this allows implementation of other discretisation methods such as boundary or spectral elements. Because the numerical implementations are independent from data structures, simulations are easily portable across desktop computers and scalable compute clusters without modifications to the program code.

In this paper we have demonstrated the use of escript at two different geophysical forward modelling examples, 3D DC resistivity and 2D magnetotellurics, with the intention to accurately compute the model responses and compare with different methods. For implementation in escript, the weak form was derived and matched against the escript template. Comparisons with analytical and other numerical solutions established confidence in the escript results. Escript is actively developed at the Centre for Geoscience Computing, University of Queensland and is released under the Open Software Licence v3.0. Up-to-date information, documentation and support are found on https://launchpad.net/escript-finley.