# On the Voltage Response of Homogeneous Earth Models in Central Loop Electromagnetic Sounding

Authors: Mauro Parise
Journal: International Journal of Antennas and Propagation
 (2022)

---

## Abstract

An accurate series-form explicit expression is derived for the voltage induced in the small receiving loop of a central loop electromagnetic sounding system positioned above a homogeneous earth model. The solution is obtained by converting the semi-infinite integral representation for the vertical magnetic field induced at the center of the two loops into a series of simpler integrals, amenable to analytical evaluation. The developed formulation is based on replacing the exponential term in the integrand of the field integral with its power series expansion with respect to the difference of the squares of the wavenumbers in air and in the conducting medium. As a result, the vertical magnetic field and the induced voltage are finally expressed as sums of the spherical Hankel functions depending on geometrical parameters and the wavenumbers. The derived solution exhibits advantages in terms of both time cost and accuracy, when compared to the conventional numerical schemes for evaluating the Sommerfeld-type integrals and to the previously published quasistatic solution to the considered problem.


---

## Body

## 1. Introduction

Central loop electromagnetic sounding is a widespread technique for the exploration of earth structures, which consists of interpreting a set of measurements of the voltage induced in a small receiving loop placed at the center of a larger transmitting loop [1–17]. The underlying concept is that the presence of buried objects or other inhomogeneities in the subsurface structures may be revealed by the discrepancy between the recorded experimental data and the theoretical response curves associated with different earth models [8, 11–13, 18]. This implies that correct data interpretation requires the a priori knowledge of highly accurate theoretical response curves, at least for the standard homogeneous earth model. Yet, despite this necessity, only under the quasistatic field assumption can we derive an accurate analytical solution for the time-harmonic voltage induced in the receiver of a transmitter-receiver loop system located on a homogeneous ground. Thus, when the quasistatic assumption fails, the induced voltage is typically through time-consuming numerical approaches or 3D electromagnetic software tools [5, 13]. An attempt to overcome the restrictions implied by the quasistatic assumption and, at the same time, providing an efficient method to compute the voltage induced in the receiver has been made in [2], where a hybrid analytical-numerical solution valid for the general case of multilayer ground has been proposed. The solution has been shown to be robust and efficient but it suffers from the limitation that its usage requires the a priori computation of a set of numerical coefficients through an expensive iterative optimization algorithm.The purpose of this work is to derive the exact series-form explicit solution describing the voltage response of a central-loop electromagnetic sounding system positioned on a homogeneous medium. The solution is obtained by converting the semi-infinite integral representation for the vertical magnetic field induced at the center of the two loops into a series of simpler integrals, amenable to analytical evaluation. The developed formulation is based on replacing the exponential term in the integrand of the field integral with its power series expansion with respect to the difference of the squares of the wavenumbers in air and in the conducting medium. As a result, the field and the voltage are expressed in terms of spherical Hankel functions of the wavenumbers and geometrical parameters. The obtained solution is not subject to the quasistatic assumption and is valid in a wide frequency range. Numerical tests are carried out to illustrate the advantages of the proposed approach. In particular, the comparison with the outcomes from the previously published quasistatic solution reveals the failure of the latter even in the low-frequency range. This evidence further confirms the usefulness of the derived expressions for the vertical magnetic field produced at the center of the two loops and the voltage induced in the small receiving loop.


## 2. Theory

