Application of 3D Electrical Resistivity Tomography in the Yeoncheon Titanomagnetite Deposit, South Korea

Abstract

The Yeoncheon titanomagnetite deposit formed by Precambrian magma differentiation is located in Gyeonggi-do, South Korea. Our team conducted an airborne magnetic survey for multiscale mineral exploration and then selected a promising survey area. An electrical resistivity survey was carried out in the potential area to image subsurface structure. Because ore minerals are mainly distributed in gabbro monzodiorite rather than quartz monzodiorite, we applied three-dimensional inversion of electrical resistivity tomography (ERT) data to identify lithology boundaries related to magma differentiation. The resistivity criterion distinguishing the lithologies of gabbro and quartz monzodiorite was determined from laboratory resistivity experimental results performed on drilling cores. The selected region for gabbro monzodiorite extends to the northeast direction, which is consistent with the geology map, magnetic anomaly, and drilling data. The inversion results of ERT can help in selecting the location of geophysical survey or drilling.

1. Introduction

The titanomagnetite (TM) deposits in South Korea, which originated from intrusion into Precambrian metasedimentary rocks, are distributed in Yeoncheon, Boremdo, and Soyeonpyeongdo. As a co-product of steel production, the rare transition metal vanadium is extracted from TM deposits [1]. The main use of vanadium is as an additive to steel alloys, and it has been applied in the automobile, aviation, and defense industries because of its high strength and resistance at high temperature. Recently, the vanadium redox flow battery has been considered a suitable secondary battery for long-term energy storage systems [2]. Therefore, exploration of TM deposits primarily enables steel production but also contributes to vanadium extraction.

The target of this study was the Yeoncheon TM deposits, found in the Gonamsan igneous layered intrusion. The Gwanin magnetite mine, currently operated by Samyang Resources, produces 200,000 tons of iron and titanium annually. These igneous deposits are orthomagmatic deposits differentiated from alkali gabbroic magmas in the middle Proterozoic [3]. Several studies, including surface geology, geochemical, and geophysical surveys [4,5], have been conducted to investigate the development of the TM deposits. However, the distribution of iron–titanium mineralization has not been clearly delineated in the vertical and horizontal directions, and an intensive survey is needed to discover unidentified ore deposits.

Ore deposit exploration requires comprehensive approaches, including geophysical, geological, and geochemical analyses. Geophysical exploration detects anomalous areas by interpreting signals recorded by exploration equipment. The choice of exploration method depends on the characteristics of the ore deposits and the geological setting of the target area. In the case of TM deposits, airborne magnetic survey works well because such deposits consist of mafic and ultramafic igneous rocks [6]. However, the interpreted magnetic anomaly data have low resolution for determining the locations of ore deposits, especially in the vertical direction. For intensive surveying in potential areas, an electrical resistivity survey can be a suitable option as a subsequent exploration method because of the low resistivity of TM deposits.

Electrical resistivity tomography (ERT) is a conventional geophysical method that images the subsurface by measuring the electrical potential between electrode pairs [7]. The constructed electrical resistivity can be used for detecting target areas that have a different range of conductivity compared to surrounding areas. Because of its ease of use and low cost, ERT is used in various fields, including ground surveying and groundwater, mineral, and tunnel exploration [8,9].

Resistivity data are acquired along parallel lines and can be interpreted using 3D inversion methods [10,11,12,13,14,15]. The 3D inversion algorithm has been developed with different types of electrode array, regularization, modeling scheme and loss function [16]. For analysis of simple structures, two-dimensional (2D) ERT is cost-effective; however, in strongly heterogeneous environments, such as mountainous areas, the distortion in the direction perpendicular to the electrodes can result in unrealistic subsurface structures. Although three-dimensional (3D) ERT requires higher computation costs, reasonable subsurface imaging can be obtained by accounting for complex structures [12]. In this study, we applied 3D ERT to the Gonamsan igneous rocks to explore the unidentified Yeoncheon TM deposits.

2. Geophysical and Geological Background of Study Area

The Yeoncheon TM deposit in South Korea (Figure 1) has been developed since 1934. Because the extension of the TM deposit has not been clearly described, despite the long operation history of the magnetite mine, it is worthwhile to search for additional TM deposits.

2.1. Selection of Survey Area

