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.
keywords:Poisson equation, steady-state diffusion equation, discontinuous coefficients, stratified media, resistivity logging, electric potential
Resistivity logging is extensively used for detecting, characterizing, and analyzing hydrocarbon-bearing zones in the subsurface earth Wait and Umashankar (1978); Moran and Gianzero (1979); Gianzero and Anderson (1982); Drahos (1984); Lovell and Chew (1990); Lovell (1993); Zhang and Zhou (2002); Doetsch et al. (2010). This sensing modality employs electrode-type devices mounted on a mandrel that inject electric currents into the surrounding earth formation Ellis and Singer (2007); Telford et al. (1990). 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 Liu and Sato (2002); Hue et al. (2005); Lee et al. (2012); Sasaki (1994); Pardo et al. (2006a, b); Nam et al. (2010); Ren and Tang (2010); Novo et al. (2007, 2010); Liu et al. (1994); Fan et al. (2000); Hue and Teixeira (2007). 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 Chew (1995). 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 arbitrary-located 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 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 Smythe (1950); Wait (1982); Carley (2013). 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 Moon et al. (2014) . 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.1 Electric potential in homogeneous media
In a homogeneous medium, the electric potential from a current electrode at the origin writes as Wait (1982)
where is the electric current flowing into the medium from the electrode, is the conductivity of the medium, and is the modified-Bessel function of the second kind of the zeroth order. For the second equality, the complete Lipschitz-Hankel integral Watson (1995) is employed. When the source is off the origin, higher order azimuthal modes appear. Using the addition theorem for , (1) is modified to
in terms of modified-Bessel functions of both first, and second, , kinds. In the above, primed coordinates () represent the source location and unprimed coordinates () represent the observation point. Also, and .
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 is a function of because 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 is a function of instead of because diverges when goes to zero. For the outgoing-potential case, the -th harmonic with dependence in layer 1 and layer 2 can be expressed, resp., as
where and are the (local) reflection and transmission coefficients at the boundary , and is an arbitrary amplitude of the primary potential. Applying the boundary conditions Wait (1982) at the interface, we obtain
For the standing-potential case, we similarly have
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 Chew (1995) 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
Applying proper constraint equations, we obtain
If another layer is added beyond layer 3, only needs to be replaced by . Therefore, the generalized reflection coefficient in cylindrically stratified media for the outgoing-potential case is
All amplitudes ’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
From (10), a new coefficient denoted by is defined as
such that . The above coefficient can be regarded as a ‘local’ transmission coefficient between two adjacent layers, as depicted in Fig. b. Generalized transmission coefficient described in Fig. 3 can be defined using the -coefficients (11) through
Note that for the outgoing-potential case. When , in (12).
For the standing-potential case depicted in Fig. 4, the potentials in three cylindrical layers can be written as
Similarly, applying proper constraint conditions, we obtain
and, more generally,
To obtain all amplitudes ’s, it is convenient to define the coefficient with decreasing subscripts . Since the relation between two successive amplitudes is
the coefficient in this case is defined as
such that . The generalized transmission coefficient depicted in Fig. 5 can then be written as
Note that for the outgoing-potential case. When , in (18).
where the two unknowns are and . Applying two constraint conditions at and , we obtain
where . When , the potential in layer is expressed as
where is the amplitude of the outgoing-potential in layer and expressed as
Therefore, the integrand factor for layer in (22) is
When , the potential in layer is expressed as
where is the amplitude of the standing-potential in layer and expressed as
The integrand of the integral for potential in layer in (25) is
In summary, the potential in cylindrically stratified media admits four different expressions depending on the relative position of and , which can be written as
|Case 1: and are in the same layer and|
|Case 2: and are in the same layer and|
|Case 3: and are in different layers and|
|Case 4: and are in different layers and|
2.3 Rescaled modified-Bessel functions
The electric potential involves products of the modified-Bessel function of the first and second kind, viz. and . Those products can involve disparate values due to the exponential behavior of the functions. For example, when , has very large value whereas has very small value. This disparity becomes progressively worse for higher order modes. On the other hand, when , has very small value while 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 Moon et al. (2014) for Bessel and Hankel functions, viz. and . In order to apply such rescaled functions, the analytical expressions for the potential need to modified accordingly, as described next.
When , and can be expressed via small argument approximations for . Noting that, for , the relationship between cylindrical functions and modified-cylindrical functions reads as
Thus, we have the following small argument approximations for and and their derivatives.
It should be noted that the multiplicative factor above is chosen so as to not depend on the radial distance, . Also, is identical for a function and its derivative, and the multiplicative factors appearing in and are reciprocal to each other. This will facilitate some computations later on.
When and , the large argument approximations for the modified-Bessel functions write as Zhang and Jin (1996)
where . 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
If the argument is neither small nor large, rescaled modified-cylindrical functions are defined, in analogy to small and large arguments, as
where the multiplicative factor is defined as in Moon et al. (2014), i.e.,
with used for double-precision arithmetic computations.
and its subscript is linked to . 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.
|small arguments||moderate arguments||large arguments|
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 Moon et al. (2014) and not repeated here
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 is associated with , whereas is associated with . Note again that the subscript refers to to the index of the radial factor in the argument. There are two important properties to note for and : ()Reciprocity: and () Boundness: . As we will see below, these two properties are important in ensuring stable numerical computations.
Recalling the expressions for the reflection and transmission coefficients obtained before, the reflection coefficient for the outgoing-potential case is modified to
Similarly, it can be shown that the reflection coefficient for the standing-potential case is modified to , and that the transmission coefficient for the outgoing-potential case and the transmission coefficient for the standing-potential case simply recover the original ones without any multiplicative factors, i.e. and .
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 and , and that, for both the outgoing-potential and standing-potential cases, .
In addition to generalized reflection and transmission coefficients, the factors and considered before are also required to compute the potential. The basic difference between the two types of coefficients is that involves two generalized reflection coefficients whereas involves only one generalized reflection coefficient. All these auxiliary coefficients can be redefined accordingly using rescaled reflection coefficients, i.e.,
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: , , , and . For convenience, we let , , , and so that , and the integrand rewrites as
where has been used. All multiplicative factors , , , have magnitudes no larger than one due to the boundness property discussed above.
Similarly, for Case 2, there are four radial parameters of interest: , , , and . For convenience, we let