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

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 sourcefree 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.


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][2][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][6][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][13][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 frequencyshifted 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-sweepuniform 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 sourcefree 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.

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 of H z , 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.

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 B z in the stretched coordinate space [28][29][30][31][32], while the other electromagnetic fields are advanced using Maxwell's curl equations ((1) and ( 2)).
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]: International Journal of Antennas and Propagation where s i = k i + σ pi / α i + jωμ , k i 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 Equation ( 5) is next transformed into the time domain to obtain the renewal equation of H z , after which the expression of S i (i = x, y, z) is inserted into (5), which can consequently be expressed as where ψ hzi is the auxiliary expression of a convolution item and is implemented as The governing equations of other fields in the stretched coordinate can be similarly derived.

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 t n+1/2 can be expressed as To obtain the discrete form of the convolution formulation, we assume that Δt n = t n+1/2 − t n−1/2 (n = 0, 1, 2, … , n), and thus, equation ( 8) can be rewritten as where Then ψ hzx n + 1/2 can be performed recursively: We can obtain ψ hzz n + 1/2 and ψ hzx n + 1/2 in the same manner: As shown in ( 6) and ( 11), ψ hzz contains H z at the current time in the discretization of Maxwell's divergence equation to obtain the iterative formula of H z at the time t n+1/2 .Therefore, the iterative formula of H z requires further derivation.Equation ( 6) can be discretized in both space and time according to the staggered Yee scheme as 3 International Journal of Antennas and Propagation Inserting ( 10)-( 12) into ( 13), we obtain where Note that H z at t n+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.

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 −7 A•m 2 , and the transmitter current is 0.7 × 10 7 A. 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.

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 Figure 1.
The FDTD-based solutions (calculated using the DBC and the CFS-PML with 8 layers) are compared with the numerical solutions, as shown in Figure 2. 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.
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 Figure 3 4 International Journal of Antennas and Propagation adopted regardless of the point in time, indicating that the CFS-PML boundary condition demonstrates a better ability to absorb electromagnetic waves.
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: where V z t represents the voltage received by the receiving coil, V zref (t) represents the numerical solutions, and V zrefmax represents the maximum of V zref (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.

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 Figure 5; 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 where σ is the conductivity, t is the observation time, and d is the diffusion depth [37].
Time slices of the anomaly responses are shown in Figure 6.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.
For further verification, an anomalous model is designed and presented in Figure 8 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   International Journal of Antennas and Propagation 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%.

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 timestep 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.

Appendix
To describe H x accurately, the integral cannot be divided as usual into [t 0 , t 1 ], [t 1 , t 2 ], … , [t n , t n+1/2 ], because the timesteps are inhomogeneous.Here, we propose a new subdivision, as shown in     International Journal of Antennas and Propagation The further calculation of (17) leads to where ξ hx t = −A x e −B x t , A i = σ pi /εk 2 i , B i = σ pi /εk i + α i /ε (i = x, y, z).

12
International Journal of Antennas and Propagation

Figure 2 :
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.

Figure 3 :Figure 4 :
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.

Figure 5 :
Figure 5: Diagrammatic drawing of model 1 in (a) x-z plane and (b) x-y plane.

Figure 6 :Figure 7 :Figure 8 :
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).

Figure 9 :
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).

Figure 10 :Figure 11 :
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.

From ( 18 ) 1 ∂H+Figure 12 :
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).

Figure 13 :
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.