PFLOTRAN-SIP: A PFLOTRAN Module for Simulating Spectral-Induced Polarization of Electrical Impedance Data

by   B. Ahmmed, et al.
Baylor University

Spectral induced polarization (SIP) is a non-intrusive geophysical method that is widely used to detect sulfide minerals, clay minerals, metallic objects, municipal wastes, hydrocarbons, and salinity intrusion. However, SIP is a static method that cannot measure the dynamics of flow and solute/species transport in the subsurface. To capture these dynamics, the data collected with the SIP technique needs to be coupled with fluid flow and reactive-transport models. To our knowledge, currently, there is no simulator in the open-source literature that couples fluid flow, solute transport, and SIP process models to analyze geoelectrical signatures in a large-scale system. A massively parallel simulation framework (PFLOTRAN-SIP) was built to couple SIP data to fluid flow and solute transport processes. This framework built on the PFLOTRAN-E4D simulator that couples PFLOTRAN and E4D, without sacrificing computational performance. PFLOTRAN solves the coupled flow and solute transport process models to estimate solute concentrations, which were used in Archie's model to compute bulk electrical conductivities at near-zero frequency. These bulk electrical conductivities were modified using the Cole-Cole model to account for frequency dependence. Using the estimated frequency-dependent bulk conductivities, E4D simulated the real and complex electrical potential signals for selected frequencies for SIP. The PFLOTRAN-SIP framework was demonstrated through a synthetic tracer-transport model simulating tracer concentration and electrical impedances for four frequencies. Later, SIP inversion estimated bulk electrical conductivities by matching electrical impedances for each specified frequency. The estimated bulk electrical conductivities were consistent with the simulated tracer concentrations from the PFLOTRAN-SIP forward model.



There are no comments yet.


page 10

page 11

page 12

page 13


Geometrically reduced modelling of pulsatile flow in perivascular networks

Flow of cerebrospinal fluid in perivascular spaces is a key mechanism un...

Numerical simulation of oxidation processes in a cross-flow around tube bundles

An oxidation process is simulated for a bundle of metal tubes in a cross...

Learned coupled inversion for carbon sequestration monitoring and forecasting with Fourier neural operators

Seismic monitoring of carbon storage sequestration is a challenging prob...

DuMu^x 3 – an open-source simulator for solving flow and transport problems in porous media with a focus on model coupling

We present version 3 of the open-source simulator for flow and transport...

A 1D-0D-3D coupled model for simulating blood flow and transport processes in breast tissue

In this work, we present mixed dimensional models for simulating blood f...

Less can be more: Insights on the role of electrode microstructure in redox flow batteries from 2D direct numerical simulations

Understanding how to structure a porous electrode to facilitate fluid, m...

PorePy: An Open-Source Simulation Tool for Flow and Transport in Deformable Fractured Rocks

Fractures are ubiquitous in the subsurface and strongly affect flow and ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.


Spectral induced polarization (SIP) is a non-intrusive geophysical method that gleans information in the form of chargeability (the ability of a material to retain charge) in the time domain or its phase shift in the frequency domain. SIP is widely used to detect sulfide minerals, clay minerals, metallic objects, municipal wastes, hydrocarbons, and salinity intrusion. However, SIP is a static method that cannot measure the dynamics of flow and solute/species transport in the subsurface. To capture these dynamics, the data collected with the SIP technique needs to be coupled with fluid flow and reactive-transport models. To our knowledge, currently, there is no simulator in the open-source literature that couples fluid flow, solute transport, and SIP process models to analyze geoelectrical signatures in a large-scale system. A massively parallel simulation framework (PFLOTRAN-SIP) was built to couple SIP data to fluid flow and solute transport processes. This framework built on the PFLOTRAN-E4D simulator that couples PFLOTRAN (a massively parallel multi-physics simulator for subsurface flow and transport) and E4D (a massively parallel geoelectrical simulator), without sacrificing computational performance. PFLOTRAN solves the coupled flow and solute transport process models to estimate solute concentrations, which were used in Archie’s model to compute bulk electrical conductivities at near-zero frequency. These bulk electrical conductivities were modified using the Cole-Cole model to account for frequency dependence. Using the estimated frequency-dependent bulk conductivities, E4D simulated the real and complex electrical potential signals for selected frequencies for SIP. These frequency-dependent bulk conductivities contain information on geochemical changes in the system. The PFLOTRAN-SIP framework was demonstrated through a synthetic tracer-transport model simulating tracer concentration and electrical impedances for four frequencies. Later, SIP inversion estimated bulk electrical conductivities by matching electrical impedances for each specified frequency. The estimated bulk electrical conductivities were consistent with the simulated tracer concentrations from the PFLOTRAN-SIP forward model. This framework is useful for practitioners of environmental hydrogeophysics and biogeophysics to monitor chemical, nuclear, and tracer transport sites as well as to detect sulphide minerals, metallic objects, municipal wastes, hydrocarbons, and salinity intrusion.

Keywords: Electrical resistivity tomography, inversion spectral induced polarization (SIP), subsurface geoelectrical signatures, hydrogeophysics, multi-physics, porous media flow.

1. Introduction

