The audio-frequency magnetotelluric (AMT) method, along with controlled-source AMT (CSAMT), has been widely used to map shallow resistivity structures with high resolution. For example, it has been employed to trace the hydrological system of thermal fluid flow in Bakreswar hot spring, India (Sinharay et al., 2010); to image high-conductive zones at mid-crustal depths in the Otjiwarongo and Katima Mulilo regions of Namibia (Share et al., 2014); to delineate the shape of resistive granitic intrusions and their conductive surrounding metamorphic zones in West Junggar, NW China (Yang et al., 2016); to investigate one of the major metallogenic belts in southern China, including iron and polymetallic deposits in the Longmen region (Hu et al., 2013); and to explore skarn-type ore belts in the Zhuxi polymetallic deposits (Shi et al., 2020). In the AMT procedure, we typically measure two orthogonal electrical field components (i.e., Ex and Ey) and two or three orthogonal magnetic field components (i.e., Hx, Hy, Hz) of the naturally induced electromagnetic wave in the frequency range of 20,000 Hz to approximately 1 Hz (i.e., audio-frequency). These measurements are used to estimate the impedance tensor and the vertical transfer functions, which are directly related to the electrical resistivity of the Earth. Modern 3D inversion techniques have significantly enhanced the resolution and reliability of resistivity models inverted from AMT data. However, to date, only a few 3D results have been reported. (e.g., Yang et al., 2016; Blake et al., 2016; Wang et al., 2017; Shi et al., 2021; Wan and Wang, 2023, etc.).

As the first porphyry Cu-Mo-Au deposit discovered in Western Junggar, NW China, Baogutu has garnered significant attention in recent years (e.g., Shen et al., 2010; Shen and Pan, 2013; Cao et al., 2014; Cao et al., 2015; Zheng et al., 2015; Cao et al., 2016; Cao et al., 2017, et al.). Previous studies on the metallogeny of Baogutu stocks and porphyries have predominantly relied on geological investigations. For example, Cao et al. (2014) suggested that Baogutu can be classified as a typical reduced porphyry Cu deposit. However, geophysical investigations have rarely been conducted, limiting our understanding of the deep characteristics of the deposit. This motivates us to utilize dense AMT surveying to evaluate the potential for further exploration.

Here, we present the results from imaging the 3D electrical resistivity structure of the V (5th) porphyry stock in the Baogutu area, using 176 AMT sites and 3D inversion techniques. The inverted 3D resistivity model has been interpreted in conjunction with the shear wave velocity model, as well as the electrical resistivity and shear wave velocity of rock samples, which show good agreement.

The Western Junggar region in NW China is situated in the southern portion of the Central Asian Orogenic Belt (CAOB, Figure 1). This area occupies a late Paleozoic convergent zone bounded by three major plates: the Tarim craton to the south, the Siberia craton to the northeast, and the Kazakhstan block to the northwest (Jahn et al., 2004; Safonova et al., 2004; Chai et al., 2009; Zhang et al., 2009; Zhao et al., 2009). The crustal growth of the CAOB, one of the most deformed areas in the Phanerozoic, has been accompanied by multi-stage and diverse crust-mantle interactions and ore-forming processes (Windley et al., 2002; Yakubchuk, 2004; Buslov et al., 2004).

The Western Junggar accretionary collage is characterized by widespread Devonian and Carboniferous island arc-related volcanic rocks, granitoids, and oceanic sediments (Buckman and Aitchison, 2004), as shown in Figure 1. The Devonian strata are mainly distributed in the northwest, while the Carboniferous strata are primarily found in the southeast (Tang et al., 2010). The most extensively distributed Lower Carboniferous layers in this area consist of three units: the Xibeikulasi Group (C1x), Baogutu Group (C1b), and Tailegula Group (C1t), from bottom to top (Zheng et al., 2009; An and Zhu, 2009; Zong et al., 2014; Zhang et al., 2015). The Late Carboniferous to Early Permian period is critical for mineral deposit formation (Yin et al., 2011). Surface exposures of ultrabasic and acidic igneous intrusions indicate that magmatism was widespread throughout the Western Junggar area during this time (Figure 1; Chen and Jahn, 2004; Chen and Arakawa, 2005; Wang and Xu, 2006).

