## Computation of potentials from current electrodes in cylindrically stratified media: A stable, rescaled semi-analytical formulation

H. Moon a, ∗ , F. L. Teixeira a , B. Donderici b

a ElectroScience Laboratory, The Ohio State University, Columbus, OH 43212, USA b Sensor Physics &amp; Technology, Halliburton Energy Services, Houston, TX 77032, USA

## Abstract

We present an efficient and robust semi-analytical formulation to compute the electric potential due to arbitrary-located point electrodes in three-dimensional cylindrically stratified media, where the radial thickness and the medium resistivity of each cylindrical layer can vary by many orders of magnitude. A basic roadblock for robust potential computations in such scenarios is the poor scaling of modified-Bessel functions used for computation of the semi-analytical solution, for extreme arguments and/or orders. To accommodate this, we construct a set of rescaled versions of modified-Bessel functions, which avoids underflows and overflows in finite precision arithmetic, and minimizes round-off errors. In addition, several extrapolation methods are applied and compared to expedite the numerical evaluation of the (otherwise slowly convergent) associated Sommerfeld-type integrals. The proposed algorithm is verified in a number of scenarios relevant to geophysical exploration, but the general formulation presented is also applicable to other problems governed by Poisson equation such as Newtonian gravity, heat flow, and potential flow in fluid mechanics, involving cylindrically stratified environments.

Keywords: Poisson equation, steady-state diffusion equation, discontinuous coefficients, stratified media, resistivity logging, electric potential

## 1. Introduction

Resistivity logging is extensively used for detecting, characterizing, and analyzing hydrocarbon-bearing zones in the subsurface earth [30, 17, 7, 4, 13, 14, 35, 3]. This sensing modality employs electrode-type devices mounted on a mandrel that inject electric currents into the surrounding earth formation [5, 29]. The ensuing electric potential is then measured at different locations to provide estimates for the surrounding resistivity. Many numerical techniques such as finite-differences, finite elements, numerical mode-matching, and finite volumes method can be used to model the response of resistivity logging tools [12, 8, 10, 26, 23, 24, 20, 25, 21, 22, 11, 6, 9]. Brute-force techniques are rather versatile and applicable to arbitrary resistivity distributions; however, at the same time, this precludes optimality in particular cases of special interest, such as when resistivity logging environment can be represented as a cylindrically stratified medium [2]. Depending on the implementation, brute-force techniques may have difficulties handling extreme sharp discontinuities in the coefficients, as is the case for the resistivity parameter for the physical scenario considered here, which can change by many orders of magnitude across adjacent layers.

In this paper, a robust semi-analytical formulation for computing the electric potential due to arbitrarylocated point electrodes in three-dimensional cylindrically stratified media is proposed. The present formulation is based on a series expansion in terms of azimuth Fourier modes and a spectral integral over the vertical wavenumber along the axial direction. The resulting problem in terms of the radial variable yields

∗ Corresponding author

Email addresses:

moon.173@osu.edu burkay.donderici@halliburton.com

(H. Moon),

(B. Donderici)

teixeira@ece.osu.edu

(F. L. Teixeira),

a set of modified Bessel equations. The present formulation removes roadblocks for numerical computations associated with the poor scaling of modified-Bessel functions for very small and/or very large arguments and/or orders [27, 31, 1]. This is done by constructing a set of rescaled, modified-Bessel functions that can be stably evaluated under double-precision arithmetic, akin to what has been done in the past for ordinary (non-modified) Bessel functions [16] . The present formulation also carefully manipulates the analytical formulae for the potential in such media to yield a set of integrand expressions can be computed in a robust manner under double-precision for a wide range of layer thicknesses, layer resistivities, and source and observation point separations. Finally, a number of acceleration techniques are implemented and compared to effect the efficient numerical integration of the Sommerfeld-type (spectral) integrals, which otherwise suffer from slow convergence. The proposed algorithm is verified in a number of practical scenarios relevant to geophysical exploration. More generally, the mathematical setting here corresponds to the classical problem of obtaining the Green's function for the steady diffusion equation (Poisson problem) with discontinuous coefficients in a separable geometry. As such, the general formalism presented here is also applicable to other problems governed by Poisson equation such as Newtonian gravity, heat flow, elasticity, neutron transport, and potential flows in fluid mechanics, in cylindrically stratified geometries.

## 2. Formulation

## 2.1. Electric potential in homogeneous media

In a homogeneous medium, the electric potential ψ from a current electrode at the origin writes as [31]

$$\psi = \frac { \mathcal { I } } { 4 \pi \sigma \sqrt { \rho ^ { 2 } + z ^ { 2 } } } = \frac { \mathcal { I } } { 2 \pi ^ { 2 } \sigma } \int _ { 0 } ^ { \infty } K _ { 0 } ( \lambda \rho ) \cos ( \lambda z ) d \lambda,$$

where I is the electric current flowing into the medium from the electrode, σ is the conductivity of the medium, and K 0 ( ) · is the modified-Bessel function of the second kind of the zeroth order. For the second equality, the complete Lipschitz-Hankel integral [32] is employed. When the source is off the origin, higher order azimuthal modes appear. Using the addition theorem for K 0 , (1) is modified to

$$\psi & = \frac { \mathcal { I } } { 2 \pi ^ { 2 } \sigma } \int _ { 0 } ^ { \infty } K _ { 0 } ( \lambda | \rho - \rho ^ { \prime } | ) \cos ( \lambda ( z - z ^ { \prime } ) ) d \lambda \\ & = \frac { \mathcal { I } } { 2 \pi ^ { 2 } \sigma } \sum _ { n = - \infty } ^ { \infty } e ^ { i n ( \phi - \phi ^ { \prime } ) } \int _ { 0 } ^ { \infty } I _ { n } ( \lambda \rho _ { < } ) K _ { n } ( \lambda \rho _ { > } ) \cos ( \lambda ( z - z ^ { \prime } ) ) d \lambda,$$

in terms of modified-Bessel functions of both first, I n ( ) and second, · K n ( ), kinds. · In the above, primed coordinates ( ρ , φ , z ′ ′ ′ ) represent the source location and unprimed coordinates ( ρ, φ, z ) represent the observation point. Also, ρ &lt; = min( ρ, ρ ′ ) and ρ &gt; = max( ρ, ρ ′ ).

## 2.2. Electric potential in cylindrically stratified media

In a cylindrically stratified medium, boundary conditions at the interfaces need to be incorporated. Let us first consider the case with two distinct cylindrical layers, as depicted in Fig. 1. When the source is embedded in layer 1, we denote it the outgoing-potential case. In this case, the primary potential ψ p is a function of K n ( λρ ) because I n ( λρ ) diverges as ρ goes to infinity. On the other hand, when the source is embedded in layer 2, we denote it the standing-potential case and ψ p is a function of I n ( λρ ) instead of K n ( λρ ) because K n ( λρ ) diverges when ρ goes to zero. For the outgoing-potential case, the n -th harmonic with e i n φ ( -φ ′ ) dependence in layer 1 and layer 2 can be expressed, resp., as

$$\psi _ { 1 } = \left [ K _ { n } ( \lambda \rho ) + R _ { 1 2 } I _ { n } ( \lambda \rho ) \right ] A _ { 0 },$$

ψ

2

=

T

12

K

n

(

λρ A ,

)

0

(3b)

and

ρ

Figure 1: Schematic description of two layers with associated coefficients in the ρz -plane.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-3280f57a-a04d-4940-a179-15fddfea0ef3.jpg)

where R 12 and T 12 are the (local) reflection and transmission coefficients at the boundary a 1 , and A 0 is an arbitrary amplitude of the primary potential. Applying the boundary conditions [31] at the interface, we obtain

$$R _ { 1 2 } = \frac { ( \sigma _ { 2 } - \sigma _ { 1 } ) K _ { n } ( \lambda a _ { 1 } ) K _ { n } ^ { \prime } ( \lambda a _ { 1 } ) } { \sigma _ { 1 } I _ { n } ^ { \prime } ( \lambda a _ { 1 } ) K _ { n } ( \lambda a _ { 1 } ) - \sigma _ { 2 } I _ { n } ( \lambda a _ { 1 } ) K _ { n } ^ { \prime } ( \lambda a _ { 1 } ) }, \quad \tilde { \sigma }.$$

$$T _ { 1 2 } = \frac { \mu _ { 1 } ^ { \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \colon } ^ { \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \colon } } { \lambda a _ { 1 } \left [ \sigma _ { 1 } I _ { n } ^ { \prime } ( \lambda a _ { 1 } ) K _ { n } ( \lambda a _ { 1 } ) - \sigma _ { 2 } I _ { n } ( \lambda a _ { 1 } ) K _ { n } ^ { \prime } ( \lambda a _ { 1 } ) \right ] }.$$

For the standing-potential case, we similarly have

$$\psi _ { 1 } & = T _ { 2 1 } I _ { n } ( \lambda \rho ) B _ { 0 },$$

$$\psi _ { 2 } = \left [ R _ { 2 1 } K _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) \right ] B _ { 0 },$$

$$R _ { 2 1 } = \frac { ( \sigma _ { 2 } - \sigma _ { 1 } ) I _ { n } ( \lambda a _ { 1 } ) I _ { n } ^ { \prime } ( \lambda a _ { 1 } ) } { \sigma _ { 1 } I _ { n } ^ { \prime } ( \lambda a _ { 1 } ) K _ { n } ( \lambda a _ { 1 } ) - \sigma _ { 2 } I _ { n } ( \lambda a _ { 1 } ) K _ { n } ^ { \prime } ( \lambda a _ { 1 } ) },$$

$$T _ { 2 1 } = \frac { \mu _ { 1 } ^ { \nu - n } \cdots \nu _ { 1 } ^ { \nu - \mu _ { 2 } ^ { \nu } } \cdots \nu _ { 2 } ^ { \nu - \mu _ { 3 } ^ { \nu - n } } \cdots \nu _ { 3 } ^ { \nu - \mu _ { 4 } ^ { \nu } } } { \lambda a _ { 1 } \left [ \sigma _ { 1 } I _ { n } ^ { \prime } ( \lambda a _ { 1 } ) K _ { n } ( \lambda a _ { 1 } ) - \sigma _ { 2 } I _ { n } ( \lambda a _ { 1 } ) K _ { n } ^ { \prime } ( \lambda a _ { 1 } ) \right ] }.$$

When more than two distinct layers are present, multiple reflections and transmissions occur. Therefore, generalized reflection and transmission coefficients should be determined. The procedure to obtain these coefficients is very similar to the one used for time-harmonic case in [2] and will not be derived in detail here. Fig. 2 depicts the relevant coefficients to the outgoing-potential case in the medium consisting of three cylindrical layers. The potentials in the three layers can be expressed as

