# Reduction of Electromagnetic Reflections in 3D Airborne Transient Electromagnetic Modeling: Application of the CFS-PML in Source-Free Media

Authors: Yanju Ji; Xuejiao Zhao; Jiayue Gu; Dongsheng Li; Shanshan Guan
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

To solve the problem of electromagnetic reflections caused by the termination of finite-difference time-domain (FDTD) grids, we apply the complex frequency-shifted perfectly matched layer (CFS-PML) to airborne transient electromagnetic (ATEM) modeling in a source-free medium. To implement the CFS-PML, two important aspects are improved. First, our method adopts the source-free Maxwell’s equations as the governing equations and introduces the divergence condition, consequently, the discrete form of Maxwell’s third equation is derived with regard to the CFS-PML form. Second, because our method adopts an inhomogeneous time-step, a recursive formula composed of convolution items based on a nonuniform time-step is proposed. The proposed approach is verified via a calculation of the electromagnetic response using homogeneous half-space models with different conductivities. The results show that the CFS-PML can reduce a 60 dB relative errors in late times. Moreover, this approach is also applied to 3D anomalous models; the results indicate that the proposed method can reduce reflections and substantially improve the identification of anomalous bodies. Consequently, the CFS-PML has good implications for ATEM modeling in a source-free medium.


---

## Body

## 1. Introduction

The airborne transient electromagnetic (ATEM) system, which is an economic alternative for acquiring electromagnetic data with highly efficient detection capabilities and a substantial depth of investigation, has been widely applied to problems associated with hydrogeological surveying, mineral exploration, and environmental monitoring [1–3]. Furthermore, high-accuracy ATEM modeling can provide a theoretical basis for subsequent data inversion and prospective instrument design. Numerous studies have been conducted using the three-dimensional (3D) TEM modeling approach developed by Oristaglio and Hohmann based on finite-difference time-domain (FDTD) grids [4]. More recently, Commer and Newman improved the accelerated simulation scheme for 3D TEM modeling using geometric multigrid concepts [5–7]. Subsequently, Guan calculated the 3D ATEM response based on the graphics processing unit (GPU) [8], and Sun et al. improved upon the 3D FDTD modeling of TEM in consideration of the ramp time [9]. However, boundary reflections still constitute one of the most important challenges for the accuracy of TEM modeling. The most widely used Dirichlet boundary condition (DBC) exhibits better effects only at earlier times before the diffusion field has arrived at the boundary [10]. With an increase in the computing time, the modeling accuracy is progressively affected by electromagnetic reflections at the boundaries. Therefore, it is necessary to develop a more effective boundary condition. Berenger proposed a perfectly matched layer (PML), which can absorb electromagnetic waves with any incident angle and frequency [11]. Subsequently, several approaches, such as the uniaxial perfectly matched layer (U-PML), multiaxial perfectly matched layer (M-PML), and complex frequency-shifted perfectly matched layer (CFS-PML), have been developed with many absorbing boundary conditions based on PML theory [12–14]. In these improvements, the CFS-PML absorbing boundary condition has shown the best performance with regard to the absorption of low-frequency induction fields and late reflections [15, 16]. Roden and Gedney more efficiently implemented the CFS-PML condition by utilizing it with the FDTD method based on recursive convolution; this condition, known as the convolutional perfectly matched layer (C-PML), is highly absorptive of evanescent modes, and it is independent of the host medium [17]. Drossaert and Giannopoulos effectively reduced the computation time and memory by using the complex frequency-shifted stretching function while applying the CFS-PML with the velocity-stress wave equations. Based on the FDTD technique [18], Giannopoulos proposed an electromagnetic modeling method using the CFS-PML with the complex frequency-shifted stretching function and analyzed the effects of the absorption coefficient in different media [19]. Gedney and Zhao applied the CFS-PML boundary condition to the simulation of two-dimensional (2D) electromagnetic waves based on an auxiliary differential equation (ADE), thereby improving the absorbing properties of the CFS-PML in their research [20]. Moreover, Li and Huang used the CFS-PML to calculate the TEM response in a source medium while discretizing the air [21]. Feng et al. reduced the memory requirements in an implementation of the CFS-PML based on the memory-minimized method (Tri-M) and adopted the discrete Zernike transform (DZT) to guarantee the optimal accuracy [22]. Yu et al. improved the C-PML parameters to absorb low-frequency electromagnetic waves in both the ground and the air and observed good absorption [23]. Furthermore, Hu et al. applied the CFS-PML within a fictitious wave domain and discussed the selection of various CFS parameters [24]. Zhao et al. improved the discretization of Maxwell’s divergence equation in the CFS-PML boundary based on a uniform time-step [25]. Feng et al. proposed a Crank–Nicolson cycle-sweep-uniform FDTD method based on CFS-PML and applied it to 3D low-frequency subsurface electromagnetic sensing problems [26].In this paper, we focus on the problem of electromagnetic reflections at boundaries in source-free media and apply the CFS-PML boundary condition to ATEM modeling based on a nonuniform time-step. Maxwell’s curl equation with a source was used as the governing equation in previous studies about CFS-PML, as the source functions are difficult to approximate and the analytical solution in the source medium is unavailable. In contrast to previous investigations, our research adopts the source-free Maxwell’s equations as the governing equations and the analytical solutions as the initial conditions to perform the iterative calculations; therefore, we can obtain modeling results with a high accuracy. Because Maxwell’s divergence equation must be incorporated into the governing equations to ensure the uniqueness of the solution [10], one of the fundamental aspects of our approach is the discretization of Maxwell’s divergence equation in the stretched coordinate space. Another crucial component of the proposed method is the derivation of a recursive convolution formula based on nonuniform time-steps in the discretized stretched coordinate formulation. Finally, the proposed method is verified via a calculation of the electromagnetic response using homogeneous half-space models and anomalous models with 3D bodies, thereby effectively demonstrating the performance of the proposed method with regard to the absorption of electromagnetic reflections during the modeling of ATEM.


