Three-Dimensional P-wave Velocity Structure of the Zhuxi Ore Deposit, South China Revealed by Control-Source First-Arrival Tomography

Abstract

The Zhuxi ore deposit, located in Jiangxi province, South China, is the largest tungsten reserve in the world. To better understand the geological structure and distribution of orebodies, we conducted a high resolution three-dimensional P-wave velocity tomography of the uppermost 0.5 km beneath the Zhuxi ore deposit and adjacent area. Our velocity model was derived from 761,653 P-wave first arrivals from 998 control-source shots, recorded by a dense array. As the first 3D P-wave velocity structure of the Zhuxi ore deposit, our model agrees with local topographic and tectonic structures and shows depth-dependent velocity similar to laboratory measurements. The Carboniferous formations hosting the proven orebodies are imaged as high velocities. The high-velocity anomalies extend to a larger area beyond the proven orebodies, and the locations of high–low velocity boundaries are in accordance with the boundaries between the Neoproterozoic formation and the Carboniferous–Triassic formation. Seismic tomography reveals that high-velocity anomalies are closely related to the mineralized areas. Our results are helpful for further evaluating the total reserves and suggest that seismic tomography can be a useful tool for mineral exploration.

1. Introduction

The Zhuxi ore deposit, located in the eastern part of the Jiangnan orogenic belt in Jiangxi Province, China, is a skarn-type mineral deposit [1]. The geological survey and borehole inspection in 2016 [2] estimated the tungsten reserves of this deposit to be more than 3.44 million tons, making it the largest proven tungsten reserve in the world [1].

The Zhuxi area experienced superimposed transformation from multi-stage tectonic–magmatic events from the Neoproterozoic to the late Mesozoic era [3,4,5]. The activity during the Yanshanian period was the most violent, which resulted in regional diagenesis and metallogenesis and was characterized by stretching, shearing, crustal thinning, granite intrusion, crustal rupture, and fluid activity [6,7]. The thrust nappe structure can cause interlaminar nappe slip of the ore-bearing stratum and form weak zones that are conducive to magma emplacement and orebody storage [8]. The Zhuxi ore deposit was formed after the multi-stage tectonic evolution of oblique intrusion of granitic magmas, skarn mineralization, hydrothermal cooling and alteration, and precipitation of metal sulfides [7]. The major orebodies beneath the Zhuxi ore deposit are located in the contacting zones between the Huanglong Formation and the Neoproterozoic slate (Figure 1), bounded by several NE-striking faults [1]. The ores are mostly hydrothermal-type copper–iron ores at shallow depths, whereas there are rich stratoid skarn scheelite (copper) ores at greater depths [7].

Until now, the reserves and distribution of the Zhuxi ore deposit have been mainly evaluated from the geological survey [2] and borehole explorations [1,7]. Two major ore belts—the Main Ore Belt (MOB) and the North Ore Belt (NOB) (Figure 1)—have been identified. As revealed by borehole inspections, most orebodies are NE-trending and dominated by a veinlet-disseminated structure [7]. Due to the limited number of inspecting boreholes, however, some geological problems are not fully understood. Such as the stratigraphic boundaries between different formations, the spatial distribution of the orebodies outside the prospecting survey area, and the shape of the concealed orebodies at depth reconstructed by nappe structure. These limited understandings restrict the study and prospecting of metallogenic regularity in the Zhuxi ore deposit and adjacent area. Imaging the subsurface structure with geophysical methods is helpful in overcoming these plights.

Recently, several geophysical explorations (e.g., gravity–magnetic, electromagnetic, and aeromagnetic survey) have been carried out in Zhuxi [9,10]. The ore zones are shown to have low Bouguer gravity anomalies, positive aeromagnetic anomalies, and high magnetic derivative norms [9,10]. However, gravity–magnetic, electrical, and electromagnetic explorations were mainly used to delineate the orebodies targeting horizontally, and it still need borehole inspection to find the depth and strata [11]. Three-dimensional fine images are still necessary to locate the orebodies and precisely estimate the reserves.

