Thin and highly conductive objects are challenging to model in 3D direct-current (DC) problems because they often require excessive mesh refinement that leads to a significant increase in computational costs. RESistor network (RESnet) is a novel algorithm that converts any 3D geo-electric simulation to solving an equivalent 3D resistor network circuit. Two features of RESnet make it an attractive choice in the DC modeling of thin and conductive objects. First, in addition to the conductivity with units of Siemens per meter (S/m) defined at the cell centers (cell conductivity), RESnet allows conductive properties defined on mesh faces and edges as face conductivity with units of S and edge conductivity with units of S·m, respectively. Face conductivity is the thickness-integrated conductivity, which preserves the electric effect of sheet-like conductors without an explicit statement in the mesh. Similarly, edge conductivity is the product of the cross-sectional area and the intrinsic conductivity of a line-like conductive object. Modeling thin objects using face and edge conductivity can avoid extremely small mesh grids if the DC problem concerns electric field responses at a much larger scale. Second, once the original simulation is transformed into an equivalent resistor network, certain types of infrastructure, similar to above-ground metallic pipes, can be conveniently modeled by directly connecting the circuit nodes, which cannot interact with each other in conventional modeling programs. Bilingually implemented in MATLAB and Python, the algorithm has been made open source to promote wide use in academia and industry. Three examples are provided to validate its numerical accuracy, demonstrate its capability in modeling steel well casings, and indicate how it can be used to simulate the effect of complex metallic infrastructure on DC resistivity data.

The direct-current (DC) resistivity method, or electric resistivity tomography (ERT), is a well-established technique in characterizing the surface using electric resistivity (or its reciprocal conductivity), especially in the areas of environmental, engineering, and near-surface investigations. However, those applications often encounter manmade metallic infrastructure that can strongly influence electric data (Holladay and West, 1984). For example, the long electrode ERT method (Rucker et al., 2011; Zhang et al., 2014) and the top-casing source method (Wilt et al., 2020) use steel well casings as extended probes to detect subtle conductivity changes in reservoirs. ERT surveys carried out in industrial or urban areas also may be interfered with by pipelines, powerlines, buried tanks, and rebar (Johnson and Wellman, 2015; Yang et al., 2021). These objects are usually thin in one or two of the three spatial dimensions but are orders of magnitude higher in conductivity compared with natural geologic units. Such objects create strong local anisotropy that can significantly distort the subsurface current distribution and appear as confusing anomalies in the DC resistivity data.

Accurate forward simulation in 3D is a critical tool in quantitatively understanding the effect of metallic infrastructure and other thin conductive objects. The most straightforward approach is to directly model reality using an extremely fine mesh (Commer et al., 2015; Hoversten et al., 2015; Heagy and Oldenburg, 2019). Those methods are particularly good at capturing the complete physics and simulating the field inside and near the conductors. Another type of approach is to approximate the effect of thin conductors at the cost of losing a certain degree of accuracy inside or near the conductors. The most common method is to assign an equivalent high conductivity to the cells at the approximate location of infrastructure (Ronczka et al., 2015; Um et al., 2015; Weiss et al., 2016; Yang and Oldenburg, 2017). Metallic objects also can be treated as equipotential, then the object’s effect can be approximated by distributing the current source over the objects (Rocroi and Koulikov, 1985; Zhu and Feng, 2011; Yang et al., 2017) or by imposing an immersive interface boundary condition (Johnson and Wellman, 2015). Steel well casings are of great significance in many applications, so many novel methods have been proposed to numerically model casing’s electric effect, e.g., the method of moment (Kohnke et al., 2018), the hierarchical finite-element method (Weiss, 2017), and the slender prismatic elements method (Yu et al., 2022).

Our algorithm, RESistor network (RESnet), takes a different route by considering the subsurface as a 3D equivalent resistor network. The most famous circuit network analog in geophysics was made by Theodore Richard Madden (1925–2013) in his technical report “the Gray Peril” (Madden, 1972; Williams et al., 2014), in which he created insightful connections between transmission network systems and geophysical forward and inverse problems (Mackie et al., 1988). The general idea of solving the DC resistivity modeling problem using an equivalent resistor network continued to be adopted by Zhang et al. (1995) and Gerami (2018), although their specific solving techniques may differ. The equivalent resistor network approach behind the RESnet program, as previously presented in Yang et al. (2016) and Li and Yang (2021), has many versatile and valuable features that supplement the existing algorithms with unique functionalities specifically designed for metallic infrastructure and other conductive objects in narrow and thin shapes. An extension of RESnet to frequency-domain electromagnetic (EM) also has been proposed to efficiently simulate EM induction in the presence of casings (Hu et al., 2022).

