Application of CSAMT Numerical Simulation in the Exploration and Electrical Structure Interpretation of Jinding Lead-Zinc Deposit

1. Introduction

The Jinding lead-zinc deposit, located in Lanping County, Nujiang Prefecture, Yunnan Province, is the largest lead-zinc deposit in China. Due to the complicated geological and topographical conditions of the "Three-River" metallogenic belt, there are many disputes about the geological genesis and metallogenic model of the Jinding lead-zinc deposit [1,2]. Regarding the geological genesis of the Jinding lead-zinc deposit, there are multiple viewpoints, including hydrothermal filling metasomatic type, stratabound epigenetic deposit, syngenetic sedimentary deposit, and exhalative (hydrothermal) sedimentary deposit [3,4,5,6,7,8]. Among them, the viewpoint of exhalative (hydrothermal) sedimentary deposit dominates. However, some studies have shown that the Jinding lead-zinc deposit is more likely to be an epigenetic filling type stratabound deposit, and this view is supported by geological evidences such as the paleoenvironment of the ore-bearing strata. In order to obtain the deep geological structure of the Jinding lead-zinc mine and better understand the geological environment for mineralization in the deep part of the deposit, it is urgent to carry out a research on a deep exploration method under the circumstance of limited geophysical exploration studies in the Jinding mining area [9]. The Controlled Source Audio frequency Magnetotellurics (CSAMT) method, renowned for its robust anti-interference capabilities, extensive detection depth, and remarkable work efficiency, has emerged as a pivotal tool in deep resource exploration. Its applications have yielded outstanding outcomes in mineral resource exploration, concealed ore body identification, and metallogenic prediction within complex geological terrains. Consequently, CSAMT holds immense significance in safeguarding mineral resource reserves and fostering the sustainable development of the national economy[10,11,12,13,14]. Numerical simulation has universal applicability and can effectively guide people's understanding and knowledge of field exploration methods. Therefore, we can verify the effectiveness of exploration methods through CSAMT numerical simulation results according to the objectives and requirements of field exploration. At the same time, we can also set the field CSAMT acquisition parameters based on the CSAMT forward modeling parameters to ensure data acquisition quality and improve exploration effectiveness.

For general geoelectric models, CSAMT do not have simple numerical solutions and often require reliable numerical solutions to be sought. Traditional methods for solving CSAMT numerical solutions can be summarized into four main algorithms: volume integral method, boundary integral method, finite difference method, and finite element method. The finite element method is adopted by many scholars due to its characteristics such as grid diversity and high calculation accuracy[15,16]. When dealing with CSAMT three-dimensional boundary conditions, the finite element method often adopts Dirichlet boundary conditions, which can lead to issues such as large discrete areas, a high number of nodes, time-consuming calculations, and significant impact on accuracy due to boundary conditions. The infinite element method is a supplement to the finite element method, and it was proposed to overcome the challenges faced by the finite element method in handling infinitely far boundary problems. The core idea of the infinite element method is to replace the traditional truncated boundary with "infinite elements." Through coordinate mapping, the "infinite elements" are infinitely extended in a certain direction, enabling the electromagnetic field to rapidly decay to zero within the infinite elements. Finally, combining the finite element discretization of the internal small area allows for rapid forward modeling[17,18,19]. So, we developed a 3D forward modeling program for CSAMT based on the coupled finite-infinite element methods, and applied it to the forward modeling of real models. The results indicate that compared to the traditional finite element method, the computational area can be reduced by 90%, and the computation speed is increased by about 30%.

Therefore, taking the CSAMT exploration of the V8 multi-functional electrical instrument as an example, we constructed a simplified geophysical model based on the geological overview of the deposit. And, we verified the effectiveness of CSAMT exploration by utilizing a coupled finite-infinite element method of 3D CSAMT forward modeling program. At the same time, we rationally designed the acquisition parameters for field data collection based on the results of the numerical simulation . Subsequently, we conducted noise testing experiments and acquisition time trials. Through meticulous analysis of signal-to-noise ratios and resistivity sounding curves, we identified the optimal transmission intensity and acquisition time parameters that significantly enhance the quality of CSAMT field data acquisition. This knowledge serves as a valuable reference for designing CSAMT fieldwork operations and ensuring the high quality of collected field data. Ultimately, leveraging the outcomes of the continuous medium inversion, we successfully delineated the electrical structure extending to a depth of 1 km beneath the mining area. This comprehensive analysis and interpretation provided a robust geophysical foundation for delving into the deep metallogenic geological backdrop and metallogenic environment of the Jinding lead-zinc deposit.

2. Geological Overview and Geophysical Model Numerical Simulation

2.1. Geological Overview of the Deposit

The Jinding Lead-Zinc Deposit is prominently situated within the densely compressed middle segment of the Sanjiang Fold System, positioned precisely at the northern terminus of the Lanping-Simao Meso-Cenozoic depression. It is strategically positioned between the north-south oriented faults tightly hemmed in by the Mishahe Fault Zone and the Lancang River Fault Zone, and is found to the western flank of the Bihe River Fault Zone, within the ancient Lanping-Yunlong Paleocene graben basin.The deposit is meticulously divided into seven distinct mining sections, including Beichang, Jiayashan, Fengzishan, Paomaping, Xipo, Nanchang, and Baicaoping. The mining area exposes a succession of strata ranging from the older to the younger, comprising the Upper Triassic Sanhedong Formation (T3s) of marine limestones (including brecciated limestones, dolomitic limestones, argillaceous limestones, and limestone mudstones), the Maichuqing Formation (T3m, the core strata of the syncline, with alternating layers of quartz sandstones and carbonaceous mudstones), the Middle Jurassic Huakaizuo Formation (J2h, interbedded layers of argillaceous siltstones, siltstones, and mudstones), the Upper Jurassic Bazhulu Formation (J2b, consisting of mudstones, argillaceous siltstones, and calcareous conglomerates), the Lower Cretaceous Jingxin Formation (K1j, with alternating layers of quartz sandstones, silty mudstones, and intercalated limestone conglomerates), the Upper Cretaceous Nanxin Formation (K2n, composed of argillaceous conglomerates, fine sandstones, siltstones, and mudstones), the Upper Cretaceous Hutousi Formation (K2h, characterized by long quartz sandstones), the Lower Tertiary Paleocene Yunlong Formation (E1y, interbedded layers of muddy gravelly rocks and argillaceous siltstones), the Guolang Formation (E2g, with alternating layers of argillaceous siltstones, silty mudstones, and intercalated fine sandstones), the Eocene Baoxiangsi Formation (Eb, featuring long quartz sandstones, conglomerates, sandstones, mudstones, and quartzites), the Upper Tertiary Pliocene Sanying Formation (N2s, with interbedded layers of conglomerates, sandy conglomerates, sandstones, and mudstones), and the Quaternary Pleistocene (Qp, sand and gravel layers) and Holocene (Qh, comprising gravels, sand grains, and sandy clays).

