A Numerical 1 . 5 D Method for the Rapid Simulation of Geophysical Resistivity Measurements

Mostafa Shahriari 1,*, Sergio Rojas 2 ID , David Pardo 1,3,4, Angel Rodríguez-Rozas 1,4 ID , Shaaban A. Bakr 1,5, Victor M. Calo 2,6,7 and Ignacio Muga 8 1 Basque Center for Applied Mathematics (BCAM), Mazarredo 14, E48009 Bilbao, Spain; dzubiaur@gmail.com (D.P.); angel.rodriguez.rozas@gmail.com (A.R.-R.); shaaban.bakr1@gmail.com (S.A.B.) 2 Applied Geology, Western Australian School of Mines, Faculty of Science and Engineering, Curtin University, Perth, WA 6845, Australia; srojash@gmail.com (S.R.); vmcalo@gmail.com (V.M.C.) 3 University of the Basque Country (UPV/EHU), Barrio Sarriena, s/n, 48940 Lejona, Spain 4 Ikerbasque (Basque Foundation for Sciences), Maria Diaz de Haro 3, 6 solairua, 48013 Bilbao, Spain 5 Department of Mathematics, Assiut University, Assiut 71515, Egypt 6 Mineral Resources, Commonwealth Scientific and Industrial Research Organisation (CSIRO), Kensington, WA 6152, Australia 7 Curtin Institute for Computation, Curtin University, Perth, WA 6845, Australia 8 Instituto de Matemáticas, Pontificia Universidad Católica de Chile, Valparaíso 2350026, Chile; ignacio.muga@gmail.com * Correspondence: m.shahriari.sh@gmail.com