In the following, we first introduce the algorithm with numerical implementations. Then, the design of the software is presented to enable users to quickly adopt the modeling technique for their specific problems. Finally, three examples of application are shown: (1) a pole-dipole survey over a half-space and comparison against the analytic solution; (2) a top-casing source configuration for reservoir monitoring and reproducing previously published results; (3) a dipole-dipole survey in an environment with complex metallic infrastructure, such as a steel cased well, a steel sheet warehouse, and an above-ground pipeline connecting them. This paper is written primarily as a tutorial to help users understand the physics behind the codes and navigate the utilization of the RESnet software for practical problems. The development of the theory and formulation of numerical approaches can be found in previous publications (Yang et al., 2016; Li and Yang, 2021).

Equivalent resistor network

RESnet uses a rectilinear mesh to build the resistor network because the nodes and edges (Figure 1a) in a mesh can be converted to the nodes and branches (Figure 1b) in a resistor network circuit. One of the apparent benefits of using the circuit analog is that when thin and elongated objects, similar to steel well casings, coincide with the circuit branches (Figure 1c), such objects can be efficiently and conveniently modeled by adding a parallel-circuit resistor of high conductance to the existing branch (Figure 1b) (Daily et al., 2004).

RESnet explicitly implements the functionality preceding by allowing the assignment of conductive properties to mesh cell centers, mesh faces, and mesh edges. Ronczka et al. (2015) refer to a similar technique as the “shunt electrode model” for casings. Based on a highly flexible, unstructured tetrahedral mesh, Weiss (2017) propose the hierarchical finite-element method, which can handle conductivities assigned to arbitrarily shaped and oriented faces and edges. With a rectilinear 3D mesh, conductive objects are always represented by cuboids. A cuboidal conductor can conduct current in one of the three orthogonal directions in the coordinate system, and its spatial dimensions are measured by the width (w), height (h), and length (l) (Figure 2). Depending on the relative size of width, height, and length, RESnet defines the following conductive properties.

  1. Cell conductivity (cellCon, σc in Siemens per meter [S/m]): defined at cell centers for volumetric objects with comparable width, height, and length (Figure 2a). CellCon describes arbitrarily shaped geologic objects similar to standard DC modeling algorithms.

  2. Face conductivity (faceCon, σf in S): defined on cell faces for sheet-like objects, similar to thin layers of sediment or fracture planes, whose width (or height) is much smaller than the length and height (or width) (Figure 2b). Face conductivity is the product of its intrinsic conductivity and the height (σf=σc×h).

  3. Edge conductivity (edgeCon, σe in S·m): defined on mesh edges for line objects whose width and height are much smaller than the length (Figure 2c). Edge conductivity is the product of its intrinsic conductivity and the cross-sectional area (σe=σc×w×h).

RESnet collapses the three types of conductivity to the conductance defined on edges. By definition, the directional conductance of a conductor with a uniform cross section is
C=σcwhl=σfwl=σe1l,
(1)
where w,h, and l are the width, height, and length of the object, respectively (Figure 2). An edge has four neighboring cells, so the edge conductance contains contributions from those four cells in parallel (Figure 2a). Because the conductivity σc and the edge length l are known, the cross-sectional area (w×h) of each of the four conductors can be calculated using the conductor’s volume, which is one-fourth of the total volume of the cell it belongs to. Using equation 1, the conductance from one neighboring cell is calculated by
Cc=V4l×σcl=Acσc,
(2)
where V=w×h×l is the volume of the corresponding cell, then the total conductance on an edge from the cell conductivity is the sum of the four conductors (the red-colored volume in Figure 2a). Similarly, an edge has four neighboring faces, and the conductance from one of the faces is
Cf=A2l×σfl=Afσf,
(3)
where A=w×l is the area of the corresponding face (the red-colored area in Figure 2b). Finally, an edge conductivity also contributes to the edge conductance according to
Ce=σel=Aeσe.
(4)
Treating σc,σf, and σe as three model vectors, the vector of total conductance assigned to the resistor on each branch of the circuit is the conductance from Cc,Cf,Ce in a parallel-circuit
C=Acσc+Afσf+Aeσe,
(5)
where Ac,Af, and Ae are the matrix representations of the calculation in equation 24, respectively.

Solutions

After the derivation of an equivalent circuit, simulating a DC resistivity survey is equivalent to connecting current sources to some nodes of the circuit and measuring the potential difference between two nodes. RESnet defines the potential φ on the nodes, then the potential difference across a branch can be computed using a differential matrix G. The current flowing along a branch is the product of the potential difference and the conductance. Finally, using Kirchhoff’s current law (KCL), we have the linear system of forward modeling as
GTdiag(C)Gφ=s,
(6)
where s indicates where the current source is connected to the circuit. Equation 6 appears similar to the ·σ equation used in most DC codes, but there is a meaningful difference. In regular DC codes, the geometric information about the mesh is embedded in the gradient and divergence matrix, so only cell conductivity can be modeled. In contrast, RESnet incorporates the mesh information into the conductance term to allow the modeling of thin objects by using face and edge conductivities. The formulation of an equivalent circuit implicitly implies a no-flux boundary condition because the current cannot escape from the circuit.