Figure 1. Simplified Geological Map of the Jinding Mining Area, Yunnan Province.

Figure 1. Simplified Geological Map of the Jinding Mining Area, Yunnan Province.

The primary geological structure of the mining area is a thrust nappe tectonic assemblage, which constitutes a significant component of the large-scale thrust nappe structure that formed during the Paleocene Yunlong period in the Lanping Basin. Klippe structures, or flysch formations, exist in numerous segments throughout the mining area. Furthermore, the thrust-sliding planes have also been incorporated into the tectonic domes, resulting in the formation of dome structures that are characterized by the co-deformation of both thrust tectonics and in-situ systems. In accordance with the preceding findings, it can indeed be concluded that the thrust tectonics predate the formation of the dome structures, with their inception likely commencing as early as the late Yunlong period. In summary, the geological structures of the Jinding mining area are complex, characterized by well-developed folds and faults that exhibit multiple stages of activity. This underscores the dynamic and protracted geological history of the region. The geological structures of the mining area are primarily characterized by the Gaoping-Laomujing syncline, the Bihe River Fault, horizontal thrust faults, transverse faults, oblique faults, and secondary fault structures.

2.2. Geological and Geophysical Model of the Ore Deposit

Based on the existing geological data, physical property measurement and 2D geological model[15] (Figure 2), we have constructed a 3D geophysical model of Jinding lead-zinc deposit, as shown in the Figure 3.

In the Figure 3, the first layer represents the surface covering layer with a thickness of 100m. Its main component is sand and gravel strata, and it has an average resistivity of 200$\Omega \mathrm{m}$. Beneath the covering layer lie domes and concealed deep fault structures with a resistivity of 100$\Omega \mathrm{m}$. On both sides of the fault are high-resistance surrounding rocks with a resistivity of 1000$\Omega \mathrm{m}$. The ore bodies are hosted within the dome structure, primarily consisting of sandstone-type and limestone breccia ore bodies. The sandstone ore body resembles a platy form, with a resistivity of 500$\Omega \mathrm{m}$. Its dimensions are 300m in length, 100m in width, and 300m in height. The limestone breccia ore bodies often appear as inclined or vertical thin layers, also with a resistivity of 500$\Omega \mathrm{m}$, and their dimensions are 100m in length, 300m in width, and 300m in height. The limestone breccia ore bodies often occur as inclined or vertical thin layers with a resistivity of 500$\Omega \mathrm{m}$. Their dimensions are typically 100m in length, 300m in width, and 300m in height.

2.3. The Coupled Finite-infinite Element Method of CSAMT

Numerical simulation research not only enhances people's understanding and knowledge about the exploration methods, but also assists in selecting the correct exploration approach and setting appropriate field acquisition parameters. Therefore, we employed a three-dimensional CSAMT forward modeling program, based on the coupled finite-infinite element method to conduct a numerical simulation on the geophysical model of Jinding lead-zinc mine. The 3D CSAMT forward modeling using the finite-element and infinite-element coupling method is an efficient numerical simulation technique for electromagnetic fields, which combines the advantages of the finite element method in handling complex geological structures and boundary conditions with the characteristics of the infinite element method in simulating infinite far-field attenuation. This method achieves rapid and high-precision forward modeling, reduces the computation domain and the number of nodes, accelerates the computation, and is of great significance to the development of exploration geophysics[20,21].

2.3.1. Fundamental Equation

The electric field generated by a horizontal electric dipole ( time factor$e^{-i\omega t}$, angular frequency$\omega$) in an isotropic conductive medium satisfies the double curl equation [22]:

$$\nabla \times \nabla \times \overrightarrow{E}-i\omega {\mu }_0\left(\sigma -i\omega {\epsilon }_0\right)\overrightarrow{E}=i\omega {\mu }_0{J}_s,$$

(1)

where $\sigma$ is the resistivity, $J_s$ is the current density vector of the electric dipole, and the magnetic permeability ${\mu }_0$ of both the air and the subsurface medium is taken as the vacuum magnetic permeability, while the dielectric constant ${\epsilon }_0$is the vacuum dielectric constant.

The electric dipole source exhibits singularity. The total electric field is decomposed into the sum of the background field${\overrightarrow{E}}^p$ (primary field ) and the induced field ${\overrightarrow{E}}^s$ (secondary field ). The background field ${\overrightarrow{E}}^p$ is directly solved using a one-dimensional geoelectric model[23], while the induced field ${\overrightarrow{E}}^s$ is solved using the finite element method to avoid the calculation of source singularity.

$$\overrightarrow{E}={\overrightarrow{E}}^p+{\overrightarrow{E}}^s,$$

(2)

The background field ${\overrightarrow{E}}^p$also satisfies the double curl equation of the electric field. Equations (1) and (2) can be used to derive the double curl equation based on the secondary field${\overrightarrow{E}}^s$:

$$\nabla \times \nabla \times {\overrightarrow{E}}^s+k^2{\overrightarrow{E}}^s=i\omega {\mu }_0\left(\sigma -{\sigma }^p\right){\overrightarrow{E}}^p,$$

(3)

