Three-dimensional inversion of inductions logs in real-time

(1)

ree-dimensional inversion of

induction logs in real-time



ter verkrijging van de graad van doctor

aan de Technische Universiteit Del,

op gezag van de Rector Magniﬁcus prof.ir. K.C.A.M. Luyben,

voorzier van het College voor Promoties,

in het openbaar te verdedigen

op donderdag 19 maart 2015 om 10.00 uur

door

Silvian BENSDORP

elektrotechnisch ingenieur,

geboren te Amsterdam

(2)

Dit proefschri is goedgekeurd door de promotoren: Prof.dr.ir. P.M. van den Berg

Prof.dr.ir. J.T. Fokkema

Samenstelling promotiecommissie:

Rector Magniﬁcus Voorzier

Prof.dr.ir. P.M. van den Berg Technische Universiteit Del, promotor

Prof.dr.ir. J.T. Fokkema Technische Universiteit Del, promotor

Prof.dr. S.A. Petersen Technische Universiteit Del, Statoil ASA

Prof.dr.ir. E.C. Slob Technische Universiteit Del

Prof.dr.ir. J. van der Kruk Forschungszentrum Jülich

Dr. P.A. Olsen Statoil ASA

Dr.ir. R.F. Remis Technische Universiteit Del

is research was ﬁnancially supported by Statoil ASA, Norway. ISBN: 978–94–6259-587–3

All rights reserved. No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or mechanical, including photocopying, recording or by any information storage and retrieval system, without prior permission of the author.

(3)

H

ydrocarbons, without any doubt, play a fundamental role in today’s society,

ei-ther as an energy source or as a building block for the chemical industry to produce the materials that are paramount to our health and well-being. Continuous advancements in geophysical prospecting made it possible to unlock resources in ever more compli-cated hydrocarbon reservoirs in an economical viable way. As current reservoirs are being depleted, there is a constant need to start operations on new reservoirs, or to improve operations in current reservoirs. New operations can be employed as soon as emerging prospecting technologies become commercially tractable, to explore and extract reserves in complicated geological seings.

Traditionally, geophysical prospecting has a strong focus on acoustic waves as the principal method for obtaining information about the subsurface; see Fokkema and van den Berg (1993), Gisolf and Verschuur (2010). e propagation characteristics of acoustic waves depend on the medium parameters, which inﬂuence how energy be-tween diﬀerent physical quantities is continuously exchanged. Consequently, the only medium parameters we can obtain information from, are the ones that partake in the energy exchange of the physical quantities we either excite or measure.

While elastodynamic medium parameters are very useful to map out the structural geology, they are ill-suited to provide information on the hydrocarbon content of a structure. e delineation of hydrocarbon reservoirs from seismic maps is done by geological interpreters. For porous subsurface rock, the content consists of water, nat-ural gas or petroleum. In an uncontaminated state, both water and almost all stable organic compounds have no electric conductivity. In nature, though, trace minerals will dissolve in the water, making the water a strong electric conductor. us the elec-tric conductivity serves as a proxy for the location hydrocarbon resources. A physical quantity naturally susceptible to changes in the electric conductivity is the electromag-netic ﬁeld (Jackson 1999). On the spatial scale of the seismic experiment, this medium parameter is probed with the controlled-source electromagnetic () method (see Gribenko and Zhdanov (2007) or overviews by Constable (2010) and Zhdanov (2010)). Although this technique gives a good indication on the location of reservoirs, the ob-tained medium descriptions have a very low spatial resolution that makes the obob-tained images ill-suited for operational planning in the hydrocarbon production process.

(6)

. . physical reality . mathematical relations . measure-ments . synthetic data . recon-struction . scaering object . modelling error . inverse modelling . forward modelling . reconstruction error

Figure 1.1: Conceptual relation between synthetic data and measurements on the le side, and the scaering object and its representation obtained through inverse modelling on the right side.

1.1 A      

More general, to obtain a description of some physical structure (in the electromag-netic case, a scaering object consisting of a distribution of the electric conductivity) we need to probe it with an electromagnetic source and record a resulting perturbation, the data, of the electromagnetic field as well. To understand how the scaering object relates to the measured data we need a forward model. e model serves as a mathe-matical description how the measured field is related to the source and the scaering object. To understand how the data relates to a scaering object for a known source we need an inverse model. e relation between these two models is illustrated in figure 1.1. Of major importance is that the characteristics of the two modelling relations are fundamentally different, as the data is in fact always incomplete.

In a forward problem, we create a virtual representation of the reality, i.e. the scaer-ing object and a source to excite an electromagnetic ﬁeld. rough our forward model this results in synthetic data. If the representation is accurate, the diﬀerence between the synthetic data and the measurement data, if we would perform the experiment in the physical world, should be small. Hence, if the modelling error falls below an acceptable error threshold we are able to simulate the physical reality.

For an inverse problem (see e.g. Colton and Kress (1992), Kirsch (1996), Vogel (2002)), the scaering object is an unknown quantity, whereas the measurements and the source characteristics are known. e mathematical representation of the inverse model is in general ill-posed. To solve the inverse problem, we create a reconstruction based on an approximate inverse operator. is reconstruction is then used to solve the forward problem. e modelling error then serves as a measure how far the reconstruction is

(7)

1.2. Induction logging 3

apart from the physical reality. Based on this error, the reconstruction is adjusted, and we restart the loop until the modelling error is reduced to a satisfactory level. While clearly there is a relation between the accuracy of the forward model and the accu-racy of the modelling error, it does not guarantee to obtain a reliable reconstruction. Furthermore, accurate forward models can be very cumbersome to process, and it is oen more beneﬁcial to iterate the loop many times with a simple approximate model rather than to perform a few iterations with a more elaborated forward model to obtain a useful reconstruction.

1.2 I 

With exploration and production drilling operations becoming more complex, there is an increased need for information to guide the drilling. Such information will typically contain rock properties such as porosity and water saturation. Predicting the distribu-tion of the properties around a drilling bit can assist in real-time operadistribu-tions planning (see ﬁgure 1.2(a)). A common method to obtain information about these properties is

low-frequency1_{electromagnetic logging with multiple sources and receivers. With the}

introduction of high-capacity telemetry systems for borehole systems, measurements are available in real-time for processing and interpretation.

Electromagnetic logging in geophysical exploration, commonly called induction log-ging, is a relevant method to discriminate between hydrocarbon-bearing and water (or shale)-bearing zones in the subsurface. eoretical principles of the induction-logging method for some relatively simple canonical configurations can be found in Kaufman and Dashevsky (2003). Many examples from a commercial-industrial viewpoint can be found in the thesis of Anderson (2001). e physical principle underlying the method is the probing of the differences in electrical conductivity between the different zones. When an induction logging instrument is lowered into a borehole, the electromagnetic field of the magnetic-dipole sources in the tool induces electrical currents in the sub-surface formation. ese induced currents contribute to the measured response in the magnetic-dipole receivers which are also located in the tool some distance apart of the magnetic-dipole sources. e interpretation of the measured response in terms of the formation conductivity gives in principle an indication for the location of the hydrocar-bon bearing zones.

1.2.1 Forward modelling

