Method of imaging the electrical conductivity distribution of a subsurface

Abstract

A method of imaging electrical conductivity distribution of a subsurface containing metallic structures with known locations and dimensions is disclosed. Current is injected into the subsurface to measure electrical potentials using multiple sets of electrodes, thus generating electrical resistivity tomography measurements. A numeric code is applied to simulate the measured potentials in the presence of the metallic structures. An inversion code is applied that utilizes the electrical resistivity tomography measurements and the simulated measured potentials to image the subsurface electrical conductivity distribution and remove effects of the subsurface metallic structures with known locations and dimensions.

Claims

1. A method for imaging electrical conductivity distribution within a subsurface containing metallic structures with known locations and dimensions, the method comprising the steps of: a. generating electrical resistivity tomography measurements by injecting electrical currents into the subsurface and measuring electrical resistivity of items within the subsurface using multiple sets of electrodes; b. generating simulated data by applying a numeric code that simulates the measured resistivity of the known metallic structures; c. calculating tomographs from the simulated and measured electrical resistivity data; and d. applying an inversion code that removes the effects of the subsurface metallic structures with known locations and dimensions on the measured data.

2. The method of claim 1 wherein the numeric code is a forward model.

3. The method of claim 2 wherein the applying the forward model to simulate the subsurface potentials comprises a solution to a Poisson equation.

4. The method of claim 3 wherein the Poisson equation includes at least one of the following: integral calculations, analytical calculations, numerical calculations, partial solutions calculations, digital calculations, and an immersed interface boundary condition solution.

5. The method of claim 3 wherein the inversion code is applied so that the simulated electrical potentials match actual measured electrical potentials in the presence of the metallic structures.

6. The method of claim 4 wherein applying the inversion code includes generating an electrical conductivity distribution that minimizes error between the simulated electrical potentials and the actual potentials measured in the presence of the metallic structures.

7. The method of claim 2 wherein applying the forward model and the inversion produces a reconstructed three-dimensional image of the subsurface, while removing the effects of the subsurface metallic structures.

8. The method of claim 1 wherein the electrodes are at least one of the following: surface electrodes, well casing electrodes, arrays of point buried electrodes, and the metallic structures with known locations and dimensions used as electrodes.

9. The method of claim 1 wherein a boundary of each metallic structure is represented by nodes, faces, or segments of unstructured tetrahedral elements.

10. A method of locating features including leaks and leaked materials beneath a subsurface containing metallic infrastructures such as pipes and tanks with known locations and dimensions, the method comprising: a. injecting current and measuring electrical potential of the subsurface with multiple sets of electrodes in the presence of the metallic infrastructures, thus generating electrical resistivity tomography measurements; b. applying a forward model to simulate the subsurface electrical potential measurement of those metallic infrastructures; and c. performing an inversion of the forward model to generate an electrical conductivity distribution that minimizes error between the simulated electrical potentials and actual measured electrical potentials, while removing the effects of the metallic structures on the measured tomography measurements.

11. The method of claim 10 wherein the applying the forward model to simulate the electrical potential measurement comprises a solution to a Poisson equation.

12. The method of claim 11 wherein the Poisson equation is solved using an immersed interface boundary condition solution to simulate subsurface potentials generated in the presence of the metallic structures.

13. The method of claim 12 wherein the Poisson equation includes at least one of the following: integral calculations, analytical calculations, numerical calculations, partial solutions calculations, digital calculations, and an immersed interface boundary condition solution.

14. The method of claim 11 wherein applying the forward model and performing the inversion produces a reconstructed three-dimensional image of the subsurface, while removing the effects of the metallic structures.

15. The method of claim 10 wherein the electrodes are at least one of the following: surface electrodes, well casing electrodes, arrays of point buried electrodes, and the metallic structures used as electrodes.

16. The method of claim 10 wherein a boundary of each metallic structure is represented by nodes, faces, or segments of unstructured tetrahedral elements.

