Scattering of electromagnetic plane waves by a buried vertical dike

Scattering of electromagnetic plane waves by a buried vertical dike

Lurimar S. Batista^I; Edson E. S. Sampaio^II

^ICentro de Ciências Formais e Tecnologia da Universidade Tiradentes Campus II, 49032-490 Aracaju, SE, Brasil

^IICentro de Pesquisa em Geofísica e Geologia da Universidade Federal da Bahia Campus Universitário de Ondina, Instituto de Geociências sala 213B 40170-290 Salvador, BA, Brasil

^{Correspondence}

ABSTRACT

The complete and exact solution of the scattering of a TE mode frequency domain electromagnetic plane wave by a vertical dike under a conductive overburden has been established. An integral representation composed of one-sided Fourier transforms describes the scattered electric field components in each one of the five media: air, overburden, dike, and the country rocks on both sides of the dike. The determination of the terms of the series that represents the spectral components of the Fourier integrals requires the numerical inversion of a sparse matrix, and the method of successive approaches. The zero-order term of the series representation for the spectral components of the overburden, for given values of the electrical and geometrical parameters of the model, has been computed. This result allowed to determine an approximate value of the variation of the electric field on the top of the overburden in the direction perpendicular to the strike of the dike. The results demonstrate the efficiency of this forward electromagnetic modeling, and are fundamental for the interpretation of VLF and Magnetotelluric data.

Key words: scattering, electromagnetic waves, vertical dike.

RESUMO

Estabelecemos a solução exata e completa do espalhamento de uma onda plana eletromagnética no domínio da freqüência e no modo elétrico transverso por um dique vertical soterrado por uma camada condutora. Uma representação integral composta de transformadas unilaterais de Fourier descreve os componentes do campo elétrico espalhado em cada um dos cinco meios: ar, cobertura, dique e as rochas encaixantes de cada lado do dique. A determinação dos termos da série que representa os componentes espectrais das integrais de Fourier requer a inversão numérica de uma matriz esparsa e o método das aproximações sucessivas. Calculamos o termo de ordem zero da série para os componentes espectrais da camada de cobertura, para valores especificados dos parâmetros geométricos e elétricos do modelo. Este resultado permitiu determinar um valor aproximado da variação do campo elétrico no contato entre o ar e a camada de cobertura em uma direção perpendicular ao traço do dique. Os resultados demonstram a eficiência desta modelagem eletromagnética direta, e são fundamentais para a interpretação de levantamentos geofísicos com os métodos VLF e Magnetotelúrico.

Palavras-chave: espalhamento, onda eletromagnética, dique vertical.

INTRODUCTION

Electromagnetic (EM) wave propagation and EM geophysical methods are employed in mineral, groundwater and petroleum exploration, in shallow geotechnical investigation, or in related subjects such as Global Positioning System (GPS). They are based on the distribution of the EM field components in the ground induced by natural or man-made EM source fields. The electromagnetic description of a medium similar to the earth's crust consists in the determination of the electromagnetic parameters of the medium from the knowledge of the distribution of the components of the EM field inside it. Maxwell's and the constitutive equations provide the means for relating those field components and the electromagnetic parameters. Among the electromagnetic parameters the conductivity is the most diagnostic for the rocks of the crust. The analysis of the overall distribution of the charges, the electric currents, and the field components is in general very difficult because of complex geological structures. Therefore there are many unsolved problems related to Maxwell's equation demanding advanced research.

The investigation of the scattering of EM plane waves caused by lateral variation of the properties of the rocks is fundamental for the success of geophysical exploration. (Sommerfeld 1896) and (Wiener and Hopf 1931) employed different techniques to solve the problem of the scattering of EM plane waves in a perfectly conductive half-plane. The need for investigating three-dimensional problems, with a relatively small computation time, led to the development of several algorithms of numerical modelling (Gupta et al. 1987), (Livelybrooks 1993), (Avdeev et al. 1997), and (Zhdanov et al. 1997). Some are flexible, and either compute the wave field at non-uniform sampling intervals, or evaluate the derivatives with unequal deggrees of precision. Others excite the system from different functions, or restrict the information at the boundaries. All of them compromise the accuracy of the final result.

