Massively parallel modeling and inversion of electrical resistivity tomography data using PFLOTRAN

Electrical resistivity tomography (ERT) is a broadly accepted geophysical method for subsurface investigations. Interpretation of field ERT data usually requires the application of computationally intensive forward modeling and inversion algorithms. For large-scale ERT data, the efficiency of these algorithms depends on the robustness, accuracy, and scalability on high-performance computing resources. In this regard, we present a robust and highly scalable implementation of forward modeling and inversion algorithms for ERT data. The implementation is publicly available and developed within the framework of PFLOTRAN, an open-source, state-of-the-art massively parallel subsurface flow and transport simulation code. The forward modeling is based on a finite-volume discretization of the governing differential equations, and the inversion uses a Gauss–Newton optimization scheme. To evaluate the accuracy of the forward modeling, two examples are first presented by considering layered (1D) and 3D earth conductivity models. The computed numerical results show good agreement with the analytical solutions for the layered earth model and results from a well-established code for the 3D model. Inversion of ERT data, simulated for a 3D model, is then performed to demonstrate the inversion capability by recovering the conductivity of the model. To demonstrate the parallel performance of PFLOTRAN's ERT process model and inversion capabilities, large-scale scalability tests are performed by using up to 131 072 processes on a leadership class supercomputer. These tests are performed for the two most computationally intensive steps of the ERT inversion: forward modeling and Jacobian computation. For the forward modeling, we consider models with up to 122 ×10⁶ degrees of freedom (DOFs) in the resulting system of linear equations and demonstrate that the code exhibits almost linear scalability on up to 10 000 DOFs per process. On the other hand, the code shows superlinear scalability for the Jacobian computation, mainly because all computations are fairly evenly distributed over each process with no parallel communication.

1 Introduction

Direct current electrical resistivity tomography (ERT) is one of the key geophysical methods for shallow-subsurface investigations, having applications in areas such as groundwater (Dahlin, 2001; Johnson et al., 2012; Meyerhoff et al., 2014; Park et al., 2016; Greggio et al., 2018; Alshehri and Abdelrahman, 2021), mineral exploration (Badmus and Olatinsu, 2009; Bery et al., 2012; Uhlemann et al., 2018; Martínez et al., 2019), environmental monitoring and remediation (Rosales et al., 2012; Rucker et al., 2013; Gabarrón et al., 2020; Rockhold et al., 2020; Kessouri et al., 2022), and engineering problems (Dahlin, 1996; Rizzo et al., 2004; Lysdahl et al., 2017). It can also be used for large-scale deep-subsurface investigations, e.g., for geothermal systems, active volcano imaging, and tectonic studies (Storz et al., 2000; Caputo et al., 2003; Johnson et al., 2010; Richards et al., 2010). Additional applications can also be found in detailed reviews by Slater (2007), Revil et al. (2012), Loke et al. (2013), Singha et al. (2022), and the references therein. In recent years, availability of multi-electrode and multichannel instrumentations permitted the acquisition of massive amounts of ERT data consisting of tens of thousands or even hundreds of thousands of observations. Efficient inversion and interpretation of such massive datasets requires fast, accurate, and highly scalable forward modeling and inversion algorithms to simulate and invert ERT data in arbitrary 3D conductivity structures. This paper presents an open-source implementation of such ERT modeling and inversion algorithms.

ERT data simulation is performed by solving the electrostatic Poisson equation. Numerical methods are needed to solve the governing equation in arbitrary 3D conductivity structures. The most used numerical methods for 3D ERT modeling are the finite-difference (FD) method (Dey and Morrison, 1979; Spitzer, 1995; Penz et al., 2013), finite-element (FE) method (Coggon, 1971; Li and Spitzer, 2002; Rücker et al., 2006; Blome et al., 2009; Johnson et al., 2010; Ren et al., 2018), and integral equation method (Lee, 1975; Schulz, 1985; Méndez-Delgado et al., 1999). Among these methods, the FD and FE methods are attractive choices for highly heterogeneous distributions of electrical conductivity in the subsurface. However, application of the FD method is mostly limited to simple model geometries due to its requirement of Cartesian grids (Rücker et al., 2006). By contrast, the FE method has proven to be effective in accounting for complex geometries, specifically by using unstructured meshes (Rücker et al., 2006; Blome et al., 2009; Johnson et al., 2010; Ren et al., 2018).

The finite-volume (FV) method can also be used effectively to account for complex geometries (Jahandari and Farquharson, 2014). The FV method is usually seen relatively close to the FD method as it inherits the simplicity similar to the FD method in its implementation. Unlike the FD method which discretizes the differential form of the governing equation, the FV method, however, directly discretizes its integral form. Nevertheless, the FV method has rarely been considered for ERT simulation problems. In our literature search, we came across only Cockett et al. (2015), which implements the method for the ERT modeling.

On the other hand, ERT data inversion is a nonlinear optimization problem to minimize a cost function that represents a measure of the difference between observed and simulated ERT data. For large-scale 3D ERT data, the inversion is often performed iteratively by linearizing the optimization problem at each iteration, calculating a model update, updating a conductivity model, and subsequently simulating ERT data for the updated model to examine the new cost function. The process continues until the cost function reduces to a predefined tolerance level. The model update is typically calculated using a gradient-based local optimization method, e.g., the steepest-descent, conjugate-gradient, or the Gauss–Newton method (Park and Van, 1991; Ellis and Oldenburg, 1994; Zhdanov and Keller, 1994; Zhang et al., 1995; Günther et al., 2006; Johnson et al., 2010). The Gauss–Newton method has usually been preferred due to its high convergence rate. Indeed, the fast convergence occurs at the cost of intensive computations because the method requires building the Jacobian (or sensitivity) and/or Hessian matrices at each inversion iteration (Nocedal and Wright, 2006).

For large-scale 3D ERT problems with tens of millions of degrees of freedom (DOFs), forward modeling and inversion are computationally demanding jobs and may need hours, if not days, of computing time if the underlying algorithms are not properly implemented to scale well on high-performance computing (HPC) resources or supercomputers. For such datasets, it is computationally infeasible to use desktop-based ERT data processing software without exploiting HPC resources. Recently, a few open-source 3D ERT modeling and inversion codes have been developed to handle moderately sized ERT problems composed of up to a few million DOFs, e.g., ResIPy (Blanchy et al., 2020), SimPEG (Cockett et al., 2015), and pyGIMLI (Rücker et al., 2017). However, these codes still lack full-scale parallel implementations to efficiently distribute computations over a supercomputer. On the other hand, another open-source code, E4D, is highly scalable on HPC resources (Johnson et al., 2010; Johnson and Wellman, 2015). It can be run on hundreds or even thousands of processes (or cores) on a supercomputer (Johnson et al., 2010). But despite this, the maximum number of processes that it can exploit is limited to the number of electrodes in a survey.