where$k^2=-i\omega {\mu }_0\left(\sigma -i\omega {\epsilon }_0\right)$. We know that on an electrical interface, the normal component of the electric field is discontinuous, while nodal finite elements require the electromagnetic field to be continuous in the normal direction. Therefore, the obtained finite element solution is often inaccurate and needs to be corrected. In the source region, the electric field solution of nodal finite elements does not satisfy condition $\nabla \cdot \left(\epsilon \overrightarrow{E}\right)=-\nabla \cdot \left(\overrightarrow{J}/i\omega \right)$; in the source-free region, it does not satisfy condition $\nabla \cdot \left(\epsilon \overrightarrow{E}\right)=0$. A divergence correction condition needs to be added to Equation (3) [24].

$$\nabla \times \nabla \times {\overrightarrow{E}}^s+k^2{\overrightarrow{E}}^s+\nabla \cdot {\overrightarrow{E}}^s=i\omega {\mu }_0\left(\sigma -{\sigma }^p\right){\overrightarrow{E}}^p,$$

(4)

2.3.2. Weak Solution Integral Form

Based on the weighted residual finite element method (Jin 2014), establish the residual formula for Equation (4):

$$r=\nabla \times \nabla \times {\overrightarrow{E}}^s+k^2{\overrightarrow{E}}^s+\nabla \cdot {\overrightarrow{E}}^s-i\omega {\mu }_0\left(\sigma -{\sigma }^p\right){\overrightarrow{E}}^p,$$

(5)

Take any test function $\overrightarrow{V}$, on the region $\Omega$:

$${\displaystyle {\iiint }_{\Omega }r}\cdot \overrightarrow{V}d\Omega =0,$$

Then we have:

$${\displaystyle \underset{\Omega }{\int }\nabla \times \nabla \times {\overrightarrow{E}}^s\cdot \overrightarrow{V}d\Omega +}{\displaystyle \underset{\Omega }{\int }{k}^2{\overrightarrow{E}}^s\cdot \overrightarrow{V}d\Omega }+{\displaystyle \underset{\Omega }{\int }\nabla \cdot {\overrightarrow{E}}^s\cdot \overrightarrow{V}d\Omega =}i\omega {\mu }_0{\displaystyle \underset{\Omega }{\int }\left(\sigma -{\sigma }^p\right){\overrightarrow{E}}^p\cdot \overrightarrow{V}d\Omega },$$

(6)

Let S be the outer boundary surface of the region $\Omega$. Instead of using traditional truncated boundary conditions, infinite elements are employed. Within the infinite elements, the electromagnetic field decays to zero, so we have:

$${\overrightarrow{E}}^s=0\mathrm{on}S,$$

(7)

By applying the first vector Green's theorem, Equation (6) can be simplified to:

$${\displaystyle \underset{\Omega }{\int }\left(\nabla \times \overrightarrow{V}\right)\cdot \left(\nabla \times {\overrightarrow{E}}^s\right)d\Omega }+{\displaystyle \underset{\Omega }{\int }{k}^2{\overrightarrow{E}}^s\cdot \overrightarrow{V}d\Omega +{\displaystyle \underset{\Omega }{\int }\nabla \cdot {\overrightarrow{E}}^s\cdot \overrightarrow{V}d\Omega }=}i\omega {\mu }_0{\displaystyle \underset{\Omega }{\int }\left(\sigma -{\sigma }^p\right){\overrightarrow{E}}^p\cdot \overrightarrow{V}d\Omega }$$

(8)

Use finite element method to discretize the internal computation domain. Assuming there are n nodes in the internal computation domain, the j-th test function can be taken as:

$${\displaystyle \underset{\Omega }{\int }\left(\nabla \times {\overrightarrow{V}}_j\right)\cdot \left(\nabla \times {\overrightarrow{E}}^s\right)d\Omega }+{\displaystyle \underset{\Omega }{\int }{k}^2{\overrightarrow{E}}^s\cdot {\overrightarrow{V}}_{j}d\Omega +{\displaystyle \underset{\Omega }{\int }\nabla \cdot {\overrightarrow{E}}^s\cdot {\overrightarrow{V}}_{j}d\Omega }=}i\omega {\mu }_0{\displaystyle \underset{\Omega }{\int }\left(\sigma -{\sigma }^p\right){\overrightarrow{E}}^p\cdot {\overrightarrow{V}}_{j}d\Omega }$$

(9)

2.3.3. The Coupled Finite-Infinite Element Method

The coupled finite-infinite element method divides the entire solution domain into a finite element region and an infinite element region, replacing the traditional outer boundary with the infinite element region. Finite element analysis and infinite element analysis are performed in the two regions respectively, and they are combined through the assembly of the overall stiffness matrix for numerical solution.

Figure 4 is a schematic diagram of the division of finite element and infinite element calculation regions. In the figure, the finite element region is the target region, which includes the field source, target body, measurement points, etc.; the infinite element region extends from the boundary of the finite element region to infinity, serving as the boundary calculation region. Infinite element analysis involves using infinite element mapping and shape functions in a certain direction within this region to map the global coordinates to local coordinates. Its principle is the same as that of finite element analysis.

• The Finite Element Analysis

When performing finite element analysis, the rectangular hexahedron is used for regional discretization, and the element node numbering and coordinates are shown in Figure 2.

In Figure 5, the corresponding relationship between the coordinates of the parent and child elements is as follows:

$$x=x_0+{\displaystyle \frac{a}{2}}\xi \hspace{1em}\hspace{1em}y=y_0+{\displaystyle \frac{b}{2}}\eta \hspace{1em}\hspace{1em}z=z_0+{\displaystyle \frac{c}{2}}\zeta ,$$

(10)

Where$x_0$, $y_0$, $z_0$ represent the midpoint of the child element, and a, b, c represent the three side lengths of the child element. The expression for the shape function of the rectangular hexahedron is as follows:

$$N_i^e={\displaystyle \frac{1}{8}}(1+{\xi }_i\xi )(1+{\eta }_i\eta )(1+{\zeta }_i\zeta ),$$

(11)