## 2. Methods

Maxwell’s divergence equation is employed as the governing equation in FDTD-based ATEM modeling in a source-free medium to advance the magnetic field ofHz, while the other fields are advanced using Maxwell’s curl equations. Following the method proposed by Wang and Hohmann [10], the initial conditions are calculated based on a homogeneous half-space. The staggered Yee grid is used to discretize the earth model inhomogeneously. In addition, the time-steps are advanced using a modified Du Fort-Frankel method [4, 27]. The ATEM model contains a ground-air boundary and a computational domain boundary; therefore, we adopt an upward continuation at the ground-air boundary and the CFS-PML as the computational domain boundary condition. However, two problems are encountered during the application of the CFS-PML: we must derive the discrete form of Maxwell’s divergence equation in the stretched coordinate, and the recursive convolution formulation based on a nonuniform time-step must be derived. Consequently, these key issues in the application of the CFS-PML boundary condition are discussed in this paper.


### 2.1. Governing Equations in the Stretched Coordinate Space

In the modeling of ATEM in a source-free medium, Maxwell’s divergence equation must be included within the governing equations during the modeling of ATEM in a source-free medium to ensure the uniqueness and stability of the solution. Therefore, Maxwell’s divergence equation (equation (3)) is adopted to advance the magnetic field Bz in the stretched coordinate space [28–32], while the other electromagnetic fields are advanced using Maxwell’s curl equations ((1) and (2)).
(1)∇s×H=εjωE+σE,(2)∇s×E=−μjωH,(3)∇s·H=0,where E is the electric field, H is the magnetic field, jω is the complex frequency variable, μ is the permeability, ε is the permittivity, and σ is the electrical conductivity.In the stretched coordinate space, the Hamilton operator (∇s) can be expressed as follows [33]:
(4)∇s=x→1sx∂∂x+y→1sy∂∂y+z→1sz∂∂z,where si=ki+σpi/αi+jωμ, ki and αi are positive and less than 1 and σpi is the conductivity of the CFS-PML layers (i=x,y,z).The difference scheme of Maxwell’s curl equations in the stretched coordinate space was developed by Roden and Gedney similar to the application of the CFS-PML [17]. Thus, in this paper, we mainly deduce the discrete form of Maxwell’s divergence equation (equation (3)).Inserting (4) into (3) leads to
(5)1sx∂Hx∂x+1sy∂Hy∂y+1sz∂Hz∂z=0.Equation (5) is next transformed into the time domain to obtain the renewal equation of Hz, after which the expression of Si (i = x, y, z) is inserted into (5), which can consequently be expressed as
(6)1kx∂Hx∂x+1ky∂Hy∂y+1kz∂Hz∂z+ψhzx+ψhzy+ψhzz=0,where ψhzi is the auxiliary expression of a convolution item and is implemented as
(7)ψhzi=−σpiεki2exp−σpiεki+αiεt∗∂Hi∂i.The governing equations of other fields in the stretched coordinate can be similarly derived.


### 2.2. The Discrete Form of Maxwell’s Divergence Equation in the Stretched Coordinate Space