In this paper, we present an implementation of massively parallel 3D ERT modeling and inversion algorithms. The forward modeling is based on the FV method, and inversion employs the Gauss–Newton method. The algorithms are implemented within PFLOTRAN (https://www.pflotran.org, last access: 27 January 2023), an open-source, massively parallel code that leverages parallel data structures and solvers within PETSc (Portable, Extensible Toolkit for Scientific Computation) (Balay et al., 2021) to simulate subsurface flow and transport processes on supercomputers (Hammond et al., 2012, 2014). Although we are presenting only ERT modeling and inversion capabilities here, the main motivation for implementing these capabilities in PFLOTRAN is to perform coupled modeling and inversion of ERT, flow, and/or transport problems in the future.

The paper is structured as follows: we first outline a detailed background theory on the ERT forward modeling including the governing Poisson equation and FV method to solve it numerically. We then describe the ERT inversion by defining it as a nonlinear optimization problem and the Gauss–Newton method to find a minimizer of the problem along with details on computing the Jacobian matrix. Thereafter, we discuss the parallel implementations of the modeling and inversion algorithms. We subsequently benchmark numerical results computed using PFLOTRAN against analytical solutions for a layered (1D) earth model and numerical solutions for a 3D earth model. Next, a 3D inversion result is presented to illustrate the inversion capability. Scalability tests are then performed to demonstrate the parallel performance of the PFLOTRAN ERT process model before discussing the results and drawing our final concluding remarks.

2 Forward ERT modeling

Forward modeling is a way of simulating ERT data for any given 1D, 2D, or 3D electrical conductivity models of the earth. For 1D conductivity models, the electrical potentials generated by an induced point current source can be obtained by using analytical or semi-analytical methods (Das, 1995; Pervago et al., 2006). However, for multi-dimensional heterogeneous conductivity models with complex geometries, one must use numerical methods to compute the electrical potential. This section describes the numerical computations of the electrical potentials in a 3D medium and their superposition to produce simulated ERT data.

2.1 Poisson's equation

The ERT forward modeling is governed by the following electrostatic Poisson equation:

$$\begin{array}{}\text{(1)}& \mathbf{\nabla }\cdot \mathit{\sigma }\left(\mathit{r}\right)\mathbf{\nabla }\mathit{\varphi }\left(\mathit{r}\right)=-I\mathit{\delta }(\mathit{r}-{\mathit{r}}_{\mathrm{s}}),\end{array}$$

where ϕ is the electrical potential at position r for a given conductivity σ(r) due to a current I injected through a point located at r_s and δ is the Dirac delta function. For brevity, hereinafter, the dependencies on r will be omitted except where it is necessary to show them.

If side and bottom boundaries ∂Γ of the 3D computational domain Γ are located at sufficiently far from the current injection location r_s, the potential and the normal components of the current density $\mathit{\sigma }\frac{\partial \mathit{\varphi }}{\partial n}$ asymptotically approach zero. On the top or surface boundary ∂Γ_S, the normal current density $\mathit{\sigma }\frac{\partial \mathit{\varphi }}{\partial n}$ is zero as no current flows through the earth surface along the outward normal vector. Consequently, we can impose zero Dirichlet or Neumann boundary conditions at the side boundaries and zero Neumann boundary at the surface boundary:

$$\begin{array}{}\text{(2)}& \mathit{\varphi }{\mathrm{|}}_{\partial \mathrm{\Gamma }}=\mathrm{0}\mathrm{or}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial n}}{\mathrm{|}}_{\partial \mathrm{\Gamma }}=\mathrm{0},\mathrm{and}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial n}}{\mathrm{|}}_{\partial {\mathrm{\Gamma }}_{\mathrm{S}}}=\mathrm{0}.\end{array}$$

To simulate ϕ in arbitrary 3D conductivity structures, Eq. (1) needs to be solved using numerical methods subject to the boundary conditions in Eq. (2).

2.2 Finite-volume method

The FV method is implemented to compute ϕ by solving Eq. (1). The 3D computational domain is first discretized into a set of control volumes known as cells (Fig. 1). These cells can be of arbitrary shape and size, but they must be bounded by planar surfaces for our implementation. The bounding discrete surfaces are known as cell faces.

Figure 1 Two adjacent FV cells i₁ and i₂ having conductivity ${\mathit{\sigma }}_{i_{\mathrm{1}}}$ and ${\mathit{\sigma }}_{i_{\mathrm{2}}}$, respectively. The distances of the common face f from the center of cells i₁ and i₁, respectively, are $d_{i_{\mathrm{1}}}$ and $d_{i_{\mathrm{2}}}$.

Let us consider the governing Poisson equation (Eq. 1) and integrate it over the ith cell of the domain Γ. This gives

$$\begin{array}{}\text{(3)}& \underset{V_i}{\int }\mathbf{\nabla }\cdot \mathit{\sigma }\mathbf{\nabla }\mathit{\varphi }\mathrm{d}V=-\underset{V_i}{\int }I\mathit{\delta }(\mathit{r}-{\mathit{r}}_{\mathrm{s}})\mathrm{d}V,\end{array}$$

where V_i is the volume of the ith cell.

By applying the Gauss divergence theorem in Eq. (3) and using the translation property of the Dirac delta function, we get

$$\begin{array}{}\text{(4)}& \underset{S_i}{\int }\mathit{\sigma }\mathbf{\nabla }\mathit{\varphi }\cdot \widehat{\mathit{n}}\mathrm{d}S=-I{\mathrm{|}}_{\mathit{r}={\mathit{r}}_{\mathrm{s}}},\end{array}$$

where S_i is the surface area of the bounding faces of the ith cell, and dS is a differential area on the bounding surface with a unit surface normal $\widehat{\mathit{n}}$ pointing outward.

If the ith cell is bounded by N_f faces, Eq. (4) can be replaced by an equivalent discrete form

$$\begin{array}{}\text{(5)}& \sum _{f=\mathrm{1}}^{N_{\mathrm{f}},i}\underset{S_{\mathrm{f}}}{\int }\mathit{\sigma }\mathbf{\nabla }\mathit{\varphi }\cdot \widehat{\mathit{n}}\mathrm{d}S=\sum _{f=\mathrm{1}}^{N_{\mathrm{f}},i}{\mathit{\sigma }}_{\mathrm{f}}(\mathbf{\nabla }\mathit{\varphi }{)}_{\mathrm{f}}\cdot {\widehat{\mathit{n}}}_{\mathrm{f}}{S}_{\mathrm{f}}=-I{\mathrm{|}}_{\mathit{r}={\mathit{r}}_{\mathrm{s}}},\end{array}$$

where f represents faces and S_f is the area of the bounding face f.

Using a two-point flux approximation for each face term in Eq. (5) by considering two adjacent cells i₁ and i₂ shown in Fig. 1, we get

$$\begin{array}{}\text{(6)}& \sum _{f=\mathrm{1}}^{N_{\mathrm{f}},i}{\mathit{\sigma }}_{\mathrm{f}}(\mathbf{\nabla }\mathit{\varphi }{)}_{\mathrm{f}}{S}_{\mathrm{nf}}=-I{\mathrm{|}}_{\mathit{r}={\mathit{r}}_{\mathrm{s}}},\end{array}$$

where S_nf is the area of the common face f projected onto the plane normal to the vector connecting centers of cells i₁ and i₂.

The conductivity at the face, σ_f, can be obtained by averaging the conductivity of cells i₁ and i₂. We use the following harmonic distance-weighted scheme to calculate σ_f:

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\mathrm{f}}={\displaystyle \frac{{\mathit{\sigma }}_{i_{\mathrm{1}}}{\mathit{\sigma }}_{i_{\mathrm{2}}}(d_{i_{\mathrm{1}}}+d_{i_{\mathrm{2}}})}{{\mathit{\sigma }}_{i_{\mathrm{1}}}{d}_{i_{\mathrm{2}}}+{\mathit{\sigma }}_{i_{\mathrm{2}}}{d}_{i_{\mathrm{1}}}}},\end{array}$$

where ${\mathit{\sigma }}_{i_{\mathrm{1}}}$ and ${\mathit{\sigma }}_{i_{\mathrm{2}}}$ are the conductivities of cells i₁ and i₂ and $d_{i_{\mathrm{1}}}$ and $d_{i_{\mathrm{2}}}$ are the distances of the projected face from the center of cells i₁ and i₁, respectively (Fig. 1).

The value of the flux at the face f, i.e., (∇ϕ)_f, is obtained by using the approximation

$$\begin{array}{}\text{(8)}& (\mathbf{\nabla }\mathit{\varphi }{)}_{\mathrm{f}}={\displaystyle \frac{{\mathit{\varphi }}_{i_{\mathrm{2}}}-{\mathit{\varphi }}_{i_{\mathrm{1}}}}{d_{i_{\mathrm{1}}}+d_{i_{\mathrm{2}}}}}.\end{array}$$