17. A method of imaging electrical conductivity distribution of features within a subsurface containing metallic structures with known locations and dimensions, the method comprising the steps of: a. injecting electrical currents into the subsurface to measure electrical potentials using multiple sets of electrodes, thus generating actual electrical resistivity tomography measurements; b. generating simulated electrical resistivity measurements by applying a Poisson equation to simulate the measured potentials in the presence of the metallic structures, wherein a boundary of each metallic structure is represented by nodes, faces, or segments of unstructured tetrahedral elements; and c. applying an inversion code to the simulated and actual electrical resistivity measurements to removes the effects of the subsurface metallic structures with known locations and dimensions on the actual data; and d. e. producing a reconstructed three-dimensional image of the subsurface electrical conductivity in the presence of the subsurface metallic structures.

18. The method of claim 17 wherein the Poisson equation includes at least one of the following: integral calculations, analytical calculations, numerical calculations, partial solutions calculations, digital calculations, and an immersed interface boundary condition solution.

19. The method of claim 17 wherein the electrodes are at least one of the following: surface electrodes, well casing electrodes, arrays of point buried electrodes, and the metallic structures used as electrodes.

20. The method of claim 17 wherein the inversion code is applied so that the simulated electrical potentials match actual measured electrical potentials in the presence of the metallic structures.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 shows a diagram of a general subsurface region containing infinitely conductive inclusions Ω.sub.1 . . . Ω.sub.n, with surface boundaries Γ.sub.1 . . . Γ.sub.n, embedded within a finite conductivity host medium Ω.sub.0 with electrical conductivity σ.sub.0.

(2) FIG. 2 shows two types of tetrahedral elements that can form a boundary, along with corresponding interior planes (P1-P6) used for local CVFM current flux computations through an infinite conductivity inclusion.

(3) FIG. 3 shows a flow diagram illustrating the parallel computation of a forward simulation.

(4) FIG. 4 is a conceptual diagram for the analytic solution to the potential generated at a point from a line source pole current injection.

(5) FIG. 5A shows a conceptual model of a system of buried metallic storage tanks, transfer pipes, and monitoring wells.

(6) FIG. 5B is a cutout of the computational mesh used to model the system of FIG. 5A.

(7) FIG. 6 shows A) computational mesh for current injection between two line sources; B) partial solution for the line source at (0,0) m; C) partial solution for the line source at (20,0) m; D) Analytic solution—or potential distribution per ampere of current injected; E) numeric solution; F) percent error distribution; and G) percent error histogram.

(8) FIG. 7 shows A) computational mesh for a pole current source in the presence of an infinite conductivity vertical line—e.g. well casing; B) partial solution for the point current source; C) partial solution for the line source; D) analytic solution—or potential distribution per ampere of current injected; E) numeric solution; F) percent error distribution; and G) percent error histogram.

(9) FIGS. 8A and 8B show two examples of potential distributions generated in the presence of the metallic subsurface infrastructure (tanks, pipes, and wells) illustrated in FIGS. 5A and 5B.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

(10) The present invention relates to methods and systems of imaging the electrical conductivity distribution of a subsurface containing metallic structures with known locations and dimensions. In one embodiment of the present invention, the method accurately accommodates the large conductivity contrast typically existing between subsurface host materials and buried metallic structures. In one embodiment, the global solution to the Poisson equation is decoupled at the metallic boundaries into a series of well-conditioned partial solutions. Each partial solution is then scaled to honor known current flux conditions based on the divergence theorem. Finally, the partial solutions are superimposed to form the global solution.

(11) In addition, other advantages of the present invention include, but are not limited to, the following: 1) the solution can be formulated on an unstructured tetrahedral mesh, enabling arbitrary shapes to be accommodated efficiently; 2) discontinuous structures that are electrically connected can be easily accommodated; 3) metallic objects can be used as current sources and/or potential measurement electrodes; 4) no special requirements or adjustments are required for the inversion problem; and 5) a method of singularity removal may be provided without the necessity of an analytic solution (Weigmann and Bube, 1998), and 6) the method or algorithm can be implemented in parallel to accommodate the computational demands of 3D imaging. With respect to computational requirements, an additional forward solution may be provided for each distinct metallic object, which is the computational equivalent of adding one point source electrode to the inversion algorithm.