Consider a current-carrying circular loop antenna and a co-axial and co-planar small receiving loop located at heightd⟶0 above a homogeneous earth structure as sketched in Figure 1. The dielectric permittivity and electric conductivity of the ground are denoted by ϵ1 and σ1 , respectively, while the magnetic permeability is assumed to be constant everywhere and equal to that of free space μ0.Figure 1
Two horizontal co-axial loops above a homogeneous ground.The current flowing in the larger transmitting loop is taken to be time-harmonic and equal toIejωt , with ω being the angular frequency. Because of the symmetry of the problem, currents in the material medium flow only horizontally. This implies that the generated electric field has no vertical component, and as a consequence, the EM field is transverse electric (TE) with respect to the vertical direction (termed as z− direction in the following) [5]. In this situation, it is well known that the magnetic field generated at the center of the two loops has only a vertical component, whose integral representation reads as [5].(1)Hz=Ib2∫0∞1+u0−u1u0+u1e−2u0dλ2u0J1λbdλ,where b denotes the radius of the emitter and Jmξ is the m th-order Bessel function, and(2)un=λ2−kn2,(3)kn2=ω2μnϵn−jωμnσn.And with the time-harmonic factorejωt suppressed for notational simplicity. The aim of this paper is to derive exact series representations for Hz. To this goal, it is first convenient to use the differentiation properties of the Bessel functions and the identity(4)u0−u1u0+u1=−1+2u0u0+u1,so as to decompose (1) into three terms as follows:(5)Hz=−Ib2∂∂b∫0∞1−e−2u0d+2u0u0+u1e−2u0d1u0J0λbλdλ.The first two integrals on the right-hand side of (5) are special cases of the well-known tabulated Sommerfeld integral [19–23](6)Snz=∫0∞e−unzunJ0λbλdλ=−jknh02knb2+z2,and as a consequence, it yields the following:(7)Hz=jIk0b2∂∂bh02k0b−h02k0r+2jk0∫0∞e−2u0du0+u1J0λbλdλ,with r=b2+4d2. Next, to evaluate the last integral within the square brackets of (7), the numerator and denominator of the fraction in the integrand are first multiplied by u0−u1. This leads to express (7) in the following form:(8)Hz=jIk0b2∂∂bh02k0b−h02k0r+2jk01k12−k02P0−P1,with(9)Pn=∫0∞une−2u0dJ0λbλdλ.The explicit expression ofP0 may be found straightforwardly, since it is easily recognized to be related to S02h in (6) by the following differential relation:(10)P0=−1b∂∂bb∂∂b+k02S02h=jk01b∂∂bb∂∂b+k02h02k0r=jk0k02h02k0r−2k0rh12k0r+k02b2r2h22k0r.On the other hand, to evaluateP1, it is convenient to rewrite it as follows:(11)P1=−limδ⟶0∂∂δ∫0∞e−u1δe−2u0dJ0λbλdλ,and then, change the variable of integration from λ to ξ=λ/k1 and express u0 as follows:(12)u0=u12+k12−k02=k1v12+α,where v1=u1/k1=ξ2−k^12 and α=k^12−k^02, being k^n=kn/k1. This makes it possible to rewrite P1 as follows:(13)P1=−k12limδ⟶0∂Q∂δ,where(14)Q=∫0∞e−v1δ^e−2d^v12+αJ0ξb^ξdξ,and being b^=k1b, d^=k1d, and δ^=k1δ. Next, since α<1, the exponential term of the integral may be replaced with its power series expansion with respect to α in Lommel’s sense, namely,(15)e−v1δ^+2d^v12+α=e−v12d^+δ^1−d^v1α+d^1+2v1d^4v13α2−d^3+6v1d^+4v12d^224v15α3+d^15+30v1d^+24v12d^2+8v13d^3192v17α4+⋯,where in the limit as δ⟶0, the quantity 2d^ within the square brackets may be confused with q^=2d^+δ^. This makes it possible to turn (15) into the following:(16)e−v1δ^+2d^v12+α=e−v1q^1−∑n=0∞−1n2n+1!!αq^v1n+1∑m=0nv1q^−m2m!!n+m!n−m!=e−v1q^+jv1q^∑n=0∞12n+1!!αq^jv1n+1hn1jv1q^,where hn1 is the n th-order spherical Hankel function of the first kind, which is defined as follows:(17)hn1jv1q^=−−jne−v1q^v1q^∑m=0nv1q^−m2m!!n+m!n−m!.Substitution of (16) into (14) provides the expression(18)Q=−αq^2+∂∂q^∫0∞e−q^v1v1J0ξb^ξdξ+jq^∑n=1∞12n+1!!−jαq^n+1∫0∞hn1jq^v1v1nJ0ξb^ξdξ,whose integrals on the right-hand side are known in the explicit form and tabulated in [19]. In particular, upon setting R^=b^2+q^2, we make use of the well-known result ([19] No. 13.47.2)(19)∫0∞Kmq^v1v1mJ0ξb^ξdξ=1q^R^jk^1q^m−1Km−1jk^1R^,with Km⋅ being the m th-order modified Bessel function of the second kind and m=n+1/2, in combination with the identities [24](20)Kn+1/2q^v1=π2jn+3/2Hn+1/21jq^v1=−jnπq^v12hn1jq^v1−π<argq^v1≤π2,(21)Kn−1/2jk^1R^=π2−jn+1/2Hn−1/22k^1R^=−jn+1/2πk^1R^2hn−12k^1R^−π2<argjk^1R^≤π,allows to obtain the expression(22)∫0∞hn1jq^v1v1nJ0ξb^ξdξ=1q^−R^jk^1q^n−1/2jk^1R^q^hn−12k^1R^,that is, equivalently,(23)∫0∞hn1jq^v1v1nJ0ξb^ξdξ=k^1q^−R^jk^1q^nhn−12k^1R^,which leads to express (18) as follows:(24)Q=Q0+∑n=1∞12n+1!!αn+1k^1n−1q^R^nhn−12k^1R^.with(25)Q0=−αq^2+∂∂q^e−jk^1R^R^=jk^1q^α2h02k^1R^−k^1R^h12k^1R^.Making use of (24) into (13) provides the following:(26)P1=−k12limδ⟶0∂Q0∂δ+∑n=1∞12n+1!!αn+1k^1n−1∂∂δq^R^nhn−12k^1R^,and, upon introducing the differentiation properties of the spherical Bessel functions, namely ([24], No. 10.1.23, 10.1.24),(27)ddR^R^−nhn2k^1R^=−k^1hn+12k^1R^R^n,(28)ddR^R^nhn−12k^1R^=k^1R^nhn−22k^1R^,it is straightforward to conclude that(29)∂Q0∂δ^=jk^1α2h02+k^1R^−1+αq^22h12+k^1q^h22,(30)∂∂δ^q^R^nhn−12=R^nhn−12+k^1q^2R^hn−22,where the argument k^1R^ of the spherical Hankel functions has been omitted for notational simplicity. Thus, after taking the limit as δ⟶0 , it turns out that R^=k1r, k^1R^=k1r, and(31)P1=k13P1′+P1″,where(32)P1′=1−τ22jh02k1r−1jk1r1+2k12−k02d2h12k1r−2k1dh22k1r,(33)P1″=−∑n=1∞1−τ2n+12n+1!!k1rnhn−12k1r+4k1d2rhn−22k1r,and with τ=k0/k1. Use of (10) in conjunction with (31)–(33) into (8) provides, after performing the b−derivatives of the spherical Hankel functions, the following expression for the Hz−field calculated at the center of the loops(34)Hz=jIk02b22r−rbh12k0b+h12k0r+P˜0+1k0rτP˜1,being(35)P˜0=2τ21−τ2h12k0r−4k0rh22k0r+b2r2h32k0r,(36)P˜1=k1rh12k1r−21−τ21+2k12−k02d2+2drh22k1r−2k1dh32k1r+2j∑n=1∞1−τ2n2n+1!!k1rn+1hn−22k1r+4k1d2rhn−32k1r.Finally, since the radiusa of the receiver has been assumed to be sufficiently small, this implies that k0a≪1 and a≪b ; it turns out that the voltage induced in the receiver is given by (34) multiplied by jωμ0πa2.