Presently there are very few options available in the literature of EM models obtained analytically. Generally the analytical solutions are related to one-dimensional models of the earth or simple structures as hemispheres or outcropping faults. The exact analytical solution of the scattering of a TE mode EM plane wave by an outcropping vertical fault in the frequency domain was done by (Sampaio and Fokkema 1992). Subsequently (Sampaio and Popov 1997) developed the correspondent analytical solution in the time domain for the zero-order term.

We present the basic formulation and the analytical solution of the scattering of a TE mode EM plane wave in the frequency domain by the following earth model: a vertical dike between two quarter-spaces and covered by an overburden. A preliminary version of this investigation has been presented by (Batista and Sampaio 1999). The analytical tools developed here will also be useful as a check of the efficiency of the numerical techniques of forward and inverse modelling.

FORMULATION OF THE PROBLEM

The Wave Equation

The application of the concept of an EM plane wave to geophysical problems was proposed originally by (Tikhonov 1950) and by (Cagniard 1953) for the study of the magnetotelluric method. We will follow this concept for the configuration of the primary field.

Let a TE mode EM plane wave, which propagates in the positive z direction, be scattered by a geological structure consisting of a buried vertical dike as depicted in Figure 1. The horizontal layer thickness is h, its electrical conductivity is , and it is in contact along the plane z = 0 with an infinitely resistive half-space (z < 0). The vertical dike presents an electrical conductivity and it is laterally limited by two vertical half-planes x = - a and x = + a. The two quarter-spaces have, respectively, conductivity for x < - a and conductivity para x > + a. The horizontal layer is in contact with both the dike and the two quarter-spaces along the horizontal plane z = + h.

For the TE mode the electric vector is always along the y direction, and for the described model the total electrical conductivity presents a constant and finite value in each medium. So the problem consists in finding the solution for the two-dimensional homogeneous Helmholtz wave equation in each medium:

where:

• k_n = , is the propagation constant in each medium; is the angular frequency; and are the free-space values of, respectively, the magnetic permeability and the dielectric permittivity; and is the conductivity.

• E_n represents the y component of the electric field vector.

E_n = E_y(x, z,) = e_y(x, z, t)e-itdt.

Representation of the Fields

The solutions of the differential equation in each medium are as follows (Figure 2):

1. For - < z < 0:

+

(3)

+

(4)

2. For 0 < z < h:

+

(5)

+

(6)

+

(7)

3. For h < z < :

+

(8)

+

(9)

+

(10)

Where: E^I = e^-ik₀z represents the incident field for z < 0; E = R_{0, j}e^ik₀z, R_{0, j} are the free-space reflection coefficients; E = T_{0, j}e^-ik₁z, T_{0, j} are the transmission coefficients from the free-space into the horizontal layer; E = R_{1, j}e^ik₁(z-h), R_{1, j} are the reflection coefficients for the horizontal layer relative to medium 2, 3 and 4 respectively; E = T_{1, j}e^-ik_j(z-h), T_{1, j} are the transmission coefficients from the horizontal layer into medium 2, 3 e 4 respectively; and u_n = , represents the wave number in the transformed space.

Boundary Conditions

To solve the problem it is necessary to find the 24 spectral components,f_{0, 2}, g_{0, 2}, f_{0, 3A}, f_{0, 3B}, g_{0, 3}, f_{0, 4}, g_{0, 4}, f_{1, 2}, g_{1, 2A}, g_{1, 2B}, f_{1, 3A}, f_{1, 3B}, g_{1, 3A}, g_{1, 3B}, f_{1, 4}, g_{1, 4A}, g_{1, 4B}, f₂, g₂, f_3A, f_3B, g₃, f₄, and g₄, employing the pertinent boundary conditions.