In the equation, ${\xi }_i$, ${\eta }_i$, ${\zeta }_i$ represent the coordinates of node i in the child element within the parent element.

2.

The Infinite Element Analysis

When performing infinite element analysis, a three-dimensional eight-node Astley-type infinite element is used[25]. Figure 3 is a schematic diagram of three-dimensional infinite element mapping, where nodes 1, 2, 3, 4, 5, 6, 7, and 8 are the basic elements of the infinite element. P is the coordinate origin, and nodes 1, 2, 3, and 4 are the four nodes of a finite element on the boundary. The distances from point P to nodes 5, 6, 7, and 8 are twice the distances from point P to nodes 1, 2, 3, and 4, respectively. Infinite element analysis involves mapping infinite coordinates to local coordinates in Figure 3(b) through infinite element mapping. In Figure 3(b), the four outermost nodes of the infinite element represent infinity, and their field values are zero.

In Figure 3(b), $\xi$represents the mapping direction of the infinite element. Within plane $\zeta -\eta$, the infinite element and the finite element adopt the same mapping form and shape function. The coordinate mapping for the infinite element is as follows:

$$y={\displaystyle \sum _{i=1}^2{L}_i{y}_i}\hspace{1em}\hspace{1em}z={\displaystyle \sum _{i=1}^2{L}_i{z}_i},$$

(12)

L_i represents the area of the quadrilateral, $L_i={\displaystyle \frac{1}{4}}(1+{\eta }_i\eta )(1+{\zeta }_i\zeta )$Combining the coordinate mapping relationship in direction $\xi$, the coordinate mapping function for the 8-node infinite element can be obtained:

$$x={\displaystyle \sum _{i=1}^8{N}_i{x}_i}\hspace{1em}\hspace{1em}y={\displaystyle \sum _{i=1}^8{N}_i{y}_i}\hspace{1em}\hspace{1em}z={\displaystyle \sum _{i=1}^8{N}_i{z}_i},$$

(12)

Where:

$$\left\{\begin{array}{l}{N}_i=L_i{\displaystyle \frac{-2\xi }{1-\xi }}\hfill \\ N_{i+4}=L_i{\displaystyle \frac{1+\xi }{1-\xi }}\hfill \end{array}\right.i=1,2,3,4,$$

(13)

quadrilateral surface element mapping and linear interpolation of area coordinates

For the shape function of the infinite element, its specific expression is as follows:

$$\left\{\begin{array}{l}{M}_i=L_i{\displaystyle \frac{1-\xi }{2}}{\displaystyle \frac{-(1-\xi )}{2}}\xi \hfill \\ M_{i+4}=L_i{\displaystyle \frac{1-\xi }{2}}(1-{\xi }^2)\hfill \end{array}\hspace{1em}\hspace{1em}\right.i=1,2,3,4,$$

(14)

This shape function is the product of the coefficient A used in the shape function adopted in the Astley Mapped Infinite Element Theory (Astley, 1994) and a second-order Lagrange interpolation polynomial.

The finite element and infinite element analysis are basically consistent, both being eight-node elements. Therefore, in numerical simulations, the infinite elements and finite elements can be perfectly combined to ensure the symmetry of the stiffness matrix, making the solution simple and convenient.

2.3.4. Solving the Equation System

In the three-dimensional numerical simulation of CSAMT, each node has three degrees of freedom in the x, y, and z directions. Therefore, we take:

$$\begin{array}{l}{\overrightarrow{V}}_j=\left({\overrightarrow{V}}_j^x,{\overrightarrow{V}}_j^y,{\overrightarrow{V}}_j^z\right)\hspace{1em}\hspace{1em}\begin{array}{c}{\overrightarrow{V}}_j^x=(N_j,0,0)\end{array}\hspace{1em}\hspace{1em}{\overrightarrow{V}}_j^y=(0,N_j,0)\hspace{1em}\hspace{1em}{\overrightarrow{V}}_j^z=(0,0,N_j)\\ {\overrightarrow{E}}^s=\left({\overrightarrow{E}}_x^s,{\overrightarrow{E}}_y^s,{\overrightarrow{E}}_z^s\right)\hspace{1em}\hspace{1em}{\overrightarrow{E}}_x^s={\displaystyle \sum _{i=1}^n{N}_i}{\overrightarrow{E}}_i^x\hspace{1em}\hspace{1em}{\overrightarrow{E}}_y^s={\displaystyle \sum _{i=n+1}^{2n}{N}_i}{\overrightarrow{E}}_i^y\hspace{1em}\hspace{1em}{\overrightarrow{E}}_z^s={\displaystyle \sum _{i=2n+1}^{3n}{N}_i}{\overrightarrow{E}}_i^z\end{array},$$

(15)

Where, $N_j$ is the shape function of the j-th measurement point, and its specific form is detailed in equations (11) and (14).

Equation (9) can be written as:

$$\begin{array}{l}{\displaystyle \underset{\Omega }{\int }\left(\nabla \times {\overrightarrow{V}}_j^x\right)\cdot \left(\nabla \times {\overrightarrow{E}}^s\right)d\Omega }+{\displaystyle \underset{\Omega }{\int }{k}^2{\overrightarrow{E}}^s\cdot {\overrightarrow{V}}_j^{x}d\Omega +{\displaystyle \underset{\Omega }{\int }\nabla \cdot {\overrightarrow{E}}^s\cdot {\overrightarrow{V}}_j^{x}d\Omega }=}i\omega {\mu }_0{\displaystyle \underset{\Omega }{\int }\left(\sigma -{\sigma }^p\right){\overrightarrow{E}}^p\cdot {\overrightarrow{V}}_j^{x}d\Omega }\\ {\displaystyle \underset{\Omega }{\int }\left(\nabla \times {\overrightarrow{V}}_j^y\right)\cdot \left(\nabla \times {\overrightarrow{E}}^s\right)d\Omega }+{\displaystyle \underset{\Omega }{\int }{k}^2{\overrightarrow{E}}^s\cdot {\overrightarrow{V}}_j^{y}d\Omega +{\displaystyle \underset{\Omega }{\int }\nabla \cdot {\overrightarrow{E}}^s\cdot {\overrightarrow{V}}_j^{y}d\Omega }=}i\omega {\mu }_0{\displaystyle \underset{\Omega }{\int }\left(\sigma -{\sigma }^p\right){\overrightarrow{E}}^p\cdot {\overrightarrow{V}}_j^{y}d\Omega }\\ {\displaystyle \underset{\Omega }{\int }\left(\nabla \times {\overrightarrow{V}}_j^z\right)\cdot \left(\nabla \times {\overrightarrow{E}}^s\right)d\Omega }+{\displaystyle \underset{\Omega }{\int }{k}^2{\overrightarrow{E}}^s\cdot {\overrightarrow{V}}_j^{z}d\Omega +{\displaystyle \underset{\Omega }{\int }\nabla \cdot {\overrightarrow{E}}^s\cdot {\overrightarrow{V}}_j^{z}d\Omega }=}i\omega {\mu }_0{\displaystyle \underset{\Omega }{\int }\left(\sigma -{\sigma }^p\right){\overrightarrow{E}}^p\cdot {\overrightarrow{V}}_j^{z}d\Omega }\end{array},$$

(16)

Using the vector curl formula and the divergence formula, Equation (16) above can be simplified and written in matrix form as follows:

$$Ax=b,$$

(17)

where:

$$\begin{array}{c}A={\displaystyle \underset{\Omega }{\int }\left[\begin{array}{c}{\displaystyle \sum _{i=1}^n\left(\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}+\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+k^2{N}_i{N}_j\right)}\hspace{1em}\hspace{1em}-{\displaystyle \sum _{i=1}^n\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial y}\hspace{1em}\hspace{1em}-{\displaystyle \sum _{i=1}^n\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial z}}}\\ -{\displaystyle \sum _{i=1}^n\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial x}\hspace{1em}\hspace{1em}{\displaystyle \sum _{i=1}^n\left(\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}+\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+k^2{N}_i{N}_j\right)\hspace{1em}\hspace{1em}-{\displaystyle \sum _{i=1}^n\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial z}}}}\\ -{\displaystyle \sum _{i=1}^n\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial x}\hspace{1em}\hspace{1em}-{\displaystyle \sum _{i=1}^n\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial y}\hspace{1em}\hspace{1em}{\displaystyle \sum _{i=1}^n\left(\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}+\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+k^2{N}_i{N}_j\right)}}}\end{array}\right]}dV\\ b=i\omega {\mu }_0{\displaystyle \underset{\Omega }{\int }\left(\sigma -{\sigma }^p\right){\left(\begin{array}{ccc}{N}_j{\overrightarrow{E}}_x^p& N_j{\overrightarrow{E}}_y^p& N_j{\overrightarrow{E}}_z^p\end{array}\right)}^{T}d\Omega }\end{array}$$