Substituting Eqs. (7) and (8) in Eq. (6) and assuming that the computational domain Γ is subdivided into N_m FV cells, i.e., $\mathrm{\Gamma }={\cup }_{i=\mathrm{1}}^{N_{\mathrm{m}}}{\mathrm{\Gamma }}^i$, where Γⁱ⊂ℝ³ is the ith cell domain, we have

$$\begin{array}{}\text{(9)}& \sum _{i=\mathrm{1}}^{N_{\mathrm{m}}}{A}_{i,k}{\mathit{\varphi }}_k=s_i,\end{array}$$

where A_i,k is the local coefficient matrix and s_i is the source vector for the ith cell.

The assembly of the local coefficient matrix into a global system, upon carrying out the summation in Eq. (9), results into the system of linear equations

$$\begin{array}{}\text{(10)}& \mathbf{A}\mathit{\varphi }=\mathit{s},\end{array}$$

where $\mathbf{A}\in {\mathbb{R}}^{N_{\mathrm{m}}\times N_{\mathrm{m}}}$ is the symmetric system matrix resulting from the FV discretization of Eq. (1), $\mathit{\varphi }\in {\mathbb{R}}^{N_{\mathrm{m}}}$ is a vector containing the unknown electric potential for all cells, and $\mathit{s}\in {\mathbb{R}}^{N_{\mathrm{m}}}$ is the source vector resulting from the right-hand side of Eq. (1). To solve this linear system, we utilize various options of iterative solvers available from PETSc. Each electrode in an ERT survey results in a different vector s. Thus, for N_e electrodes, we need to solve Eq. (10) for N_e times with a different right-hand side. The solutions provide electrical potential at each cell center of the discrete computational domain and for each current injection location. Finally, the potential distribution for any combination of source/sink electrodes can be obtained by subtracting the potential associated with the sink electrodes from the source electrode potential (Johnson et al., 2010).

3 Inversion

ERT inversion is a procedure of estimating electrical conductivity of the earth from a set of observed ERT data and/or some imposed constraints. The ERT inverse problem is a nonlinear optimization problem. For large-scale 3D ERT data, the inversion is often performed iteratively by linearizing the problem at each iteration. This section details our approach of performing inversion of ERT data.

3.1 Optimization problem

We pose the inverse problem as a nonlinear minimization problem of finding a conductivity model

$$\begin{array}{}\text{(11)}& {\mathit{m}}_{\mathrm{inv}}=\arg \underset{\mathit{m}\in \mathbf{M}}{min}\mathit{\psi }\left(\mathit{m}\right),\end{array}$$

where $\mathit{m}=(m_{\mathrm{1}},m_{\mathrm{2}},\mathrm{\dots },m_{N_{\mathrm{m}}}{)}^{\mathrm{T}}$ is an unknown model parameter vector in the model set M. Since the conductivity within subsurface can vary over several order of magnitudes, the inversion operations are performed for logarithmic conductivity, i.e., we assume m=ln σ. This also enforces positivity of conductivity values. The quantity ψ is the cost function, defined as

$$\begin{array}{}\text{(12)}& \mathit{\psi }\left(\mathit{m}\right)={\mathit{\psi }}_{\mathrm{d}}\left(\mathit{m}\right)+\mathit{\beta }{\mathit{\psi }}_{\mathrm{m}}\left(\mathit{m}\right),\end{array}$$

where ψ_d is the data cost function that represents a measure of the misfit between the observed d^obs and synthetic d^syn ERT data, ψ_m is the regularization cost function representing a measure of the variation in the model parameters, and β is a regularization parameter.

A least-squares method is used to minimize ψ (Nocedal and Wright, 2006). Thus, ψ_d and ψ_m can be expressed as

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathit{\psi }}_{\mathrm{d}}=[{\mathit{d}}^{\mathrm{obs}}-{\mathit{d}}^{\mathrm{syn}}(\mathit{m}\left){]}^{\mathrm{T}}{\mathbf{W}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{W}}_{\mathrm{d}}\right[{\mathit{d}}^{\mathrm{obs}}-{\mathit{d}}^{\mathrm{syn}}\left(\mathit{m}\right)],\text{(14)}& {\displaystyle }{\mathit{\psi }}_{\mathrm{m}}=(\mathit{m}-{\mathit{m}}_{\mathrm{ref}}{)}^{\mathrm{T}}{\mathbf{W}}_{\mathrm{m}}^{\mathrm{T}}{\mathbf{W}}_{\mathrm{m}}(\mathit{m}-{\mathit{m}}_{\mathrm{ref}}),\end{array}$$

where vectors d^obs and d^syn are N_d dimensional vectors containing the observed and synthetic ERT data. ${\mathbf{W}}_{\mathrm{d}}\in {\mathbb{R}}^{N_{\mathrm{d}}\times N_{\mathrm{d}}}$ is the data weighting matrix, usually a diagonal matrix, whose elements are estimated based on the standard deviation of the noise; ${\mathbf{W}}_{\mathrm{m}}\in {\mathbb{R}}^{N_{\mathrm{m}}\times N_{\mathrm{m}}}$ is the regularization matrix, and m_ref is a N_m dimensional vector containing the logarithmic of reference conductivity model parameters.

3.2 Gauss–Newton method

To find a minimizer m_inv of the cost function in Eq. (12), we use a Gauss–Newton minimization approach. This is an iterative approach where after linearizing ψ(m_k) in the vicinity of a model m_k for a small model perturbation δm_k at the kth iteration, we obtain a set of linear equations or the normal equation

$$\begin{array}{}\text{(15)}& {\mathcal{H}}_k\mathit{\delta }{\mathit{m}}_k=-{\mathit{g}}_k,\end{array}$$

where the Gauss–Newton Hessian ${\mathcal{H}}_k\in {\mathbb{R}}^{N_{\mathrm{m}}\times N_{\mathrm{m}}}$ matrix is given by

$$\begin{array}{}\text{(16)}& {\mathcal{H}}_k=\mathfrak{R}\mathfrak{e}\mathit{\left\{}{\mathcal{J}}_k^{\mathrm{T}}{\mathbf{W}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{W}}_{\mathrm{d}}{\mathcal{J}}_k\mathit{\right\}}+\mathit{\beta }{\mathbf{W}}_{\mathrm{m}}^{\mathrm{T}}{\mathbf{W}}_{\mathrm{m}}\end{array}$$

and the gradient vector by

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\mathit{g}}_k& =-\mathfrak{R}\mathfrak{e}\left[{\mathcal{J}}_k^{\mathrm{T}}{\mathbf{W}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{W}}_{\mathrm{d}}\mathit{\right\{}{\mathit{d}}^{\mathrm{obs}}-{\mathit{d}}^{\mathrm{syn}}\mathit{\left\}}\right]\\ & +\mathit{\beta }{\mathbf{W}}_{\mathrm{m}}^{\mathrm{T}}{\mathbf{W}}_{\mathrm{m}}\mathit{\{}{\mathit{m}}_k-{\mathit{m}}_{\mathrm{ref}}\mathit{\}}.\end{array}\end{array}$$

The matrix 𝓙_k is the Jacobian matrix also known as the sensitivity or the Fréchet derivative matrix. Following Johnson et al. (2010), we implement a parallel conjugate-gradient least-square (CGLS) solver (Hestenes and Stiefel, 1952) to solve the normal Eq. (15). The normal equation is first reformulated into an equivalent least-squares problem and then solved for the model update δm_k using only the Jacobian matrix 𝓙_k without explicitly forming the Hessian matrix 𝓗_k. As N_m≫N_d for typical ERT surveys, this approach avoids the massive computational cost associated with the explicit formation of 𝓗_k. We therefore only need to compute 𝓙_k at each iteration, which is also a computationally intensive step. The effective computation of 𝓙_k is explained in the next subsection (Sect. 3.3).

The solution of Eq. (15) gives a model update vector δm_k at the kth iteration such that a new model

