3D Numerical Simulation-Based Targeting of Skarn Type Mineralization within the Xuancheng-Magushan Orefield, Middle-Lower Yangtze Metallogenic Belt, China


 Recent exploration has identified a series of Cu-Mo skarn deposits within the Xuancheng-Magushan orefield. The orefield forms part of the Nanling-Xuancheng mining district, which is located within the Middle-Lower Yangtze River Metallogenic Belt (MLYRMB) of central-eastern China. However, this area contains thick and widespread unprospective sedimentary cover sequences that have impeded traditional approaches to mineral exploration. This study presents the results of 3D numerical simulation modeling that identifies possible mineral exploration targets within the entire Xuancheng-Magushan orefield. This modeling enables the identification of unexplored areas with significant exploration potential that are covered by thick sedimentary sequences that cannot be easily explored using traditional exploration approaches. This study outlines the practical value of 3D numerical simulation-based targeting in areas with thick sedimentary cover sequences and uses the Flac3D software package to couple processes involved in ore formation such as stress, pressure, and heat transfer. Here, we use volumetric strain increments calculated during numerical modeling as the thermodynamic representation of the generation of space during prograde skarn formation, with this space filed by sulfides either penecontemporaneously or soon after magmatism. This process occurred during retrograde hydrothermal ore formation and the genesis of the skarn-type mineralization in this area. The results of the volumetric strain increment calculated during this numerical modeling study matches the distribution of known mineralization as well as delineating eight potential targets that have not yet been explored but represent areas of significant exploration potential within the Xuancheng-Magushan orefield, indicating these targets should be considered prospective for future mineral exploration. One of these targets was also identified during our previous Comsol-based numerical modeling of the formation of the Magushan Cu-Mo skarn deposit. The fact that this area has been identified as prospective using two different numerical modeling methods indicates that this area should be prioritized for future exploration and also validates the numerical modeling approaches used here and in our previous research that more specifically focused on the Magushan skarn deposit. Overall, our study indicates that prospectivity modeling using 3D numerical simulation-based approaches can be both effective and economical and should be considered an additional tool for future mineral exploration to reduce exploration risks when targeting mineralization in areas with thick and unprospective sedimentary cover sequences.