In this thesis we develop a method approximating the system response of an induction logging tool for the purpose of generating logging responses for realistic earth conduc-tivity conﬁgurations. e method aims to approximate in a reliable and computational fast way the response of a logging tool along an arbitrarily prescribed borehole tra-jectory in a full three-dimensional (3) earth model, such that diﬀerent realizations of

1_{With low-frequency we specifically refer to electromagnetic fields showing dominant diffusive}

(8)

... borehole . earth depth .. . .. . .

earth surface or seabed

(a) Global illustration of borehole measurement systems.

... . . source coil . receiver coil . drilling tool . borehole

(b) Illustration the configuration of a induction logging system inside a borehole.

Figure 1.2: Schematic overview of the logging configuration in a geological seing (a). e field amplitude contour lines are based on a time-harmonic magnetic field originating from a magnetic dipole source polarized in the borehole direction (b).

both borehole trajectories and earth models can be evaluated eﬀectively. From a physi-cal point of view, current logging tools consist of a number of magnetic-dipole sources (implemented as source coil antennas) located at the tool axis in the direction of the tool axis and a number of magnetic-dipole receivers (implemented as receiver coil antennas) located at the tool axis in an arbitrary orientation, as illustrated in ﬁgure 1.2(b).

e computation of the response of a logging tool in a 3 inhomogeneous medium requires a code in full 3, based on the Maxwell equations. ese codes can be clas-siﬁed according to their numerical implementation. In induction logging three types are commonly employed: contrast-type of integral-equation methods (amongst

(9)

oth-1.2. Induction logging 5

ers Zwamborn and van den Berg (1992), Avdeev et al. (2002), Zhang and Liu (2003), Avdeev and Knizhnik (2009), Nie et al. (2013)); numerical representations based on finite-element methods, e.g. Badea et al. (2001), Nam et al. (2013); and finite-difference methods (e.g. Newman and Alumbaugh (2002), Weiss and Newman (2002), Davydy-cheva et al. (2003 2009), Wang et al. (2009) amongst others). Other approaches such based on finite-dimensional integration are described by Bossavit and Mayergoyz (1998) and Hiptmair (2002), yet to catch the aention of the geophysical community. While these forward modelling programs are nowadays available or becoming available, the computational burden is too large to carry out computations for different realizations of borehole trajectory and realistic earth models in a time-efficient maer. Hence, if we are to rely on an effective approximate model it must include all the necessary physics. Direct linearization of the problem, assuming that the actual electric field in the domain of observation is equal to the background field, thus neglecting the scaered field, the so-called Born approximation, is too crude. Several variations to the Born approxima-tion have been proposed to alleviate the crudeness of this approximaapproxima-tion while trying to retain some of the computational simplicity of the approximation (e.g. Zhdanov et al. (2000), Torres-Verdín and Habashy (2001), Tseng et al. (2003), Abubakar and Habashy (2005), Guozhong and Torres-Verdín (2006)).

e logging-curve response encompasses measurements of an induced magnetic field along a trajectory through a 3 inhomogeneous medium for one or more pre-scribed frequencies of operation. e method we develop allows for the definition of an arbitrarily curved logging trajectory dictated by the borehole position, along which the electromagnetic response is computed. In our analysis, the borehole trajectory is replaced by locally straight line segments. Along each line segment, the electromag-netic field strength is only significant within a 3 volumetric window of limited spatial dimensions. During the computation the window can move and rotate as it follows the trajectory. e size of this window of observation depends both on the frequency of operation and the local electrical conductivity of the earth formation around the tool. In principle only the earth formation has control over the induced electromagnetic field: for simplicity we choose not to compensate for secondary borehole effects such as the effect of steel casing of the borehole (Kim and Lee 2006), or induction tool eccentricity (Lovell and Chew 1990, Wang et al. 2003, Tompkins et al. 2004).

To numerically describe the electromagnetic field in this thesis we use the contrast-type of integral equation as mathematical tool to formulate the logging problem at hand. However, in each reduced window, a homogeneous background medium is cho-sen, where the electromagnetic field in this background medium is denoted as the pri-mary electromagnetic field. In each local window, this pripri-mary electromagnetic field is obtained directly from a simple closed-form expression. Within each window, a back-ground medium is chosen in such a way that the variety of the actual conductivity with respect to the value of the chosen background is as small as possible. One way to obtain such conductivity background is to average the either synthesized or a-priori defined conductivity around the tool domain. e background conductivity is determined from

(10)

the measured data by assuming that the measured ﬁelds are approximated by the pri-mary ﬁelds only.

Subsequently, the electric currents due to the differences in the electrical conductiv-ity with respect to the effective conductivconductiv-ity of the background medium in the window under investigation, are seen as contrast currents that generate a scaered (secondary) field. In principle, this approach leads to contrast-type of integral-equation formulation (see e.g. Ward and Hohmann (1988) or van Bladel (2007)). e computational solution of the integral equation over the reduced 3 window requires too much computational time for real-time processing of the well logs. However, in view of the reduced size of each local window and the relative small changes of the contrast in electrical con-ductivity with respect to the concon-ductivity of the matched homogeneous background, the interaction between different regions within the local window can be neglected and each contrasting region may be seen as a single spherical scaerer (Slob 1994). is so-called Single Spherical Scaerer (SSS) approximation is used advantageously in this thesis to provide a simple and effective forward model for the true disturbance of the electromagnetic field resulting from the contrasting conductivity in the reduced win-dow2_.

1.2.2 Inverse modelling

e inversion of induction logging data in the exploration industry has been performed for over 35 years. Some early contributions include Oristaglio and Worthington (1980), Freedman and Minerbo (1991), Chew and Liu (1994), van den Berg and van der Horst (1995), Tabarovsky and Rabinovich (1998). All these approaches rely on a 1 conductiv-ity model. Contemporary approaches focus on 2.5 medium reconstructions (Abubakar et al. 2006b) and parametrized models using conductivity type reconstructions (Abuba-kar et al. 2006a). Lately, using finite-difference or volume integral equation methods, direct petrophysical parameter inversion for parameters such as fluid permeability and porosity (Alpak et al. 2006 2008) are employed. Although finite-element simulations in 2 domains are becoming available (Nam et al. 2009, Chen 2010, El Bouajaji et al. 2011), such methods require continuous grid adoption to function effectively (Böhm et al. 2000). eir 3 counterparts are currently too computationally burdensome to be used in practice in inverse reconstruction algorithms. To overcome the 3 complexity, in Haber et al. (2004) the authors construct a linear system in which they both include the solution constraints and the forward problem, solving both demands simultane-ously at the cost of an increased system matrix. Other approaches that extract a more structural medium description include Mendoza et al. (2012), where bed layers and dip angles are estimated.

Using a different source-receiver geometry, another common data-acquiring con-figuration is cross-well induction logging, where sources and receivers are in different

(11)

1.3. esis outline 7

wells with the scaering object between them — see e.g. Daily and Owen (1991), Abu-bakar and van den Berg (1998) and AbuAbu-bakar et al. (2005).

An approach using constraining equations for reconstruction images algorithms can be found in Afonso et al. (2011). Reconstruction for logging problems where both source and receiver are located on the same axis are in general more diﬃcult to resolve than for cross-well problems, where the receiver placement is much more favourable with respect to the source placement (Alumbaugh and Morrison 1995, Li et al. 2010). All these deterministic approaches have in common that they are based on the minimization of some cost functional. An entirely diﬀerent approach to these kind of problems is a statistical approach. Some nice introductions into this topic can be found in Evans and Stark (2002) and Tenorio (2001).