Engineered subsurface systems are dynamic due to natural and anthropogenic activities that alter porosity, permeability, fluid saturation, and geochemical properties over time [national2000research]. Geophysical techniques such as seismic (deep or near-surface seismic) and potential-based methods (electromagnetic, magnetic, electrical resistivity tomography (ERT), spectral induced polarization (SIP)) can characterize changes in the subsurface [Revil2011, kearey2013, snieder2007advanced]. Among these geophysical methods, ERT and SIP map the distribution of bulk electrical conductivity (i. e., the reciprocal of resistivity) that arises due to changes in subsurface fluid flow, temperature, deformation, and reactive transport [carpenter2017stress, kaselow2004stress, gresse2017, Robinson2015, Byrdina2013]. Because structural, topological, and geochemical properties (e.g., pore structures, fracture networks, electron donor, etc.) influence bulk electrical conductivity [Revil2011, snieder2007advanced], ERT and SIP have a wide application in environmental and energy industries to characterize subsurface interactions. Hence, coupling ERT and/or SIP process models to flow and reactive-transport process models can enhance the interrogation of engineered subsurface systems.

ERT’s data-collection component measures the electric potentials resulting from an applied direct current (DC), while the data-processing component inverts these measured potentials to map the spatial distribution of bulk electrical conductivity [Byrdina2013, Revil2011, Revil2004]. Because ERT injects DC (which is at near zero frequency), it cannot interrogate the polarization features of geologic materials, heavy metals, and induced polarization minerals (e. g., clay minerals, hydrothermal-alteration products, pyrite, finely disseminated sulfide minerals, etc.) [yanetal2014, he2005hydrocarbon, Revil2011]. However, by injecting oscillatory/alternating currents (AC), the induced polarization (IP) method can measure “chargeability” in the time domain or “phase shift” in the frequency domain, which is the phase angle (alternatively phase lag) between the applied current and induced voltage of polarized geologic materials [sumner2012principles, vaulet2011]. The IP method measures the energy storage capacity of certain minerals and can be used to detect hydrocarbons [luo1988], contaminant plumes [olhoeft1986direct, morgan1999induced, vanhala1997mapping], municipal waste, green waste (agricultural and biodegradable wastes) [landfill], sulfide minerals [yanetal2014, butler2005], and hydrothermal products [yanetal2014, butler2005]. IP is a single- or double-frequency method that generally fails to distinguish between a true IP response (the response from polarized geologic materials) and noise (the response from the surroundings or electromagnetic interference) [butler2005, luo1988]. The SIP method extends IP to measure geoelectrical signatures across a wide range of user-specified frequencies to facilitate comprehensive data collection to identify the sources of a true IP response. Also, the SIP method has the potential to characterize subsurface structures and processes [luo1988].

SIP is a comprehensive but static method to extract subsurface polarization signals. It cannot measure evolving dynamics of subsurface processes (e. g., flow, deformation, geochemical reactive-transport, etc.) [JOHNSON2017]. Moreover, it cannot detect all contaminants or chemical reactions occurring in a system [vaulet2011]. However, these shortcomings can be overcome by coupling the SIP method with fluid flow and reactive-transport models. To summarize the coupling, the process models related to fluid flow and reactive-transport simulate the spatio-temporal distribution of conducting fluids in rocks, fluid content in the pores, and fluid chemistry. The concentration of chemical species is transferred to the SIP method, which simulates the electrical potential signals due to polarization. The electrical conductivity in the SIP method contains information on the spatial distribution of conducting fluids and fluid chemistry at different times. In addition, the SIP method inverts for frequency-dependent electrical conductivity based on measured/simulated-electrical impedance and phase-shift data. This coupling facilitates detection, extraction, and understanding of the evolution of hydrogeophysical and biogeophysical signatures at both lab and field scales [KEMNA2012, Atekwana2009, mellage2018linking, SLATER2002].

Software to model and invert geoelectrical data include: Res2Dinv [loke1996rapid, loke2003comparison, res3dinv], Aarhusinv [Aarhusinv2013], BERT [rucker2006three, gunther2006three], EarthImager3D [agi2008], E4D [johnsonetal2010], pyGIMLi [Ruecker2017], and ZondRes3D [zond3res]. The preceding software packages can also image frequency-dependent electrical conductivity of the subsurface but cannot capture dynamic subsurface processes. Moreover, they are not coupled with flow and reactive-transport models, which decreases image quality when interrogating/imaging fluid flow in the subsurface. To overcome these problems, Johnson et al. [JOHNSON2017] developed a massively parallel simulator called PFLOTRAN-E4D, which couples PFLOTRAN [Hammondetal2014], an open-source, state-of-the-art, massively parallel subsurface flow and reactive-transport code to an open-source ERT model in E4D, which is a massively parallel finite element code for simulating and inverting geoelectrical data. However, the PFLOTRAN-E4D framework does not account for induced polarization. To capture dynamics of subsurface processes and the true sources of induced polarization, a computationally efficient framework is needed to couple fluid flow and solute transport with the SIP process model. This work extended the capabilities of PFLOTRAN-E4D to include SIP in a framework called PFLOTRAN-SIP. Here, PFLOTRAN-SIP capabilities are demonstrated with a representative tracer transport model to provide more information than the ERT model, if the medium has polarization properties.

1.1. Main contributions and outline of the paper

In this work, we built a coupled framework PFLOTRAN-SIP to simulate flow, solute transport, and SIP processes for a large-scale system. This framework adds to the recently developed PFLOTRAN-E4D code and takes advantage of the PFLORAN-E4D massively parallel capabilities. The paper is organized as follows: Sec. 1 discusses the limitations of ERT along with the advantages of the SIP method. The importance of coupling between fluid flow, reactive-transport, and SIP process models is discussed. Discussion on the state-of-the-art simulators are also provided. Sec. 2 introduces the PFLOTRAN-SIP framework. A methodology for coupling PFLOTRAN and E4D simulators is described. Process models related to SIP, fluid flow, and solute transport in PFLOTRAN and E4D

