DCIP3D manual: 
Updates at version 2.0

1.0 Introduction

DCIP3D Version 1.0 has received considerable use since its release in 1998. The code requires that the electrodes be on the mesh nodes and this generally posed little difficulty in flat terrains. However, in regions of significant topography, the field electrode geometry can be highly irregular and the meshes require a large number of cells, many of which have a high aspect ratio. The primary motivation for Version 2.0 has been to overcome this problem and to provide an inversion algorithm that can easily accommodate any field geometry. Other improvements, suggested by Placer Dome, Falconbridge, and Newmont, and arising from our own experiences with the code, have also been made. In particular some numerical aspects have been improved to increase efficiency, models at intermediate points in the inversion are saved, potential problems with signs of input voltages are flagged, and an easy procedure to include surface smoothing is incorporated.

The following material outlines the major changes and provides examples of the improvements. This document is an addendum to the main document that outlines technical material for forward modelling and inversion for DCIP3D Version 1.0.

2.0 Forward modelling

For forward modelling the DC resistivity equations, DCIP3D uses a nodal-based finite volume technique in which the current is input on a node. This is illustrated in Figure 1.

Figure 1: Mesh design with DCIP3D. a.) Version 1.0 required electrodes be placed on cell nodes. b.) Version 2.0 allows the electrodes to be placed anywhere.

When inverting field data, the usual procedure is to generate a mesh whose nodes coincide with the location of the current electrodes. The difficulty with accomplishing this task is illustrated in Figure 1a. The left panel is an attempt to design a mesh that permits each electrode to be on a node. The number of cells required to accomplish this is large and the aspect ratios are undesirable. High aspect ratios of cells reduces the numerical accuracy and also reduces the speed at which the forward modelling equations can be solved. This problem is greatly exacerbated when field surveys are carried out in regions of considerable topography. It would be preferable to have a uniform gridding in which the cell size is dictated by the resolving power of the data rather than by small details regarding exact placement of electrodes. A desired grid is shown in Figure 1b.

Figure 2: Current electrode can be placed at an arbitrary position (xs, ys, zs) within a cell, or on a node. The currents are distributed to each node of the cell through linear-interpolation. (a) Current at arbitrary position (b) Current on a cell node

To handle a current electrode that is at an arbitrary position (x_s, y_s, z_s) in the cell we have two choices: (1) alter the structure of DCIP3D to change from a nodal based to cellcentered solution (this was the approach taken for EH3D), (2) keep DCIP3D intact but make a modification to distribute any current amongst the 8 nodes of the cell. The latter approach is shown in Figure 2, where a current I is distributed onto nodes P1 through P8. Effectively, we write

      (1a)

where r_s = (x_s, y_s, z_s) is the position of the current electrode, r_i = (x_i, y_i, z_i) is the position of the i′th node, and w_i is the linear-interpolation weighting for the i′th node.

Figure 3: Forward modelled apparent resistivities of a halfspace showing differences between placing the current electrode at cell centers versus cell nodal points. (a) Current electrode at cell center (b) Current electrode on a cell node.

      (2)

so the total current distributed among the 8 nodes is equal to I.With the linear interpolation we note that if the source electrode is on one of the faces of the cell, then only 4 nodes will be activated. If the source electrode is along an edge then the current is distributed between two nodes, and if the source electrode is at a corner of the prism, then only one node is activated. If potential electrodes are not on the nodes then fields are linearly interpolated.

As an example we model a halfspace in which the cell size is 50m and a current is injected at the surface in the center of the cell. Potentials are obtained on a 25m grid. Apparent resistivities should equal 100Ohm−m, which is the true halfspace value. The results are shown in Figure 3a. Errors up to 25% are observed at locations that are 25m (1/2 cell) from each of the four corners where the distributed current is input. At distances 50m (one cell width) the error has dropped to about 7%. These are expected results and are in accordance with testing using Version 1.0. In Figure 3b the apparent resistivities for a current on a nodal location are plotted. At a distance of 50m from the current, the error is less than 5% if we’re on a node. The errors increase somewhat between distances of one and two cells. We conclude that for numerical accuracies of about 5% or less, the observation should be at least one cell width away from the location of a current electrode.

Figure 4: Pyramid model used for forward modeling and inversion result. (a) Pyramid model mesh. (b) Pyramid model cross-section.