There are two groups of intermediate-acidic intrusions in the region. The first group consists of barren granite batholiths that intruded into Devonian to Lower Carboniferous volcanic rocks and were dominant around 300 Ma (Su et al., 2006). The second group includes intermediate-acidic stocks that intruded into the Carboniferous Darbut volcanic belt, primarily between 322 and 305 Ma (Shen and Jin, 1993). Separated by the Darbut fault, alkali-feldspar granite intruded as batholiths to the north of the fault (e.g., Akebasitao and Miaoergou granitic intrusions), while the southern side of the fault is dominated by intermediate-acidic small intrusions (e.g., Baogutu and Karamay granitic intrusions) (Yin et al., 2011). Frequent and large-scale magmatism, along with simultaneous and/or late tectonic deformation and particularly multi-stage superimposed magmas, contributed to the development of the Hatu, Saertuohai, and Baogutu Cu-Au deposits in West Junggar (Tang et al., 2010; Ma et al., 2012). Based on the geological structure and evolutionary history, ore deposits and mineralization have been linked to hydrothermal activity, ultrabasic rocks, and intermediate-acidic small intrusions (Cheng and Zhang, 2006). However, the nature of the granitoids and the genesis of these deposits remain subjects of debate (Cao et al., 2016).

The faults and folds are well developed in the study area (Figure 2A), with Stock V located on the eastern flank of the Baogutu anticline (Figure 2B; Yang, 2008; Zhang et al., 2010). There are three types of faults identified: NS-trending faults from the early stage, NE-trending faults from the middle stage, and sub-parallel NS-trending faults from the later stage. The early extension-shear faults are larger, extending over 10 km, and control the distribution of intermediate-acidic rocks. The middle NE-trending faults cut through the earlier NS-trending faults and intrusions, partially affecting porphyry formations. These NE-trending faults are parallel to most of the intermediate-acidic dikes and ore-bearing quartz-sulfide veins in the Baogutu Cu-Au deposits, serving as critical ore-controlling structures in the area (Wei et al., 2014). Previous studies suggest that the early NS-trending faults provide space for Cu-Au mineralization, while the NE-trending faults generally control the development of small stocks. The late-stage NS-trending faults are smaller and manifest as small cracks, which serve as the main ore-hosting structures for mineralization, often partially filled with dikes (Yang, 2008).

The Baogutu porphyry copper belt is located south of the Darbut fault in the western Junggar, approximately 40 km from Karamay city. This belt represents the eastern extension of the Tuoli Cu-Au belt in Xinjiang (Shen et al., 2009; Zhang et al., 2006). The deposits are primarily distributed in the contact zones of Hercynian intermediate-acidic intrusions and their associated dikes (Cheng and Zhang, 2006). Porphyry copper deposits develop within intermediate-acidic stocks, resulting in mineralization with thicknesses ranging from tens to hundreds of meters. While copper (Cu) is the major ore-forming element, gold (Au) and molybdenum (Mo) are also present (Cheng et al., 2009). Extensive studies have identified more than 20 ore-bearing stocks (marked with Roman numerals) (e.g., Shen et al., 2009; Shen et al., 2012; Shen and Pan, 2015; Tang et al., 2010; Zheng et al., 2015; Cao et al., 2016; Ma et al., 2012, etc.). However, limited geophysical data have been collected to delineate the deposits in Baogutu. Motivated by the successful applications of AMT data for various geological targets, we deployed a dense AMT array to image the Baogutu porphyry Cu-Au deposit within Stock V, which we present here as a case study in mineral exploration.

We conducted 11 AMT profiles, with each profile consisting of 16 stations, resulting in a total of 176 sites covering a region of 3 km by 5 km (Figure 2B). The profile spacing is approximately 300 m, and the site spacing is also about 300 m in the central part of the array, increasing to around 400 m for the last 2 and 3 sites at the northern and southern ends of each profile. To minimize noise from a power line that crosses the study area from east to west, Site #10 of each profile was relocated 100 m north or south from its designated position (Figure 2B).

