COMSOL 6.4 - A Geoelectrical Forward Problem

Electrical resistivity tomography (ERT) is a geophysical method used for imaging the resistivity of soil and rocks in the subsurface, based on measurements of electrical quantities by means of electrodes at the surface or inside boreholes. The method has its roots in the vertical electrical sounding techniques (VES) (Ref. 1) and exploits the fact that resistivity varies over several orders of magnitude in different types of rock. In particular, it depends on the conductivity and quantity of fluids and minerals contained in the rocks.

ERT is today routinely applied in order to investigate targets from diverse applications such as mining, oil and gas exploration, groundwater, waste disposal, landslides, dams and embankments, geological faults, volcanoes, and archaeology. The method is closely related to electrical impedance tomography (EIT) as it is used in medical imaging or nondestructive material testing. One peculiarity of geophysical ERT is, however, that the boundaries of the target region are typically far away from the measurement system, virtually at infinity.

Typically, the physical domain in which the geoelectrical survey takes place is so large that it can be considered unbounded to the sides and the bottom. This property must be properly taken into account during modeling.

A single geoelectrical measurement involves the injection of a current I through two point-like electrodes, typically called C1 and C2, and the measurements of the potential difference V12 between two other electrodes (P1 and P2) resulting in a measured resistance R = V12/I. This application describes the 3D ERT forward problem for 25 electrodes in a homogeneous ground of resistivity 100 Ω⋅m, and compares the result with the analytical solution.

Although technically a low-frequency alternating current is used, inductive and capacitive effects are typically ignored and the governing equation is

(1)

Here, V represents the electric potential, σ is the subsurface electric conductivity, and

are point current sources from the electrodes located at positions rC1 and rC2 (actually, one electrode acts as current source, and the other as current sink).

The goal of ERT is to determine the unknown subsurface resistivity from measurements with a multitude of electrode configurations, R(C1,C2,P1,P2)i, i = 1,…, N. The procedure typically involves two steps:

•

The prediction of measurement responses on a given model (forward step).

•

The correction of the model to achieve a better fit (inverse step).

The procedure is repeated iteratively until an acceptable fit between predicted and measured data, and possibly an agreement with additional constraints (for example, smoothness or a priori information), are found.

While the forward step is well defined by Equation 1, many algorithms have been proposed for the inverse step. Most commonly, gradient-based inverse methods are used (Ref. 2), which require the calculation of a parameter sensitivity matrix S that describes how the measured resistance would be affected by the change of a parameter (typically the resistivity in a certain subsurface region).

It can be shown (Ref. 3) that the measured resistance can also be expressed as

with the sensitivity function Si defined as

(2)

which is given by the reciprocity theorem (Ref. 4). Here, rC1 and rC2 designate the positions of the electrodes injecting the electric current into the ground, and rP1 and rP2 are the position for the electrodes measuring the difference in electric potential V12, for a given configuration “i”.

Model Definition

The model illustrates the forward step and the sensitivity calculation for a typical situation of near-surface ERT with 25 point electrodes spaced one meter from each other. The model is geometrically parameterized to be easily adaptable to various sizes.

At the sides and at the bottom of the computational region, the infinite element domains account for the fact that the domain is not bounded.

Figure 1: Computational domain. A line of 25 electrodes is placed at the top of a 50-by-50-by-20 meters box bounded by infinite element regions.

As a material, a homogeneous half-space that often serves as a reference model for forward algorithms is used. The electric conductivity is σ = 0.01 S/m, which corresponds to a resistivity of 100 Ωm.

The solution is compared to the analytical solution

(3)

To achieve sufficient accuracy, the mesh density is refined directly at the electrodes and in the region of particular interest underneath them. A swept mesh is used in infinite element regions; see Figure 2.

Figure 2: Refined mesh around the electrodes and swept mesh in infinite element regions.

Use the continuation solver to switch between different electrode configurations.

Results and Discussion

Figure 3 shows the potential distribution around the active pair for the second of the two electrode configurations that the model solves for.

Figure 4 compares the simulation results for the potential at the electrode positions for both electrode configurations to the analytical solution given by Equation 3.

Figure 5 shows the relative error of the computed solutions with respect to the analytical solutions given by Equation 3 for the two electrode configurations. The plot is taken along the line joining the point sources.

Figure 3: Electric potential around the active electrode pair.

Figure 4: Comparison between analytical (lines) and modeled (markers) results.

Figure 5: The relative error of the computed solutions with respect to the analytical solutions for the two electrode configurations. The plot is taken along the line joining the point sources.

The model solves for two current dipole situations with a common midpoint. Together they allow the calculation of the sensitivity of a Wenner-α configuration according to Equation 2.

S = with(1,ec.Jx)*with(2,ec.Jx) + with(1,ec.Jy)*with(2,ec.Jy) +

with(1,ec.Jz)*with(2,ec.Jz))

