Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://geophys.geol.msu.ru/shevnin/publ/180.pdf Äàòà èçìåíåíèÿ: Wed Oct 17 15:56:41 2012 Äàòà èíäåêñèðîâàíèÿ: Sat Feb 2 23:48:58 2013 Êîäèðîâêà:

Geophysical Prospecting, 2006, 54, 651­661

Analytical solution for the electric potential in arbitrary anisotropic layered media applying the set of Hankel transforms of integer order
E. Pervago1 , A. Mousatov1 and V. Shevnin1
1

´ ´ Instituto Mexicano del Petroleo, Mexico

Received November 2004, revision accepted March 2006

ABSTRACT The analytical solution and algorithm for simulating the electric potential in an arbitrarily anisotropic multilayered medium produced by a point DC source is here proposed. The solution is presented as a combination of Hankel transforms of integer order and Fourier transforms based on the analytical recurrent equations obtained for the potential spectrum. For the conversion of the potential spectrum into the space domain, we have applied the algorithm of the Fast Fourier Transform for logarithmically spaced points. A comparison of the modelling results with the power-series solution for two-layered anisotropic structures demonstrated the high accuracy and computing-time efficiency of the method proposed. The results of the apparent-resistivity calculation for both traditional pole-pole and tensor arrays above three-layered sequence with an azimuthally anisotropic second layer are presented. The numerical simulations show that both arrays have the same sensitivity to the anisotropy parameters. This sensitivity depends significantly on the resistivity ratio between anisotropic and adjacent layers and increases for the models with a conductive second layer.

INTRODUCTION The geoelectrical study of heterogeneous anisotropic media provides important additional information about geological structures and petrophysical rock properties. The various kinds of resistivity anisotropy in layered formations are due to the presence of fractured systems or cross-bedded sediments. The determination of anisotropy parameters in heterogeneous anisotropic media is a complex problem that requires fast and accurate methods for simulating the electric potential in such an environment. The general method for solving the electromagnetic problem in an anisotropic layered medium uses 2D Fourier transforms (2D FT) in the boundary plane. Gurevich (1975) proposed this approach for simulating vertical electrical sounding on the surface of a half-space consisting of layers with arbitrarily orientated anisotropy axes. By applying 2D FTs, Le Masne and Vasseur (1981) computed the electromagnetic field of sources on the surface of a homogeneous conductive half-space with a horizontal anisotropy axis. The scheme for the numerical computation of the electric potential from a point source in a layered earth with arbitrary anisotropy was developed by Das (1995) and extended by Li (1996) for the case of arbitrary sources. The 2D FT approach was used by Li and Pedersen (1991) to determine the electromagnetic response of a horizontal electric dipole on the surface of an azimuthally anisotropic medium. The current density and magnetic field for a direct current source in a layered conductor with an arbitrary anisotropy was calculated by Yin and Weidelt (1999). Yin and Maurer (2001) considered electromagnetic induction in a layered earth with arbitrary anisotropy.



E-mail: epervago@imp.mx

C

2006 European Association of Geoscientists & Engineers

651


652 E. Pervago, A. Mousatov and V. Shevnin

Another approach to modelling an azimuthal resistivity sounding over a two-layered structure with arbitrarily orientated transverse anisotropy was applied by Pervago (1998) and Bolshakov et al. (1998). In this case, the solution of the equation for an azimuthally anisotropic medium was represented by a set of Hankel transforms of integer order. Such a solution corresponds to the angular harmonic decomposition of the azimuthal resistivity diagrams. The analysis of the harmonic coefficients allows the heterogeneity and anisotropy effects to be separated (Pervago et al. 2001). Additionally, by using the decomposition characteristics, the sensitivity of electrical surface arrays and well-logging sondes to anisotropy was estimated by Bolshakov et al. (1997) and Mousatov and Pervago (2002). Here, we propose an analytical solution and computation technique for calculating the electric potential of a point DC source in a medium composed of layers with three-axial, arbitrarily orientated anisotropy. The solution obtained is expressed by a set of Hankel transforms of integer order. We have found recurrent expressions for the coefficients of the electric-potential spectrum in each layer. In order to convert the potential spectrum to the space domain, we applied the algorithm of the Fast Fourier­Bessel Transform for logarithmically spaced points (FFTLog) designed by Hamilton (2000). For the solution and algorithm verifications, the calculation results were compared with the data published by Bolshakov et al. (1998). Using the technique developed, we presented the responses of the traditional azimuthal resistivity surveys and advanced tensor array (Mousatov et al. 2000, 2002, 2003) above a three-layered anisotropic model.