$$\psi _ { 1 } & = \left [ K _ { n } ( \lambda \rho ) + \widetilde { R } _ { 1 2 } I _ { n } ( \lambda \rho ) \right ] A _ { 1 },$$

$$\psi _ { 2 } = \left [ K _ { n } ( \lambda \rho ) + R _ { 2 3 } I _ { n } ( \lambda \rho ) \right ] A _ { 2 },$$

$$\psi _ { 3 } & = T _ { 2 3 } K _ { n } ( \lambda \rho ) A _ { 3 }.$$

Applying proper constraint equations, we obtain

$$\widetilde { R } _ { 1 2 } = R _ { 1 2 } + T _ { 2 1 } R _ { 2 3 } ( 1 - R _ { 2 1 } R _ { 2 3 } ) ^ { - 1 } T _ { 1 2 }.$$

If another layer is added beyond layer 3, only R 23 needs to be replaced by ˜ R 23 . Therefore, the generalized reflection coefficient in cylindrically stratified media for the outgoing-potential case is

$$\widetilde { R } _ { i, i + 1 } & = R _ { i, i + 1 } + T _ { i + 1, i } \widetilde { R } _ { i + 1, i + 2 } ( 1 - R _ { i + 1, i } \widetilde { R } _ { i + 1, i + 2 } ) ^ { - 1 } T _ { i, i + 1 }.$$

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-bcb15dfa-576a-4aa6-87ea-839a01d9d61a.jpg)

Layer 2

R

21

23

R

Figure 2: Schematic description of three layers with associated coefficients for the outgoing-potential case: (a) R 12 and (b) S 12 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-cebe9681-fa17-4ac1-9025-be4772740fc3.jpg)



Figure 3: Generalized transmission coefficient for the outgoing-potential case.

All amplitudes A i 's as well as generalized reflection coefficients should be determined in order to obtain the potential everywhere. The relationship between two successive amplitudes in cylindrically stratified media can be written as

$$A _ { i + 1 } & = T _ { i, i + 1 } A _ { i } + R _ { i + 1, i } R _ { i + 1, i + 2 } A _ { i + 1 }. & ( 1 0 )$$

From (10), a new coefficient denoted by S is defined as

$$S _ { i, i + 1 } = ( 1 - R _ { i + 1, i } \widetilde { R } _ { i + 1, i + 2 } ) ^ { - 1 } T _ { i, i + 1 } & & ( 1 1 )$$

such that A i +1 = S i,i +1 A i . The above coefficient can be regarded as a 'local' transmission coefficient between two adjacent layers, as depicted in Fig. 2b. Generalized transmission coefficient described in Fig. 3 can be defined using the S -coefficients (11) through

$$\widetilde { T } _ { j i } = T _ { i - 1, i } S _ { i - 2, i - 1 } \cdots S _ { j, j + 1 } = T _ { i - 1, i } \prod _ { k = j } ^ { i - 2 } S _ { k, k + 1 } = T _ { i - 1, i } X _ { j, i - 1 }.$$

Note that i &gt; j for the outgoing-potential case. When i = j +1, X j,i -1 = 1 in (12). For the standing-potential case depicted in Fig. 4, the potentials in three cylindrical layers can be written a 2

La

a 1

Layer 2

T

12

R

21

R

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-d616795d-64fc-4914-89b7-a7f138373cc8.jpg)

23

Figure 4: Schematic description of three layers with associated coefficients for the standing-potential case: (a) R 32 and (b) S 32 .

as

$$\psi _ { 1 } & = I _ { n } ( \lambda \rho ) B _ { 1 }, & ( 1 3 \tt a )$$

$$\psi _ { 2 } & = \left [ R _ { 2 1 } K _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) \right ] B _ { 2 }, & ( 1 3 b )$$

$$\psi _ { 3 } & = \left [ \widetilde { R } _ { 3 2 } K _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) \right ] B _ { 3 }. & ( 1 3 c )$$

Similarly, applying proper constraint conditions, we obtain

$$\widetilde { R } _ { 3 2 } = R _ { 3 2 } + T _ { 2 3 } R _ { 2 1 } ( 1 - R _ { 2 3 } R _ { 2 1 } ) ^ { - 1 } T _ { 3 2 }.$$

and, more generally,

$$\widetilde { R } _ { i, i - 1 } & = R _ { i, i - 1 } + T _ { i - 1, i } \widetilde { R } _ { i - 1, i - 2 } ( 1 - R _ { i - 1, i } \widetilde { R } _ { i - 1, i - 2 } ) ^ { - 1 } T _ { i, i - 1 }.$$

To obtain all amplitudes B i 's, it is convenient to define the S coefficient with decreasing subscripts S i,i -1 . Since the relation between two successive amplitudes is

$$B _ { i - 1 } = T _ { i, i - 1 } B _ { i } + R _ { i - 1, i } R _ { i - 1, i - 2 } B _ { i - 1 },$$

the S coefficient in this case is defined as

$$S _ { i, i - 1 } = ( 1 - R _ { i - 1, i } \widetilde { R } _ { i - 1, i - 2 } ) ^ { - 1 } T _ { i, i - 1 }$$

such that B i -1 = S i,i -1 B i . The generalized transmission coefficient depicted in Fig. 5 can then be written as

$$\widetilde { T } _ { j i } = T _ { i + 1, i } S _ { i + 2, i + 1 } \cdots S _ { j, j - 1 } = T _ { i + 1, i } \prod _ { k = j } ^ { i + 2 } S _ { k, k + 1 } = T _ { i + 1, i } X _ { j, i + 1 }.$$

Note that i &lt; j for the outgoing-potential case. When j = i +1, X j,i +1 = 1 in (18).

Using generalized reflection and transmission coefficients, we can extend (2) to incorporate multiple reflections and transmissions. The integral part of (2) can be modified to

$$\int _ { 0 } ^ { \infty } [ I _ { n } ( \lambda \rho _ { < } ) K _ { n } ( \lambda \rho _ { > } ) + K _ { n } ( \lambda \rho ) A _ { j n } ( \rho ^ { \prime } ) + I _ { n } ( \lambda \rho ) B _ { j n } ( \rho ^ { \prime } ) ] \cos ( \lambda ( z - z ^ { \prime } ) ) d \lambda,$$

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-1e4b82c0-7cfc-49b0-b3bb-47f04ddee97c.jpg)



Figure 5: Generalized transmission coefficient for the standing-potential case.

where the two unknowns are A jn ( ρ ′ ) and B jn ( ρ ′ ). Applying two constraint conditions at ρ = a j -1 and ρ = a j , we obtain

$$A _ { j n } ( \rho ^ { \prime } ) = ( 1 - \widetilde { R } _ { j, j - 1 } \widetilde { R } _ { j, j + 1 } ) ^ { - 1 } \widetilde { R } _ { j, j - 1 } \left [ K _ { n } ( \lambda \rho ^ { \prime } ) + I _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j + 1 } \right ],$$

$$B _ { j n } ( \rho ^ { \prime } ) = ( 1 - \widetilde { R } _ { j, j + 1 } \widetilde { R } _ { j, j - 1 } ) ^ { - 1 } \widetilde { R } _ { j, j + 1 } \left [ I _ { n } ( \lambda \rho ^ { \prime } ) + K _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j - 1 } \right ].$$

Substituting (20a) and (20b) into (19) and rearranging the integrand excluding the cosine factor gives

[

I

n

(

λρ

&lt;

)

K

n

(

λρ

&gt;

) +

K

n

(

λρ A

)

jn

(

ρ

′

) +

I

n

(

λρ B

)

jn

(

ρ

′

)]

$$& \stackrel { \iota \nolimits \nu \nu ^ { \prime } \nolimits \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nolimits \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu\nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu\nu \nu \nolimits \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nolimits \nu \nu\nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nu \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nulimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimim\nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimim\nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimim\nolimits \nolimits \nolimits \nolimits \nolimits \nolimim\nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimim\nolimits \nolimits \nolimits \nolimim\nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimim\nolimits \nolimits \nolimits \nolimim\nolimits \nolimits \nolimits \nolimimits \nolimits \nolimits \nolimim\nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimits \nolimim\nolimits \nolimits \nolimits \nolimim\nolimits \nolimim\nolimits \nolimits \nolimits \nolimim\nolimits \nolimits \nolimim\nolimits \nolimits \nolimim\nolimits \nolimits \nolimim\nolimits \nolimits \nolimim\nolimits \nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimits \nolimits \nolimim\nolimits \nolimim\nolimits \nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimits \nolimim\nolimits \nolim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimim\nolimits \nolim\nolim\nolimits \nolimim\nolim\nolimits \nolim\nolim\nolim\nolim\nolimits \nolim\nolim\nolim\nolimits \nolim\nolim\nolim\nolimits \nolimim\nolimits \nolim\nolim\nolimits \nolimim\nolim\nolimits \nolim\nolimits \nolim\nolimits \nolimim\nolimits \nolim\nolim\nolim\nolim\nolimits \nolimim\nolim\nolimits \nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolimits \nolim\nolim\nolim\nolim\nolimits \nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolimits \nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolimits \nolim\nolimits \nolimim\nolimits \nolimits \nolimim\nolimits \nolimim\nolim\nolim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimits \nolimim\nolimim\nolimits \nolimim\nolimits \nolimim\nolimim\nolimits \nolimim\nolimits \nolimim\nolim\nolimits \nolimim\nolimits \nolimim\nolim\nolimim\nolimits \nolimim\nolim\nolimits \nolimim\nolimim\nolimits \nolimim\nolimits \nolim\nolim\nolim\nolimits \nolimim\nolimim\nolim\nolimits \nolimim\nolimits \nolimim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolimits \nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolimits \nolimim\nolimits \nolim\nolim\nolim\nolim\nolimits \nolimim\nolimim\nolim\nolim\nolim\nolim\nolimits \nolimim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolimits \nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim\nolim
& \quad = \begin{cases} \\ & \quad \left [ K _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) ) \end{cases} \right ] \left [ I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n, j - 1 } ( \lambda \rho ) + I _ { n, j - 1 } ) \end{cases} \right ] \left [ K _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) + I _ { n, j - 1 } ) \end{cases} \right ] \left [ K _ { n } ( \lambda \rho ) + I _ { n, j - 1 } \end{cases$$

where M j = ( 1 -˜ R j,j -1 ˜ R j,j +1 ) -1 . When i &gt; j , the potential in layer i is expressed as

$$\psi _ { i } = \frac { I } { 2 \pi ^ { 2 } \sigma _ { j } } \sum _ { n = - \infty } ^ { \infty } e ^ { i n ( \phi - \phi ^ { \prime } ) } \int _ { 0 } ^ { \infty } \left [ K _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) \widetilde { R } _ { i, i + 1 } \right ] A _ { i n } \cos ( \lambda ( z - z ^ { \prime } ) ) d \lambda,$$

where A in is the amplitude of the outgoing-potential in layer i and expressed as

$$A _ { i n } = \left ( 1 - R _ { i, i - 1 } \widetilde { R } _ { i, i + 1 } \right ) \widetilde { T } _ { j i } A _ { j n } = N _ { i + } \widetilde { T } _ { j i } A _ { j n } = N _ { i + } \widetilde { T } _ { j i } \left [ I _ { n } ( \lambda \rho ^ { \prime } ) + K _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j - 1 } \right ] M _ { j }. \quad \ ( 2 3 )$$