Equation 6 does not have a unique solution for the potential. As a result, we modify the coefficient matrix of equation 6 by adding one to the first entry of the first row, which forces the potential at the first node to be zero. The modified coefficient matrix is positive definite, so RESnet factorizes the original sparse matrix to a lower triangular matrix, then obtains the solution of φ through quick back substitution. When there are many different right sides representing different sources, which is common in DC resistivity surveys, the direct method is more efficient because one matrix factorization can be used for the parallel solutions of different sources. Once the potentials are computed at all the nodes, the potential differences and current intensities on all the branches also can be derived. Interested readers and potential users are referred to the appendices to learn more about the design of software and the explanation of codes in detail.

RESnet is designed to conveniently, efficiently, and accurately simulate the effect of thin conductors. One of the motivating applications is the steel well casing in oilfields. This section quantitatively demonstrates the superior performance of RESnet in comparison to the conventional finite-difference (FD) method in solving geo-electric problems involving steel casings.

The formulation of RESnet in the discretized form is equivalent to the FD method using rectilinear or structured meshes, which has been established as the standard of 3D DC resistivity modeling. However, RESnet can be considered as an extension of the regular FD with the additional capability of modeling the casing using edge conductivity. In contrast, when only the regular FD code is available, it was a common practice to approximate the casing using a string of very conductive volumetric cells, colloquially called the “conductive pillar” method. Here, the edge conductivity in RESnet is compared with the conductive pillars of different thicknesses and the FD solution using extremely fine cells down to a few millimeters.

The geo-electric model under investigation is a 1000 m long vertical steel casing buried in a 0.01 S/m uniform half-space. The casing pipe has an outer diameter of 10 cm and an inner diameter of 8 cm, resulting in a cross-sectional area of approximately 0.0028 m2; the steel’s intrinsic conductivity is assumed to be 5 × 106 S/m; the edge conductivity assigned to the edges representing the well path is 14,137 S·m. The synthetic DC resistivity survey deploys one end of the DC current supply to the top of casing at (0, 0, 0) and places the return electrode 2000 m away at (−2000, 0, 0). Ten potential electrodes are aligned along a survey line from (10, 0, 0) to (100, 0, 0) at a uniform spacing of 10 m. Such an electrode array produces a set of pole-dipole data of potential difference.

The base mesh contains three levels of refined core volumes depending on the distance to the well head. The first level for the range of X = [−1220, 1220], and Y = [−1220, 1220] is discretized by 40 m cells; the second level for the range of X = [−500, 500], and Y = [−500, 500] has 20 m cells; the third level for the range of X = [−100, 100], and Y = [−100, 100] uses 10 m cells. The core volume is vertically discretized by 50 m cells for the top 1 km. Outside of the core volume, the mesh cells expand by a factor of 1.4 until the modeling domain is truncated at approximately 10 km on the side or at the bottom. The background model is first simulated to obtain the pole-dipole data for the scenario without a casing, which serves as a reference for a better understanding of the casing effect (Figure 3a). Then, the following forward simulations are carried out:

  1. Edge conductivity using RESnet. The vertical edges along the 1000 m wellbore are assigned an edgeCon of 14,137 S·m. No further mesh refinement is needed.

  2. Thick conductivity pillar using the regular FD. The cell column of X = [−10, 10], Y = [−10, 10], Z = [−1000, 0] is filled with a conductivity of 35.34 S/m, which retains the same cross-sectional area-conductivity product. No further mesh refinement is needed.

  3. Thin conductivity pillar using the regular FD. An extra level of refinement representing the casing (1 m horizontal cell size) is implemented for the range of X = [−1, 1] and Y = [−1, 1]. This refined volume is filled with a conductivity of 3534 S/m, which still retains the same cross-sectional area-conductivity product as in (1) and (2).

  4. Exact pipe geometry using the regular FD. Additional refinement involving extremely fine cells of 8.86 mm is applied to exactly depict the actual geometry and the hollow structure of casing. The casing cells are assigned the steel’s intrinsic conductivity of 5 × 106 S/m; the noncasing cells in the refined region have the background conductivity.

Half-space

For a test of the numerical accuracy, the first example simulates a pole-dipole survey over a uniform half-space, whose analytic solution exists. The script for the first example is Example_Halfspace.m or Example_Halfspace.py. The survey places the positive current source (A electrode) at the origin; ten measurement electrodes are 10 m equally spaced along a line in the x-direction from X = 10 m to 100 m (Figure 4).