simulators are presented. Also, this section discusses the mesh interpolation to transfer state variables between

PFLOTRAN and E4D meshes. Sec. 3 describes a reservoir-scale model setup to perform high-fidelity numerical simulations. It describes the model domain, related meshes, initial conditions, boundary conditions, and various input parameters related to PFLOTRAN-SIP framework. Additionally, this section provides details on the inversion procedure for electrical conductivity at different frequencies. This exercise was performed with the intention of comparing the simulated conductivity/tracer distribution from PFLOTRAN-SIP framework and inverted electrical conductivity from SIP inversion. Sec. 4 explains simulated electrical potentials for a measurement and compares the true/simulated and estimated conductivities. Also, it shows that estimated frequency-dependent conductivities from the SIP inversion are consistent with the tracer concentration/simulated conductivity from the PFLOTRAN-SIP framework. Sec. 5 provides a discussion on the numerical results, limitations of the proposed framework, and how geoscientists can use the PFLOTRAN-SIP tool for their applications. Finally, conclusions are drawn in Sec. 6.

2. Pflotran-Sip: Process Models and Coupling Framework

The PFLOTRAN-SIP framework couples flow and reactive-transport process models in PFLOTRAN [hammond2012pflotran, Hammondetal2014, pflotran-web-page, pflotran-user-ref] with SIP process models in E4D [e4dinweb, johnsonetal2010, e4dmanual] to characterize fluid-driven electrical impedance signatures across multiple frequencies. At each time-step in the framework, the simulation outputs from PFLOTRAN (fluid saturation, tracer concentration, porosity, tortuosity etc.) were supplied to Archie’s model [archie1942] to calculate fluid-dependent bulk electrical conductivities for E4D simulations. The estimated bulk electrical conductivities was decomposed into real and imaginary components for each user-defined frequency using the Cole-Cole model [colecole1941, tarasov2013use], which is an empirical description of frequency-dependent behaviour of bulk electrical conductivities [tarasov2013use]. The above mentioned processes were repeated until the entire transient simulation was completed.

2.1. E4D Process Model

E4D is an open-source, massively parallel, finite element code for simulating and inverting three-dimensional time-lapse ERT and SIP data [e4dinweb, e4dmanual, johnsonetal2010, johnson20173]. The process models in E4D for ERT and SIP assume that displacement currents are negligible and current density can be described by Ohm’s constitutive model [johnson20173]. These assumptions result in a Poisson equation relating induced current to the potential field that determines the electrical potential field in the domain:



[m] is the position vector,

[S ] is the effective electrical conductivities at position , [] is the current injected at position , [] is the electrical potential at , and is the Dirac delta function.

Equation (2.1) models the DC effect, which is required in ERT forward/inverse modeling; however, it does not account for induced polarization under alternating currents. Induced polarization under alternating current results in a secondary potential that needs to be accounted for in the SIP forward/inverse modeling. This requires modification of the preceding equation to solve for the total electrical potential field under IP effects:


where [] is the total electrical potential field, which includes IP effects from a polarized material with chargeability distribution [milliradians] [seigel1959]. The secondary potential resulting from the IP effect is [oldenburg1994]:


and the apparent chargeability [–] is [seigel1959]:


The secondary potential and apparent chargeability can be computed by solving Eqs. (2.1) and (2.2). These potentials , , and are time-domain signatures of induced polarization. Eq. (2.3) is in the time domain and is transformed into the frequency domain in E4D


where [Hz] is the frequency. [S/m] and [V] are the frequency-dependent electrical conductivities and electrical potential, respectively. consists of real and complex electrical potentials corresponding to induced polarization. Zero potentials are enforced on the boundary of the domain [johnson20173, Section 3] to solve Eq. (2.5).

E4D simulates four-electrode configurations (e.g., Wenner array, dipole-dipole array) [johnsonetal2010]. Current is injected from source to sink electrodes while measurements are recorded between the other two electrodes [johnsonetal2010, kearey2013, Robinson2015]. For ERT, the measured response is the potential difference (Voltage) between the two electrodes while SIP also includes the phase shift (radians). Based on the user-defined survey design, E4D simulates hundreds to tens of thousands of ERT/SIP measurements to compute the electrical potential distribution or to invert the lab/field/simulated data to identify the best fitted bulk electrical conductivities distribution throughout the model domain. As the governing equations are linear in and , E4D solves Eq. (2.5) by superimposing pole solutions with different current sources that makes ERT or SIP forward modeling highly scalable [johnsonetal2010, johnson20173].

E4D solves the ERT and SIP process models in the frequency domain using a low-order finite element method (FEM). The output of the FEM solution for the ERT process model is electrical potential throughout the domain, which is real valued and frequency independent. Because the SIP process model is frequency dependent, the corresponding output of the FEM solution has both real and imaginary components of electrical potential. The complex-valued electrical potential or equivalently the phase-shift distribution in the model domain provides new information on IP in the subsurface. This new information (i.e., phase shift) is not captured by ERT.

The finite element methodology [johnsonetal2010] in E4D is based on the standard Galerkin weak formulation [hughes2012finite] on a unstructured low-order tetrahedral finite element mesh [tetgen2015]. E4D iteratively computes the total electrical potential field due to IP effects [johnson20173, Section 3]. Equations for computing the real and imaginary components of the complex-valued electrical potential are decoupled and the finite element analysis is performed in the real number domain. First, E4D solves for the real component without considering IP effects. Second, the current source for the imaginary component is computed from the real component. Third, the imaginary component of the total electrical potential is calculated based on this computed current source. Fourth, the secondary current source arising from the imaginary component is computed. This secondary current source takes IP effects into account. Later, the real component is calculated based on this secondary current source. These steps are repeated until a convergence criterion is satisfied; change in the Euclidean norm of the real potential field between subsequent iterations is less than a user-defined tolerance.