ELECTRIC POTENTIAL IN AN ANISOTROPIC LAYERED MEDIUM The model consists of parallel homogeneous layers with three-axial, arbitrarily orientated anisotropy. We assume that the z-axis of the Cartesian coordinate system is perpendicular to layer boundaries and is orientated downwards. The electrical conductivity of the mth layer is defined by the tensor m . In an anisotropic homogeneous region, the electric potential Um (x, y, z) produced by a point source of direct current I with coordinates (xA , xA , yA ) satisfies the divergence equation, div( mgradUm) = I (x - xA, y - yA, z - z A) (1)

The solution of an inhomogeneous equation is the sum of a general solution for a homogeneous equation and the source function. Applying the 2D Fourier transform in the horizontal xy-plane to the homogeneous equation corresponding to (1), we obtain the homogeneous equation in the spectral domain with spatial frequencies, k x and k y :
m zz

~ 2U m m + 2i kx xz + ky z2
- -

m yz

~ Um 2 - kx z

m xx

m + k2 yy + 2kx ky y

m xy

~ Um = 0

(2)

~ Here, the electric potential Um (x, y, z) and its spectrum Um(kx , ky , z) are the Fourier transform pair, 1 ~ U m(kx , ky , z) = 2 Um(x, y, z) = 1 2
- -

U (x, y, z)e-

i(kx x+ky y)

dxdy,

(3)

~ Um(kx , ky , z)ei(

kx x+ky y)

dkx dky

(4)

~ In order to find the solution of (2) in the form, U (z) = evz , we first obtain the quadratic equation for the parameter v,
m m zz v2 + 2i kx xz + ky m yz 2 v - kx m xx m + k2 yy + 2kx ky y m xy

= 0,

(5)

which has the following roots: v1
,m

=

1 + m zz

2 mm kx xx zz - m yz

m xz

2

mm m + 2kx ky xy zz - xz

m yz

mm + k2 yy zz - y

m yz

2

m -i kx xz + ky

,
m xz 2 mm m + 2kx ky xy zz - xz m yz mm + k2 yy zz - y m yz 2

v2

,m

=

1 - m zz

2 mm kx xx zz - m yz

m -i kx xz + ky

.

(6)

C

2006 European Association of Geoscientists & Engineers, Geophysical Prospecting, 54, 651­661


Electric potential in anisotropic layered media 653

Since v2,m(kx , ky ) = -v1,m(kx , ky ) where v is the complex conjugate of v, the common solution of (2) with complex coefficients can be written as the sum of refracted and reflected parts, ~ U m(kx , ky , z) = Am(kx , ky )evm z + Bm(kx , ky )e
-vm z



,

(7)

This solution can also be written in the cylindrical coordinate system (r, , z) as a series of Hankel (Fourier­Bessel) transforms of integer order. Transforming to cylindrical coordinates by substituting, x = r cos , y = r sin , kx = -kr sin k , ky = kr cos k ,

(4) can be rewritten in the form, Um(r, , z) = 1 2
- 0 n

~ kr U m(-kr sin k , kr cos k , z)e-

ikr r sin(k - )

dkr dk

(8)

Then, using the Jacobi­Anger expansion for Bessel functions of the first kind of integer order J e
it sin p

