- Email: [email protected]
EMLCLLER—a program for computing the EM response of a large loop source over a layered earth model
ARTICLE IN PRESS
Computers & Geosciences 29 (2003) 1301–1307
EMLCLLER—a program for computing the EM response of a large loop source over a layered earth model$ N.P. Singh*, Toru Mogi Institute of Seismology and Volcanology, Hokkaido University, N10W8 Kita-ku, Sapporo 060-0810, Japan Received 15 August 2002; received in revised form 24 July 2003; accepted 5 August 2003
Abstract A simple program EMLCLLER coded in FORTRAN 77 is presented in this paper for modeling the electromagnetic (EM) response of a large circular loop source over a layered earth model in quasi-static as well as non-quasi-static (general frequency) regions. The program is based on a semi-analytical numerical approach for evaluating the improper integrals occurring in the expressions of EM field components. It takes into account both conduction as well as displacement current factors. The program is formulated in such a way that it is well suited for computing the EM response for any arbitrary position of the source loop in air or on the surface of the model, in contrast to the earlier methods which face convergence problem for source position on the surface of the model. The program has wide application and is capable of computing the EM response at any arbitrary point either inside or outside the source loop. The validity and accuracy of the program is demonstrated by computing the EM response of a large loop source over the homogeneous and multi-layer earth models. Response curves depict their characteristic variations. The computed results are in coincidence with the published results for the quasi-static region (fp50 kHz) and are extension of their characteristic variation in the non-quasi-static (50 kHzpfp1000 kHz) region. Matching of computed results with the published results demonstrate the validity of the program. r 2003 Elsevier Ltd. All rights reserved. Keywords: EM response; Large loop source; Layer earth model; Conduction current; Displacement current
1. Introduction Computation of theoretical geophysical response of a model is a routine step in the problems dealing with modeling and inversion schemes of geophysical data interpretation. The problem of computation of EM response over a layered earth structure due to various EM sources has attracted the attention of scientists since the early inception of EM prospecting methods. There $ Code available from http://www.iamg.org/CGEditor. index.htm. *Corresponding author. Department of Geophysics, Banaras Hindu University, Varanasi 221005, India. Tel.: +91-5422315878; fax: +91-542-2368174. E-mail addresses: [email protected], singhnp [email protected] (N.P. Singh), [email protected] (T. Mogi).
have been many studies related to the small loop of infinitesimal size (i.e. magnetic dipole) (e.g. Frischknecht, 1967; Vanyan, 1967; Ward, 1967; Wait, 1982; Kaufman and Keller, 1983; Ward and Hohmann, 1988; Kaufman and Eaton, 2001), whereas those related to the large loop of finite size are limited (Morrison et al., 1969; Ryu et al., 1970; Patra and Mallick, 1980; Poddar, 1982,1983; Ward and Hohmann, 1988; Kaufman and Eaton, 2001). Ryu et al. (1970) presented the EM response of a large circular loop source over a layered earth model on the basis of Morrison et al. (1969) formulation. Poddar (1982, 1983) has presented the EM response of a rectangular loop source over the layered earth models. In their studies, Ryu et al. (1970) were restricted in presenting the results for the quasi-static approximation pertaining to low frequencies, and only for the measurement points outside the source loop, whereas,
0098-3004/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2003.08.008
ARTICLE IN PRESS 1302
N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307
Poddar (1982, 1983) results lack for the higher frequency region, neglecting the displacement current and manually setting a non-zero (close to zero) value for the zero loop height, i.e. when source loop lies on the surface of the model. Neglect of the displacement current in the ground and/or in air may not be valid for certain geological environments and/or frequencies higher than 50 kHz: For example, with a multi-coil DIGHEMIV system, which employs three frequencies ð900; 7200 Hz and 56 kHzÞ; the neglect of displacement currents at the highest frequency ð56 kHzÞ becomes a problem when dealing with highly resistive regions. Further, in highly resistive environments, such as in certain regions of northern Canada, which are characterized by Precambrian igneous rocks overlain by electrically thin overburden and permafrost, the displacement current factor cannot be neglected (Fraser et al., 1992). Thus, it is pertinent to develop a suitable computational technique, which can compute, without limitations, the EM response comprising of both conduction as well as displacement currents over the wide frequency range and at any arbitrary receiver position either inside or outside the source loop. Therefore, in view of practical applications and filling the gap in existing literature, we developed a 1-D forward modeling program, EMLCLLER, which is complete and capable of computing the EM response of a large loop source over the layered earth model. Moreover, the development of such an algorithm is important and useful for the interpretation and inversion of large loop EM data, which may be preferable because of low transmitter power requirements, and convenience, economy and time saving in measuring over the entire ground inside as well as outside the source loop.
2. Theoretical considerations The geometry of the earth model, the loop source and co-ordinate system under consideration are shown in Fig. 1. The earth model comprises a sequence of nhomogeneous layers, and the co-ordinate system (x; y; z) is taken with its z-axis directed vertically downward. For taking the advantage of circular symmetry, cylindrical co-ordinate system ðr; f; zÞ has been also used frequently. The integral expressions of EM field components due to a finite horizontal circular loop of radius a; carrying a current Ieiot and placed at height z ¼ h above the surface of an n-layer earth can be given by integrating the transverse electric (TE) expression of potential at a point due to a vertical magnetic dipole, over the entire loop (Ward and Hohmann, 1988). The expression of various field components, namely Hz ; Hr and Ef at a measurement point on the surface
I(ω) ω
z= -h
aa
σ0, µ0, ε0 z= 0
( 0,0,0 )
z= z1 Y z= zi
h1
σ1, µ1, ε1
hi
σi,
µi,
X
εi
z= zi+1 z= zn-1
σn, µn, εn Z
Fig. 1. Earth model and measurement system under consideration.
of earth ðz ¼ 0Þ can be written as Hz ¼
Hr ¼
Ia 2 Ia 2
Z
N
½eu0 h ð1 þ rTE Þ 0
Z
l2 J1 ðlaÞJ0 ðlrÞ dl; u0
ð1Þ
N
½eu0 h ð1 rTE ÞlJ1 ðlaÞJ1 ðlrÞ dl
ð2Þ
0
and iom0 Ia Ef ¼ 2
Z
N
½eu0 z ð1 þ rTE Þ 0
l J1 ðlaÞJ1 ðlrÞ dl; ð3Þ u0
where rTE ¼ ðY0 Y# 1 Þ=ðY0 þ Y# 1 Þ with Y0 ¼ u0 =iom0 (intrinsic admittance of free space) and Y# 1 ¼ HyTE =ExTE ¼ ðHxTE =EyTE Þ (surface admittance at z ¼ 0). For an n-layer case, the surface admittance are given by the recurrence relation: Y# 2 þ Y1 tanhðu1 h1 Þ ; Y# 1 ¼ Y1 Y1 þ Y# 2 tanhðu1 h1 Þ Y# nþ1 þ Yn tanhðun hn Þ Y# n ¼ Yn Yn þ Y# nþ1 tanhðun hn Þ and Y# n ¼ Yn with Yn ¼ un =iomn ; un ¼ ðkx2 þ ky2 kn2 Þ1=2 ¼ ðl2 kn2 Þ1=2 ; and kn2 ¼ o2 mn en iomn sn : Here, r is the distance of the measurement point from the center of the loop. The harmonic time factor eiot is implied, and SI units are used everywhere. In numerical calculations, tanhðun hn Þ is used in its exponential form (Knight and Raiche, 1982) for stability reasons. The magnetic permeability for the rock formations in general is taken equal to the free space value, i.e. mi ði ¼ 1; 2; y; nÞ ¼ m0 :
ARTICLE IN PRESS N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307
3. Computational scheme The computation of improper integrals occurring in expressions of various EM field components, Eqs. (1)– (3), is complex because the product of two Bessel functions makes the integrand have an oscillatory nature, and this problem becomes more severe when h ¼ 0; i.e. when source loop lies on the surface of the earth model. Ryu et al. (1970), however, tried to overcome this problem by subtracting a known term inside the integral sign and adjusting the same or its equivalent analytical term outside the integral sign. But overall they only managed to adjust with this problem and accomplished their computation with certain constraints and limitations, such as (a) it is valid only for z ¼ h > 0; i.e. when source lies in air, and for the case when source lies on the surface of the model, they treated it as a limiting case, (b) it is valid for r > a; i.e. only for the points outside the loop. Moreover, in presentation of results, they were restricted only to the quasi-static approximations pertaining to the low frequencies (fp50 kHz) and for the receiver points outside the source loop. Therefore, in order to overcome these limitations, we have developed a new program using some special functional relationships and the Hankel transform method as described next. For computing the improper integrals involving the product of Bessel functions of orders 0 and 1 and occurring in the expressions of EM field components, these integrals are converted into the corresponding Hankel transforms of orders 0 and 1. Accordingly, the Hankel transform of the kernel KðlÞ of integer order n is defined as (Anderson, 1979; Guptasarma and Singh, 1997):
f ðrÞ ¼
Z
N
KðlÞJn ðrlÞ dl; for r > 0;
ð4Þ
0
where Jn is the Bessel function of first kind and of order n: The transform argument r > 0 is real but KðlÞ and therefore f ðrÞ may be complex functions of real variables. The kernel function used in the Hankel transform must be a continuous bounded function as the argument increases, ði:e: KðlÞ-0; as l-NÞ: Therefore, in order to use the Hankel transform method for evaluation of these integrals, convergence of integrals are tested. In case of slow convergence or divergent nature, they are made to converge rapidly by subtracting or adding a convergent integral expression under the integral sign and subsequently adjusting it or its equivalent analytic expression outside the integral sign. For example, the expression of Hz (given by Eq. (1)) is
1303
adjusted as follows: Z Ia N u0 h l2 Hz ¼ e ð1 þ rTE Þ l elh 2 0 u0 J1 ðlaÞJ0 ðlrÞ dl Z Ia N lh þ e lJ1 ðlaÞJ0 ðlrÞ dl: 2 0
ð5Þ
In Eq. (5), for h > 0 both integrals on the right-hand side are convergent, whereas for h ¼ 0; the first integral on the right-hand side remains convergent, whereas the second becomes divergent. Therefore, for the case when source lies in the air at a height, z ¼ h; the computation of Hz component given by Eq. (5), is carried out using the digital linear filter algorithm for the Hankel transform developed using the filter coefficients of Guptasarma and Singh (1997). However, for the case when source lies on the surface of the model h ¼ 0; the first integral on the right-hand side is computed using the digital linear filter algorithm, whereas the second integral, whose integrand becomes divergent, is computed using the unique functional relationship given by Z N lJ1 ðlaÞJ0 ðlrÞ dl 0 8 a 2 1 1 a2 > for ar2 o1; > 2rðr2 a2 Þ F 2; 2; 2; r2 > > > < 1 2 2 ¼ ð6Þ for ar 2 o1; F 12; 12; 1; ar 2 2 2 > ða r Þ > > > > 2 : N for ar2 ¼ 1; where F ð:Þ is the hypergeometric function or Gauss hypergeometric series. This relationship enables us to compute the Hz field response at any arbitrary point inside or outside the source loop. In a similar manner, the expressions of other field components can be modified to a form suitable for numerical computation for any position of the source loop and for any arbitrary position of the measurement point either inside or outside the source loop. More explicitly, the expressions of Hr (at z ¼ 0) and Ef (at z ¼ 0) fields can be arranged and computed as follows: Z Ia N u0 h Hr ¼ ½e ð1 rTE Þ elh 2 0 lJ1 ðlaÞJ1 ðlrÞ dl Z Ia N lh þ e lJ1 ðlaÞJ1 ðlrÞ dl; ð7Þ 2 0 Z iom0 Ia N u0 h l e ð1 þ rTE Þ elh Ef ¼ 2 u0 0 J1 ðlaÞJ1 ðlrÞ dl Z iom0 Ia N lh e J1 ðlaÞJ1 ðlrÞ dl: 2 0
ð8Þ
ARTICLE IN PRESS N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307
In Eqs. (7) and (8), when the source lies in air at a height z ¼ h; both the integrals on the right-hand side are convergent and can be computed using the previous digital linear filter algorithm for the Hankel transform. For the case, when the source lies on the surface of the model h ¼ 0; the first integral in both the equations can be computed using the digital filter algorithm for Hankel transform, whereas the second integrals which becomes divergent can be computed using the following relationship given by Eqs. (9) and (10), respectively. Z N lJ1 ðlaÞJ1 ðlrÞ dl 0 8 2i a2 a2 > 3 5 a2 > for F ; 0; ; r2 o1; > 2 2 2 r > p r2 ðr2 a2 Þ > > < r2 2 2 ¼ 2i ð9Þ for ar 2 o1; F 32; 0; 52; ar 2 > 2 2 2 > a ða r Þ p > > > > 2 : N for r ¼ 1 a2
and Z N 0
J1 ðlaÞJ1 ðlrÞ dl 8 a 2 > F 12; 32; 2; ar2 > 2 > 2r > < r 2 F 12; 32; 2; ar 2 ¼ 2 > 2a > > > : N
2 for ar2 o1; 2 for ar 2 o1; 2 for ar2 o1:
ð10Þ
4. Computer program A FORTRAN 77 program EMLCLLER to compute the forward EM response of a large circular loop source over a layered earth model at any arbitrary point inside or out side the source loop is developed and presented in this paper. The program uses two input files. The input file ‘INPUTR.IN’ comprises of survey parameters, viz. loop radius, source-receiver offset, source loop height, frequency limits and current circulating in the loop, and the file ‘MODELR.IN’ consists of model parameters, such as resistivity and thickness of the layers. Thus, the required survey parameters are supplied through INPUTR.IN, whereas, the desired model parameters are adjusted via MODELR.IN. The output file in this program RLDR.DAT, AMPR.DAT and CHEKR.DAT contain the desired output results. In their present form, these files are written to print out the amplitude and phase, normalized amplitude and phase, and real and imaginary part of the vertical magnetic field, respectively. The main program consists of a number of function subroutines and subprograms. It is developed and compiled in Compaq Visual Fortran (Professional Edition V6.5, Intel Version) on a PC (Pentium III, Toshiba Satellite model). Moreover, the
program takes just few seconds more for computation of EM response over the homogeneous earth ð0:01 s=mÞ with displacement current factor in comparison to those without displacement current factor.
5. Examples and discussions 5.1. Check for validity of the program For demonstrating the validity of the program, we have applied it for computing the EM response of a large circular loop source over the layered earth models, and the results computed using this program are compared with those of published results in literature. 5.1.1. EM response over a homogeneous earth model As the first check for validity and accuracy of the program, we have applied it for computing the EM response at the center of a circular loop source over a homogeneous earth model (as shown in inset of Fig. 2) and the results showing the variation of real and imaginary component of Hz field with frequency are presented in Fig. 2. From Fig. 2, it is noticed that our results are in agreement with those given by Ward and Hohmann (1988). Moreover, for low frequencies (corresponding to quasi-static region) both results are in coincidence with each other. This indicates the validity of the program. 5.1.2. EM response over a 2-layer earth model As another check for the validity and accuracy of the program, we have applied it for computing the EM response over a published 2-layer earth model (Ryu et al., 1970), and the results showing the variation of amplitude and phase with induction Plot of Hz with frequency at center of a circular loop source over a 0.01 S/m homogeneous earth model 1.E-01 1.E-02
Real (+)
1.E-03 Hz (A/m)
1304
1.E-04 1.E-05
Imaginary (-) Real (-)
1.E-06 1.E-07
Present study Ward and Hohmann (1988)
1.E-08 1.E-01 1.E+00 1.E+01 1.E+02 1.E+03 1.E+04 1.E+05 1.E+06 Frequency (Hz) a = 50 m a r = 0.0 m σ = 0.01 S/m
Fig. 2. Comparison with results of Ward and Hohmann (1988) for Hz field at center of circular loop of radius 50 m (carrying a unit current) over 0.01 s/m homogeneous earth model.
ARTICLE IN PRESS N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307 Amplitude vs induction number Normalized amplitude (IHz/Hz0I)
1.E-04 Amplitude (A/m)
h= 100 m 12.5 m
1.E-05 r
a
1. 56 m
a = 10 m r = 100 m
1.E-06 h σ = 0.01 S/m 1
1.E-07 1.E-02 (A)
σ
1.E-01
1.E+00
1.E+01
Normalized amplitude vs frequency for a 2-layer model 1.5 h/r = 0.025 0.05
1 1.0
0.5
0.2
σ1= 0.001 S/m
0 1.E+00
1.E+02
(A)
1.E+04
Frequency (Hz)
a = 200 m r = 2a
Phase vs frequency
h = 100 m
h/r = 0.025 0.05 0.1
0
12.5 m 1.56 m
Present study Ryu etal (1970)
1.E-01 1.E+00 Induction number, B
1.E+06 r
a
Induction number, B
Phase (in degree)
Phase (in degree)
0.1
σ1= 0.01 S/m
h
Phase vs induction number 20 0 -20 -40 -60 -80 -100 -120 -140 -160 1.E-02 (B)
1305
-60 1.0
-180
0.2
-240 -300 -360 1.E+00
1.E+01
Fig. 3. Comparison with results of Ryu et al. (1970) for circular loop (radius 10 m) over 2-layer model at distance of 100 m from the loop center: (A) amplitude of Hz field versus induction number; and (B) phase versus induction number.
-120
1.E+02 1.E+04 Frequency (Hz)
(B)
1.E+06
Fig. 4. EM response of loop source over 2-layer model at offset point outside loop at r ¼ 2a: (A) normalized amplitude versus frequency curves; and (B) phase versus frequency curves.
5.2. Examples illustrating applicability of the program For illustrating the applicability and utility of the program, we have applied it for computing the EM response over different 2- and 3-layer models with varying conductivity contrast and thickness of the layers for arbitrary in loop and offset loop measurement points. Computations are performed for the normalized amplitude (absolute amplitude of the field component over the model divided by its free space value) and phase of Hz field (per unit current) over the various 2- and 3layer models and the results are presented as a function of frequency in Figs. 4–6. Fig. 4A shows computed results for the normalized amplitude of Hz field with frequency, and Fig. 4B shows corresponding results for the phase of Hz field with frequency over a 2-layer earth model (as shown in Fig. 4) for the receiver position outside the source loop. Fig. 5A and B show normalized amplitude versus frequency and phase versus frequency curves, respectively, for a 2-layer model (as shown in
1.5
h/r = 4.0
1.0 0.4 0.2
0.1 0.05
1 0.5 0 1.E+00
h
σ1=0.01 S/m
(A)
σ1=0.01 S/m
1.E+02 1.E+04 Frequency (Hz) a r
1.E+06 a = 200 m r = a/2
Phase vs frequency 90 Phase (in degree)
number, B ¼ ðom0 s1 =2Þ1=2 r are presented in Fig. 3. From the Fig. 3, it is observed that our results are in coincidence with Ryu et al. (1970) for the quasi-static region ðf o50 kHzÞ and are extending to the non-quasistatic region. This justifies the accuracy and validity of the program.
Normalized amplitude (IHZ/HZ0I)
Normalized amplitude vs frequency for a 2-layer model 2
60 30
h/r = 4.0 1.0 0.4
0 -30 -60 -90 1.E+00
(B)
0.2
0.1
0.05
1.E+02 1.E+04 Frequency (Hz)
1.E+06
Fig. 5. EM response of loop source over 2-layer model at an inloop point at r ¼ a=2: (A) normalized amplitude versus frequency curves; and (B) phase versus frequency curves.
Fig. 5) for the receiver position inside the source loop. From Figs. 4 and 5, it is observed that all the curves show their characteristic variations. The results
ARTICLE IN PRESS N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307
1306
Normalized amplitude vs frequency for a 3-layer model Normalized amplitude (IHz/Hz0I)
1.5
Phase (in degree)
(A)
h2/h1= 10
1.0 0.5 0.1
1
0.5
0 1.E+00
1.E+02
1.E+04
1.E+06
there may be some decrease in accuracy of Bessel functions. Moreover, the digital filter algorithm for computation of Hankel transform uses filter coefficients of Guptasarma and Singh (1997), which have maximum accuracy up to 9-digits for the values of r in the range (0.01–100). Therefore, the program, EMLCLLER has maximum accuracy for the values of r in the range 0.01– 100, the value of xo8; and for the measurement points neither very close to the loop nor too far (i.e. at infinite distance) away from the source loop ðraa and NÞ:
Frequency (Hz) Phase vs frequency for a 3-layer model 20 h2/h1= 0.1 0 0.5 -20 -40 10 -60 -80 h = 10 m σ =0.01 S/m 1 1 r -100 h a σ2=0.3 S/m 2 -120 a = 50 m r = 2a σ3=0.001 S/m -140 -160 1.E+00 1.E+02 1.E+04 1.E+06
(B)
Frequency (Hz)
Fig. 6. EM response of loop source over 3-layer model at offset point at r ¼ 2a: (A) normalized amplitude versus frequency curves; and (B) phase versus frequency curves.
demonstrate the capability of the program to compute the EM response at any arbitrary measurement point either inside or outside the source loop. Fig. 6A presents normalized amplitude of Hz field versus frequency curves for a 3-layer earth model (as shown in Fig. 6), and Fig. 6B presents corresponding phase versus frequency curves, respectively, for the receiver position outside the source loop. From Fig. 6, it is noticed that program is capable of computing the EM response over the 3-layer model and, thus, can be used for a multi-layer model. 5.3. Criteria for optimum accuracy of the program Like any other program, the performance and accuracy of EMLCLLER is also dependent on the accuracy of its several intermediate subroutines and subprograms. For instance, the program uses a function subprogram hypgeoða; b; c; zÞ for the computation of Gauss hypergeometric series/function F ða; b; c; zÞ; where a; b; c and z are the coefficients of the series. The subprogram hypgeo converges for jzjo1 and gives better results for the moderate values of a; b and c (Press et al., 1992). Further, the hypgeo and main program itself involves Bessel subprograms bessj0ðxÞ and bessj1ðxÞ (Press et al., 1992), which make use of an approximate expression for large values of arguments ðx ¼ la > 8Þ: Therefore, for large values of Bessel arguments ðx > 8Þ;
6. Conclusions The present paper illustrates a new program for computing the forward EM response of a large circular loop source on or above the surface of a layered earth in both quasi-static as well as non-quasi-static (50 kHzpfp1000 kHz) regions. It makes use of both conduction as well as displacement currents and is capable of computing the EM response at any arbitrary point inside or outside the large source loop. Moreover, it is equally suitable for computing the EM response for any position of source loop either in air or on the surface of the model. Numerical examples illustrating the coincidence of computed EM response of a large loop source over the homogeneous earth and 2-layer models with those of the published response, justify the validity and accuracy of the program. The overall results demonstrate the capability of the program to compute the EM response for the higher frequency range, for an arbitrary measurement point inside or outside the source loop, and for any position of source loop in air or on the surface of the model, where as other available algorithms fail due to various reasons. The program can be applied for interpretation and inversion of large loop EM data. It can be also used as an aide to the EM depth sounding for determining the normal correction to the Turam observations, where a circular loop of finite size is generally used as transmitter. Moreover, the program can be used for approximate computation of EM response of a large square loop source, if the receiver lies at the center of the loop or far away from the source loop.
Acknowledgements The authors are grateful to the anonymous reviewers for their constructive suggestions, which certainly improved the manuscript. One of the author (NPS) gratefully acknowledges the financial assistance provided by Japan Society for Promotion of Sciences (JSPS) for the completion of this work in the form of JSPS Post-doctoral fellowship for foreign researchers.
ARTICLE IN PRESS N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307
References Anderson, W.L., 1979. Numerical integration of related Hankel transform of orders 0 and 1 by adaptive digital filtering. Geophysics 44, 1287–1305. Fraser, D.C., Stodt, J.A., Ward, S.H., 1992. The effect of displacement current on the response of a high frequency electromagnetic system. In: Ward, S.H. (Ed.), Geo-technical and Environmental Geophysics, Vol. II, Environmental and Groundwater. Society of Exploration Geophysicists, Tulsa, Oklahoma, pp. 89–95. Frischknecht, F.C., 1967. Fields about an oscillating magnetic dipole over a 2-layer earth. Colorado School of Mines Quarterly 62 (1), 370. Guptasarma, D., Singh, B., 1997. New digital linear filters for Hankel J0 and J1 transforms. Geophysical Prospecting 45, 745–762. Kaufman, A.A., Eaton, P.A., 2001. The Theory of Inductive Prospecting. Elsevier, Amsterdam, 681pp. Kaufman, A.A., Keller, G.V., 1983. Frequency and Transient Soundings. Elsevier, Amsterdam, 685pp. Knight, J.H., Raiche, A.P., 1982. Transient electromagnetic calculations using the Gaver–Stehfest inverse Laplace transform method. Geophysics 47, 47–50. Morrison, H.F., Phillips, R.J., O’Brien, D.P., 1969. Quantitative interpretation of transient electromagnetic
1307
fields over a layered half space. Geophysical Prospecting 21, 1–20. Patra, H.P., Mallick, K., 1980. Geo-sounding Principles, Vol. 2. Elsevier, Amsterdam, 419pp. Poddar, M., 1982. A rectangular loop source of current on a two-layered earth. Geophysical Prospecting 30, 101–114. Poddar, M., 1983. A rectangular loop source of current on a multi-layered earth. Geophysics 48, 107–109. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in Fortran, 2nd Edition. Cambridge University Press, Cambridge, 934pp. Ryu, J., Morrison, H.F., Ward, S.H., 1970. Electromagnetic field about a loop source of current. Geophysics 35, 862–896. Vanyan, L.L., 1967. Electromagnetic Depth Sounding. Science Consultant Bureau, New York, 312pp (selected and translated by G.V. Keller). Wait, J.R., 1982. Geo-electromagnetism. Academic Press, New York, 268pp. Ward, S.H., 1967. The electromagnetic methods. In: Ward, S.H. (Ed.), Mining Geophysics, Vol. 2. Society of Exploration Geophysicists, Tulsa, OK, pp. 224–372. Ward, S.H., Hohmann, G.W., 1988. Electromagnetic theory for geophysical applications. In: Nabighian, M.N. (Ed.), Electromagnetic Methods in Applied Geophysics, Theory, Vol. 1. Society of Exploration Geophysicists, Tulsa, OK, pp. 131–311.
Computers & Geosciences 29 (2003) 1301–1307
EMLCLLER—a program for computing the EM response of a large loop source over a layered earth model$ N.P. Singh*, Toru Mogi Institute of Seismology and Volcanology, Hokkaido University, N10W8 Kita-ku, Sapporo 060-0810, Japan Received 15 August 2002; received in revised form 24 July 2003; accepted 5 August 2003
Abstract A simple program EMLCLLER coded in FORTRAN 77 is presented in this paper for modeling the electromagnetic (EM) response of a large circular loop source over a layered earth model in quasi-static as well as non-quasi-static (general frequency) regions. The program is based on a semi-analytical numerical approach for evaluating the improper integrals occurring in the expressions of EM field components. It takes into account both conduction as well as displacement current factors. The program is formulated in such a way that it is well suited for computing the EM response for any arbitrary position of the source loop in air or on the surface of the model, in contrast to the earlier methods which face convergence problem for source position on the surface of the model. The program has wide application and is capable of computing the EM response at any arbitrary point either inside or outside the source loop. The validity and accuracy of the program is demonstrated by computing the EM response of a large loop source over the homogeneous and multi-layer earth models. Response curves depict their characteristic variations. The computed results are in coincidence with the published results for the quasi-static region (fp50 kHz) and are extension of their characteristic variation in the non-quasi-static (50 kHzpfp1000 kHz) region. Matching of computed results with the published results demonstrate the validity of the program. r 2003 Elsevier Ltd. All rights reserved. Keywords: EM response; Large loop source; Layer earth model; Conduction current; Displacement current
1. Introduction Computation of theoretical geophysical response of a model is a routine step in the problems dealing with modeling and inversion schemes of geophysical data interpretation. The problem of computation of EM response over a layered earth structure due to various EM sources has attracted the attention of scientists since the early inception of EM prospecting methods. There $ Code available from http://www.iamg.org/CGEditor. index.htm. *Corresponding author. Department of Geophysics, Banaras Hindu University, Varanasi 221005, India. Tel.: +91-5422315878; fax: +91-542-2368174. E-mail addresses: [email protected], singhnp [email protected] (N.P. Singh), [email protected] (T. Mogi).
have been many studies related to the small loop of infinitesimal size (i.e. magnetic dipole) (e.g. Frischknecht, 1967; Vanyan, 1967; Ward, 1967; Wait, 1982; Kaufman and Keller, 1983; Ward and Hohmann, 1988; Kaufman and Eaton, 2001), whereas those related to the large loop of finite size are limited (Morrison et al., 1969; Ryu et al., 1970; Patra and Mallick, 1980; Poddar, 1982,1983; Ward and Hohmann, 1988; Kaufman and Eaton, 2001). Ryu et al. (1970) presented the EM response of a large circular loop source over a layered earth model on the basis of Morrison et al. (1969) formulation. Poddar (1982, 1983) has presented the EM response of a rectangular loop source over the layered earth models. In their studies, Ryu et al. (1970) were restricted in presenting the results for the quasi-static approximation pertaining to low frequencies, and only for the measurement points outside the source loop, whereas,
0098-3004/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2003.08.008
ARTICLE IN PRESS 1302
N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307
Poddar (1982, 1983) results lack for the higher frequency region, neglecting the displacement current and manually setting a non-zero (close to zero) value for the zero loop height, i.e. when source loop lies on the surface of the model. Neglect of the displacement current in the ground and/or in air may not be valid for certain geological environments and/or frequencies higher than 50 kHz: For example, with a multi-coil DIGHEMIV system, which employs three frequencies ð900; 7200 Hz and 56 kHzÞ; the neglect of displacement currents at the highest frequency ð56 kHzÞ becomes a problem when dealing with highly resistive regions. Further, in highly resistive environments, such as in certain regions of northern Canada, which are characterized by Precambrian igneous rocks overlain by electrically thin overburden and permafrost, the displacement current factor cannot be neglected (Fraser et al., 1992). Thus, it is pertinent to develop a suitable computational technique, which can compute, without limitations, the EM response comprising of both conduction as well as displacement currents over the wide frequency range and at any arbitrary receiver position either inside or outside the source loop. Therefore, in view of practical applications and filling the gap in existing literature, we developed a 1-D forward modeling program, EMLCLLER, which is complete and capable of computing the EM response of a large loop source over the layered earth model. Moreover, the development of such an algorithm is important and useful for the interpretation and inversion of large loop EM data, which may be preferable because of low transmitter power requirements, and convenience, economy and time saving in measuring over the entire ground inside as well as outside the source loop.
2. Theoretical considerations The geometry of the earth model, the loop source and co-ordinate system under consideration are shown in Fig. 1. The earth model comprises a sequence of nhomogeneous layers, and the co-ordinate system (x; y; z) is taken with its z-axis directed vertically downward. For taking the advantage of circular symmetry, cylindrical co-ordinate system ðr; f; zÞ has been also used frequently. The integral expressions of EM field components due to a finite horizontal circular loop of radius a; carrying a current Ieiot and placed at height z ¼ h above the surface of an n-layer earth can be given by integrating the transverse electric (TE) expression of potential at a point due to a vertical magnetic dipole, over the entire loop (Ward and Hohmann, 1988). The expression of various field components, namely Hz ; Hr and Ef at a measurement point on the surface
I(ω) ω
z= -h
aa
σ0, µ0, ε0 z= 0
( 0,0,0 )
z= z1 Y z= zi
h1
σ1, µ1, ε1
hi
σi,
µi,
X
εi
z= zi+1 z= zn-1
σn, µn, εn Z
Fig. 1. Earth model and measurement system under consideration.
of earth ðz ¼ 0Þ can be written as Hz ¼
Hr ¼
Ia 2 Ia 2
Z
N
½eu0 h ð1 þ rTE Þ 0
Z
l2 J1 ðlaÞJ0 ðlrÞ dl; u0
ð1Þ
N
½eu0 h ð1 rTE ÞlJ1 ðlaÞJ1 ðlrÞ dl
ð2Þ
0
and iom0 Ia Ef ¼ 2
Z
N
½eu0 z ð1 þ rTE Þ 0
l J1 ðlaÞJ1 ðlrÞ dl; ð3Þ u0
where rTE ¼ ðY0 Y# 1 Þ=ðY0 þ Y# 1 Þ with Y0 ¼ u0 =iom0 (intrinsic admittance of free space) and Y# 1 ¼ HyTE =ExTE ¼ ðHxTE =EyTE Þ (surface admittance at z ¼ 0). For an n-layer case, the surface admittance are given by the recurrence relation: Y# 2 þ Y1 tanhðu1 h1 Þ ; Y# 1 ¼ Y1 Y1 þ Y# 2 tanhðu1 h1 Þ Y# nþ1 þ Yn tanhðun hn Þ Y# n ¼ Yn Yn þ Y# nþ1 tanhðun hn Þ and Y# n ¼ Yn with Yn ¼ un =iomn ; un ¼ ðkx2 þ ky2 kn2 Þ1=2 ¼ ðl2 kn2 Þ1=2 ; and kn2 ¼ o2 mn en iomn sn : Here, r is the distance of the measurement point from the center of the loop. The harmonic time factor eiot is implied, and SI units are used everywhere. In numerical calculations, tanhðun hn Þ is used in its exponential form (Knight and Raiche, 1982) for stability reasons. The magnetic permeability for the rock formations in general is taken equal to the free space value, i.e. mi ði ¼ 1; 2; y; nÞ ¼ m0 :
ARTICLE IN PRESS N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307
3. Computational scheme The computation of improper integrals occurring in expressions of various EM field components, Eqs. (1)– (3), is complex because the product of two Bessel functions makes the integrand have an oscillatory nature, and this problem becomes more severe when h ¼ 0; i.e. when source loop lies on the surface of the earth model. Ryu et al. (1970), however, tried to overcome this problem by subtracting a known term inside the integral sign and adjusting the same or its equivalent analytical term outside the integral sign. But overall they only managed to adjust with this problem and accomplished their computation with certain constraints and limitations, such as (a) it is valid only for z ¼ h > 0; i.e. when source lies in air, and for the case when source lies on the surface of the model, they treated it as a limiting case, (b) it is valid for r > a; i.e. only for the points outside the loop. Moreover, in presentation of results, they were restricted only to the quasi-static approximations pertaining to the low frequencies (fp50 kHz) and for the receiver points outside the source loop. Therefore, in order to overcome these limitations, we have developed a new program using some special functional relationships and the Hankel transform method as described next. For computing the improper integrals involving the product of Bessel functions of orders 0 and 1 and occurring in the expressions of EM field components, these integrals are converted into the corresponding Hankel transforms of orders 0 and 1. Accordingly, the Hankel transform of the kernel KðlÞ of integer order n is defined as (Anderson, 1979; Guptasarma and Singh, 1997):
f ðrÞ ¼
Z
N
KðlÞJn ðrlÞ dl; for r > 0;
ð4Þ
0
where Jn is the Bessel function of first kind and of order n: The transform argument r > 0 is real but KðlÞ and therefore f ðrÞ may be complex functions of real variables. The kernel function used in the Hankel transform must be a continuous bounded function as the argument increases, ði:e: KðlÞ-0; as l-NÞ: Therefore, in order to use the Hankel transform method for evaluation of these integrals, convergence of integrals are tested. In case of slow convergence or divergent nature, they are made to converge rapidly by subtracting or adding a convergent integral expression under the integral sign and subsequently adjusting it or its equivalent analytic expression outside the integral sign. For example, the expression of Hz (given by Eq. (1)) is
1303
adjusted as follows: Z Ia N u0 h l2 Hz ¼ e ð1 þ rTE Þ l elh 2 0 u0 J1 ðlaÞJ0 ðlrÞ dl Z Ia N lh þ e lJ1 ðlaÞJ0 ðlrÞ dl: 2 0
ð5Þ
In Eq. (5), for h > 0 both integrals on the right-hand side are convergent, whereas for h ¼ 0; the first integral on the right-hand side remains convergent, whereas the second becomes divergent. Therefore, for the case when source lies in the air at a height, z ¼ h; the computation of Hz component given by Eq. (5), is carried out using the digital linear filter algorithm for the Hankel transform developed using the filter coefficients of Guptasarma and Singh (1997). However, for the case when source lies on the surface of the model h ¼ 0; the first integral on the right-hand side is computed using the digital linear filter algorithm, whereas the second integral, whose integrand becomes divergent, is computed using the unique functional relationship given by Z N lJ1 ðlaÞJ0 ðlrÞ dl 0 8 a 2 1 1 a2 > for ar2 o1; > 2rðr2 a2 Þ F 2; 2; 2; r2 > > > < 1 2 2 ¼ ð6Þ for ar 2 o1; F 12; 12; 1; ar 2 2 2 > ða r Þ > > > > 2 : N for ar2 ¼ 1; where F ð:Þ is the hypergeometric function or Gauss hypergeometric series. This relationship enables us to compute the Hz field response at any arbitrary point inside or outside the source loop. In a similar manner, the expressions of other field components can be modified to a form suitable for numerical computation for any position of the source loop and for any arbitrary position of the measurement point either inside or outside the source loop. More explicitly, the expressions of Hr (at z ¼ 0) and Ef (at z ¼ 0) fields can be arranged and computed as follows: Z Ia N u0 h Hr ¼ ½e ð1 rTE Þ elh 2 0 lJ1 ðlaÞJ1 ðlrÞ dl Z Ia N lh þ e lJ1 ðlaÞJ1 ðlrÞ dl; ð7Þ 2 0 Z iom0 Ia N u0 h l e ð1 þ rTE Þ elh Ef ¼ 2 u0 0 J1 ðlaÞJ1 ðlrÞ dl Z iom0 Ia N lh e J1 ðlaÞJ1 ðlrÞ dl: 2 0
ð8Þ
ARTICLE IN PRESS N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307
In Eqs. (7) and (8), when the source lies in air at a height z ¼ h; both the integrals on the right-hand side are convergent and can be computed using the previous digital linear filter algorithm for the Hankel transform. For the case, when the source lies on the surface of the model h ¼ 0; the first integral in both the equations can be computed using the digital filter algorithm for Hankel transform, whereas the second integrals which becomes divergent can be computed using the following relationship given by Eqs. (9) and (10), respectively. Z N lJ1 ðlaÞJ1 ðlrÞ dl 0 8 2i a2 a2 > 3 5 a2 > for F ; 0; ; r2 o1; > 2 2 2 r > p r2 ðr2 a2 Þ > > < r2 2 2 ¼ 2i ð9Þ for ar 2 o1; F 32; 0; 52; ar 2 > 2 2 2 > a ða r Þ p > > > > 2 : N for r ¼ 1 a2
and Z N 0
J1 ðlaÞJ1 ðlrÞ dl 8 a 2 > F 12; 32; 2; ar2 > 2 > 2r > < r 2 F 12; 32; 2; ar 2 ¼ 2 > 2a > > > : N
2 for ar2 o1; 2 for ar 2 o1; 2 for ar2 o1:
ð10Þ
4. Computer program A FORTRAN 77 program EMLCLLER to compute the forward EM response of a large circular loop source over a layered earth model at any arbitrary point inside or out side the source loop is developed and presented in this paper. The program uses two input files. The input file ‘INPUTR.IN’ comprises of survey parameters, viz. loop radius, source-receiver offset, source loop height, frequency limits and current circulating in the loop, and the file ‘MODELR.IN’ consists of model parameters, such as resistivity and thickness of the layers. Thus, the required survey parameters are supplied through INPUTR.IN, whereas, the desired model parameters are adjusted via MODELR.IN. The output file in this program RLDR.DAT, AMPR.DAT and CHEKR.DAT contain the desired output results. In their present form, these files are written to print out the amplitude and phase, normalized amplitude and phase, and real and imaginary part of the vertical magnetic field, respectively. The main program consists of a number of function subroutines and subprograms. It is developed and compiled in Compaq Visual Fortran (Professional Edition V6.5, Intel Version) on a PC (Pentium III, Toshiba Satellite model). Moreover, the
program takes just few seconds more for computation of EM response over the homogeneous earth ð0:01 s=mÞ with displacement current factor in comparison to those without displacement current factor.
5. Examples and discussions 5.1. Check for validity of the program For demonstrating the validity of the program, we have applied it for computing the EM response of a large circular loop source over the layered earth models, and the results computed using this program are compared with those of published results in literature. 5.1.1. EM response over a homogeneous earth model As the first check for validity and accuracy of the program, we have applied it for computing the EM response at the center of a circular loop source over a homogeneous earth model (as shown in inset of Fig. 2) and the results showing the variation of real and imaginary component of Hz field with frequency are presented in Fig. 2. From Fig. 2, it is noticed that our results are in agreement with those given by Ward and Hohmann (1988). Moreover, for low frequencies (corresponding to quasi-static region) both results are in coincidence with each other. This indicates the validity of the program. 5.1.2. EM response over a 2-layer earth model As another check for the validity and accuracy of the program, we have applied it for computing the EM response over a published 2-layer earth model (Ryu et al., 1970), and the results showing the variation of amplitude and phase with induction Plot of Hz with frequency at center of a circular loop source over a 0.01 S/m homogeneous earth model 1.E-01 1.E-02
Real (+)
1.E-03 Hz (A/m)
1304
1.E-04 1.E-05
Imaginary (-) Real (-)
1.E-06 1.E-07
Present study Ward and Hohmann (1988)
1.E-08 1.E-01 1.E+00 1.E+01 1.E+02 1.E+03 1.E+04 1.E+05 1.E+06 Frequency (Hz) a = 50 m a r = 0.0 m σ = 0.01 S/m
Fig. 2. Comparison with results of Ward and Hohmann (1988) for Hz field at center of circular loop of radius 50 m (carrying a unit current) over 0.01 s/m homogeneous earth model.
ARTICLE IN PRESS N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307 Amplitude vs induction number Normalized amplitude (IHz/Hz0I)
1.E-04 Amplitude (A/m)
h= 100 m 12.5 m
1.E-05 r
a
1. 56 m
a = 10 m r = 100 m
1.E-06 h σ = 0.01 S/m 1
1.E-07 1.E-02 (A)
σ
1.E-01
1.E+00
1.E+01
Normalized amplitude vs frequency for a 2-layer model 1.5 h/r = 0.025 0.05
1 1.0
0.5
0.2
σ1= 0.001 S/m
0 1.E+00
1.E+02
(A)
1.E+04
Frequency (Hz)
a = 200 m r = 2a
Phase vs frequency
h = 100 m
h/r = 0.025 0.05 0.1
0
12.5 m 1.56 m
Present study Ryu etal (1970)
1.E-01 1.E+00 Induction number, B
1.E+06 r
a
Induction number, B
Phase (in degree)
Phase (in degree)
0.1
σ1= 0.01 S/m
h
Phase vs induction number 20 0 -20 -40 -60 -80 -100 -120 -140 -160 1.E-02 (B)
1305
-60 1.0
-180
0.2
-240 -300 -360 1.E+00
1.E+01
Fig. 3. Comparison with results of Ryu et al. (1970) for circular loop (radius 10 m) over 2-layer model at distance of 100 m from the loop center: (A) amplitude of Hz field versus induction number; and (B) phase versus induction number.
-120
1.E+02 1.E+04 Frequency (Hz)
(B)
1.E+06
Fig. 4. EM response of loop source over 2-layer model at offset point outside loop at r ¼ 2a: (A) normalized amplitude versus frequency curves; and (B) phase versus frequency curves.
5.2. Examples illustrating applicability of the program For illustrating the applicability and utility of the program, we have applied it for computing the EM response over different 2- and 3-layer models with varying conductivity contrast and thickness of the layers for arbitrary in loop and offset loop measurement points. Computations are performed for the normalized amplitude (absolute amplitude of the field component over the model divided by its free space value) and phase of Hz field (per unit current) over the various 2- and 3layer models and the results are presented as a function of frequency in Figs. 4–6. Fig. 4A shows computed results for the normalized amplitude of Hz field with frequency, and Fig. 4B shows corresponding results for the phase of Hz field with frequency over a 2-layer earth model (as shown in Fig. 4) for the receiver position outside the source loop. Fig. 5A and B show normalized amplitude versus frequency and phase versus frequency curves, respectively, for a 2-layer model (as shown in
1.5
h/r = 4.0
1.0 0.4 0.2
0.1 0.05
1 0.5 0 1.E+00
h
σ1=0.01 S/m
(A)
σ1=0.01 S/m
1.E+02 1.E+04 Frequency (Hz) a r
1.E+06 a = 200 m r = a/2
Phase vs frequency 90 Phase (in degree)
number, B ¼ ðom0 s1 =2Þ1=2 r are presented in Fig. 3. From the Fig. 3, it is observed that our results are in coincidence with Ryu et al. (1970) for the quasi-static region ðf o50 kHzÞ and are extending to the non-quasistatic region. This justifies the accuracy and validity of the program.
Normalized amplitude (IHZ/HZ0I)
Normalized amplitude vs frequency for a 2-layer model 2
60 30
h/r = 4.0 1.0 0.4
0 -30 -60 -90 1.E+00
(B)
0.2
0.1
0.05
1.E+02 1.E+04 Frequency (Hz)
1.E+06
Fig. 5. EM response of loop source over 2-layer model at an inloop point at r ¼ a=2: (A) normalized amplitude versus frequency curves; and (B) phase versus frequency curves.
Fig. 5) for the receiver position inside the source loop. From Figs. 4 and 5, it is observed that all the curves show their characteristic variations. The results
ARTICLE IN PRESS N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307
1306
Normalized amplitude vs frequency for a 3-layer model Normalized amplitude (IHz/Hz0I)
1.5
Phase (in degree)
(A)
h2/h1= 10
1.0 0.5 0.1
1
0.5
0 1.E+00
1.E+02
1.E+04
1.E+06
there may be some decrease in accuracy of Bessel functions. Moreover, the digital filter algorithm for computation of Hankel transform uses filter coefficients of Guptasarma and Singh (1997), which have maximum accuracy up to 9-digits for the values of r in the range (0.01–100). Therefore, the program, EMLCLLER has maximum accuracy for the values of r in the range 0.01– 100, the value of xo8; and for the measurement points neither very close to the loop nor too far (i.e. at infinite distance) away from the source loop ðraa and NÞ:
Frequency (Hz) Phase vs frequency for a 3-layer model 20 h2/h1= 0.1 0 0.5 -20 -40 10 -60 -80 h = 10 m σ =0.01 S/m 1 1 r -100 h a σ2=0.3 S/m 2 -120 a = 50 m r = 2a σ3=0.001 S/m -140 -160 1.E+00 1.E+02 1.E+04 1.E+06
(B)
Frequency (Hz)
Fig. 6. EM response of loop source over 3-layer model at offset point at r ¼ 2a: (A) normalized amplitude versus frequency curves; and (B) phase versus frequency curves.
demonstrate the capability of the program to compute the EM response at any arbitrary measurement point either inside or outside the source loop. Fig. 6A presents normalized amplitude of Hz field versus frequency curves for a 3-layer earth model (as shown in Fig. 6), and Fig. 6B presents corresponding phase versus frequency curves, respectively, for the receiver position outside the source loop. From Fig. 6, it is noticed that program is capable of computing the EM response over the 3-layer model and, thus, can be used for a multi-layer model. 5.3. Criteria for optimum accuracy of the program Like any other program, the performance and accuracy of EMLCLLER is also dependent on the accuracy of its several intermediate subroutines and subprograms. For instance, the program uses a function subprogram hypgeoða; b; c; zÞ for the computation of Gauss hypergeometric series/function F ða; b; c; zÞ; where a; b; c and z are the coefficients of the series. The subprogram hypgeo converges for jzjo1 and gives better results for the moderate values of a; b and c (Press et al., 1992). Further, the hypgeo and main program itself involves Bessel subprograms bessj0ðxÞ and bessj1ðxÞ (Press et al., 1992), which make use of an approximate expression for large values of arguments ðx ¼ la > 8Þ: Therefore, for large values of Bessel arguments ðx > 8Þ;
6. Conclusions The present paper illustrates a new program for computing the forward EM response of a large circular loop source on or above the surface of a layered earth in both quasi-static as well as non-quasi-static (50 kHzpfp1000 kHz) regions. It makes use of both conduction as well as displacement currents and is capable of computing the EM response at any arbitrary point inside or outside the large source loop. Moreover, it is equally suitable for computing the EM response for any position of source loop either in air or on the surface of the model. Numerical examples illustrating the coincidence of computed EM response of a large loop source over the homogeneous earth and 2-layer models with those of the published response, justify the validity and accuracy of the program. The overall results demonstrate the capability of the program to compute the EM response for the higher frequency range, for an arbitrary measurement point inside or outside the source loop, and for any position of source loop in air or on the surface of the model, where as other available algorithms fail due to various reasons. The program can be applied for interpretation and inversion of large loop EM data. It can be also used as an aide to the EM depth sounding for determining the normal correction to the Turam observations, where a circular loop of finite size is generally used as transmitter. Moreover, the program can be used for approximate computation of EM response of a large square loop source, if the receiver lies at the center of the loop or far away from the source loop.
Acknowledgements The authors are grateful to the anonymous reviewers for their constructive suggestions, which certainly improved the manuscript. One of the author (NPS) gratefully acknowledges the financial assistance provided by Japan Society for Promotion of Sciences (JSPS) for the completion of this work in the form of JSPS Post-doctoral fellowship for foreign researchers.
ARTICLE IN PRESS N.P. Singh, T. Mogi / Computers & Geosciences 29 (2003) 1301–1307
References Anderson, W.L., 1979. Numerical integration of related Hankel transform of orders 0 and 1 by adaptive digital filtering. Geophysics 44, 1287–1305. Fraser, D.C., Stodt, J.A., Ward, S.H., 1992. The effect of displacement current on the response of a high frequency electromagnetic system. In: Ward, S.H. (Ed.), Geo-technical and Environmental Geophysics, Vol. II, Environmental and Groundwater. Society of Exploration Geophysicists, Tulsa, Oklahoma, pp. 89–95. Frischknecht, F.C., 1967. Fields about an oscillating magnetic dipole over a 2-layer earth. Colorado School of Mines Quarterly 62 (1), 370. Guptasarma, D., Singh, B., 1997. New digital linear filters for Hankel J0 and J1 transforms. Geophysical Prospecting 45, 745–762. Kaufman, A.A., Eaton, P.A., 2001. The Theory of Inductive Prospecting. Elsevier, Amsterdam, 681pp. Kaufman, A.A., Keller, G.V., 1983. Frequency and Transient Soundings. Elsevier, Amsterdam, 685pp. Knight, J.H., Raiche, A.P., 1982. Transient electromagnetic calculations using the Gaver–Stehfest inverse Laplace transform method. Geophysics 47, 47–50. Morrison, H.F., Phillips, R.J., O’Brien, D.P., 1969. Quantitative interpretation of transient electromagnetic
1307
fields over a layered half space. Geophysical Prospecting 21, 1–20. Patra, H.P., Mallick, K., 1980. Geo-sounding Principles, Vol. 2. Elsevier, Amsterdam, 419pp. Poddar, M., 1982. A rectangular loop source of current on a two-layered earth. Geophysical Prospecting 30, 101–114. Poddar, M., 1983. A rectangular loop source of current on a multi-layered earth. Geophysics 48, 107–109. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in Fortran, 2nd Edition. Cambridge University Press, Cambridge, 934pp. Ryu, J., Morrison, H.F., Ward, S.H., 1970. Electromagnetic field about a loop source of current. Geophysics 35, 862–896. Vanyan, L.L., 1967. Electromagnetic Depth Sounding. Science Consultant Bureau, New York, 312pp (selected and translated by G.V. Keller). Wait, J.R., 1982. Geo-electromagnetism. Academic Press, New York, 268pp. Ward, S.H., 1967. The electromagnetic methods. In: Ward, S.H. (Ed.), Mining Geophysics, Vol. 2. Society of Exploration Geophysicists, Tulsa, OK, pp. 224–372. Ward, S.H., Hohmann, G.W., 1988. Electromagnetic theory for geophysical applications. In: Nabighian, M.N. (Ed.), Electromagnetic Methods in Applied Geophysics, Theory, Vol. 1. Society of Exploration Geophysicists, Tulsa, OK, pp. 131–311.