2.2. PFLOTRAN Process Models


solves a system of nonlinear partial differential equations describing multiphase, multicomponent, and multiscale reactive flow and transport using the finite volume method (FVM)

[hammond2012pflotran, Hammondetal2014, lichtner2015pflotran]. In this paper, we restrict to solving the governing equations for single-phase fluid flow and solute transport towards predicting the spatio-temporal distribution of solute concentrations. The governing mass conservation equation for single-phase variably saturated flow is:


where [kg ] is the fluid density, [–] is the porosity, [–] is the saturation, t [] is time, is the Darcy flux, and is the volumetric source/sink term. Darcy flux is:


where [m2] is the intrinsic permeability, [–] is the relative permeability, is the dynamic viscosity, [Pa] is pressure, [m s–1] is gravity, and is the vertical component of the position vector . The source/sink term is:


where is the mass flow rate, is the formula weight of water, and denotes the location of the source/sink. The governing equation for tracer transport:


where [molality] is the solute concentration, [m2 s–1] is the diffusion/dispersion coefficient, [–] is tortuosity (related to the path length of the fluid flow), and [molality ] is the solute source/sink term. PFLOTRAN admits different types of boundary conditions (e.g., Dirichlet, Neumann, Robin) when solving Eqs. (2.6)–(2.9).

Figure 1. PFLOTRAN Process Modeling: Peer and child process model class of PFLOTRAN (redrawn from reference [JOHNSON2017]).

The coupled governing equations given by Eqs. (2.6)–(2.9) are solved using a two-point flux FVM for spatial discretization and a fully implicit backward Euler method for time discretization. The resulting nonlinear algebraic equations are solved using a Newton–Krylov solver [hammond2012pflotran, balay2017petsc]. The workflow of PFLOTRAN is divided into a tree structure (see Fig. 1), which allows hierarchical coupling of flow and solute-transport process models to numerical methods (e.g., FVM, backward Euler, Newton-Raphson, Newton-Krylov, logarithmic formulation etc.) [hammond2012pflotran, JOHNSON2017]. A process model tree as shown in Fig. 1 has a master-process model and two pointers – child process model and peer process model . Here, the flow model is master process model while and are for the solute transport and E4D/SIP models, respectively. The time step for solving the flow model may be different from solute-transport model. Transfer of information between (e.g., flow) and (e.g., solute transport) takes place before and after each of ’s time steps. Synchronization of and (e.g., ERT or SIP) occurs at specified times. Execution starts with the master-process model , which can take as many adaptive time steps as needed to reach the synchronization point. ’s time step may be the same as ’s or smaller. and march forward in time according to their own time steps to reach the synchronization point. When , , and all reach the synchronization point, key variables and parameters (e.g., saturation, solute concentration, porosity etc.) are updated between and .

2.3. PFLOTRAN-SIP Coupling

The coupling involves six steps: (1) the flow process model in PFLOTRAN is called to calculate fluid pressure and velocity; (2) the solute transport process model is solved using simulation outputs (e.g., Darcy flux and fluid saturation) from the flow process model to obtain solute concentrations; (3) solute concentrations in each PFLOTRAN mesh cell are used to calculate the DC electrical conductivities for E4D based on Archie’s model; (4) the Cole-Cole model calculates the frequency-dependent electrical conductivities; (5) real and complex electrical conductivities are interpolated onto the E4D mesh; and (6) the SIP process model in E4D is called to solve the forward problem to calculate electrical impedance and phase shift.

PFLOTRAN and E4D use message passing interface calls for inter-process communication. Based on user specification, PFLOTRAN divides the computing resources between PFLOTRAN and E4D at the initial step. PFLOTRAN and E4D read their corresponding input files and complete pre-simulation steps. These include setup of the flow model, the solute transport model, the SIP model, and mesh interpolation matrix. Mesh interpolation is needed for two reasons: (1) the meshes of PFLOTRAN and E4D are different and (2) the solution procedure of PFLOTRAN is based on the FVM while E4D’s solution procedure is based on the FVM. As a result, the state variables (e.g., solute concentration, fluid saturation) computed at the cell center by PFLOTRAN need to be accurately transferred from the PFLOTRAN mesh to the E4D mesh to calculate electrical conductivities. Generation of the mesh interpolation matrix is described in Sec. 2.5. Algorithm 1 and Fig. 2 summarize the coupling of PFLOTRAN and E4D SIP models.

Figure 2. Coupling PFLOTRAN and SIP Process Models: Steps involved in coupling fluid flow, solute transport, and SIP process models in the PFLOTRAN-SIP framework; further details are available [johnson20173, lichtner2015pflotran].

2.4. Petrophysical Transformation

To simulate SIP signals during fluid flow and solute transport, a mathematical relationship linking fluid flow state variables and bulk electrical conductivities is required. Archie’s model [archie1942, glover2010generalized, shah2005generalized] is a petrophysical transformation relating state variables simulated by PFLOTRAN to bulk electrical conductivities:


