*4.1. Forward Modeling Results of 3D Chargeable Body*

%QW

In order to analyze the influence of system parameters on the IP effect, we build a three-dimensional model with the consideration of a chargeable body, as shown in Figure 8. The ground conductivity was <sup>σ</sup>∞<sup>=</sup> 0.001 S/m, the air layer conductivity was <sup>σ</sup>air<sup>=</sup> <sup>10</sup>−<sup>6</sup> S/m, Tx was the long wire source, the length was LTx= 400 m, Rx represents the receiving point, the height from the ground was h = 25 m, and the offset distance was r = 200 m. The dimensions of the chargeable body were <sup>400</sup> <sup>×</sup> <sup>400</sup> <sup>×</sup> 300 m3, which was 100 m below the surface. The model parameters of the chargeable body were selected as follows: chargeability η= 0.5, frequency dependence c = 0.5, and characteristic time constant τ= 0.01 s. The conductivities at infinite frequency were σ∞= 0.01 S/m, σ<sup>∞</sup> = 0.5 S/m, and σ<sup>∞</sup> = 0.1 S/m. The calculation results are shown in Figure 9. The snapshots of the induced current system in the fictitious wave field are shown in Figure 10a,b.

**Figure 8.** Three-dimensional (3D) model with IP effect.

K P Ș IJ V F

ı 6P ı 6P ı 6P

7LPHV

**Figure 9.** Response curves under different conductivities.

**Figure 10.** The snapshots of the induced current system in the fictitious wave field (**a**) early time; (**b**) late time.

It can be seen from Figure 9 that, with the increase in conductivity at the infinite frequency, the attenuation of the early induction field gradually decreased, the attenuation of the late IP response gradually quickened, and the occurrence time of the inverse sign gradually delayed. This is because the influence of the induction field of the low-resistance body was greater than that of the high-resistance body. Therefore, the second field decayed slowly, and the negative response appeared later. Thus, the polarization phenomenon was more easily observed on the high-resistance body. The chargeable body position can be clearly observed from Figure 10, and the method of wave field transformation could effectively identify the polarizable body. Since we set the air layer as a finite high-resistance body, the conductivity was not zero, and the propagation of induced current in the air layer was not zero. In the fictitious wave field, due to the influence of the chargeable body, the propagation mode of the electromagnetic field changed, and it no longer propagated in the form of symmetry.

After setting σ<sup>∞</sup> = 0.1 S/m, as well as the ground conductivity σground = 0.001 S/m, σground = 0.01 S/m, and σground = 0.1 S/m, which are associated with high-, medium-, and low-resistance areas, the IP effects were calculated as shown in Figure 11.

 With the increase in ground conductivity, the amplitude of the early induction field and the late polarization field increased. In the high- and middle-resistance regions, the negative responses occurred at 1.7 ms and 3.3 ms. On the other hand, in the low-resistance region, the negative response could not be observed. When the negative value reached the minimum, the curve decayed approximately exponentially. The GREATEM had better resolution in the low-resistance region, but it is more suitable for the high-resistance region when measuring the IP effect.

ıJURXQG 6P

K PȘ IJ VF ı 6P

%QW

7LPHV

 

]P

**Figure 11.** IP effects under different background conductivities.

## *4.2. The Influence of System Parameters on IP Response*

ORJ

[P

E

 

 

]P

ORJ