Introduction
There exist a wide variety of geophysical prospection methods.In this work, we focus on resistivity methods.We categorize these resistivity prospection methods according to their acquisition location as (a) on surface, such as the ones obtained using controlled sources electromagnetics (CSEM) [1][2][3] and magnetotelluric [4], and (b) in the borehole, such as the ones obtained using logging-while-drilling (LWD) devices.LWD is a technology that conveys borehole logging tools (e.g., gamma ray, resistivity, density, and sonic) downhole, record measurements, and transmit them to the subsurface for real time interpretation while the hole is being drilled [5][6][7][8][9][10][11][12].These tools provide two pieces of information: (a) real-time data, which is processed on the field while drilling, and (b) data that is stored in the device to process after pulling it out from the hole.We use real-time data to evaluate the formation for geosteering, which is the act of adjusting inclination and azimuth angles of the borehole to reach a geological target [5][6][7][8][9][10][11].
The first commercial LWD device appeared in the 1970s.They were commercially used for formation evaluation, especially in high-angle wells.Nowadays, LWD devices are used both for reservoir characterization [6][7][8] and geosteering applications [9][10][11].Modern borehole resistivity instruments can measure all nine couplings of the magnetic field, namely xx, xy, xz, yx, yy, yz, zx, zy and zz couplings (the first letter indicates the orientation of the transmitter and the second one indicates the receiver orientation) [13,14].
Since the depth of investigation of LWD resistivity measurements is limited compared to the assumed thickness of the geological layers, it is common to approximate subsurface models in the proximity of the logging instrument with a sequence of 1D models [5].In a 1D model, we reduce the dimension of the problem via a Hankel or a 2D Fourier transform along the directions over which we assume the material properties to be invariant [15][16][17].The resulting ordinary differential equations (ODEs) can be solved (a) analytically, leading to a semi-analytic approach after the numerical inversion of a Hankel transform [15,[17][18][19][20], or (b) numerically, leading to a numerical approach [16].
Solving the ODEs analytically has some major limitations, for example, (a) the analytical solution can only account for piecewise constant material properties and other resistivity distributions cannot be solved analytically, which prevents to properly model, for example, an oil-to-water transitioning zone when fluids are considered to be immiscible; (b) a specific set of cumbersome formulas has to be derived for each physical process (e.g., electromagnetism, elasticity, etc.) anisotropy type, etc.; (c) analytical derivatives of certain models (e.g., cross-bedded formations, or derivatives with respect to the bed boundary positions) are often difficult to obtain and have not been published to the best of our knowledge [13].
Solving the resulting ODEs numerically is also possible.The resulting method also exhibits a linear cost with respect to the discretization size, since it consists of a sequence of independent 1D problems.An example of this approach can be found in [16], where Davydycheva et al. use a 2D Fourier transform to reduce the dimension of the problem.Then, they employ a highly accurate 1D finite difference method (FDM) to solve the resulting ODEs.This method is relatively simple to code.However, this combined methodology requires a computational cost that is over 1000 times larger than that observed in semi-analytical methods.This occurs due to the elevated number of unknowns required to properly discretize the ODEs.If a common factorization of the system matrix based on a lower and an upper triangular matrix (the so-called LU factorization) would be precomputed for all source positions, the situation would only worsen: one would need a refined grid to properly model all sources, and since the cost of backward substitution is also proportional to the discretization size (as it occurs with the LU factorization), the large computational cost required to perform backward substitution would significantly increase the total cost of the solver.The use of other traditional techniques such as a finite element method (FEM) (see, e.g., [21]) to solve the resulting ODEs would not alleviate those problems.Additionally, in [16], the use of a 2D Fourier transform presents another burden on the performance of the solver, since the number of ODEs (and the total solver cost) increases quadratically with respect to the number of integration points of the 1D Inverse Fourier Transform (IFT), while in the case of a Hankel Transform, the cost only grows linearly with respect to the number of integration points on the 1D Inverse Hankel Transform (IHT).
The main objective of this work is to overcome the above limitations of both semi-analytic and existing numerical methods by solving each 1D problem (associated with a Hankel mode) using an efficient multi-scale FEM.Additionally, we seek a method with the following features: (a) we can consider arbitrary resistivity distributions along the 1D direction, and (b) we can easily and rapidly construct derivatives with respect to the material properties and position of the bed boundaries by using an adjoint formulation, which allows us to compute numerically the derivatives forming the Jacobian matrix needed by the Gauss-Newton inversion method at (almost) no additional cost.Despite these advances, presently our proposed multi-scale method is slower than the semi-analytic one.However, it is approximately two orders of magnitude faster than a traditional 1D FEM or a 1D FDM, like the one presented in [16].This speedup is essential for practical applications since one often needs to simulate thousands of logging positions and estimate millions of derivatives to solve a real-time inversion problems in the field for geosteering operations.Additionally, we employ a 1D Hankel transform rather than a 2D Fourier transform.This leads to a more complex mathematical formulation, but the resulting method exhibits a superior performance due to the lower number of ODEs that need to be solved.
In this work, we consider two traditional symmetric LWD tools available from oil service companies.Each of them contains two transmitters and two receivers.The receivers are symmetrically located at different sides of the transmitters [22].

3D Formulation
In this section, we describe the strong and weak formulations for our three dimensional (3D) electromagnetic (EM) problem.

Strong Formulation
Let σ(x, y, z) be the real-valued conductivity tensor with positive determinant.Let J be a volumetric current source density and M a magnetic volumetric current source.Then, the EM fields satisfy Maxwell's equations in 3D space [23]: ∇ where E and H are the time-harmonic (sinusoidal) complex-valued electric and magnetic fields, respectively, ω = 2π f is the angular frequency, where f > 0 is the frequency of the transmitter, ε and µ are the permittivity and magnetic permeability tensors of the media, respectively, and i is the imaginary unit, i 2 = −1.The problem domain is Ω = R 3 .By pre-multiplying (1) by σ−1 = (σ − iωε) −1 , applying the curl operator, and substituting (2) into the result, we arrive at the following reduced wave equation for the magnetic field: where the right-hand-side is: In order to ensure the uniqueness of the magnetic field, we use the Silver-Müller radiation condition [21,[24][25][26].

Weak Formulation
Let F be an arbitrary test function and F * its conjugate transpose.Pre-multiplying Equation (3) by F * and integrating over the domain Ω, we obtain the following equation: In the above, as a sufficient condition to ensure integrability, we select F ∈ H(curl; Ω), where: Using integration by parts assuming that the solution and its flux are continuous, and considering the radiation condition, we obtain the following variational formulation: (5)

1.5D Formulation
In our model problem, we consider the 3D Maxwell's equations in a 1D transversely isotropic (TI) layered formation.That is, the formation conductivity is constant along the x and y directions, and the conductivity tensor is defined as: where σ h > 0 is the conductivity of the media along the x and y directions, and σ v > 0 is the conductivity along z direction.Our formulation allows for parameter variations in the parameters along the z-axis.Analogously, ε(z) and µ(z) are considered to be transversely isotropic tensors.Since material properties are uniform in the xy-plane, it is convenient to use a Hankel transform to represent the magnetic field along x and y directions.

Hankel Transform
We consider H to be the 2D Fourier transform of H along x and y directions, where the material properties are homogeneous.We have: where x = (x, y) and k = (k x , k y ) (see Figure 1).We switch from the Cartesian system of coordinates to a cylindrical one according to the following transformations: Substituting ( 8) into (7) and applying the change of coordinates under the integral signs, we obtain: where ρ = (ρ, φ, z).Using the trigonometric identity: x We obtain: We now use the following relation between exponentials and Bessel functions: to obtain where We compute the cylindrical components of the magnetic field as follows: where Similarly for H φ , we have: By substituting (15) and ( 17) into (13), the Hankel representation of the magnetic field becomes: The curl of the magnetic field in cylindrical coordinates is: By substituting (18) into (19), the first component of the curl becomes: By using the property of Bessel functions given by Equation (A6) of Appendix A, we obtain: Using the formula of the derivative of the Bessel function given by Equation (A4) of Appendix A, we have: For an arbitrary function g(ξ, z) = (g − (ξ, z), g + (ξ, z), g z (ξ, z)) in the spectral domain, we introduce the following notation to simplify computations: Using (23), we obtain: For the second component of (19), using (18), we have: Using the property of the Bessel functions given by Equation (A4) of Appendix A and ( 23), we obtain: The third component of ( 19), becomes after using ( 18): For the derivative of J m+1 , we use Equation (A6) of Appendix A, and for the derivative of J m−1 , we use Equation (A4) of the same appendix.Then, We use Equation (A1) of Appendix A to simplify J m−2 .As a result, we obtain: 3.2.Hankel Finite Element (HFE) Full Field Formulation L 2 -orthogonality holds for Bessel functions of with the same order (see Equation (A3) of Appendix A).Hence, in order to simplify the terms of the variational formulation containing Bessel functions, we introduce the following matrix: Q is a unitary matrix, since: Hence, the change of coordinates implied by Q preserves the inner product.In particular, for arbitrary vector-valued functions U and V, we have: By using the above matrix, we obtain the following equalities: where σh = σ h − iωε h and σv = σ v − iωε v .For the L 2 terms, we obtain: For a specific Hankel mode ξ q > 0 and an exponential order t, we select a mono-modal test function of the form: where: Using (31) and the test functions defined in (34), we have: where Separating the integrals according to each variable and using the L 2 -orthogonality property of exponential functions, we have: By the orthogonality property of the Bessel functions given by Equation (A3) of Appendix A, we obtain: Similarly, for the other components of the curl, we have: For the L 2 terms, using (32) and the test functions defined in (34), and using L 2 -orthogonality property of exponentials and orthogonality property of the Bessel functions given by (A3) of Appendix A, we obtain: Using (37)-( 40), for each Hankel mode ξ q > 0 and exponential order t, the stiffness matrix becomes: where (42) Symbol < ., .> L 2 represents the L 2 inner product given by: A sufficient condition to guarantee that the above integrals are finite is to require In ( 43), the weak derivative of the function is considered.

Load vector
In 1D layered medium, we consider (0, 0, z Tx ) to be the general representation of a point source location.We use the following identities to describe the right-hand-side vector in cylindrical coordinates: where x, ŷ and ẑ are the unitary vectors in Cartesian coordinates.The right-hand-side of (3) in cylindrical coordinates for a z-oriented point source is: where δ is the Dirac delta distribution.We consider l to be the right-hand-side of (5).Using F q,t as our test function and separating the integrals according to each variable, we obtain: By L 2 -orthogonality of the exponentials, the load vector is non-zero when t = 0. Since J 0 (0) = 1, for a z-oriented point source, the right-hand-side becomes: Hence, we obtain the field by solving the following variational formulation: Lets consider R x and R y to be the right-hand-sides of (3) for x and y-oriented sources, respectively.By using (44), we obtain: Similarly to (46), the right-hand-side of the variational formulation is non-zero only when t = −1, 1.For those values of t, and for each Hankel mode, we have: where, Based on (18), we define the following notation: We have H t = (H t ρ , H t φ , H t z ).Therefore, the magnetic field for the x-oriented source is: For a y-oriented source, the field is computed as:

Multi-Scale Hankel Finite Element Method (Ms-HFEM)
We now describe our multi-scale FE method in the Hankel domain.In order to make the computational problem tractable, we truncate our domain along the z direction.We consider Ω z = (z 0 , z N ) to be our problem domain along z direction and we have −∞ < z 0 and z N < ∞.Moreover, we consider our solution to satisfy a zero Dirichlet boundary condition at both ends, since the waves amplitude rapidly decreases as we move away from the source.Thus, we have , with: In the following, for simplicity, we shall remove symbols ξ q and t from the notation.For each Hankel mode, we need to solve three problems associated with t = −1, 0, 1.The curl operator is the one defined in ( 24), ( 26) and (29).Similarly, Π + , Π − , and Π z are the symbols defined in Equation ( 23), and l is the right-hand-side of the variational formulation described in Equations ( 46) and (50) for t = −1, 0, 1.Our multi-scale approach consists of the following steps for each Hankel mode: 1.
Divide the domain into a finite number of sub-domains.We consider z 0 , 2).We call them decomposition points.We define each sub-domain as Ω i = (z i−1 , z i ), and we have: are the sub-domains associated to the decomposition points.

