Computational Modeling of Geoelectrical Soundings using PML-FDTD

The Finite-Difference Time-Domain (FDTD) method was applied in order to analyze the transient responses of geoelectrical soundings that use circular electric dipole (CED) as source over stratified formations. The model was developed in cylindrical coordinates and a perfectly matched layer (PML) was incorporated to the domain to absorb wave reflections at computational grid boundaries. Numeric results are validated with analytic solutions. Comparisons between the transient response of two different type of soundings are performed and results indicate that the transient response of soundings that excite purely TM mode are more sensitive to the variation of electrical characteristics of the medium.


INTRODUCTION
Fast and efficient numerical simulations of geoelectrical and geoelectromagnetic sounding techniques are of great importance for Oil & Gas industry.Besides assisting the process of exploration of hydrocarbon reservoirs, they eliminate significant costs with new prototypes development and field tests execution.
Due to the enhance of computational capability and new numerical techniques development, different types of soundings have been computationally modeled, and several have excelled.Some of them are the Marine Controlled Source Electromagnetic (mCSEM) [1]- [7], which has been used successfully to complement deep and ultradeep reservoir exploration seismic activities, and Loggingwhile-drilling (LWD) tools modeling, that is very useful for direction and/or horizontal wells geodirectional drilling.[8]- [19].
At the end of last decade, a new kind of geoelectrical sounding, that applies a Circular Electric Dipole (CED) as excitation source, was proposed [20].The authors developed a unidimensional analytical model (1D) and presented an amount of field tests conducted by the measurement tool patent-owner company.CED is a sort of sounding that excites transverse magnetic (TM) fields in the ground, which results in a transient response to a pulse excitation strongly dependent on medium constitutive parameters.Recently, new field tests have been held in order to evaluate the potential of using this type of sounding to identify resistive layers that represent hydrocarbon saturated formations [21]- [23].
In this work, a computational model that applies Finite-Difference Time Domain (FDTD) method is where and are the electric and magnetic field intensity vectors, respectively; and are the electric and magnetic current density vectors of the source, respectively; ϵ, σ and μ are medium's electric permittivity, magnetic permeability and electric conductivity, respectively.Equations ( 1) are subject to the following boundary condition: ( where is the domain boundary.

A. Discretization
The proposed computational modeling applies the Finite-Difference Time-Domain (FDTD) method and an interlaced spatial grid scheme in cylindrical coordinates, decomposing the domain in two groups of cells, known as primary and dual cells.Fig. 1 illustrates spatial positioning of electric and magnetic fields in a three-dimensional (3D) cylindrical cell.A cylindrical coordinates system was chosen to avoid approximation errors while modeling the sources.Primary and dual cells grid vertices are defined in (ρ i , ϕ j , z k ) and (ρ i +Δρ/2, ϕ j +Δϕ/2, z k +Δz/2), respectively.The nodal primary grid indexes are i, j and k, while Δρ, Δϕ and Δz are the spatial increments in radial, azimuthal and longitudinal directions, respectively.At the time domain, electric and magnetic field components are interleaved at intervals Δt/2, which is the time increment.and components are located at the cylindrical grid so that field continuity conditions at the interface between media with different materials can be guaranteed.

1) Circular Electric Spiral Source
Applying a circular electric spiral (loop) as excitation source of the system and taking into account that the investigated medium is stratified, due to the symmetry of the scenario, there is no change in field intensity over ϕ-axis.Therefore, the three-dimensional (3D) model can be reduced to a bidimensional (2D) one.In this case, the cylindrical cell considered for the problem discretization can be illustrated as shown in Fig. 2. Therewithal, this type of antenna creates purely TE fields, with the components E ρ , E z and H ϕ equal to zero, and the equations that define the problem are: Centered finite-difference approximation to the time and space partial derivatives can be applied in equations ( 3), so: It is important to note that components are defined in integer time instants, n.

2) Circular Electric Dipole Source
As introduced earlier, the circular electric dipole (CED) will be represented in this work by a circular magnetic spiral (not feasible in practice) that, in the same way as CED, creates purely TM fields, but also simplifies computational effort to bi-dimensional (2D) modelling, due to source symmetry.In this case, it is possible to apply the Duality Theorem to Maxwell's equations and modelling process is equivalent to the one done in last topic.Fig. 3 shows the cylindrical cell considered for FDTD discretization.It is important to observe that, at this point, magnetic field components are defined in primary grid of FDTD, while electric field components are in dual grid.Since this kind of source creates purely TM fields, the components H ρ , H z and E ϕ are equal to zero and the equations that define the problem are: Thus, the following discretized equations are obtained: In equations ( 6), components are defined in integer time instants, n.