(12) In one example, the present invention allows for the imaging of contaminated soil beneath a leaking tank using surface electrodes, well casing electrodes, and/or vertical arrays of point source electrodes. The capability to accurately model infrastructure with known dimensions and location—as described in the present invention—enables ERT imaging to be accurately assessed and executed.

(13) Subsurface Electrical Potential

(14) The direct-current electrical potential arising from a given current source is given by the Poisson equation:
∇.Math.σ(r)∇φ(r)=J(r)  eq. (1)

(15) where r is a position vector, σ is electrical conductivity, φ is electrical potential, and J is the source current density. Given σ, J and appropriate boundary conditions, the objective of the forward model is to solve equation 1 for φ, which is measured at known electrode locations during an ERT experiment. Consider the general subsurface region Ω which is comprised of a background sub-region Ω.sub.0 with electrical conductivity σ.sub.0(r), and n sub-regions Ω.sub.1 . . . Ω.sub.n with infinite conductivity, such that the potential in each sub-region has a unique, constant potential, designate φ.sub.1 . . . φ.sub.n as shown in FIG. 1. In this example, the surface boundary (Γ.sub.s) condition is Neumann, whereby the current flux normal to the surface boundary is zero. In this example, the subsurface boundary conditions are Dirichlet (Γ.sub.ss), and specify zero potential on each subsurface boundary. In practice, the subsurface boundaries of the computational mesh are placed far from the region of interest to approximate zero potential at infinite distance. The sub-regions Ω.sub.1 . . . Ω.sub.n represent (for example) electrically conductive well casings, pipes, or tanks with known location and dimension that are assumed to have infinite conductivity. The error introduced by the infinite conductivity assumption is insignificant under most circumstances, as discussed below.

(16) The solution to equation 1 under the conditions shown in FIG. 1 is computed herein using, as one example, a specific application of the more general framework described by Klapper and Shaw (2006), which provides a generalization of the immersed interface boundary condition method for accommodating large coefficient discontinuities in the numerical solutions of elliptic and parabolic equations. In one embodiment, the global solution φ is generated by 1) solving a series of well-conditioned partial solutions that are decoupled at the discontinuous conductivity boundaries (i.e., Γ.sub.1 . . . Γ.sub.n), 2) scaling the partial solutions to honor a set of known current flux conditions, and 3) superimposing the partial solutions to provide the global solution. Below is shown the development for the global solution arising from a single current source with a current sink at infinite distance from the source, termed herein the pole solution. To simulate an ERT data set, one such solution is generated for every electrode used as a current source or sink. The total potential for any source-sink electrode pair may then be generated by subtracting the pole solution for the current sink electrode from the pole solution for the current source electrode. This approach provides a significant computational savings over explicitly computing the global solution for each measurement, the number of which is typically much larger than the number of current source electrodes.

(17) The pole solution for a point electrode located in Ω.sub.0 comprises n+1 partial solutions, one for the electrode and one for each of the n infinite conductivity sub-regions. The partial solution for region Ω.sub.0 (φ.sub.0) is given by solving:
∇.Math.σ(r)∇φ.sub.0(r)rεΩ.sub.0
φ(r)=0 rεΓ.sub.1 . . . Γ.sub.n
σ(r)∇φ(r).Math.{circumflex over (n)}=0 rεΓ.sub.s
φ(r)=0 rεΓ.sub.ss  eq. (2)

(18) where {circumflex over (n)} is the unit normal vector to the corresponding surface. In equation 2, the partial solution in region Ω.sub.0 is computed for given current source in Ω.sub.0 (e.g. an electrode point source), with the potential at each infinite-conductivity boundary set to zero, and the outer boundary conditions as previously specified.