The core mesh region of X = [−4, 104], Y = [−4, 4], and Z = [−4, 0] under and around the survey line is uniformly discretized with the smallest cell size 2 m; the cells gradually coarsen at a rate of 1.4 from the core mesh to the side and bottom boundary approximately 4 km away. The expansion rate and distance to the boundary condition are two empirical parameters in mesh design; the former is usually between 1.1 and 1.5, depending on the complexity of the model and the trade-off between accuracy and efficiency; the latter must be chosen so that the boundary effect in the long offset data is negligible. The entire mesh is uniformly filled with a cellCon of 0.01 S/m; no faceCon or edgeCon is applied. The script explicitly populates the return (B) electrode at X = −inf, but inside the code, the B electrode is automatically snapped to the node with the smallest possible x-coordinate. When the modeling domain is sufficiently large, the approximation to the pole source is deemed accurate sufficient.

The script runs for less than a second on a desktop computer with a 3.2 GHz Intel Core i7-8700 processor. Nine pole-dipole data points are computed using the 10 electrode locations (the white dots in Figure 4). The simulated potential difference data are compared with the analytic solutions with a satisfactory agreement (Figure 5a); the relative error, calculated as the difference between the analytic and RESnet normalized by the analytic, is less than 2% for the entire survey line (Figure 5b). Further reduction of the modeling error can be achieved by refining the mesh or replacing the trilinear interpolation with higher order formulas. Additional tests on layered earth models also are carried out to ensure RESnet has competent performance on conventional 3D DC modeling tasks that only involve the cell conductivity.

Steel casing

RESnet can efficiently model thin and elongated objects, similar to steel well casings, pipelines, and powerlines, with the concept of edgeCon. In the second example, Example_Casing.m or Example_Casing.py, we reproduce some simulation exercises in previous publications so that the results can be benchmarked. The example, taken from Figure 5 in Heagy and Oldenburg (2019), uses a top-casing source configuration, in which the positive current (A) electrode is in direct contact with the head of a steel casing. Such configuration has the advantage of channeling a substantial amount of current down to greater depth for enhanced sensitivity to subtle changes in conductivity due to injected fluid or other events (Li and Yang, 2021). In this particular example, a 1000 m long vertical steel casing is buried in a uniform background subsurface. The A electrode is fixed to the top of the casing at the origin, but the return (B) electrode is placed at four different offsets, X = −2000 m, −750 m, −500 m, and −250 m (Figure 6). The strength of the surface electric field in the area of X = [−1250, 1250] and Y = [−1250, 1250] is visualized to investigate the effect of buried steel casing.

Our 3D rectilinear mesh contains the smallest cell size of 50 m in the core mesh region specified by X = [−1300, 1300], Y = [−1300, 1300], and Z = [−1000, 0]; the discretization of the core mesh region is uniform, although our experience indicates that fine grids at depth may not be necessary if only the surface field is of interest. Outside of the core mesh region, the cells gradually coarsen at a rate of 1.4 until the presumed boundary condition is satisfied. The subsurface is first filled with a uniform cellCon of 0.01 S/m; then, an edgeCon is assigned to the edges representing the well path from the surface to Z = −1000 m. According to Heagy and Oldenburg (2019), the casing has an outer diameter of 0.1 m, a thickness of 0.01 m, and a conductivity of 5 × 106 S/m. So, the edgeCon of the casing can be calculated as the cross-sectional area times the intrinsic conductivity.

Four sets of source excitations are set up for the return electrode offset 2000 m, 750 m, 500 m, and 250 m. To visualize the electric field on the surface, we deploy receiver stations over a uniform data grid in the area of X = [−1250, 1250] and Y = [−1250, 1250] at a 20 m spacing. Each station measures the Ex and Ey fields using a pair of x-oriented and y-oriented M-N electrodes, respectively. The same receiver array is repeated for the four sets of sources.

The script runs for approximately 4 s on the same computer. The coefficient matrix is only factorized once because the conductivity models remain unchanged for the four simulations, and the solutions are obtained by quick back-substitution for different right sides (sources). The calculated potential difference data are first normalized by the electrode spacing, and then the total strength is computed using the Ex and Ey components at every station. The results of RESnet are in good agreement with the top row in Figure 5 of Heagy and Oldenburg (2019). For the return electrode offset of 2000 m, the field strength is almost symmetrically single-peaked at the wellhead location; but the finite-distance effect of the return electrode becomes more evident at receivers of X < −1000 m (Figure 7a). Other short-offset cases feature a double-peaked pattern in the field strength when the return electrode is within the area of investigation (Figure 7b7d).