In our research we focus on a novel way of reconstructing the 3 electric conduc-tivity distribution around a drilling tool. To obtain such distribution, the reconstruction problem must be cast into a mathematical formulation for which we can compute a met-ric such that we can quantify amongst different reconstructions which one is the best. e process to find the smallest/largest metric is called optimization. A good introduc-tion on the wealth of this topic for convex metrics can be found in Kelley (1999) or Boyd and Vandenberghe (2004). Typically this metric takes into account the modelling error (cf. figure 1.1) in combination with extra conditions imposed on the structure of the reconstructed scaering object. We use a formulation based on the contrast-source inversion algorithm by van den Berg and Kleinman (1997) to define and solve our op-timization problem. is formulation is extended to force a solution which minimizes the relative total variation3_{at each iteration (van den Berg et al. 2003).}

1.3 T 

In this thesis we develop a fast methodology to reconstruct the geology around a bore-hole using synthetic induction logging data generated by a generic tool conﬁguration. e physics of the induction logging system is described by electromagnetic wave the-ory. is theory is ubiquitous throughout this thesis.

In apter 2, we derive a forward model based on the integral-equation representa-tion of the Maxwell equarepresenta-tions. e integral equarepresenta-tion is derived specifically for diffusive wave-fields, where the displacement currents are ignored. e derivation of the numer-ical representation is achieved using the finite-dimensional representation on a regular grid, giving rise to a convolution-type matrix structure that can be evaluated rapidly using the fast Fourier transform. Because we are primarily interested in the magnetic fields at specific points in space, a finite-dimensional evaluation of the magnetic field pertained to the electric field is discussed as well.

In apter 3, we derive the Single Spherical Scaerer approximation, an approxima-tion to the diffusive electromagnetic field that is free of any differential or integral-type

(12)

structure. For different source- and receiver configurations, we derive expressions for the magnetic fields at the receiver positions based on the Single Spherical Scaerer (SSS) approximation. Furthermore, we derive an estimate for the penetration depth of the tool based on the operating frequency and the medium surrounding the induction logging instrument. In a statistical analysis, we compute figures of merit for different induction logging systems parameters using a geological model. In this analysis we directly compare the SSS approximation with the well-known Born approximation.

In apter 4, we describe a reconstruction procedure for the electric conductivity around a borehole, based on the contrast-source inversion algorithm. Typically, com-mercial induction logging systems have only few receivers. To obtain enough data to make a meaningful interpretation of the geology surrounding the borehole, we use many different tool positions, up to a point were there is no overlap in the area that influence the pertaining measured data. We modify the contrast-source inversion algo-rithm to accommodate for the different probing areas in a single unified reconstruction procedure. Moreover, a regularization procedure based on total variation minimization is described to enhance the resolution of the reconstructed image.

As the electric conductivity can in principle only be recovered from the scaered electric or magnetic field, in apter 5, we describe how we can separate the incident magnetic field component from a measured magnetic field. is procedure yields a background electric conductivity for a local window that is specific for a single tool location and operating frequency. Also, we derive a simple decision algorithm for the effective gradient direction of a medium surrounding a tool at a specific position.

To show the results of our reconstruction procedure, in apter 6, we compute syn-thetic induction logs for numerous artiﬁcial and realistic media. We show the per-formance of our procedure where both the data generation and the inverse model are based on the SSS approximation, as well as the performance for data generated using the accurate model described in chapter 2. We show convergence plots of the complete reconstruction, as well as the performance of the individual components partaking in the reconstruction. is serves as an analytical aid to determine the reliability of the reconstruction in the diﬀerent parts of the reconstructed object.

Last, in apter 7, we show some practical improvements on the basic ideas as out-lined in chapter 4. We perform multiple reconstructions based on a synthetic well-log stretching over 180 meters. To enhance the image quality, we separate the well-logs in a low spatial and a high spatial frequency part (k-domain ﬁltering), and use these new logs in a joint reconstruction procedure. To optimize the reconstruction procedure for deployment in a real-time environment, we modify our reconstruction system such that new data is included as it becomes available, while we keep parts of the reconstruction that have already been performed. In this way, instead of starting a new reconstruction procedure from scratch each time, we can produce a reconstruction using new data at only a fraction of the processing time.

(13)

1.3. esis outline 9

our ﬁnal thoughts, and oﬀer discussion material for future development perspectives in induction logging inversion.

(14)

(15)

Chapter Two

e electromagnetic ﬁeld equations and their

numerical solution

Introduction

Induction logging is a technique of probing the surroundings of a source with an elec-tromagnetic field, and measuring the resulting magnetic field at different locations and polarizations along the measurement device. Due to differences in the medium prop-erties the electromagnetic field radiated by the source scaers continuously over the entire domain. In case of induction logging, the medium properties are the electric permiivity, the electric conductivity and the magnetic permeability. We assume a negligible electric permiivity and constant magnetic permeability. In this case, the magnetic field in any point of space follows directly from the electric vector field.

In this chapter, we first derive the electromagnetic-field integral equations for both the magnetic and the electric fields, Secondly, we derive a finite-dimensional approxi-mation for the electromagnetic field, which allows us to compute the electromagnetic fields over arbitrary medium distributions.

2.1 F  

Maxwell’s equations

All classical (i.e. non-quantum, non-relativistic) electrodynamic effects in nature can be described by the set of Maxwell equations. e electric-field strength vector E and the magnetic-field strength vector H satisfy a pair of vector differential equations with spatial coordinates x = {x1, x2, x3} ∈ R3, expressed using Cartesian unit vectors (i1, i2, i3)and temporal coordinate t∈ R+of form

−∇ × H(x, t) + ∂tD(x, t) + Jcon(x, t) =−Jext(x, t), (2.1)

∇ × E(x, t) + ∂tB(x, t) =−Kext(x, t), (2.2)

where vector ﬁelds D(x, t) and B(x, t) represent the electric and magnetic ﬂux density, respectively, Jext_{(x, t)}_{denotes the external electric current density and J}con_{is the}

electric current density from conduction eﬀects. For mathematical convenience, the magnetic current source density Kext_{(x, t)}_{is introduced. e operator ∂}

(16)

temporal derivative, the operator∇ = {∂1, ∂2, ∂3} denotes the spatial diﬀerentiation with respect to x1, x2and x3.

e tangential components of E(x, t) and H(x, t) are continuous throughout the entire space R3_{. is is represented by the boundary conditions}

n(x)× H(x, t)|_Ω1 = n(x)× H(x, t)|_Ω2, (2.3)

n(x)× E(x, t)|_Ω

1 = n(x)× E(x, t)|Ω2, (2.4)

where n(x) is the normal vector along any source-free interface between two spatial domains Ω1, Ω2, both in R3_.

Compatibility relations

We take the divergence∇· on both sides of the Maxwell equations (2.1) and (2.2)) to

obtain the compatibility relations

∂t∇ · D(x, t) + ∇ · Jcon(x, t) =−∇ · Jext(x, t), (2.5)

∂t∇ · B(x, t) = −∇ · Kext(x, t). (2.6)