During ATEM modeling, the time-step of an FDTD method usually increases gradually according to the attenuation characteristics of TEM responses (i.e., the TEM field exhibits relatively sharp variations at earlier times and gradually becomes smooth thereafter). Here, the time-step is advanced using a modified Du Fort-Frankel method. The traditional recursive convolution formula is deduced based on a uniform time-step; therefore, we need to modify the derivation procedure based on inhomogeneous time-steps which is shown in the appendix. In the stretched coordinate space,ψhzx at tn+1/2 can be expressed as
(8)ψhzxn+12=−σpxεkx2exp−σpxεkx+αxεtn+1/2∗∂Hxn+1/2∂x=ξhxtn+1/2∗∂Hxn+1/2∂x.To obtain the discrete form of the convolution formulation, we assume thatΔtn=tn+1/2−tn−1/2 (n=0,1,2,…,n), and thus, equation (8) can be rewritten as
(9)ψhzxn+12=ξhxtn+1/2∗∂Hxn+1/2∂x=∫0tn+1/2ξhxτ∂Hxtn+1/2−τ∂xdτ,where ξhxt=−Axe−Bxt, Ai=σpi/εki2, Bi=σpi/εki+αi/ε (i=x,y,z).Thenψhzxn+1/2 can be performed recursively:
(10)ψhzxn+12=12AxBxe−BΔtn−AxBx∂Hxn+1/2∂x+∂Hxn−1/2∂x+e−BxΔtnψhzxn−12.We can obtainψhzzn+1/2 and ψhzxn+1/2 in the same manner:
(11)ψhzzn+12=12AzBze−BzΔtn−AzBz∂Hzn+1/2∂z+∂Hzn−1/2∂z+e−BzΔtnψhzzn−12,(12)ψhzxn+12=12AxBxe−BxΔtn−AxBx∂Hxn+1/2∂x+∂Hxn−1/2∂x+e−BxΔtnψhzxn−12.As shown in (6) and (11), ψhzz contains Hz at the current time in the discretization of Maxwell’s divergence equation to obtain the iterative formula of Hz at the time tn+1/2. Therefore, the iterative formula of Hz requires further derivation. Equation (6) can be discretized in both space and time according to the staggered Yee scheme as
(13)1kxHxn+1/2i+1,j+1/2,k+1/2−Hxn+1/2i,j+1/2,k+1/2Δx+1kyHyn+1/2i+1/2,j+1,k+1/2−Hyn+1/2i+1/2,j,k+1/2Δy+1kzHzn−1+1/2i+1/2,j+1/2,k+1−Hzn−1+1/2i+1/2,j+1/2,kΔz+ψhzxn+1/2+ψhzyn+1/2+ψhzzn+1/2=0.Inserting (10)–(12) into (13), we obtain
(14)Hzn+1/2i+12,j+12,k=Hzn+1/2i+12,j+12,k+1+CzDzΔzi+1/2,j+1/2,kHzn−1+1/2i+12,j+12,k+1−Hzn−1+1/2i+12,j+12,k+DzkyΔyi+1/2,j,k+1/2Hyn+1/2i+12,j+1,k+12−Hyn+1/2i+12,j,k+12+DzkxΔxi,j+1/2,k+1/2Hxn+1/2i+1,j+12,k+12−Hxn+1/2i,j+12,k+12+Dzψhzxn+1/2+Dzψhzyn+1/2+e−BzΔtnDzψhzzn−1/2,where Cz=1/2Az/Bze−BzΔtn−Az/Bz, Dz=kz/1+kzCzΔzi+1/2,j+1/2,k.Note thatHz at tn+1/2 is related to the value of ψhzz of the previous time-step, and thus, we must refresh ψhzz at the end of the calculation for the next iteration.


## 2.1. Governing Equations in the Stretched Coordinate Space

In the modeling of ATEM in a source-free medium, Maxwell’s divergence equation must be included within the governing equations during the modeling of ATEM in a source-free medium to ensure the uniqueness and stability of the solution. Therefore, Maxwell’s divergence equation (equation (3)) is adopted to advance the magnetic field Bz in the stretched coordinate space [28–32], while the other electromagnetic fields are advanced using Maxwell’s curl equations ((1) and (2)).
(1)∇s×H=εjωE+σE,(2)∇s×E=−μjωH,(3)∇s·H=0,where E is the electric field, H is the magnetic field, jω is the complex frequency variable, μ is the permeability, ε is the permittivity, and σ is the electrical conductivity.In the stretched coordinate space, the Hamilton operator (∇s) can be expressed as follows [33]:
(4)∇s=x→1sx∂∂x+y→1sy∂∂y+z→1sz∂∂z,where si=ki+σpi/αi+jωμ, ki and αi are positive and less than 1 and σpi is the conductivity of the CFS-PML layers (i=x,y,z).The difference scheme of Maxwell’s curl equations in the stretched coordinate space was developed by Roden and Gedney similar to the application of the CFS-PML [17]. Thus, in this paper, we mainly deduce the discrete form of Maxwell’s divergence equation (equation (3)).Inserting (4) into (3) leads to
(5)1sx∂Hx∂x+1sy∂Hy∂y+1sz∂Hz∂z=0.Equation (5) is next transformed into the time domain to obtain the renewal equation of Hz, after which the expression of Si (i = x, y, z) is inserted into (5), which can consequently be expressed as
(6)1kx∂Hx∂x+1ky∂Hy∂y+1kz∂Hz∂z+ψhzx+ψhzy+ψhzz=0,where ψhzi is the auxiliary expression of a convolution item and is implemented as
(7)ψhzi=−σpiεki2exp−σpiεki+αiεt∗∂Hi∂i.The governing equations of other fields in the stretched coordinate can be similarly derived.