Seismic reflection and refraction are effective tools for exploring the subsurface structure and detecting reservoir related reflective interfaces [12,13]. However, seismic reflection is rarely used in orebody explorations because most metal mines have complex structures and seismic records in metal mine areas are subject to low signal-to-noise ratios, weak reflections, and strong scattering [11]. Fortunately, laboratory measurements indicate that strong seismic velocity differences exist between the underlying and overlying formations around the Zhuxi ore veins. Therefore, the fine-scale 3D velocity structure images around the Zhuxi ore deposit will be helpful in revealing the distribution of the metamorphic belt and assessing mineral reserves.

To image the fine-scale structure around the Zhuxi ore deposit, we conducted a geophysical experiment in Zhuxi and the adjacent region in November and December of 2017 (hereafter, referred to as the Zhuxi experiment). In this study, we present the 3D velocity structure around the Zhuxi ore deposit obtained from the control source seismic experiment data.

2. Zhuxi Experiment

From November to December 2017, a multi-disciplinary experiment was conducted around the Zhuxi ore deposit, including seismic, gravity, and electric surveys [14]. During this experiment, we deployed one 2D array (black triangles in Figure 1) and four dense seismic profiles (blue triangles in Figure 1) for the seismic survey. The 2D array covered an area of 30 km × 40 km with 178 portable three-component short period (5 s–150 Hz) seismometers (manufactured by the Chongqing Geological Instrument Co., Ltd, Chongqing, China) with station spacing of 1–5 km. The four profiles were 26.5 km long in total, equipped with 2311 vertical component 10 Hz geophones (manufactured by the Sercel, Nantes, France) with an average station spacing of 11.5 m.

Four types of control sources, namely, vibrator (28 tons force; supplied by the BGP Inc., China National Petroleum Corporation, Zhuozhou, Hebei, China), electric hammer (800 kg force; supplied by the Institute of Geophysics, China Earthquake Administration, Beijing, China), airgun (2000 in³; supplied by the Institute of Geophysics, China Earthquake Administration, Beijing, China), and methane detonation source (supplied by the China Academy of Engineering Physics, Mianyang, Sichuan, China) [19], were used to generate seismic waves. The vibrator is an electronically controlled, hydro mechanical system driven by a servo-valve assembly that exerts ground force through its baseplate against the Earth’s surface [20]. During the Zhuxi experiment, the vibrator was fired along four profiles at 738 points, with an average spatial interval of 40 m. At each point, we lunched one 12 s long sweep with frequency 1.5–96 Hz at 28 tons force. The hammer generates seismic waves by hitting a square iron plate fixed to the ground [21]. The airgun source can rapidly release compressed air to generate a pressure pulse and relatively long period air bubble oscillation that is then transmitted as seismic waves [22,23]. After being fired in a sealed container, methane detonation produces seismic waves when the high pressure air is quickly released to impact the surroundings. The methane detonation source was used for the first time in this experiment [19].

The vibrator and hammer sources were fired along profiles 1, 2, and 3 (Figure 1b) and were synchronized with the receivers (Sercel 428XL) with a wire connection. The seismometers on profile 4 and the 2D array (three-component, short-period seismographs) were continuously recorded and synchronized with the Global Positioning System (GPS; manufactured by the Trimble, Sunnyvale, CA, USA) timing. In addition, the origin time of the source was recorded using the GPS timing system with millisecond accuracy. The database using in the paper comes from the observations supported by the Institute of Geophysics, China Earthquake Administration, which can be requested from the authors.

3. First Arrival Picking

To remove contamination from the source time functions and improve the signal-to-noise (SNR), we first preprocessed the waveform data. For the vibrator, we reconstructed the record from the continuous signals using a cross-correlation method [24]. The airgun signals were linearly stacked to enhance the SNR [25,26], exploiting the high repeatability of the airgun source [23,27]. We then manually picked first P-wave arrivals from raw waveform for different sources and the P phases were clearly identified for all sources (Figure 2).