The with operator makes reference to the modeled configuration (1 or 2 in the first argument of the operator). Figure 3 shows a plot of this expression.

Figure 6: Sensitivity plot.

References

1. F. Wenner, “A Method of Measuring Earth Resistivity,” National Bureau of Standards, Bull. 12(4) 258, pp. 478–496, 1915.

2. M.H. Loke and R.D. Barker, “Practical techniques for 3D resistivity surveys and data inversion,” Geophysical Prospecting, vol. 44, pp. 499–523, 1996.

3. S. Friedel, “Resolution, Stability and Efficiency of Resistivity Tomography Estimated from a Generalized Inverse Approach,” Geophys. J. Int., vol. 153, pp. 305–316, 2003.

4. D.B. Geselowitz, “An application of electrocardiographic lead theory to impedance plethysmography,” IEEE Trans. Biomed. Eng. BME-18, pp. 38–41, 1971.

ACDC_Module/Devices,_Resistive/geoelectrics

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Geometry 1

Add parameters that are useful for defining physical and geometric quantities.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
L	50[m]	50 m	Domain length
W	50[m]	50 m	Domain width
H	20[m]	20 m	Domain height
WI	4[m]	4 m	Infinite layer thickness
N	25	25	Number of electrodes
a	1[m]	1 m	Electrode separation
x0	L/2-N/2*a	12.5 m	First electrode x-position
y0	W/2	25 m	First electrode y-position
rho0	100[m/S]	100 Ω·m	Resistivity

Add two extra parameters which are going to be useful for sweeping among excitations applied to different entities. For further details, look at coming comments in these instructions, where the parameters are used.

In the table, enter the following settings:


Name	Expression	Value	Description
dom_C1	39	39	Positive electrode domain number
dom_C2	44	44	Negative electrode domain number

Geometry 1

Block 1 (blk1)

In the toolbar, click

In the window for , locate the section.

In the text field, type L.

In the text field, type W.

In the text field, type H.

Locate the section. In the text field, type -H.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	WI

Find the subsection. Select the checkbox.

Select the checkbox.

Click

Polygon 1 (pol1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type range(x0,a,x0+N*a).

In the text field, type y0+0*range(0,1,N).

In the text field, type 0*range(0,1,N).

Click

Block 2 (blk2)

In the toolbar, click

In the window for , locate the section.

In the text field, type (N+4)*a.

In the text field, type N*a/3.

Locate the section. In the text field, type x0-2*a.

In the text field, type y0-(N+4)*a/2.

In the text field, type -N*a/3.

Click

Form Union (fin)

In the window, click .

In the window for , click

Click the

button in the toolbar.

Click the

button in the toolbar to see the interior of the geometry.

The completed geometry is shown in Figure 1.

Click the

button in the toolbar to return to the default state.

Definitions

Analytic 1 (an1)

In the toolbar, click

In the window for , type V_ref in the text field.

Locate the section. In the text field, type 1[A]*rho0/(2*pi)*(1/abs(x-x1)-1/abs(x-x2)).

In the text field, type x, x1, x2.

The analytical solution given by V_ref will be used to validate the computed solution.

Ground Boundaries

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select all the exterior boundaries located in the ground. This operation can be performed quickly by selecting the checkbox, then clicking on one boundary from each of the sides and the bottom of the box.

In the text field, type Ground Boundaries.

Infinite Element Domains

In the toolbar, click

Select Domains 1–9, 11, 12, and 14–19 only.

In the window for , type Infinite Element Domains in the text field.

Infinite Elements Boundaries

In the toolbar, click

Select Domains 1–8, 11, 12, and 14–19 only.

In the window for , locate the section.

From the list, choose .

Select the checkbox.

In the text field, type Infinite Elements Boundaries.

Infinite Element Domain 1 (ie1)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Set the physical size of the infinite elements domain to be 100 times the size of the geometry. Larger sizes can give even more accurate results but will require a more refined mesh to be resolved.

Locate the section. In the text field, type 1e2*dGeomChar.

Materials

100 Ohmm Homogeneous

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	1/rho0	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic

In the text field, type 100 Ohmm Homogeneous.

Electric Currents (ec)

Point Current Source 1

In the toolbar, click

and choose .

Select Points 29–54 only.

In the window for , locate the section.

In the Qj,p text field, type 1*(dom==dom_C1)-1*(dom==dom_C2).

The dom variable is a built-in variable that represents the domain number. It uniquely identifies each entity within the specified dimension. In this model, it is used to select different points, since on point 1 it has a value of 1, on point 2 it has a value of 2, and so on. The == operator, together with a sweep over a parameter, allows for scanning over different sources, identifying them by their domain number. Point number dom_C1 injects a unit charge. Point number dom_C2 extracts a unit charge. See the COMSOL Multiphysics Reference Manual for more information about built-in variables and operators.

Locate the section. Click