2.
Decompose the magnetic field into primary and secondary fields.For each Hankel mode, we decompose our magnetic field as follows: where H P and H S are primary and secondary fields, respectively.

3.
Find a local primary field.Lets assume that z Tx ∈ Ω p (see Figure 3).We define our local primary field H P ∈ V 0 (Ω p ) as the one that satisfies: Extending the local primary field to the entire domain with zero, we have H P ∈ V 0 (Ω z ).The local primary field has a discontinuous flux at z p−1 and z p .For the special case when the source is located at one decomposition point, we consider z Tx = z p , where The 10 −5 displacement is an arbitrary choice.We use a small number compatible with our grid size, but other values are also valid.

4.
Solve (N − 1) pairs of local problems.We consider . For each sub-domain Ω M i , we solve a pair of local problems which correspond to a discontinuous flux at the node z = z i .Specifically, the flux of the first local problem has a jump equal to 1 1 = (1, 0, 0), and the flux of the second local problem has a jump equal to 1 2 = (0, 1, 0).The local functions ψ i,k ∈ V 0 (Ω M i ) solve the following variational problems: where , and l M,k corresponds to the jump of the flux of the solutions at z = z i .Specifically: Similarly to the local primary field, we consider the extension by zero of the local solutions on Ω z , and we have ψ i,k ∈ V 0 (Ω z ).We denote the solutions of these local problems as multi-scale basis functions.We define the following space of multi-scale basis functions: Solve the secondary field formulation using the multi-scale basis functions.Since the flux components of the local primary field are discontinuous, we need our secondary field to balance these artificial discontinuities.Thus, by combining the primary and secondary fields, we recover a continuous flux for the full field.From (56), we obtain: We describe our secondary field as follows: where H S = (H S − , H S + , H S z ), and ψ i,k ∈ V M .By the definition of the multi-scale basis functions, (63) satisfies the reduced wave equation.Moreover, we consider F ∈ V M .Therefore, by substituting (63) into (62), we have: Finally, we add the local primary field and the secondary field to evaluate the full field.
In the next sections, we further describe the formulation for each step.