We used the MTU-NET system manufactured by Phoenix Geophysics to acquire the AMT data. For logistical reasons, we only recorded the horizontal electric fields ( E x , E y ) and the horizontal magnetic fields ( H x , H y ), meaning that the tipper data is not available for inversion. Two pairs of non-polarizing electrodes with 50 m arms were employed to measure the electrical field components. The study area is situated in a remote region, far from sources of cultural noise. Field tests showed that a half hour time series was sufficient to estimate the AMT transfer functions. Thus, we recorded both the electrical and magnetic field components for approximately 30 min of each site. We processed the recorded time series using the standard robust remote reference approach (open source package named EMTF, Egbert and Booker, 1986), resulting in reliable estimations of impedance from 10,000 Hz–1 Hz. The data quality is generally good for most sites, characterized by small error bars and low scattering. However, a few sites experienced noisy time series due to poor electrode-grounding contact and interference from the power line. For some sites, the responses in the frequency band of 10,000-1,000 Hz were discarded due to a poor signal-to-noise ratio in the magnetic field between 4,000 and 1,000 Hz, a phenomenon known as the “dead band” of AMT (Garcia and Jones, 2002). Figure 3 shows the observed apparent resistivity and phase of the off-diagonal impedance tensor for three representative stations. Sites 03-13, 06-11, and 09-07, marked by blue squares in Figure 2B, are situated atop the Xibeikulasi Group, at the boundary between the Xibeikulasi Group and the Baogutu Group, and within the Baogutu Group, respectively. These sites exhibit distinct features and good data quality.

We employed the latest parallel 3D MT data inversion scheme, the Modular system for Electromagnetic Inversion (ModEM) (Egbert and Kelbert, 2012; Kelbert et al., 2014), for 3D modeling and inversion. ModEM is widely used for interpreting MT data from continental-scale arrays (e.g., Meqbel et al., 2014; Yang et al., 2015; Bedrosian, 2016; Murphy and Egbert, 2017; Egbert et al., 2022) as well as local profiles (e.g., Xu et al., 2016b; Xu et al., 2016c), due to its power, flexibility, and availability to academic communities. The preferred 3D resistivity model presented in this paper (Figures 6, 7) was obtained by inverting the full impedance tensor from all 176 sites across 59 periods, ranging from 10,400 Hz–0.6 Hz. The topography of the study area is relatively flat (Figure 2C), with elevation variations of less than 100 m throughout the entire region; therefore, we ignored the topographic effect in the inversion process. We tested various a priori models, data error floor settings, model gridding, and other inversion parameters by running the inversion multiple times to recover a reliable inverted model. For the preferred model presented here, we assigned error floors of 5% of | Z x y | for Z x x and Z x y , and 5% of | Z y x | for Z y x and Z y y . Since Z x x and Z x y are estimated from the same electrical-channel, they may be contaminated by similar noise sources in the electric field. Therefore, it makes sense to apply the same error floor to both components. The same rationale applies to the other row of impedance components, Z y x and Z y y . We believe this method of setting error floors is more reasonable, as it accounts for the inherent differences in SNR between diagonal and off-diagonal components while considering potential common noise sources.

We discretized the study area horizontally with 75 m grids in the core, padding with seven grids on all four edges, with increasing widths by a factor of 1.2 outward to the boundaries exponentially, to satisfy the boundary conditions. The model was discretized vertically with 30 layers, starting from 5 m for the first layer, and then the thickness of the layers increased logarithmically to ∼ 20 km for the last layer. We relied upon a very fine parameterization of the uppermost layers to simulate any near-surface electrical resistivity heterogeneities, to accommodate the static shift effects caused by galvanic distortion in the observed data. Static shifts in MT data are primarily caused by static charges concentrating at the edges of shallow electrical heterogeneities. These static charges produce an electric field that does not vary with time and superposes on the varying electric field across the entire frequency band. When an MT site is located near the edge of a shallow anomaly—typically small-scale but not exclusively—the static electric field can either enhance or weaken the original time-varying E-field. This interaction results in changes to the amplitude of the MT impedance, leading to shifts in the apparent resistivity curves. Regarding the handling of static shifts, while it is theoretically possible to invert for the static shift coefficient, our practical experience suggests that parameterizing the near-surface region finely enough to model these effects directly might be a more effective approach. This method allows for a more accurate representation of the complex near-surface structures that contribute to static shifts. This discretization resulted in a 112 × 86 × 30 grid, with seven additional air layers implicitly included in the modeling code.