(18)

Equation (17) is a large, sparse, symmetric complex coefficient linear equation system. In this system, matrix A represents the overall stiffness matrix of the coupled finite element-infinite element method, which is an $3\times N_x\times N_y\times N_z$order square matrix. Here, $N_x$, $N_y$, $N_z$ enote the total number of nodes in the x, y, and z directions, respectively. Each row in matrix A has a maximum of 81 non-zero elements. x represents the electric field values at each node that need to be solved; b is the right-hand side term, which includes the anomalous body and the primary field. For solving large-scale sparse linear equation systems, this paper adopts Pardiso, an open-source solver with excellent performance and high parallelization. Pardiso uses LU decomposition for direct solving, making it particularly suitable for handling multi-source CSAMT problems.

2.4. The Forward Modeling of the Jinding Lead-Zinc Geophysical Model

We use the coupled finiment-infiniment element method, the forward modeling calculation parameters are as follows:

In the x direction, the subdivision area ranges from -3000 to 15000. Within the survey line area, the grid size is set to 50 meters, resulting in a total of 121 nodes.

In the y direction, the subdivision area ranges from -3000 to 3000. Within the survey line area, the grid size is set to 50 meters, resulting in a total of 63 nodes. Near the field source and the target body, the grid subdivision is relatively dense.

In the z direction, the subdivision area ranges from -3200 to 3000, with a total of 36 nodes. The air region spans from -3200 to 0, and is subdivided into 6 layers. Due to the rapid attenuation of the field in underground media, the grid subdivision size is smaller, the grid size is set to 25 meters.

Frequencies: 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1 Hz, totaling 14 frequencies. The average calculation time for each frequency point is approximately 240 seconds. The numerical simulation results are presented in Figure 4.

In Figure 7, from left to right, the vertical source-receiver distances of the five survey lines are 9600m, 9800m, 10000m, 10200m, and 10400m, respectively. From the apparent resistivity pseudo-section map , it can be observed that the CSAMT 3D forward modeling results exhibit a significant response to the low-resistivity dome structure of the Jinding lead-zinc deposit, and provides a clear reflection to the structure's morphology. However, due to the boundary effects, the boundary between the dome structure and the high-resistivity limestone is relatively blurred. In addition, due to the low-resistivity shielding effect of electromagnetic methods, the relatively high-resistivity lead-zinc ore bodies show little response in the simulation. Based on the numerical simulation results, it can be seen that when exploring for high-resistivity target bodies in low-resistivity surrounding rocks, the CSAMT electromagnetic response is weak, and the numerical effect is not significant. However, when exploring for low-resistivity structures such as faults and fracture zones in high-resistivity surrounding rocks, the effect is very pronounced. Therefore, it is feasible to conduct CSAMT exploration research in the Jinding lead-zinc mining area to identify low-resistivity faults, dome structures, and other features within the ore deposit.

Additionally, the forward modeling parameters used in the 3D CSAMT numerical simulation, such as grid size, vertical source-receiver distance, frequency band range, survey line length, as well as the three-dimensional size and location of the target body, can assist us in designing the field CSAMT acquisition parameters to achieve good exploration results.

3. Acquisition Design for CSAMT