Local Primary Field
We consider the local primary field defined in Equation (57).We further decompose our primary field as follows: where H F is the fundamental solution of the electromagnetic reduced wave equation, and H C is a correction field.Since the fundamental field has no boundary condition at the boundaries of Ω p , the correction field is intended to enforce the zero Dirichlet (tangential) boundary condition on H P (see Figure 4).We define our correction field as follows: where H C,i,k (i = p − 1, p and k = 1, 2) are the correction basis functions.To define them, we perform the following decomposition: where H C,i,k 0 ∈ V 0 (Ω p ), and H C,i,k 1 ∈ V j (Ω p ), for j = i, is a lift of the correction field at z i to impose the non-zero Dirichlet boundary condition.V i (Ω p ) ⊂ V(Ω p ), for i = p − 1, p, is the space of all vector-valued functions F ∈ V(Ω p ) satisfying a zero Dirichlet (tangential) boundary condition at z = z i .By substituting (67) into variational formulation (41), we arrive at: (68) In order to impose a zero Dirichlet (tangential) boundary condition for our local primary field, we enforce the following conditions: where n p−1 = − ẑ and n p = ẑ are the outward unit normal vectors at z p−1 and z p , respectively.Using (66), we obtain:

Secondary Field Formulation
We define the secondary field to be the difference between the full field and the local primary field.Therefore, we have: Since the flux of the local primary field may be discontinuous on the boundaries of its domain, the flux of the secondary field should be discontinuous on the boundaries of the primary field's domain as follows: [ where n p−1 = − ẑ and n p = ẑ.Thus, the full field has a continuous flux on Ω z , since the secondary field satisfies (72).Hence, the secondary field formulation is: We have: Using the fact that H P (z − p−1 ) = 0, and considering (74), we obtain: Using the first component of the curl (24), we have: Similarly, since H P (z + p ) = 0, the jump on the right boundary of primary field's domain is equal to: By using the orthogonality of the Bessel functions and the exponentials, similar to the computations for (41), we obtain: where Similarly, by considering the jump condition of the flux of multi-scale basis functions, they are computed using (60).