To detect the magnetic anomaly including TM deposits, an airborne magnetic survey was conducted by the Korea Institute of Geoscience and Mineral Resources (KIGAM) in 2019 [17]. To eliminate the distortion due to the magnetic field inclination, we applied the reduction to the pole (RTP) technique to the total-field magnetic anomaly data. RTP technology helps the interpretation by making the anomaly more symmetrical and better-placed on the target. The magnetic anomaly data after RTP was overlaid onto the location map in Figure 1. Several locations, including the operating mine marked by the black dashed line in Figure 1, show higher amplitude than surrounding areas. The airborne magnetic survey covers a wide range of regions, but the inversion results are insufficient to describe subsurface structure. By applying the electrical resistivity method to highly anomalous areas, along with using drilling data, vertical information can be supplemented.

Among the candidate TM deposits, we selected the area with the GPS coordinates of 38.12° N and 127.22° E, which has good accessibility and high magnetic amplitude (Figure 1). We conducted the electrical resistivity survey along three lines to determine the distribution of resistivity.

2.2. Geological Setting

A schematic geologic map of the Yeoncheon TM deposit is shown in Figure 2. Iron–titanium oxide mineralization is developed in the Gonamsan igneous rocks, which intrude a Precambrian metamorphic complex. At the surface, this igneous intrusion extends about 3 km × 6 km with elongation in the north–south direction. It consists of quartz monzodiorite (QMD) on the west and monzogabbro–monzodiorite (GMD) on the east. The fractionation of iron–titanium oxide is closely related to the GMD, which differentiated from magma earlier than the QMD and contains more mafic rock. Most areas with high magnetic anomalies (Figure 1) are located in the GMD.

The Gonamsan igneous intrusion can be classified by lithofacies into lower, middle, and upper zones [18]. The upper zone consists mostly of silicate minerals, and the lower and middle zones contain iron-titanium ore bodies with alternations of silicate minerals. Therefore, it is necessary to investigate the subsurface structure in the lower and middle zones using geophysical data. Because the electrical resistivity survey was performed around the boundary of the QMD and GMD (Figure 2), 3D ERT could contribute to dividing the two types of monzodiorite and facilitate the discovery of TM deposits.

3. Electrical Resistivity Survey of the Yeoncheon TM Deposits

3.1. Data Acquisition

An electrical resistivity survey along three parallel profile lines was conducted to obtain precise information about the shape and extent of the anomaly in the airborne magnetic survey (Figure 3). The resistivity data were acquired with a SuperSting R8/IP (Advanced Geosciences, Cedar Park, TX, USA) resistivity meter. A cross-line survey was not conducted because of difficulties caused by the rough mountainous terrain and because of the short cross distance of the survey area. Although it is not clearly shown in Figure 1, the terrain near survey line 3 is very steep and the outcrop is exposed along the mountain slope. Therefore, the design of this survey could not cover the entire anomaly.

The survey lines had a length of about 680 m and included 35 electrodes, which were spaced at 20 m intervals. The dipole-dipole electrode array [19] was used for data acquisition. We chose this array because of its high resolution and the ease of deployment in the field. The low signal-to-noise ratio, a problem with this array, does not have a significant impact on data quality due to the high resistivity of near-surface in this area. Furthermore, to enhance the data quality and depth of interpretation, we measured electrical potential by changing the dipole spacing to 20, 40, and 60 m. In general, the measured potential and depth of interpretation increased with dipole spacing, but the spatial resolution decreased. The distance between survey lines was 120 m. For 3D analysis, the distance between survey lines is not appropriate compared to the electrode interval, but it is supplemented using the data with 40 and 60 m of dipole spacing. After all electrodes were installed, their locations were measured by high-precision GPS, and this information was used for inversion by modifying the node position to the measured electrode position.

3.2. Methods

We performed 2D and 3D inversion of the acquired resistivity data. In the 2D interpretation, DC2DPro version 1.01 [20] was used for processing. Data editing and 2D inversion with topography were carried out using this software. The software can perform both forward modeling and inversion using the finite element method with various kinds of constraints. We used the L2 norm for inversion with the second-order smoothness constraint [21] for processing of the acquired data.

For 3D resistivity inversion, we used a complex resistivity inversion algorithm to consider induced polarization effects [22]. This algorithm was constructed simply by extending the inversion variable in the conventional resistivity inversion to the complex value. Except for the type of data, all the procedures are the same as conventional methods, although they require more memory and computation. The complex version of Poisson’s equation can be defined as [23]:

$$-\nabla ·({\sigma }^{*}(\omega )\nabla {\varphi }_{\mathsf{\sigma }}^{*}(\omega ))=I\delta (r-r_{\mathrm{s}}),$$

(1)

where σ* and ϕ_σ* are the complex conductivity and electrical potential, respectively. ω is the angular frequency, I is the amplitude of current source, and r_s is the spatial coordinates of electrode. The objective function for resistivity inversion can be represented with the Lagrange multiplier (λ) as follows:

$$E(m^{*})=\left|\right|W_{\mathrm{d}}^{*}[d^{*}-f^{*}(m^{*})]|{|}^2+\lambda ||W_{\mathrm{m}}\Delta m^{*}|{|}^2,$$

(2)

where W_d* is a weighting matrix for data fitting, W_m is the roughness matrix and f* is a modeling operator based on finite element method. d* is the observed data, m* is the complex resistivity, and Δm* is the perturbation model. Taking the derivative with respect to the model parameter for the second-order Taylor expansion of misfit function results in:

$$(J^{*\mathrm{H}}{J}^{*}+\lambda W_{\mathrm{m}}^{\mathrm{T}}{W}_{\mathrm{m}})\Delta m^{*}=W_{\mathrm{d}}^{*}[d^{*}-f^{*}(m^{*})],$$

(3)

where J* is the Jacobian matrix and the superscript H denotes the complex conjugate. The Lagrange multiplier is the smoothness constraint, which helps model parameters update robustly in the iterative inversion process. The smoothness constraint is suitable for regions where the subsurface structure does not change rapidly. Therefore, this regularization can be applied to the exploration area in this study, which does not contain a discontinuous area such as a fault zone.

3.3. Results

Three sets of data were available for the single survey because we obtained resistivity data with different dipole spacing. We separately processed each set of resistivity data and compared the results for quality control. After combining the processed resistivity data for the three dipole spacing, 2D inversion was carried out. Figure 4 shows the final 2D inversion results for each survey line. In survey line 1, high resistivity zones are distributed at the upper layers than in other survey lines, and low resistivity anomaly is observed at bottom. The trend of electrical resistivity in survey line 2 is contrary to that in survey line 1, with higher resistivity in the south than in the north. The inversion result of survey line 3 increases with depth. Additionally, the range of resistivity is higher in survey line 3 than in the other survey lines. The 2D inversion results for each survey line can be interpreted independently, but these can also be used as the bases for 3D inversion.

For a comprehensive interpretation of the subsurface structure, we performed a 3D electrical resistivity inversion. We used only the data with a dipole spacing of 40 and 60 m considering the distance between survey lines. The coordinate transform of topography was applied, which was extracted from the digital map. The electrode locations for the 3D process were properly transformed from the 2D resistivity data.

The 3D geomodeling software GOCAD (Emerson, Houston, TX, USA) was used for visualization of the inversion results and the integration of multi-geophysical data. Figure 5 shows the 3D inversion results of electrical resistivity for the entire survey area obtained using SKUA-GOCAD version 2019. In the 3D resistivity images, no region corresponds to the low resistivity of ore deposits. However, due to the different resistivity ranges of the two monzodiorites, the inverted resistivity can be used to locate the boundaries of the GMD, where the iron-titanium oxide deposits are mainly developed.

Figure 6 shows pseudo-sections of the observed and predicted data of survey line 1 for 2D and 3D inversion. Because pseudo-sections of each survey show similar behavior, we include only the results of survey line 1. Because of the smoothness constraint, the predicted data shows smoother than the observed data. The predicted data is close to the observed data for both 2D and 3D inversion. We could verify that the inversion algorithm we used works well for the survey data in the exploration area. For 2D inversion, final root mean square (RMS) errors after 4 iterations for the data of survey lines 1, 2 and 3 were 0.1962, 0.0883 and 0.1656, starting from 0.7196, 0.4889 and 1.0024, respectively. For 3D inversion, final RMS errors after 10 iterations were 0.3130, where initial error was 1.8364.

4. Discussion

To interpret the interior of the 3D electrical resistivity inversion results, we used the drilling data obtained by KIGAM in 2019 [17]. A total of 4 drilling up to 300 m depth were carried out with azimuth angle of 12.6, 42.0, 70.4 and 149.6 degrees. The dip of drilling cores was about −60 degree. In this study, 2 drilling cores corresponding to azimuth angle of 42.0 and 70.4 degree were used, which contain both GMD and QMD. Figure 7 shows the lithology of drilling cores and their locations. The lithology was determined by the measurement of X-ray fluorescence. According to drilling cores, GMD is located above QMD, but it is a spatially limited interpretation.