where [S m-1] is the bulk electrical conductivity at near zero frequency (), [–] is the cementation exponent (an empirical constant, e.g., ranges from 1.8 to 2.0 for sandstone), [–] is the solute concentration simulated by PFLOTRAN, [–] is the saturation exponent (also an empirical constant that is usually close to 2.0), and [S ] is the fluid electrical conductivity. Note that Archie’s model assumes that the matrix of the porous medium is non-conductive. However, for porous materials with clay minerals, this assumption is not valid due to the structure of the clay and its cation exchange capacity. For such materials, a modified version of Archie’s model (e.g., Waxman–Smits model [vinegar1984induced]) can be used instead of Eq. (2.10).

The saturation exponent is related to the wettability of the rock and reflects the presence of non-conductive fluids (hydrocarbons) in the pore space. The cementation exponent quantifies how much the pore structure and the pore network change the electrical resistivity (or decreases the electrical conductivity) that is usually assumed independent of temperature. Additionally, increasing permeability decreases the cementation exponent. Common values for are between 1.0 and 4.1. For unconsolidated sands, and for consolidated sandstones,

. In carbonate rocks, the cementation exponent shows higher variance due to complex pore structures and networks with values of

between 1.7 and 4.1 [glover2009cementation].

To account for frequency dependence, Eq. (2.10) was modified using the Cole-Cole model [colecole1941, cole1942dispersion, dias2000developments, tarasov2013use, revil2014spectral]:


where is a complex number such that , [–] is a shape parameter (an empirical constant), and [s] is the characteristic relaxation time constant (required time for imaginary electrical component to reach equilibrium after perturbation) related to characteristic pore or grain size.

2.5. Mesh Interpolation

Once the frequency-dependent real and imaginary components of bulk electrical conductivities were calculated on PFLOTRAN mesh, they were interpolated onto the E4D mesh. The conductivity at any intermediate point in a PFLOTRAN mesh cell was approximated using tri-linear interpolation. These approximated values were computed using values at the PFLOTRAN cell centers surrounding the point (see Fig. 3) using [JOHNSON2017]:


where [m3] is the volume of the th element of E4D mesh, [S ] is the bulk electrical conductivity in the PFLOTRAN mesh element, [–] is the number of subdivisions by which th E4D element is divided for integral approximation, and [S ] is the bulk conductivity within the th element and subdivision . The value of is:


where [–] is the total number of mesh elements in the PFLOTRAN mesh, is the bulk conductivity of the th element in the PFLOTRAN mesh, and is the linearly interpolated weight for to determine the value of sub-element in the th E4D mesh element.

The preceding equations interpolated data onto the E4D mesh based on the computed values of in the PFLOTRAN mesh:


The matrix form of the preceding equation is:


where is a vector of length (total number of PFLOTRAN mesh elements) that contains the bulk conductivities computed on the PFLOTRAN mesh, is a vector of length that contains the bulk conductivities values in the E4D mesh, is the total number of E4D mesh elements, and () is a sparse matrix containing weights . A parallel approach to compute is available [JOHNSON2017]. Please, see the work of Johnson et. al., 2017 [johnson20173] for more details.

Note that there are various mesh transfer/interpolation algorithms (e.g., seminal works by Michael Heath, Xiangmin Jiao, Charbel Farhat, Patrick Farrell, Pavel Bochev, Mikhail Shashkov, etc. [jiao2004common, jiao2004overlaying1, jiao2004overlaying2, farhat1998load, farrell2011conservative, farrell2009conservative, bochev2011formulation, bochev2013fast]) that provide robust and reliable transfer (e.g., preserving mass, chemical species, and energy) of state variables across non-matching meshes. Extension of this proposed approach using conservative interpolations will be considered in future efforts.

Figure 3. Mesh Interpolation from PFLOTRAN to E4D: Schematic of interpolation of state variables (e.g., solute concentration) on the PFLOTRAN mesh (cube) on to E4D mesh (tetrahedron), redrawn from [JOHNSON2017].
1:  INPUT: Initial and boundary conditions for fluid flow and solute transport models in PFLOTRAN, fluid density, porosity, saturation, volumetric source/sink with its location, intrinsic and relative permeabilities, dynamic viscosity, mass flow rate, diffusion/dispersion coefficients, tortuosity, solute source/sink with its location, Archie’s and Cole-Cole model parameters, total simulation time, time-step for PFLOTRAN, interrogation frequencies, electrode locations and measurement configuration, number of processors for PFLOTRAN and E4D, and meshes for PFLOTRAN and E4D.
2:  Solve Eqs. (2.6)–(2.8) for fluid pressure, fluid saturation, and fluid velocity.
3:  Solve Eq. (2.9) to calculate the spatio-temporal distribution of solute concentration.
4:  Transfer solute concentration from PFLOTRAN to the E4D master processor to perform SIP simulations at specific times.
5:  Receive numerical model setup information from PFLOTRAN input files to perform mesh interpolation for SIP simulations.
6:  Broadcast run information and distribute mesh assignments to E4D slave processors.
7:  Calculate the mesh interpolation matrix Eq. (2.14) to interpolate PFLOTRAN simulation outputs (e.g., solute concentrations) onto the E4D mesh for SIP simulations.
8:  Calculate electrical conductivities using Archie’s model Eq. (2.10).
9:  Calculate frequency-dependent electrical conductivities using the Cole-Cole model Eq. (2.11).
10:  Pass real and complex conductivities calculated at different frequencies to the E4D master processor to perform SIP simulations.
11:  Broadcast real and complex conductivities to E4D slave processors to compute pole solutions for electrode configurations.
12:  Solve Eq. (2.5) to compute real and complex electrical potential at different frequencies and solute concentrations at specified times.
Algorithm 1 Overview of the proposed PFLOTRAN-SIP framework for simulating electrical impedance data

3. Numerical Model Setup