Global Problem
Using (63), (79), and the definition of the multi-scale basis functions, we obtain: where, for an arbitrary test function F and a multi-scale basis function ψ i,k : We consider test functions F ∈ V M .Hence, by using (80), we have: By using (23), we can further simplify m(., .)to:

Implementation
For simplicity and to compare our numerical method directly with the state-of-the-art analytic implementations [15], we assume µ = µ 0 I 3 and ε = ε 0 I 3 (I 3 is the 3D identity matrix) to be constant in each layer.ε 0 is set to 8.85 × 10 −12 (F/m), which corresponds to the free-space permittivity, while µ 0 is set to 4π × 10 −7 (H/m), i.e., the magnetic permeability constant.
We consider each layer as a sub-domain.Therefore, the decomposition points are the boundaries of our layers.By doing so, we can evaluate the local primary fields in sub-domains which have smoothly varying materials.In particular, if we assume that the layer properties are homogeneous, the fundamental solution in (65) is independent of the tool position.Moreover, the correction basis functions are independent of the tool position.Consequently, instead of solving one primary field for each tool position, we find one fundamental field and four correction basis functions per layer.This simplification allows us to increase the speed of the method almost by a factor equal to the number of tool positions.
In our model problem, we consider two different Cartesian coordinate systems: (a) a system of coordinates related to the Earth, and (b) a system of coordinates related to the logging device, which consists of a rotation of the Earth system of coordinates in a way that the logging device extends along the z direction.We denote the angle between the logging instrument and the z direction of the Earth system of coordinates as α (relative dip angle).β is the azimuthal angle (see Figure 5).Therefore, the transformation between the systems of coordinates of the Earth and the tool is given by the following rotation matrix: where rotation R is defined by the following composition of rotations: In this notation, subscripts e and b denote the Earth system of coordinates and the logging instrument system of coordinates, respectively.If the logging instrument is perpendicular to the layering of the medium, we have α = 0.
In this work, for simplicity, we consider β = 0. Hence, the possibly non-zero couplings of the magnetic field can only be xx, xz, yy, zx and zz, where the first and the second letters in the subscript indicate the transmitter and receiver directions, respectively.
In this work, we consider LWD tools equipped with magnetic dipoles sources.We consider a traditional symmetric LWD instrument, which is a standard tool similar to those offered by oil service companies containing two transmitters and two receivers.In the aforementioned instrument, the receivers are located symmetrically around the transmitters (see Figures 6 and 10) [22].
To analyze the result of our experiments and compare them against those typically obtained in borehole resistivity applications, we further postprocess the values of the magnetic field.First, we evaluate the zz coupling (H zz ) of the magnetic field at two different receivers.We denote these values as H zz (Rx 1 ) and H zz (Rx 2 ), which correspond to the first and the second receiver, respectively.Then, we compute: where ph denotes the phase of a complex number.Often, the attenuation is defined in the decibel scale (which corresponds to the above number multiplied by 20).Subsequently, we compute the relation between attenuation and resistivity in a homogeneous media.This transformation, when applied to a heterogeneous media, delivers the apparent resistivity based on attenuation (see [5]).We similarly define the apparent resistivity based on the phase difference.
For the inverse Hankel transform, we use a fast Hankel transform algorithm based on digital filters (see [19] for details).
The entire Hankel Finite Element method has been implemented using FORTRAN 90.We consider the logging instrument described in Figure 6. Figure 7 through Figure 9 show the apparent resistivities (logs) for different dip angles.The distance between two consecutive logging points is half a foot (0.1524 m), as this is the resolution provided by most commercial logging-while-drilling (LWD) devices.In all cases, our simulation framework delivers a perfect agreement with the semi-analytic solutions.The large apparent resistivity values observed in Figures 8 and 9 form the so-called horns.They are a typical artifact that appears due to the employed post-processing method [10,16].Moreover, they do not affect significantly to the interpretation of the results.