$$\begin{array}{}\text{(18)}& {\mathit{m}}_{k+\mathrm{1}}:={\mathit{m}}_k+\mathit{\alpha }\mathit{\delta }{\mathit{m}}_k\end{array}$$

decreases the cost function ψ, i.e., $\mathit{\psi }\left({\mathit{m}}_{k+\mathrm{1}}\right)<\mathit{\psi }\left({\mathit{m}}_k\right)$. There is a common practice to find step length α using a line search method (Nocedal and Wright, 2006); however, from our experience, a full step length, i.e., α=1, works well for the ERT data inversion.

In our implementation, the Gauss–Newton iteration usually starts with an average apparent conductivity model as an initial model and continues until the data misfit drops below a predefined tolerance level (defined by a chi-squared data misfit ${\mathit{\chi }}^{\mathrm{2}}={\mathit{\varphi }}_{\mathrm{d}}/N_{\mathrm{d}}\le \mathrm{1}$) or it reaches a predefined maximum number of iterations. Furthermore, a predefined input conductivity model can also be used as an initial model if it can be estimated from the existing geological and/or geophysical information. After the inversion converges, it provides the minimizer m_inv which represents the inverted conductivity model σ_inv=exp (m_inv).

3.3 Jacobian computation

Calculating 𝓙 at each iteration is one of the most computationally intensive steps in Gauss–Newton inversion. Therefore, an efficient 𝓙 drives the robustness of the developed inversion scheme and maximizes its scalability. The Jacobian matrix characterizes the change in the synthetic ERT data d^syn relative to a change in the model parameters m and is defined as the partial derivatives of d^syn with respect to m=ln σ; i.e.,

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\mathcal{J}}_{j,k}& ={\displaystyle \frac{\partial d_j^{\mathrm{syn}}}{\partial m_k}}={\mathit{\sigma }}_k{\displaystyle \frac{\partial d_j^{\mathrm{syn}}}{\partial {\mathit{\sigma }}_k}};j=\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{d}}\\ & \mathrm{and}k=\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{m}}.\end{array}\end{array}$$

An adjoint state method is used to effectively compute the Jacobian matrix. Let A and B be the source and sink electrodes located at r_A and r_B and M and N be the two measuring or receiving electrodes located at r_M and r_N. The potential difference between M and N can be obtained using the potential computed due to the source A at M and N, i.e., ϕ_AM and ϕ_AN, and the sink B at M and N, i.e., ϕ_BM and ϕ_BN. A simplified four-electrodes Wenner configuration is shown in Fig. 2. Note that, contrary to this configuration, electrodes can be placed in any other configurations, e.g., Schlumberger, pole–pole, pole–dipole, dipole–dipole, and gradient arrays (Dahlin and Zhou, 2004).

Figure 2 Wenner electrode configuration: four electrodes A, B, M, and N are deployed in-line with an equal spacing a between the two neighboring electrodes. Electrodes A and B act as the source and sink electrodes where current I is injected, while electrodes M and N act as the receiving electrodes measuring a potential difference Δϕ.

Using the expression derived in Appendix A (Eq. A9), an element of the Jacobian matrix is obtained by

$$\begin{array}{}\text{(20)}& {\mathcal{J}}_{j,k}={\displaystyle \frac{\partial d_j^{\mathrm{syn}}}{\partial m_k}}=({\mathit{\varphi }}_{\mathrm{A}}-{\mathit{\varphi }}_{\mathrm{B}}{)}^{\mathrm{T}}{\displaystyle \frac{\partial \mathbf{A}}{\partial m_k}}({\mathit{\varphi }}_{\mathrm{M}}-{\mathit{\varphi }}_{\mathrm{N}}),\end{array}$$

where ϕ_A and ϕ_B are potential distributions due to the source A and sink B and ϕ_M and ϕ_N are potential distributions due to pseudo-source and pseudo-sink, respectively, at the receiving electrodes M and N.

Assuming ${\mathit{\varphi }}_{\mathrm{S}}={\mathit{\varphi }}_{\mathrm{A}}-{\mathit{\varphi }}_{\mathrm{B}}$ and ${\mathit{\varphi }}_{\mathrm{R}}={\mathit{\varphi }}_{\mathrm{M}}-{\mathit{\varphi }}_{\mathrm{N}}$, where ϕ_S is the net potential due to the source and sink electrodes A and B, and ϕ_R is the net potential due to source and sink assumed to be located, respectively, at the receiving electrodes M and N, this results in

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\partial d_j^{\mathrm{syn}}}{\partial m_k}}={\mathit{\varphi }}_{\mathrm{S}}^{\mathrm{T}}{\displaystyle \frac{\partial \mathbf{A}}{\partial m_k}}{\mathit{\varphi }}_{\mathrm{R}}.\end{array}$$

Equation (21) is valid for any numbers of the source/sink electrodes and the receiving electrodes as well as any of their combinations. One only needs to compute the net potentials due to the source/sink electrodes ϕ_S and receiving electrodes ϕ_R using the superposition principle.

4 Parallel implementation within PFLOTRAN

The parallel implementation strategy for the ERT forward modeling and inversion capabilities follows the existing PFLOTRAN’s strategy for flow and transport process models. Forward ERT modeling is implemented by extending PFLOTRAN's hierarchy of Fortran classes and data structures to create a new ERT process model that constructs and solves the linear systems of equations described in Sect. 2. The implementation of inversion, as described in Sect. 3, is completely new to PFLOTRAN and requires all the necessary steps (or calculations) for ERT inversion, such as the cost-function computation, regularization, Jacobian computation, solution to the dense Gauss–Newton normal equation, and model update.

To execute simulations in parallel, domain decomposition is employed to distribute the discrete computational domain across (computer) processes using a division along the three principal axes for structured grids (calculated by PETSc or explicitly defined by the user), or in the case of unstructured grids, the ParMETIS parallel graph partitioner (Karypis and Schloegel, 2013) is employed to divide the domain. All interprocess communication is carried out with MPI (Message Passing Interface), and PETSc's parallel data structures facilitate much of the MPI communication.

Efficient partitioning and the overlapping of communication with computation (e.g., calculating floating point operations while data are being transferred between processes) helps improve parallel performance. The latter helps mask communication costs. The goal is to maximize the time spent in computation and minimize interprocess communication and enforce load balance, the even distribution of computation across processes employed. Efficient partitioning can be accomplished by dividing the domain evenly and minimizing the number and size of (MPI) messages being passed between processes. However, an even distribution does not necessarily guarantee load balance as certain regions of the domain with time-consuming operations (e.g., localized physics, intense input/output or I/O) may require more computation. Thus, efficient partitioning is not trivial.

PFLOTRAN's parallel file I/O is through HDF5 (Hierarchical Data Format version 5), a platform-independent binary file format that leverages MPI-IO operations to read/write a single file using multiple (potentially thousands of) processes (The HDF Group, 2022). PFLOTRAN writes ASCII-formatted observation files (files storing ERT data at specific locations in space) locally on the processes owning the observation points. Thus, the writing of observation files requires no parallel communication. For more comprehensive details on the parallel implementation strategy employed within PFLOTRAN, interested readers are referred to Hammond et al. (2012, 2014).

5 Modeling benchmarking results

To examine the accuracy of the PFLOTRAN ERT process model, simulation results are compared to solutions obtained using well-established analytical and numerical methods. Two cases are considered for the comparison: three-layer earth and 3-D earth resistivity models. All benchmarking simulations were performed on the Deception supercomputer housed at the Pacific Northwest National Laboratory. It is composed of 96 compute nodes where each node has 64-core AMD EPYC 7502 processors running at 2.5 GHz boost to 3.35 GHz with 256 GB of memory.

5.1 Layered earth model