3.1. PFLOTRAN Model Setup

The domain was  m3 and consisted of three layers as shown in Fig. 4. The PFLOTRAN mesh had a total of 125,000 finite volume cells. The upper layer was  m3 and extended from = 0 to  m as a highly conductive material with a permeability of  m2. Fluid was water and these rock properties (e.g., permeability, porosity, diffusion coefficient etc.) are representative of sandstone. The middle layer was less permeable (permeability was  m2), with size  m3 extending from = to  m. This permeability is representative of rocks such as shale or granite. The low-permeability layer, however, included a small-volume, highly permeable ( m2) material between and  m, and  m, and and  m. The bottom layer was also a highly conductive region with a permeability of  m2 and dimensions of  m3.

A solute (conservative tracer) at 10 mol was placed below the low permeable zone as shown in Fig. 4 as the purple  m3 block. The initial and boundary conditions for the model included: pressure of 1 atm at the top with a hydrostatic pressure gradient from top to bottom. The right face () was assigned 2 atm to drive fluid from right to left. For solute transport, the boundary conditions were zero-concentration Dirichlet inflow at the right face and zero diffusive gradient outflow at the left face that allowed only advective outflow. The remaining faces were specified as zero-solute flux boundaries.

For low- and high-permeability zones, tortuosity was 1.0 while porosities were 0.3 and 0.25, respectively. Solute diffusivity was  m2 s–1. The Newton solver (20-iteration maximum) was applied for flow and solute transport. For the flow solver, relative and absolute tolerances were with a relative update tolerance of while for solute transport solver, relative and absolute tolerances were with a relative update tolerance of . The simulation was run for one year with an initial time step of  year, which was allowed to accelerate by a factor of 8.

Figure 4. PFLOTRAN model domain: Schematic of the initial boundary value problem.

3.2. SIP Model Setup

Although the domain dimensions for SIP simulations were identical to the PFLOTRAN simulation, there was only a single layer. The corresponding E4D mesh for the simulation had 86,780 nodes and 609,562 tetrahedral elements. Zero potentials were enforced on the boundaries. A total of 80 point electrodes were placed in the domain, all located at  m arranged in 5 lines along the -axis, with each line comprising 16 electrodes. The electrode coordinates started at and ended at with a  m separation between lines. Electrode measurement configurations included a combination of Wenner and dipole-dipole arrays. A current of 1A was injected and received at a pair of electrodes and the potential difference was measured at another pair of electrodes. There are various advantages of injecting and receiving the current through a pair of electrodes. For example, such a measurement system can eliminate any inaccuracies caused by the injecting circuit impedance (the contact impedance between the probe and the medium, which can be high). Using the prescribed measurement configuration, a total of 1,062 simulated measurements were collected to capture electrical impedance and phase shift.

The electrical conductivity of the fluid at  Hz was  S/m. The cementation and saturation exponents were 0.564 and 0.576, respectively. The characteristics-relaxation time was 0.061 s, all representative of sandstone [titov2002theoretical]. SIP analysis was performed for five different frequencies: 0.1, 1.0, 10, 100, and 1,000 Hz. Forward model simulations were performed using 61 processors, where 20 processors were assigned for PFLOTRAN and 41 for E4D. Out of those 41 processors, 40 performed SIP simulations for different measurement configurations and the remaining processor gathered the simulated data.

3.3. SIP Inversion of electrical conductivity

For verification, E4D’s inversion module was used to estimate frequency-dependent electrical conductivity based on the simulated electrical impedance and phase-shift data. This estimated conductivity was compared with the simulated conductivity generated by the PFLOTRAN-SIP framework. The employed inversion process was blind (i.e., we did not provide prior constraints on the conductivity). This can be improved by providing detailed conductivity information to E4D’s inversion module. The SIP inversion employs an unstructured mesh, which consisted of 51,124 nodes with 316,183 mesh elements. Low-order mesh elements were generated to make inversion process simple and computationally efficient because high-order mesh elements did not improve the outcome [e4dmanual]. The simulated measurements by PFLOTRAN-SIP were the data supplied to the inversion process as observations.

The E4D inversion is based on minimizing the following objective function to estimate the frequency-dependent electrical conductivity distribution, :


where is an operator that provides a scalar measure of the misfit between observed and simulated data (e.g., electrical impedance and phase shift) based on the user-specified norm (e.g., Euclidean norm), is another operator that provides the scalar measure of the difference between  [S ] and constraints placed upon the structure of [S ], is the regularization parameter, is the data-weighting matrix, and is the model-weighting matrix. and are the estimated and reference frequency-dependent electrical conductivities. The user specified bounds on the frequency-dependent conductivity in each mesh cell were 0.00001 and 1.0. The and were the observed and simulated data, respectively. Eq. (3.1) is solved using the iteratively reweighted least square method [scales1988robust]. Further details on the parallel inverse modeling algorithm and its implementation in E4D are available [johnsonetal2010].

The value was 100 at the beginning of the inversion and decreased as the non-linear iteration progressed. Before was reduced, the minimum fractional decrease in the objective function, , between iterations had to be less than 0.25 upon which was reduced to 0.5. The convergence of the SIP inversion procedure was based on the value of the current iteration after data culling, computed as:


where the data residual is the difference between observed and estimated values divided by the standard deviation for that measurement.

is the total number of survey measurements and is the number of measurements selected from the total number of measurements during the current iteration. SIP inversion converged after 48 iterations when reached 60 using the example model from Sec. 4.

4. Results