(19) The partial solution (φ.sub.j, j=1 . . . n) for region Ω.sub.0, is given by solving:
∇.Math.σ(r)∇φ.sub.j(r)=0 rεΩ.sub.0
φ(r)=0 rεΓ.sub.1 . . . Γ.sub.n,r.Math.Γ.sub.j
φ(r)=1 rεΓ.sub.j
σ(r)∇φ(r).Math.{circumflex over (n)}=0 rεΓ.sub.s
φ(r)=0 rεΓ.sub.ss  eq. (3)

(20) The solution to equation 3 provides the potential in Ω.sub.0 given a unit potential on the infinite conductivity boundary Γ.sub.j, zero potential on the remaining infinite conductivity boundaries, and outer boundary conditions as previously specified.

(21) The partial solutions φ.sub.0 . . . φ.sub.n provide the basis desired to construct the global solution. However, direct superposition of φ.sub.0 . . . φ.sub.n will provide a global solution that violates the known flux conditions through each infinite conductor. Those conditions are given by the divergence theorem, which specifies:

(22) $\begin{array}{cc}{\int }_{{\Gamma }_j}\sigma \left(r\right)\nabla \varphi \left(r\right).Math.\hat{n}\mathrm{dS}={\int }_{{\Omega }_j}J\left(r\right)\mathrm{dV}& \mathrm{eq}.\left(4\right)\end{array}$

(23) In this work there are two possible source conditions for equation 4. In the first case there is no current source within Ω.sub.j, and the right hand side of equation 4 is zero. This represents the case were Ω.sub.j is not being used as a current source electrode. When region Ω.sub.j is the current source electrode, the right hand side of equation 4 is set to one in order to simulate ERT potential measurements which are typically normalized to one unit of injected current.

(24) The partial solutions are assembled to form the global solution by the equation:
φ(r)=φ.sub.0(r)+Σ.sub.j=1.sup.nC.sub.jφ.sub.j(r)  eq. (5)

(25) where the scaling factors C.sub.j are chosen to satisfy the flux conditions specified by equation 4. The scaling factors are determined by solving the matrix equation:

(26) $\begin{array}{cc}\left[\begin{array}{cccc}{f}_{1,1}& f_{1,2}& \mathrm{.Math.}& f_{1,n}\\ f_{2,1}& & & \\ \mathrm{.Math.}& & \ddots & \\ f_{2,n}& & & f_{n,n}\end{array}\right]\left[\begin{array}{c}{C}_1\\ C_2\\ \mathrm{.Math.}\\ C_n\end{array}\right]=\left[\begin{array}{c}{F}_1\\ F_2\\ \mathrm{.Math.}\\ F_n\end{array}\right]& \mathrm{eq}.\left(6\right)\end{array}$

(27) where f.sub.i,j is the current flux through Ω.sub.i arising from φ.sub.j(r), computed as described in the forthcoming sections. F.sub.i is chosen to satisfy equation 4. For example, if the current source is a point electrode within Ω.sub.0, then:
F.sub.i=−f.sub.i,0, i=1 . . . n  eq. (7)

(28) so that the net current flux out of Ω.sub.i, i=1 . . . n arising from φ.sub.0 is negated by the cumulative current flux out of Ω.sub.i, i=1 . . . n arising from C.sub.jφ.sub.j(r), j=1 . . . n.

(29) If region Ω.sub.k, k=1 . . . n is a non-point current injection electrode, then:
F.sub.i=0, if i≠k
F.sub.i=1, if i=k  eq. (8)

(30) so that the net current flux from Ω.sub.k is one unit, and the current flowing out of the remaining infinite conductivity regions is zero. Note also that in this case, φ.sub.0(r) is zero because there is no current source within Ω.sub.0.

(31) Let np be the number of point source electrodes and nnp be the number of non-point source electrodes used to generate an ERT survey. Inspection of equations 2 through 6 reveals that the forward simulation for an ERT survey requires np solves of equation 2 (one for each point source electrode), n solves of equation 3 (one for each infinite conductivity region), and np+nnp solves of equation 6 (one for each non-point current source).

(32) Numerical Solution