The objective of the CSAMT exploration is to obtain the electrical structure within a depth of 1km below the surface of the mining area, and help us to understand the deep geological background and mineralization environment, at the same time, it can provide geophysical evidence for geological researchers to study the genesis of the Jinding lead-zinc deposit, and if it is possible, which can predict and delineate the mineralization target areas.

3.1. Line Layout

For the CSAMT field data acquisition, a total of 7 survey lines were laid out. The specific locations of the survey lines are shown in Figure 5. From bottom to top in the figure, they are Line 700, Line 1100, Line 1900, Line 2100, Line 2400, Line 3600, and Line 4500. The total length of the survey lines is 13.68 kilometers. Among them, Lines 1900, 2100, and 2400 cross the main ore body and are the key research objects.

3.2. Acquisition Parameter

The parameters for CSAMT field data acquisition mainly consist of two parts: transmitter parameters and receiver parameters. The transmitter operating parameters primarily include the transmitter position (transmitter-receiver distance), dipole length : AB, ground resistance: RAB, transmitting current : I, and transmitting frequency: f. The receiver operating parameters mainly include the measurement mode, receiver MN length, and acquisition time.

For this CSAMT field measurement, the transmitter dipole AB= 2.25km, ground resistance RAB,=20$\Omega \mathrm{m}$. To ensure an exploration depth exceeding 1km, based on the previous forward modeling parameters, the transmitter-receiver distance is set between 10.8km and 14.2km. The transmitting frequency ranges from 7680Hz to 0.9375Hz, with a total of 52 frequency points. The dipole spacing MN =40 m. The CSAET method is adopted for the field CSAMT measurement, with 6 channels of electric field data and 1 channel of magnetic field data collected each time. Two sets of V8 multifunctional electrical instruments are used for field data acquisition. The construction layout is shown in Figure 9.

3.3. Quality Control and Analysis for the CSAMT Data

In the field data acquisition, the strength of the received signal and the duration of data collection are crucial to the quality of CSAMT data acquisition. Therefore, in this CSAMT field measurement, on the one hand, through white noise measurement experiments, the signal-to-noise ratio was calculated to judge the strength of the received signal, and the transmission signal strength was adjusted accordingly. On the other hand, through data collection time experiments, the quality of measurement data under different collection times was analyzed to find the optimal data collection time scheme.

3.3.1. CSAMT Signal-to-Noise Ratio Experiment

To calculate the signal-to-noise ratio (SNR) of CSAMT data, both the noise signal and the data signal need to be obtained. Firstly, the open-pit mine at the center of the work area was selected as the experimental site. Data acquisition was conducted when the artificial field source was not transmitting, and the measured signals were considered as noise signals. Secondly, under the same experimental setup conditions, when the artificial field source was transmitting, the collected signals were considered as data signals. Finally, the signal-to-noise ratio (SNR) formula was used to calculate the SNR. By comparing the SNR values, the quality of the collected data was evaluated, which further indicated the strength of the received signal. The experimental results of the signal-to-noise ratio were shown in Figure 10.

In Figure 10, there are four sub-figures: the electric field noise signal amplitude diagram, the magnetic field noise signal amplitude diagram, the electric field data signal amplitude diagram, and the magnetic field data signal amplitude diagram. From the experimental results, it can be observed that the noise signals have a large range of amplitude value fluctuations and significant error bars. On the other hand, the data signals exhibit smooth and continuous amplitude curves with small error bars, and their values are more than 10 times those of the noise signals. Therefore, the signal-to-noise ratio experiment we designed is reasonable and feasible, and can be used as a criterion for evaluating the strength of the received signal.

The calculation formula for signal-to-noise ratio is as follows:

$$SNR=20*Log{\displaystyle \frac{signal}{noise}}$$

(19)

The calculation results of the signal-to-noise ratio for the electric field and magnetic field at the test point are shown in Figure 11.

In the Figure 11, within the frequency band of 5000-10 Hz, the signal-to-noise ratio of both the electric field and the magnetic field is greater than 20 Db, and the noise percentage is less than 10%. This indicates that the signal strength is sufficiently large, and the quality of the measured CSAMT field data is high. In the high-frequency and low-frequency bands, the SNR is less than 10 Db, and the noise percentage is greater than 30%, indicating low quality of the collected data. This is due to two main reasons: on the one hand, the V8 instrument and magnetic probe exhibit response distortion in both high and low frequency bands, resulting in low measurement accuracy; on the other hand, when transmitting high frequencies, the transmission current is small, making the signal more susceptible to interference. In the case of low frequencies, factors such as shorter data collection time and stronger interfering signals contribute to a relatively lower signal-to-noise ratio in the data.

3.3.2. The Experiment of CSAMT Data Acquisition Time

During the field work, we also conducted experiments on acquisition time, and obtained apparent resistivity sounding curves for different acquisition times. Through comparative analysis, we determined the optimal acquisition time parameters.

Experimental conditions: With a single acquisition cycle of 40 minutes, three acquisition times of 40, 80, and 120 minutes were selected. Under the same survey line, survey point, equipment, and field source, three apparent resistivity sounding curves were obtained. The experimental results for different acquisition times are shown in Figure 9.

As can be observed in Figure 12, the apparent resistivity sounding curve with an acquisition time of 40 minutes exhibits poor continuity of resistivity values, especially below 30 Hz. The apparent resistivity curve appears jagged, and the characteristics of the transition zone and near zone are not pronounced, indicating poor data quality. Compared to the apparent resistivity curve of 40-minute, the apparent resistivity sounding curve with an 80-minute is smoother and more continuous, with smaller jumps. Below 30 Hz, it starts to enter the near zone, and the curve rises at an angle exceeding 45°. Meanwhile, there are jump points, and the transition zone characteristics are not pronounced, indicating intermediate data quality.The apparent resistivity sounding curve with 120-minute acquisition time is smooth and exhibits good continuity. The characteristics of the far zone, transition zone, and near zone are distinct. A minimum value appears at 10 Hz, which is a characteristic of the transition zone. Below 10 Hz, the sounding curve rises at nearly a 45-degree angle, which is a typical characteristic of a CSAMT sounding curve, indicating high data quality.