As a further example of forward modelling, we consider the pyramid model which was one of the test examples used for Version 1.0. Figure 4 shows a conductive prism buried beneath a topographic high. A pole-dipole survey with 50m dipoles was simulated and the apparent resistivity for a line data across the center of the pyramid is shown in Figure 5. Figure 5a shows the result with current electrodes in the centers of cells. It is to be compared with 5b where the currents are at nodes. The two pseudosections are very similar. Differences occur primarily because the experiments are not identical. The pole-dipole data in Figure 5a correspond to n = 1.5,6.5 while in Figure 5b the data are n = 1,6. The two pseudosections should be slightly shifted vertically with respect to each other. This is observed.

Figure 5: Forward model of the pyramid example with both current and potential electrodes on cell nodes. (a) Forward model with electrodes on nodes (b) Forward model with electrodes on cell centers

With the above tests, and others that we have carried out, we conclude that the alterations to the forward modelling to distribute currents among the nodes associated with the cell in which the true current resides, is working properly. We next look at an example for the inverse problem.

2.1 Inversion with DCIP3D Version 2.0

The inversion methodology presented in Version 1.0 has not been altered. The important changes for Version 2.0 are concerned with optimizing numerical computations. The major items include: no longer having to store the sensitivity matrix dcinv3d.mtx, and optimizing the way in which reciprocity is used in developing the sensitivity matrix. Generally, if the same electrode location is used more than once in a data file, its solution is stored for future use. These changes have no impact upon the final solution obtained by the inversion algorithm. Discussion of these, and some additional minor changes, are incorporated in the revised manual. In what follows, we focus upon the major implications of having the electrodes off of the nodal positions.

The test example for inversion will be the pyramid example shown in Figure 4. This example has 15,548 cells, with the smallest cell dimension of (50,50,25). 10 lines of poledipole data were collected above the prism, with 50m dipoles, and n = 1,6. The currents were on the surface but at the mid-point of the cell. Gaussian random errors of 5% + 0.0001 were added to the data. The model, recovered after 15 iterations, fit the data to the expected misfit. A plot of the resistivities of the surface cells, and a cross-section of the recovered model, are shown in Figure 6.

Figure 6: Inversion of the pyramid model. Both current and potential electrodes are placed at cell centers. (a) Surface cells of pyramid inversion. (b) Cross section of recovered model.

To compare, we carry out the same inversion but with data generated by having currents at the nodal locations. Gaussian noise (5% + 0.0001) was added and the data inverted. The surface cells and cross-section through the model are shown in Figure 7.

Figure 7: Inversion of the pyramid model. Both current and potential electrodes are placed on cell nodes. (a) Surface cells of pyramid inversion. (b) Cross section of recovered model.

One characteristic of the inversion, irrespective of whether currents are on or off the nodes is that the model exhibits rough structure at the surface (or more generally, at current and potential electrode locations). The reason for this lies in the sensitivity values. The sensitivities for the DC resistivity (or IP problem) are highest at the location of an electrode (either a current or potential electrode). Changing the conductivity at a location closest to the electrode makes the greatest change in the measured datum and hence that is the preferential location to add conductivity. This effect is well understood and a number of GIF software programs have incorporated ways to ameliorate this "electrode" artefact.

In Figure 8a we plot the sensitivity due to a current at a nodal location and for a potential field value at a distance of 300m from the current. This is a rough function which exhibits sign changes at the location of both electrodes. In the inverse problem the conductivity perturbation is made up of a linear combination of these sensitivities. This is one of the fundamental reasons why the first layer of surface cells often has considerable "chatter". In Version 2.0 we distribute the true current amongst neighboring electrodes and hence there is a superposition of the sensitivities shown in Figure 8a. In the case that the true current is in the center of the cell the sensitivity due to the modelled distributed currents, is that shown in Figure 8b.

Figure 8: Sensitivities calculated by placing the current electrode at a cell node, and at a cell center. (a) Cell node. (b) Cell center.

One way to ameliorate the structure caused by electrode sensitivities is to introduce a weighting that penalizes structure very near to the electrode. Weighting schemes ([1], [2], [3]) can be used. However, because we anticipate that this weighting will be used on a regular basis, we have provided a utility program, make_wdat.exe, to generate W_x and W_y (see reference manual) which are the weighting matrices for the horizontal gradients. The first p layers are assigned weights by the user. In the following example, the top four layers of the model were weighted with values of (64,16,4,2). The inversion result is displayed in Figure 9.