(33) Equations 2 and 3 are discretized and solved using the weak form finite element solution on an unstructured tetrahedral mesh (Gunther et al., 2006). The unstructured mesh enables arbitrary shapes to be accommodated, and facilitates efficient modeling of complex structures (e.g. piping systems, buried tanks etc.). Potentials are computed on the nodes of the mesh, and conductivities are specified on the mesh elements, each element being formed by four bounding nodes. As will be shown, linear structures (e.g. pipes, wellbore casings) can be accurately modeled as a series of connected nodes (e.g. a line source), and 2D structures (e.g. sheet piles or other planar metallic barriers) can be modeled as a mesh of connected nodes (e.g. a planar source). In addition, because the potential of an infinite conductivity feature is everywhere equal, it is only necessary to model the boundary of 3D structures, with boundary conditions as described in the previous section. Thus only the shells of 3D structures are incorporated into the computational mesh, the internal structure being hollow. FIG. 5B shows an example of a computational mesh incorporating buried tanks, pipes, and well casings that will be used in the forthcoming synthetic modeling and inversion demonstration.

(34) Flux Computations

(35) To solve equation 6 for the flux scaling factors (C.sub.j), the current flux terms f.sub.i,j comprising the left hand side matrix is computed. It is well known that finite element flux computations are globally conservative, but locally non-conservative, eliminating the possibility of computing f.sub.i,j by the weak form finite element solution. Instead, localized control volume finite element (CVFM) computations (Voller; 2009) are used to compute the flux through each infinite conductivity inclusion.

(36) To illustrate the CVFM computations, consider the neighboring tetrahedral elements A and B shown in FIG. 2, which illustrate the two types of elements that can bound a point, line, plane, or volume representing an infinite conductivity feature. Element A illustrates the case where only two nodes of the element lie on an infinite conductivity boundary (N3 and N4, or N1 and N5). Element B illustrates the case where either three nodes (N2, N3,or N4), or 1 node (N1), lies on the boundary. Note that it is not possible for all four nodes to lie on the boundary with a non-degenerate tetrahedron. Typically many elements of types A and B will be combined to bound a given infinite conductivity structure.

(37) In the CVFM method, six sub-planes (P1-P6) are formed within each bounding element by connecting segment midpoints, element face centroids, and the nodal centroid in the case of element B as shown in FIG. 2. These subplanes are defined on every element having nodes that bound an infinite conductivity surface Γ.sub.i. As a whole, they form a closed boundary around Ω.sub.i by which the left hand side of equation 4 may be approximated to compute f.sub.i,j. To illustrate, let e represent an element bounding Ω.sub.i, and let ne be the total number of elements bounding Ω.sub.i. Equation 4 is approximated by:

(38) $\begin{array}{cc}{\int }_{{\Gamma }_i}\sigma \left(r\right)\nabla \varphi \left(r\right).Math.\hat{n}\mathrm{dS}\approx \underset{e=1}{\overset{\mathrm{ne}}{.Math.}}{\sigma }_e\underset{p=1}{\overset{6}{.Math.}}{A}_{e,p}.Math.\nabla {\varphi }_{e,j},{\hat{n}}_{e,p}.Math.=f_{i,j}& \mathrm{eq}.\left(9\right)\end{array}$

(39) where σ.sub.e is the conductivity of element e, A.sub.e,p is the area of plane p in element e, φ.sub.e,j is the potential within element e arising from source j, {circumflex over (n)}.sub.e,p is the outward pointing normal vector of plane p in element e, and represents the inner product. Note that φ.sub.e,j is a solution of equation 2 or 3 as described above, and ∇φ.sub.e,j is constant within each element, computed using the linear finite element shape functions for each element. In words, equation 9 computes and sums the outward flux flowing through every subplane to provide the total net flux flowing through Ω.sub.i due to the potential field φ.sub.j(r). For computational efficiency, equation 9 is reduced to the form:

(40) $\begin{array}{cc}{f}_{i,j}=\underset{e=1}{\overset{\mathrm{ne}}{.Math.}}{\sigma }_e\underset{k=1}{\overset{4}{.Math.}}{K}_{e,k,i}{\varphi }_{e,k,j}& \mathrm{eq}.\left(10\right)\end{array}$