Model Problem 2: Multilayered Formation
Figure 10 describes the logging instrument employed for this example.As in our previous examples, the distance between two consecutive logging points is half a foot (0.1524 m).The main result in terms of apparent resistivities for this model problem for a vertical well is shown in Figure 11.The attenuations and phase differences for this case are shown in Figures 12 and 13, respectively.As before, the numerical solutions coincide with the semi-analytic ones.Figure 14 compares the average computational time needed to solve one tool position using: (a) a traditional FE method, (b) a traditional FE method when we reuse the LU factorization, and (c) our proposed multi-scale method.As shown in the figure, reusing the LU factorization worsens the situation since we need a more refined grid to model for all source positions accurately.Table 1 shows a time comparison between a semi-analytic method and the proposed numerical method.For one tool position, computing the local primary field and pre-computing the multi-scale basis functions is computationally expensive and we observe a significant discrepancy between numerical and semi-analytic solutions.However, as expected, by increasing the number of tool positions, the ratio between the time of the proposed numerical method and the semi-analytic one is decreasing.Figure 15 shows the average time used to solve one tool position.Using a multi-scale method, the average time per position rapidly decreases as we augment the total number of tool positions.This occurs because the pre-computed multi-scale basis functions, fundamental fields, and correction basis functions only need to be computed once for any number of tool positions.The results in terms of apparent resistivities for this model problem for 60 • and 89 • deviated wells are shown in Figures 16 and 17, respectively.Again, the numerical and semi-analytic solutions coincide.).
Table 1.A time comparison (in seconds) for a vertical well using a semi-analytic method and the proposed multi-scale FEM as a function of the number of tool positions.
Positions Semi-Analytic Numerical Ratio We consider the logging instrument described in Figure 6 and a model problem that exhibits a sub-domain with a linearly varying resistivity, which is the case of an oil-to-water saturated transitioning (OWT) zone corresponding to a two phase-flow of immiscible fluids.To compare our solution to the semi-analytic one in the later case, we approximate our linearly varying resistivity distribution using multiple piecewise constant resistivity distributions.Figure 18 shows that the apparent resistivity of the semi-analytic solution is converging to the multi-scale solution as the number of layers with piecewise constant resistivity distribution increases.With our numerical method, we can model those arbitrarily varying conductivities without the need of approximating them.Figure 19 compares the attenuation and phase difference of the ultimate semi-analytic solution and the multi-scale one for the zz coupling.