It is important to note that the field around the return electrode appears stronger than around the wellhead, as the casing, effectively acting as a long electrode, spreads the current source over the entire length of the well, resulting in a reduced current density on the surface near the casing. The same phenomena have been observed and carefully examined in the previous section. The results of this modeling exercise suggest that the return electrode should be placed as far away as possible to maximize the magnitude of measured anomaly signals in a top-casing source survey. The original example in Heagy and Oldenburg (2019) is modeled using a 3D cylindrical-mesh finite-volume method, which simulates the exact geometry of a hollow casing with a small cell size of 0.25 cm. The SimPEG codes and the script used in their example are open source (Cockett et al., 2015); readers can refer to their paper for details. This example of casing demonstrates the effectiveness, efficiency, and convenience of RESnet in modeling thin and conductive objects of practical significance.

Complex infrastructure

In the third example, we show a modeling scenario involving complex metallic infrastructure that cannot be easily handled by any other algorithms. We design a synthetic ERT survey using a standard dipole-dipole array for the detection of subsurface contamination in an industrial environment with a lot of metallic infrastructure (Figure 8). The measured data, plotted as the apparent resistivity on pseudosections, are examined to study the effect of manmade infrastructure on ERT data.

The survey site is within a core mesh area defined by X = [−45, 45] and Y = [−12, 12] and discretized by the smallest cell size of 1 m. The top of the vertical mesh grid is at Z = 4 m, although the surface is still assumed to be at Z = 0 m because we must model buildings and pipelines above the ground. The dipole-dipole array uses 21 electrodes in the range of X = [−40, 40] at a constant 4 m spacing along the x-direction (the white dots in Figure 8). There are 18 sets of source excitations in each simulation.

Four different electric models are simulated to demonstrate the effect of complex infrastructure. The first is a uniform half-space, implemented by filling the region of Z > 0 m with a very resistive cell conductivity of 10−5 S/m, and Z < 0 m with a typical earth conductivity of 10−2 S/m. The script for the first model runs for approximately 3.7 s on the same computer. The resultant pseudosection shows a consistent uniform pattern in apparent resistivity within a 1% error compared to the true value of 100 Ω·m, except for the shortest source-receiver separation with an error not exceeding 5% (Figure 9a). Anomaly-free modeling is another test of numerical accuracy and the overall correctness in computation.

In the second model, two anomalous blocks, one conductive (0.1 S/m) and the other resistive (0.001 S/m), are added using cellCon to the uniform subsurface as the target of the ERT survey (the blue and yellow blocks in Figure 8). The blocks vary in size, shape, and depth of burial and are distributed at different sections of the survey line. This model runs for about the same time as the first model. The resultant pseudosection clearly shows the well-known anomalous data patterns associated with the deep conductive block and the shallow resistive block (Figure 9b).

The third model adds a vertical steel casing and a metallic warehouse to the model. The casing is 100 m long and is located 6 m away from the survey line (the vertical red line in Figure 8); the casing is modeled by the same edgeCon parameters as in the previous subsection and is not connected to any other object. The warehouse is 10 m long, 6 m wide, and 6 m tall and is built with five pieces of steel sheets 5 mm thick for the four side walls and the roof; the warehouse walls also have a 2 m section buried in the earth for better mechanical stability (the red box in Figure 8). The warehouse is modeled by assigning faceCon parameters to the mesh faces representing the warehouse walls and roof. Because the addition of nonzero edgeCon and faceCon does not increase the number of circuit nodes, the computing time remains about the same. On the resultant pseudosection, the indication of the two sought blocks is still visible, but there also is another significant presence of anomalies from the well casing and warehouse. The interaction among the targets, the infrastructure, and the measurement system has become complicated because that resistive anomaly in apparent resistivity cannot be straightforwardly explained by the metallic warehouse. Such a distorted data pattern may confuse practitioners in the field, and our algorithm provides a handy tool for a quantitative understanding of the effect of infrastructure.

In the fourth model, we connect the head of steel casing to the southwestern corner of the warehouse using a pipeline suspended 4 m above the surface (the horizontal red line in Figure 8). Such an object is difficult for conventional methods to simulate if not impossible. In RESnet, this type of above-ground pipeline is implemented by indexing and connecting the corresponding nodes in the equivalent resistor network, regardless of how far the nodes are apart. The newly established branch has a total conductance calculated using the same concept of edge conductivity. The fourth model runs, again, for approximately 3.7 s because adding a new branch to the circuit only slightly alters the sparsity pattern of the coefficient matrix, and the number of unknowns (nodes) remains the same. The simulation result shows a pseudosection drastically different from the one from the third model (compare Figure 9d with 9c), although the only change to the model is a pipeline 4 m above the surface, which often is ignored in practice. This simulation demonstrates the importance of considering seemingly irrelevant infrastructure when carrying out ERT surveys, as well as the necessity of having advanced modeling tools to study the coupling relationship between the geologic targets and infrastructure.