We used the nonlinear conjugate gradient (NLCG) scheme built into ModEM to fit the data iteratively. Based on the laboratory measurements of the rock samples in the study area (Table 1), the regional background resistivity is relatively large (more than 1,000 Ω ⋅ m). Meanwhile, the median value of apparent resistivity of observed AMT sites is about 600 Ω ⋅ m. Therefore, for the preferred model we started from a 600 Ω ⋅ m half-space a priori model, fitting the data to a normalized root mean square (nRMS) misfit of 2.05 by 72 iterations (Figure 4). Figure 2C shows the nRMS of each site for the preferred model. Generally, the data fit is good for most of the sites, with only a few sites contaminated by noise or with poor data quality in the “dead band” showing slightly higher nRMS (orange and red dots in Figure 2C). We presented the apparent resistivity and phase maps of all four impedance tensor components at two representative frequencies, 320 Hz and 9.4 Hz, in Figure 5. Comparing the observed data and the predicted data (i.e., responses generated by the inverted model), almost all the main anomalies in the observed data have been reproduced in the responses of the inverted model, especially for the off-diagonal components of the impedance tensor. These three representative sites in Figure 3 also show a good fit between the predicted and observed apparent resistivities and impedance phases.

Note that the variation of nRMS with respect to the number of iterations shows fluctuations (Figure 4), particularly in the early stages. Based on our understanding, occasional small increases in nRMS are likely due to the line search scheme used in the NLCG module of ModEM. This scheme evaluates only two additional values of the objective function beyond the current nRMS value. It then uses these three points to interpolate a quadratic parabolic function. The algorithm subsequently selects the minimum of this parabola as the result of the line search. However, if the actual function along the search direction is not well-represented by a parabolic function, the step size determined by this method may be imprecise. This can lead to an increase in nRMS. Specifically, when the true function deviates from a parabolic shape, the chosen step size might not be optimal, resulting in temporary increases in nRMS.

As presented in the selected depth slices (Figure 6) and the vertical cross-sections (Figure 7), the inverted 3D electrical resistivity model clearly reveals the detailed structure of the study area, ranging from dozens of meters to more than 2 km in depth. For convenience, we plotted the outlines of the outcrops V1 to V3 of Baogutu stock V from south to north in Figure 6G, along with the local faults and geological boundaries of the Xibeikulasi Group and Baogutu Group in all subplots of Figure 6. The profiles A-A’ to E-E’ in Figure 6H indicate the locations of the vertical cross-sections shown in Figure 7.

The most pronounced feature in the shallow depths of the model is the east-west conductive zone in the central part of the study area, which is related to the mineralized stock V. We named these conductors C1-C6 (Figure 6C). The shallow layers in the model (Figures 6A, B) indicate that C1, C3, and C4 connect to each other, while C2 does not appear until a depth of nearly 250 m (Figure 6C). In the deeper layers, as shown in the profiles A-A’ and B-B’ (Figure 7), C1 and C2 have merged together and connected to C3 along the east-west direction in the southern area of V2. The C4 (Figures 7D, D’) and C3 (Figures 7E, E’) are clearly separated from C5 by resistive intrusions beneath the outcrops V1 and V2. We have noted that two conductive spots exist in the middle and northern parts of stock V3, as shown in Figure 6D; they appear to be two independent branches. Since the current AMT array cannot constrain these two anomalies precisely, we cannot confidently describe them in detail. The conductive zone in the V1 outcrop area, which coincides with the outcrop boundary (Figures 6A–E), extends vertically from the surface to approximately 700 m in depth.