Figure 9: Inversion of the pyramid model. Both current and potential electrodes are placed at cell centers. Horizontal smoothing in the top 4 layers of cells has been incorporated in the inversion. (a) Surface cells of pyramid inversion. (b) Cross-section of the recovered resistivity model.

We stress that the above weighting should be applied with care. A strong horizontal smoothing can often eliminate horizontal changes in the conductivity, even if the earth model truly had these. The inclusion, and details, of the surface weighting therefore involves subjective decisions by the user. Without weighting, the model at the surface can be quite rough (due to electrode effects). On the other hand, local geology is also sometimes quite rough and hence the true earth model should exhibit a large heterogeneity.

Under the assumption that the top layers of the earth should be laterally smooth, then we suggest the following procedure for designing the weighting: Given the electrode locations, forward model data using a halfspace resistivity. Add noise that is representative of the field data. Carry out an inversion beginning with an incorrect value of the halfspace. This provides insight into the severity of the inversion artifacts. In fact, this is a reasonable procedure even if electrode artifacts are not of concern, because this inversion will provide first order information about how well the data might constrain the earth model. Next, experiment with your surface weighting parameters (number of layers to be smoothed and the relative strengths of weighting for each layer of cells) and repeat the inversion. Ideally, the minimum amount of weighting should be sought so that when the algorithm is applied to field data, there is no excessive penalty for true surface structure.

As a final synthetic example we invert both surface and borehole data for the pyramid model. The previous data set was augmented with surface potentials along one line measured from a pole current in a borehole. Five current locations were used and the currents were not at cell nodes. In this case the true current is distributed over 8 nodes. The location of the drill hole is shown in Figure 10. A cross-section of the recovered model is shown in Figure 11. Not only has the inversion algorithm performed well, but the additional information provided by the borehole data has improved the result. The bottom and sides of the block are more clearly defined. For this example, the final target misfit was achieved in 9 iterations.

Figure 10: DCIP3D 2.0 Inversion. There were several pole-dipole lines done along the surface. Borehole data was also included by using 5 current electrode positions down the borehole, and a single line of dipoles collected over the center of the pyramid.

Figure 11: DCIP3D 2.0 Inversion. Surface pole-dipole and borehole data generated by 5 current electrode positions down the borehole, and a single line of dipoles over the center of the pyramid, were inverted.

The above work has illustrated the essential alteration of DCIP3D Version 2.0, and shown it’s effect on forward modelling and inversion of DC data. Similar conclusions apply to the IP data. A suite of test results for the pyramid, 5-prism, and a field data set is provided with the documentation for the manual. As mentioned in the introduction, there were other modifications incorporated into Version 2.0. These are presented in point form in the following Appendix. The items for which the user has some control, are further elaborated upon in the Technical Manual. Overall, we believe that the improvements in Version 2.0 will greatly increase the practical usability of the code.

3.0 Summary of Major Changes Incorporated into DCIP3D Version 2.0

The following is a summary of the changes made to DCIP3D in version 2.0.

• All input files that ran on the previous version of DCIP3D will work with this new version with no changes.
• Batch-files or scripts used to run DCIP3D will need minor changes to comply with changes described below.
• In addition there is a new utility for building weighting files that will help produce smoother models in the top rows of the mesh. This utility is described at the end of this file.
• The remainder of this file describes changes to components of the DCIP3D library
  of programs. Components of the library that are not mentioned have no changes.

Changes to input/output file formats and DCIP3D 2.0 operation.

1. Electrodes can be located anywhere and are not restricted to mesh node locations.
2. New input format for dcipf3d. This allows IP data to be computed by evaluating J_n where n is a chargeability model (in msec, mrad, etc). That is, data are simulated using the linearized representation of IP data. Instead of command line parameters, dcipf3d now requires an input file.

Usage: dcipf3d dcipf3d.inp

Format of dcipf3d.inp:


   ip          ! dc | ip | ipL
   mesh.txt    ! mesh file
   obs_loc.txt ! electrode locations file
   model.con   ! conductivity model
   model.chg   ! chargeability model
   null        ! topography | null
   0           ! output potentials 0|1
   1.e-5       ! tolerance (optional)