On a source-free interface between Ω1and Ω2, the continuity relations

n(x)· [∂tD(x, t) + Jcon(x, t)] = 0, (2.7)

n(x)· ∂tB(x, t) = 0, (2.8)

hold.

e temporal Fourier transform

e Fourier transform of a function f(t) on R+_{is deﬁned as}

˜ f (ω) =

∫ R+

f (t) exp(iωt) dt, (2.9)

where i2₌_{−1 is the imaginary unit. For the transformation of the temporal dimension,} the space of all ω is commonly known as the frequency domain. A particular benefit of the Fourier transform is that it works on the dual space of functions, the so-called generalized functions, which includes functions that are not continuously differentiable. With the temporal Fourier transform defined, we set the restriction that we will only use monochromatic sources, and assume that the relation between the field quantities is linear. is implies that the resulting electromagnetic field is harmonic with respect to time.

To transform the time-domain electromagnetic equations to the frequency domain, we transform the original Maxwell equations (2.1–2.2) with the Fourier transform, we

(17)

2.1. Field integral equations 13

map the set of square integrable functions on R3_{× R}+ _{to a set of functions in the}

complex space C3_{, using the transformation signature exp(}_{−iωt). is results in} −∇ × ˜H(x, ω)− iω ˜D(x, ω) + ˜Jcon(x, ω) =− ˜Jext(x, ω), (2.10)

∇ × ˜E(x, ω)− iω ˜B(x, ω) =− ˜Kext(x, ω), (2.11) where the tilde indicates a dependency on frequency rather than time. e substituted continuity relations are transformed as

−iω ∇ · ˜D(x, ω) +∇ · ˜Jcon(x, ω) =−∇ · ˜Jext(x, ω), (2.12) −iω ∇ · ˜B(x, ω) =−∇ · ˜Kext(x, ω). (2.13) In all subsequent treatment, we will only work with the representation in the frequency domain, and hence the tilde symbol will be omied for brevity.

Constitutive relations

While the microscopic interaction between fields and maer can be very complex, we are usually only interested in an approximation of the collective behaviour of such ap-proximations on a scale proportional to the temporal change of the field; these are called the macroscopic approximations. Such approximations will greatly simplify the descrip-tion of electromagnetic interacdescrip-tion, while for most applicadescrip-tions, including inducdescrip-tion logging, providing sufficient representation of the underlying physics.

For the electric ﬁeld the two useful relations are Ohm’s law for the conduction cur-rents in the le-hand side of (2.10),

Jcon(x, ω) = σ(x, ω)E(x, ω), (2.14)

and the electric-ﬂux density

D(x, ω) = ε(x, ω)E(x, ω), (2.15)

where σ(x, ω) is the electric conductivity and ε(x, ω) is the electric permiivity of the medium under consideration. For induction logging configurations (i.e. low frequency wave fields) it is assumed that the electric-flux density currents play a negligible role.

is implies we have ωε(x, ω)≪ σ(x, ω), and therefore we obtain the simpliﬁcation

iωD(x, ω) + Jcon(x, ω) = σ(x, ω)E(x, ω) (2.16)

For non-magnetizable materials we may apply the simple constitutive relation for the magnetic ﬂux

(18)

where µ0is the frequency independent and constant permeability in R3_{. Using the} constitutive relations and grouping together all current sources we obtain the diﬀusive Maxwell equations

−∇ × H(x, ω) + σ(x, ω)E(x, ω) = −Jext_{(x, ω),} _(2.18)

∇ × E(x, ω) − iωµ0H(x, ω) =−Kext(x, ω). (2.19)

ese will be the electromagnetic equations from which we will derive integral repre-sentations for further analysis. ese are linear equations, meaning that any superpo-sition of electromagnetic ﬁelds will result in a valid electromagnetic ﬁeld as well.

Field integral representations in a homogeneous medium

In a homogeneous space the conductivity becomes independent of the spatial location x. We deﬁne a background conductivity σb such that σ(x, ω) = σb(ω). en, by

substituting (2.18) in (2.19) and vice versa to eliminate either E(x, ω) or H(x, ω), we obtain

−∇ × ∇ × E(x, ω) + iωσb(ω)µ0E(x, ω) =

−iωµ0Jext(x, ω) +∇ × Kext(x, ω), (2.20)

−∇ × ∇ × H(x, ω) + iωσb(ω)µ0H(x, ω) =

σb(ω)Kext(x, ω)− ∇ × Jext(x, ω), (2.21)

respectively. With the compatibility relations (2.12–2.13) and some rearrangements we arrive at the vector Helmholtz equations

(

∇2_{+ k}2 b )

E(x, ω) =

−iωµ0Jext(x, ω)− (σb)−1∇∇ · Jext(x, ω) +∇ × Kext(x, ω), (2.22) (

∇2 + kb2

)

H(x, ω) =

σbKext(x, ω) + (iωµ0)−1∇∇ · Kext(x, ω)− ∇ × Jext(x, ω), (2.23) where we use the background wave number

k2_b= iωµ0σb. (2.24)

e explicit solution of these equations is obtained via the Green’s function method, viz.

E(x, ω) = iωµ0A(x, ω) + σ_b−1∇∇ · A(x, ω) − ∇ × F (x, ω), (2.25) H(x, ω) =−σbF (x, ω)− (iωµ0)−1∇∇ · F (x, ω) + ∇ × A(x, ω). (2.26)

(19)

with vector potentials A(x, ω), F (x, ω) deﬁned as

A(x, ω) = ∫ Ω g(x− x′, ω)Jext(x, ω) dV, (2.27) F (x, ω) = ∫ Ω g(x− x′, ω)Kext(x, ω) dV, (2.28)

over our domain of interest Ω⊂ R3_{that contains the external sources, as is illustrated} in ﬁgure 2.1. e Green’s function for the 3 Helmholtz vector equation on a Cartesian coordinate system, in a homogeneous background, is the solution of

(

∇2_{+ k}2 b )

g(x− x′, ω) =−δ(x − x′), (2.29)

where δ(·) is the Dirac delta distribution. e solution is obtained using the spatial Fourier transform technique as

g(x− x′, ω) = exp(ikb∥x − x

′_∥)

4π∥x − x′∥ , (2.30)

where we take the branch of kbsuch that the ﬁeld decays for increasing∥x − x′_{∥, i.e.}

kb= 1 + i √ 2 √ µ0σbω. (2.31)

e scattering problem

e fields defined in a homogeneous medium are called the incident fields. Due to changes in the medium parameters (in our case the conductivity σ), the total fields E(x, ω)and H(x, ω) will deviate from the incident field. As any linear combination of vector fields that are a solution to the Maxwell equations, there must exist a solution, called the scaered field, such that

Esct(x, ω) = E(x, ω)− Einc(x, ω),

Hsct(x, ω) = H(x, ω)− Hinc(x, ω). (2.32)

e incident ﬁeld

e incident ﬁeld is the electromagnetic ﬁeld that is generated by all external current sources, both electric and magnetic, in a homogeneous background. We only consider

magnetic-dipole sources. Assuming the source has an eﬀective current density I0

po-larized in the M direction, we can denote the external magnetic dipole located at xS

as

(20)

.. xS .. xR .. borehole . R3

(a) Incident field problem configuration on an arbitrary borehole.

..

Ω .