## 2.2. The Discrete Form of Maxwell’s Divergence Equation in the Stretched Coordinate Space

During ATEM modeling, the time-step of an FDTD method usually increases gradually according to the attenuation characteristics of TEM responses (i.e., the TEM field exhibits relatively sharp variations at earlier times and gradually becomes smooth thereafter). Here, the time-step is advanced using a modified Du Fort-Frankel method. The traditional recursive convolution formula is deduced based on a uniform time-step; therefore, we need to modify the derivation procedure based on inhomogeneous time-steps which is shown in the appendix. In the stretched coordinate space,ψhzx at tn+1/2 can be expressed as
(8)ψhzxn+12=−σpxεkx2exp−σpxεkx+αxεtn+1/2∗∂Hxn+1/2∂x=ξhxtn+1/2∗∂Hxn+1/2∂x.To obtain the discrete form of the convolution formulation, we assume thatΔtn=tn+1/2−tn−1/2 (n=0,1,2,…,n), and thus, equation (8) can be rewritten as
(9)ψhzxn+12=ξhxtn+1/2∗∂Hxn+1/2∂x=∫0tn+1/2ξhxτ∂Hxtn+1/2−τ∂xdτ,where ξhxt=−Axe−Bxt, Ai=σpi/εki2, Bi=σpi/εki+αi/ε (i=x,y,z).Thenψhzxn+1/2 can be performed recursively:
(10)ψhzxn+12=12AxBxe−BΔtn−AxBx∂Hxn+1/2∂x+∂Hxn−1/2∂x+e−BxΔtnψhzxn−12.We can obtainψhzzn+1/2 and ψhzxn+1/2 in the same manner:
(11)ψhzzn+12=12AzBze−BzΔtn−AzBz∂Hzn+1/2∂z+∂Hzn−1/2∂z+e−BzΔtnψhzzn−12,(12)ψhzxn+12=12AxBxe−BxΔtn−AxBx∂Hxn+1/2∂x+∂Hxn−1/2∂x+e−BxΔtnψhzxn−12.As shown in (6) and (11), ψhzz contains Hz at the current time in the discretization of Maxwell’s divergence equation to obtain the iterative formula of Hz at the time tn+1/2. Therefore, the iterative formula of Hz requires further derivation. Equation (6) can be discretized in both space and time according to the staggered Yee scheme as
(13)1kxHxn+1/2i+1,j+1/2,k+1/2−Hxn+1/2i,j+1/2,k+1/2Δx+1kyHyn+1/2i+1/2,j+1,k+1/2−Hyn+1/2i+1/2,j,k+1/2Δy+1kzHzn−1+1/2i+1/2,j+1/2,k+1−Hzn−1+1/2i+1/2,j+1/2,kΔz+ψhzxn+1/2+ψhzyn+1/2+ψhzzn+1/2=0.Inserting (10)–(12) into (13), we obtain
(14)Hzn+1/2i+12,j+12,k=Hzn+1/2i+12,j+12,k+1+CzDzΔzi+1/2,j+1/2,kHzn−1+1/2i+12,j+12,k+1−Hzn−1+1/2i+12,j+12,k+DzkyΔyi+1/2,j,k+1/2Hyn+1/2i+12,j+1,k+12−Hyn+1/2i+12,j,k+12+DzkxΔxi,j+1/2,k+1/2Hxn+1/2i+1,j+12,k+12−Hxn+1/2i,j+12,k+12+Dzψhzxn+1/2+Dzψhzyn+1/2+e−BzΔtnDzψhzzn−1/2,where Cz=1/2Az/Bze−BzΔtn−Az/Bz, Dz=kz/1+kzCzΔzi+1/2,j+1/2,k.Note thatHz at tn+1/2 is related to the value of ψhzz of the previous time-step, and thus, we must refresh ψhzz at the end of the calculation for the next iteration.


## 3. Numerical Tests

To investigate the feasibility of the proposed method, we first compare the solutions of our method with the numerical solutions calculated by the integral method in homogeneous half-space models [8, 34, 35]. Then, the proposed method is applied to several models with 3D conductors. With the exception of the slanted model, every model has 101 × 101 × 50 grids [36]. The grids are nonuniform with a minimum spacing of 10 m and a maximum spacing of 120 m. The spacing of the CFS-PML boundary is 120 m, and the boundary is not contained within any of the abovementioned grids. A transmitting coil is located at the center of the earth model with a height of 120 m. The magnetic moment is 4π × 10−7A·m2, and the transmitter current is 0.7 × 107A. The receiving coil is situated at a height of 60 m, and it is 130 m away from the transmitting coil in the x-direction.


### 3.1. Validation with Homogeneous Half-Space Models