Therefore, in the Jinding lead-zinc mining area, by extending the data acquisition time and increasing the number of stacking, the characteristics of the apparent resistivity sounding curve at the measurement points can be significantly improved, thereby enhancing the quality of CSAMT single-point data.

4. CSAMT Data Processing and Inversion

4.1. CSAMT Data Processing

Before inversion, CSAMT data processing mainly includes the elimination of flying points and smoothing of the Cagniard apparent resistivity and phase data at each measurement point, as well as static correction and spatial filtering of apparent resistivity and phase along the entire survey line. Figure 13 shows the apparent resistivity frequency sounding curve obtained after data processing for some measurement points. From Figure 10, it can be seen that the apparent resistivity curve is smooth, continuous, and has a clear morphology. The characteristics of the far zone, transition zone, and near zone of the curve are obvious, which conform to the standards of CSAMT apparent resistivity sounding curves.

Figure 14 (a) shows the original pseudo-section of frequency-apparent resistivity for survey line 2100. It can be seen from the figure that due to the influence of uneven geological bodies and topographic fluctuations on the surface, the Cagniard apparent resistivity varies greatly and exhibits “noodle-like” characteristics in the vertical direction, indicating severe static effects that require processing. Figure 14 (b) is the result diagram of the frequency apparent resistivity pseudo-section after static correction and spatial filtering. It can be seen from the figure that after processing, the apparent resistivity data changes smoothly, the vertical “noodle-like” phenomenon disappears, and the static effects are well suppressed.

From the sounding curves of apparent resistivity at the measurement points and the pseudo-section diagram of frequency-apparent resistivity along the survey line, it can be seen that compared with the original data, the processed data has achieved significant suppression in terms of static effects. The data is authentic, reliable, and of high quality, making it suitable for further inversion.

4.2. CSAMT Inversion

Because the direction of the survey lines is perpendicular to the trend of the ore body, we can use two-dimensional continuous medium inversion software to invert CSAMT data. The principle of continuous medium inversion is based on the assumption of continuous variation of electrical properties in the subsurface. It discretizes the subsurface space using numerical methods and adjusts the resistivity values of each unit to achieve the best fit between the forward electromagnetic field data and the observed data. The method is closer to the real subsurface conditions and is suitable for exploration in complex geological settings, providing higher precision information on the subsurface electrical structure[26,27].

Based on the statistics of the work area, the frequency of measurement points entering the near zone is around 10Hz, and the surface resistivity ranges between 100-200 ohm-meters. According to the skin depth formula calculation, the exploration depth is approximately 1.5Km, which meets the exploration requirements. Therefore, we selected a total of 40 frequency points in the 8192-10Hz frequency band for inversion. The initial inversion model is the result of one-dimensional Bostick inversion. After 34 iterations of inversion, the data misfit is 2.7. The inversion results of each survey line are detailed in Figure 15.

5. The Electrical Structure Interpretation of the Jinding Lead-Zinc Mine

5.1. Physical Property Analysis

To interpret the electrical structure detected by CSAMT, this study used the two-dimensional inversion results of Line-1100 as a reference section. The main reason for selecting Line-1100 is that it has abundant borehole data, covering most of the stratigraphic lithologic units (including ore bodies) and major structural interfaces (including thrust nappes) in the mining area. The interpretation method we adopted is to combine the shallow subsurface electrical structure with borehole data, and then conduct a comparative analysis with the inversion results.

As shown in Figure 16, the profile exhibits a low-resistivity background. The central low-resistivity body primarily consists of fine sandstone, siltstone, limestone breccia containing lead-zinc ore and gypsum, with localized fault breccia zones.The center of the high-resistivity body is entirely composed of continuous fine sandstone and siltstone strata containing gypsum. Limestone, limestone breccia siltstone, as well as mudstone, sandstone, and siltstone containing gypsum exhibit moderate or high resistivity. From the above electrical patterns, it can be observed that the resistivity of geological bodies is controlled by the material composition and transport properties of the rocks themselves, with a more pronounced influence from the transport properties. Therefore, for the geological bodies in the Jinding mining area, due to the influence of large-scale tectonic and fluid activities, the distribution of their electrical conductivity is a response to the tectonic state and fluid type of the geological bodies.For example, rocks cemented by gypsum with low permeability exhibit high resistivity(Figure 17a-b ), while breccia cemented by calcium or sand with high permeability exhibit low resistivity(Figure 17c-d ). Relatively intact limestone and sandstone have moderate to high resistivity and occur as either independent stratigraphic systems or tectonically entrained rock slices (Figure 17e). Therefore, the electrical structures of the subsequent survey line profiles are interpreted based on the electrical characteristics of the aforementioned geological bodies.

5.2. Interpretation of the Electrical Structure

The relative positions of the survey lines on the plane are shown in Figure 4. Line-4500 and Line-3600 are located on the north side of the PaoMaping mining area; Line-2400, Line-2100, and Line-1900 pass through the Jiayashan mining area; Line 1100 and Line-700 are located on the west side of the Nanchang mining area. The north-south span of the survey lines is approximately 3.8km, and all survey line profiles are oriented east-west, which provides assistance for comparing the three-dimensional spatial positions of the underground electrical structures. Due to the significant differences and distinct contrasts in resistivity presented in the two-dimensional continuous medium inversion results, as well as the abundant details in the electrical structures, the interpretation of the electrical structures primarily relies on the two-dimensional inversion results.

3.

Line-2400