• ’dc’ and ’ip’ options on line 1 are the same as before.
• ’ipL’ is a new option that calculates the ip data by multiplying the sensitivity matrix by chargeability.
• When ’dc’ is chosen, the chargeability model line is ignored.
• When ’1’ is entered on line 7, the potentials for every cell are output in file potentials_???.txt where ??? is the transmitter number. This file can then be viewed using Meshtools3d.

3. Tolerance: dcipf3d, dcinv3d, and ipsen3d have an additional optional line at the end of the input file for tolerance. This value indicates how well the forward modelled system is solved. If this line is missing, a default of 1.e-5 will be used.
4. dcinv3d now saves amodel after each iteration. Themodels are ordered: dcinv3d_01.con, dcinv3d_02.con, dcinv3d_03.con, etc.
5. dcinv3d.out has a header line.
6. The input files can have either length scales or alpha values. When three values are entered, it is assumed the numbers are length scales: (L_e,L_n,L_z). When four values are entered, it is assumed the numbers are alpha values: (_s,_e,_n,_x). The .log file would have both sets of values printed out.
7. The dc data are checked to see if some datamight have a wrong sign. A file check_sign.txt is created. This is just a warning about the potential of an incorrect sign.
8. The predicted data files have an extra column for true data.
9. Improved line search when running dcinv3d with mode=1.

Internal Changes to DCIP3D 2.0:

1. All floating point arithmetic is done in double precision. More accurate results are obtained, however more memory is required.
2. In dcinv3d, when the sensitivities are calculated, the dcinv3d.mtx file is not created. This results in a big performance gain.
3. In dcinv3d, the diagonal of G^TG is calculated during the sensitivity calculation, not after. A big improvement in performance is gained.
4. Code got cleaned and reorganized. Large working arrays are only allocated and used when needed. This results in neater code and less memory required.

A new utility, make_wdat.exe, is included. It is used to make a w.dat file that has smoothing in the X and Y directions for the first few layers that underly the topography surface. This suppresses the tendency of the algorithm to make highly variable structure in these top layers of the model.

Usage: make_wdat make_wdat.inp

Format of make_wdat.inp:

 mesh.txt  ! mesh file
 topo.dat  ! topography file | null
 3         !# of layers where extra smoothing will be applied
 10. 5. 2. ! weighting values to place in top layers of ’w.dat’

In the above example, the top 3 layers below air will have extra smoothing in the X and Y directions. The first layer will have a smoothing value of 10, the second layer down will have weighting value 5, etc. The first "layer" is defined as the first cell below topography, the second "layer" is defined as the next cell down, etc. No weighting is applied to any cells that are in the air above topography. Weighting values in all other cells of the mesh will remain at 1.0.

Sensitivity Calculation: (used in dcinv3d, ipsen3d, dcipf3d with ipL option)

When calculating the sensitivity matrix G, the number of times a forward system must be solved is equal to the number of transmitters plus the number of receivers.
To speed up the process of calculating G, if the same electrode location appears more than once in the data file, its solution is stored in memory for future use. The program the maximum number of solution vectors that need to be stored in memory and this number is displayed on the screen when the program starts. For example:

Maximum number of solution vectors to be stored for
reciprocity calculation: 37
Memory needed: 27.56 Mb

If the memory needed is too large, the number of vectors to be stored can be changed in the last line of the input file. If -1 is entered, then there is no limit. If 10 is entered for the example above, the following is output:

Maximum number of solution vectors to be stored for
reciprocity calculation: 37
Memory needed: 27.56 Mb
Only 10 will actually be stored.
Memory needed: 7.45 Mb

To have the maximum number of vectors stored be as small as possible, make sure that all the identical electrode locations are close together in the data file.

References

[1] Jiuping Chen, Eldad Haber, and Douglas W. Oldenburg. Three-dimensional numerical modelling and inversion of magnetometric resistivity data. Geophys. J. Int, 149:679 – 697, 2002.

[2] Yaoguo Li and Douglas W. Oldenburg. 3-d inversion of magnetic data. Geophysics, 61:394–408, 1996.

[3] Yaoguo Li and Douglas W. Oldenburg. 3-d inversion of magnetic data. Geophysics, 63:109–119, 1998.