=

n=

Jn (t )einp ,

(9)

and introducing the modified spectrum, ~ ~ ~ U m(kr , k , z) = kr Um(-kr sin k , kr cos k , z), we can rewrite (8) as series of the Hankel transform, Um(r, , z) =


(10)

e
n=

in 0



Jn (kr , r )

1 2

-

~ ~ U m(kr , k , z)e

-ink

dk dkr .

(11)

Finally, this equation can be written in the compact form, Um(r, , z) =
n=- m Cn (r, z)e in

.

(12)

m The complex harmonic coefficient Cn (r, z) is defined by the Fourier and Hankel transforms with spatial frequencies k r , k as m Cn (r, z) = 0

Jn (kr , r )

1 2

-

~ ~ U m(kr , k , z)e

-ink

dk dkr .

(13)

The application of the Hankel transform allows us to avoid the singularity of the potential spectrum of a point DC source at spatial frequency zero, if the 2D Fourier method is used. Additionally, in this case, the potential spectrum is discrete and its harmonic coefficients decrease rapidly when the harmonic number increases. If the medium is isotropic, the modified spectrum does not depend on the spatial frequency k , i.e. ~ ~ ~ ~ U m(kr , k , z) = U m(kr , z),
m and the harmonic coefficients Cn (r, z) become equal to zero for all n = 0 due to the property of the Fourier integral,

1 2

-

~ ~ U m(kr , z)e

-ink