Therefore, the integrand factor for layer i in (22) is

$$\left [ K _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) \widetilde { R } _ { i, i + 1 } \right ] \left [ I _ { n } ( \lambda \rho ^ { \prime } ) + K _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j - 1 } \right ] N _ { i + } \widetilde { T } _ { j i } M _ { j }.$$

When i &lt; j , the potential in layer i is expressed as

$$\psi _ { i } = \frac { I } { 2 \pi ^ { 2 } \sigma _ { j } } \sum _ { n = - \infty } ^ { \infty } e ^ { i n ( \phi - \phi ^ { \prime } ) } \int _ { 0 } ^ { \infty } \left [ I _ { n } ( \lambda \rho ) + K _ { n } ( \lambda \rho ) \widetilde { R } _ { i, i - 1 } \right ] A _ { i n } \cos ( \lambda ( z - z ^ { \prime } ) ) d \lambda,$$

where A in is the amplitude of the standing-potential in layer i and expressed as

$$A _ { i n } = \left ( 1 - R _ { i, i + 1 } \widetilde { R } _ { i, i - 1 } \right ) \widetilde { T } _ { j i } A _ { j n } = N _ { i - } \widetilde { T } _ { j i } A _ { j n } \left [ K _ { n } ( \lambda \rho ^ { \prime } ) + I _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j + 1 } \right ] M _ { j }$$

The integrand of the integral for potential in layer i in (25) is

$$\left [ I _ { n } ( \lambda \rho ) + K _ { n } ( \lambda \rho ) \widetilde { R } _ { i, i - 1 } \right ] \left [ K _ { n } ( \lambda \rho ^ { \prime } ) + I _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j + 1 } \right ] N _ { i - } \widetilde { T } _ { j i } M _ { j }.$$

In summary, the potential in cylindrically stratified media admits four different expressions depending on the relative position of ρ and ρ ′ , which can be written as

$$\psi _ { i } = \frac { \mathcal { I } } { 2 \pi ^ { 2 } \sigma _ { j } } \sum _ { n = - \infty } ^ { \infty } e ^ { i n ( \phi - \phi ^ { \prime } ) } \int _ { 0 } ^ { \infty } F _ { n } ( \rho, \rho ^ { \prime } ) \cos ( \lambda ( z - z ^ { \prime } ) ) d \lambda,$$

where

Case 1 : ρ and ρ ′ are in the same layer and ρ ≥ ρ ′

$$F _ { n } ( \rho, \rho ^ { \prime } ) = \left [ K _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) \widetilde { R } _ { j, j + 1 } \right ] \left [ I _ { n } ( \lambda \rho ^ { \prime } ) + K _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j - 1 } \right ] M _ { j },$$

Case 2 : ρ and ρ ′ are in the same layer and ρ &lt; ρ ′

$$F _ { n } ( \rho, \rho ^ { \prime } ) = \left [ I _ { n } ( \lambda \rho ) + K _ { n } ( \lambda \rho ) \widetilde { R } _ { j, j - 1 } \right ] \left [ K _ { n } ( \lambda \rho ^ { \prime } ) + I _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j + 1 } \right ] M _ { j },$$

Case 3 : ρ and ρ ′ are in different layers and ρ &gt; ρ ′

$$F _ { n } ( \rho, \rho ^ { \prime } ) = \left \lceil K _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) \widetilde { R } _ { i, i + 1 } \right \rceil \left \lceil I _ { n } ( \lambda \rho ^ { \prime } ) + K _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j - 1 } \right \rceil N _ { i + } \widetilde { T } _ { j i } M _ { j },$$

Case 4 : ρ and ρ ′ are in different layers and ρ &lt; ρ ′

$$F _ { n } ( \rho, \rho ^ { \prime } ) = \left [ I _ { n } ( \lambda \rho ) + K _ { n } ( \lambda \rho ) \widetilde { R } _ { i, i - 1 } \right ] \left [ K _ { n } ( \lambda \rho ^ { \prime } ) + I _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j + 1 } \right ] N _ { i - } \widetilde { T } _ { j i } M _ { j }.$$

## 2.3. Rescaled modified-Bessel functions

The electric potential involves products of the modified-Bessel function of the first and second kind, viz. I n and K n . Those products can involve disparate values due to the exponential behavior of the functions. For example, when | z | &lt;&lt; 1, K n ( z ) has very large value whereas I n ( z ) has very small value. This disparity becomes progressively worse for higher order modes. On the other hand, when glyph[Rfractur] e z [ ] &gt;&gt; 1, K n ( z ) has very small value while I n ( z ) has very large value. This leads to unreliable results under double-precision computations. To eliminate this problem, a new set of rescaled modified-Bessel functions are defined in a similar fashion to what has been done in [16] for Bessel and Hankel functions, viz. J n and H (1) n . In order to apply such rescaled functions, the analytical expressions for the potential need to modified accordingly, as described next.

When | z | glyph[lessmuch] 1, I n ( z ) and K n ( z ) can be expressed via small argument approximations for n &gt; 0. Noting that, for -π &lt; arg z ( ) &lt; π/ 2, the relationship between cylindrical functions and modified-cylindrical functions reads as

$$I _ { n } ( z ) = \mathfrak { i } ^ { - n } J _ { n } ( i z ),$$

$$K _ { n } ( z ) & = \frac { \pi } { 2 } \mathfrak { i } ^ { n + 1 } H _ { n } ^ { ( 1 ) } ( i z ). & ( 3 0 b )$$

Thus, we have the following small argument approximations for I n and K n and their derivatives.

$$I _ { n } ( \lambda a _ { i } ) \approx \frac { 1 } { n! } \left ( \frac { \lambda a _ { i } } { 2 } \right ) ^ { n } = \frac { 1 } { n! } \left ( \frac { \lambda } { 2 } \right ) ^ { n } \cdot a _ { i } ^ { n } \cdot 1 = G a _ { i } ^ { n } \hat { I } _ { n } ( \lambda a _ { i } ),$$

$$I ^ { \prime } _ { n } ( \lambda a _ { i } ) \approx \frac { \stackrel { \cdot } { 1 } } { 2 ( n - 1 )! } \left ( \frac { \lambda a _ { i } } { 2 } \Big ) ^ { n - 1 } = \frac { 1 } { n! } \left ( \frac { \lambda } { 2 } \right ) ^ { n } \cdot a _ { i } ^ { n } \cdot \frac { n } { \lambda a _ { i } } = G a _ { i } ^ { n } \hat { I } ^ { \prime } _ { n } ( \lambda a _ { i } ),$$

$$K _ { n } ( \lambda a _ { i } ) \approx \frac { \theta ^ { \prime } - 1 )! } { 2 } \left ( \frac { \theta ^ { \prime } } { \lambda a _ { i } } \right ) ^ { n } = n! \left ( \frac { 2 } { \lambda } \right ) ^ { n } \cdot a _ { i } ^ { - n } \cdot \left ( \frac { 1 } { 2 n } \right ) = G ^ { - 1 } a _ { i } ^ { - n } \hat { K } _ { n } ( \lambda a _ { i } ),$$

$$K _ { n } ^ { \prime } ( \lambda a _ { i } ) \approx - \frac { n! } { 4 } \left ( \frac { 2 } { \lambda a _ { i } } \right ) ^ { n + i } = n! \left ( \frac { 2 } { \lambda } \right ) ^ { n } \cdot a _ { i } ^ { - n } \cdot \left ( - \frac { 1 } { 2 \lambda a _ { i } } \right ) = G ^ { - 1 } a _ { i } ^ { - n } \hat { K } _ { n } ^ { \prime } ( \lambda a _ { i } ).$$

It should be noted that the multiplicative factor G above is chosen so as to not depend on the radial distance, a i . Also, G is identical for a function and its derivative, and the multiplicative factors appearing in I n and K n are reciprocal to each other. This will facilitate some computations later on.

When glyph[Rfractur] e z [ ] glyph[greatermuch] 1 and z = λa i = ( λ ′ +i λ ′′ ) a i , the large argument approximations for the modified-Bessel functions write as [34]

$$I _ { n } ( \lambda a _ { i } ) & = \frac { e ^ { \lambda a _ { i } } } { \sqrt { 2 \pi \lambda a _ { i } } } \left [ 1 - \frac { ( \mu - 1 ) } { 1! ( 8 \lambda a _ { i } ) } + \frac { ( \mu - 1 ) ( \mu - 9 ) } { 2! ( 8 \lambda a _ { i } ) ^ { 2 } } + \cdots \right ] \\ & = e ^ { \lambda ^ { \prime } a _ { i } } \hat { I } _ { n } ( \lambda a _ { i } ), \\ & \longmapsto \frac { r } { r } \quad \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \colon$$

$$K _ { n } ( \lambda a _ { i } ) & = \sqrt { \frac { \pi } { 2 \lambda a _ { i } } } \stackrel { \dots } { e } ^ { - \lambda a _ { i } } \left [ 1 + \frac { ( \mu - 1 ) } { 1! ( 8 \lambda a _ { i } ) } + \frac { ( \mu - 1 ) ( \mu - 9 ) } { 2! ( 8 \lambda a _ { i } ) ^ { 2 } } + \cdots \right ] \\ & = e ^ { - \lambda ^ { \prime } a _ { i } } \hat { K } _ { n } ( \lambda a _ { i } ),$$

where µ = 4 n 2 . Again, the associated multiplicative factors are reciprocal to each other. The derivatives of rescaled modified-cylindrical functions for large arguments can be derived through the recursive formulas such that