l_i

= l_m

and E_{i, j} = E_{m, j+1}, in the x direction;

l_i

= l_m

and E_{i, j} = E_{m, j+1}, in the z direction.

The weighing factor l_n is the inverse of impedance per unit length in each medium; = ±a; and = 0 or h.

Substituting Equations (2)-(10) in the boundary conditions and putting in the left hand side only the terms that are function either of cos(

x) or of cos(z), we obtain the following system of twenty four integral equations in twenty four unknowns - the spectral components:

(11)

[

(12)

(13)

[

(14)

(15)

[

(16)

(17)

[

(18)

[

(19)

[

(20)

[-

(21)

[

(22)

(23)

[

(24)

(25)

[

(26)

[

(27)

[

(28)

[-

(29)

[

(30)

[

(31)

[

(32)

[-

(33)

[

(34)

It will prove to be useful to organize this system in a matrix form. Therefore, Equations (11)-(34) can be rewritten in the following form:

(35)

Where: b= z, in Equations (11) to (22), or b= x, in Equations (23) to (34); () is a pentadiagonal square matrix ( 24×24), which contains the constant coefficients constituted by six sparse submatrixes,

=

j = 0, 1,

=

,

=

i = 2, 3, 4,

(;) is a sparse square matrix ( 24×24), formed by six sparse submatrixes,

=

,

=

,

=

,

=

.

Where: S^j = sin(j) and C^j = cos(j), j = a, h or 0; f () is a vector with 24 positions that represents the spectral components; and Y(b) is a sparse vector, which represents the independent constants.

Determination of the Spectral Components

Applying the inverse cosine Fourier transform

(

) in Equation (35), we obtain:

• Left Hand Side,

  

  (36)

• Right Hand Side,

  

  (37)

Where:

=

,

=

,

=

,

=

,

=

,

=

.

For

• S = sin(j);

• C = cos(j).

For j = a, h or 0 and n = 1, 2, 3 and 4.

() =

.

From Equations (36) and (37) we obtain,

(

(38)

Multiplying Equation (38) by W^-1(), we obtain,

(

(39)

Notice that () represents the inverse matrix, obtained by a straightforward operation; y () = () ^. y () represents the vector of the independent constants; and (;) = () ^. (;) is a sparse square matrix, denominated nucleus matrix.

The integral equations of the system (39) are classified as Fredholm singular integral equations of the second kind. The regular part,

(;)()d is a Riemann improper integral, possessing a finite value. So, the solution of the equation (39) can be obtained employing the method of successive approximations (Kondo 1991):

(40)

because

(

(41)

and

(42)

(43)

Where:

ⁿ(;) = ¹(;) o ^n-1(;), for n > 1, and

¹(;) = (;).

The elements of the vector () are the spectral components used to calculate the scattered electric field in the nine domains of the model under study.

It is sufficient the series (40) be convergent in order for

(44)

A necessary but not sufficient condition for the convergence of the series is that

() = 0. (Sampaio and Fokkema 1992) verified numerically that the terms of the series decrease and the series remains bounded for up to five terms for the model of a vertical fault, and we expect the model of the buried dike to present a similar behavior. However the rigorous proof of the convergence of the series is a difficult mathematical problem and we lack such a proof.

NUMERICAL RESULTS

Basic Concepts

The values of frequency and the values of the earth's crust conductivity employed in geophysics are such that the modulus of the constant of propagation of the subsurface is always much larger than the modulus of the constant of propagation of the air. In the present case this implies that | k_n| | k₀|, for n = 1, 2, 3, 4. Therefore the reflection and transmission coefficients assume the following values:

r_j = ;

s_j = ;

t_j = ;

where h is the overburden thickness.

Computation of the Spectral Components

Due to the fast convergence of the terms of the series, we will compute only the first term of the series (zero-order term), employing Equation (41). With the simplification of the series the original system reduces to a sparse system of 12 linear equations in 12 unknowns.

(

(45)

The solution of the system of equations (45) is relatively simple, because its matrix is sparse, pentadiagonal, and contains only three 4 × 4 independent submatrixes.

Graphical Representation of the Spectral Components

The convergence of the spectral components can be verified through a graphical analysis from the behaviour of each element of the vector f as a function of the integration variable . In order to do that we applied the zero-order solution to the following model of the earth: a vertical dike of width 2a = 0.5| k₃| and conductivity = 0.1 S/m, immersed between two quarter-spaces of conductivities = 0.01 S/m and = 0.005 S/m, and covered by a horizontal layer of conductivity = 0.05 S/m and thickness h = 0.01| k₃|.

Figure (3) displays the variation of the real and the imaginary parts of the spectral components on the surface of the earth. They are the kernels of the equations that represent the electric field on the top of the overburden.

Figure (3) displays the convergence of the spectral components relative to the integration variable. All the components of the vector y() have a similar behaviour: they start at a finite value for = 0, as increases they pass by a maximum or a minimum value, and rapidly decay to zero as + . This demonstrates the convergence of the sequence y(). Therefore the integrals converge and the computation is performed with a good accuracy.

Computation and Representation of the Normalized Electric Field

After the determination of the spectral components and the check of the convergence of the integrals we computed the normalized electric field, E , at the surface of the earth for the proposed dike model, employing Equations (5), (6) and (7). We employed different geoelectric parameters to simulate distinct models.

Figure (4) displays the variation of the real and the imaginary parts of the normalized electric field on z = 0 as a function of | k₃| x, for the vertical dike model with the following parameters: variable dike width 2a and overburden thickness h; = 0.05 S/m; = = 0.005 S/m; and = 0.1 S/m. The electric field has been computed for the following values of p₁ = 2a/h: 1; 5; 10; 50; and 200.

Figure (4) shows that as ( p₁

) the contacts between the dike and the surrounding rocks are very well defined at | x| = a. In the limit the model becomes that of an outcropping dike. On the other hand, as p₁ 0, which in the limit means either an infinite thickness of the overburden or a zero width of the dike, the same contacts become, as expected, progressively not well defined.

Figure (5) displays the variation of the real and the imaginary parts of the normalized electric field on z = 0 as a function of | k₃| x, for the vertical dike model with the following parameters: variable conductivity of the surrounding rocks and the dike; = 0.05 S/m; overburden thickness such that h = 0.02| k₃|; and width of the dike such that 2a = 0.4| k₃|. The electric field has been computed for the following values of p₂ = /: 0.1; 1; 10; 100; and 1000.

We observe in Figure (5) that, as expected, the behaviour of the curves for p₂ = 0.1 is the opposite to the curves for p₂ > 1. The real part curves for p₂ > 1 reach a maximum value at x = 0 and have an inflection point at | x| = a. For p₂ = 0.1 these curves reach a minimum at x = 0 and they also have an inflection point at | x| = a. The behaviour of the imaginary curves is exactly the opposed.

Figure (6) displays the variation of the real and the imaginary parts of the normalized electric field on z = 0 as a function of | k₃| x, for the vertical dike model with the following parameters: variable conductivity of the overburden and the dike; overburden thickness such that h = 0.02| k₃|; dike width such that 2a = 0.4| k₃|; and = = 0.0005 S/m. The electric field has been computed for the following values of p₃ = /: 0.1; 0.4; 1; 10; and 100.

We observe that in the real part curves the contacts between the dike and the surrounding rocks are well defined for p₃ > 1, whereas for an overburden more conductive than the dike they are not. However this distinction does not happen in the imaginary part curves. For them the contacts are always well defined. The maximum value is observed at x = 0 for p₃ 1, whereas for p₃ < 0.1 the position of the maximum is not well defined.

CONCLUSION

The complete and exact algebraic solution of the scattering of a monochromatic EM plane wave, for the case of a vertical dike immersed in two quarter-spaces and overlaid by a horizontal layer, was determined. As a first step, zero-order terms of the series representation of the spectral components were selected to compute an approximate value of the electric field above the vertical dike.

The results of the modeling show that: 1) for the outcropping dike, h 0, the contacts between the dike and the surrounding rocks are very well defined; 2) for a large thickness of the overburden, h , the contacts between the dike and the quarter-spaces are not defined; 3) when the conductivity of the overburden is much larger than the conductivity of the dike, > > , it is not possible to define the contacts between the dike and the surrounding rocks employing the real component of the electric field, but the contacts are defined through the imaginary component.

The results show the precision of the zero-order terms in the calculation of the secondary electric field. So they can be used in the computation to substitute other techniques with advantage. Therefore the expressions of the analytical solution of the electric field are fundamental for the interpretation of magnetotelluric, or VLF (Very Low Frequency) data, associated to the geophysical exploration.

ACKNOWLEDGMENTS

We acknowledge our fellowships and the grant from CNPq (Brazilian National Science Foundation). We thank the fruitful comments from Prof. M. Popov. We also thank the help from Mr. J. Lago in the preparation of the manuscript.

Manuscript received on May 3, 2002; accepted for publication on November 11, 2002;

presented by DIOGENES A. CAMPOS

• AVDEEV DB, KUVSHINOV AV, PANKRATOV OV AND NEWMAN GA. 1997. High-performance three-dimensional electromagnetic modelling using modified Neumann series. Wide-band numerical solution and examples. J Geomagn Geoeletr 49: 1519-1539.
• BATISTA LS AND SAMPAIO EES. 1999. Scattering of Electromagnetic Plane Waves by a Buried Vertical Dike. In: 6^th International Congress of the Sociedade Brasileira de Geofísica, Rio de Janeiro. Expanded Abstracts in CD-ROM, Rio de Janeiro.
• CAGNIARD L. 1953. Basic Theory of the Magneto-Telluric Method of Geophysical Prospecting. Geophysics 18: 605-635.
• GUPTA PK, BENNETT LA AND RAICHE AP. 1987. Hybrid calculations of the three-dimensional electromagnetic response of buried conductors. Geophysics 52: 301-306.
• KONDO J. 1991. Integral Equations. Oxford applied mathematics and computing science series. Oxford: Oxford University Press, 440 p.
• LIVELYBROOKS D. 1993. Program 3D FEEM: a multidimensional electromagnetic finite element model. Geophys J Int 114: 443-458.
• SAMPAIO EES AND FOKKEMA JT. 1992. Scattering of Monochromatic Acoustic and Electromagnetic Plane Waves By Two Quarter Spaces. J Geoph Res 97: 1953-1963.
• SAMPAIO EES AND POPOV MM. 1997. Zero-Order Time Domain Scattering of Electromagnetic Plane Waves by Two Quarter Spaces. Radio Science 32: 305-315.
• SOMMERFELD A. 1896. Mathematical Theory of Diffraction. Math Ann 47: 317-374.
• TIKHONOV AN. 1950. Determination of the Electrical Characteristics of the Deep Strata of the Earth's Crust. Doklady 73: 295-297.
• WIENER N AND HOPF E. 1931. One Class of Singular Integrals. Proc Prussian Acad Math Phys Sec 696-796.
• ZHDANOV MS, VARENTSOV IM, WEAVER JT, GOLUBEV NG AND KRYLOV VA. 1997. Methods for modeling electromagnetic fields: Results from COMMEMI. The international project on the Comparison of Modeling Methods for Electromagnetic Induction. Appl Geophys 37: 133-271.

+

+

(1)

Publication Dates

• Publication in this collection
  31 July 2003
• Date of issue
  June 2003

History

• Received
  03 May 2002
• Accepted
  11 Nov 2002