The electromagnetic responses of the homogeneous half-space models are calculated, and the FDTD solutions are compared with the numerical solutions to verify the accuracy of our method. To demonstrate the absorbing effect of the proposed method in different host media, the conductivities are set as 0.1 S/m, 0.01 S/m, and 0.005 S/m; the diagrammatic drawing of the homogeneous half-space model is shown in Figure1.Figure 1
Diagrammatic drawing of homogeneous half-space models.The FDTD-based solutions (calculated using the DBC and the CFS-PML with 8 layers) are compared with the numerical solutions, as shown in Figure2. In the models with different conductivities, we can discern that these solutions are in good agreement at earlier times because the diffusion fields have not reached the boundaries yet. With an increase in the time; however, the effectiveness of the CFS-PML gradually becomes clearer, especially in the models with a smaller conductivity (as presented in Figure 2 b~c). The average relative errors of the FDTD solutions using the CFS-PML and DBC are shown in Figure 2(d). In the model with a large conductivity (0.1 S/m) and a slower propagation speed, these boundary conditions both demonstrate a better absorption of electromagnetic waves with average relative errors of 1.61% and 0.79% (corresponding to the DBC and CFS-PML, respectively). When the model conductivity is 0.01 S/m, the average relative errors of the DBC and CFS-PML are 21.47% and 2.63%, respectively, effectively revealing the better ability of the CFS-PML to absorb electromagnetic waves. In addition, when the ground conductivity is 0.005 S/m, the strong electromagnetic wave reflected from the DBC leads to an enormous deviation in the numerical solution, whereas the CFS-PML still exhibits good absorption; in this case, the average relative errors of the DBC and CFS-PML are 42.9% and 4.17%, respectively.Figure 2
Comparison of the induced voltage at the center of the receiving coil at 10 ms when the conductivity is (a) 0.1 S/m, (b) 0.01 S/m, and (c) 0.005 S/m; (d) shows a comparison of the average relative errors obtained at these conductivities.
(a)
(b)
(c)
(d)Except for the decay curves at the center of the receiving coil, we also investigate the electromagnetic reflections at the boundaries. Therefore, time slices of the electromagnetic responses from a homogeneous half-space model with a conductivity of 0.01 S/m are presented in Figure3. As illustrated, the electromagnetic responses calculated by the DBC and CFS-PML boundaries are compared at three separate moments. At earlier times, the electromagnetic responses of the DBC and CFS-PML boundaries display the same variations; meanwhile, at later times, serious distortions are observed at both 4.5 ms and 9 ms when the DBC is used. In contrast, no distortions are evident when the CFS-PML is adopted regardless of the point in time, indicating that the CFS-PML boundary condition demonstrates a better ability to absorb electromagnetic waves.Figure 3
Time slices of the electromagnetic responses in the air using the DBC at (a) 1 ms, (b) 4.5 ms, and (c) 10 ms in addition to those using the CFS-PML at (d) 1 ms, (e) 4.5 ms, and (f) 10 ms. The ground conductivity is 0.01 S/m.
(a)
(b)
(c)
(d)
(e)
(f)To study the reflection errors, the time-dependent field at the center of the receiving coil is recorded when the conductivity of the host medium is 0.01 S/m. The error relative to the numerical solution (in dB) is computed using:(15)ErrordB=20log10Vzt−VzreftVzrefmax,where Vzt represents the voltage received by the receiving coil, Vzref(t) represents the numerical solutions, and Vzrefmax represents the maximum of Vzref(t). The relative errors (in dB) computed via (13) are recorded in Figure 4, from which we can see that the relative errors are similar at earlier times. Moreover, compared with the DBC, the CFS-PML can reduce the relative error by nearly 60 dB at later times.Figure 4
Relative errors at the center of the receiving coil in dB using the CFS-PML and DBC at 10 ms when the conductivity is 0.01 S/m.


### 3.2. Validation with Anomalous Models