We also filtered the seismic records to respective predominate frequency bands (determined from the spectrum characteristics analysis) for each seismic source to compare their waveform characteristics and propagation capabilities (Figure 2). The longest distance the airgun can propagate 20 km with relatively lower frequency. The vibrator source has clear P-wave first arrival phases with propagation distances up to 6 km. The methane detonation source can realize a detection of 10–5 km with higher frequency than the airgun source, which is suitable for small-scale fine-structure exploration [19]. The hammer signal has the highest frequency with the shortest propagation distance, and the SNR of the waveform is relatively lower compared with the signals from other sources (Figure 2). We did not use the low SNR data (Figure 2d) from the hammer sources in this study.

The manually picked travel times are distributed around a linear trend (black line in Figure 3a) with an average velocity of 5.62 km/s (obtained from the least-square fitting). There are some outliers in the phase picking, mainly due to poor SNRs (Figure 3a). We only kept the picks with deviations less than 0.3 s (discarding 1397 phases). In total, 761,653 reliable P-wave arrivals (Figure 3a) were used for further tomographic inversion.

4. Method and Input Velocity Model

4.1. Inversion Method

In this study, the 3D shallow P-wave velocity structures beneath the Zhuxi ore deposit and the surrounding areas were obtained by using the simul2000 program [28,29,30,31], which has been widely used in local seismic imaging [32,33,34].

In the simul2000 inversion, the arrival time residuals are represented as functions of seismic velocities:

$$r_{ij}^{}=t_{ij}^{obs}-t_{ij}^{calc}$$

(1)

$$r_{ij}^{}=\frac{\partial T_{ij}}{\partial x}\mathsf{\Delta }x+\frac{\partial T_{ij}}{\partial y}\mathsf{\Delta }y+\frac{\partial T_{ij}}{\partial z}\mathsf{\Delta }z+\mathsf{\Delta }{\tau }_i+{{\displaystyle \sum }}_{k=1}^{k=total}\frac{\partial T_{ij}}{\partial m_k}\mathsf{\Delta }{m}_k$$

(2)

where $t_{ij}^{obs}$, $t_{ij}^{calc}$, $r_{ij}$, and $T_{ij}$ are the P-wave observed arrival time, calculated arrival time, arrival time residual, and travel time from event $i$ to station $j_{,}$ respectively. The earthquake location coordinate is $\left(x,y,z\right),$ and $m_k$ is the seismic velocity of the $k$ th grid node. The perturbation of a parameter from its initial value is denoted by $\mathsf{\Delta }$. Since we used control sources in this paper, the first four items in Equation (2) vanished. The simul2000 algorithm uses damped least-squares inversion to solve for the model perturbations and can calculate the covariance matrix to estimate the resolution of the final model.

4.2. Input Velocity Model

Seismic tomography is generally a linearized approximation of a nonlinear problem [31]; therefore, the initial velocity model is of key importance in the inversion process. The optimal one-dimensional (1D) P-velocity model can be obtained using the VELEST package [35,36]. This algorithm seeks the optimal 1D (layered) velocity model that minimizes the difference between the calculated and observed travel times [36] and is widely used for initial reference models in 3D local earthquake tomography (e.g., Matrullo et al. [37]).

We selected 19 starting 1D P-wave velocity models with different surface (0 km) velocities and velocity gradients (Figure 4) for the VELEST inversion procedure. Starting from these different initial models, all inversions converged to similar structures (red lines in Figure 4) after 20 iterations. The 19 inverted velocity models were further averaged and smoothed as the final 1D velocity model (green line in Figure 4) for 3D tomography.

5. Three-Dimensional P-Wave Velocity Structure

5.1. Model Setup and Parameter Selection

In inverting Equation (2), we set horizontal grid spacing of 0.5 km and 0.1 km spacing in depth (pink dots in Figure 3b). Considering that the highest elevation station in our study was 0.45 km above sea level, we placed an unmodeled layer at 0.5 km above sea level.