~ ~ dk = {U (kr , z)

n = 0,

0, n = 0.

Thus, the potential Um (r, z) is equal to the zeroth harmonic coefficient Cm (r, z) (equation (12), and 11 is converted into the 0 well-known solution in the form of the Hankel integral, Um(r, z) =
0

~ ~ J0 (kr r )U m(kr , z) dkr .

C

2006 European Association of Geoscientists & Engineers, Geophysical Prospecting, 54, 651­661


654 E. Pervago, A. Mousatov and V. Shevnin

~ ~ Combining (7) and (10), the modified spectrum function U m(kr , k , z) in the mth layer can be written as ~ ~ ~ U m(kr , k , z) = A(kr , k )e where wm(k ) = m - im,
wm(k ) = m + im, 2 2 wm kr z

~ + Bm(kr , k )e

-wm kr z

,

(14)

(15) (16)
m xz mm m sin2 k + zz xy - xz m yz mm sin 2k + yy zz - m yz

m(k ) =

mm xx zz -

cos2 k


m m xz sin k + yz cos k

m zz

,

(17)

m(k ) =



m zz

.

(18)

The vertical component of the current density vector j tensor m and the potential U m as j~z.m(x, y, z) =
m xz

z,m

in the mth layer can be expressed in terms of the conductivity

m m Um(x, y, z) + yz Um(x, y, z) + zz Um(x, y, z). x y z ~ Um(kx , ky , z). z

(19)

After a Fourier transformation of (19), the spectrum j~z,m(kx , ky , z) has the following form:
m m ~ j~z,m(kx , ky , z) = xz kx + yz ky Um(kx , ky , z) + m zz

(20)

Using the modified potential spectrum (10) in cylindrical coordinates, we obtain the Hankel spectrum for the vertical component of the current density tensor, ~ ~ m m ~ j~z,m(kr , k , z) = kr - xz sin k + yz cos k U m(kr , k , z) +
m zz

~ ~ U m(kr , k , z). z

(21)

~ ~ ~ Substituting (14) for U m(kr , k , z) into (21), the spectrum j~z,m(kr , k , z) can be written as ~ ~ j~z,m(kr , k , z) = kr cm Am(kr , k )e where
m cm(k ) = m(k ) + zz wm(k ), m dm(k ) = m(k ) - zz wm(k ), m m m(k ) = -xz sin k + yz cos k . wm kr z

~ + dm Bm(kr , k )e

-wm kr z

,

(22)

(23) (24) (25)

~ ~ The unknown coefficients, Am(kr , k ) and Bm(kr , k ), can be found using the boundary conditions, ~ ~ ~ ~ U m(kr , k , zm) = U
m+1

(kr , k , zm), (kr , k , zm).

(26) (27)

~ ~ j~z,m(kr , k , zm) = j~z,

m+1

Using the conditions far from the source point, i.e. ~ ~ U 0 (kr , k , z) 0, ~ ~ U N(kr , k , z) 0, z -, z +, (28) (29)

C

2006 European Association of Geoscientists & Engineers, Geophysical Prospecting, 54, 651­661


Electric potential in anisotropic layered media 655

~ ~ ~ as well as the U M (kr , k , z) continuity and the j~z, M (kr , k , z) discontinuity on the plane z A of a source placed in the Mth layer when tends to zero ( 0), we have ~ ~ ~ ~ ~ ~M ~M ~ U M (kr , k , z A + ) - U M (kr , k , z A - ) = U 0 (kr , k , z A + ) - U 0 (kr , k , z A - ), ~ ~z ~z ~ j~z, M (kr , k , z A + ) - j~z, M (kr , k , z A - ) = j~0, M (kr , k , z A + ) - j~0, M (kr , k , z A - ). (30) (31)

~ ~M The source potential spectrum U 0 (kr , k , z) and the spectrum of the vertical component of the source current density ~0 (k , k , z) are j~z, M r ~ ~M U 0 (kr , k , z) = and ~ j~0, M (kr , k , z) = z Ikr 4 ikr -
M

I e M 4 M zz

-kr (i M z+ M |z|)

(32)



M zz M xz

sin k +



M yz M zz

cos k

+ i M + M sign(z) e-

kr (i M z+ M |z|)

.

(33)

By substituting (32) and (33) into (30) and (31), we obtain the following conditions on the source plane: ~ ~ ~ ~ U M (kr , k , z A + ) - U M (kr , k , z A - ) = 0, Ikr ~ ~ . j~z, M (kr , k , z A + ) - j~z, M (kr , k , z A - ) = 2 (34) (35)

~ On solving this system of linear equations, we obtain the analytical recurrent equations for the coefficients, Am(kr , k ) and ~ m(kr , k ) (see Appendix). B ~ ~ After determining the quantities, Am(kr , k ) and Bm(kr , k ), we calculate the spectrum of the electric potential (14) and then the complex harmonic coefficients using (13). Because the internal integral in (13) corresponds to the Fourier transform,
m Dn (kr , z) =

1 2

-

~ ~ U m(kr , k , z) e

-ink

dk ,

(36)

it is convenient to perform its evaluation using the FFT method. The external integral of (12) represents the Hankel transform of integer order,
m Cn (r, z) = 0 m Dn (kr , z)Jn (kr r ) dkr .

(37)

In geophysical applications, the calculation of the Hankel integral is presented with a specific problem as a limited number of points has to cover a relative distance range of 104 ­105 . In the case of linearly spaced points (normal scale), the matrices of transforms are large (104 â 104 ­ 105 â 105 ). To avoid this problem, we applied the algorithm of the Fast Fourier Transform in logarithmically spaced points (FFTLog, developed by Hamilton (2000). Similarly to the normal FFT, the FFTlog method gives the exact Fourier transform of a discrete function that is periodic and uniformly spaced in logarithmic space (Talman 1978). The FFTLog calculates the fast Fourier and Hankel (Fourier­Bessel) transform of arbitrary order for a relative distance range of several orders with a limited number of points. The algorithm proposed by Hamilton (2000) converts the computation of the Hankel transform into two FFTs performed sequentially, and it improves the accuracy and speed of computation.

SIMULA TION EXAMPLE Based on the analytical solution and the calculation algorithms presented above, we simulated apparent-resistivity curves for the vertical electrical soundings in media composed of transversely anisotropic layers. These media are of interest for practical applications and correspond to geological structures formed by fractured rocks or intercalation of electrically contrasting thin beds. In this case, we assume that each layer is azimuthally anisotropic with the anisotropy axis perpendicular to the vertical

C

2006 European Association of Geoscientists & Engineers, Geophysical Prospecting, 54, 651­661


656 E. Pervago, A. Mousatov and V. Shevnin

axis z.In the local coordinate system that coincides with the anisotropy axes, the conductivity tensor of the mth layer is given by
m

m = 0 ^ 0

0 m 0

0 m

0

(38)

If the Cartesian coordinate system is used throughout the model, the conductivity tensor m can be written as m = RT mRm, m^ where cos m sin m cos m 0 0 is the rotation matrix, (40) (39)

Rm = - sin m 0

0 1

denotes the rotation angle around the z-axis. All simulations were performed for the relative distance of 1012 ­1013 with a total of 128­256 points for the electric potential spectra. To verify the technique described above, we calculated the harmonics of the apparent resistivity using our approach and the analytic solution given by power series obtained for the special case of the two-layered azimuthally anisotropic model (Bolshakov et al. 1998). In this model, the anisotropy axes in each layer are arbitrarily orientated in the boundary plane and are described by the conductivity tensor (39). Each layer m is characterized by the anisotropy coefficient m and the mean resistivity m = ( m m)-1/2 . The result of observations above an anisotropic medium is usually given as the apparent resistivity a , which ¯ can also be written as the sum of harmonics,
m n=- Cn (r, z)ein ,

a (r, , z) =

(41)

where C (r, z) is the complex amplitude of the nth harmonic. n It is convenient to compare the coefficients C (r, z) of the apparent-resistivity harmonics because the zeroth harmonic n amplitude is much higher than the amplitude of any other harmonic and a comparison of the apparent resistivities (harmonic sums) will not be relevant (Fig. 1). The C (r, z) graphs are plotted as a function of the normalized distance R/h where R is the n array spacing and h is the upper layer thickness. These graphs demonstrate good agreement between the two analytical solutions and indicate the high accuracy of simulation. The relative error in this case is less than 10-8 for all harmonics shown. We simulated the apparent resistivity a and apparent anisotropy coefficient a for a traditional azimuthal survey with a pole-pole array on the surface of three-layered structures for different model parameters and observation azimuths (Fig. 2). In this model, the second layer is transversely anisotropic with the anisotropy axis perpendicular to the vertical axis z (azimuthally anisotropic layer). The apparent-anisotropy coefficient a was calculated as the ratio of the apparent resistivities measured perpendicular to and along the anisotropy axis. When the mean resistivity of the anisotropic second layer is higher than the resistivity of adjacent layers ( 1 = 3 = 1, 2 = 20), the apparent-resistivity coefficient does not exceed the value 1.1, even for the true value 2 = 4 and thickness h 2 = 2h 1 . For the conductive anisotropic layer ( 1 = 3 = 20, 2 = 1), the apparent-anisotropy coefficient takes high enough values, even for a small relative thickness of the second layer, h 2 = 0.2h 1 .In this case, the behaviour of the a and a curves shows the feasibility of thin anisotropic-layer detection by performing measurements on various azimuths (the traditional technology of an array rotation around a reference point).

C

2006 European Association of Geoscientists & Engineers, Geophysical Prospecting, 54, 651­661


Electric potential in anisotropic layered media 657

Figure 1 Comparison of the coefficientsCn of the apparent-resistivity harmonics calculated for the two-layered azimuthally anisotropic structure by the power-series (solid line) and by the technique presented (circles). The curve index is the harmonic number. Model parameters are (A): = 1, 100 m, = 1, 2; (B): = 100, 1 m, = 1, 2. R/hours is the normalized spacing.



Figure 2 The apparent resistivity a (A) (B) and the apparent-anisotropy coefficients a (C) (D) for a pole-pole array above the three-layered structures with an azimuthally anisotropic second layer. Model parameters are (A) (C): = 1, 20, 1 m, h 1 = h; (B) (D): = 20, 1, 20 m, h 1 = h. 1: = 1, 2, 1, h 2 = 0.2 h; 2: = 1, 2, 1, h 2 = 2 h; 3: = 1, 4, 1, h 2 = 0.2 h; 4: = 1, 4, 1, h 2 = 2 h. AM/h is the normalized spacing.

C

2006 European Association of Geoscientists & Engineers, Geophysical Prospecting, 54, 651­661


658 E. Pervago, A. Mousatov and V. Shevnin

Figure 3 Scheme of the current and measuring electrode distribution for the tensor array, TEN

FD

.

Figure 4 The apparent mean resistivity a (A) (B) and apparent-anisotropy coefficients a (C) (D) for the tensor array, TENFD. , above three¯ layered structures with an azimuthally anisotropic second layer. The model parameters correspond to those given in Figure 2. AO/h is the normalized spacing.

Additionally, for the same models, we calculated the response of the tensor arrays, proposed by Mousatov et al. (2000, 2002, 2003), for the anisotropy parameter determination in heterogeneous media. This array (TENFD ) is based on the measurements of the first ( Ux = U 2 ­U 1 , Uy = U 3 ­U 4 ) and second ( Uyy = 2U 0 ­U 3 ­U 4 ) differences of the electric potential from a point source (Fig. 3). The tensor array provides the determination of all anisotropy parameters (mean resistivity, anisotropy coefficient and azimuth of anisotropy axis) by carrying out the observation in one arbitrary direction only, and it does not require the array rotation. The apparent mean resistivities and apparent-anisotropy coefficients a obtained by the TENFD (Fig. 4) are similar to the parameters measured by the conventional azimuthal survey.

CONCLUSIONS The analytical solution for media composed of arbitrarily anisotropic layers was obtained. Based on this solution, we developed an algorithm for the simulation of the electric field produced by a point DC source in an anisotropic multilayered medium.

C

2006 European Association of Geoscientists & Engineers, Geophysical Prospecting, 54, 651­661


Electric potential in anisotropic layered media 659

Comparison of the modelling results with the power series solution published for two-layered anisotropic structures demonstrated the high accuracy and computing-time efficiency of the method proposed. The analytical solution obtained is expressed as a sum of discrete angular harmonics with rapidly decreasing coefficients when the harmonic number increases. These harmonic coefficients are given by a set of Hankel transforms of integer order. Application of the Hankel transform allows us to avoid the singularity of the electric-potential spectrum at the spatial frequency zero that appears in the 2D Fourier transform. To calculate the Hankel transform we applied the fast Fourier-Bessel transform for logarithmically distributed points, thus providing high accuracy of calculation and significantly reducing the computation time. This calculation algorithm can be used to model the electric field for surface and borehole arrays, designed to determine the anisotropy parameters in heterogeneous media with plane or cylindrical layers. ACKNOWLEDGEMENTS ´ The authors are grateful to Instituto Mexicano del Petroleo, where this study was carried out within the framework of the Yacimientos Naturalmente Fracturados program. We thank Dr Pedro Anguiano Rojas for his valuable help and support. REFERENCES
Bolshakov D.K., Modin I.N., Pervago E.V. and Shevnin V.A. 1997. Separation of anisotropy and inhomogeneity influence by the spectral analysis of azimuthal resistivity diagrams. 3rd EEGS-ES Meeting, Aarhus, Denmark, Expanded Abstracts, 147­150. Bolshakov D.K., Modin I.N., Pervago E.V. and Shevnin V.A. 1998. Modelling and interpretation of azimuthal resistivity sounding over twolayered model with arbitrary-oriented anisotropy in each layer. 60th EAGE Conference, Leipzig, Germany, Extended Abstracts, P110. Das U.C. 1995. Direct current electric potential computation in an inhomogeneous and arbitrary anisotropic layered earth. Geophysical Prospecting 43, 417­432. Gurevich Y.M. 1975. On the theory of vertical electrical sounding in anisotropic media. Physics of the Earth 7, 102­105, (in Russian). Hamilton A.J.S. 2000. Uncorrelated modes of the nonlinear power spectrum. Monthly Notices of the Royal Astronomical Society 312, 257­284. Le Masne D. and Vasseur G. 1981. Electromagnetic field of sources at the surface of a homogeneous conducting half-space horizontal anisotropy: application to fissured media. Geophysical Prospecting 29, 803­821. Li P. 1996. A new solution for the electric potential in an arbitrary anisotropic and inhomogeneous layered earth. 66th SEG Meeting, Denver, USA, Expanded Abstracts, 962­965. Li X. and Pedersen L.B. 1991. The electromagnetic response of an azimuthally anisotropic half-space. Geophysics 56, 1462­1473. Mousatov A. and Pervago E. 2002. Azimuthally anisotropy estimation with electrical sondes applying harmonic decomposition. 72nd SEG Meeting, Salt Lake City, USA, Expanded Abstracts, 700­703. Mousatov A., Pervago E. and Shevnin V. 2000. New approach to resistivity anisotropic measurements. 70th SEG Meeting, Calgary, Canada, Expanded Abstracts, 1381­1384. Mousatov A., Pervago E. and Shevnin V. 2002. Anisotropy determination in heterogeneous media by tensor measurements of the electric field. 72nd SEG Meeting, Salt Lake City, USA, Expanded Abstracts, 1444­1447. Mousatov A., Pervago E. and Shevnin V. 2003. Arrays for tensor measurements of the electric field. Proceedings of SAGEEP, San Antonio, USA, Pp. 502­515. Pervago E. 1998. Anisotropy and inhomogeneity influence on the results of vertical electrical soundings. PhD Thesis, Moscow State University. Pervago E., Mousatov A. and Shevnin V. 2001. Joint influence of resistivity anisotropy and inhomogeneity on example of a single dipping interface between isotropic overburden and anisotropic basement. Proceedings of the SAGEEP, Denver, USA, ERP 7. Talman J.D. 1978. Numerical Fourier and Bessel transforms in logarithmic variables. Journal of Computational Physics 29 (1, October), 35­48. Yin C. and Maurer H.-M. 2001. Electromagnetic induction in layered earth with arbitrary anisotropy. Geophysics 66, 1405­1416. Yin C. and Weidelt P. 1999. Geoelectrical fields in layered earth with arbitrary anisotropy. Geophysics 64, 426­434.

APPENDIX For computing convenience (14) and (22) can be rewritten in the following form (Fig. 5): for the layers below the source point (z zA ,m = M,..., N): ~ + + ~m U + (kr , k , z) = bm am e
wm kr (z-zm )

+e

-wm kr (z-zm )

,

(A1)

C

2006 European Association of Geoscientists & Engineers, Geophysical Prospecting, 54, 651­661


660 E. Pervago, A. Mousatov and V. Shevnin

Figure 5 Scheme for calculating the coefficients a- , a+ , b- , b+ ,inan anisotropic mulmmmm tilayered medium.

~ J~

+ z,m

+ + (kr , k , z) = kr bm cmam ewm

kr (z-zm )

+ dme

-wm kr (z-zm )

,

(A2)

and for the layers above the source point (z zA ,m = 0, ... , M): ~ - - ~m U - (kr , k , z) = bm am e
-wm kr (z-zm-1 )

+e

wm kr (z-zm-1 )

,

(A3)

~ J~

- z,m

- - (kr , k , z) = kr bm dmam e-

wm kr (z-zm-1 )

+ cm e

wm kr (z-zm-1

),

(A4)

where
+ ~~ am = Am/ Bme + ~ bm = Bme 2Re(wm )kr zm

,

(A5) (A6) , (A7) (A8)

-wm kr zm

,

- ~ ~ am = Bm/ Ame - ~ bm = Amewm

-2Re(wm )kr zm-1

kr zm-1

.

In these equations, z m -1 and z m are the top and bottom of the mth layer. Using expressions (A1)­(A4), the system of equations (26­29) and (34,35) can be rewritten as follows:
- a 0 = 0,

(A9)

m = 1..,M,
- bm- 1

a

-wm-1 kr hm-1 - m-1 e

+e

wm-1 kr hm-1

- - = bm am + 1 ,

(A10) ,
wM kr (z A-z M-1 )

- bm-

1

dm-1 a

-wm-1 kr hm-1 - m-1 e

+c

m-1

ewm-1

kr hm-1

- - = bm dmam + c

m

(A11) = 0, (A12)

+ b+ a M e M

wM kr (z A-z M )

+e

-w kr (z A-z M ) M

- - b- a M e M

-w kr (z A-z M-1 ) M

+e

C

2006 European Association of Geoscientists & Engineers, Geophysical Prospecting, 54, 651­661


Electric potential in anisotropic layered media 661

+ b+ c M a M e M

wM kr (z A-z M )

+ dM e-

w kr (z A-z M ) M -w kr (z A-z M-1 ) M

- -b- dM a M e M

+ cMe

wM kr (z A-z M-1 )

=

I , 2

(A13)

and m = M, ... , N ­ 1,
+ + + bm am + 1 = bm+ 1

a

+ -wm+1 kr hm+1 m+1 e

+e

wm+1 kr hm+1

,
wm+1 kr hm+1

(A14) , (A15) (A16)

+ + + bm cmam + dm = bm+

1

cm+1 a

+ -wm+1 kr hm+1 m+1 e

+ dm+1 e

a+ = 0 N

where hm = zm - z m-1 is the mth layer thickness. Starting from the known coefficients, a- and a+ ,we can find the coefficients a+ , and a- , sequentially towards the Mth layer, m m N 0 which contains the source (Fig. 5):
- am =

(cm-1 - cm) + a (cm-1 - dm) + a

- m-1 - m-1

(dm-1 - cm)e-2Re(wm-1 )kr (dm-1 - dm)e-2Re(wm-1 )kr (dm - cm+1 )e-2Re(wm+1 )kr (cm - cm+1 )e-2Re(wm+1 )kr

hm-1 hm-1

,

m = 1,..., M,

(A17)

+ am =

(dm - dm+1 ) + a (cm - dm+1 ) + a

+ m+1 + m+1

hm+1 hm+1

,

m = N - 1,..., M.

(A18)

The coefficients, b- and b+ , in the Mth layer, which contains the source, can be found based on (A12) and (A13): M M b- = M
+ 1 + am e-2Re(wM )kr (zM -z A) I e +- 2 (c M - dM )(a M a M e-2Re(wM )kr h M - 1) - 1 + am e-2Rc(wM )kr (z A-zM-1 ) I e +- 2 (c M - dM )(a M a M e-2Rc(wM )kr h M - 1) -wM kr (z A-z M-1 )

,

(A19)

b+ = M

-w kr (z M -z A) M

.

(A20)

Finally, the coefficients, b+ and b- , are expressed through the coefficients previously obtained from the layers above and M M below the layer containing the source:
- - bm = bm+

1+a
1

a

- -wm kr hm me

- m+1

+e

wm kr hm

,

m = M - 1,..., 0,

(A21)

+ + bm = bm-

1+a
1

+ m-1

a

+ -wm kr hm me

+e

wm kr hm

,

m = M + 1,..., N.

(A22)

C

2006 European Association of Geoscientists & Engineers, Geophysical Prospecting, 54, 651­661