Discussion and Conclusions
We propose a multi-scale Hankel FEM for solving Maxwell's equations in a 1D Transversely Isotropic media excited by a 3D arbitrarily oriented point dipole.The multi-scale FEM pre-computes the fundamental fields, and correction and multi-scale basis functions.As a result, this computation is expensive if only a single logging position is studied, but it becomes competitive as the number of logging position grows.
The numerical method produces highly accurate solutions, as our numerical validations experiments show.Additionally, computation of the parametrization derivatives is straightforward by simply considering the adjoint formulation.By using this method, it is possible to consider arbitrary resistivity distributions along the z direction, while semi-analytic methods only allow for piecewise constant material coefficients.
The method we propose is still slower than the semi-analytic one.The most time consuming part of our method is to compute the primary field.However, in the case of piecewise constant resistivity distribution, we can obtain the fundamental field analytically, which decreases the cost of computing the primary field considerably.Nonetheless, this work employs a full FEM implementation to preserve the generality of the method.This enables us to consider arbitrary resistivity distributions (other than just piecewise constant ones).Because of the lack of existence of a semi-analytic method that can model the aforementioned case, our methodology is a viable alternative.In order to use a semi-analytic method, we can approximate the non-constant resistivity distribution using multiple constant resistivity distributions.However, this piecewise-constant representation leads to an extensive error and a bothersome implementation when computing the derivatives to form the Jacobian matrix to perform the inversion.
As future work, we plan to extend our method to other multi-physics problems e.g.elasto-acoustic problems and to account for other material parameter distributions, such as, cross-bedded formations.Also, we will investigate the optimal selection of the primary field.Moreover, we shall also integrate this simulator on an inversion software platform.A trivial parallelization of this software through operating frequencies and transmitters of a logging device and logging positions is straightforward.Additionally, the parallelization of the method through multi-scale basis functions, fundamental fields, and correction basis functions is under development, which will boost the speed of the method.

Figure 1 .
Figure 1.Cartesian and cylindrical systems of coordinates.

Figure 4 .
Figure 4. Fundamental, correction and local primary fields.

Figure 5 .6. Numerical Examples 6 . 1 .
Figure 5. x e , y e and z e are the axes of the Cartesian coordinates.x b , y b and z b are the axes of the borehole coordinates.α and β are the dip and azimuthal angles, respectively.

Figure 6 .Figure 7 .
Figure 6.Logging instrument for model problem 1. Tx 1 and Tx 2 are the induction transmitters, and Rx 1 and Rx 2 are the receivers.

Figure 8 .
Figure 8. Model problem 1. Apparent resistivity for the zz coupling for a 70 • deviated well.

Figure 9 .
Figure 9. Model problem 1. Apparent resistivity for the zz coupling for a 89 • deviated well.

Figure 10 .Figure 11 .
Figure 10.Logging instrument for model problem 2. Tx 1 and Tx 2 are the induction transmitters, and Rx 1 and Rx 2 are the receivers.

Figure 12 .Figure 13 .
Figure 12.Model problem 2. Attenuations for the non-zero couplings (xx and yy couplings) of the magnetic field for a vertical well (dip angle = 0 • ).

Figure 14 .Figure 15 .
Figure 14.Model problem 2. Average time (in seconds) to solve for one tool position as the number of logging positions varies using (a) a traditional FE approach, (b) a FE approach where we reuse the LU factorization, and (c) our proposed multi-scale method ( time(seconds) number of tool positions).

Figure 18 .
Figure18.Model problem 3. Convergence of the zz coupling apparent resistivity of the semi-analytic solution to the multi-scale one for a vertical well (dip angle = 0 • ).l is the number of piecewise constant resistivity distributions that we select to approximate the linearly varying resistivity.

Figure 19 .
Figure 19.Model problem 3. Attenuation and Phase difference of the magnetic field for a vertical well (dip angle = 0 • ).