$$I _ { n } ^ { \prime } ( \lambda a _ { i } ) = I _ { n - 1 } ( \lambda a _ { i } ) - \frac { n } { \lambda a _ { i } } I _ { n } ( \lambda a _ { i } ) = e ^ { \lambda ^ { \prime } a _ { i } } \hat { I } _ { n - 1 } ( \lambda a _ { i } ) - e ^ { \lambda ^ { \prime } a _ { i } } \frac { n } { \lambda a _ { i } } \hat { I } _ { n } ( \lambda a _ { i } ) = e ^ { \lambda ^ { \prime } a _ { i } } \hat { I } _ { n } ^ { \prime } ( \lambda a _ { i } ),$$

$$K _ { n } ^ { \prime } ( \lambda a _ { i } ) = - K _ { n - 1 } ( \lambda a _ { i } ) - \frac { \ddots \vdots } { \lambda a _ { i } } K _ { n } ( \lambda a _ { i } ) = - e ^ { - \lambda ^ { \prime } a _ { i } } \hat { K } _ { n - 1 } ( \lambda a _ { i } ) - e ^ { - \lambda ^ { \prime } a _ { i } } \frac { n } { \lambda a _ { i } } \hat { K } _ { n } ( \lambda a _ { i } ) = e ^ { - \lambda ^ { \prime } a _ { i } } \hat { K } _ { n } ^ { \prime } ( \lambda a _ { i } ).$$

If the argument is neither small nor large, rescaled modified-cylindrical functions are defined, in analogy to small and large arguments, as

$$I _ { n } ( \lambda a _ { i } ) = P _ { i } \hat { I } _ { n } ( \lambda a _ { i } ),$$

$$I _ { n } ^ { \prime } ( \lambda a _ { i } ) = P _ { i } \hat { I } _ { n } ^ { \prime } ( \lambda a _ { i } ),$$

$$K _ { n } ( \lambda a _ { i } ) & = P _ { i } ^ { - 1 } \hat { K } _ { n } ( \lambda a _ { i } ), & ( 3 5 c )$$

$$K _ { n } ^ { \prime } ( \lambda a _ { i } ) = P _ { i } ^ { - 1 } \hat { K } _ { n } ^ { \prime } ( \lambda a _ { i } ),$$

where the multiplicative factor P i is defined as in [16], i.e.,

$$\text{If } | I _ { n } ( \lambda a _ { i } ) | ^ { - 1 } < T _ { m }, \ \ P _ { i } = 1.$$

$$\text{If } | I _ { n } ( \lambda a _ { i } ) | ^ { - 1 } \geq T _ { m }, \ \ P _ { i } = | I _ { n } ( \lambda a _ { i } ) |.$$

with T m = 10 00 used for double-precision arithmetic computations. 1

and its subscript is linked to a i . As Table 1 shows, the argument for rescaled modified-cylindrical functions can be categorized into small, moderate, and large, with different and appropriate multiplicative factors defined accordingly.

The numerical threshold values used here to define small, moderate, and large argument ranges are identical to those used for the time-harmonic case involving cylindrical Bessel and Hankel functions detailed in [16] and not repeated here 1 .

## 2.4. Rescaled reflection and transmission coefficients

We can classify the multiplicative factors associated with the rescaled modified-cylindrical functions into two types, denoted as α and β , and shown in Table 2. The factor α i is associated with K n ( λa i ), whereas β i is associated with I n ( λa i ). Note again that the subscript i refers to to the index of the radial factor a i in the argument. There are two important properties to note for α and β : ( ) i Reciprocity : α i = 1 /β i and ( ii ) Boundness : | β α i j | ≤ 1 , for i &lt; j . As we will see below, these two properties are important in ensuring stable numerical computations.

1 A third type of threshold, used in connection with modified-Bessel functions magnitudes (not arguments) is also necessary within the moderate argument range. Again, this threshold value is identical to the one applied to ordinary Bessel functions in [16]

Table 1: Definition of rescaled modified-cylindrical functions for all types of arguments.

|                | small arguments                | moderate arguments       | large arguments              |
|----------------|--------------------------------|--------------------------|------------------------------|
| I n ( λa i )   | Ga n i ˆ I n ( λa i )          | P i ˆ I n ( λa i )       | e λ ′ a i ˆ I n ( λa i )     |
| I ′ n ( λa i ) | Ga n i ˆ I ′ n ( λa i )        | P i ˆ I ′ n ( λa i )     | e λ ′ a i ˆ I ′ n ( λa i )   |
| K n ( λa i )   | G - 1 a - n i ˆ K n ( λa i )   | P - 1 i ˆ K n ( λa i )   | e - λ ′ a i ˆ K n ( λa i )   |
| K ′ n ( λa i ) | G - 1 a - n i ˆ K ′ n ( λa i ) | P - 1 i ˆ K ′ n ( λa i ) | e - λ ′ a i ˆ K ′ n ( λa i ) |

Table 2: Definition of α i and β i .

| argument type   | α i           | β i       |
|-----------------|---------------|-----------|
| small           | G - 1 a - n i | Ga n i    |
| moderate        | P - 1 i       | P i       |
| large           | e - λ ′ a i   | e λ ′ a i |

Recalling the expressions for the reflection and transmission coefficients obtained before, the reflection coefficient R 12 for the outgoing-potential case is modified to

$$\text{or use outgoing-posteria case is motionive } \\ R _ { 1 2 } & = \frac { ( \sigma _ { 2 } - \sigma _ { 1 } ) K _ { n } ( \lambda a _ { 1 } ) K _ { n } ^ { \prime } ( \lambda a _ { 1 } ) } { \sigma _ { 1 } I _ { n } ^ { \prime } ( \lambda a _ { 1 } ) K _ { n } ( \lambda a _ { 1 } ) - \sigma _ { 2 } I _ { n } ( \lambda a _ { 1 } ) K _ { n } ^ { \prime } ( \lambda a _ { 1 } ) } \\ & = \alpha _ { 1 } ^ { 2 } \frac { ( \sigma _ { 2 } - \sigma _ { 1 } ) \hat { K } _ { n } ( \lambda a _ { 1 } ) \hat { K } _ { n } ^ { \prime } ( \lambda a _ { 1 } ) } { \sigma _ { 1 } \hat { I } _ { n } ^ { \prime } ( \lambda a _ { 1 } ) \hat { K } _ { n } ( \lambda a _ { 1 } ) - \sigma _ { 2 } \hat { I } _ { n } ( \lambda a _ { 1 } ) \hat { K } _ { n } ^ { \prime } ( \lambda a _ { 1 } ) } \\ & = \alpha _ { 1 } ^ { 2 } \hat { R } _ { 1 2 }.$$

Similarly, it can be shown that the reflection coefficient R 21 for the standing-potential case is modified to R 21 = β R 2 1 ˆ 21 , and that the transmission coefficient T 12 for the outgoing-potential case and the transmission coefficient T 21 for the standing-potential case simply recover the original ones without any multiplicative factors, i.e. T 12 = ˆ T 12 and T 21 = ˆ T 21 .

Based on the above modifications for the reflection and transmission coefficients, generalized reflection and transmission coefficients for thre or more layers can be similarly modified. After some algebra, it can be shown that ˜ R i,i +1 = α R 2 i ˆ ˜ i,i +1 and ˜ R i +1 ,i = β R 2 i ˆ ˜ i +1 ,i , and that, for both the outgoing-potential and standing-potential cases, ˜ T ij = ˆ ˜ T ij .

In addition to generalized reflection and transmission coefficients, the factors M j and N i ± considered before are also required to compute the potential. The basic difference between the two types of coefficients is that M j involves two generalized reflection coefficients whereas N i ± involves only one generalized reflection coefficient. All these auxiliary coefficients can be redefined accordingly using rescaled reflection coefficients, i.e.,

$$M _ { j } = \left [ 1 - \widetilde { R } _ { j, j - 1 } \widetilde { R } _ { j, j + 1 } \right ] ^ { - 1 } = \left [ 1 - \beta _ { j - 1 } ^ { 2 } \alpha _ { j } ^ { 2 } \hat { \widetilde { R } } _ { j, j - 1 } \hat { \widetilde { R } } _ { j, j + 1 } \right ] ^ { - 1 },$$

$$N _ { i + } & = \left [ 1 - R _ { i, i - 1 } \widetilde { R } _ { i, i + 1 } \right ] ^ { - 1 } = \left [ 1 - \beta _ { i - 1 } ^ { 2 } \alpha _ { i } ^ { 2 } \hat { R } _ { i, i - 1 } \hat { \widetilde { R } } _ { i, i + 1 } \right ] ^ { - 1 },$$

$$N _ { i - } = \left [ 1 - R _ { i, i + 1 } \widetilde { R } _ { i, i - 1 } \right ] ^ { - 1 } = \left [ 1 - \beta _ { i - 1 } ^ { 2 } \alpha _ { i } ^ { 2 } \hat { R } _ { i, i + 1 } \hat { \widetilde { R } } _ { i, i - 1 } \right ] ^ { - 1 }.$$

## 2.5. Rescaled integrand

Nest step is to modify the full integrand using rescaled modified-cylindrical functions. Since there are four integrand expressions, depending of the relative position of ρ and ρ ′ , each case is considered separately. For Case 1 , there are four radial parameters of interest: a j -1 , ρ ′ , ρ , and a j . For convenience, we let a j -1 = a 1 , ρ ′ = a 2 , ρ = a 3 , and a j = a 4 so that a 1 &lt; a 2 &lt; a 3 &lt; a 4 , and the integrand rewrites as

$$\begin{array} { c c c c c c } \, \cdot & \, \cdot & \, \cdot & \, \cdot & \, \cdot & \, \cdot \\ F _ { n } ( \rho, \rho ^ { \prime } ) & = \left [ K _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) \widetilde { R } _ { j, j + 1 } \right ] \left [ I _ { n } ( \lambda \rho ^ { \prime } ) + K _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j - 1 } \right ] M _ { j } \\ & = \left [ \beta _ { 2 } \alpha _ { 3 } \hat { K } _ { n } ( \lambda \rho ) + ( \beta _ { 2 } \alpha _ { 4 } ) ( \beta _ { 3 } \alpha _ { 4 } ) \hat { I } _ { n } ( \lambda \rho ) \hat { \widetilde { R } } _ { j, j + 1 } \right ] \left [ \hat { I } _ { n } ( \lambda \rho ^ { \prime } ) + ( \beta _ { 1 } \alpha _ { 2 } ) ^ { 2 } \hat { K } _ { n } ( \lambda \rho ^ { \prime } ) \hat { \widetilde { R } } _ { j, j - 1 } \right ] M _ { j } \\ & = \left [ A _ { 1 } \hat { K } _ { n } ( \lambda \rho ) + A _ { 2 } \hat { I } _ { n } ( \lambda \rho ) \hat { \widetilde { R } } _ { j, j + 1 } \right ] \left [ A _ { 3 } \hat { I } _ { n } ( \lambda \rho ^ { \prime } ) + A _ { 4 } \hat { K } _ { n } ( \lambda \rho ^ { \prime } ) \hat { \widetilde { R } } _ { j, j - 1 } \right ] M _ { j }, \end{array}$$