3) Degenerated Cells
It is important to notice that, at the origin of the coordinate system (ρ=0), there is a singularity in equations (4b) and (6b).In order to solve this issue, the procedure adopted was to degenerate ρ-axisadjacent cells.A sketch of the degenerated cells with its field components is shown in Fig. 4. Because of the impossibility to calculate spatial derivatives at the origin, in this the integral form of Maxwell's equations are discretized.Thus, Ampère's law and Faraday's law can be expressed as: It should be noted that, as the electric field was defined in primary grid in the case of circular electric spiral source, it is necessary to calculate H z field only in i = 1/2.E ϕ and H ρ components are not defined in i = 0. Similarly, in case of using a circular magnetic spiral as source, the only component which value changes in i = 0 is E z .
Therefore, assuming that fields do not change inside the cell, the surface and line integrals in Ampère's and Faraday's laws can be easily discretized.Following the same discretization procedure applied previously, the discretized equation for the degenerated cells are: Cylindrical Perfectly Matched Layers A perfectly matched layer (PML) is incorporated in the outermost cells of the grid in order to absorb outgoing waves.This is done by modifying the constitutive parameters inside the PML region, and hence it does not require any modification in Maxwell's equations themselves.In cylindrical coordinates, the PML constitutive tensors that are matched to a homogeneous nondispersive medium characterized by constitutive parameters ϵ and µ, are given by and , with [24]: (9) where is the analytic extension of the radial coordinate to a complex variable domain, and s ρ and s z are complex stretching variables.These variables are defined as: and (11) where a ρ (ρ), a z (z), σ ρ (ρ) and σ z (z) are functions of position only [24].Inside the PML region, longitudinal eigenmodes turn into , where ; and, for radial eigenmodes, the same occurs in terms of Hankel functions.Hence, the transformed eigenmodes exhibit exponential decay inside the PML region so as to reduce spurious reflections from the grid terminations.
An effective PML design requires balancing the theoretical reflection error and the numerical discretization error.Several profiles have been suggested for grading in the context of PML.The most successful implementations use a polynomial or geometric variation of the PML loss.
Here, the following polynomial grading is adopted [25]: where , d is the PML thickness, m is the polynomial grading factor, is the PML conductivity at the outer boundary, θ is the incidence angle over PML and R(θ) is the theoretical reflection error.

III. NUMERICAL RESULTS
The basic configuration of the sounding operations scenario is shown in Fig. 5. Circular spiral sources (electric and magnetic) have 500 m radius and the excitation pulse width is 0.05 s.Unless mentioned otherwise, the ground layers conductivity is 0.2 S/m and fields are sampled in ρ = 250 m.For all simulations performed in this work, relative electric permittivity and relative magnetic permeability are equal to 1.The PML is set to have 4 cells and a quadratic polynomial profile for imaginary term of the stretching variable (s ρ ).The real terms of s ρ and s z were fixed at 1.
To validate the PML-FDTD algorithm, the results of a homogeneous ground simulation were compared to analytic solutions.For the validation, a circular electric spiral excitation source was adopted.The computational domain was discretized in a (N ρ , N z ) = (25, 25) standardized grid, where Δ ρ = Δ z = 80 m and the field amplitude is sampled in ρ = 800 m.Fig. 6 shows an acceptable similarity between PML-FDTD and analytic results.To compare the results obtained for circular electric spiral and CED (modeled as a circular magnetic spiral) sensors, the two types of excitation sources were implemented in a nonhomogeneous ground scenario.Fig. 5 illustrates the simulation domain, where there is air (σ = 0 S/m) above the ground and two formation layers with different conductivities.
At first, a convergence analysis of the PML-FDTD algorithm was executed through variation of cells size, maintaining medium total size.For this simulation, upper and lower formation layers conductivities were set to 0.2 S/m and 0.1 S/m respectively.Fig. 7 and Fig. 8 show transient results of the two sensor types, electric and magnetic spirals, respectively.It is important to note that the algorithm shows adequate convergence with spatial increment reduction.
In order to evaluate sensors sensibility, maintaining the upper layer conductivity constant at 0.2 S/m and varying the lower layer conductivity from 0.01 to 0.1 S/m, new simulations were performed.
Transient responses of the electric and magnetic spiral sensor types can be observed in Fig. 9 and Fig. 10, respectively.It is important to note that fields produced in formations by a magnetic spiral source (purely TM fields) are more sensitive to the variation of electrical characteristics of the lower layer than the ones created by the electric spiral (purely TE fields).

Fig. 5 .
Fig. 5. Basic configuration of the modeled medium, with excitation source and PML technique implemented.

Fig. 7 .Fig. 8 .
Fig. 7. Results variation due to cells size change, using a circular electrical spiral as excitation source.

Fig. 9 .
Fig. 9. Results obtained for the circular electrical spiral, when σ = 0 S/m above the ground, σ = 0.2 S/m at upper layer and lower layer conductivity varying from 0.01 to 0.1 S/m.

Fig. 10 .
Fig. 10.Results obtained for the circular magnetic spiral, when σ = 0 S/m above the ground, σ = 0.2 S/m at upper layer and lower layer conductivity varying from 0.01 to 0.1 S/m.
Ribeiro 1 and Marcela S.Novo 2developed to analyze the transient responses of geoelectrical soundings that use circular electric dipole (CED) as source over stratified formations.Knowing that CED excites purely transverse magnetic (TM) fields and that its model would be three-dimensional (3D), for this pioneer study proposed, CED will be represented by a circular magnetic spiral (not feasible in practice), which also that the algorithm is precise, robust and stable.As expected, results indicate that transient response of soundings that excite purely TM fields in the ground are more sensitive to the variation of electrical characteristics of the medium.