(41) where k spans the four nodes in element e, φ.sub.e,k,j is the simulated potential in node k of element e, arising from source j, and K.sub.e,k,j is row e and column k of a flux computation matrix K.sub.i formed to account for the geometric terms of equation 9 for region Ω.sub.i. Since K.sub.i depends only on the mesh geometry, it is computed only once, and provides a means to efficiently compute the net flux through Ω.sub.i given any potential field φ.sub.j(r) and conductivity distribution σ(r). As shown below, computation of K.sub.i facilitates efficient computation of the forward problem in the parallel algorithm. This is particularly helpful in the inverse problem, where the forward problem must be solved multiple times.

(42) Parallel Implementation

(43) To summarize, steps to compute the forward solution include:

(44) 1. Compute the np point source partial solutions φ.sub.0(r), one for each point electrode (equation 2).

(45) 2. Compute the n non-point partial solutions φ.sub.j(r), j=1 . . . n, one for each infinite conductivity region (equation 3).

(46) 3. Using the partial solutions computed in parts 1 and 2, compute the np+npp solutions to equation 6, one for each np point source/sink electrode and each npp non-point current source/sink electrode.

(47) 4. Construct the np+npp pole solutions using equation 5.

(48) 5. Superimpose appropriate pole solutions and extract potentials to construct simulated measurements.

(49) Steps 1 through 5 are amenable to parallel computation. The parallel algorithm described herein uses a master/slave configuration (Johnson et al., 2010), whereby a single master process handles input/output and orchestrates computations. One or more slave processes handle the primary computational burdens, message passing, and memory requirements.

(50) For the sake of simplicity, assume one slave process is assigned to each of the np point source electrodes, and to each of the n infinite conductivity regions. FIG. 3 shows a flow diagram illustrating the parallel computation of a forward simulation. Rectangles represent the master process and circles represent slave processes. Each slave step is perfectly parallel—no inter-process communications are required—except for step 3, which represents a small portion of the overall computational burden, enabling the algorithm to remain highly scalable. In the setup phase the master process reads input information (i.e. computational mesh, electrode locations, measurement sequence etc.) and broadcasts pertinent information to each of the slave processors. The slave processors construct the finite element matrix. The slave processes responsible for infinite conductivity inclusions (slaves np+1 through np+n) also construct the flux computation matrices K.sub.i for regions Ω.sub.i, i=1 . . . n . Steps 1 and 2, whereby the partial solutions for each point electrode source and each infinite conductivity source are computed, are executed at the same time. In step three, message passing between slaves is required in order to efficiently compute the left and right hand sides of equation 6. In step 4, each slave uses equation 5 to compute the pole solution for its corresponding electrode. Finally, in step 5, each slave extracts those potentials necessary to simulate the ERT measurement sequence and sends them to the master process where they are assembled.

(51) In this work the partial solutions (i.e. equations 2 and 3) are not parallelized, enabling the forward computation to remain highly scalable, except for the message passing required in step 3. The matrix equation required to solve equations 2 and 3 is sparse, such that relatively large problems (many millions of elements) can easily be computed stored on the memory typically allocated to a single processor. The message passing requirements to formulate and solve equation 6 are typically small such that the computation burden associated with steps 3 is much smaller than that of steps 1 and 2, even for large numbers of infinite conductivity inclusions. Thus, message passing requirements for the forward simulation are minimal, enabling the algorithm to remain highly scalable.

(52) Inverse Problem

(53) Analytic Solutions

(54) In the results section numeric forward modeling results using the procedure described above are validated through comparison to analytic solutions for a line source. Referring to FIG. 4, the potential at a point within an infinite medium of homogeneous conductivity caused by a vertical line source whose axis passes through the origin is given by (Wait, 1982, Bhagwan and Trofimenkoff, 1982):