Line-2400, Line-2100, and Line-1900 are located in the Jiayashan mining segment of the Jinding mining area. As shown in Figure 18, the entire profile of Line-2400 can be divided into three parts, with the western and eastern segments being high-conductivity areas, and the middle segment (Region I) being a high-resistivity area.Region I is the mineral occurrence segment, which exhibits alternating high and low resistivity from shallow to deep, reflecting a complex structural-rock system. The eastern segment is located below the outcropping area of the Guolang Formation (E2g), which is mainly composed of alternating layers of argillaceous siltstone and siltstone-bearing mudstone, with a high conductor C1 underlying it. The strata exposed in the western segment belong to the Yunlong Formation (E2y), with normal lithology consisting of alternating layers of muddy conglomerate and argillaceous siltstone, and there is a high conductor C2 present in the deep part. Based on the electrical characteristics of the rocks, the presence of high conductors C1 and C2 indicates that both sets of strata have undergone varying degrees of disruption and modification, resulting in increased conductivity compared to intact lithologies. From the electrical gradient zones, it can be observed that exogenous rock slices have been thrusted from west to east, overlying the Yunlong Formation strata. The mineral bodies are controlled by the thrusted base, while the Guolang Formation contacts the exogenous rock slices through high-angle strike-slip faults.

4.

line-2100

Line-2100 is the longest survey line, located to the south of Line-2400. As shown in Figure 19, the electrical structure of the Line-2100 survey line is complex, being separated by five faults. Similar to Line-2400, Region II exhibits a high-resistivity background with alternating high and low resistivity from shallow to deep. However, the exposed area of the Yunlong Formation in Region IV on the west side is different from the western segment of Line-2400, showing high resistivity, which reflects a relatively intact J-K stratigraphic system. The difference between Region III and Line-2400 is that, although both are exposed areas of the Yunlong Formation, their electrical characteristics are completely opposite. The high-resistivity background indicates that the deep part of this region may have structural features similar to Region II. However, the presence of high conductor C5 corresponds to C2 in Line-2400. Similarly, high conductor C4 corresponds to C1 in Line-2400, reflecting the Guolang Formation. Additionally, there is a high conductor C3 on the west side of the Fengzishan mineral segment, which requires further investigation.

5.

line-1900

Line-1900 is located to the south of Line-2100 and also crosses the Jiayashan mineral segment. The Line-1900 profile is divided into three parts by faults. The low-resistivity zone C6 in the eastern segment, compared with C1 in Line-2400 and C4 in Line-2100, indicates the stable electrical characteristics of the Guolang Formation in the north-south direction. However, the presence of the Bijiang Fault enhances the extent of the high conductor C6. The moderate resistivity in the middle mineral segment (Region V) shows slight differences from the electrical structures of the previous two survey lines. The western segment (Region VI) can be compared with Region III of Line-2100, with the difference being that Region VI has a wider distribution, and the underlying high conductor C7 can be connected to C5 in Line-2100.

6.

Line-4500 and Line-3600

Line-4500 and Line-3600 are located to the north of the Paomaping mineral segment, as shown in Figure 21. In the electrical profile of Line-4500, Region VII reflects the J-K stratigraphic system. Similar to the previously mentioned survey lines, the Guolang Formation exhibits low resistivity (C8) and extends deeply downward. Line-3600 crosses three stratigraphic zones from west to east, which are T3, E2, and J2. The lateral electrical gradient zones also separate the fault occurrences in the three exposed stratigraphic zones. The high conductor C9 represents an extension of the Guolang Formation, and is connected to the deep part of Region VIII. The upper part may reflect an exotic thrust sheet, while the lower part belongs to the in-situ system. The vertical electrical gradient zone in Region VIII indicates the depth of the thrust surface.

7.

Line-1100 and line-700

Line-1100 and Line-700 are located to the south of the Jinding mining area and to the west of the Nanchang mineral segment. As shown in Figure 22, the shallow high-resistivity bodies in Line-1100 outline the distribution of the exotic thrust bodies, while the in-situ stratigraphic system may be located in Regions XI and XII. Line-700 is situated to the south of Line-1100 and can be divided into eastern and western segments. The low-resistivity background in Region XIII indicates that the exotic thrust system is relatively thin, with the in-situ system dominating at depth. In Region XIV, the surface exposes strata of the Yunlong Formation, and the alternating high and low resistivity features may indicate that the stratigraphic system has undergone a certain degree of structural deformation or experienced large-scale fluid activity. Its electrical characteristics are similar to the mineralization system observed in the previously mentioned profiles.

Based on the interpretation of the electrical structure from the seven surveyed lines, it can be understood that the lithological system of the Jinding lead-zinc mining area is mainly divided into three parts: the exotic stratigraphic system, the middle lithological system, and the in-situ system. The electrical characteristics of the three lithological systems are as follows: the shallow exotic stratigraphic system is mainly characterized by a high resistivity background; the middle lithological system, influenced by tectonic and fluid cementation, displays alternating high and low resistivity; and the in-situ system is relatively complex, with the Yunlong Formation and Guolang Formation exhibiting different electrical characteristics. The Yunlong Formation shows significant electrical variations, while the Guolang Formation maintains low resistivity and high conductivity, and high conductors(C1, C4, C6, C8, C9) are prevalent in the underlying geological bodies of the exposed areas of the Guolang Formation.

6. Conclusion

We have developed a new exploration mode on the CSAMT. Based on the geological overview of the mineral deposit, we established a geophysical model, and performed forward modeling using the coupled finite-infinite element method. On the one hand, we verified the effectiveness of CSAMT exploration based on the forward results, on the another hand, we designed the field acquisition parameters for CSAMT by referencing the forward modeling parameters

we proposed a new method for evaluating CSAMT data. Through the noise experiment, we quantitatively evaluated the strength of the transmitted signal by analyzing the signal-to-noise ratio, thereby ensuring the quality of data acquisition. At the same time, through the acquisition time experiment, we qualitatively evaluated the quality of the sounding curve, and obtained the optimal acquisition time parameters, and ensured the complete characteristic of the apparent resistivity curve

Using two-dimensional continuous medium inversion, the electrical distribution at a depth of 1km below each survey line in the Jinding lead-zinc mining area was obtained. At the same time, through the interpretation of the electrical structure of the survey lines, the electrical characteristics of the lithologic system in the Jinding lead-zinc mining area were revealed.

The CSAMT exploration mode adopted in this paper can provide some reference for CSAMT practitioners in designing field work and ensuring data accuracy.