The one-year PFLOTRAN-SIP model simulations were completed in two minutes. The computation was performed on 61 Intel® Xeon® CPU E5-2695 V4 @ 2.1GHz processors. Fig. 5 shows the tracer concentration at the end of simulation. In one year, the pressure gradient drove the tracer about 100 m from its initial location in the -direction and also moved it upward about 20 m.

Figure 5. PFLOTRAN Simulation: Spatial distribution of tracer concentrations after one year.

The SIP module in PFLOTRAN-E4D simulated real and complex electrical impedances at 0.1, 1, 10, 100, and 1,000 Hz. Because of minimal differences between 1 and 10 Hz, only results for 0.1, 10, 100, and 1,000 Hz are discussed. This indicated that some frequencies may be redundant because they yielded similar impedances. Sensitivity analysis can be performed to identify redundant frequencies; however, this was beyond the scope of this paper. Fig. 6 shows the real and complex potentials due to changes in tracer concentration for the various frequencies. Also, this figure provides information on the change in electrical potential at different frequencies for a single measurement configuration, indicating the maximum tracer concentration. This 80-node measurement configuration was selected because tracer concentration was most visible. The response clearly shows the polarization feature of the tracer. The gradient of the real electrical potential was high near  m (top row of Fig. 6) because the tracer concentration was maximum. From Fig. 6, it is evident that the real potential response for 0.1 Hz is different from the responses at 10, 100, and 1,000 Hz The root-mean-square error (RMSE) between these responses was approximately 15% of the maximum real potential value indicating that frequency has an impact on the real potential distribution.

Figure 6. Slices of simulated real (top) and complex (bottom) electrical potentials/impedances at = 250 m for a single electrode-measurement configuration after one year.

The bottom row of Fig. 6 shows complex electrical potential responses where the polarity was switched (colors interchanged). Unlike the real electrical potential, each complex electrical potential was notably different suggesting that it was frequency dependent. The corresponding RMSE between responses was 85% of the maximum complex potential value. Such high variation was expected as the complex electrical potential depends on frequency, chargeability, and relaxation time although the last two were constant in this study. Because the response of the complex potential was clearly visible in the simulation, this indicated that the PFLOTRAN-SIP framework can effectively simulate polarized geologic materials.

Figure 7. Simulated and estimated frequency-dependent electrical conductivities at = 250 m after one year (a)-(d) True-electrical conductivities from the PFLOTRAN-SIP framework, (e)-(h) estimated bulk-real conductivities from SIP inversion, and (i)-(l) estimated bulk complex electrical conductivities from SIP inversion.

Fig. 7 shows the simulated and estimated frequency-dependent electrical conductivities using the PFLOTRAN-SIP framework with the SIP inversion module in E4D. The true (PFLOTRAN simulated) and estimated (inversion of survey data) real electrical conductivities are plotted in Fig. 7(a)-(d) and Fig. 7(e)-(h), respectively. SIP inversion was performed using the simulated electrical impedance and phase-shift data obtained from PFLOTRAN-SIP model runs after one year. The computational time required to perform SIP inversion was approximately two hours on 41 Intel® Xeon®d CPU E5-2695 V4 @ 2.1GHz processors. Estimated electrical conductivities showed high contrast around the high tracer distribution/simulated conductivities although they were more diffuse than the true/simulated distribution (Fig. 7(a-h)). However, estimated conductivities at 1,000 Hz were more accurate than frequencies <1,000 Hz with the same constraints. Later, real conductivity values were used in Eq. 2.11 to provide initial guesses for complex conductivities for SIP inversion. Estimated complex electrical conductivity distributions are shown in Fig. 7(i)-(l). Similar to estimated real conductivities, complex conductivities computed from SIP inversion were also diffuse. The inversion process can be improved by providing prior information and structural constraints on electrical conductivities. However, both estimated conductivities were generally consistent with the tracer distribution, which showed that the SIP inversion module can simulate electrical impedance and phase-shift data. To summarize, SIP provides a major benefit, which ERT lacks. SIP provides greater information content than ERT. This is because the SIP survey yields multiple datasets at different frequencies, which help to overcome false positives. For example, from Fig. 7 it is evident that the SIP inversion analyses at different frequencies are indicating the same tracer region, (not a false positive). With an ERT survey, it may be difficult to delineate a false positive from a true positive because ERT generates a single dataset.

Fig. 8(a)-(c) show simulated outputs of tracer concentrations, real potentials, and complex potentials for the 80-electrode measurement configuration at frequencies of 0.1, 10, 100, and 1,000Hz. The location of maximum tracer concentration was around  m. The locations of current and potential measurement electrodes were at ( and  m,  m, and  m) (Fig. 8(a)). Note that electrodes were not placed at the location of maximum concentration because of the lack of a prior knowledge of tracer fate (which is the case in real-world applications). Because the measurement electrodes were offset from the maximum tracer concentration, this resulted in an offset of peaks between tracer distribution and geoelectrical signals. However, the measured potentials provided meaningful information on the bounds of the tracer distribution as well as revealing the significance of higher frequencies obtained from a combination of the electrical impedance and phase shift.

Fig. 8 (b) and (c) showed that the absolute real potential and complex potential decreased as frequency increased. As noted in Sec. 2.1 and in Reference [johnson20173], for SIP simulations, E4D first solves the real potential . That is, , where is the real component of and the real potential, is inversely proportional to . Also, increases as increases; hence the absolute value of the real potential distribution (as shown in Fig. 8 (b)) decreases as increases. After is evaluated, E4D computed the complex potential by solving where is the imaginary part of and is the imaginary potential. Thus, is proportional to . Also, decreases as increases; hence the absolute value of the complex potential distribution (as shown in Fig. 8 (c)) decreased as increased.