## 3. Numerical Results

In this section, the results arising from (34) are compared with those provided by the previously published quasistatic solution to the considered problem, namely [5],(37)Hz=−Ik12b33−3+3jk1b−k12b2e−jk1b,which is obtained by assuming that both the emitter and the receiver lie on the top surface of the conducting medium to be explored (that is by assumingd=0). For the sake of validation of the proposed formulation, the results generated by (34) are also compared with the data arising from numerical integration of (1), performed by using an adaptive Gauss–Kronrod G7-K15 quadrature rule, resulting from the combination of a 7-point Gauss rule with a 15-point Kronrod rule. In particular, after mapping the semi-infinite integration interval onto the finite interval [0 : 1] through a suitable change of variable, the latter interval has been recursively decomposed until the estimates of the relative errors associated with all the sub-integrals, each one given by the relative difference between the Gaussian and Kronrod rule [25], do not exceed a fixed tolerance level.First, the three different approaches are applied to the computation of the amplitude-frequency spectrum of theHz−field produced at the center of a loop carrying 1 A of current and placed on a clay soil. The index N at which the infinite sum in (36) is terminated is taken to be equal to 15, while two distinct values for the radius b of the source loop are considered, that is, b=10 m and b=100 m. Finally, the electrical conductivity and dielectric permittivity of the clay soil are taken to be σ1=10 mS/m and ϵ1=10ϵ0 [8], respectively. The obtained curves, illustrated in Figure 2, show how the numerical data agree well with the outcomes of (34) all over the considered frequency interval. Thus, the results from (34) may be referred to as exact curves. Conversely, discrepancy is observed between the quasistatic and exact trends even in the low-frequency range, to a greater extent the larger the radius of the source loop. Moreover, the quasistatic curves start to deviate from the exact ones at the frequencies of about 3 MHz when b=10 m, and 300 kHz for b=100 m. This means that the extension of the quasistatic frequency range decreases as b is increased.Figure 2
Amplitude-frequency spectrum of theHz−field at the center of the two loops.The decrease in accuracy of the quasistatic solution as frequency grows up is further illustrated in Figure3, which shows the relative error originating from substituting (37) in place of (34)–(36) with N=15, which is again plotted against frequency. Here, the electromagnetic parameters of the material medium are given the same values as in the previous example. From the analysis of the plotted curves, it is found that the relative error increases almost linearly in the low-frequency range, up to the limits of about 100% for b=10 m and b=100 m, and 30% for b=1 m. This implies that the quasistatic solution exhibits a higher accuracy if the radius of the source coil is reduced. Once the upper limit is reached, the error generated by (37) does not change any longer with increasing frequency.Figure 3
Relative error of the outcomes of (37) as compared to the results from (34), plotted against frequency.One would ask how quickly (34) converges to the exact solution for different values of the radius b of the source loop. This aspect is investigated by Figure 4, which depicts the relative truncation error that results from using (34) with a truncated sum in place of the infinite sum, plotted against b with the truncation index N taken as a parameter. The operating frequency is assumed to be equal to 100 kHz, and the estimate of the infinite sum is obtained by terminating the sequence of partial sums when the relative difference between the two last partial sums is smaller than 10−14. The same material medium as in the preceding examples is considered. A glance at Figure 4 reveals that convergence is in general slower for small but especially for large values of b. For instance, if b=100 m it suffices to use a sum made up of N=15 terms to get an acceptable relative truncation error (1.3⋅10−8), but as b is further increased the accuracy of the result of the computation rapidly drops, and it can be enhanced only by increasing the truncation order (up to N=30, as an example). However, the need to resort to longer sums does not imply that the proposed approach is more computationally expensive than numerical integration. In fact, with accuracy being equal, the time cost of the analytical approach is always smaller than that of numerical integration. The validity of this statement is confirmed by the entries of Table 1, which shows the time taken by (34) and the Gauss–Kronrod scheme, on a standard 1.8 GHz processor, to determine the Hz−field produced for b=100 m and b=500 m. The estimate of the relative error generated by Gauss-Kronrod rule and the relative truncation error associated with (34) are also indicated in the third column of Table 1. As can be observed, for b=100 m, the use of (34) with N=15 ensures approximately the same level of accuracy offered by numerical quadrature, but with significant timesaving with respect to the latter. On the other hand, when b=500 m, the relative error associated with N=15 is no longer small, and longer sums must be used in order to improve accuracy.Figure 4
Relative truncation error of (34) against the radius of the source loop, with the truncation index N taken as a parameter.Table 1
CPU time comparisons for the computation of theHz−field.
Loop radius (m)ApproachRelative errorCPU time (s)100G7-K15 (Tol =10−7)0.48⋅10−78.3100Proposed solution (N = 15)0.13⋅10−70.336500G7-K15 (Tol =10−12)0.67⋅10−1317.4500Proposed solution (N = 15)0.35⋅10−40.372500Proposed solution (N = 30)0.4⋅10−131.19The entries of Table1 confirm that increasing N up to 30 permits breaking down the truncation error and, at the same time, tell us that the improvement in accuracy is not achieved at the cost of significantly worsening the computation time, which is still far smaller than that implied by a numerical scheme which offers the same degree of precision.It has been shown how increasing the radiusb of the source loop reduces the frequency range of validity of the quasistatic assumption, and that at higher frequencies, the quasistatic solution fails regardless of the value of b. The only exception to this rule occurs when the loops lie on a poorly conducting ground. This point is clarified by Figure 5, which illustrates profiles of the voltage induced in the small receiving loop against the conductivity σ1, calculated by assuming ϵ1=10ϵ0 and a=5 cm. The operating frequency is taken to be 100 MHz, and the same values for b as in the previous example are considered.Figure 5
Voltage induced in the small loop as a function of the ground conductivity.As is seen, for larger source loops (b=100 m), there is no chance to obtain accurate results by using the quasistatic solution for the voltage, and this is expected because of the relatively high operating frequency. Conversely, for b=10 m, the quasistatic assumption still holds at 100 MHz, provided that the ground conductivity σ1 does not exceed few mS/m.


## 4. Conclusion

The purpose of the present work has been to derive the exact explicit expression describing the time-harmonic voltage response of a central loop electromagnetic sounding system placed on a homogeneous earth structure. The solution has been obtained by turning the semi-infinite integral representation for the vertical magnetic field induced at the center of the emitter into a series of simpler integrals amenable to analytical evaluation. The series expansion of the field integral has been made possible by the replacement of the exponential term of the integrand with its power series expansion with respect to the difference of the squares of the wavenumbers in air and in the conducting medium. As a result, the vertical magnetic field and the induced voltage have been expressed in terms of the spherical Hankel functions of the wavenumbers and geometrical parameters. The results of numerical tests have shown the validity of the proposed approach. In particular, it has been seen how the derived solution exhibits advantages in terms of time cost and accuracy, when compared, respectively, to the numerical quadrature rules devoted to the evaluation of the Sommerfeld-type integrals, and to the previously published quasistatic solution to the problem.



---

*Source: 8294000-2022-06-16.xml*