Let us consider the layered earth or 1D model shown in Fig. 3. The model is composed of three layers, whose resistivity (ρ) and thickness (t), i.e., [ρ, t] are [100 Ωm, 30 m], [300 Ωm, 30 m], and [10 Ωm, half space] from the top to bottom. The dimension of the model is $\mathrm{400}\times \mathrm{400}\times \mathrm{100}$ m³. The model boundaries are extended to ±3000 m in the x and y directions and 3000 m in the z direction to accommodate zero Dirichlet boundary conditions at the side and bottom boundaries. These boundary extensions are not shown in Fig. 3. Vertical electrical sounding (VES) data are simulated at the center of the model considering a Wenner configuration. This configuration consists of four electrodes A, B, M, and N, which are deployed in-line with an equal spacing a between the two neighboring electrodes (Fig. 2). Electrodes A and B again act as the current source and sink electrodes, while electrodes M and N act as the potential receiving electrodes. Starting with a smaller value of a for a shallow-depth investigation, the value is progressively increased to examine deeper depths. In our case, we start with a minimum value a_min=8 m and progressively increase it to a maximum a_max=132 m following $a=a_{\min }+i\mathit{\delta }a$, where δa=4 m and $i=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{31}$. This results in N_d=32 VES sounding recordings using N_e=128 electrodes.

Figure 3 Vertical cross section of a three-layer earth resistivity model used for benchmarking PFLOTRAN results against the analytical solutions.

Figure 4 Apparent resistivity ρ^a responses computed for the layered earth resistivity model in Fig. 3 (upper plot). They are calculated using the analytical method (solid black lines) and PFLOTRAN (filled red circles). The lower plot (blue line) shows the relative difference (in percentage) between the analytical and PFLOTRAN results.

To simulate the VES data, the model is discretized using a mix of uniform and non-uniform cell sizes in the x, y, and z directions. The main computational domain is discretized with a uniform grid spacings of 2 m in the x and y directions and 1 m in the z direction. The boundary regions are discretized with severely stretched non-uniform grid spacings following Jaysaval et al. (2014). The discretization results in a grid with $\mathrm{240}\times \mathrm{240}\times \mathrm{120}$ cells. Therefore, using the FV method to simulate electric potentials leads to a system of linear equations with 6 912 000 DOFs. PFLOTRAN simulated the VES data using 512 processes on Deception. The total simulation time was 45 s.

The simulated VES results are compared to analytical solutions obtained from SimPEG (Cockett et al., 2015). Figure 4 shows the comparison of the apparent resistivity ρ^a responses. Analytical and PFLOTRAN results are shown by the solid black line and filled red circles, respectively. We observe that both responses agree very well. To make the comparison more quantitative, we calculate the relative difference between the two responses and show the difference in Fig. 4 (blue line). The relative difference varies from 0.02 % to 1.3 %. We also calculate the average relative difference using

$$\begin{array}{}\text{(22)}& \mathit{ϵ}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{d}}}}\sum _i^{N_{\mathrm{d}}}{\displaystyle \frac{{{\mathit{\rho }}_{\mathrm{1}}^{\mathrm{a}}}_i-{{\mathit{\rho }}_{\mathrm{2}}^{\mathrm{a}}}_i}{{{\mathit{\rho }}_{\mathrm{1}}^{\mathrm{a}}}_i}},\end{array}$$

where ρ^a with subscripts 1 and 2, respectively, represents the apparent resistivity computed using PFLOTRAN and the analytical method from SimPEG. The average relative difference ϵ is 0.18 %. These relatively small differences demonstrate the high degree of accuracy of the PFLOTRAN ERT process model for the 1D model. Note that the apparent resistivity was computed using

$$\begin{array}{}\text{(23)}& {\mathit{\rho }}^{\mathrm{a}}={\displaystyle \frac{\mathrm{\Delta }\mathit{\varphi }}{IG}},\end{array}$$

where $\mathrm{\Delta }\mathit{\varphi }={\mathit{\varphi }}_{\mathrm{AM}}-{\mathit{\varphi }}_{\mathrm{AN}}-{\mathit{\varphi }}_{\mathrm{BM}}+{\mathit{\varphi }}_{\mathrm{BN}}$ and G is a geometric factor given as

$$\begin{array}{}\text{(24)}& G={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\left({\displaystyle \frac{\mathrm{1}}{r_{\mathrm{AM}}}}-{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{AN}}}}-{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{BM}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{BN}}}}\right),\end{array}$$

where distances between the source electrode A and the receiving electrodes M and N are, respectively, r_AM and r_AN, while those between the sink electrode B and the receiving electrodes M and N are, respectively, r_BM and r_BN. For the Wenner configuration, r_AM=a, r_AN=2a, r_BM=2a, and r_BN=a. Therefore, geometric factor $G=\mathrm{1}/\mathrm{2}\mathit{\pi }a$ was used for calculating ρ^a using Eq. (23).

5.2 Three-dimensional earth model

In the previous example, PFLOTRAN accuracy for simulating ERT responses was validated for a simplistic 1D earth model. In the real world, however, the earth models are 3-D with varying degrees of heterogeneity. Therefore, it is important to validate the accuracy of numerical results computed using PFLOTRAN against results computed using well-established methods for 3D models.

Figure 5 Three-dimensional earth resistivity model composed of two anomalous blocks embedded in a homogeneous background. ERT profiling data are simulated along a survey profile at $y=z=\mathrm{0}$ m from $x=-\mathrm{16}$ m to x=16 m using the Wenner configuration with increasing electrode spacing a from 1 to 6 m.

We consider the 3D model shown in Fig. 5. The model dimension is $\mathrm{32}\times \mathrm{16}\times \mathrm{8}$ m³. It has a homogeneous background of conductivity 0.002 Sm⁻¹ and includes two anomalous blocks: a conductive block (0.2 Sm⁻¹) and a resistive block (0.0002 Sm⁻¹). Both anomalous blocks have a dimension of $\mathrm{2}\times \mathrm{2}\times \mathrm{2}$ m³ and are separated by 4 m in the x direction. As in the previous case, the model boundaries are extended to ±1000 m in the x and y directions and 1000 m in the z direction to accommodate zero Dirichlet boundary conditions at the side and bottom boundaries. ERT profiling data are simulated along a survey profile (yellow line) at $y=z=\mathrm{0}$ m from $x=-\mathrm{16}$ m to x=16 m using the Wenner configuration with increasing electrode spacing a from 1 to 6 m. Unlike the VES performed for the 1D model where data were simulated at a single observation point at the center of the model, the ERT profiling, performed for the 3D model, records data at multiple observation points along the survey profile. The profiling, therefore, helps in investigating lateral conductivity variations in the model along the x direction. In this case, it results in N_d=129 ERT data recording by using N_e=32 electrodes placed uniformly at every 1 m from $x=-\mathrm{15.5}$ m to x=15.5 m.

Figure 6 Apparent resistivity ρ^a responses computed for the 3D earth resistivity model in Fig. 5 using the Wenner configuration with spacings: (a) a=1 m, (b) a=2 m, (c) a=3 m, (d) a=4 m, (e) a=5 m, and (f) a=6 m (upper plots). They are calculated using PFLOTRAN (solid black lines) and E4D (filled red circles) for an ERT profiling along a survey profile at $y=z=\mathrm{0}$ from $x=-\mathrm{16}$ m to x=16 m. The lower plots (blue lines) show the relative difference (in percentage) between PFLOTRAN and E4D results.

The main computational domain is discretized with a uniform grid spacing of 0.2 m in all the three directions, while the boundary domain is again discretized with severely stretched non-uniform grid spacings. The model discretization leads to a grid with $\mathrm{200}\times \mathrm{120}\times \mathrm{60}$ FV cells and, hence, 1 440 000 DOFs in the corresponding system of linear equations. PFLOTRAN simulated the ERT profiling data using 128 processes on Deception. The total simulation time was 9 s.