where β -1 2 = α 2 has been used. All multiplicative factors A 1 , A 2 , A 3 , A 4 have magnitudes no larger than one due to the boundness property discussed above.

Similarly, for Case 2 , there are four radial parameters of interest: a j -1 , ρ , ρ ′ , and a j . For convenience, we let a j -1 = a 1 , ρ = a 2 , ρ ′ = a 3 , and a j = a 4 so that a 1 &lt; a 2 &lt; a 3 &lt; a 4 , and the integrand rewrites as

$$\L _ { n } ^ { \ } ( \rho, \rho ^ { \prime } ) & = \left [ I _ { n } ( \lambda \rho ) + K _ { n } ( \lambda \rho ) \widetilde { R } _ { j, j - 1 } \right ] \left [ K _ { n } ( \lambda \rho ^ { \prime } ) + I _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j + 1 } \right ] M _ { j } \\ & = \left [ \beta _ { 2 } \alpha _ { 3 } \hat { I } _ { n } ( \lambda \rho ) + ( \beta _ { 1 } \alpha _ { 2 } ) ( \beta _ { 1 } \alpha _ { 3 } ) \hat { K } _ { n } ( \lambda \rho ) \hat { \widetilde { R } } _ { j, j - 1 } \right ] \left [ \hat { K } _ { n } ( \lambda \rho ^ { \prime } ) + ( \beta _ { 3 } \alpha _ { 4 } ) ^ { 2 } \hat { I } _ { n } ( \lambda \rho ^ { \prime } ) \hat { \widetilde { R } } _ { j, j + 1 } \right ] M _ { j } \\ & = \left [ B _ { 1 } \hat { I } _ { n } ( \lambda \rho ) + B _ { 2 } \hat { K } _ { n } ( \lambda \rho ) \hat { \widetilde { R } } _ { j, j - 1 } \right ] \left [ B _ { 3 } \hat { K } _ { n } ( \lambda \rho ^ { \prime } ) + B _ { 4 } \hat { I } _ { n } ( \lambda \rho ^ { \prime } ) \hat { \widetilde { R } } _ { j, j + 1 } \right ] M _ { j }.$$

Again, all multiplicative factors B 1 , B 2 , B 3 , B 4 magnitudes are bounded by one.

For Case 3 , there are six radial parameters of interest: a j -1 , ρ ′ , a j , a i -1 , ρ , and a i . For convenience, we let a j -1 = a 1 , ρ ′ = a 2 , a j = a 3 , a i -1 = a 4 , ρ = a 5 , and a i = a 6 so that so that a 1 &lt; a 2 &lt; a 3 ≤ a 4 &lt; a 5 &lt; a 6 , and the integrand rewrites as

$$& \Lambda _ { n } ( \rho, \rho ^ { \prime } ) = \left [ K _ { n } ( \lambda \rho ) + I _ { n } ( \lambda \rho ) \widetilde { R } _ { i, i + 1 } \right ] \left [ I _ { n } ( \lambda \rho ^ { \prime } ) + K _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j - 1 } \right ] N _ { i + \widetilde { T } _ { j i } M _ { j } } \\ & \quad = \left [ \beta _ { 2 } \alpha _ { 5 } \hat { K } _ { n } ( \lambda \rho ) + ( \beta _ { 2 } \alpha _ { 6 } ) ( \beta _ { 5 } \alpha _ { 6 } ) \hat { I } _ { n } ( \lambda \rho ) \hat { \widetilde { R } } _ { i, i + 1 } \right ] \left [ \hat { I } _ { n } ( \lambda \rho ^ { \prime } ) + ( \beta _ { 1 } \alpha _ { 2 } ) ^ { 2 } \hat { K } _ { n } ( \lambda \rho ^ { \prime } ) \hat { \widetilde { R } } _ { j, j - 1 } \right ] N _ { i + \widetilde { T } _ { j i } M _ { j } } \\ & \quad = \left [ C _ { 1 } \hat { K } _ { n } ( \lambda \rho ) + C _ { 2 } \hat { I } _ { n } ( \lambda \rho ) \hat { \widetilde { R } } _ { i, i + 1 } \right ] \left [ C _ { 3 } \hat { I } _ { n } ( \lambda \rho ^ { \prime } ) + C _ { 4 } \hat { K } _ { n } ( \lambda \rho ^ { \prime } ) \hat { \widetilde { R } } _ { j, j - 1 } \right ] N _ { i + \widetilde { T } _ { j i } M _ { j } }.$$

All multiplicative factors C 1 , C 2 , C 3 , C 4 have magnitudes never greater than unity.

Finally, for Case 4 , there are again six radial parameters of interest: a i -1 , ρ , a i , a j -1 , ρ ′ , and a j . For convenience, we let a i -1 = a 1 , ρ = a 2 , a i = a 3 , a j -1 = a 4 , ρ ′ = a 5 , and a j = a 6 so that a 1 &lt; a 2 &lt; a 3 ≤ a 4 &lt; a 5 &lt; a 6 , and the integrand rewrites as

$$\begin{array} { c c c } \cdot & \quad \cdot & \quad \cdot \\ F _ { n } ( \rho, \rho ^ { \prime } ) & = \left [ I _ { n } ( \lambda \rho ) + K _ { n } ( \lambda \rho ) \widetilde { R } _ { i, i - 1 } \right ] \left [ K _ { n } ( \lambda \rho ^ { \prime } ) + I _ { n } ( \lambda \rho ^ { \prime } ) \widetilde { R } _ { j, j + 1 } \right ] N _ { i - } \widetilde { T } _ { j i } M _ { j } \\ & = \left [ \beta _ { 2 } \alpha _ { 5 } \hat { I } _ { n } ( \lambda \rho ) + ( \beta _ { 1 } \alpha _ { 2 } ) ( \beta _ { 1 } \alpha _ { 5 } ) \hat { K } _ { n } ( \lambda \rho ) \hat { \widetilde { R } } _ { i, i - 1 } \right ] \left [ \hat { K } _ { n } ( \lambda \rho ^ { \prime } ) + ( \beta _ { 5 } \alpha _ { 6 } ) ^ { 2 } \hat { I } _ { n } ( \lambda \rho ^ { \prime } ) \hat { \widetilde { R } } _ { j, j + 1 } \right ] N _ { i - } \widetilde { T } _ { j i } M _ { j } \\ & = \left [ D _ { 1 } \hat { I } _ { n } ( \lambda \rho ) + D _ { 2 } \hat { K } _ { n } ( \lambda \rho ) \hat { \widetilde { R } } _ { i, i - 1 } \right ] \left [ D _ { 3 } \hat { K } _ { n } ( \lambda \rho ^ { \prime } ) + D _ { 4 } \hat { I } _ { n } ( \lambda \rho ^ { \prime } ) \hat { \widetilde { R } } _ { j, j + 1 } \right ] N _ { i - } \widetilde { T } _ { j i } M _ { j }. \end{array}$$

Once more, all multiplicative factors D 1 , D 2 , D 3 , D 4 have magnitudes never greater than one.

## 2.6. Numerical integration

The electric potential expression includes a semi-infinite integral and an infinite series summation, as shown in (28). Therefore, truncation errors are inevitable and an error analysis should be made to ensure reliable results. In addition, the involved Sommerfeld-type integrals can be notoriously slowly convergent. A variety of extrapolation methods are applied here in order to accelerate the numerical integration, as described in more detail in Appendix B.

## 3. Results

A number of cases of practical interest for geophysical exploration are considered next. The results below are obtained on 2.6 GHz Opteron with 8 cores and 32 GB memory. For all cases below, the relative permittivity glyph[epsilon1] r and relative permeability µ r are set to one. The innermost layer has 6 ′′ radius and represents a mud-filled borehole. The outer layers represent the adjacent Earth formations, invasion zones, and/or casing layers. Each layer assumes different resistivity values R , where R = 1 /σ . As Fig. 6 illustrates, a current source emitting a DC current of 1 A is located at the Survey Electrode position. The electric potential is measured at two points: the Measurement Electrode 1 ( V 16 ′′ ) positioned 16 inches away from the source along the vertical direction and the Measurement Electrode 2 ( V 32 ′′ ) positioned 32 inches away from the source along the vertical direction. Both the source and two measurement electrodes are spaced 5 inches away from the z -axis and their azimuthal positions are the same ( φ = φ ′ = 0). The relevant error parameters as defined in the Appendix and used in the following are e tol = 10 -4 and e thr = 10 -4 . The electric potential and resistivity units used here are [ V ] and [Ω · m ], respectively. The results of the present algorithm are compared against those from the finite element method (FEM).

## 3.1. Logging simulation results

Case 1 corresponds to a homogeneous problem, where analytical solutions are available. As Table 3 shows, the new algorithm produces very accurate results with very fast computing time. This is very important for making feasible the solution of the inverse problem (i.e., determining the unknown resistivity profile from known electric potential at the received electrodes) using iterative methods predicated on repeated forward solves. Homogeneous Formation Results

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-ddc65e2e-ef0c-4b30-9a12-48d982319f6e.jpg)

(Rm =  1 Ω-m , Rf =  1 Ω-m)









Cross-section of the formation in ρ-z plane is shown.

Survey electrode transmits a DC

current of 1A.

Both mud resistivity and formation resistivity are 1 Ω-m.

Voltages at measure electrodes and voltage difference between measure

electrodes for this case are:

V\_ME1 [V]

0.17878839

V\_ME2 [V]

∆V [V]

1

Figure 6: Case 1 in the ρz -plane with R 1 = R 2 = 1.

Table 3: Comparison of electric potential for Case 1.

|                 | Analytical                      | FEM                             | Present Algorithm                                 |
|-----------------|---------------------------------|---------------------------------|---------------------------------------------------|
| V 16 ′′ V 32 ′′ | 1.9581 × 10 - 1 9.7905 × 10 - 2 | 1.7878 × 10 - 1 9.4416 × 10 - 2 | 1.9580 × 10 - 1 (2 sec.) 9.7900 × 10 - 2 (2 sec.) |
| ∆ V             | 9.7905 × 10 - 2                 | 8.4371 × 10 - 2                 |                                                   |
|                 |                                 |                                 | 9.7902 × 10 - 2                                   |

0.094416797