To stabilize the inversion of Equation (2), a damping factor is generally used, which plays a role in the inversion process. Damping factors are often selected by applying a trade-off analysis between data misfit and model variance [33,38]. We explored a wide range of damping factors to determine the optimal damping parameter (Figure 5). According to the trade-off curve, we chose 300 as the damping value for our tomography, producing a good compromise between data misfit and model variance.

5.2. Model Quality and Resolution Test

After 12 iterations, the tomographic inversion converged with the root–mean–square of arrival time residual reduced from 0.03 s to 0.01 s (Figure 6). Before the inversion, the residual varied from −0.06 s to 0.09 s with an average value of approximately 0.03 s, suggesting a systematic deviation in the initial velocity model. After the inversion, the residual followed a Gaussian distribution centered at 0 s.

To assess resolving power of the data and parametrization, we performed a checkerboard resolution test with the same inversion parameters as for the real data. We computed synthetic travel times through the 1D velocity model (green line in Figure 4) with alternative ±5% velocity perturbations across two adjacent grid nodes. The input and inverted-P-wave velocity structures are shown in Figure 7. The resolution matrix which is a stricter criterion than ray density (Lin et al.), [32,33,34] was calculated to estimate the resolution of the model and the uncertainties in the model parameters. The green contours enclose the well-resolved area with the diagonal element of the resolution matrix greater than 0.1 (1.0 represents the best resolution and 0.0 is not resolved at all [31,33]). As shown in Figure 7, the velocity structure of the uppermost 0.5 km is well resolved at the central part of the survey area.

5.3. Three-Dimensional P-Wave Velocity Structure

At the surface (0 km), the velocity structure shows a clear banded distribution, consistent with the geologic faults (NE-trending faults) and topographic distribution, as shown in Figure 1 (i.e., the high-velocity anomalies correlate with the high topography). From 0 to 0.5 km, the most notable features in our model are the high-velocity anomalies beneath the MOB and the NOB (Figure 8) (including Zhuxi, Tanling, and Hongmeiling). The MOB for the Zhuxi ore deposit corresponds to the high magnetic gradients, low Bouguer gravity anomalies [9,10], and low resistivity anomalies [14]. In our model, the MOB is imaged as a high-velocity anomaly from the surface to a depth of 0.5 km, and the majority of the high-velocity anomalies exist in Permian and Carboniferous formations (Figure 8). The NOB is a small copper–molybdenum deposit in the cretaceous granitic porphyry [7,39,40]. Our velocity structure presents some sporadic high-velocity anomalies (black rectangles in Figure 8) beneath the NOB.

6. Discussion

6.1. Absolute Velocities of the Orebodies

As previously mentioned, many boreholes were drilled down to 2.2 km [41] to inspect the distribution of orebodies. These inspections provided an unprecedented opportunity to calibrate our velocity model. We collected 78 rock samples recovered from different strata along five prospecting lines (green lines in Figure 9a). The depth of rock samples ranged from the surface to 2.2 km (with 11 samples shallower than 0.6 km). P-wave velocities (black dash line in Figure 9b) were measured under atmospheric pressure and room temperature using the method proposed by Xu et al. [42].

We then compared our velocity profiles with the laboratory measurements. Three velocity profiles (initial 1D, inverted average of whole study area, and inverted average around the prospecting area) were used for comparison. The average velocity around the prospecting area (the magenta line in Figure 9b) is higher than the initial 1D model and the average velocity from 3D tomography, which reveals that the orebodies have higher velocities than other areas. The velocity profiles from the laboratory measurements (black dash line in Figure 9b) and our inverted model (magenta line in Figure 9b) have similar depth dependencies. The velocities measured in the laboratory are generally lower than the tomographic results (especially at 0.4 km), probably because the laboratory measurements were affected by pores and micro-cracks introduced during the core recovery and sample preparation. Furthermore, the laboratory measurements were performed under atmospheric pressure much lower than the in situ pressure. At 0 km depth, the tomographic velocity is lower than the laboratory measurements partly because of the wide spreading of surface sediments. At 0.5 km depth, P-wave velocity are an estimated 6 km/s from both tomography and laboratory.