σ(x, ω)_. ̸= σb

σ(x, ω) = σb

(b) Generic setup for the scaering problem.

Figure 2.1: Schematic overview of the incident field problem and the scaering prob-lem. e solution to the electromagnetic field problem is the superposition of both fields.

where iωµ0M is the dipole moment. Computing the vector potential F (x, ω), the

expressions for the incident ﬁelds in the points x̸= xS_are

Einc(x, ω) =−iωµ0I0∇ × g(x − xS, ω)M , (2.34)

Hinc(x, ω) =−I0 [

k_b2g(x− xS, ω)M − ∇∇ · g(x − xS, ω)M]. (2.35)

e scattered ﬁeld

e scaered fields are the electromagnetic fields generated by the electric conduc-tion currents. As we will work with non-magnetizable maer only, there are no in-duced magnetic current sources. Along the same path of the derivation of the homoge-neous electromagnetic fields, asserting by definition that for the scaered fields we have Jext(x, ω) = 0and Kext(x, ω) = 0in (2.18) and (2.19), we obtain for the scaered

(21)

ﬁelds

−∇ × Hsct

(x, ω) + σbEsct(x, ω) =− [σ(x, ω) − σb] E(x, ω), (2.36)

∇ × Esct_{(x, ω)}_{− iωµ}_0Hsct_{(x, ω) = 0,} _(2.37)

where in view of (2.18) we interpret the right-hand side of (2.36) as an electric current density source.

In analogue procedure to the derivation for the electromagnetic ﬁelds in a homoge-neous medium, we introduce the vector potential equations such that

Esct(x, ω) = σ−1_b (k2_bAsct+∇∇ · Asct(x, ω)), (2.38)

Hsct(x, ω) =∇ × Asct. (2.39)

e vector potential is deﬁned similarly as in the homogeneous case by

Asct(x, ω) = ∫

Ω

g(x− x′, ω) [σ(x′, ω)− σb(ω)] E(x′, ω) dx′. (2.40)

e electric-ﬁeld integral equation

To define the field equations, we first introduce the electric conductivity contrast

χ(x, ω) = σ(x, ω)− σb(ω)

σb(ω) , (2.41)

which is zero wherever the anomalous conductivity σ(x, ω) is equivalent to the

back-ground conductivity σb(ω)outside Ω, as in ﬁgure 2.1(b). For convenience the

normal-ized vector potential for the scaering problem is usually deﬁned as

ˆ

Asct(x, ω) = ∫

Ω

g(x− x′, ω)χ(x′, ω)E(x′, ω) dx′. (2.42) e term χ(x, ω)E(x, ω) is called the contrast source. Once the incident electric field is evaluated for the external source distribution, such as (2.34), we add this field to (2.38) and obtain the total field as

E(x, ω) = Einc(x, ω) + k2_bAˆsct(x, ω) +∇∇ · Âsct(x, ω). (2.43) From (2.42) we observe that the vector potential ˆAsctonly depends on the electric field. For a given medium χ(x, ω), the electric field E(x, ω) is the only unknown in (2.43). In terms of the electric field we write the integral equation explicitly as

Einc(x, ω) = E(x, ω)−(kb2+∇∇· )∫

Ω

(22)

Once an (approximate) solution to (2.44) is obtained, the electric ﬁeld E(x, ω) is substituted in the right-hand side of (2.40). e resulting vector potential is then used by (2.39) to obtain the expression for the scaered magnetic ﬁelds Hsct_(xR_{, ω)}_{, expressed}

through the normalized scaering vector potential ˆAsct_{(x, ω)}_as

Hsct(xR, ω) = σb(ω)∇ × ˆAsct(xR, ω). (2.45)

2.2 F 

To simplify notation, we introduce the electric scaering operator Gχ operating on a

vector ﬁeld u(x, ω) as

Gχu(x, ω) =

(

k_b2+∇∇·)∫ x′∈Ω

g(x− x′, ω)χ(x′, ω)u(x′, ω) dV (2.46) and the identity operator I with Iu = u. en we can write the integral equation (2.44) in operator notation as

(I− Gχ)E(x, ω) = Einc(x, ω), x∈ Ω. (2.47)

To find the electric field E(x, ω) that (uniquely) satisfies this equation, we have to resort to numerical techniques, as analytic solutions are generally not available. A numeri-cal method essentially is to find an approximation in the finite-dimensional subspace VN ⊂ V on Ω that has minimal distance to the true solution in V for some metric in

V. Let us deﬁne at set of functions{ϕ1,k, . . . , ϕN,k}, k = 1, 2, 3, each associated with

a particular centre coordinate xnsuch that Vnwill consist of local functions that form

a basis in Vn. Each expansion function ϕn,k(x) = ϕk(x; xn)has a vector component

subscript k to indicate on which component the function is operating. Moreover, where appropriate, we use the vector valued function

ϕ_n,k(x) = ϕ_k(x; xn) = ϕk(x; xn)ik. (2.48)

We deﬁne the centre coordinates xn ∈ Ω associated with ϕn,kon a Cartesian grid

with mesh size h (see ﬁgure 2.2). For the local functions to span a volume equivalent to

Ω, we choose basis functions ϕn,kwith property

∫ Ω

ϕk(x; xn) dx = h3, k = 1, 2, 3. (2.49)

With these basis functions any function EN(x)∈ Vncan be wrien as

EN(x, ω) = 3 ∑ k=1 N ∑ n=1 en,k(ω)ϕn,k(x)ik. (2.50)

(23)

2.2. Finite-dimensional representation 19 ..xn . h . h . Ω

Figure 2.2: Two dimensional cross section of the grid{x1, . . . , xN} in Ω. e mesh

size h and an arbitrary coordinate xnare indicated for illustrative purpose.

To obtain a representation for E(x, ω) in Vn, we need to ﬁnd the set of coeﬃcients

{e1,k, . . . , eN,k} that minimizes the residual

R(x, ω) = Einc(x, ω)− 3 ∑ k=1 N ∑ i=n [ en,k(ω)ϕn,k(x)− Gχen,k(ω)ϕn,k(x) ] , (2.51)

where (2.47) is used. A general way to solve this problem is to introduce the method of mean weighted residuals (see Boyd (2002), chapter 3)

⟨wn,k, R(x, ω)⟩_Ω= 0, ∀ wn,k ∈ Wn; k = 1, 2, 3, (2.52)

where the N vector-valued weighting functions are deﬁned

wn,k(x) = wk(x; xn) = wk(x; xn)ik, k = 1, 2, 3, (2.53)

and⟨·, ·⟩Ω : V∗× V → R denotes the inner product on V , Wn ⊂ V∗ is the set of

weighting functions of size N and V∗_{is the dual space of V . is procedure is known}

as a weak formulation. e ﬁnite-dimensional representation for GχE(x, ω)consists

of three consecutive diﬀerent stages that need to be applied which we will detail next. e basis functions we shall use are

ϕk(x; xn) =    1 if|xj− xn,j| < 1 2h, ∀j = 1, 2, 3, 0 otherwise. (2.54)

(24)

(i) e contrast source

In the ﬁrst step, we require a ﬁnite-dimensional representation for the contrast source density

ˆ

J (x, ω) = χ(x, ω)E(x, ω). (2.55)

