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

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.


INTRODUCTION
Engineered subsurface systems are dynamic due to natural and anthropogenic activities that alter porosity, permeability, fluid saturation, and geochemical properties over time [1]. 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 [2][3][4]. 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 [5][6][7][8][9]. Because structural, topological, and geochemical properties (e.g., pore structures, fracture networks, electron donor, etc.) influence bulk electrical conductivity [2,4], 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 [2,9,10]. 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.) [2,11,12]. 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 [13,14]. The IP method measures the energy storage capacity of certain minerals and can be used to detect hydrocarbons [15], contaminant plumes [16][17][18], municipal waste, green waste (agricultural and biodegradable wastes) [19], sulfide minerals [11,20], and hydrothermal products [11,20]. 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) [15,20]. 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 [15].
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.) [21]. Moreover, it cannot detect all contaminants or chemical reactions occurring in a system [14]. 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 [22][23][24][25].
Software to model and invert geoelectrical data include: Res2Dinv [26][27][28], Aarhusinv [29], BERT [30,31], EarthImager3D [32], E4D [33], pyGIMLi [34], and ZondRes3D [35]. 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. [21] developed a massively parallel simulator called PFLOTRAN-E4D, which couples PFLOTRAN [36], 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.

PFLOTRAN-SIP: PROCESS MODELS AND COUPLING FRAMEWORK
The PFLOTRAN-SIP framework couples flow and reactive-transport process models in PFLOTRAN [36][37][38][39] with SIP process models in E4D [33,40,41] 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 [42] 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 [43,44], which is an empirical description of frequency-dependent behaviour of bulk electrical conductivities [44]. The above mentioned processes were repeated until the entire transient simulation was completed.

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 [33,40,41,45]. 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 [45]. These assumptions result in a Poisson equation relating induced current to the potential field that determines the electrical potential field in the domain: is the position vector, σ [S m −1 ] is the effective electrical conductivities at position x, I [A] is the current injected at position x 0 , Φ σ (x) [V] is the electrical potential at x, 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: is the total electrical potential field, which includes IP effects from a polarized material with chargeability distribution η (r) [milliradians] [46]. The secondary potential resulting from the IP effect is [47]: and the apparent chargeability [-] is [46]: The secondary potential Φ s and apparent chargeability η a can be computed by solving Eqs.  E4D simulates four-electrode configurations (e.g., Wenner array, dipole-dipole array) [33]. Current is injected from source to sink electrodes while measurements are recorded between the other two electrodes [3,8,33]. 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 [33,45].
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 [33] in E4D is based on the standard Galerkin weak formulation [48] on a unstructured low-order tetrahedral finite element mesh [49]. E4D iteratively computes the total electrical potential field due to IP effects [45,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. PFLOTRAN solves a system of nonlinear partial differential equations describing multiphase, multicomponent, and multiscale reactive flow and transport using the finite volume method (FVM) [36,37,50]. 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: is the volumetric source/sink term. Darcy flux is: is the vertical component of the position vector x. The source/sink term Q w is: is the formula weight of water, and x ss [m] denotes the location of the source/sink. The governing equation for tracer transport: 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 [37,51]. 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.) [21,37]. A process model tree as shown in Fig. 1 has a master-process model A and two pointers -child process model B and peer process model C. Here, the flow model is master process model A while B and C 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 A (e.g., flow) and B (e.g., solute transport) takes place before and after each of A's time steps. Synchronization of A and C (e.g., ERT or SIP) occurs at specified times. Execution starts with the master-process model A, which can take as many adaptive time steps as needed to reach the synchronization point. B's time step may be the same as A's or smaller. B and C march forward in time according to their own time steps to reach the synchronization point. When A, B, and C all reach the synchronization point, key variables and parameters (e.g., saturation, solute concentration, porosity etc.) are updated between A and C.  (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.

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 [42,52,53] is a petrophysical transformation relating state variables simulated by PFLOTRAN to bulk electrical conductivities: The saturation exponent β is related to the wettability of the rock and reflects the presence of nonconductive 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, α ≈ 1.3 and for consolidated sandstones, 1.8 < α < 2.0. 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 [55].
To account for frequency dependence, Eq. (2.10) was modified using the Cole-Cole model [43,44,[56][57][58]: where i is a complex number such that i 2 = −1, γ [-] 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.

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 [21]: where V i [m 3 ] is the volume of the i th element of E4D mesh, σ c [S m −1 ] is the bulk electrical conductivity in the PFLOTRAN mesh element, n k [-] is the number of subdivisions by which i th E4D element is divided for integral approximation, and σ c i,k [S m −1 ] is the bulk conductivity within the i th element and subdivision k. The value of σ c i,k is: where n c [-] is the total number of mesh elements in the PFLOTRAN mesh, σ c j is the bulk conductivity of the j th element in the PFLOTRAN mesh, and W i,k,j is the linearly interpolated weight for σ c j to determine the value of sub-element k in the i th E4D mesh element.
The preceding equations interpolated data onto the E4D mesh based on the computed values of σ(x, ω) in the PFLOTRAN mesh: 14) The matrix form of the preceding equation is:  [59][60][61][62][63][64][65][66]) 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.

