Pageoph, Vol. 112, 1974

Birkh~iuser Verlag, Basel

A Numerical Method to Calculate the Electromagnetic Field of a Horizontal Current Dipole ~) By H. SCRIBA2)

Summary - A compuler program is described to calculate the field of an oscillating eurrem dipole over a horizontally stralified medium for the quasi-static case. Multi-layer master curve~ fearfrequency stranding can be calculaled wilb Ibis program. Examples of these curves are presented 1. Inn'oducgior Various electromagnetic controlled-source methods are currently used (KELLER [ [ ]). Compared to DC methods, they offer the advantage of a better signal-to-noise ratio (e.g, FR6Hr_tr [2], 8CRmA [3]) and the possibility of deep sounding with fixed electrode separation and variable frequency (frequency sounding). Considerable work has been done on the evaluation of the field of a vertical magnetic dipole (see e.gJ VAr~YAN [5]). In this case the vector potential has only one component A* and one obtains relatively simple expressions for the field components. The field of a horizontal current dipole is more complicated, since the vector potential has not only a primary component A~ (dipole oriented along x-axis) but also a secondary component Az. The only field component depending on A~ alone is Bz (VA~Y~,~ et aL [4]), all other components depend on Az too. The purpose of this study is the description of a procedure to calculate the field of a horizontal current dipole.

2. Fundamental equations Let u~ consider a homogeneous isotropic conducting medium. We shall use the Giorgi system of unils in our calculations. We start from Maxwell's equations

V•

(1)

v • ,q- b =3

(2)

v. B = o

(3)

v. z3=,,

(4)

t) Contribution No. 90 of the Geophysical Institute, Swiss Federal Institute of Technology, Ziirich. ~) Geophysical Institute of the Swiss Federal Institute of Technology, P.O.B. 266, CH-8049 Ziirich, Switzerland.

802

H. Scriba