[P

D

In the field measurements of GREATEM, due to the influence of airflow and environment, the flight height was not a constant. Therefore, it was necessary to analyze the influence of flight altitude on the IP effect. We take the 3D model (Figure 8) as an example. The air layer conductivity was set to σair= 10−<sup>6</sup> S/m and the earth conductivity was set to 0.001 S/m, while the size of the chargeable body was 400 × 400 × 300 m. The model parameters of the chargeable body were as follows: chargeability η= 0.5, frequency dependence c = 0.5, characteristic time constant τ= 0.01 s, conductivity at the infinite frequency σ∞= 0.1 S/m, offset distance r = 200 m, and flight altitude h = 25 m, h = 50 m, h = 75 m, and h = 100 m. The induced polarization (IP) effect at different flight altitudes was calculated as depicted in Figure 12. In order to understand the influence of the parameters of GREATEM on the IP response, other parameters remained unchanged, while the flight height was set to h = 25 m, and the emission source length was LTx= 50 m, LTx= 100 m, and LTx= 200 m. The responses of different emission source lengths are shown in Figure 13. Taking LTx= 400 m and the offset distance r = 100 m, r = 200 m, and r = 400 m, the responses of different offsets are shown in Figure 14.

Based on Figure 12, flight altitude mainly affected the amplitude of the response. The effect on the time of sign reversal was not obvious, and the sign reversal phenomenon was mainly concentrated in the range 1.5–2 ms. With the increase in flight altitude, the received electromagnetic response decreased gradually. It can be seen from Figure 13 that the length of the emission source mainly affected the amplitude of the response. The length of emission source increased gradually, and the amplitude of response increased gradually, which also shows that a longer emission source resulted in a higher resolution of the polarizable body. In Figure 14, with the decrease in offset distance, the time of the sign reversal appeared earlier, and the polarization response became easier to observe at a short offset distance.

Next, we set the air layer conductivity to σ = 10−<sup>6</sup> S/m and earth conductivity to 0.001 S/m. The size of the chargeable body was 400 <sup>×</sup> <sup>400</sup> <sup>×</sup> 100 m3. The model parameters of the chargeable body were chosen as η= 0.5, c = 0.5, τ= 0.01 s, σ∞= 0.1 S/m, offset distance r = 200 m, emission source length LTx= 400 m, and flight altitude h = 25 m. The chargeable body was buried 25 m and 50 m below the surface. The response curves of the chargeable body at different depths are shown in Figure 15. Moreover, we changed the size of the chargeable body to 300 <sup>×</sup> <sup>300</sup> <sup>×</sup> 200 m<sup>3</sup> and the chargeable body depth to 100 m, while keeping other parameters the same, and setting the distance from the left boundary of the chargeable body to the emission source as 0 m, 50 m, and 100 m. The response curves at different source positions are shown in Figure 16.

**Figure 12.** IP effects at different flight altitudes.

**Figure 13.** The response curves under different emission source lengths.

**Figure 14.** The response curves under different offset distances.

From Figures 15 and 16, it can be seen that all the amplitudes decreased, while the depth of chargeable body increased. However, the time of the reversal sign did not change. In the early stage, as the distance between the emitter source and the chargeable body decreased, the amplitude of the induced field increased, the attenuation of the curve was delayed. In the later stage, the time of sign reversal appeared later when the distance between the emitter and the chargeable body decreased.

7LPHV

'HSWKP 'HSWKP

Ș IJ V F ı 6P

%QW

7LPHV

U P U P U P

K P Ș IJ V F ı 6P

%QW

**Figure 15.** The influence of chargeable body depth on response curve.

**Figure 16.** The influence of different emitter locations on IP effect.

#### **5. Conclusions**

Three-dimensional numerical simulations are the basis of data interpretation and imaging. In this paper, a three-dimensional numerical simulation of underground polarized media was realized, which achieved simulation results closer to the actual mineral resources. In the actual measurement process with GREATEM, the geodetic model could be established using this method, and the optimal configuration structure of the equipment was obtained using the numerical simulation method. This not only provides the basis for the exploration of underground mineral resources, but also provides theoretical guidance for the design of detection instruments. Wave field transformation was proposed to calculate the response of the GREATEM system using the Cole–Cole model. The electromagnetic wave was absorbed by the boundary condition of CFS-PML in the fictitious wave field. Compared with the analytical solution and the transformed FD method, the relative error was less than 10%, satisfying the standard calculation requirements. Compared with the existing method, the calculation time was reduced five-fold. According to the simulations, the flight altitude, electrical source length, electrical source position, and the offset distance all affected the electromagnetic response monotonically. In the homogeneous half-space model, when the characteristic time constant (τ > 1 s) and conductivity at the infinite frequency (σ<sup>∞</sup> > 1 S/m) were too large, the negative response was likely overwhelmed by the noise; thus, the negative response could not be observed. Therefore, the occurrence time of the negative response could be enhanced by reducing the offset distance. Low-altitude flights can facilitate the recognition ability of the chargeable body, but they are often limited by geological conditions. The response of the GREATEM system had good recognition ability in the low-resistance area; however, the

chargeable body could hardly be detected. Therefore, it is necessary to improve the detection accuracy of the GREATEM in low-resistance areas.

This paper realized a three-dimensional numerical simulation of polarized media based on the Cole–Cole model, and improved the calculation efficiency. However, in the actual measurement process, the underground medium is diverse. It is not comprehensive to only use the Cole–Cole model to describe the characteristics of polarized medium; the existence of the GEMTIP model in particular is a major problem in the field of geophysics. In future work, we need to consider the GEMTIP model in the numerical simulation, so as to improve the universality of this method.

**Author Contributions:** Conceptualization, X.M. and Y.J.; methodology, X.M.; software, X.M.; validation, X.M., Y.J. and G.L.; formal analysis, Y.W.; investigation, Y.W.; resources, W.H.; data curation, X.M.; writing—original draft preparation, X.M.; writing—review and editing, W.H.; supervision, G.L.; funding acquisition, Y.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (Grant Number 41674109).

**Acknowledgments:** We thank the editors and reviewers for comment. We thank the National Natural Science Foundation of China (Grant Number 41674109) for support.

**Conflicts of Interest:** The authors declare no conflicts of interest.