The drastic change of pseudosections from Figure 9b9d can be surprising because some resistive features turn to conductive in the next simulation. We note that the minimum value of apparent resistivity decreases from Figure 9b9d, as a result of adding more and more conductive infrastructure. The pseudosection can be conceptually considered as a superposition of resistive and conductive anomalies all around the electrode array, so the presence of highly conductive objects significantly stretches the lower limit of the apparent resistivity, overshadows other moderate anomalies, and may eventually confuse practitioners in characterizing the resistive or conductive nature of the actual targets.

Using an equivalent circuit also has other benefits and conveniences not shown in our examples. Instead of assigning very resistive values to the mesh cells in the air, one can model the air by removing the nodes and resistors in the air from the circuit, leading to a numerically smaller and easier problem. Similarly, when extremely resistive objects, like plastic sheets, exist in the subsurface, one can model their effect by removing the corresponding branches and resistors, effectively creating open circuits. Some nodes and branches at the corners or boundaries of a modeling domain also can be safely dropped if the current there is deemed sufficiently small, depending on the geo-electric model and survey configuration.

Quantitative modeling of thin and conductive infrastructure is known to be numerically complicated because the spatial dimension of metallic objects is much smaller than the sizes of regular geologic units. We have developed software called RESnet based on the concept of approximating the subsurface conductivity structure with an equivalent resistor network. RESnet allows the definition of conductivity properties not only at the cell centers as in the conventional methods but also on mesh faces and edges. As a result, a 3D rectilinear mesh can contain three conductivity models — cellCon, faceCon, and edgeCon. Such innovation brings users great convenience and flexibility in dealing with thin metallic objects in 3D DC resistivity modeling.

In this paper, we have made the RESnet software in MATLAB and Python freely available to the public. We also have introduced the physical principles, mathematical tools, and computing devices involved in developing the software. A numerical study is carried out to demonstrate the convenience, efficiency, and accuracy of RESnet for the modeling of steel casings in comparison with the conventional conductive pillar approach using the regular FD method. In addition, three simulation examples are provided to demonstrate the working of RESnet in modeling common geologic problems, top-casing sources for reservoir monitoring, and ERT surveys in an industrial environment. In the last example, our result has shown that above-ground pipelines, when electrically coupled with other infrastructure, can significantly distort the measured data in an ERT survey. The source codes have been explained in such a detailed way that users in academia and industry can quickly adopt RESnet in their research or production.

The development of RESnet program opens many new opportunities. The new formulation of conductivity model that includes edge, face, and cell conductivity values makes it possible to recover the edge and face conductivities in the inversion, if thin ore bodies or metallic objects are the targets to be imaged. The concepts of edge and face conductivity also introduce a new way of modeling anisotropy without invoking the tensor conductivity because they are by definition directionally conductive. More applications of RESnet are to be explored after the public release of this software.

This work was supported by the National Natural Science Foundation of China (grant no. 41974087) and Guangdong Provincial Key Laboratory of Geophysical High-Resolution Imaging Technology (grant no. 2022B1212010002). The computational work was supported by the Center for Computational Science and Engineering at Southern University of Science and Technology. The author thanks S. Chen, K. Wang, L. Yang, Y. Li, and M. Ravasi for their help in the development, maintenance, and release of the RESnet software. The proofreading of A. Malcolm has improved the quality of this manuscript.

The entire package of the RESnet software, along with the examples presented in this paper, is open source and can be freely accessed through the MATLAB and Python versions (Yang, 2018b, 2023). Feedback and issues can be reported on the repository website.

APPENDIX A DESIGN OF SOFTWARE

RESnet follows a linear and straightforward flow of operation that consists of three major steps (green, purple, and blue boxes in Figure A-1).

  1. Problem setup: users specify the mesh, model, and survey (electrode arrangement).

  2. Solving: the original simulation task is transformed into an equivalent resistor network circuit problem and solved for electric potentials at mesh nodes.

  3. Postprocessing: users obtain the simulated potential difference data as demanded by the setup parameters and carry out further analysis or plotting.

For most applications, users can find it straightforward to accomplish their specific modeling goals by combining different components from the template scripts. In general, end users mostly interact with the problem setup and postprocessing, whereas the solving part has been standardized and does not expect changes from users.