Figure 6a–f (upper curves) show the apparent resistivity ρ^a responses computed for the 3D model in Fig. 5 using the ERT profiling with the Wenner configuration with electrode spacings of 1, 2, 3, 4, 5, and, 6 m. The solid black lines show the responses computed using PFLOTRAN. All responses exhibit the characteristics of the conductive and resistive anomalous blocks embedded in a homogeneous background of conductivity 0.002 Sm⁻¹ or resistivity 500 Ωm; see, e.g., the low apparent resistivity on the left and high apparent resistivity on the right sides of the plots. These responses are compared to responses computed using E4D (Johnson et al., 2010; Johnson and Wellman, 2015). In E4D, an unstructured-mesh FE method is implemented to perform the forward modeling. The filled red circles show the responses computed using E4D. The relative differences between the two codes are shown by the lower blue lines, which vary from 0.02 % to 2.5 %. The average relative difference ϵ computed using Eq. (22) is 0.77 %. These relatively small differences imply a good agreement between responses computed using both codes.

6 Inversion results

To examine the efficiency of PFLOTRAN to invert ERT data, we consider a model modified from Fig. 5 by cropping the model for $\left|x\right|>\mathrm{8}$ m. The modified model, therefore, has the dimension of $\mathrm{16}\times \mathrm{16}\times \mathrm{8}$ m³. The model is discretized using the same grid spacings as in the previous 3D modeling example, which resulted in $\mathrm{120}\times \mathrm{120}\times \mathrm{60}=\mathrm{864}\mathrm{000}$ FV cells. We computed synthetic ERT data for three survey lines located at the ground surface at $y=-\mathrm{5}$ m, y=0 m, and y=5 m from $x=-\mathrm{8}$ m to x=8 m. Electrodes are placed uniformly with an electrode spacing of 1 m along the three survey lines from $x=-\mathrm{7.5}$ m to x=7.5 m, i.e., 16 electrodes along each survey line. Using these N_e=48 electrodes, ERT data were simulated considering different combinations of a set of 4 electrodes, where 2 electrodes acted as the current and the remaining 2 as the potential receiving electrodes. These electrode combinations result in N_d=2070 simulated ERT data. PFLOTRAN simulated these ERT data using 128 processes on Deception. The total simulation time was 6 s. In this case, ERT data were simulated in the form of the resistance obtained using $\frac{\mathrm{\Delta }\mathit{\varphi }}{I}$, which basically equals to the potential difference Δϕ for unit current (I=1 A).

The simulated ERT data were then inverted using PFLOTRAN starting with a model having a half-space conductivity of 0.0022 Sm⁻¹. The starting model conductivity was obtained by averaging the apparent conductivity computed from the given ERT data. The previously described modeling grid also acted as the inversion grid; i.e., each cell in the model was used as an inversion parameter. Therefore, we needed to invert for N_m=864 000 model parameters. PFLOTRAN inverted the ERT data using 512 processes on Deception. The total inversion time was 74 s.

Figure 7 Inversion result of the simulated ERT data for a model modified from Fig. 5 by cropping the model for $\left|x\right|>\mathrm{8}$ m.

Figure 8 A measure of data misfit as the inversion progresses: χ² versus number of iteration.

Figure 9 Histogram of data-fit error, which is defined as ratio of the residual resistance (difference in the observed and predicted resistances) and observed resistance expressed as a percentage.

Figure 7 shows the inverted conductivity model. The inverted model agrees well with the true model: the inversion recovered the conductive and resistive anomalies in the model. The inverted resistivity, however, varies smoothly in the model because a smooth regularization constraint is applied by minimizing the difference in the cell conductivity and the conductivity of its neighboring cells. Further discussion on the implementation of our regularization schemes is beyond the scope of the current paper.

The initial chi-squared data misfit χ² was 14.8. As the inversion progressed, χ² decreased as shown in Fig. 8. The inversion took 11 iterations to satisfy the convergence criteria defined by χ²≤0.9. To further examine the efficiency of PFLOTRAN to fit the individual ERT data, a histogram of data-fit error is also shown in Fig. 9. The data-fit error is defined as ratio of the residual resistance (difference in the observed and predicted resistances) and observed resistance expressed as a percentage. The histogram implies that the ERT data were fitted comparatively well.

7 Scalability tests

Scalability tests demonstrate parallel performance or how efficiently a simulator utilizes additional processes to either speed up a simulation of fixed problem size (strong scaling) or maintain speed as the problem size and number of processes increase proportionally (weak scaling). For ideal strong scaling, the speedup is proportional to the increase in the number of processes (e.g., doubling the number of processes cuts the run time in half), while ideal weak scaling results in an unaltered or constant run time. In either (ideal) case, the additional computing capacity is efficiently utilized. Ideal strong scaling is often termed linear speedup, whereas superlinear speedup is better than ideal (e.g., doubling the number of processes reduces the run time to less than half). Superlinear performance is not uncommon, as access to an increased cache (efficient, high-performance memory on the processor core) can generate superlinear speedup.