(55) $\begin{array}{cc}\varphi \left(r,z\right)=\frac{I_s}{8\pi l\sigma }\ln \left[\frac{\left(z_2-z\right)+\sqrt{{\left(z_2-z\right)}^2+r^2}}{\left(z_1-z\right)+\sqrt{{\left(z_1-z\right)}^2+r^2}}\right]& \mathrm{eq}.\left(11\right)\end{array}$

(56) where z,z.sub.1,z.sub.2, l and r are distances as outlined in FIG. 4, I.sub.s is the current source strength of the line source, and σ homogeneous conductivity. Using the method of images, the corresponding half space potential is given by:

(57) $\begin{array}{cc}\varphi \left(r,z\right)=\frac{I_s}{8\pi l\sigma }\ln \left[\frac{\left(z_2-z\right)+\sqrt{{\left(z_2-z\right)}^2+r^2}}{\left(z_1-z\right)+\sqrt{{\left(z_1-z\right)}^2+r^2}}+\frac{\left(z_2+z\right)+\sqrt{{\left(z_2+z\right)}^2+r^2}}{\left(z_1+z\right)+\sqrt{{\left(z_1+z\right)}^2+r^2}}\right]& \mathrm{eq}.\left(12\right)\end{array}$

(58) where the half-space boundary is located at z=0.

(59) Using equations 11 and 12, a semi-analytic solution for the halfspace potential arising from a point source of current in the presence of a vertical conductor may be derived. The solution is similar to the numeric procedure described above as outlined in detail by Johnston et al 1987. The vertical conductor is divided into a number of small segments, and the current strength of each segment is adjusted so that the total current flux emanating from the line source is perfectly negated by the current flux through the line source that arises from the point source potential, thereby honoring conservation of charge within the conductor. Once the current strength of each segment is known, the potential distribution arising from that segment is computed with equation 12, and the global potential is produced by summing the contributions from each segment with the point source potential.

(60) Synthetic Analysis: Imaging Vadose Zone Contamination in the Presence of Buried Conductive Pipes, Tanks, and Well Casings

(61) To demonstrate a real-world application, consider the conceptual model shown in FIG. 5A, which consists of a series of six metallic storage tanks connected together by buried metallic piping, such that the tanks and pipes have the same electrical potential. The pipes are tanks connected physically and electrically to form a single ‘electrode’ of constant potential. The conductivity of soils surrounding the tanks and pipes is 0.001 S/m, except for the regions where leaks have occurred. The three pipe-leak regions have a conductivity of 0.1 S/m, and the tank leak region has a conductivity of 1.0 S/m. The system also includes 14 metal-cased monitoring wells to a depth of 40 m, distributed as shown in FIG. 5A. The objective of the analysis is to investigate the capability to image the leak zones using 1) surface electrodes deployed directly above the transfer pipes, 2) well casings as long electrodes, and 3) borehole electrodes installed beneath the leaking tank. The surface and borehole electrodes are shown as dots, and well casings are used as electrodes for ERT imaging. The surface array, in this example, consists of 43 electrodes over the main transfer line, and 39 electrodes over each of the three lines connecting to the six tanks, for a total of 160 electrodes at 2m spacing. Each of the four boreholes arrays contains 16 electrodes at 2 m spacing, for a total of 64 electrodes. Wenner array surveys were simulated along each of the four surface lines, for a total of 815 measurements. The borehole survey included a series of 3995 cross-hole dipole-dipole measurements. The well-casing survey comprised every possible unique measurement, including potential measurements on the tank/pipe ‘electrode’, for a total of 4095 measurements using 15 electrodes (14 wells and 1 tank/pipe electrode).

(62) A cutout of the computational mesh is shown superimposed on the tanks in FIG. 5B. As discussed previously, only the conductive boundary of the tank is included in the mesh, the interior tank regions being hollow. The pipes and wells are modeled as line sources with 0.5 m spacing between each node on the line, which produces a refined mesh around pipe and well locations. The mesh is also refined near point electrodes for modeling accuracy in the high potential gradient regions near electrodes. Mesh boundaries extend far beyond the imaging region (˜5000 m) to reduce boundary effects (not shown), with expanding elements to reduce computational demand. In total, the mesh consists of 1,567,844 elements and 254,409 nodes.