To test the effectiveness of the proposed method, we design three anomalous models. In each model, the ground conductivity is 0.01 S/m, and the thickness of the CFS-PML is 10 cells. The first model is shown in Figure5; two anomalous bodies with different sizes, depths, and conductivities are selected. An analytical solution exists only in the homogeneous half-space. Therefore, to study the reflection errors due to the boundaries, the mesh is extended by an additional 40 cells in each direction to obtain the reference solution, and the CFS-PML is used to minimize the reflections. To ensure that the electromagnetic waves cannot reach the boundaries during the observation time (10 ms), the diffusion depth is set as
(16)d=2tσμ0,where σ is the conductivity, t is the observation time, and d is the diffusion depth [37].Figure 5
Diagrammatic drawing of model 1 in (a) x-z plane and (b) x-y plane.
(a)
(b)Time slices of the anomaly responses are shown in Figure6. Both at earlier times (3 ms) and at later times (10 ms); reflected electromagnetic waves exist in the model when the DBC is adopted. Moreover, especially at 10 ms, the locations of anomalous bodies cannot be identified accurately. In contrast, the reflections are absorbed well in the model with the CFS-PML, and the results are in agreement with the reference solution. A comparison of ∂B/∂t at the center of the lattice along the x-direction between 3 ms and 10 ms is shown in Figure 7. In earlier times, the DBC solution is in good agreement with the reference solution at the center of the model; however, such good agreement is not observed near the boundary. Furthermore, at 10 ms, the DBC solution severely deviates from the reference solution along the whole boundary. Meanwhile, the responses in the model using the CFS-PML are consistent with the reference solution regardless of the time. The relative errors are shown in Figures 7(c) and 6(d). The maximum relative error is less than 3% at earlier times and no more than 4.6% at later times when the CFS-PML is adopted. However, the results obtained with the DBC are consistent with the reference solution only at the center of the model area at earlier times.Figure 6
Comparison of the responses in model 1 in the plane of receiving coil at 3 ms (a) (b) (c) and 10 ms (d) (e) (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 7
Comparison of the induced voltage at the line cross the center of receiving coil along x-direction at (a) 3 ms and (b) 10 ms and the relative errors using (c) the CFS-PML and (d) DBC of model 1.
(a)
(b)
(c)
(d)For further verification, an anomalous model is designed and presented in Figure8. In this model, the anomalous body adjoins the boundary. The dimensions of the anomalous body are 480 m × 120 m × 180 m, and the depth of the body is 190 m. Time slices of the anomaly responses are shown in Figure 9. Intense reflections of electromagnetic waves appear when the DBC is utilized, and the information of the anomalous body cannot be retrieved accurately. In contrast, through the implementation of the CFS-PML boundary condition, the reflections are effectively suppressed, and the results are approximately equivalent to the reference solution. Consequently, the information of the anomalous body is presented clearly. The responses at the boundary at the center of the model along the x-direction are shown in Figure 10. The DBC solution evidently coincides with the reference solution at the center of the boundary, although the DBC solution is different near the boundary at earlier times and performs progressively worse with increasing time. Meanwhile, the responses obtained using the CFS-PML are still in good agreement. The maximum relative error of the results acquired using the CFS-PML is less than 5% at later times. Furthermore, the DBC solutions exhibit serious reflections, especially near the anomalous body. From Figures 8 and 9, we can conclude that the CFS-PML is very stable and can be established very close to the anomalous body.Figure 8
Diagrammatic drawing of model 2 in (a) x-z plane and (b) x-y plane.
(a)
(b)Figure 9
Comparison of the responses in model 2 in the plane of receiving coil at 3 ms (a) (b) (c) and 10 ms (d) (e) (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 10
Comparison of the induced voltage at the line cross the center of receiving coil along x-direction at (a) 3 ms and (b) 10 ms and the relative errors using (c) the CFS-PML and (d) DBC of model 2.
(a)
(b)
(c)
(d)Finally, we design an oblique model. The anomalous body is oblique with an angle of 45°, and it is established at the center of the x-y plane, as shown in Figure11. The dimensions of the anomalous body are 400 m × 610 m × 90 m, and the depth is 350 m. Time slices of the anomaly responses are shown in Figure 12. Reflections of electromagnetic waves evidently occur along the boundary at earlier times when the DBC is adopted, while the responses at the center are basically equivalent to those in the reference solution. At 10 ms, the reflections at the DBC boundaries are so intense that the diffusion is completely distorted. In contrast, the responses with the CFS-PML boundary condition are highly consistent with the reference solution, even at later times. The received ∂B/∂t at the boundary across the center of the receiving coil in the x-direction is shown in Figure 13. It is clear that the DBC solution coincides with the reference solution at the center of the model area at earlier times and seriously deviates from the reference solution at later times. However, the CFS-PML can effectively absorb the reflections from the boundary, even at later times, and the maximum relative error is less than 5%.Figure 11
Diagrammatic drawing of model 3 in (a) x-z plane and (b) x-y plane.
(a)
(b)Figure 12
Comparison of the responses in model 3 in the plane of receiving coil at 3 ms (a) (b) (c) and 10 ms (d) (e) (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 13
Comparison of the induced voltage at the line cross the center of receiving coil along x-direction at (a) 3 ms and (b) 10 ms and the relative errors using (c) the CFS-PML and (d) DBC of model 3.
(a)
(b)
(c)
(d)


## 3.1. Validation with Homogeneous Half-Space Models

The electromagnetic responses of the homogeneous half-space models are calculated, and the FDTD solutions are compared with the numerical solutions to verify the accuracy of our method. To demonstrate the absorbing effect of the proposed method in different host media, the conductivities are set as 0.1 S/m, 0.01 S/m, and 0.005 S/m; the diagrammatic drawing of the homogeneous half-space model is shown in Figure1.Figure 1
Diagrammatic drawing of homogeneous half-space models.The FDTD-based solutions (calculated using the DBC and the CFS-PML with 8 layers) are compared with the numerical solutions, as shown in Figure2. In the models with different conductivities, we can discern that these solutions are in good agreement at earlier times because the diffusion fields have not reached the boundaries yet. With an increase in the time; however, the effectiveness of the CFS-PML gradually becomes clearer, especially in the models with a smaller conductivity (as presented in Figure 2 b~c). The average relative errors of the FDTD solutions using the CFS-PML and DBC are shown in Figure 2(d). In the model with a large conductivity (0.1 S/m) and a slower propagation speed, these boundary conditions both demonstrate a better absorption of electromagnetic waves with average relative errors of 1.61% and 0.79% (corresponding to the DBC and CFS-PML, respectively). When the model conductivity is 0.01 S/m, the average relative errors of the DBC and CFS-PML are 21.47% and 2.63%, respectively, effectively revealing the better ability of the CFS-PML to absorb electromagnetic waves. In addition, when the ground conductivity is 0.005 S/m, the strong electromagnetic wave reflected from the DBC leads to an enormous deviation in the numerical solution, whereas the CFS-PML still exhibits good absorption; in this case, the average relative errors of the DBC and CFS-PML are 42.9% and 4.17%, respectively.Figure 2
Comparison of the induced voltage at the center of the receiving coil at 10 ms when the conductivity is (a) 0.1 S/m, (b) 0.01 S/m, and (c) 0.005 S/m; (d) shows a comparison of the average relative errors obtained at these conductivities.
(a)
(b)
(c)
(d)Except for the decay curves at the center of the receiving coil, we also investigate the electromagnetic reflections at the boundaries. Therefore, time slices of the electromagnetic responses from a homogeneous half-space model with a conductivity of 0.01 S/m are presented in Figure3. As illustrated, the electromagnetic responses calculated by the DBC and CFS-PML boundaries are compared at three separate moments. At earlier times, the electromagnetic responses of the DBC and CFS-PML boundaries display the same variations; meanwhile, at later times, serious distortions are observed at both 4.5 ms and 9 ms when the DBC is used. In contrast, no distortions are evident when the CFS-PML is adopted regardless of the point in time, indicating that the CFS-PML boundary condition demonstrates a better ability to absorb electromagnetic waves.Figure 3
Time slices of the electromagnetic responses in the air using the DBC at (a) 1 ms, (b) 4.5 ms, and (c) 10 ms in addition to those using the CFS-PML at (d) 1 ms, (e) 4.5 ms, and (f) 10 ms. The ground conductivity is 0.01 S/m.
(a)
(b)
(c)
(d)
(e)
(f)To study the reflection errors, the time-dependent field at the center of the receiving coil is recorded when the conductivity of the host medium is 0.01 S/m. The error relative to the numerical solution (in dB) is computed using:(15)ErrordB=20log10Vzt−VzreftVzrefmax,where Vzt represents the voltage received by the receiving coil, Vzref(t) represents the numerical solutions, and Vzrefmax represents the maximum of Vzref(t). The relative errors (in dB) computed via (13) are recorded in Figure 4, from which we can see that the relative errors are similar at earlier times. Moreover, compared with the DBC, the CFS-PML can reduce the relative error by nearly 60 dB at later times.Figure 4
Relative errors at the center of the receiving coil in dB using the CFS-PML and DBC at 10 ms when the conductivity is 0.01 S/m.


## 3.2. Validation with Anomalous Models

To test the effectiveness of the proposed method, we design three anomalous models. In each model, the ground conductivity is 0.01 S/m, and the thickness of the CFS-PML is 10 cells. The first model is shown in Figure5; two anomalous bodies with different sizes, depths, and conductivities are selected. An analytical solution exists only in the homogeneous half-space. Therefore, to study the reflection errors due to the boundaries, the mesh is extended by an additional 40 cells in each direction to obtain the reference solution, and the CFS-PML is used to minimize the reflections. To ensure that the electromagnetic waves cannot reach the boundaries during the observation time (10 ms), the diffusion depth is set as
(16)d=2tσμ0,where σ is the conductivity, t is the observation time, and d is the diffusion depth [37].Figure 5
Diagrammatic drawing of model 1 in (a) x-z plane and (b) x-y plane.
(a)
(b)Time slices of the anomaly responses are shown in Figure6. Both at earlier times (3 ms) and at later times (10 ms); reflected electromagnetic waves exist in the model when the DBC is adopted. Moreover, especially at 10 ms, the locations of anomalous bodies cannot be identified accurately. In contrast, the reflections are absorbed well in the model with the CFS-PML, and the results are in agreement with the reference solution. A comparison of ∂B/∂t at the center of the lattice along the x-direction between 3 ms and 10 ms is shown in Figure 7. In earlier times, the DBC solution is in good agreement with the reference solution at the center of the model; however, such good agreement is not observed near the boundary. Furthermore, at 10 ms, the DBC solution severely deviates from the reference solution along the whole boundary. Meanwhile, the responses in the model using the CFS-PML are consistent with the reference solution regardless of the time. The relative errors are shown in Figures 7(c) and 6(d). The maximum relative error is less than 3% at earlier times and no more than 4.6% at later times when the CFS-PML is adopted. However, the results obtained with the DBC are consistent with the reference solution only at the center of the model area at earlier times.Figure 6
Comparison of the responses in model 1 in the plane of receiving coil at 3 ms (a) (b) (c) and 10 ms (d) (e) (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 7
Comparison of the induced voltage at the line cross the center of receiving coil along x-direction at (a) 3 ms and (b) 10 ms and the relative errors using (c) the CFS-PML and (d) DBC of model 1.
(a)
(b)
(c)
(d)For further verification, an anomalous model is designed and presented in Figure8. In this model, the anomalous body adjoins the boundary. The dimensions of the anomalous body are 480 m × 120 m × 180 m, and the depth of the body is 190 m. Time slices of the anomaly responses are shown in Figure 9. Intense reflections of electromagnetic waves appear when the DBC is utilized, and the information of the anomalous body cannot be retrieved accurately. In contrast, through the implementation of the CFS-PML boundary condition, the reflections are effectively suppressed, and the results are approximately equivalent to the reference solution. Consequently, the information of the anomalous body is presented clearly. The responses at the boundary at the center of the model along the x-direction are shown in Figure 10. The DBC solution evidently coincides with the reference solution at the center of the boundary, although the DBC solution is different near the boundary at earlier times and performs progressively worse with increasing time. Meanwhile, the responses obtained using the CFS-PML are still in good agreement. The maximum relative error of the results acquired using the CFS-PML is less than 5% at later times. Furthermore, the DBC solutions exhibit serious reflections, especially near the anomalous body. From Figures 8 and 9, we can conclude that the CFS-PML is very stable and can be established very close to the anomalous body.Figure 8
Diagrammatic drawing of model 2 in (a) x-z plane and (b) x-y plane.
(a)
(b)Figure 9
Comparison of the responses in model 2 in the plane of receiving coil at 3 ms (a) (b) (c) and 10 ms (d) (e) (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 10
Comparison of the induced voltage at the line cross the center of receiving coil along x-direction at (a) 3 ms and (b) 10 ms and the relative errors using (c) the CFS-PML and (d) DBC of model 2.
(a)
(b)
(c)
(d)Finally, we design an oblique model. The anomalous body is oblique with an angle of 45°, and it is established at the center of the x-y plane, as shown in Figure11. The dimensions of the anomalous body are 400 m × 610 m × 90 m, and the depth is 350 m. Time slices of the anomaly responses are shown in Figure 12. Reflections of electromagnetic waves evidently occur along the boundary at earlier times when the DBC is adopted, while the responses at the center are basically equivalent to those in the reference solution. At 10 ms, the reflections at the DBC boundaries are so intense that the diffusion is completely distorted. In contrast, the responses with the CFS-PML boundary condition are highly consistent with the reference solution, even at later times. The received ∂B/∂t at the boundary across the center of the receiving coil in the x-direction is shown in Figure 13. It is clear that the DBC solution coincides with the reference solution at the center of the model area at earlier times and seriously deviates from the reference solution at later times. However, the CFS-PML can effectively absorb the reflections from the boundary, even at later times, and the maximum relative error is less than 5%.Figure 11
Diagrammatic drawing of model 3 in (a) x-z plane and (b) x-y plane.
(a)
(b)Figure 12
Comparison of the responses in model 3 in the plane of receiving coil at 3 ms (a) (b) (c) and 10 ms (d) (e) (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 13
Comparison of the induced voltage at the line cross the center of receiving coil along x-direction at (a) 3 ms and (b) 10 ms and the relative errors using (c) the CFS-PML and (d) DBC of model 3.
(a)
(b)
(c)
(d)


## 4. Conclusions and Perspectives

We applied the CFS-PML boundary condition to the modeling of ATEM data in a source-free medium and introduced the source-free Maxwell’s equations as the governing equations. Furthermore, we provided discrete forms of the governing equations in a CFS-PML formulation and deduced the recursive convolution formula with a nonuniform time-step in the stretched coordinate space. The validated results of a homogeneous half-space model demonstrate that the proposed method can absorb incident wave better, and it can reduce the reflections at later times. The CFS-PML is highly absorptive and can considerably reduce the relative error with a 60 dB improvement compared with the traditional boundary condition (i.e., the DBC) at later times.Finally, we employed the CFS-PML in airborne anomalous models and compared the electromagnetic responses with the reference solution. The results show that the CFS-PML demonstrates better absorptive capabilities and can accurately retrieve the information of the anomalous body regardless of whether the anomalous body is set obliquely or near the boundary of the model area. Moreover, the maximum relative error of the improved method is less than 5% at later times (10 ms). The proposed method also displays a better recognition precision in the presence of multiple or complex 3D conductive bodies, thereby verifying the effectiveness of this method for ATEM modeling. Consequently, the proposed method can provide higher-precision simulation results for deep geological exploration.



---

*Source: 1846427-2018-09-12.xml*