RESnet is a stand-alone scientific program for the purpose of carrying out a 3D DC simulation in geophysics. The following summarizes our ideas of software engineering when developing the codes and can aid recycling of the codes for other purposes.

  1. Procedural-oriented programming (POP). Small and mid-sized software projects can benefit from the flexibility and agility of POP if the computation and operations are encapsulated in functions. Our Python implementation, RESnet-py, also is developed with the programming paradigm of POP. It is our intention that, in addition to people who use the software as a whole, some users may integrate the functionalities of RESnet in their programs. POP also is good at documenting algorithms with the least coding effort and learning curve.

  2. Minimum external dependency. Our codes are designed to run on computers with only the simplest installation of MATLAB and Python. The MATLAB version of RESnet does not require any add-on toolboxes. The Python version only imports from numpy, scipy, matplotlib, and ctypes; the latter is used to call MKL Pardiso from Python. Users may find RESnet lacks many supporting functionalities, for example, the mesh generation and data visualization. However, we have decided only to include the core functionalities of numerical simulation in the software because users have the freedom of using other utilities and integrating RESnet into their established routines.

  3. Templated user scripts. The user interface of RESnet is the command line or MATLAB/Python integrated development environment. We encourage users to develop their scripts by following the example templates. By using scripts to run the software, we find a trade-off between hiding tedious details and revealing the machinery of code because the software is intended to be used as a packaged production tool and a sandbox for research. Users also are supposed to develop graphic user interfaces and other postprocessing applications to work with RESnet if necessary.

APPENDIX B EXPLANATION OF CODES

RESnet comes with a minimum working example, RUNME.m, which solves a nominal 3D DC resistivity problem with a miniature mesh. This example involves all the necessary steps required by simulating an electric survey in RESnet and can be a good starting point to learn how the code works. We recommend that users run the test code RUNME.m first to ensure the software is installed correctly. The following description of the code is based on the MATLAB version RUNME.m. The Python version of RUNME.py is a one-to-one translation of RUNME.m. The components of the software as listed subsequently, along with the logical flow connecting them, are shown in Figure A-1.

Mesh design

RESnet adopts a right-handed coordinate system, in which the x-direction represents easting and the direction from left to right across the page, the y-direction represents northing and the direction from the front to the back into the paper, and the z-direction represents elevation and the direction from the bottom to the top.

In the current implementation, a 3D rectilinear mesh is used to describe the geo-electric model and to build the resistor network with its nodes and edges (branches). RESnet specifies a rectilinear mesh using three vectors: nodeX, nodeY, and nodeZ, which list the coordinates of nodes in the x, y, and z-direction. Note that the elements in nodeX and nodeY are in ascending order, but those in nodeZ are descending to accommodate the custom of counting from the top to bottom in geophysics. An open source repository rectmesh3d-m provides convenient utilities for mesh/model creation, manipulation, and visualization (Yang, 2018a).

Model design

In RESnet, the model design is to assign conductivity values to the cells, faces, and edges of a given mesh. RESnet uses the concept of blocks to specify the conductivity model. The first variable, “blkLoc,” defines the 3D ranges of all the intended blocks to be designated to certain positions in the mesh; each row of blkLoc defines an individual block by specifying the coordinate limits in the format of [xmin, xmax, ymin, ymax, zmax, zmin]. If the minimum and the maximum coordinate limits are identical for a particular dimension, that dimension vanishes, and a face conductivity or edge conductivity value is implied. The second variable, “blkCon,” is a vector storing cell, face, or edge conductivity values of the corresponding blocks in blkLoc. If a face or edge conductivity is intended, users are responsible for calculating the spatially integrated conductivity values in advance.

Survey design

In this step, users have to decide the electrode arrangement by completing the source parameter variable tx and the receiver parameter variable rx. In MATLAB, tx is a cell array; each element of tx defines a separate set of source excitation by specifying the location and current intensity of the current electrodes. Typically, the A and B electrodes, representing the positive and negative current electrodes, respectively, can be described as [Ax, Ay, Az, 1; Bx, By, Bz, −1], in which the subscript x, y, and z indicate the coordinates. However, RESnet is flexible in defining current sources with a wide distribution. In RUNME.m, the second set of source consists of three electrodes because the positive current source is split into two parts at two distant locations. RESnet requires that the total current intensity within the same set of source (tx’s cell array element) must sum to zero, or the computed result is erroneous.

The variable rx also is a cell array, whose number of elements is the same as that of tx. Each element of rx specifies all the requested potential electrode (M, N) locations under the excitation of the corresponding set of source; the information about the potential electrode (receiver) is organized in the format of [Mx, My, Mz, Nx, Ny, Nz] for each potential difference datum.

formRectMeshConnectivity