(63) Forward modeling and inversion are executed on the same mesh using 161 processors for the surface array simulation, 16 processors for the well-electrode array simulation, and 65 processors for the borehole array simulation. To simulate field conditions, 5% normally distributed noise is added to each measurement.

(64) Results

(65) FIG. 6 shows a comparison of the analytic and numeric solutions computed for a current injection of 1 Ampere between two line sources (e.g. well casings). FIG. 6A shows a cutout of the computational mesh and the locations of the current source and sink lines. FIGS. 6B and 6C show the two partial solutions to equation 2 need to form the pole solutions, one for each line source. FIG. 6B shows the partial solution for the line source at (0,0) m. FIG. 6C shows the partial solution for the line source at (20,0) m. For each partial solution the potential on the source line is one volt and the potential on the remaining line is zero as required. As described previously, the partial solutions are scaled and summed so that the net current flux emanating from the pole line source is equal to one ampere and net flux through the remaining line source is zero, thereby providing the pole solution for the given line source. The pole solution for the line source at x=20 m is subtracted from the pole solution for the line source at x=0 m to provide the global solution shown in FIG. 6E, which is a numeric solution. For comparison, FIG. 6D shows the corresponding analytic solution— i.e., potential distribution per ampere of current injected—computed using equation 12. FIG. 6F shows the numeric solution error, which is given by:

(66) $\begin{array}{cc}\%{\mathrm{Err}}_i=100\left(\frac{{\varphi }_i-{\varphi }_{t,i}}{{\varphi }_{t,i}}\right)& \mathrm{eq}.\left(13\right)\end{array}$

(67) where φ.sub.i,j and φ.sub.i are respectively the analytic and numeric potential solutions for node i. FIG. 6G shows the error histogram.

(68) FIG. 7 compares the semi-analytic solution for a point current pole source in the presence of a vertical line of infinite conductivity. As with FIG. 6, the FIG. 7 shows: A) the computational mesh for a pole current source in the presence of an infinite conductivity vertical line (e.g., a well casing), B) the partial solution for the point electrode source (see equation 2), C) the partial solution for the line source (see equation 3), D) the semi-analytic solution (i.e., potential distribution per ampere of current injected), and E) the numeric solution, F) the percent error distribution or percent difference between the numeric and semi-analytic solution, and G) the percent difference histogram. Both FIGS. 6 and 7 display good agreement between analytic and numeric solutions, and demonstrate the accuracy of the numeric solution.

(69) FIG. 8 shows two examples of potential distributions generated in the presence of the metallic subsurface infrastructure (tanks, pipes, and wells) illustrated in FIG. 5. FIG. 8A shows the potential generated from a 1 Ampere current injection between point source electrodes on the surface, and 8B shows the potential generated from a 1 Ampere current injection between two well casing line sources. Current flows perpendicular to constant potential isosurfaces, so the potentials are presented as isosurfaces to illustrate the nature of current flow through the subsurface. Careful inspection of 8A and 8B shows how the pipes and tanks provide primary current flow paths between source and sink electrodes, dramatically altering the potential field that would exist in the absence of the infrastructure. Isosurfaces are shown at +/−(10.0, 1.0, 0.1, and 0.001 V/A). Pipes and tanks are electrically connected. As current flows perpendicular to the isosurfaces, the isosurfaces indicate how current flows into and out of metallic features, and how those features influence the potential field. Noting the dominant influence the infrastructure exerts on the potential field, it is valid to conclude that erroneous modeling of infrastructure effects in this case and similar cases, such as assuming tanks are electrically isolated, will likely result in significant inversion artifacts as the inversion algorithm attempts to compensate for modelling error.

(70) In compliance with the statute, embodiments of the invention have been described in language more or less specific as to structural and methodical features. It is to be understood, however, that the entire invention is not limited to the specific features and/or embodiments shown and/or described, since the disclosed embodiments comprise forms of putting the invention into effect.