We deﬁne the ﬁnite-dimensional representation ˆJN(x, ω)in a similar way as the

ﬁnite-dimensional electric ﬁeld (2.50) such that

ˆ JN(x, ω) = 3 ∑ k=1 N ∑ n=1 jn,k(ω)ϕn,k(x)ik. (2.56)

Similar to (2.52), we obtain the coeﬃcients jn,kby minimizing the diﬀerence between

the vector ﬁeld ˆJ (x, ω)and the approximation, viz., ⟨ δ(x− xn)ik, χ(x, ω)EN(x, ω)− 3 ∑ k′=1 N ∑ n′=1 jn′,k′(ω)ϕn′,k′(x) ⟩ = 0, ∀ n = 1, . . . , N; k = 1, 2, 3. (2.57) As weighting function we have used the Dirac delta distribution, since the contrast is usually provided as a piece-wise constant distribution on Ω. Substituting with (2.50) we can evaluate the contrast source coeﬃcients jnexplicitly as

jn,k(ω) = χ(xn, ω)en,k(ω). (2.58)

(ii) e vector potential

e ﬁnite-dimensional representation for the vector potential is the approximation of the vector potential ˆA(x, ω)generated by a continuous convolution. It is approximated by a ﬁnite sum of basis functions

ˆ AN(x, ω) = 3 ∑ k=1 N ∑ n=1 an,k(ω)ϕn,k(x)ik ≈ ∫ Ω g(x− x′, ω) ˆJN(x′, ω) dx′, (2.59)

that is obtained through a similar procedure as for the contrast source ˆJN, though we

require some extra care due to the weak singularity in the Green’s function (2.30) at x = x′. As a weighting function, we consider the ball function around xn

Bn(x) =