Dozens of boreholes have been drilled in stock V-2 (e.g., Shen et al., 2009; Shen et al., 2010). Two lithologic profiles have been compiled from the drilling cores of these boreholes (Shen et al., 2010). To verify the conductive anomalies in the preferred model, we plotted the vertical cross-sections alongside the lithologic and mineral profiles revealed by the borehole results (Shen et al., 2010). The lithologic outlines were adapted from Shen et al. (2010). The green polygons in Figure 8B, C indicate that the orebodies are primarily distributed in the upper 600 m of stock V-2, where a highly conductive anomaly (20-40 Ω ⋅ m) corresponds with the copper mineralization. To enhance the interpretation of the inverted model, we also collected rock samples from the Baogutu area and measured their electrical resistivity and shear wave velocity in the laboratory. As shown in Table 1, all samples, whether from volcanic rocks or sediments, are relatively resistive (1,000-4,000 Ω ⋅ m), while the copper orebodies are significantly more conductive (approximately 30 Ω ⋅ m). The two borehole profiles (Figures 8B, C) demonstrate that a dioritic intrusion extends from depth to 400 m, coinciding with the AMT imaging model.

Considering the lack of high-frequency data (10,000–1,000 Hz), which is primarily sensitive to depths of 100–300 m in a 600 Ω ⋅ m half-space, and the limited number of AMT sites (only 4-5) available to constrain the structures near the borehole profiles, the inverted model aligns surprisingly well with the borehole data. This consistency suggests that 3D inversion of AMT data from an array of regularly distributed sites can yield reliable and high-resolution imaging results.

The porphyry Cu-Au-Mo mineralization in West Junggar is associated with intermediate intrusive rocks. Typically, ore-bearing mineralization occurs in the upper stocks and in the contact zones between the stocks and their country rocks (Shen et al., 2009). In our inverted model, several high-conductivity zones can be identified as potential mineralized zones. As noted earlier, all major conductive anomalies are located on the flanks or in the uppermost zones of the stocks.

For stock V2, the orebodies are likely distributed along the southern boundary, exhibiting a shallow east-west extension and a deeper southwest trend. Drilling results in stock V2 indicate that the ore-bearing porphyry extends to approximately 600 m in depth (Shen et al., 2009), which aligns well with our model. Furthermore, C3 connects with C2 at depth and extends more than 1 km, suggesting significant mineralization occurs at depths of 200 m to 1 km along the southern edge of stock V2. The conductor C4 may represent an orebody associated with small faults. As suggested by Richards (2003), large porphyry deposits are likely formed where magma ascents are concentrated at structural intersections. However, conductor C1 might also represent an isolated orebody developed at shallower depths within the Xibeikulasi Group. In all depth slices shallower than 1 km, we can clearly identify an independent high-conductivity area, C6, at the northeast boundary of stock V1. This area is likely an ore-bearing stock extending vertically to 700 m, associated with the intermediate intrusion. For stock V3, conductor C5 likely highlights an orebody that may develop independently on the northeast flank of the stock at depths of 100 m–800 m (Figure 7D).

Shen and Pan (2015) suggested that the mineralized intrusive phases at the Baogutu porphyry belt comprise the main-stage diorite stock and minor late-stage diorite porphyry dikes. The main-stage host experienced significant hydrothermal alteration and contains the majority of the Cu-Mo-Au mineralization. This implies that the presence of high-conductivity zones near stock V is likely associated with the Cu-Au deposit.

Xu H. et al. (2016) employed the eikonal tomography method to map the shallow 3D velocity structure of the Baogutu V stock belt. They noted that the Xibeikulasi Group (C1x) in the west exhibits a higher shear wave velocity compared to the Baogutu Group (C1b) in the east. However, our model does not reveal any significant resistivity differences between these two groups, suggesting that they likely have similar conductivity. To facilitate comparison with the shear wave velocity model, we plotted two depth slices (at 450 m and 700 m) in Figures 9A, B, using the same map range as in Figure 6. In Figure 9A, the low shear wave velocity anomaly at the southern edge of stock V2 aligns closely with conductor C3. Additionally, measurements of rock samples indicate that ore-bearing intrusive rocks (e.g., gabbro, dioritic-porphyrite) possess low shear wave velocities and relatively low resistivity, especially when mineralized, suggesting that C3 could be a promising mineralized zone. To examine the overall correlation between the shear wave velocity (Vs.) model and the resistivity model, we computed their Pearson’s correlation coefficient by first interpolating them onto the same grid and then calculating the correlation coefficient point by point with respect to the depth variations of the models using following Equation 1 R σ , V = 1 N − 1 ∑ i = 1 N σ i − μ σ δ σ V i − μ V δ V (1) where R is the correlation coefficient, σ i , μ σ and δ σ are the electrical conductivity, its mean and standard deviation respectively, V i , μ V and δ V are the shear wave velocity, its mean and standard deviation.