6.2. P-Wave Velocity Beneath Prospecting Lines

The borehole inspections along several lines also reveal the spatial variation of the lithology, and lines 54 and 42 (Figure 9) are the most well-studied [7]. In this section, we further compare our velocity model with these two sections across the Zhuxi ore deposit (Figure 10).

High-velocity anomalies are observed beneath the MOB for both prospecting lines 54 and 42 (Figure 10b,d). According to the borehole explorations, the Carboniferous Huanglong Formation and the Permian Qixia Formation beneath the MOB are the dominant reserve layers of the copper–tungsten ore from surface to deeper (~2 km), containing abundant Cu, Zn, Fe, W, Mo, and other elements [1,10,39,43]. These high-velocity anomalies (with red arrows, most in the Carboniferous formation) correspond to the hydrothermal vein-type copper and iron ores beneath the MOB in the shallow part (<0.2 km) and a small part of tungsten orebody beneath 0.2 km [1,7]. Our results also reveal that the orebodies of prospecting line 42 are shallower than those of prospecting line 54, which is in accordance with the results from borehole inspection [7,8].

Our velocity results show that there are obvious high–low velocity boundaries between the Neoproterozoic (Pt3 in Figure 10b,d, low velocity) and the Carboniferous–Triassic (C–T in Figure 10b,d, high-velocity), and the locations of velocity boundaries are consistent with stratigraphic classification (black lines in Figure 10) from geological survey and borehole inspection [8]. These velocity contrasts are in accordance with lithologic settings, where the Neoproterozoic is mostly phyllite with argilloarenaceous rock and lightly metamorphosed rock, and the Carboniferous–Triassic formation is mainly composed of carbonates, limestone, and dolomite [1,7]. Shi et al. [14] also revealed that the average resistivities for the Neoproterozoic are higher than those of Carboniferous–Permian. The observed velocity contrasts extend NW and deep beyond the MOB (Figure 10b,d), and this trend may depict the depth extension of stratigraphic boundaries (thick dashed solid lines in Figure 10).

6.3. Cross Sections of the Whole Zhuxi Ore Deposit

The correlation between the high-velocity structure and mineralized zone shown in the last section suggests that velocity structure is a good indicator for stratigraphic classification. To determine the spatial distribution of the high-velocity zone, we present seven more cross sections (CC’ to II’ parallel to AA’ and BB’).

Similar to the prospecting lines 54 and 42, high-velocity anomalies are observed beneath the MOB and NOB for all seven cross sections (Figure 11). The high-velocity anomalies along the CC’ to FF’ beneath the MOB (Figure 11b–e) are considerably higher than those along other cross sections. There is a high-velocity zone (red circle in Figure 8e) on the south edge of the MOB at 0.4 km (Figure 8e), where a low electrical resistivity structure is observed to be extending deeper than 2 km [14]. This is attributed to the Tanling deposit and may be a large concealed magmatic rock [39]. The range of the high-velocity zone beneath the NOB is significantly smaller than the MOB, which is consistent with the orebody reserve estimations. Other geophysical results also indicate that the MOB has lower Bouguer gravity anomalies and a higher magnetic field gradient than the NOB [9,10]. In addition, there are some high-velocity anomalies in the southwest of the NOB (black rectangles in Figure 8), which may correspond to the sporadic distribution of the cretaceous granitic porphyry. The low-velocity zones (0–0.2 km) beneath the MOB (Figure 11b,e,f) are likely the Quaternary sediment [1,7].