{

1 if ∥x − xn∥ < 1₂h,

(25)

2.2. Finite-dimensional representation 21

Note that, like the weighting function ϕn,k, the value of Bn does not depend on k.

We write the residual problem for (2.59) with Bn(x)as weighting function, where the

objective is to minimize between the continuous expression and the ﬁnite dimensional approximation similar to (2.57), as ⟨ Bnik, ∫ Ω g(x− x′, ω) ˆJN(x′, ω) dx′− 3 ∑ k′=1 N ∑ n′=1 an′,k′(ω)ϕn′,k′(x) ⟩ Ω= 0, ∀ n = 1, . . . , N; k = 1, 2, 3. (2.61) If we restrict the support of ϕn,k(x)to completely contain the support of Bn(x), i.e.

Supp(Bn(x))⊆ Supp(ϕk(x; xn)), (2.62)

we immediately observe that aer the interchange of integration variables (2.61) can be wrien as an,k(ω) ∫ Bn(x) ϕk(x; xn) dx = ∫ x′_∈Ω N ∑ n′=1 jn′,k(ω)ϕn′,k(x′) ∫ x_∈Bn(x) g(x− x′, ω) dx dx′, (2.63)

where we have used the equality g(x− x′_{, ω) = g(x}′_{− x, ω). e integral on the}

le-hand side can be evaluated explicitly resulting in the constant value

η = ∫ Bn(x) ϕk(x; xn) dx = 4 3πh 3 _(2.64)

e integral (2.64) only diﬀers from (2.49) by the domain of integration. To carry out the integral over x on the right-hand side of (2.63) we use a spherical coordinate system. en we have

∥x − x′_{∥ =}(_{∥x − x}

n|2+∥x′− xn∥2− 2 cos(θ) ∥x′− xn∥∥x − xn∥

)1

2_, (2.65)

where∥x′_{− x}_n_{∥ is constant and ∥x − x}_n_{∥ = r is the variable of integration. e}

evaluation of the Green’s function integral over Bn(x) around the centre point xn,

where r <1 2h, yields ∫ Bn g(x− x′, ω) dx = ( 1− ikb1₂h ) exp(ikb1₂h)− 1 k2 b , x = x′. (2.66)

For integration over the sphere Bn(x)with∥x − x′∥ > 1₂h, we obtain

∫

Bn

g(x− x′, ω) dx = g(xn− x′, ω)ζ(kb, 1

(26)

in which ζ is given by

ζ(kb, a) = 2π

(iakb+ 1) exp(−iakb) + (iakb− 1) exp(iakb)

(ikb)3 , (2.68) with lim h→0 ζ(kb, h) h3 = 4 3π. (2.69)

Combining the results of (2.66–2.67) we introduce the spherical average Green’s func-tion, where we have included the volumetric scaling parameter η, by

˜ g(x− x′, ω) =        η−1 ( 1− ikb1₂h ) exp(ikb1₂h)− 1 k2 b , if x = x′, η−1g(x− x′, ω)ζ(kb, 1 2h), ∥x − x ′_{∥ >} 1 2h. (2.70) When∥x − x′_{∥ >} 1

2hand the discretization scale kbhis ﬁne enough such that lim kbh→0 η−1ζ ( kb, 1 2h ) = 1, (2.71)

the weak Green’s function ˜g(x−x′₎_{becomes equivalent to its strong form g(x}_−x′_{, ω)}_.

Combining the diﬀerent results we can rewrite (2.63) as

an,k = N ∑ n′=1 jn′,k ∫ x′∈Ω ˜ g(xn′− x′)ϕk(x′; xn′) dx′. (2.72)

is integral has no analytic solution. However, if we assume limited variation of ˜

g(xn − x′), we move the Green’s function out of the integral, approximating it by

a constant being the contribution at the midpoint x′_{= x}_n_{. Employing (2.49) we obtain}

the discrete convolutional expression

an,k≈ h3 N ∑ n′=1 ˜ g(xn− xn′, ω)jn′,k. (2.73)

Note that correctness of the limited variation assumption relates to the mesh size h. With all points xnare spaced equidistantly on an aﬃne grid, the sum (2.73) can be

evaluated using a discrete Fourier routine which generates all coeﬃcients {a1,k, . . . , aN,k}, k = 1, 2, 3 simultaneously per k. is speeds up computations considerably if

(27)

(iii) Spatial derivatives

e last step generates the ﬁnite-dimensional representation of the estimate for the scaered ﬁeld inside the domain Ω of observation

GχEN(x, ω) = 3 ∑ k=1 N ∑ n=1 esct_n,k(ω)ϕn,k(x)ik ≈ kb2Aˆ sct N (x, ω) +∇∇ · ˆA sct N (x, ω), (2.74) associated with the electric ﬁeld EN(x, ω). is can be evaluated by point-wise

evalu-ation on{xn} using a finite-difference scheme to obtain ∇∇ · ÂsctN (x, ω)(Abramowitz

and Stegun 1964), or we can use a weakening procedure by minimizing the difference between the scaered electric field obtained from the vector potential ˆAsct_N (x, ω)via (2.38) and the coefficient-based representation (2.74). is results in the residual equa-tion (Zwamborn and van den Berg 1992)

⟨ wn,k, 3 ∑ k′=1 N ∑ n′=1 [ esct_n′,k′ϕn′,k′− an′,k′ ( k_b2ϕ_n′,k′+∇∇ · ϕn′,k′ )] ⟩ = 0, ∀ wn,k∈ WN; k = 1, 2, 3. (2.75)

To ﬁnd a non-trivial solution, in this case the orthogonal basis functions ϕn,kand wn,k

need to basis functions ˜ϕn,k linear varying in space, with (2.49) still valid, such that

the intersect of adjacent cells on a mesh is non-empty. e details of this procedure for

∇∇·ϕn,kcan be found in the appendix. We note here that the used expansion functions

˜

ϕn,k will give the same result as the conventional ﬁnite-diﬀerence approximation (see

appendix A). e resulting scaering coeﬃcients can be formally wrien by

esct_n,k(ω) = ˆGχ({en′,k′(ω)}) , n, n′= 1, . . . , N ; k, k′= 1, 2, 3, (2.76)

where ˆGχ is the ﬁnite-dimensional counterpart of the operator embodying the three

steps described above.

Discrete linear system of equations

To find a solution for the discrete field equation (2.51), we consider the set of equations as weighted system of equations. We define its residual equation (2.52) as

⟨ wn,k, Einc(x, ω)− [ EN(x, ω)− ˆGχEN(x, ω) ] ⟩ Ω = 0, ∀ wn,k∈ WN. (2.77)

Equivalent to (2.57), we use the Dirac delta function as the weighting function. With this weighting function, we define the finite-dimensional representation of the incident electric field as

eincn,k(ω) =

∫ Ω

(28)

If we place the coeﬃcients en,k, eincn,k in vectors e and einc, respectively, we can write

the solution of (2.77) as the linear system of equations

( ˆI− ˆGχ)e = einc, (2.79)

where Î is the identity matrix and the scaered field coefficients are obtained from

(2.76).

Solving a system of coeﬃcients

While we have no need to determine the exact form of this matrix representation, by the Cayley-Hamilton theorem (Decell 1965) we know the solution to (2.79) can be wrien

as some M-th order matrix polynomial with coeﬃcients qmwith M ≤ N as

e≈ eM = M ∑ m=0 qm ( ˆ I− ˆGχ )m einc, (2.80)

where the relative importance of the higher order components will control the con-vergence. In case the coeﬃcients qm = 1, for all m, then the representation of (2.80)

reduces to the Neumann series. If the norm of ˆGχis less than 1 this series converges

for M → ∞.

By applying the three steps described in this section we can compute the product of ˆI− ˆGχ with a coeﬃcient vector. An iterative application of this operator creates

the polynomial. is is the basis of the Krylov subspace methods. An eﬃcient

imple-mentation to ﬁnd the polynomial coeﬃcients qmfor an orthogonalized set of vectors

( ˆI− ˆGχ)meincis the  algorithm by Saad (2003).

Retaining only the zero-order part (M = 0) we obtain the so-called Born approxi-mation. Although in practice the Born approximation is computationally advantageous the electric-field solution contains higher order terms. In the next chapter we show that we can make a more efficient approximation without sacrificing the computational simplicity.

Magnetic-ﬁeld evaluation

Once we have obtained a ﬁnite-dimensional representation EN(x, ω)that makes the

residual (2.51) suﬃciently small, the corresponding contrast source ˆJN(x, ω)is

com-puted via (2.58). Using (2.40) and (2.39) we can evaluate the ﬁnite-dimensional repre-sentation of the magnetic ﬁeld directly as

Hsct(x, ω) = σb∇ × ∫

Ω

g(x− x′, ω) ˆJN(x′, ω) dx′, x = xR. (2.81)

Note that if we treat xR_{as a parameter rather than a variable, the magnetic ﬁeld does not}

(29)

on the vector ﬁeld under the integral sign, and the integrand is evaluated as ∇ × g(x − x′_{, ω) ˆ}_J N(x′, ω) = g(x− x′, ω) ( ikb− 1 ∥x − x′_∥ ) Θ(x− x′)× ˆJN(x′, ω), (2.82)

with unit vector

Θk(x− x′) =

xk− x′k

∥x − x′_∥. (2.83)

As (2.81) is a weakly singular integral, we have

∂k

∫ x′∈Bϵ

1

4π∥xR_{− x}_′_∥dx′= 0 (2.84)

for a vanishing ball Bϵaround xR. If xR ∈ {xn}, there is no contribution from the

point xn = xR. Using (2.49) and taking the midpoint xn for all terms, we obtain

analogous to (2.72) and (2.73) the scaered magnetic ﬁeld as

Hsct(xR, ω) = σbh3 N ∑ n=1 xR̸=xn g(xR− xn, ω) × ( ikb− 1 ∥xR_{− x} n∥ )∑3 k=1 Θ(xR− xn)× jn,k(ω)ik. (2.85)

is completes the numerical description of the electromagnetic problem.

Summary

Induction logging systems use electromagnetic fields to probe their environment. Based on Maxwell’s equations we have described the pertaining physics of the classical elec-tromagnetic field and we have given specific equations to compute both the electric field and the magnetic field at any point in space for a conductive medium. As equations based on continuous vector spaces can only be solved for some very specific circum-stances, we have derived a finite-dimensional approximation to the integral form of the electromagnetic field equations.

e ﬁnite-dimensional approximation is obtained by iteratively applying a number

of computational steps. First, we compute the contrast sources ˆJN(x, ω)from some

estimate of the electromagnetic ﬁeld E(x, ω). en, we compute an approximate vec-tor potential ˆAsct

N (x, ω)from a convolution with the scalar Green’s function and the

contrast source. Subsequently we apply the gradient-divergence operator to this vec-tor potential. From the resulting diﬀerence in the linear matrix relation (2.79) we form

(30)

a new estimate for the electric field E(x, ω). Based on the previous estimate for the electric field a new estimate for the contrast source ˆJ (x, ω)is constructed, and the pro-cedure is repeated until the linear matrix relation agrees up to a predefined tolerance. From the approximated contrast source density we compute the magnetic field at the receiver positions xR_{by summing over all contrast source densities.}

e modelling of the electromagnetic ﬁeld will be required to analyse the behaviour of induction logging tools in realistic and complex media, the so-called forward problem, and are indispensable if from a physical measurement we want to reconstruct what kind of medium caused the measurements as they were observed, the so-called inverse problem.

(31)

Chapter ree

Induction logging and the Single Spherical Scaerer

approximation

Introduction

e computation of the electromagnetic response of an induction tool in an unbounded medium can be a cumbersome process. Measurements are taken multiple times per metre in the borehole, while induction logging tools have multiple external sources, each generating a new electromagnetic ﬁeld. Boreholes can have lengths of hundred of metres inside a target of interest, thus requiring thousands of controlled current-source locations in a single log, each generating a new electromagnetic ﬁeld problem consequently.

is chapter, we start with the derivation of a simpliﬁcation of the vector potential ˆ

Asct(x, ω)deﬁned by (2.42) such that we can circumvent the convolutional nature of

its computation in favour of a point-wise multiplication.

In view of the exponential decay of the diﬀusive electromagnetic ﬁeld, we analyse the replacement of the true conductivity medium on R3_{for a spatially restricted} win-dow on Ω⊂ R3_{. is local window of observation around the logging tool will decrease} the rank of the system of equations to a tractable amount.

We will use all aforementioned procedures to synthesize induction logging mea-surements on a reservoir medium, and compare results of our simpliﬁed model with the solution obtained from the integral equation.

3.1 T    

e iterative solution of the discretized version of integral equation (2.80), i.e. the so-lution of the electric field for a given medium and incident field, can be quite tedious. In particular, the number of operations required to obtain the electric field solution de-pends on the condition number of the matrix ˆI− ˆGχ (or the spectral radius for the

continuous equation of the system operator I− Gχ, see Arveson (2002)), which is

de-pendent on the contrast distribution χ(x, ω). For a multitude of evaluations of elec-tromagnetic ﬁelds extensive computer processing power is then required. Hence, an approximation is oen needed to capture the essential physics while keeping computa-tional times and/or system cost to a manageable extent.

(32)

From the diffusive electromagnetic field equations (2.18), (2.19) we have established that all sources, external and induced, generate fields with exponential amplitude damp-ing (cf. (2.30) and (2.31)). e exponential behaviour of the electromagnetic field results in high-order scaering terms existing primarily in a local seing. If we take the limit to this behaviour,the electric scaered field Esct _{at point x is generated only by the}

scaering source ˆJat the same position x. Starting with the expression in terms of the vector potential,

Esct(x, ω) = k2_bAˆsct(x, ω) +∇∇ · ˆAsct(x, ω), (3.1)

we consider the response from a restricted domain around x. erefore, we change the domain of integration Ω in (2.42) to that of a vanishing ball Bϵcentred at xϵwith radius

ϵ. e consequence is that the vector potential is approximated by ˆ Asct(x, ω)≈ ∫ Bϵ exp(ikb∥x − x′∥) 4π∥x − x′∥ ˆ J (x′, ω) dx′. (3.2)

In the vanishing limit ϵ→ 0, we may assume that ˆJ (x, ω) = χ(x, ω)E(x, ω)as well

as the exponential part of the Green’s function are constant over a vanishing ball. e integral then simpliﬁes to

ˆ Asct(x, ω)≈ [∫ Bϵ 1 4π∥x − x′∥dx ′]_exp(ikb_{∥x − x} ϵ∥) ˆJ (xϵ, ω). (3.3)

e integral can be evaluated with the aid of spherical coordinates to arrive at

ˆ Asct(x, ω)≈ ( 1 2ϵ 2₋1 6∥x − xϵ∥ 2 ) exp(ikb∥x − xϵ∥) ˆJ (xϵ, ω). (3.4)

It is readily observed that ˆAsct(x, ω)

x→xϵ disappears when ϵ → 0, i.e. for the ﬁrst

part of (3.1) we obtain ˆ Asct(x, ω) x_→xϵ =−1 6∥x − xϵ∥ 2_{J (x}_ˆ ε, ω)≈ 0. (3.5)

From (3.5) we easily compute the divergence in the limit x→ xϵ, viz.

∇ · ˆAsct(x, ω) x→xϵ ≈ − 3 ∑ k=1 1 3(xk− xk,ϵ)Jk(xϵ, ω)≈ 0, (3.6)

thereaer the gradient of (3.6) in the limit x→ xϵis expressed by

∇∇ · ˆAsct(x, ω) x→xϵ ≈ −1 3 ˆ J (xϵ, ω), (3.7)

(33)

3.2. Simpliﬁed expressions and sensitivity 29

Note that the result (3.7) can also be obtained by taking the integral of the divergence

instead of the other way around (see Lee et al. (1980)). Taking the limit ϵ → 0 (and

consequently with it x→ xϵ), we ﬁnally obtain from (3.1), (3.5) and (3.7) the

approxi-mation

Esct(x, ω)≈ −1 3

ˆ

J (x, ω). (3.8)

Using the explicit limit approximation (3.1) in the integral equation (2.43) for the electric ﬁeld E(x, ω) results in

E(x, ω)≈ Einc(x, ω)−1

3χ(x, ω)E(x, ω). (3.9)

From this approximation we state the Single Spherical Scaerer (SSS) approximation explicitly as

E(x, ω)≈ 3

3 + χ(x, ω)E

inc

(x, ω). (3.10)

is approximation is an improved version of the Born approximation (Bensdorp et al.

2014ba). If χ(x, ω)→ 0 we arrive at the standard Born approximation. e SSS

ap-proximation is still a ﬁrst order scaering apap-proximation.

3.2 S   

From a physical point of view, all electromagnetic phenomena take place in an un-bounded space. is unun-bounded space implies that the simulation of logging measure-ments would require the discretization of the complete space R3_{. e nature of diffusive} fields, where the electromagnetic fields decay exponentially in space, provides a rem-edy to this problem, as diffusive fields damp out exponentially fast. us the simulated

response on a subspace Ω⊂ R3_{converges exponentially to the response on R}3_.

3.2.1 Magnetic ﬁeld evaluation

Induction logging measures the magnetic ﬁeld H(x, ω) along multiple positions on the tool axis. In the forthcoming analysis, we shall use the SSS approximation (3.10) to obtain an approximate expression for the scaered magnetic ﬁeld Hsct_{(x, ω). If we}

substitute the scaered vector potential equation (2.40) in the deﬁnition for the scaered magnetic ﬁeld (2.39), we obtain

Hsct(xR, ω) =∇R× ∫

(x∈Ω)

(34)

where∇R ₌ _{∂R

1, ∂2R, ∂R3} denotes the spatial diﬀerentiation with respect to the re-ceiver coordinates. Substituting the SSS approximation for the electric ﬁeld E(x, ω) in (3.11) results into

Hsct(xR, ω) =∇R× ∫

x_∈Ω

g(xR− x, ω)κ(x, ω)Einc(x, ω) dx, (3.12) where, for convenience, we have introduced the modiﬁed contrast

κ(x, ω) = 3 χ(x, ω)

3 + χ(x, ω). (3.13)

To clarify in the relation between diﬀerent source-receiver polarizations, we substitute the incident ﬁeld Einc

(x, ω)in (3.11) with the magnetic-source formulation (2.33), i.e.

Einc(x, ω) =−iωµ0I0∇S× g(x − xS, ω)M , (3.14)

where∇S ₌_{∂S

1, ∂2S, ∂3S} denotes the spatial diﬀerentiation with respect to the source coordinates.

For computational purposes we shall assume I0 = (iωµ0)−1. en we carry out the

curl operations∇× operating on xS_{and x}R_{, respectively, to arrive at}

Hsct(xR, ω) = ∫

Ω

ξ(xR, xS, x, ω)Q(xR, xS, x)κ(x, ω) dx M . (3.15) All scalar parameters are collected by

ξ(xR, xS, x, ω) = ( ikb− 1 ∥x − xR∥ ) g(x− xR, ω) × ( ikb− 1 ∥x − xS_∥ ) g(x− xS, ω), (3.16)

while all terms dependent on their geometric angle are put in a single tensor

Q(xS, xR, x) = 1 ∆S_(x)∆R_(x) ×  −∆ S 2∆R2 − ∆S3∆R3 ∆S1∆R2 ∆S1∆R3 ∆S 2∆R1 −∆S1∆R1 − ∆S3∆R3 ∆S2∆R3 ∆S 3∆R1 ∆S3∆R2 −∆S1∆R1 − ∆S2∆R2   . (3.17) e ∆ functions indicating the distance between the position vector x and the pertain-ing source or receiver position vector. ey are given by

∆S(x) =∥x − xS∥, ∆R(x) =∥x − xR∥, (3.18)

∆S_k = xk− xSk, ∆

R