(Pageoph,

where E represents the electric field, /3 the dielectric displacement, H the magnetic field,/~ the magnetic induction, J the electric current density and P the electric charge density. The material equations are 6 = u~

(5)

/5 = eE

(6)

3= ~

(7)

where e is the dielectric constant (5 v a c u u m ~" 8 0 ) , /2 is the permeability (/~ vacuum= #o) and a is the electric conductivity. The field vectors can be derived from the electromagnetic potentials A and U = v x ~

(8)

and f: = - J

- v u.

(9)

Let us assume that all quantities are harmonic functions of time, proportional to e ''~ where ~ois the angular frequency. By splitting offthe factoir e u~ we obtain functions of ~o alone. We shall use capital letters for functions of time and small letters for functions of frequency. Observing the Lorentz convention

~2 V. a + - - u = 0

i~o

(10)

where = ~ / - i ~ 2 + i#ao~

(l 1)

is the complex wave number, we get from equation (9) i(o = --ico~ - ~'i V(V" ~).

(12)

Thus all the field vectors can be derived from the vector potential ~ which fulfils the inhomogeneous Helmholtz equation v 2 ~ _ ~2~ = _ _ ~ f

(13)

This equation must be solved for the boundary conditions given, 3. Solution of the boundary value problem 3.1. Theory We assume that the earth is horizontally stratified and has N infinitely extended isotropic layers. We let the surface of the earth coincide with the plane z = 0 (see Fig. I). The layers have been numbered starting at the top. The uppermost layer is the upper half-space (air = vacuum) and has the index 0. Let the dielectric constant of the layer

Vol. 112, 1974)

The Electromagnetic Field of a Horizontal Current Dipole

V

E~,Bx

Io'dV

~

803

EzV,Bz

//////////////////////-//////////////////////~

layer

1 Pl,dl

layer j pj,dj

layer

N-1 Pu-l,dN-1

layer

N ON,dN==r

Z I

Figure 1 Horizontally stratified medium with N layers. The horizontal current dipole is placed at the origin of the coordinate system

j be e~ and the permeability be #j. To a good approximation we can put #j = #o for allj. We assume that the displacement currents are small compared to the conduction currents and can be neglected (quasi-static approximation). This is no serious limitation except for electromagnetic prospecting at very shallow depths (cf. VANYAN [5]). Let the thickness of the l a y e r j be dj with the conductivity trj. The lowest layer is a halfspace. The boundary conditions may all be satisfied if the vector potential has only x and z components = (ax, az).

(14)

The tangential components of ~"and ~ must be continuous at the boundaries, from which follows the continuity of a~, a,, Oax/Oz and (1/~2)V.~ (SOMMERFELD [6]). A solution for the vector potential at the earth's surface is given by VANYAN [5]"

ax -

Iodlpo f me -mh 2--~ m + nl/rl Jo(mr)dm

a~

Iodlpox f me -"h 2~ r m + nl/rlJl(mr)dm

(15)

oo

0 oo

(16)

804

H. Scriba

(Pageoph,

where Io is the amplitude of the current in the dipole, dl the length of the dipole, h the height of the dipole above the earth, Jo and Jx the Bessel function of order zero and one, m is the integration parameter and n~ = q - m - ~ : ' ~ .

The rl and fl take into account the layered structure and can be calculated recursively 1 - k l e -2nldl r1

1 + k l e -2nla' n 2 - - n I r2

k I ~-~

n 2 + n l r2

(18) 1 -- kN_ 1 e -2n~-~dN-~

1 + kN_1 e - 2 n N - l d l v - I nN -- nN-l kN_ 1

nN W nN-i

~1 e-2ntdl 1 + f~l e - 2 m a l 1

-

k i = n2P2 - n l P l n2 P2 + nl pl

(19) 1 - kN-1 e -2nN-laN-1

fN-1 = 1 + ION-1 e -2nN-td~v-I nspN - ns-lps-1 feN-1

nNpN + nN-lPN-a "

The coefficients kj and ~ may be interpreted as the reflection coefficients for a plane wave with polarization perpendicular and parallel, respectively, to the plane of incidence (WAIT [7]). The field components are given by equations (8) and (12). Because of the complicated structure of the kernel functions, an analytical solution of the integrals can be found only for very simple cases such as the homogeneous half-space (WAIT [8], BANNISTER [9]). For a layered medium the integrals have to be evaluated numerically. 3.2. C o m p u t e r p r o g r a m 3.2.1. Infinitesimal dipole. Numerical solutions of the integrals (15), (16) and (17) exist only if the kernel functions converge to zero with increasing m. For h = 0 this is

Vol. 112, 1974)

The ElectromagneticField of a HorizontalCurrent Dipole

805

no problem and the integrals can be directly evaluated by stepwise integration. For h = 0 convergence of the kernels can be achieved by subtracting their asymptotic values. This gives a good rate of convergence for small values of 7- The integrals with the asymptotic kernels can generally be integrated analytically (SCRmA [3]). For large values of 7 the rate of convergence is improved by splitting up the field into two parts, one produced by a homogeneous half-space with the wave number 71, and a residual field caused by the layered structure of the medium (VANYAN[5])

/__?o +7l.

o0)

The fields at the surface of a homogeneous half-space are well known (FOSTER[10], WAIT [8], BANNISTER [9]) and thus we obtain for the field components at the surface of a stratified medium

ex = -cko

mT~rJo(mr)dm +

m[T~ - T~r]Jo(mr)dm

~0

+

0

r2-2x 2 _1 r2-3y 2 r-----q~ [Td - T#,lJl(rnr)dm 7• rS +

+ 71r e_,,,] ra ]

(21)

0

e, =--cko -xy 7

m[rl

-

[rl - rSlJdmr)dm

T ~ , ] J o ( m r ) d m - --7-

0

~'~ r' ) (22)

0

e, = -cio9 - -

nl T~ Jl(mr) d m - -~ 11

K1

(23)

0

I;

2xy toni T~ Jo(mr)dm + - ra

bx = r - ~ 2

0

by = c

;

nl T~ Jdmr)dm

0

{(x2 ); -~ - 1

ran1 T~ Jo(mr) dm +

r22;nlJimr,m rs

0

11 T

+;;

0

K1 T

(3Y2- x2)

72ry2(Io(~)Kl(~)-ll(~)Ko(~))]}

bz=c