0.084371591

Cases 2 and 3 correspond to a two-layer problem with 1:5 contrast between the adjacent borehole and formation resistivities.

with One Radial Layer with One Radial Lay

(Rm =  1 Ω-m , Rf =  5 Ω-m)

Mud resistivity is Rm = 1 Ω-m

Formation resistivity is Rf = 5 Ω- m.

Voltages at measure electrodes and voltage difference between

measure electrodes for this case are:

V\_ME1 [V]

V\_ME2 [V]

∆V [V]

(Rm =  5 Ω-m , Rf =  1 Ω-

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-36927dc6-aa73-4081-bb4e-e811c289a8a9.jpg)

Mud resistivity is Rm = 5 Ω-m

Formation resistivity is Rf = 1 Ω- m.

Voltages at measure electrodes and voltage difference between

measure electrodes for this case are:

V\_ME1 [V]

V\_ME2 [V]

∆V [V]

0.16945594

0.090310919

0.079145021

4

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-52cb85dd-1064-49d9-ba10-c07a0fdfb883.jpg)

5

Figure 7: Cases 2 and 3 in the ρz -plane: (a) Case 2 with R 1 = 1, R 2 = 5 and (b) Case 3 with R 1 = 5, R 2 = 1.

Table 4: Comparison of electric potential for Case 2.

|         | FEM             | Present Algorithm        |
|---------|-----------------|--------------------------|
| V 16 ′′ | 9.5544 × 10 - 1 | 9.7802 × 10 - 1 (2 sec.) |
| V 32 ′′ | 5.4241 × 10 - 1 | 5.4981 × 10 - 1 (2 sec.) |
| ∆ V     | 4.1303 × 10 - 1 | 4.2822 × 10 - 1          |

Table 5: Comparison of electric potential for Case 3.

|                 | FEM                         | Present Algorithm        |
|-----------------|-----------------------------|--------------------------|
| V 16 ′′ V 32 ′′ | 1.6945 × 10 - 1 - 2         | 2.0533 × 10 - 1 (2 sec.) |
| ∆ V             | 9.0310 × 10 7.9145 × 10 - 2 | 9.7677 × 10 - 2 (2 sec.) |
|                 |                             | 1.0766 × 10 - 1          |









Cases 4 and 5 include a highly conductive casing. As Table 7 shows, there is disagreement between the semi-analytical and FEM results w.r.t. the absolute value of the electric potentials, although the computed potential differences (voltage drop) between the two electrodes are very similar. This disagreement is easy to explain as a FEM mesh truncation effect. The electric current in this case flows primarily along the thin, highly conductive casing, which does not produce sufficient decay on the electric potential before before it reaches the mesh boundary. Consequently, the FEM result has a spurious potential offset. The FEM has difficulty in simulating this problem unless a mesh truncation treatment is included on both top and bottom boundaries and/or a very long mesh is used along the z direction. Notwithstanding such discrepancy, note that, in a resistivity logging context, the primary quantity of interest is the potential difference between electrodes.

Simulation Results with Casing

Simulation Results with Casing and Formation

(Rm =  1 Ω-m , Rc =  10 -8  Ω-m)



In this case, formation is replaced with (an infinite) casing.

Mud resistivity is Rm = 1 Ω-m

Casing resistivity is Rc = 10 -8  Ω- Electrode 2

m.

Voltages at measure electrodes  Electrode 1

and voltage difference between 16'

measure electrodes for this case are:

V\_ME1 [V]

V\_ME2 [V]

∆V [V]

(Rm =  1 Ω-m , Rc =  10 -8  Ω-m, Rf  = 2 Ω-m )

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-0b3d4544-0789-42d9-a2fd-af6b1b680e65.jpg)

In this simulation, casing width is shrinked to an eight of an inch and a

formation layer that fills the rest of the computational domain is added.

Mud resistivity is Rm = 1 Ω-m

Casing resistivity is Rc = 10 -8  Ω-m.

Formation resistivity is Rf = 2 Ω- m.

Voltages at measure electrodes and voltage difference between measure

electrodes for this case are:

V\_ME1 [V]

0.22658744E-03

V\_ME2 [V]

∆V [V]

7

6

Figure 8: Cases 4 and 5 in the ρz -plane: (a) Case 4 with R 1 = 1, R 2 = 10 -8 and (b) Case 5 with R 1 = 1, R 2 = 10 -8 , R 3 = 2.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-961d631d-22f2-4db0-a853-0cfe9ca6bbbb.jpg)

Table 6: Comparison of electric potential for Case 4.

|         | FEM             | Present Algorithm        |
|---------|-----------------|--------------------------|
| V 16 ′′ | 1.2980 × 10 - 4 | 1.3873 × 10 - 4 (2 sec.) |
| V 32 ′′ | 2.1372 × 10 - 7 | 2.1415 × 10 - 7 (2 sec.) |
| ∆ V     | 1.2959 × 10 - 4 | 1.3852 × 10 - 4          |

Table 7: Comparison of electric potential for Case 5.

|         | FEM             | Present Algorithm        |
|---------|-----------------|--------------------------|
| V 16 ′′ | 2.2658 × 10 - 4 | 1.7241 × 10 - 3 (3 sec.) |
| V 32 ′′ | 9.6331 × 10 - 5 | 1.5885 × 10 - 3 (2 sec.) |
| ∆ V     | 1.3025 × 10 - 4 | 1.3562 × 10 - 4          |

0.96331022E-04

0.13025642E-03











As Fig. 9 illustrates, Cases 6 and 7 feature two layers besides the borehole, where the mud layer can represent an invasion zone with a resistivity between those of the borehole and the outer formation.

Simulation Results with Two Formation Layers

Simulation Results with Two Formation Laye

(Rm =  1 Ω-m , Rf1 = 1.5 Ω-m, Rf2  = 2 Ω-m )

(Rm =  1 Ω-m , Rf1 = 1.5 Ω-m, Rf2  = 2 Ω-

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-789d65a0-3093-4b80-ad44-686a1ba7e49d.jpg)

In this simulation, casing in the previous slide is replaced with a resistive layer of

resistivity Rf1 = 1.5 Ω-m.

Mud resistivity is Rm = 1 Ω-m         Measure

Resistivity of second layer is Rf1 = 1.5 Ω-  e 2

m.                                                   16'

Resistivity of second layer is Rf2 = 2

Ω- m.

Voltages at measure electrodes and voltage difference between measure  urvey

electrodes for this case are:

V\_ME1 [V]

V\_ME2 [V]

∆V [V]

8

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-1c5fccfe-f69a-4ba3-8f31-1d3b5a3800fe.jpg)

9

Figure 9: Cases 6 and 7 in the ρz -plane: (a) Case 6 with R 1 = 1, R 2 = 1 5, . R 3 = 2 and (b) Case 7 with R 1 = 1, R 2 = 1 5, . R 3 = 2.

Table 8: Comparison of electric potential for Case 6.

|         | FEM             | Present Algorithm        |
|---------|-----------------|--------------------------|
| V 16 ′′ | 3.7097 × 10 - 1 | 3.9061 × 10 - 1 (2 sec.) |
| V 32 ′′ | 1.9796 × 10 - 1 | 2.0245 × 10 - 1 (2 sec.) |
| ∆ V     | 1.7301 × 10 - 1 | 1.8816 × 10 - 1          |

Table 9: Comparison of electric potential for Case 7.

|         | FEM             | Present Algorithm        |
|---------|-----------------|--------------------------|
| V 16 ′′ | 3.6077 × 10 - 1 | 3.8160 × 10 - 1 (2 sec.) |
| V 32 ′′ | 1.9694 × 10 - 1 | 2.0156 × 10 - 1 (2 sec.) |
| ∆ V     | 1.6383 × 10 - 1 | 1.8003 × 10 - 1          |

## 3.2. Electric potential maps

Plots of the spatial distributions of the electric potential of Cases 2-5 are provided next. Since the potential varies by many orders of magnitude near the electrodes, a log-scale is used; i.e. we plot the quantity 10 log 10 | V | . Fig. 10 depicts the potential distributions using a three-dimensional view. Fig. 11 and Fig. 12 show the cross-sections of the potential distributions at the z = 16 ′′ and z = 32 ′′ planes, respectively, while Fig. 13 shows the cross-sections at the y = 0 ′′ plane. In all figures, thicker black lines represent interfaces between cylindrical layers and thinner black lines represent equipotential contours. Note that the third layer of Case 5 is too thin for visualization in these plots.



In this simulation, width of the first formation layer is expanded to 3'.

Mud resistivity is Rm = 1 Ω-m

Resistivity of second layer is Rf1 = 1.5 Ω- m.

Resistivity of second layer is Rf2 = 2

Ω- m.

Voltages at measure electrodes and voltage difference between measure

electrodes for this case are:

V\_ME1 [V]

0.36077656

V\_ME2 [V]

∆V [V]

0.19694043

0.16383612











![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-1c611ee3-add4-48e9-855f-c528c6bd06cf.jpg)

(a)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-ab8fa43d-a74a-41cd-8b19-e9beac7f8060.jpg)

(b)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-98d4ffc1-1d96-4854-817f-e0c901edde9b.jpg)

Figure 10: Electric potential around the current electrode: (a) Case 2, (b) Case 3, (c) Case 4, and (d) Case 5.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-f056dc48-63f5-4efe-9f3e-c551d796f9c3.jpg)

Case5

Top : 0.

Botto

Left : 0

Right :

Case7-

Top : 0.

Botto

Left : 0

Right :

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-ca3a5a73-f26c-4e07-9340-9d951f36b566.jpg)

Figure 11: Electric potential at z = 16 ′′ plane around the current electrode: (a) Case 2, (b) Case 3, (c) Case 4, and (d) Case 5.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-88332121-9ac4-48c5-a146-6f25e9a8937a.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-e8cad186-7184-41aa-bdff-8061e0c2f8d9.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-17f51bc4-9ed2-4393-b6a1-44feba960c57.jpg)

Case5 - z

Top : 0.2

Bottom :

Left : 0.0

Right : 4.

Case7 - z16

Top : 0.2

Bottom : 3.

Left : 0.0

Right : 4.6 o

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-c78f02e0-02e4-452b-9e2c-8f0e349bd4af.jpg)

Figure 12: Electric potential at z = 32 ′′ plane around the current electrode: (a) Case 2, (b) Case 3, (c) Case 4, and (d) Case 5.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-8da9ecdd-9f29-4dce-a188-c4722f4f41e6.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-c41af8a1-3688-4562-a269-90b8ae68bfd2.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-c3e0c711-b3fb-4c62-8777-4f2be30955af.jpg)