Introduction
Skarn deposits are widely distributed in both space and time (e.g., [1]) and are important sources of copper, lead, zinc, molybdenum, and other metals [2][3][4]. These mineral deposits are dominated by calc-silicate minerals, such as garnet and pyroxene, which are generated during the prograde development of skarns, leading to the presence of these minerals and being key in terms of the definition and identification of skarns (e.g., [1]). The Middle-Lower Yangtze River Metallogenic Belt (MLYRMB) is located in central-eastern China and contains numerous porphyry-skarn copper-gold and magnetite-apatite iron deposits that formed in an intracontinental setting [5][6][7][8]. The Xuancheng-Magushan orefield is located in the southeastern part of the Nanling-Xuancheng mining district, a newly defined mining camp within the MLYRMB [7,[9][10][11]. Most of the orefield is covered by a thick and barren sedimentary cover sequence that means very few of the prospective lithological units in this area are exposed, making exploration using traditional approaches (e.g., surface geochemistry or potential field geophysics) difficult.
Rapid recent advancements in computing hardware and theoretical and applied computational science have enabled the development of three-dimensional prospectivity modeling (e.g., [12][13][14][15][16][17][18]) and the complex coupled numerical simulation of physical processes [19][20][21][22][23], including those involved in ore formation and mineralizing systems. However, traditional three-dimensional prospectivity modeling approaches (e.g., weights-of-evidence) frequently encounter issues related to conditional independence, where datasets are biased by relationships where given exploration criteria can generate responses in different datasets, particular mineralizing processes can generate more than one exploration criteria, or responses present within one dataset can be conditioned by responses in another dataset [15]. This can be overcome by using the numerical simulation approaches utilized in this study, providing another approach for the identification of areas that are prospective for deeper exploration in regions with thick sedimentary cover sequences. This approach therefore acts as a complement to more traditional exploration methods and could potentially remove risks related to the verification of conditional independence involved in other types of modeling (e.g., Weights of Evidence; [21,[24][25][26]). This study uses Itasca Flac 3D , a professional numerical simulation software package that is typically used for the modeling of multiple physical processes in rock and soil engineering. This means that even though this software is not specifically designed for geological research, it allows the coupling of multiple physical fields (e.g., temperature and pressure) that enables the easier formulation of coupled simulation models. This approach allows the numerical simulation of the hydrothermal and mineralizing processes that generated the skarn deposits in the study area (as well as other types of mineralization). In addition, this type of modeling and simulation can highlight unexplored areas where hydrothermal and mineralizing processes are likely to have occurred but where mineralization has not yet been identified or exploration has not as yet taken place, yielding targets for future exploration.
This study utilizes 14 geophysical data-based inversion cross-sections constructed using gravity and magnetic data for the Xuancheng-Magushan orefield to build a 3D model. This model is then combined with a numerical simulation approach to generate a 3D simulated model that can be used to identify areas for future exploration within the Xuancheng-Magushan orefield. We also address the following questions using numerical simulations, namely, (1) Can numerical simulations enable the identification or verification of existing geological data? and (2) Can numerical simulations be used to identify prospective areas for exploration targeting in areas with thick and unprospective sedimentary cover sequences, such as the Xuancheng-Magushan orefield, as well as providing guides for future mineral exploration? Both of these points are assessed during this study, demonstrating that numerical simulation can be used in both fundamental geological and economic geology research as well as in exploration targeting in regions of orefield scale with thick and unprospective cover sequences.
The Nanling-Xuancheng basin is located within the eastern MLYRMB ( Figure 1) and is divided into upper and deeper sections by an unconformity, changes in deformation and magmatic histories, and the sedimentary stratigraphy of the basin. The deeper part of the basin contains conformable or disconformable Silurian (O) to lower and middle Triassic (T 1-2 ) sedimentary units. The Indosinian orogeny associated with the collision between India and Eurasia generated intense NE-SW-oriented folding in these units, which are also cross-cut by NW-SE-and NE-SW-oriented faults. The upper part of the basin contains Jurassic (J) to Quaternary (Q 4 ) units, including the Jurassic volcanic Zhongfencun Formation (J 3 z). The majority of these formations have unconformable contacts, are variably folded, contain interbedded volcanic units, and have been intruded by magmatic units. The majority of the mining camp area is covered by thick Quaternary sediments although this has not prevented the discovery of a number of polymetallic Cu, Au, and Mo deposits as evident in Figure 2. These include the Chating porphyry Cu-Au and the Shizishan Cu, Qiaomaishan Cu-S, Fenghuangshan Cu-Mo, and Magushan Cu-Mo skarn deposits (e.g., [34]), with the majority of the deposits in this 2 Lithosphere area located within the Liqiao-Tongshan and Xuancheng-Magushan orefields ( Figure 2; [35]).

Orefield
Geology. The Xuancheng-Magushan orefield is located in the southeast part of the Nanling-Xuancheng basin ( Figure 2) and contains NW-SE-and NE-SW-oriented groups of faults that form a tectonic framework when combined with the intense folding that is present across the northeast and southeast parts of the basin. The whole area is dominated by Triassic and Cretaceous units that are associated with the mineralization in this area (Figures 3 and 4) and form the host rocks for the skarn-related intrusions that in turn host the orebodies that define the skarn-type Fenghuangshan and Magushan Cu-Mo deposits.
2.3. Deposit Geology. The Magushan Cu-Mo deposit is a skarn copper-molybdenum deposit that is located in the northeast part of the Xuancheng-Magushan orefield (Figures 3 and 5(a)). It hosts resources and reserves containing an estimated 78,000 t of contained Cu, 11,000 t of contained Mo, and 1800 t of contained Au [36]. This area contains Cretaceous to Quaternary sediments, with the latter cropping out over more than 75% of the region (Figure 3).
Drillholes in this area indicate that the Quaternary sediments in this region cover Triassic, Jurassic, Carboniferous, and Silurian units that are hosted by two NE-SW synclines located within the northeastern and southwestern parts of this region (Table 1; Figure 4). The intrusions in this area appear to be spatially controlled by the presence of the Permian units within the southeast limb of the synclines as these units are structurally weak, meaning that the magmas that formed these intrusions exploited these weak points during their emplacement ( Figure 5(b)). The majority of the Cu and Mo mineralization within the Magushan deposit is hosted at the contact between a porphyritic granodiorite intrusion with mixed mantle-crustal origins and surrounding sedimentary units, especially the Permian and Carboniferous limestones in this area that appear to have reacted with exsolved magmatic fluids to form the skarn mineralization within the deposit [37]. The area also contains significant contact metamorphism associated with the intrusions present in this region. The Magushan deposit itself consists of one main orebody, one subsidiary orebody, and a number of smaller orebodies. The molybdenite and chalcopyrite mineralization in this area is associated with a skarn containing concentric zones of diopside and garnet bearing assemblages. The majority of the orefield is covered by thick and unprospective Cretaceous, Eocene, and Quaternary sediments (Table 1) meaning that traditional geochemical and geological exploration techniques have limited use with the exception of drilling. This led to the use of geophysical exploration using gravity and magnetic potential fields approaches [37] and the generation of a joint inversion constrained by a geological map and drillholes data for this area that yielded a number of interpreted geological cross-sections (Figures 3 and 4). All 14 sections are constrained by geological mapping of the Xuancheng-Magushan orefield ( Figure 3) and associated logging of drillcores in the area around the Magushan deposit ( Figure 5).

Modeling of the Xuancheng-Magushan Orefield
A total of 14 cross-sections that cover the entirety of the Xuancheng-Magushan orefield were constructed based on the constrained inversion of existing airborne gravity and magnetic data (Figures 3 and 4). These sections are used to define the basic 3D form of the model for the entire study area with Intrepid Geomodeller software used during entirety of the 3D modeling process. The workflow for this paper reflects all of the 3D modeling and numerical simula-tion processes used during this study and is shown in Figure 6.
3.1. 3D Modeling of the Research Area. There are three steps involved in the 3D modeling using the Intrepid Geomodel-ler™ software package. The first of these is the importing of source data (e.g., interpreted cross-sections, geological maps showing the locations of different geological units, structures, and mineralization, and other related data) into the software, ensuring real coordinates match with the registered images within the software system. The second step is the digitization of all of the geological elements within the model (e.g., boundaries and contacts between sedimentary units, between country rocks and intrusions, and the locations and shapes of faults). The final step is to define the relative temporal and spatial relationships between all of the different elements within the model (e.g., the relationship between different faults, sediments, and intrusions within the study area). The only remaining step is to undertake the computation and the initiation of the automatic computing and building process. The resulting model constructed using the Intrepid Geomodeller software package is shown in Figure 7, representing the initial 3D model generated during this study. Each geological unit constructed during this modeling was then assigned their associated gravity and magnetic parameters determined using the geophysical data available for the  4 Lithosphere study area. This quantification allowed forward modeling calculations that enable the modification of the initial 3D model to minimize the differences between the calculated 3D geological model and the real geological structure in the study area, a standard approach for forward modeling ( Figure 8). The final models for the individual geological units within the study area are shown in Figure 9.
3.2. Data Transformation. The Flac 3D software package incorporates an approach to the generation of Flacformat 3D block models that differs from the voxel format 3D block models used in the Geomodeller software. This means that the voxel models exported from the latter cannot be used in the former without applying a transformation despite the fact that both software packages use basic components in the form of blocks. This transformation requires the coding outlined in Figure 10 and Tables 2 and 3.
If we assume the model contains a block m that has a line number of R m in the voxel data file, then the location of this block can easily be calculated in 3D space along axes i, j, and k (i.e., the block is the number i block along the x-axis, number j block along the y-axis, and number k block along the z -axis). Thus, the m0, m1, m2, m3, m4, m5, m6, and m7 values are used to represent the 8 vertices of the block yielding point numbers of n m0 , n m1 , n m2 , n m3 , n m4 , n m5 , n m6 , and n m7 . The values of i, j, and k can be directly calculated using the following equations: and the coordinates of all vertices can be therefore calculated using the following equations: Finally, vertex numbers can be calculated using the following equations: These equations allow the transformation of the voxelbased Geomodeller block model for the Xuancheng-Magushan orefield into the block model system used by the Flac 3D software, with the resulting transformed model shown in Figure 11 (where each block is a 500 × 500 × 500 m cube).

Conditional Settings and Mathematic
Modeling of the Magushan Area. Previous research undertaken within the   Lithosphere study area [24] suggests that the formation of the Magushan skarn deposit took place at a deeper location than the current depth of mineralization. This change in elevation is addressed in our model by the addition of a 1.3 km layer (actually 1.5 km due to the transformation of data) that represents Cretaceous sediments that were denuded and lost after the formation of the mineralization in this region ( Figure 12). The initial temperature of the intrusion is set as 700°C with the starting temperature of the sediments in this region changing with depth as determined using a typical geothermal gradient (25°C/km; [38,39]) and a surface temperature that is set at room temperature (normally 20°C at 1 atmospheric pressure, i.e., 0.1013 MPa). The pressure within the Flac 3D model is automatically calculated using density data and the volume of the imported model, indicating there is no requirement to set the pressure gradient for this model. The numerical simulation undertaken during this study involves coupling fluid flow and mechanical and heat transfer processes. The intrusion in the study area and the surrounding sedimentary units are set as homogenous and isotropic material with porous viscoelastic characteristics. The mechanical behavior of these units is represented using a Mohr-Coulomb model with the main associated processes and their formulas and equations given below [21,40].
(1) Fluid-flow Fluid-flow in this model is governed by Darcy's law [41], where q W is the seepage velocity of fluid, k is the coefficient of the flow, P is the pore pressure, ρ w is the density of fluid, and g is the gravity constant.
(2) Heat-transfer In Equation (5), q T i is the thermal flux vector, k T is the effective thermal conductivity, and T i is temperature.
(3) Conservation of mass In Equation (6), ξ is the change of fluid capacity, q i,i is the   (4) Conservation of energy In Equation (7), C T is the effective heat ratio and c w is the heat capacity of the fluid.

(5) Conservation of momentum
In Equation (8), v i is the velocity component along the x -axis, and σ ij,j is the stress tensor.
In Equations (9) and (10), ε ij is the thermal strain tensor, α t is the linear coefficient of thermal expansion, δ ij is the Kronecker operator, M is the Biot modulus, α is the Biot coefficient, β is the coefficient of thermal expansion of porous medium, and ε is the volumetric strain increment.
The related mechanical parameters of the rock units within the Xuancheng-Magushan orefield are given in Table 4 with their thermal parameters given in Table 5.

Results
Volumetric strain increment (VSI) values can be thought of as the thermodynamic representation of ore-forming processes of skarn deposit. The VSI values identified in numerical simulations represent the modeling of rock fractures caused by physical and chemical processes such as thermal expansion, thermal contraction, stress, strain, and changes in mass and energy [21,26,42]; (Zhao et al., 2011). The skarn alteration in the study area involved the initial prograde formation of garnet and other calc-silicate minerals from the calcite within the hosting limestones, increasing the density and decreasing the volume of the resulting skarn [43]. This process also generated increased volumetric strain values that can be identified during our modeling as well as increasing the porosity and permeability of the skarn. This increase in porosity and permeability allowed fluid flow during the retrograde phase of skarn formation (e.g., [43]), leading to the filling of the resulting voids by the precipitation of spacefilling chalcopyrite and molybdenite during the genesis of these Cu-Mo skarns.
Previous research on skarn deposits indicates that areas with a certain range of VSI values (generally between 0.1% and 1%) represent areas prospective for hydrothermal mineralization, meaning that the distribution of VSI values can be used to identify targets for future mineral exploration [26,42]; (Zhao et al., 2011). This approach has been successfully used in other 3D prospectivity modeling and associated economic geology-focused research (e.g., [26,42]; Zhao et al., 2011). This link reflects volumetric changes associated with variations in internal inner stress and temperature; both of which are thought to be associated with the mineralizing  9 Lithosphere processes that operate in skarn systems [43]. Combining this with the fact that the garnet skarns within the study area are the dominant hosts for the known Cu-Mo mineralization in this area [36] means that the distribution of VSI values within our model can be used to identify areas prospective for exploration for skarn-type mineralization. Here, VSI values above 11 Lithosphere zero represent increases in density, porosity, and permeability; each of which provided both conditions suitable for the flow of ore-forming fluids as well as space that could be filled by sulfide mineralization. In comparison, VSI values below zero represents a decrease in density, porosity, and permeability, reducing both fluid flow and mineralizing potential, indicating these areas should be considered unprospective for mineral exploration.    Figure 12: The transformed 3D model before (a) and after (b) the addition of the Cretaceous sediments that were denuded and lost between mineral deposit formation and the present day.

Lithosphere
The results of our modeling are visualized by the generation of depth slices; essentially a series of plan-view maps constructed at 500 m depth intervals for the entirety of the model (Figures 13 and 14). The map for the -500 m depth interval indicates a cluster of high VSI values that are proxi-mal to the locations of the Fenghuangshan deposit, indicating a relationship between high VSI values and areas of known mineralization (Figure 15). The VSI data at this depth also indicate the unexplored area surrounding the Fenghuangshan deposit should be considered prospective. There  14 Lithosphere is no known mineralization at depths of -1000, -1500, and -2000 m in the study area, meaning that the high VSI regions at these depths may represent concealed deep-seated skarn mineralization that is yet to be discovered. Finally, there are no anomalously high VSI values at depths between -2500 and -4500 m, reflecting either a lack of mineralization or poor data at these extreme depths.
The VSI data outlined above enabled the identification of eight possible targets for future exploration for skarn-type mineralization in the study area ( Figure 16). Of these, we 15 Lithosphere consider target 2 to be the highest priority as this target is located in the same area and at the same depth as a target identified during the Comsol-based numerical modeling of ore-forming processes that focused on the area immediately surrounding the Magushan deposit [24]. This identification of the same target using very different mathematical methods and numerical simulation approaches strongly supports both the modeling approaches used in this study and the high prospectivity of this target. In addition, targets 3, 4, and 5 are considered highly prospective as a result of the colocation of these targets and gravity anomalies that may reflect the location of mineralization, intrusions, or alteration that change the density and gravity characteristics of this area ( Figure 15). This study also identified targets 3, 4, and 5; all of which appear to have formed at depths of -1.8 km (present day depth of -500 m combined with 1.3 km of denuded sediments, consistent with the known depths of skarn formation (e.g., [1]), further validating the approaches undertaken during this study and the prospectivity of these targets. It is important to note that although these targets originally formed at a depth of -1.8 km, their present day depth of -500 m is much more amenable for both exploration and economic exploitation.

Discussion
Mineral exploration in regions containing thick and unprospective cover sequences requiring the targeting of deepseated mineralization remains challenging, especially over time as near-surface targets are identified and brought into production, a process that requires exploration to target deeper levels over time. The deep nature of this exploration means that only certain types of exploration tools can be used, such as geophysics and drilling based on structural and geophysical targeting [44][45][46][47][48]. Recent advances in 3D prospectivity modeling have highlighted the potential use of prospectivity [12][13][14][15][16] and numerical modeling [21] in exploration targeting. However, the simulation undertaken during this study provides an additional approach to identifying areas prospective for deep-seated exploration purely based on 3D numerical simulation, supplementing (although not replacing given the reliance of this approach on, e.g., geophysical data) more traditional exploration approaches and removing some of the risks conditional dependence-related risks associated with other types of prospectivity modeling.
The modeling undertaken during this study has identified areas of known mineralization as prospective for mineral  16 Lithosphere exploration despite this modeling being independent of any training dataset or inputs indicating areas of previously known mineralization. This independence indicates that simulation-based approaches could provide independent verification of areas thought to be highly prospective targets for future exploration as well as having potential use in its own right for exploration for deep-seated mineralization both in the study area and elsewhere. The results of this study are consistent with research by Li et al. [21], who used a similar numerical simulation approach to undertake a 3D prospectivity analysis of the skarn-dominated Yueshan orefield in China. Both the research undertaken during this study and the research presented by Li et al. [21] indicate that the incorporation of numerical simulation can enhance the results of prospectivity modeling, again supporting the usefulness of this approach in mineral exploration, especially in areas with significant cover or where mineralization is located at depth.
This study used numerical simulation to identify areas for future exploration in a large area containing two known skarn deposits. The results of this study indicate that this approach, ideally when combined with modeling of mineralizing processes on a more localized scale [24], can be an effective and economical addition to traditional exploration approaches, ensuring that most can be made of the data available for an area. The Xuancheng-Magushan orefield is located within the southern part of the newly discovered Nanling-Xuancheng mining district, an area with high exploration potential. This significant potential is exemplified by the recent discovery and exploitation of the Magushan and Fenghuangshan skarn deposits [36]. This approach is far more economical than the significant amount of grid-type drilling that would be needed to test the entire area, which covers some 23 × 21 km, with an approach combining our modeling with selected drilling potentially proving a much more effective approach to exploration. The situation in the 17 Lithosphere study area emphasizes the increasing challenges facing exploration globally, not just within this region, and reflects a gradual increase of exploration depth as the majority of outcropping and near-surface mineralization has been identified and exploited over time. This in turn suggests that mathematical/numerical approaches to mineral exploration will become increasingly important in the near future, with this type of numerical simulation is likely to be more widely applied during future mineral exploration.
However, there are still some minor shortcomings in our research. The relationship between volumetric strain increment and ore-forming potential values is nonlinear, which means it is impossible to evaluate differences in potential between targets by merely comparing VSI values. Locations with maximum VSI values are obviously more prospective than the surrounding nonanomalous areas with more background-like VSI values, but it is not possible to compare areas with anomalously high VSI values to determine which individual target is more prospective. For example, saying a target with a VSI of 1% is 10 times more prospective than a target with a VSI value of 0.1% is incorrect as the relationship between these VSI values is not linear. In addition, our numerical simulation approach is somewhat simplified in that it only involves mechanics and heat-transfer processes whereas real mineralizing systems are far more complex and involve a variety of different interacting processes and variables such as mechanics, heat transfer, fluid flow, chemical reaction, and material migration. This means that simple simulations combining just mechanics and heat transferal cannot model all of the details of a mineralizing system, although equally the more complex a model, the more uncertainty is involved, and the more computing power is needed, meaning a trade-off between simplicity and essentially cost must take place. It will be increasingly possible to develop more complex and fully coupled 3D numerical models as computing power develops and as exploration sites become more data rich in terms of the generation of higher resolution geophysical and geological data as well as at-rig geochemistry. This means that future models may well provide more insight into the geology and mineralization in areas covered by thick sedimentary sequences, enabling the delineation of more accurate targets during future mineral exploration. However, for now, the modeling presented in this study clearly indicates that the concept of this numerical modeling approach to mineral exploration is valid.

Conclusions
This study uses a 3D numerical simulation method to analyze the Xuancheng-Magushan orefield, an area where prospective units are frequently covered by barren and thick sedimentary sequences, and has produced the following key findings: (1) The simulated distribution of volumetric strain increment values matches areas of a known skarn deposit identified during drilling and exploitation of the deposits in this area to date, indicating that that our modeling could potentially be used to predict areas of hitherto unknown mineralization in areas with thick and barren sedimentary cover (2) Our modeling also identified eight prospective exploration targets within the entirety of the Xuancheng-Magushan orefield. One of these targets (target 2) matches the location of the results of a more locally focused simulation model constructed using Comsol software and utilizing a finite element method. The other three targets (targets 3, 4, and 5) identified during this study match the locations of gravity anomalies that may indicate intrusions or alterations associated with mineralization. The gravity data for the area also includes a large low gravity anomaly that cannot effectively be used for exploration targeting (Figure 8). This indicates the value provided by our numerical modeling, where the much smaller area outlined as prospective around target 8 in Figure 8 providing a more selective basis for future mineral exploration relative to the larger area covered by the associated gravity low. This indicates that 3D numerical simulation could be an effective approach to 3D prospectivity modeling and exploration targeting in areas with thick and barren sedimentary cover sequences where traditional geological and geochemical techniques are less effective or have significant costs, with this modeling representing a cost-effective add on to ensure the maximum use of the available data

Data Availability
Joint inversion sections and parameters for numerical modeling can be found in the main manuscript. The gravity and magnetic anomaly geophysical source data used in this study are confidential. As such, coordinates and absolute anomaly values have been removed but relative values are illustrated within the manuscript. Lithosphere