m2T~,Jl(mr)dm + - ~ [ 3 - ( 3 0

(25)

+ 37, r + 7~r2)e-,,q

(26)

806

H. Scriba

(Pageoph,

where lz and K~ are the modified Bessel functions of order I and Tt =

1 -

-

mr I

r#, =

T~ =

1

h- nl

m + nl

rl

1

mr1 + n~

m + n~

p~

1(1) - I

(27)

(28)

(29)

/o d/~o C ~

2~

Figure 2 shows the flow chart of the corresponding computer program. The listing will not be given here but it is available on request. After reading the data, the program evaluates the integrals of equations (21) to (26). This is done by several subprograms: KERN computes the complex kernel functions. The rj and fj are calculated recursively according to equations (18) and (19). Bessel simultaneously computes the Bessel functions Jo and Jt by polynomial approximations (ALLEN [11]). The results are combined in FCTS to give the integrand. The integration is carried out by ROMBIN, a library subprogram that uses the Romberg method of integration (RoMBERG [12]). ROMBIN permits the simultaneous integration of several integrands. The integration procedure is controlled by INTEG. The length of the integration step is usually 1/r. The integrals are summed up until the required precision is reached. The program has been written in ALGOL-60 except for the subprogram KERN which is written in FORTRAN IV, a language that allows a simple formulation of complex calculus. The expressions for the homogeneous half-space are computed separately and the results are fed to the program as data.

3.2.2. Finhe dipole. For the computation of the field of a finite dipote, we have restricted ourselves to the y-component of the electric field so far, because this is the only component which is independent of frequency for the homogeneous half-space. If the length of the dipole is 2.a and is assumed short compared to the wavelength, the current Io can be assumed constant over the whole length of the dipole. The integration is then straightforward and we get

eyti.

The

Io l~oie~ [ [ y ; 2~

o

( Tl _ T~,)jl(rnr) drn ] i+i + [[ y_ ]~X +j ax_ .

integral in this equation contains the same kernels as the integral in equation (22).

(30)

Vo1.112,1974)

The Electromagnetic Field of a Horizontal Current Dipole

807

? INPUT OF X,Y,R, D(J),RHO(J), EXO,EYO,EZO, BXO,BYO,BZO

COMPUTATION OF GAM(J)

t INTEG: COMPUTATION OF THE INTEGRALS (21) TO (27)

[

~

ROMBIN: S IMULTANEOUS INTEGRATI01'I

t

COMPUTATION OF RE,IM,AMP,FI

OUTPUT OF X,Y, R,D(J),RHO(J), EX,EY,EZ,BX,BY, BZ

ICOMPUTATION OF ~COMPUTATION OF [THE INTEGRAND _J ]THE COMPLEX ~ |KERNEL FUNCTION B~SSEL: " L (FORTRAN) COMPUTATION OF [Jo AND J,

I

Is:: -s] NO

Figure 2 Simplified flow chart of the computer program for the horizontal N-layer case 3.3. E x a m p l e s

Some examples of the calculated fields at the surface of a horizontally stratified medium will be shown with the dipole situated at the origin and oriented in the xdirection (Fig. 1). The fields are presented as functions o f d ~ I. v -~/2 and are normalized such that f x j , v=o = 1

where v is the frequency. Thus the curves become independent of the resistivity and thickness of the first layer. Some results have been compared with CECCHINI and ANDRIEUX [13] and show good agreement. Also, the calculations of this paper appear to be consistent with those of ANDERSON (private communication). The first example (Fig. 3) shows the amplitude of the y-component of the electric field of an infinitesimal dipole for a three-layer case. The resistivities of the upper and

808

H. Seriba

(Pageoph,

eyl

5

10

3 1 0.3

1

0.1

0.1

i

i

i

10-5

10-4

10-3

i

i

[m-lsl/2]

lJii

10-2

d~1 v-V2

Figure 3 A m p l i t u d e of the y - c o m p o n e n t of the electric field (er) o f a n infinitesimal dipole at the surface of a threelayer earth. T h e resistivity ratios are 1 :x : 1, where x is the p a r a m e t e r indicated at the curves. T h e thickness ratios are 1 : 3 : ~o with r/d~ = 30 a n d ~0= 45 ~

lower layer have been normalized to 1 whereas the resistivity of the second layer is variable. The distance-to-thickness ratios are r/d1 --- 30 and r/d2 = 10 and the angle between the dipole and the radius vector is 45 ~ For dla/~ > 3.104 m.s ~/2 we obtain N--1 the high-frequency asymptote which depends on layer 1 alone. For ~J=l d j v ~ < 103 m.s ~/2 we have the D C value of the apparent resistivity. For P2 < P3 we get a clear

10

3 1 0.3

1

0.1

0.1

I

10-4

I

!

I

I

I Ill

I

10-3

I

I

I

I

I

I II

10-2

I

I

I

[m-Is1/2]

I

I

I I II

10-1

d~'l v-1/2

Figure 4 A m p l i t u d e o f the y - c o m p o n e n t o f the electric field (e,) of a finite dipole at the surface o f a four-layer earth. T h e resistivity ratios are 1 : 3 : x = I, where x is the p a r a m e t e r indicated at the curves. T h e thickness ratios are 1 : 3 : 1 0 : ~ with a/d~ = 100, r/a = 1/3 a n d r = 45 ~ T h e length of the dipole is 2 " a

Vol. 112, 1974)

The Electromagnetic Field of a Horizontal Current Dipole

809

extremum o f the field function whereas for P2 >/03 the extremum gradually disappears and the third layer is hidden by the screening effect o f the high resistivity o f the second layer. In Fig. 4 the amplitude o f ey is shown for a finite dipole over a four-layer earth. The thickness ratios are dl :d2:da = 3: 10:30. The length o f the dipole is 2 . a = 20-d3 with r/a = 1/3. The angle between the dipole and the radius vector is 45 ~ The resistivity ratios in this case are pl:p2:pa:p4 = 1 : 3 : x : 1 where x is variable. Again we get pronounced extrema for the different layers.

Acknowledgement

I would like to thank Prof. St. MUELLER for his continuous interest and encouragement and Dr. E. WIELANDT for reading the manuscript. I also wish to thank Drs. A. CECCmM and P. ANDRmUX (Compagnie G6n6rale de G6ophysique, Massy, France) for a copy o f their paper presented at the 1973 E A E G meeting, and Dr. W. L. ANDERSON (US Geological Survey, Denver) for exchanging c o m p u t e d examples. REFERENCES [1] G. V. KELLER,Natural-fieM and controlled-source methods in electromagnetic exploration, Geoexploration 9 (1971), 99-147. [2] R. K. FROHLICH, Geoelektrisehe Tiefensondierungen mit grossen Eindriagtiefen bei hohen Stfrfeldern, Dissertation, Technische Universitfit Clausthal, 1966, 79 pp. [3] H. SCRmA,Geoelektrische Widerstandsmessung mit tieffrequenten Rechteckstr6men, Dissertation No. 5180, ETH Ziirich 1973, 62 pp. [4] L. L. VANYAN,G. M. MOROZOVA,V. L. LOSHENITZINA,E. I. TEREKHINand A. I. SHTIMMER, Four-layer master curvesfor frequency electromagnetic sounding, in L. L. VANYANet al., Electromagnetic Depth Soundings (Moscow 1964), selected and translated from the Russian by G. V. KELLER,Consultants Bureau, New York, 1967. [5] L. L. VANYAN,Fundamentals of Electromagnetic sounding, in L. L. VANYANet al., Electromagnetic Depth Soundings (Moscow 1965), selected and translated from the Russian by G. V. KELLER, Consultants Bureau, New York, 1967. [6l A. N. SOMMERFELD,Ueber die Ausbreitung der Wellen in der drahtlosen Telegraphic, Ann. Phys. 4 (1926), Folge, 81, 1135-1153. 17] J. R. WAIT,Fields of a horizontal dipole over a stratified anisotropic half-space, IEEE Trans. Ant. Prop., AP 14 (1966), 790-792. [8] J. R. WAIT, The electromagnetic fields of a horizontal dipole in the presence of a conducting half: space, Can. J. Phys. 39 (1961), 1017-1028. /9] P. R. BANNISTER,Quasi-static fields of dipole antennas at the earth's surface, Radio Science 1 (1966), 1321-1330. [10l R. M. FOSTER,Mutual impedance of grounded wires lying on the surface of the earth, Bell Syst. Tech. J. 10 (1931), 408-419. 111] E. E. ALLEN,Analytical approximations, M.T.A.C. 8 (1954), 240-241. [12] W. ROMBERG, Vereinfachte numerische Integration, Det Kong. Norske Videnskabers Selskab Forhandlinger, Trondheim, 28 (1955), Nr. 7. [13] A. CECCHINIand P. ANDRIEUX,The horizontal alternating electric dipole at the surface of a layered earth application to electromagnetic soundings, abstract of a paper presented at the 35th Meeting of the European Association of Exploration Geophysicists at Brighton (Great Britain), 1973, Geophys. Prospecting 21 (1973), 599. (Received 19th April 1974)