As discussed earlier, there are two computationally intensive steps to invert ERT data: forward modeling and Jacobian computation. To demonstrate the parallel performance of PFLOTRAN for these two steps, large-scale strong scaling tests are reported in this section by using up to 1024 processes on the Deception supercomputer or up to 131 072 processes on the Theta supercomputer. Configuration details for Deception were previously given in the “Modeling benchmarking results” section (Sect. 5). On the other hand, Theta is a leadership-class supercomputer, ranked within the top 100 fastest computers in the world and housed at the Argonne Leadership Computing Facility (ALCF) at Argonne National Laboratory (https://www.alcf.anl.gov/alcf-resources/theta, last access: 27 January 2023). With a peak performance of 11.7 PFLOPS, it is composed of 4392 compute nodes connected through a Cray Aries network. Each node has 64-core Intel Xeon Phi 7230 processors running at 1.3 GHz with 192 GB of memory. Therefore, there are 281 088 processes available on Theta. The PFLOTRAN compilations on Deception and Theta were linked to Intel's Math Kernel Library (MKL). The code leverages solely the distributed-memory parallel capability on the machines and not the shared-memory multi-threading.

7.1 Forward modeling scalability

The three-layer earth model in Fig. 3 is first rediscretized to generate several test matrices for the analysis by considering different cell sizes in the x, y, and z directions. These test matrices result in six different modeling scenarios from coarser to finer grids. Table 1 completely specifies details of the test matrices M2, M7, M19, M35, M64, and M122 along with the number of DOFs in the governing matrix equations that ranges between 2 ×10⁶ and 122 ×10⁶. The naming of these matrices is based on the number of DOFs they represent; e.g., M19 represents a matrix equation with about 19 ×10⁶ DOFs. To the best of our knowledge, the modeling with 122 ×10⁶ DOFs stands as the largest ERT modeling reported in the geophysical literature. The next largest 3D ERT modeling problem reported in the literature was composed of about 10 ×10⁶ DOFs (Ren et al., 2018), while most others were often composed of less than 1 ×10⁶ DOFs (Günther et al., 2006; Rücker et al., 2006; Johnson et al., 2010; Johnson and Wellman, 2015). Having the capability to solve for up to a 100 ×10⁶ DOFs efficiently in parallel is a crucial requirement for large-scale ERT modeling and inversion problems in the future if model dimensions are significantly increased and/or finer grids are employed.

Table 1 Test matrices obtained by rediscretizing the three-layer earth model in Fig. 3 using different cell sizes dx, dy, and dz in the x, y, and z directions. N_x, N_y, and N_z are, respectively, the number of FV cells in the x, y, and z directions. Note that the cell sizes are given only for the main computational domain; the boundary domain is discretized using severely stretched non-uniform cell sizes.

Figure 10 PFLOTRAN scalability for simulating ERT data for the model in Fig. 3 for N_e=128 electrodes. The test is performed on the Theta supercomputer considering six test matrices M2, M7, M19, M35, M64, and M122 with the respective number of DOFs of about 2, 7, 19, 35, 64, and 122 ×10⁶ (see Table 1 for details).

Considering the six test matrices above, we performed the forward modeling scalability test on Theta for simulating ERT data for N_e=128 electrodes, the same number of electrodes as in the three-layer benchmarking example. Figure 10 illustrates the PFLOTRAN ERT process model strong scaling for each test case (wall-clock time versus the number of processes employed for a fixed problem size). The wall-clock time is shown on the vertical axis, while the horizontal axis shows the number of processes; both axes are shown in logarithmic base-2 scale. The line labeled “M19 Ideal” is the plot of ideal performance based on the 64 process run time for the M19 simulation:

$$\begin{array}{}\text{(25)}& T_n={\displaystyle \frac{T_{\mathrm{64}}\times \mathrm{64}}{n}},\end{array}$$

where T_n is the time on n processes and T₆₄ is the time on 64 processes.

For all problem sizes, the wall-clock times flatten at higher numbers of processes. This phenomenon is expected as the amount of work per process (DOFs per process) decreases at higher process counts and the cost of communication between processes begins to dominate. The M19 performance begins to deviate from ideal at ∼2048 processes, which equates to ∼9000 DOFs per process. Similar behavior is observed for many of the other problem sizes. Hammond et al. (2014) demonstrate that for groundwater flow and reactive transport simulation PFLOTRAN runs most efficiently when at least 10 000 DOFs are assigned to each process. These results show that this 10 000 DOF threshold also applies for PFLOTRAN ERT simulations. That is, PFLOTRAN's scalability on the supercomputer should be more ideal when the number of DOFs per process is greater than 10 000.

Another factor limiting parallel performance is the limited scalability of linear Krylov solvers at high process counts. Mills et al. (2009) demonstrates that the vector inner products performed within Krylov solvers become increasingly expensive as tens of thousands of processes are employed to solve a problem. These inner products (or MPI global reduction operations) require repeated global synchronization across all processes involved. As more processes are employed to solve the problems (requiring communication across a growing network of compute nodes), the transfer of information for each vector inner product is more expensive time-wise (i.e., increased number of hops within the high-speed network, longer distance traveled). Figure 10 reveals a similar phenomenon.

7.2 Jacobian computation scalability

To test the scalability of PFLOTRAN for the Jacobian computation, we consider a Jacobian matrix computed for inverting the ERT data in the “Inversion result” section (Sect. 6). Recalling that N_d=2070 ERT data and N_m=864 000 model parameters for the inversion example, the dimension of the Jacobian matrix 𝓙 is therefore 2070×864 000. The Jacobian matrix is computed at each inversion iteration after the forward modeling run and used in the normal equation (Eq. 15) to get a model update δm.

Figure 11 PFLOTRAN strong scaling for computing the Jacobian for the ERT inversion example in Fig. 7. The dimension of the Jacobian matrix 𝓙 is 2070×864 000 and is distributed uniformly over each computing process. PFLOTRAN exhibits superlinear performance for larger numbers of processes.

The Jacobian scalability test is performed on the Deception supercomputer. Figure 11 illustrates the strong scalability of PFLOTRAN wall-clock time for computing 𝓙. It shows that PFLOTRAN exhibits linear (and sometimes superlinear) scalability for computing 𝓙 with increasing numbers of processes. Consequently, the total wall-clock time for computing 𝓙 using 1024 processes reduces to 0.5 s compared to 887 s when using a single process on Deception. This ideal performance in computing 𝓙 is most likely due to even load balance across processes and the lack of MPI communication in its construction. Moreover, as noted before, the superlinear performance is most likely due to access to increased cache with an increasing number of processes.

For large-scale ERT inversion problems, the linear scaling for computing 𝓙 can be exploited to reduce the overall inversion wall-clock time even if PFLOTRAN exhibits deteriorating scalability for the forward modeling on a very large number of processes. For example, for an ERT inversion problem where the 𝓙 calculation time is much higher than the forward modeling time, increasing the number of processes to reduce the overall inversion wall-clock time is recommended, even when forward modeling exhibits a significant departure from the linear scalability. One has to find a balance between the two, and the most time-consuming operation should govern.

8 Discussions

In the past 20 years, a few open-source 3D ERT modeling and inversion codes were developed to handle moderate-scale ERT problems (e.g., ResIPy, SimPEG, pyGIMLI, and E4D). Although some of these codes can practically tackle most present ERT modeling and inversion problems having less than 10 ×10⁶ DOFs, they still lack full-scale parallel implementations to efficiently distribute numerical computations over massive supercomputers. The limitation can also lead to computational challenges for large-scale ERT problems with 100 ×10⁶ or more DOFs, which may not be solved efficiently with the existing codes. The 3D ERT modeling and inversion capability developed within PFLOTRAN removes this limitation. To examine the accuracy of the forward modeling capability, numerical results from PFLOTRAN were compared to analytical solutions for a layered earth model and numerical results obtained from E4D for a 3D earth model. The average relative differences stayed within 1 %, which demonstrated the high degree of accuracy of the PFLOTRAN modeling capability. Furthermore, large-scale strong scalability tests, using up to 131 072 processes, were performed on a leadership class supercomputer to demonstrate efficient distribution of computational workloads over massive supercomputers. The scalability tests demonstrate that computational times for the Jacobian computation scale linearly with increasing number of processes. For the forward modeling, PFLOTRAN showed a linear to near-linear scaling behavior up to 10 000 DOFs per process, beyond which performance tends to degrade.

Due to the efficient parallel performance of PFLOTRAN observed above, a large-scale ERT simulation problem with 122 ×10⁶ DOFs and 122 electrodes could be solved in about 2 min by using 32 768 processes as compared to taking about 1 h with 512 processes. Combining this forward modeling scalability along with the ideal scaling of the Jacobian computation, we can aim for inverting moderate-scale ERT problems within minutes and large-scale problems within hours with the access to massive supercomputers. The moderate-scale inversion example of Sect. 6 demonstrated that inversion of N_d=2070 data and N_m=864 000 model parameters could be performed in about 1 min by using 512 processes.

Finally, there were two key motivations for implementing ERT capabilities in PFLOTRAN. The first motivation was to utilize the massive parallel infrastructure implemented within PFLOTRAN for the flow and transport simulations as illustrated by Hammond et al. (2012, 2014) for simulating subsurface problems composed of billions of DOFs utilizing hundreds of thousands of computing processors. This infrastructure made PFLOTRAN a natural fit to achieve our goal of massive parallel implementation for the ERT forward modeling and inversion capabilities. Researchers familiar with PFLOTRAN can now have immediate access to the added geophysical ERT capabilities with minimal learning efforts.

The second motivation was to have a native implementation of ERT modeling and inversion capabilities along with flow and transport simulations. The new capabilities will allow us to perform fully coupled hydro-geophysical simulations and inversion. The idea behind the joint hydro-geophysical simulation is to use known flow and transport parameters, such as porosity and intrinsic permeability, and simulate for hydrological variables, e.g., liquid saturation and solute concentration. Using petrophysical relationships (e.g., Archie's law or its modifications; Slater, 2007), these parameters and variables can be mapped into electrical conductivity, which subsequently can be used for simulating ERT responses. In our previous work, a joint hydro-geophysical simulation capability was developed by combining the capabilities of PFLOTRAN and E4D (Johnson et al., 2017). However, it was not straightforward to use PFLOTRAN–E4D for joint hydro-geophysical simulations. This is because the implementation was based on obtaining flow and transport parameters and variables from PFLOTRAN using a structured mesh, which then needed to be interpolated on an unstructured tetrahedral mesh for ERT simulations by E4D. The process was cumbersome and had the potential for errors. PFLOTRAN’s native implementation of ERT avoids employing any kind of interpolation because both flow and transport and ERT simulations use the same modeling grid.

Furthermore, there are growing interests in joint hydro-geophysical inversion, i.e., estimating the distribution of subsurface flow and transport parameters, from time-lapse ERT data; see, e.g., Mboh et al. (2012), Camporese et al. (2015), Rücker et al. (2017), González-Quirós and Comte (2021), and Pleasants et al. (2022). The inclusion of the ERT simulation capability in PFLOTRAN fulfills the requirement of the forward modeling necessary for native implementation of hydro-geophysical inversion. Fully coupled hydro-geophysical simulation and inversion capabilities are currently under development within PFLOTRAN and will be a focus of our forthcoming papers.

9 Conclusions

Two new capabilities – forward modeling and inversion of geophysical ERT data – are developed within PFLOTRAN, a massively parallel open-source, state-of-the-art, multi-phase, multi-component subsurface flow and transport simulation code. The capabilities are accurate, robust, and highly scalable on HPC platforms. The accuracy of the forward modeling capability was demonstrated for layered earth and 3D earth conductivity models by comparing our numerical results against well-established analytical and numerical results. The average relative differences stayed within 1 %, demonstrating the high degree of accuracy of the PFLOTRAN ERT process model. The inversion capability was demonstrated for a 3D synthetic case by recovering the conductivity of the model. Large-scale scalability tests illustrated that the forward modeling for models with 100 ×10⁶ DOFs can be performed in a few minutes; e.g., ERT modeling for a model with 122 ×10⁶ DOFs and 122 electrodes was performed in 128 s by using 32 768 processes on a supercomputer. The Jacobian computation, a key ingredient for the inversion, exhibited ideal or linear scaling, and as a result a Jacobian of 2070×864 000 dimension was computed in 0.5 s using 1024 processes compared to 887 s using a single process. Moreover, integration of geophysical ERT capabilities within PFLOTRAN opens the door to perform coupled hydro-geophysical modeling and inversion, which is a focus of our future research.

Appendix A: Jacobian computation

Let us assume that p_M and p_N are two interpolation vectors and ϕ_A and ϕ_B are potential distributions due to the source A and sink B (Fig. 2) such that

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\mathit{\varphi }}_{\mathrm{AM}}& ={\mathit{p}}_{\mathrm{M}}^{\mathrm{T}}{\mathit{\varphi }}_{\mathrm{A}},{\mathit{\varphi }}_{\mathrm{AN}}={\mathit{p}}_{\mathrm{N}}^{\mathrm{T}}{\mathit{\varphi }}_{\mathrm{A}},{\mathit{\varphi }}_{\mathrm{BM}}={\mathit{p}}_{\mathrm{M}}^{\mathrm{T}}{\mathit{\varphi }}_{\mathrm{B}},\\ & \mathrm{and}{\mathit{\varphi }}_{\mathrm{BN}}={\mathit{p}}_{\mathrm{N}}^{\mathrm{T}}{\mathit{\varphi }}_{\mathrm{B}}.\end{array}\end{array}$$

Therefore, the jth simulated potential difference between the receiving electrodes M and N can be obtained as

$$\begin{array}{}\text{(A2)}& d_j^{\mathrm{syn}}={\mathit{\varphi }}_{\mathrm{AM}}-{\mathit{\varphi }}_{\mathrm{AN}}-{\mathit{\varphi }}_{\mathrm{BM}}+{\mathit{\varphi }}_{\mathrm{BN}}.\end{array}$$

Using Eq. (A1) in Eq. (A2), we get

$$\begin{array}{}\text{(A3)}& d_j^{\mathrm{syn}}={\mathit{p}}_{\mathrm{M}}^{\mathrm{T}}({\mathit{\varphi }}_{\mathrm{A}}-{\mathit{\varphi }}_{\mathrm{B}})-{\mathit{p}}_{\mathrm{N}}^{\mathrm{T}}({\mathit{\varphi }}_{\mathrm{A}}-{\mathit{\varphi }}_{\mathrm{B}}).\end{array}$$

Taking the derivative of Eq. (A3) with respect to conductivity m_k of kth cell yields

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\partial d_j^{\mathrm{syn}}}{\partial m_k}}={\mathit{p}}_{\mathrm{M}}^{\mathrm{T}}{\displaystyle \frac{\partial }{\partial m_k}}({\mathit{\varphi }}_{\mathrm{A}}-{\mathit{\varphi }}_{\mathrm{B}})-{\mathit{p}}_{\mathrm{N}}^{\mathrm{T}}{\displaystyle \frac{\partial }{\partial m_k}}({\mathit{\varphi }}_{\mathrm{A}}-{\mathit{\varphi }}_{\mathrm{B}}).\end{array}$$