Case5 - z

Top : 0.2

Bottom :

Left : 0.0

Right : 4.

Case7 - z32

Top : 0.2

Bottom : 3.

Left : 0.0

Right : 4.6 o

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-5df0496d-1e79-4f05-be7d-a63e83618e74.jpg)

Figure 13: Electric potential at y = 0 ′′ plane around the current electrode: (a) Case 2, (b) Case 3, (c) Case 4, and (d) Case 5.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-835cc03c-0e4a-46c9-934d-61f54178e79a.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-65105b39-33a9-4c88-99ac-d115588fbadf.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-e376050b-1f5e-4f4a-a996-69c4b1a8a1aa.jpg)

Case5 - y

Top : 0.2

Bottom :

Left : 0.0

Right : 4.

Case7 - y

Top : 0.2

Bottom :

Left : 0.0

Right : 4.

## 4. Conclusion

We have introduced a stable semianalytical formulation to compute the potential due to arbitrarylocated electrodes in a cylindrically layered media, where the physical parameters (viz., layer resistivity) can vary by many order of magnitude. Stability was achieved by a rescaling of the modified Bessel functions and subsequent manipulation of the integrand expressions to avoid underflow and overflow problems under double-precision arithmetic. Several extrapolation methods to compute the involved Sommerfeld integrals were also considered and evaluated. The resulting algorithm was verified in a number of cases relevant to borehole geophysics. The fast speed and robustness of the algorithm makes it also quite suited as a forward solver engine for the inverse roblem, where the formation resistivity needs to be estimated from a limited numberof computed (measured) potential data. The same basic algorithm can be applied to other steadystate diffusion problems obeying Poisson's equation with discontinous coefficients in cylindrically layered geometries.

## Acknowledgments

The authors are grateful to Halliburton Energy Services for the permission to publish this work and to Dr. Baris Guner for generating validation data.

## Appendix A: Folding azimuth summation

The computation of potential in (28) involves infinite series over the azimuth mode n . For an arbitrary order ν ≥ 0 [34], we have

$$I _ { - \nu } ( z ) & = I _ { \nu } ( z ) + \frac { 2 } { \pi } \sin ( \nu \pi ) K _ { \nu } ( z ), & ( 4 3 \tt a ) \\ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

$$K _ { - \nu } ( z ) & = K _ { \nu } ( z ).$$

If the order is a positive integer, ν = n , the above reduces to I -n ( z ) = I n ( z ) and K -n ( z ) = K n ( z ). Consequently, we can write

$$\cdot \quad \vec { \psi } _ { i } & = \frac { \mathcal { I } } { 2 \pi ^ { 2 } \sigma _ { j } } \sum _ { n = - \infty } ^ { \infty } e ^ { i n ( \phi - \phi ^ { \prime } ) } \int _ { 0 } ^ { \infty } F _ { n } ( \rho, \rho ^ { \prime } ) \cos ( \lambda ( z - z ^ { \prime } ) ) d \lambda \\ & = \frac { \mathcal { I } } { 2 \pi ^ { 2 } \sigma _ { j } } \int _ { 0 } ^ { \infty } \left [ F _ { 0 } ( \rho, \rho ^ { \prime } ) + 2 \sum _ { n = 1 } ^ { \infty } F _ { n } ( \rho, \rho ^ { \prime } ) \cos ( n ( \phi - \phi ^ { \prime } ) ) \right ] \cos ( \lambda ( z - z ^ { \prime } ) ) d \lambda.$$

## Appendix B: Extrapolation methods and convergence study

Among numerous extrapolation methods, some popular ones for Sommerfeld-type integrals are very briefly revisited here. For more details, the reader can refer, e.g., to [15, 33]. Before a given sequence is extrapolated, a Sommerfeld-type integral can be divided into a number of subintervals as

$$S = \int _ { a } ^ { \infty } g ( \lambda ) p ( \lambda ) d \lambda = \sum _ { i = 0 } ^ { \infty } \int _ { \lambda _ { i - 1 } } ^ { \lambda _ { i } } g ( \lambda ) p ( \lambda ) d \lambda = \sum _ { i = 0 } ^ { \infty } u _ { i },$$

where g λ ( ) is an exponentially decaying part and p λ ( ) is an oscillatory part. This is called the partitionextrapolation approach [28]. In general, the remainders are defined as

$$r _ { n } = S _ { n } - S = - \int _ { \lambda _ { n } } ^ { \infty } g ( \lambda ) p ( \lambda ) d \lambda.$$

Furthermore, the remainders are assumed to feature Poincar´-type asymptotic expansions [15] of the form e

$$r _ { n } \sim \omega _ { n } \sum _ { i = 0 } ^ { \infty } a _ { i } \lambda _ { n } ^ { - i }, \quad n \to \infty,$$

where ω n is the remainder estimate and a i are associated coefficients. The estimates ω n play an important role in the extrapolation and can be obtained analytically or numerically. The coefficients a i are unknowns but they are not necessary for the extrapolation itself. In our case, g λ ( ) in (45) can be asymptotically expressed as

$$g ( \lambda ) \sim \frac { e ^ { - \lambda | \rho - \rho ^ { \prime } | } } { \lambda } \sum _ { i = 0 } ^ { \infty } \frac { c _ { i } } { \lambda ^ { i } },$$

where c i are arbitrary constants. Furthermore, p λ ( ) = cos( λ z ( -z ′ )) with half-period equal to π/ ρ | -ρ ′ | . After some algebra, it can be shown that the remainder estimates write as

$$\omega _ { n } = \frac { ( - 1 ) ^ { n + 1 } } { \lambda _ { n } } e ^ { - \frac { n \pi | \rho - \rho ^ { \prime } | } { | \varepsilon - \varepsilon ^ { \prime } | } }.$$

We list next three popular extrapolation methods for a given sequence, { S n } = S , S 0 1 , S 2 , · · · , S n , where we consider S n = S (0) n .

## (i) Euler transformation

$$S _ { n } ^ { ( k + 1 ) } = \frac { 1 } { 2 } \left ( S _ { n } ^ { ( k ) } + S _ { n + 1 } ^ { ( k ) } \right ), \quad n, k \geq 0$$

The best approximation in this case is S ( k ) 0 | k = n , and this choice is most effective for logarithmic alternating sequences.

## (ii) Iterative Aitken transformation

$$S _ { n } ^ { ( k + 1 ) } = S _ { n } ^ { ( k ) } - \frac { [ \Delta S _ { n } ^ { ( k ) } ] ^ { 2 } } { \Delta ^ { 2 } S _ { n } ^ { ( k ) } }, \quad n, k \geq 0$$

Note that ∆ S ( k ) n = S ( k ) n +1 -S ( k ) n and ∆ 2 S ( k ) n = S ( k ) n +2 -2 S ( k ) n +1 + S ( k ) n . Obviously, this is a nonlinear transformation and it can be applied to both linear monotone and alternating sequences. When an odd number of sequences is given, the best approximation is S ( k ) 0 | k = n . For an even number, the best approximation is S ( k ) 1 | k = n .

## (iii) Weighted-averages method

$$S _ { n } ^ { ( k + 1 ) } = \frac { S _ { n } ^ { ( k ) } + \eta _ { n } ^ { ( k ) } S _ { n + 1 } ^ { ( k ) } } { 1 + \eta _ { n } ^ { ( k ) } }, \quad n, k \geq 0,$$

where η ( k ) is the weight and defined as

$$\eta ^ { ( k ) } = - \frac { \omega _ { n } } { \omega _ { n + 1 } } = \frac { \lambda _ { n + 1 } } { \lambda _ { n } } e ^ { \frac { \pi | \rho - \rho ^ { \prime } | } { | s - z ^ { \prime } | } }.$$

The weighted-averages method [18, 19] can be regarded as generalized Euler transformation since η = 1 recovers the Euler transformation. The relative power of this method, compared to the two other methods, comes from using remainder estimates. As pointed out in [15], the weighted-averages method is quited suited for extrapolating Sommerfeld-type integrals.

Figure 14: Two scenarios of the integrand: (a) Type 1 with | ρ -ρ ′ | = 0 05, . | z -z ′ | = 0 01 and (b) Type 2 . with | ρ -ρ ′ | = 0 05, . | z -z ′ | = 0 5. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-bef5710b-e114-4d83-bba4-f64988a4967e.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-9eb4daa9-6a64-40be-9b37-ac6f80033434.jpg)

To evaluate the integral (45), the subinterval length q should be first determined. As suggested in [15], the half-period of the oscillating part of the integrand is a good choice for q because it makes the sequences alternating; i.e., q = π/ z | -z ′ | . However, this is not appropriate here in some circumstances. Let us consider the two scenarios depicted in Fig. 14. The two figures show the behavior of g λ ( ) and p λ ( ) for different combinations of | ρ -ρ ′ | and | z -z ′ | as λ increases. For better visualization, the functions are normalized by their respective maxima. Type 1 integrand occurs when | ρ -ρ ′ | &gt; z | -z ′ | and shows linear monotone convergence. On the other hand, Type 2 integrand occurs when | ρ -ρ ′ | &lt; z | -z ′ | and shows logarithmic alternating convergence. Therefore, such an expression for q is not appropriate for Type 1 because the integrand is near zero before the first half-period comes. For the Type 1 integrand, a different subinterval length can be defined as q = π/ ρ | -ρ ′ | , so that the subinterval length can be written in general as

$$q = \frac { \pi } { \max ( | \rho - \rho ^ { \prime } |, | z - z ^ { \prime } | ) }.$$

Once q is determined, it becomes necessary to determine how many subintervals are required to reach sufficient convergence. To do so, the relative error below is defined

$$e _ { i } = \frac { | T ( S _ { i + 1 } ) - T ( S _ { i } ) | } { | T ( S _ { i + 1 } ) | },$$

where T S ( i ) is the extrapolated (transformed) value for given sequences, { S i } = S , S 0 1 , · · · , S i . If e i is less than a given error tolerance e tol , the sequence is stopped and the number of subintervals is determined. Next, the number of quadrature points per the subinterval is increased until the relative error between two adjacent iterations meets the desired criterion. To distinguish it from e tol , the criterion in this step is denoted error threshold e thr . The same step is repeated for determining the number of orders. The overall procedure is schematically depicted in Fig. 15.

Some convergence tests are performed next to evaluate the numerical integration procedures above. The domain is assumed homogeneous because exact (analytical) solutions are available as reference. Three cases of source/observation distances are considered.

Case 1 : | ρ -ρ ′ | = 0 001 m, . | z -z ′ | = 0 1 m .