This function converts the topological structure of mesh to the connectivity parameters of an equivalent circuit. When a rectilinear mesh is given, all of its nodes are extracted to form a list called “nodes,” in which each row is the [x, y, z] coordinates of each individual node. By indexing the edges connecting the nodes, a list called “edges” is generated; each row of edges contains the index numbers of the starting and ending nodes. The vector “lengths” is for the length of each edge. The list of “faces” has all the mesh faces, and each row of faces is the index number of the edges surrounding that face. The vector “areas” is for the area of each face. The list of “cells” has all the mesh cells, and each row of cells is the index number of the faces surrounding that cell. The vector “volumes” is for the volume of each cell.

makeRectMeshModelBlocks

This function converts the information about the blocks in blkLoc and blkCon to three full-length conductivity model vectors cellCon, faceCon and edgeCon, whose values are defined at cell centers, faces and edges, respectively. The elements in cellCon are the intrinsic conductivity of each cell, and the counting order is from +z to −z, then from −x to +x, and finally from −y to +y. Faces (normal direction) and edges are directional objects, so faceCon and edgeCon vectors are composed of x, y, and z-oriented faces and edges, in that order. Within each orientation type, elements are ordered by following the rule same as in cellCon. The elements in cellCon, faceCon, and edgeCon are with units of S/m, S, and S·m, respectively. The calculation of face and edge conductivities is the user’s responsibility outside of the RESnet program. For the cells, faces and edges not covered by the blocks in blkLoc, the corresponding conductivity elements are assumed zero.

formXXXX2EdgeMatrix and total conductance

These three functions calculate the weighting and summing matrix Cell2Edge, Face2Edge, and Edge2Edge. When applied to cellCon, faceCon, or edgeCon, the matrix can transform the conductivity models on a mesh to the conductance on the branches of the equivalent circuit. Formulating the transformation in the format of matrix multiplication facilitates the derivation of sensitivity for the purpose of inversion. Multiplying Cell2Edge, Face2Edge, and Edge2Edge matrix with cellCon, faceCon, and edgeCon yields the conductance on branches Cc, Cf, and Ce from the contribution of cell, face, and edge conductivities, respectively. The total conductance on branches of the equivalent resistor network circuit is the sum of the three.

calcTrilinearInterpWeights

This utility function essentially calculates the weighting coefficients for the trilinear interpolation. Such interpolation is required in handling the sources and receivers that can be arbitrarily located in the 3D space. When a point source is specified as a row in the variable tx, the interpolation finds the surrounding eight nodes in a lattice grid structure and redistributes the intended current intensity to those nodes using the trilinear weights. This discretization of source is necessary because we can only apply current sources at the nodes of a circuit. The superposition of all the current injections specified in tx for a particular set of sources yields a vector of external current intensities sources at all the nodes, which acts as a right-hand side (RHS) in solving the circuit problem. When there are multiple sets of sources, the variable sources is a multicolumn matrix, and each column represents an individual set of sources.

Potential (M, N) electrodes may not necessarily sit at the nodes, so the same trilinear interpolation finds the values at the query points by the weighted averaging of calculated potentials at the surrounding eight nodes. The calculated weights are for the mapping from field points to the nodes, so by applying the potential vector to the transpose of the weighting matrix (Mw, Nw), the interpolated potentials at field points can be obtained.

solveRESnet and output

Because the equivalent resistor network has been established, solving the circuit problem requires only the variable edges as the connectivity information, C is the total conductance on all the branches and sources for the external current intensities at all the nodes. The function solveRESnet first creates a potential difference operator matrix that maps the potentials on the nodes to the potential differences on the edges (branches). Such a matrix is the circuit analog to the gradient operator in the FD method. Then, the potential differences are left-multiplied with a diagonal matrix of branch conductance to obtain the currents on the edges. Application of KCL requires calculating the net current flowing in and out of each node, and this operation turns out to be the transpose of the potential difference matrix with the sign flipped. As a result, we can build the linear system of equations with the difference matrix, the diagonal matrix of conductance, and the sources as the RHS. The solution of the linear system is obtained using direct solvers. In MATLAB, a Cholesky decomposition is used to factorize the coefficient matrix, followed by a back-slash operation that runs the back substitution for quick solutions of possibly multiple RHS. In Python, matrix factorization is implemented by calling the Pardiso solver in MKL through the wrapper code PyPardiso.py. Advanced users can replace the solvers with their own choices because the coefficient matrix and the RHS have been explicitly given.

In the end, solveRESnet outputs potentials for the electric potential data at all the nodes, and potentialDiffs and currents for the potential differences and currents, respectively, on all the branches. By applying the transpose of Mw and Nw to potentials, we obtain the electric potentials at the locations of M and N electrodes. Differentiating the potentials at M and N electrodes yields data requested by the receiver parameter variable rx.

A biography and photograph of the authors are not available.

Freely available online through the SEG open-access option.

Supplementary data