It is important to note that the Vs. model is best resolved only in the central part of the array (see Xu et al., 2016a, for details), so we limited our correlation analysis to this region. As shown in Figure 9C, the resistivity model strongly positively correlates with Vs. at the southern edge of stock V2, while it shows a strong negative correlation with Vs. at stock V1. This suggests that stock V1 may have different lithological properties compared to stock V2. The strong positive correlation between S-wave velocity (Vs.) and conductivity observed in Stock V1 and V2 likely indicates the presence of a fracture zone associated with a local fault. Additionally, a prominent negative correlation belt surrounding Stock V1 is one of the most notable features. This area’s low conductivity and high Vs. values may suggest the formation of a fracture contact zone during the intrusion of the magmatic stock.

Regarding the metallogenesis of the Baogutu porphyry copper belt, Zhang et al. (2010) suggested that the Baogutu porphyry copper deposit may have formed in the Late Carboniferous, with its mineralization closely linked to the intrusion of intermediate-acidic porphyry bodies. They posited that the ore-forming fluid originates from deep magmatic waters. Tang et al. (2010) and Yin et al. (2010) proposed that the Baogutu Cu-Au deposits could occur above a slab window during ocean ridge subduction, where the interaction between high-oxygen fugacity slab melt and mantle peridotite leads to the decomposition of metal sulfides, subsequently promoting Cu and Au mineralization. Furthermore, Shen et al. (2009) suggested that the ore-bearing porphyry system derives from the partial melting of multiple sources, including oceanic crust and a subduction-modified mantle wedge, with significant crystal fractionation occurring during the convergence between the paleo-Junggar ocean and the Darbut arc. Given the limitations of our AMT data, which lacks long-period information and covers only a relatively small aperture, we are unable to resolve any structures deeper than 2 km. This restriction complicates discussions surrounding the metallogenesis of the Baogutu porphyry copper deposit in our study area. However, based on the tectonic evolution inferred from the regional-scale resistivity model (Xu et al., 2016b), we support the view that the Baogutu porphyry copper belt is likely related to an intra-oceanic subduction event beneath the Darbut arc during the Late Paleozoic.

We imaged the high-resolution resistivity structure of stock V of the Baogutu intermediate-acidic intrusions using an AMT dataset consisting of 176 sites. The Baogutu porphyry copper belt is distinctly represented as conductive zones along the flanks of outcrops V1 to V3, particularly at the northeast flanks of V1 and V3, and extending in an east-west direction between V1 and V2. Our model, inverted from the AMT data, reveals that all significant high electrical conductivity zones extend from the surface to approximately 600 m in depth. When integrated with previous geological, geochemical, and geophysical studies, our model provides insights into the 3D structure of the Cu-Mo-Au mineralization zones associated with the magmatic activities of the Baogutu intrusions. Most mineralization appears to develop irregularly within the ore-bearing stocks and at the contact zones between the stocks and their country rocks, suggesting that the porphyry copper deposits may be influenced by small intermediate-acidic intrusions rather than showing significant correlations with existing faults. Comparing our results with borehole profiles, the AMT imaging aligns closely with the lithological profiles compiled from drilling cores. The resistivity model also demonstrates a strong correlation with the shear wave velocity model. We recommend that further drilling be conducted at other promising conductive anomalies. Additionally, to enhance our understanding of metallogenesis, the current array should be expanded to include broadband MT stations.