The high-low velocity boundaries are visible and can be tracked in the cross sections (Figure 11), and these velocity boundaries may indicate the stratigraphic boundaries (black dotted lines in Figure 11) between the Neoproterozoic and the Carboniferous–Triassic. From all cross sections (including prospecting lines 54 and 42), we suggest that the MOB can span 2 km in accordance with the results of borehole explorations [7]. Even though the high-velocity anomalies between the MOB and the NOB become weaker from northeast to southwest, all of these belts present a north-dipping structure. Our results also indicate that the velocity beneath the Tanling is higher than Zhuxi (Figure 11). Shi et al. [14] also revealed that Tanling has lower resistivity than Zhuxi and suggested that the deep part of Tanling is probably of copper mineralization rather than tungsten mineralization.

6.4. Formation of the Zhuxi Ore Deposit

The borehole explorations show that the high-velocity bodies in the shallow part of the Zhuxi ore deposit (Figure 8, Figure 10 and Figure 11) are mainly chalcopyrite and pyrite [40]. There is a large amount of pyrite distribution (especially in the metamorphic rocks). These copper and iron orebodies in the shallow part are mainly lenticular and vein-type. The lenticular-type ore is caused by the difference in the ore-bearing space formed by tectonism and the vein-type ore is mainly originated from hydrothermal activity. Mineralization mainly occurred in the Carboniferous formation and the shallow part of the Permian Qixia Formation, whereas little metallogenic activity is observed beneath the Triassic formation. Whether the other high-velocity anomaly zones in the Permian have been mineralized needs further analysis from geological, geochemical, and geophysical aspects. Two mineralization mechanisms exist for the copper–tungsten and iron deposits; the hydrothermal activity is the dominant mechanism, and the sedimentation is the secondary one. The hydrothermal activity is described in detail below (Figure 12).

The granitic magma activities during the Yanshanian period [43,44,45,46,47] bring ore-bearing hydrothermal fluid for mineralization [1], and there are a series of faults (especially the NE-trending thrust faults) forming possible upwelling channels. In the early stage of mineralization, magma and ore-bearing hydrothermal fluids originated from the high-temperature and high-pressure environment. During the upward migration, hydrothermal fluids continuously extract metal elements (mainly tungsten) from carbonate rocks and then alter into scheelite during the latest stage of magma crystallization [48]. As the ore-bearing hydrothermal fluids continue to rise, the temperature and pressure will gradually decrease. A large portion of copper–iron–zinc minerals will be separated as sulfides, and then copper–tungsten (zinc) orebody and iron and copper (zinc) orebody are formed by metasomatism [1,49].

In short, the MOB and NOB were formed under the superimposed effects of tensile structure, magma upwelling, hydrothermal metasomatism, and mineralization [1,7,8,50]. Our 3D body wave tomography results suggest the extension of stratigraphic boundaries outside the prospecting survey area and reveal that these potential mineralized areas may be gradually connected in the deep area beneath the MOB and NOB. Therefore, the MOB and NOB may have originated from the same hydrothermal system.

Subject to the low power of seismic sources, our first-arrival P-wave tomography only presents the velocity structure of the uppermost 0.5 km. In a future study, we will use the continuous records from the 2D array to perform an ambient seismic noise tomography. The ambient seismic noise tomography is expected to provide S-wave velocity structure coarser but deeper than the current P-wave velocity model presented here. In addition, the S-wave velocity structure will provide more constraints on the deep distribution of tungsten orebodies.

7. Conclusions

In this study, we imaged the high-resolution 3D velocity structure beneath the Zhuxi ore deposit using the data of the Zhuxi experiment conducted in 2017. Our velocity model shows a clear banded distribution, consistent with the local topographic and tectonic structures, and a similar depth-dependent velocity trend with the laboratory measurements of rock samples. The proven orebodies, mainly related to magmatic hydrothermal activities during the Yanshanian period, are shown as high-velocity anomalies beneath the MOB and NOB, corresponding to widespread copper–iron and a few tungsten–molybdenum orebodies. According to the distribution of high-velocity anomalies, we suggest that the stratigraphic boundaries extending outside the prospecting survey area (presenting north-dipping structure) and these potential mineralized areas are probably connected at depth. Our results can help in further evaluating the total reserves, suggesting that seismic tomography can be a useful tool for mineral exploration.