Fig. 8 (d) shows phase-shift data distribution along the same line for tracer distribution, real and imaginary potential distribution. Mathematically, the phase shift is the inverse tangent of the ratio between complex and real potential responses. Physically, it is the shift between voltage and current signals that is largely governed by the polarization characteristics of the subsurface. In this study, phase shift leveraged signals from both real and complex potential responses to improve interpretation of survey data. From Fig. 8, there was a change in phase shift where tracer transport was predominant. Moreover, the 1,000 Hz frequency bounded the tracer zone better than lower frequencies that cannot be distinguished with ERT. This change in phase helped constrain the polarized region or bound the interface between tracer-filled and tracer-free fluids. Obtaining the region of interest using these constraints, further geoelectrical interrogation can be performed with this volume. Hence, through phase-shift signatures across multiple frequencies, the PFLOTRAN-SIP framework will facilitate identification of polarized or geochemically altered zones.

Figure 8. Distribution of (a) tracer concentration where I and V represents current and potential electrodes, respectively, (b) real potential, (c) complex potential, and (d) phase shift along the -axis at  m and  m.

5. Discussion

IP arises from solute transport and accumulation of ions/electrons in polarized materials (e.g., materials that have different grain types, colloids, biological materials, phase-separated polymers, blends, and crystalline minerals, etc.) that are subject to an external electric field. A total of five different mechanisms govern the IP phenomena at frequencies  MHz (which are of interest in this paper). These include: (1) Maxwell-Wagner polarization, which is a phenomenon that occurs at high frequency [alvarez1973complex, chelidze1999electrical, lesmes2001dielectric, chen2006geometrical]; (2) polarization at the inner part of the interface between minerals and water [de1992generalized, leroy2009mechanistic, vaudelet2011induced, revil2012spectral]; (3) polarization at the outer part of the interface between minerals and water [dukhin1974dielectric, de1992generalized]; (4) membrane polarization for multi-phase systems [marshall1959induced, vinegar1984induced, titov2002theoretical]; and (5) electrode polarization observed in the presence of disseminated conductive minerals such as sulfide minerals and pyrite [wong1979electrochemical, merriam2007induced, seigel2007early].

Our PFLOTRAN-SIP simulations were geared toward IP mechanisms that fall in (1), (4), and (5). To simulate the mechanisms mentioned in (2) and (3), Eq. (2.11) must be replaced with conductivity models that account for interface polarization with consideration of effective pore size, electrical formation factor, distribution of relaxation times, and sorption mechanisms [revilandflorsch2010, vaulet2011]. Note that the Cole-Cole model given by Eq. (2.11) neglects the effects of polarization at interfaces or sorption onto mineral surfaces. The PFLOTRAN-SIP framework can easily account for such modifications in conductivity, but this is a future endeavor. Additionally, the proposed framework considers conductivity models, which provide better information on porosity, permeability, tortuosity, and fluid phase distribution.

6. Conclusions

This work developed a computational framework to couple PFLOTRAN and E4D to model electrical-impedance and phase-shift data for SIP due to changes in subsurface processes. PFLOTRAN and E4D are massively parallel codes that simulate processes related to fluid flow, reactive transport, and SIP. A mathematical relationship based on Archie’s and the Cole-Cole models linked flow and solute-transport state variables at various frequencies. A reservoir-scale tracer transport model demonstrated the proposed PFLOTRAN-SIP framework where fluid flow and tracer concentration evolution were simulated over one year. Then, we simulated electrical potentials for various electrode configurations at different frequencies. These simulations showed that contrast in real potential was minimal even as the frequency varied. However, there was a significant change in contrast of complex potentials across frequencies. Phase shift (combination of real and complex potentials) helped identify the region where tracer concentration was high. This analysis showed that SIP has two major advantages over ERT. First, SIP provided frequency-dependent electrical impedance data. Second, phase-shift signatures obtained from SIP analysis identified and constrained geochemically altered zones. Combining frequency-dependent real potential, complex potential, and phase responses from a SIP survey/simulation paints a more detailed picture of the subsurface with an enhanced ability to detect contaminants/tracers. Moreover, coupling fluid flow, reactive transport, and SIP models can better detect contaminants compared to either the ERT or SIP method alone. For instance, through our numerical example, solute transport simulations provided insight into the tracer distribution. This information was used to customize SIP inversion to estimate frequency-dependent electrical conductivities, which yielded an improved image of tracer concentration at different frequencies. In summary, we developed an open-source computational framework that allows practitioners to better monitor sites with polarized materials. Although this work focused on simulating tracer transport, it could also be applied to detect hydrocarbon flow, changes in the subsurface due to geochemical reactions, sulphide minerals, metallic objects, municipal wastes, and salinity intrusion.

Data and Resources

PFLOTRAN source code can be downloaded at E4D source code can be downloaded at The PFLOTRAN-SIP simulation input files used for this manuscript are available in the public Github repository Additional information regarding the simulation datasets can be obtained from Bulbul Ahmmed (Email: and Maruti Kumar Mudunuru (Email:


This research was funded by the U.S. Department of Energy (DOE) Basic Energy Sciences (BES) and Fossil Energy (FE) programs. MKM, SK, and BA also thank the support from Center for Space and Earth Sciences (CSES) Emerging Ideas R&D Program. The authors thank Glenn Hammond (Sandia National Laboratories) and Tim Johnson (Pacific Northwest National Laboratory) for the coupled framework PFLOTRAN-E4D upon which PFLOTRAN-SIP is built. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).