There are two distinct trends of electrical resistivity, one with respect to elevation and another with respect to northeast direction, in Figure 5. Even when considering the low resistivity of the weathered zone in the shallow part, the electrical resistivity shows an overall increase with depth. Thus, most of the GMD region is located in the upper layers, corresponding with the borehole data, with QMD at the bottom and GMD at the middle and top. Another feature of the inversion result is that the electrical resistivity is lower in the northeast direction than in other regions. This tendency of resistivity toward the east is consistent with the geologic map shown in Figure 2, and the information in the north direction is revealed by the 3D resistivity structures.

For quantitative analysis of drilling data, the measured electrical resistivity results were obtained from laboratory experiments. Figure 8 shows box plots of the ranges of measured electrical resistivity according to rock type. Although the TM deposits were not included in the drilling data, the resistivity of the GMD differs obviously from the values of other groups. We separated the GMD from the entire domain based on the electrical resistivity value of 7000 Ωm, which is the minimum value of the QMD in the drilling data.

Figure 9 illustrates the regions where electrical resistivity is lower than 7000 Ωm. The selected area represents the candidate region of the GMD. We validated our results by matching the boundaries of the selected region with the lithologies of the drilling data (Figure 9). Furthermore, the selected region extends to the northeast direction, which is consistent with the geology map in Figure 2 and the magnetic anomaly in Figure 1. By referring to this result, we can determine the locations of the following geophysical survey and the drilling. Additionally, because it is difficult to obtain geological information in the subsurface, the classified region can contribute to analysis of the magma differentiation around the target area.

3D effects of inversion results can be quantitatively evaluated by analyzing the differences between 2D and 3D data. Figure 10 represents 2D sections of 3D inversion results along the acquisition profile. For survey lines 2 and 3, the overall trends in Figure 8 are similar to those of 2D inversion results in Figure 4. However, for survey line 1, 2D and 3D inversion results differ significantly. The low resistivity anomalies at the bottom in 2D inversion results are not observed in 3D inversion results. Survey line 1 corresponds to QMD areas with relatively high resistivity on the geological map in Figure 2. In addition, compared to the magnetic anomaly in Figure 1, the low resistivity zones are less likely to be observed. This can be seen as an error that can occur when interpreting a single-line survey.

The reliability of the subsurface structure in the 3D inversion is affected by the number of 2D survey lines. The numerical inversion in this study was conducted using only three survey lines, without a cross-line survey. Although we accounted for geometric effects by changing the dipole spacing for each line survey, relatively limited survey data are a limitation of our inversion results. Furthermore, because data are not acquired by a 3D coordinate system, they are not true 3D ERT, and there is a fundamental problem with lateral resolution. A sensitivity test can be a way of verifying inversion results, and in this study we supplemented this issue by changing dipole spacing. Two boreholes are adjacent, and the number of drillings is insufficient to represent the target area. Additional surveying and drilling would improve the resistivity images, and multi-geophysical data, including induced polarization and local magnetic survey data, can be used for cross-validation.

5. Conclusions

We applied 3D ERT to explore the Yeoncheon TM deposit in the Gonamsan igneous intrusion. This intrusion is classified as QMD and GMD according to the ratio of mafic rocks. The location of the survey area, which was determined by an airborne magnetic survey, is near the boundary between the two monzodiorites. Because iron–titanium mineralization is mainly developed in the GMD, separating the two monzodiorites provides helpful information for identifying ore deposits. The electrical resistivity survey was carried out along three parallel survey lines using dipole spacings of 20, 40, and 60 m. We performed 3D electrical resistivity inversion from 2D resistivity models to investigate the complex mountainous terrain. From the 3D resistivity inversion results, we were able to separate the region of the GMD based on the measured resistivity of the drilling data. There were two distinct trends in the resistivity model. First, electrical resistivity increases with depth, which suggests that the GMD is located above the QMD, consistent with the lithology borehole data. Second, the GMD is distributed in the northeastern part of the survey area, which corresponds to the location of monzodiorite in the geologic map. The separated 3D resistivity model for the GMD facilitates determining the locations and directions of subsequent geophysical surveys for exploration of the ore deposit.