NUMERICAL MODEL SETUP
3.1. PFLOTRAN Model Setup. The domain was 500 × 500 × 500 m 3 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 500 × 500 × 350 m 3 and extended from z = 0 to −350 m as a highly conductive material with a permeability of 7.38 × 10 −13 m 2 . Fluid was water and these rock properties (e.g., permeability, porosity, Algorithm 1 Overview of the proposed PFLOTRAN-SIP framework for simulating electrical impedance data 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. For low-and high-permeability zones, tortuosity was 1.0 while porosities were 0.3 and 0.25, respectively. Solute diffusivity was 10 −9 m 2 s -1 . The Newton solver (20-iteration maximum) was applied for flow and solute transport. For the flow solver, relative and absolute tolerances were 10 −50 with a relative update tolerance of 10 −60 while for solute transport solver, relative and absolute tolerances were 10 −4 with a relative update tolerance of 10 −60 . The simulation was run for one year with an initial time step of 10 −8 year, which was allowed to accelerate by a factor of 8.  x-axis, with each line comprising 16 electrodes. The electrode coordinates started at (40, 50, −425) and ended at (460, 450, −425) with a 28 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 ω = 0 Hz was 2 × 10 −3 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 [67]. 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.

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 [41]. 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 frequencydependent electrical conductivity distribution, σ est :

1)
where Φ d 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), Φ m is another operator that provides the scalar measure of the difference between σ est [S m −1 ] and constraints placed upon the structure of σ ref [S m −1 ], ζ is the regularization parameter, W d is the data-weighting matrix, and W m is the model-weighting matrix. σ est and σ ref are the estimated and reference frequencydependent electrical conductivities. The user specified bounds on the frequency-dependent conductivity in each mesh cell were 0.00001 and 1.0. The Φ obs and Φ pred were the observed and simulated data, respectively. Eq. (3.1) is solved using the iteratively reweighted least square method [68]. Further details on the parallel inverse modeling algorithm and its implementation in E4D are available [33]. 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 χ 2 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. n d is the total number of survey measurements and n c is the number of measurements selected from the total number of measurements during the current iteration. SIP inversion converged after 48 iterations when χ 2 reached 60 using the example model from Sec. 4.

RESULTS
The one-year PFLOTRAN-SIP model simulations were completed in two minutes. The computation was performed on 61 Intel R Xeon R 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 y-direction and also moved it upward about 20 m. 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 x = 300 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-meansquare 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. 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. 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 R Xeon R 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 x = 300 m. The locations of current and potential measurement electrodes were at (x = 208, 236, 264, and 292 m, y = 250 m, and z = −420 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 [45], for SIP simulations, E4D first solves the real potential Φ r . That is, −div [σ r grad [Φ r ]] = I, where σ r is the real component of σ * (x, ω) and the real potential, Φ r is inversely proportional to σ r . Also, σ r increases as ω increases; hence the absolute value of the real potential distribution (as shown in Fig. 8 (b)) decreases as ω increases. After σ r is evaluated, E4D computed the complex potential by solving div [σ r grad [ where σ c is the imaginary part of σ * (x, ω) and Φ c is the imaginary potential. Thus, Φ c is proportional to σ c . Also, σ c 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.

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 < 1 MHz (which are of interest in this paper). These include: (1) Maxwell-Wagner polarization, which is a phenomenon that occurs at high frequency [69][70][71][72]; (2) polarization at the inner part of the interface between minerals and water [73][74][75][76]; (3) polarization at the outer part of the interface between minerals and water [73,77]; (4) membrane polarization for multi-phase systems [54,67,78]; and (5) electrode polarization observed in the presence of disseminated conductive minerals such as sulfide minerals and pyrite [79][80][81].
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 [14,82]. 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.

CONCLUSIONS
This work developed a computational framework to couple PFLOTRAN and E4D to model electricalimpedance 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. 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 opensource 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.