$$\text{Case 2} \, \colon \, | \rho - \rho ^ { \prime } | = \ 0. 1 \, m, \ | z - z ^ { \prime } | = 0. 1 \, m$$

$$\text{Case } 3 \, \colon \, \stackrel { \dots } { | \rho - \dots } \rho ^ { \prime } | = \ \ 1 0 \, m, \ \ | z - z ^ { \prime } | = 0. 1 \, m$$

Figure 15: Flowchart of the procedure for computation of electric potential using extrapolation methods.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0316v1-figure-f34c6003-07f2-4d62-854a-e6983e8faf4c.jpg)

In each case, the three extrapolation methods mentioned before are compared. The excitation current magnitude and the medium resistivity are both set to one. The relevant error parameters are e tol = 10 -6 and e thr = 10 -4 . The smaller e tol is chosen to examine the effect of the methods on the number of subintervals. Tables 10, 11, and 12 compare the results, where the first rows represent the number of subintervals needed to achieve convergence in terms of e tol and the second rows represent the relative error against the analytical solution. As Table 10 shows, the iterative Aitken method does not work for Case 1 because it corresponds to a Type 2 integrand with logarithmic alternating convergence. The Euler transformation works well for all cases, but it can be seen that the weighted-averages method provides best results, corroborating the conclusions stated in [15].

Table 10: Comparison of the three extrapolation methods for Case 1.

|                   | Euler           | Aitken   | weighted-averages   |
|-------------------|-----------------|----------|---------------------|
| # of subintervals | 16              | > 100    | 10                  |
| Relative error    | 1.9990 × 10 - 6 | N. A.    | 3.5863 × 10 - 7     |

Table 11: Comparison of the three extrapolation methods for Case 2.

|                   | Euler           | Aitken          | weighted-averages   |
|-------------------|-----------------|-----------------|---------------------|
| # of subintervals | 13              | 13              | 4                   |
| Relative error    | 1.4131 × 10 - 6 | 1.1921 × 10 - 6 | 1.4058 × 10 - 6     |

Table 12: Comparison of the three extrapolation methods for Case 3.

|                   | Euler           | Aitken          | weighted-averages   |
|-------------------|-----------------|-----------------|---------------------|
| # of subintervals | 13              | 13              | 5                   |
| Relative error    | 4.9227 × 10 - 7 | 4.9227 × 10 - 7 | 4.9799 × 10 - 7     |

## References

- [1] Michael Carley. Numerical solution of the modified Bessel equation. IMA Journal of Numerical Analysis , 33(3):1048-1062, Jul. 2013.
- [2] Weng Cho Chew. Waves and Fields in Inhomogeneous Media . IEEE Press Series on Electromagnetic Waves. IEEE Press, New York, 1995.
- [3] J. A. Doetsch, I. Coscia, S. Greenhalgh, N. Linde, A. Green, and T. G¨nther. u The borehole-fluid effect in electrical resistivity imaging. Geophysics , 75(4):F107-F114, Jul.-Aug. 2010.
- [4] D. Drahos. Electrical modeling of the inhomogeneous invaded zone. Geophysics , 49(10):1580-1585, Oct. 1984.
- [5] Darwin V. Ellis and Julian M. Singer. Well Logging for Earth Scientists . Springer, Dordrecht, The Netherlands, second edition, 2007.
- [6] Guo-Xin Fan, Qing Huo Liu, and Sean P. Blanchard. 3-D numerical mode-matching (NMM) method for resistivity well-logging tools. IEEE Transactions on Antennas and Propagation , 48(10):1544-1552, Oct. 2000.
- [7] S. Gianzero and B. Anderson. An integral transform solution to the fundamental problem in resistivity logging. Geophysics , 47(6):946-956, 1982.
- [8] Y. K. Hue, F. L. Teixeira, L. San Martin, and M. S. Bittar. Three-dimensional simulation of eccentric LWD tool response in boreholes through dipping formations. IEEE Transactions on Geoscience and Remote Sensing , 43(2):257-268, Feb. 2005.
- [9] Yik-Kiong Hue and F. L. Teixeira. Numerical mode-matching method for tilted-coil antennas in cylindrically layered anisotropic media with multiple horizontal beds. IEEE Transactions on Geoscience and Remote Sensing , 45(8):24512462, Aug. 2007.
- [10] Hwa Ok Lee, Fernando L Teixeira, Luis E San Martin, and Michael S Bittar. Numerical modeling of eccentered LWD borehole sensors in dipping and fully anisotropic earth formations. IEEE Transactions on Geoscience and Remote Sensing , 50(3):727-735, Mar. 2012.
- [11] Qing-Huo Liu, Barbara Anderson, and Weng Cho Chew. Modeling low-frequency electrode-type resistivity tools in invaded thin beds. IEEE Transactions on Geoscience and Remote Sensing , 32(3):494-498, May 1994.
- [12] Sixin Liu and Motoyuki Sato. Electromagnetic logging technique based on borehole radar. IEEE Transactions on Geoscience and Remote Sensing , 40(9):2083-2092, Sep 2002.
- [13] J. R. Lovell and W. C. Chew. Effect of tool eccentricity on some electrical well-logging tools. IEEE Transactions on Geoscience and Remote Sensing , 28(1):127-136, Jan. 1990.
- [14] John Richard Lovell. Finite element methods in resistivity logging . Ph.d. dissertation, Delft University, Delft, The Netherlands, 1993.
- [15] K. A. Michalski. Extrapolation methods for Sommerfeld integral tails. IEEE Transactions on Antennas and Propagation , 46(10):1405-1418, Oct. 1998.
- [16] H. Moon, F. L. Teixeira, and B. Donderici. Stable pseudoanalytical computation of electromagnetic fields from arbitrarilyoriented dipoles in cylindrically stratified media. Journal of Computational Physics , 273:118-142, 2014.
- [17] J. H. Moran and S. Gianzero. Effects of formation anisotropy on resistivity-logging measurements. Geophysics , 44(7): 1266-1286, Jul. 1979.
- [18] JR Mosig and FE Gardiol. Analytical and numerical techniques in the Green's function treatment of microstrip antennas and scatterers. In IEE Proceedings H (Microwaves, Optics and Antennas) , volume 130, pages 175-182. IET, 1983.
- [19] Juan Mosig. The weighted averages algorithm revisited. IEEE Transactions on Antennas and Propagation , 60(4):20112018, Apr. 2012.
- [20] M. J. Nam, D. Pardo, C. Torres-Verd´ ın, S. Hwang, K. G. Park, and C. Lee. Simulation of eccentricity effects on shortand long-normal logging measurements using a Fourier-hp-finite-element method. Exploration Geophysics , 41(1):118-127, 2010.
- [21] Marcela S. Novo, Luiz C. da Silva, and F. L. Teixeira. Finite volume modeling of borehole electromagnetic logging in

- 3-D anisotropic formations using coupled scalar-vector potentials. IEEE Antennas and Wireless Propagation Letters , 6: 549-552, 2007.
- [22] Marcela S. Novo, Luiz C. da Silva, and Fernando L. Teixeira. Three-dimensional finite-volume analysis of directional resistivity logging sensors. IEEE Transactions on Geoscience and Remote Sensing , 48(3):1151-1158, Mar. 2010.
- [23] D. Pardo, L. Demkowicz, C. Torres-Verd´ ın, and M. Paszynski. Two-dimensional high-accuracy simulation of resistivity logging-while-drilling (LWD) measurements using a self-adaptive goal-oriented hp finite element method. SIAM Journal on Applied Mathematics , 66(6):2085-2106, Oct. 2006.
- [24] D. Pardo, C. Torres-Verd´ ın, and L. F. Demkowicz. Simulation of multifrequency borehole resistivity measurements through metal casing using a goal-oriented hp finite-element method. IEEE Transactions on Geoscience and Remote Sensing , 44 (8):2125-2134, Aug. 2006.
- [25] Zhengyong Ren and Jingtian Tang. 3D direct current resistivity modeling with unstructured mesh by adaptive finiteelement method. Geophysics , 75(1):H7-H17, Jan.-Feb. 2010.
- [26] Yutaka Sasaki. 3-D resistivity inversion using the finite-element method. Geophysics , 59(12):1839-1848, Dec. 1994.
- [27] William Ralph Smythe. Static and Dynamic Electricity . International Series in Pure and Applied Physics. McGraw-Hill, New York, second edition, 1950.
- [28] W. Squire. Partition-extrapolation methods for numerical quadratures. International Journal of Computer Mathematics , 5(1):81-91, 1975.
- [29] W. M. Telford, L. P. Geldart, and Robert E. Sheriff. Applied Geophysics . Cambridge University Press, Cambridge England and New York, second edition, 1990.
- [30] J. R. Wait and K. R. Umashankar. Analysis of the earth resistivity response of buried cables. Pure and Applied Geophysics , 117(4):711-742, 1978.
- [31] James R. Wait. Geo-Electromagnetism . Academic Press, New York, 1982.
- [32] G. N. Watson. A Treatise on the Theory of Bessel Functions . Cambridge Mathematical Library. Cambridge University Press, New York, second edition, 1995.
- [33] Ernst Joachim Weniger. Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series. Computer Physics Reports , 10:189-371, 1989.
- [34] Shanjie Zhang and Jian-Ming Jin. Computation of Special Functions . John Wiley, New York, 1996.
- [35] Z. Zhang and Z. Zhou. Real-time quasi-2-D inversion of array resistivity logging data using neural network. Geophysics , 67(2):517-524, Mar.-Apr. 2002.
Haksu Moon Fernando L. Teixeira Burkay Donderici
Computation of potentials from current electrodes in cylindrically stratified media: A stable, rescaled semi-analytical formulation
We present an efficient and robust semi-analytical formulation to compute the
electric potential due to arbitrary-located point electrodes in
three-dimensional cylindrically stratified media, where the radial thickness
and the medium resistivity of each cylindrical layer can vary by many orders of
magnitude. A basic roadblock for robust potential computations in such
scenarios is the poor scaling of modified-Bessel functions used for computation
of the semi-analytical solution, for extreme arguments and/or orders. To
accommodate this, we construct a set of rescaled versions of modified-Bessel
functions, which avoids underflows and overflows in finite precision
arithmetic, and minimizes round-off errors. In addition, several extrapolation
methods are applied and compared to expedite the numerical evaluation of the
(otherwise slowly convergent) associated Sommerfeld-type integrals. The
proposed algorithm is verified in a number of scenarios relevant to geophysical
exploration, but the general formulation presented is also applicable to other
problems governed by Poisson equation such as Newtonian gravity, heat flow, and
potential flow in fluid mechanics, involving cylindrically stratified
environments.