Furthermore, taking the derivative of Eq. (10) with respect to m_k yields

$$\begin{array}{}\text{(A5)}& {\displaystyle \frac{\partial \mathbf{A}}{\partial m_k}}\mathit{\varphi }+\mathbf{A}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial m_k}}=\mathbf{0},\end{array}$$

where we used $\partial \mathit{s}/\partial m_k=\mathrm{0}$ as the source does not depend on the conductivity.

If A is invertible, rearranging Eq. (A5) gives

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial m_k}}=-{\mathbf{A}}^{-\mathrm{1}}{\displaystyle \frac{\partial \mathbf{A}}{\partial m_k}}\mathit{\varphi }.\end{array}$$

Using Eq. (A6) in Eq. (A4), we get

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}{\displaystyle \frac{\partial d_j^{\mathrm{syn}}}{\partial m_k}}& =-{\mathit{p}}_{\mathrm{M}}^{\mathrm{T}}{\mathbf{A}}^{-\mathrm{1}}{\displaystyle \frac{\partial \mathbf{A}}{\partial m_k}}({\mathit{\varphi }}_{\mathrm{A}}-{\mathit{\varphi }}_{\mathrm{B}})\\ & +{\mathit{p}}_{\mathrm{N}}^{\mathrm{T}}{\mathbf{A}}^{-\mathrm{1}}{\displaystyle \frac{\partial \mathbf{A}}{\partial m_k}}({\mathit{\varphi }}_{\mathrm{A}}-{\mathit{\varphi }}_{\mathrm{B}}).\end{array}\end{array}$$

As $\partial d_j^{\mathrm{syn}}/\partial m_k$ is a scalar quantity, using the property x^Ty=y^Tx for two vectors x and y, Eq. (A7) yields

$$\begin{array}{}\text{(A8)}& {\displaystyle \frac{\partial d_j^{\mathrm{syn}}}{\partial m_k}}=-({\mathit{\varphi }}_{\mathrm{A}}-{\mathit{\varphi }}_{\mathrm{B}}{)}^{\mathrm{T}}{\displaystyle \frac{\partial \mathbf{A}}{\partial m_k}}({\mathbf{A}}^{-\mathrm{1}}{\mathit{p}}_{\mathrm{M}}-{\mathbf{A}}^{-\mathrm{1}}{\mathit{p}}_{\mathrm{N}}),\end{array}$$

where A^T=A is used due to its symmetry.

If the same interpolation scheme is used for placing the current electrodes within the discretized cells and recording the potential at electrodes from the cell potentials, we will have ${\mathbf{A}}^{-\mathrm{1}}{\mathit{p}}_{\mathrm{M}}=-{\mathit{\varphi }}_{\mathrm{M}}$ and ${\mathbf{A}}^{-\mathrm{1}}{\mathit{p}}_{\mathrm{N}}=-{\mathit{\varphi }}_{\mathrm{N}}$. Therefore,

$$\begin{array}{}\text{(A9)}& {\displaystyle \frac{\partial d_j^{\mathrm{syn}}}{\partial m_k}}=({\mathit{\varphi }}_{\mathrm{A}}-{\mathit{\varphi }}_{\mathrm{B}}{)}^{\mathrm{T}}{\displaystyle \frac{\partial \mathbf{A}}{\partial m_k}}({\mathit{\varphi }}_{\mathrm{M}}-{\mathit{\varphi }}_{\mathrm{N}}).\end{array}$$