Numerical Simulation Based Targeting of the Magushan Skarn Cu–Mo Deposit, Middle-Lower Yangtze Metallogenic Belt, China

The Magushan Cu–Mo deposit is a skarn deposit within the Nanling–Xuancheng mining district of the Middle-Lower Yangtze River Metallogenic Belt (MLYRMB), China. This study presents the results of a new numerical simulation that models the ore-forming processes that generated the Magushan deposit and enables the identification of unexplored areas that have significant exploration potential under areas covered by thick sedimentary sequences that cannot be easily explored using traditional methods. This study outlines the practical value of numerical simulation in determining the processes that operate during mineral deposit formation and how this knowledge can be used to enhance exploration targeting in areas of known mineralization. Our simulation also links multiple subdisciplines such as heat transfer, pressure, fluid flow, chemical reactions, and material migration. Our simulation allows the modeling of the formation and distribution of garnet, a gangue mineral commonly found within skarn deposits (including within the Magushan deposit). The modeled distribution of garnet matches the distribution of known mineralization as well as delineating areas that may well contain high garnet abundances within and around a concealed intrusion, indicating this area should be considered a prospective target during future mineral exploration. Overall, our study indicates that this type of numerical simulation-based approach to prospectivity modeling is 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, and represent important sources of copper, lead, zinc, molybdenum, and other metals [1][2][3].Skarns are dominated by calc-silicate minerals such as garnet and pyroxene, and the presences of these minerals are also key in terms of the definition and identification of skarns [4].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 Magushan Cu-Mo deposit is a skarn deposit in the southeast of the Nanling-Xuancheng area, one of the eight mining camps within the MLYRMB [8][9][10][11].This area is covered by a thick and barren sedimentary cover sequence with little exposure of prospective lithologies, making traditional (e.g., geochemical, geological mapping, some geophysics) exploration difficult.
Rapid recent advancements in computing hardware and theoretical and applied computational science have enabled the development of prospectivity modeling [12][13][14][15][16] and the complex coupled numerical simulation of geological processes [17][18][19][20][21].However, typical three-dimensional prospectivity modeling workflows (e.g., weights-of-evidence approaches) frequently encounter issues of conditional dependence, 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 [12].This can be overcome by using the simulation approaches utilized during this study, providing another method of identifying areas that are prospective for deep exploration in regions with barren cover sequences, providing a supplement to more traditional exploration approaches and potentially removing risks related to the verification of conditional independence.This study combines this approach with Comsol Multiphysics, a software package that allows researchers to customize distributed ordinary differential equations and formulas to calculate variables that allow the simulation of complex and coupled processes that occur within hydrothermal systems.This means that although this software is not specifically designed for geological research, it allows the coupling of multiple physical and chemical fields, allowing the easier formulation of coupled simulation models.This approach allows the numerical simulation of the hydrothermal and mineralizing processes that generated the Magushan Cu-Mo deposit (as well as other types of mineralization) in addition to highlighting unexplored areas where these processes are likely to have occurred, yielding targets for future exploration and additional data generation.
The Magushan skarn Cu-Mo deposit was discovered in the 1960s and mining of the deposit began in the 1990s before the mine shut shortly after as a result of a lack of exploration in this region and the somewhat poor understanding of the system.A review of resources and reserves in 2003 indicated that this deposit still contained a significant amount of copper and molybdenum, leading to the recommencement of mining in 2010 [22].The deposit has current resources and reserves containing an estimated 78,000 t of contained Cu, 11,000 t of contained Mo, and 1800 t of contained Au, and the acquisition of new geophysical data in 2016 also suggested the presence of potentially mineralized intrusions beneath the thick sedimentary cover sequences that surround the Magushan ore deposit [23].Here, we use these and other data for the Magushan deposit and the surrounding region to present a new conceptual model for the mineralization in this area.We also address the following questions using new numerical simulations combined with the results of previous geological and geophysical research, namely: (1) Can numerical simulations enable the identification or verification of existing geological data?; (2) Can numerical simulations be used to identify prospective areas for exploration targeting in areas with thick and unprospective sedimentary cover sequences, such as those around the known Magushan deposit, and provide a guide 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 with thick and unprospective cover sequences.

Regional Geology
The MLYRMB (Figure 1) hosts world-class Cu-Fe polymetallic mineralization and is one of the most important areas of mineralization in China [5,6,11,[24][25][26].The area is cross-cut by a series of major faults and hosts eight large mining camps, namely (from west to east) the southeast Hubei, Jiurui, Anqing-Guichi, Luzong, Tongling, Nanling-Xuancheng, Ningwu, and Ningzhen mining districts Minerals 2019, 9, 588 3 of 19 (Figure 1), [5,11,25,[27][28][29].The mineral deposits in this region host 13 Mt of contained Cu and 800 t of contained Ag, the majority of which are present in porphyry, skarn, and epithermal deposits.All of these deposits are genetically associated with Mesozoic magmatism spread over three separate stages at 149-135, 133-125, and 123-105 Ma [7] and the evolution of this area during the associated Yanshanian tectonic event [30].and 800 t of contained Ag, the majority of which are present in porphyry, skarn, and epithermal deposits.All of these deposits are genetically associated with Mesozoic magmatism spread over three separate stages at 149-135, 133-125, and 123-105 Ma [7] and the evolution of this area during the associated Yanshanian tectonic event [30].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 (T1-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 (Q4) units, including the Jurassic volcanic Zhongfencun Formation (J3z).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.These include the Chating porphyry Cu-Au and the Shizishan Cu, Qiaomaishan Cu-S, and Magushan Cu-Mo skarn deposits, with the majority of the deposits in this area located within the Chating and Magushan orefields [31].

Orefield Geology
The 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 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.These include the Chating porphyry Cu-Au and the Shizishan Cu, Qiaomaishan Cu-S, and Magushan Cu-Mo skarn deposits, with the majority of the deposits in this area located within the Chating and Magushan orefields [31].

Orefield Geology
The 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 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 Magushan Cu-Mo deposit (Figure 5).
basin.The 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 Magushan Cu-Mo deposit (Figure 5).basin.The 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 Magushan Cu-Mo deposit (Figure 5).

Deposit Geology
The Magushan Cu-Mo deposit is a skarn copper-molybdenum deposit located in the northeast part of the Magushan orefield (Figures 3 and 5a) and 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 [22].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 5b).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 the surrounding sedimentary units, especially the Permian and Carboniferous limestones in this area that reacted with exsolved magmatic fluids to form the skarn mineralization within the deposit [23].The area also contains significant contact metamorphism associated with the intrusions present in this region.The Magushan deposit consists of one main orebody, one subsidiary orebody, and a number of smaller orebodies.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 beyond e.g.drilling.This led to geophysical exploration using gravity and magnetic potential fields approaches [23] and the generation of a constrained joint inversion for this area that yielded a number of interpreted geological cross-sections (e.g.A-A' in Figure 3).This section passes through the southwest part of the Magushan deposit zone as well as the NE-SW trending syncline, the southeastern limb of which hosts the Magushan skarn.This interpreted section also indicates that another intrusive apophyse is likely to be present within the northwest limb, with this intrusion also potentially being mineralized (Figure 4).

Deposit Geology
The Magushan Cu-Mo deposit is a skarn copper-molybdenum deposit located in the northeast part of the Magushan orefield (Figures 3 and 5a) and 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 [22].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 5b).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 the surrounding sedimentary units, especially the Permian and Carboniferous limestones in this area that reacted with exsolved magmatic fluids to form the skarn mineralization within the deposit [23].The area also contains significant contact metamorphism associated with the intrusions present in this region.The Magushan deposit consists of one main orebody, one subsidiary orebody, and a number of smaller orebodies.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 beyond e.g., drilling.This led to geophysical exploration using gravity and magnetic potential fields approaches [23] and the generation of a constrained joint inversion for this area that yielded a number of interpreted geological cross-sections (e.g., A-A' in Figure 3).This section passes through the southwest part of the Magushan deposit zone as well as the NE-SW trending syncline, the southeastern limb of which hosts the Magushan skarn.This interpreted section also indicates that another intrusive apophyse is likely to be present within the northwest limb, with this intrusion also potentially being mineralized (Figure 4).

Modeling of the Magushan Skarn Deposit
There are three steps in the computation of geoscience workflows, namely geological conceptual, simulation, and mathematical modeling [32][33][34].Here we construct conceptual geological, mathematical, and simulation models for the Magushan Cu-Mo deposit to both advance our knowledge of the processes that control the formation of skarn deposits as well to determine whether numerical simulations can be used to identify prospective areas in regions with thick and barren cover sequences (such as those around the known Magushan deposit) and as a guide for future exploration targeting.The workflow used in this paper to reflect the processes involved in ore formation is shown in Figure 6.

Modeling of the Magushan Skarn Deposit
There are three steps in the computation of geoscience workflows, namely geological conceptual, simulation, and mathematical modeling [32][33][34].Here we construct conceptual geological, mathematical, and simulation models for the Magushan Cu-Mo deposit to both advance our knowledge of the processes that control the formation of skarn deposits as well to determine whether numerical simulations can be used to identify prospective areas in regions with thick and barren cover sequences (such as those around the known Magushan deposit) and as a guide for future exploration targeting.The workflow used in this paper to reflect the processes involved in ore formation is shown in Figure 6.

Conceptual Model and Conditions
Previous research on skarn deposits [4] and the energy release that occurs in the subvolcanic environment [35] provide the basis for our conceptual geological model, with a geological model for the formation of the Magushan skarn deposit and the associated generation of lobate skarn boundaries shown in Figure 7.

Conceptual Model and Conditions
Previous research on skarn deposits [4] and the energy release that occurs in the subvolcanic environment [35] provide the basis for our conceptual geological model, with a geological model for the formation of the Magushan skarn deposit and the associated generation of lobate skarn boundaries shown in Figure 7.

Modeling of the Magushan Skarn Deposit
There are three steps in the computation of geoscience workflows, namely geological conceptual, simulation, and mathematical modeling [32][33][34].Here we construct conceptual geological, mathematical, and simulation models for the Magushan Cu-Mo deposit to both advance our knowledge of the processes that control the formation of skarn deposits as well to determine whether numerical simulations can be used to identify prospective areas in regions with thick and barren cover sequences (such as those around the known Magushan deposit) and as a guide for future exploration targeting.The workflow used in this paper to reflect the processes involved in ore formation is shown in Figure 6.

Conceptual Model and Conditions
Previous research on skarn deposits [4] and the energy release that occurs in the subvolcanic environment [35] provide the basis for our conceptual geological model, with a geological model for the formation of the Magushan skarn deposit and the associated generation of lobate skarn boundaries shown in Figure 7.The conceptual model developed during this study for the Magushan deposit is based on previous research on skarn mineralizing environments [4] combined with the geological characteristics of the deposit (Figures 4 and 5b).The entirety of the syncline is not shown in Figure 5b but the data in this figure was combined with the interpreted section shown in Figure 4, with the results shown in Figure 8.No research has been undertaken on the diagenetic processes that have affected these sediments, meaning there is a lack of constraints that would enable the entirely accurate restoration of denuded layers in our model; as an alternative, we have estimated values based on the stratigraphic column shown in Table 1.The emplacement of the Early Cretaceous intrusions in this area indicates that this syncline formed prior to this magmatic event.This in turn indicates that the Cretaceous Zhongfencun and Gecun formations have been denuded in this area and the pre-denudation restoration should include additional layers thicker than 1020 m (the minimum thickness of the Zhongfencun and Gecun formations).The mean thickness of the Cretaceous sediments in this area is 1300 m, indicating that a 1.3 km layer should be added onto the syncline model based on the interpreted geophysical data to account for denudation (Figure 8).Previous research indicated that early skarn-related metamorphism within the Nanling-Xuancheng mining district occurred at temperatures of 270-410 • C [36] and we used a temperature of 700 • C for the intrusion in our model with a constant source at 700 • C, reflecting the later stages of granodiorite magma crystallization and solidification.Other parts of the model have a temperature gradient of 25 • C/km that increases from a surface temperature of 20 • C. Our model also incorporates a pressure gradient of 26.5 MPa/km that reflects the average density of the rocks in this area and a surface pressure that is equal to atmospheric pressure.The model has open boundaries and boundary conditions that are the same as the geothermal and lithostatic gradients outlined above.Our model had a run time of 50,000 years with individual time steps of 25 years and yielded transient results.Our simulation model represents the timing of events after magma upwelling and intrusion but before skarn generation, reflecting the relative timing of mineralization within this area enabling the identification of areas that might contain hitherto unknown areas of mineralization that represent targets for future exploration.

Simulation Modeling
Our simulation model is based on geophysical interpretations and cross-sections for the study area (Figures 4 and 5b).The intrusion associated with the Magushan deposit is located along the contact between Permian and Carboniferous sedimentary units with these contacts between the intrusion and the surrounding sediments also concentrating the formation of skarn mineralization in this area (Figure 8).This model involves fluid flow, heat transfer, chemical reaction, and material migration and the conditions and boundaries conditions are as stated above with the exception of the setting of the reactant concentrations within the hydrothermal fluids involved in this model.The characteristics of the different lithological units used in this model are outlined in Table 2.
Skarn-type ore-forming processes usually follow intrusive events, with this skarn-type metamorphism causing calcite within the limestone to metamorphose to garnet, increasing the density and decreasing the volume of the resulting skarn.These processes increase porosity and permeability, enabling the penetration of ore-forming fluids usually during retrograde skarn formation rather than prograde contact metamorphism [37].Garnet is the main ore-bearing gangue mineral in the Magushan deposit and as such is the main focus of the chemical reaction component of our model.The main chemical reactions during the ore-forming processes that generated the Magushan deposit occurred at temperatures of 270 • C-410 • C [36] and can be simplified as shown in Equation [4]: where X is calcium within calcite in the limestone and Y and SiO 4 are material derived from the magma.We used Y and SiO 4 concentrations of 3 × 10 −3 mol/m 3 and a concentration of 3 × 10 −2 mol/m 3 for X; no fluid inclusions microthermometry or analysis has been undertaken in the area so these represent estimates of the concentrations of these reactants within the ore-forming fluids.In addition, the equilibrium coefficient of this reaction is unknown, meaning that rather than actual concentration values our model provide relative results, meaning that these customized values are suitable for the identification of areas where this reaction (and therefore potentially mineralization) has been concentrated (in relative terms).The computed concentrations and the spatial distribution of garnet generation that resulted from this modeling indicates that the northwest limb of the syncline is most likely a prospective part of the Magushan area and should be the focus of future exploration (Figures 4  and 5).

Mathematical Modeling
The main driving forces that enable the release and subsequent evolution of exsolved magmatic-hydrothermal fluids from a magma body are dynamic porosity and permeability, fluid-flow driven by temperature and pressure differences, and chemical reactions within pore spaces.The entire magmatic-hydrothermal system involves a number of physical and chemical processes that are controlled by fluid-flow, heat transfer, mechanics, chemical reactions, and material migration.These factors affect each other and combine in various ways in the overall process of ore formation (Figure 9; note that our research on the Magushan deposit does not involve either stress fields or deformation).
The main associated processes and their formulas and equations are given in the following sections and all of the key symbols used in these equations are outlined in Table S1.The main associated processes and their formulas and equations are given in the following sections and all of the key symbols used in these equations are outlined in Table S1.

Heat Transfer and the Body Force
Temperature and pressure gradients are calculated using: where  is temperature in °C,  is room temperature (normally 20°C at 1 atmospheric pressure, i.e., 0.1013 MPa),  is depth in m,  is the temperature gradient within the geological units in the model (25°C/km) [38,39],  is pressure in MPa,  is the atmospheric pressure, and  is the pressure gradient within the geological units in the model, reflecting the average density of these units (26.5 MPa/km; Table 3) [40][41][42].The equations that describe the conservation of energy are given below:

Heat Transfer and the Body Force
Temperature and pressure gradients are calculated using: where T is temperature in • C, T 0 is room temperature (normally 20 • C at 1 atmospheric pressure, i.e., 0.1013 MPa), y is depth in m, G T is the temperature gradient within the geological units in the model (25 • C/km) [38,39], P is pressure in MPa, P 0 is the atmospheric pressure, and G P is the pressure gradient within the geological units in the model, reflecting the average density of these units (26.5 MPa/km; Table 3) [40][41][42].
Table 3. Pressure gradients within different layers in the lithosphere (where dp is pressure increment, dh is depth increment, g is the gravity constant, and ρ is the density of the object).Adapted from [42].The equations that describe the conservation of energy are given below:

Formula
where d z is the thickness of the section (set as 1200 m), ρ is fluid density, C p is the heat capacity of the mixture of fluid and chemical reactants (for details see Section 3), T is temperature, t is time, ν is the speed of fluid flow, q is heat dissipation, q 0 is the generalized inward heat flux, Q is reaction heat, Q vd is the heat transferred from the surrounding environment, ρ p is the density of porous rock, θ p is the volume fraction of porous rock, C p,p is the specific heat capacity, k e f f is the effective heat conductivity, C pr is the heat conductivity of porous rock, C m is the heat conductivity of the mixture of fluid and chemical reactants (for details see Section 3), k disp is the coefficient of the initial heat dissipation to the surrounding environment, and ε is porosity.

Fluid-Flow Driven by Pressure Differentiation
Fluid-flow in our model is governed by Darcy's law [43][44][45]: where ν is the speed of fluid flow, k is the permeability, µ is viscosity, ε is porosity, Q m is the fluid source term (fluids are modeled as vectors in the COMSOL software package, so they have a direction and a value but no entity), Q m is used to describe the total of the fluid flow in/through the model system (this term would be a sink term when fluids are flowing out of the model; [30]), t is time, ∇P is the pressure gradient, and gρ l is the body force (gravity and liquid density) where g is the gravitational constant and ρ l is the density of the liquid in question.

Chemical Reactions and Coupling Equations with Heat Transfer
The chemical reaction component of our modeling used a typical garnet formation equation to describe chemical reactions during ore-forming processes at temperatures of 270-410 • C [36] as follows: C m = 0.5 where R Garnet is the reaction rate during the simulated process, w i is the ratio fraction of material i, C p,i is the heat capacity of material i, M i is the molar weight of material i, c i is the concentration of material i; x i is the concentration fraction of material i, and C c,i is the heat conductivity of material i, where i can be X, Y, or SiO 4 .

Transport of Dilute Minerals in Porous Rocks
The transportation of the fluids that precipitate minerals through rocks using Darcy's law is implicitly related to the nature of the porous media, as reflected by the fluid-flow equation used for porous rocks: where D i is the diffusion coefficient of reactant i in fluid, ν is the speed of fluid flow, R i is the reaction rate of reactant i, and N i is the amount of reactant i that takes part in the chemical reaction.

Simulated Ore-Forming Processes
Our modeling predicts the generation of garnet within the contact zone between the intrusion and the surrounding sedimentary units as well as the surrounding regions within both southeast and northwest limbs of the syncline (Figure 10a).This is consistent with the presence of known mineralization within the southeast limb of this fold (i.e., the Magushan deposit) but also suggests that the other limb of this syncline represents a target for mineral exploration.As mentioned above, the equilibrium coefficient of the reaction equation is unknown, meaning that the entire simulation has a low positive equilibrium status and yields garnet concentrations with values of 10 −21 , in practice, and as also mentioned above, this means that the concentration values obtained during this modeling are relative, and the actual values have no practical significance.
where  is the diffusion coefficient of reactant  in fluid, ν is the speed of fluid flow,  is the reaction rate of reactant , and  is the amount of reactant  that takes part in the chemical reaction.

Simulated Ore-Forming Processes
Our modeling predicts the generation of garnet within the contact zone between the intrusion and the surrounding sedimentary units as well as the surrounding regions within both southeast and northwest limbs of the syncline (Figure 10a).This is consistent with the presence of known mineralization within the southeast limb of this fold (i.e. the Magushan deposit) but also suggests that the other limb of this syncline represents a target for mineral exploration.As mentioned above, the equilibrium coefficient of the reaction equation is unknown, meaning that the entire simulation has a low positive equilibrium status and yields garnet concentrations with values of 10 −21 , in practice, and as also mentioned above, this means that the concentration values obtained during this modeling are relative, and the actual values have no practical significance.The mineralized potential of the northwest limb of the syncline was determined by calculating garnet area ratio values (i.e. the proportion of the area of the limb containing predicted garnet) for the data shown in Figure 10b-d, with the results of this calculation shown in Figure 11 and given in Table 4.The northwest limb has a higher garnet area ratio (2.2%) than the southeast limb (1.98%) but at a low concentration (>2 × 10 −21 mol/m 3 ).However, these area ratio values decrease rapidly in both limbs when considering areas of moderate or greater garnet concentration only (>4 × 10 −21 mol/m 3 ; 0.96% in the northwest limb and 0.69% in the southeast limb) before decreasing even further in both limbs (0.24% in the northwest limb and 0.31% in the southeast limb) for areas with only high concentrations of garnet (>6 × 10 −21 mol/m 3 ).It is also important to note that the distribution of The mineralized potential of the northwest limb of the syncline was determined by calculating garnet area ratio values (i.e., the proportion of the area of the limb containing predicted garnet) for the data shown in Figure 10b-d, with the results of this calculation shown in Figure 11 and given in Table 4.The northwest limb has a higher garnet area ratio (2.2%) than the southeast limb (1.98%) but at a low concentration (>2 × 10 −21 mol/m 3 ).However, these area ratio values decrease rapidly in both limbs when considering areas of moderate or greater garnet concentration only (>4 × 10 −21 mol/m 3 ; 0.96% in the northwest limb and 0.69% in the southeast limb) before decreasing even further in both limbs (0.24% in the northwest limb and 0.31% in the southeast limb) for areas with only high concentrations of garnet (>6 × 10 −21 mol/m 3 ).It is also important to note that the distribution of predicted garnet formation matches the locations of known mineralization within the southeast limb despite the fact that this modeling had no inputs (e.g., training points) derived from the known mineralization within the Magushan deposit.This indicates that the northwest limb of the syncline is also likely to host mineralization, potentially over a larger area than that associated with the known Magushan deposit in the southeast limb, despite the fact that the former area contains little known mineralization (primarily as a result of a low level of exploration).We have taken this a step further and have identified a target within the northwest limb of the syncline based on the results of our simulation and the previously constructed geophysics-based cross-sections for this region (Figure 12), providing targets for future exploration in the region around the present Magushan deposit.
predicted garnet formation matches the locations of known mineralization within the southeast limb despite the fact that this modeling had no inputs (e.g.training points) derived from the known mineralization within the Magushan deposit.This indicates that the northwest limb of the syncline is also likely to host mineralization, potentially over a larger area than that associated with the known Magushan deposit in the southeast limb, despite the fact that the former area contains little known mineralization (primarily as a result of a low level of exploration).We have taken this a step further and have identified a target within the northwest limb of the syncline based on the results of our simulation and the previously constructed geophysics-based cross-sections for this region (Figure 12), providing targets for future exploration in the region around the present Magushan deposit.predicted garnet formation matches the locations of known mineralization within the southeast limb despite the fact that this modeling had no inputs (e.g.training points) derived from the known mineralization within the Magushan deposit.This indicates that the northwest limb of the syncline is also likely to host mineralization, potentially over a larger area than that associated with the known Magushan deposit in the southeast limb, despite the fact that the former area contains little known mineralization (primarily as a result of a low level of exploration).We have taken this a step further and have identified a target within the northwest limb of the syncline based on the results of our simulation and the previously constructed geophysics-based cross-sections for this region (Figure 12), providing targets for future exploration in the region around the present Magushan deposit.Variations in garnet area ratios reflect the proportion of these regions that contains garnet at different relative concentrations.Note that the concentration values obtained during this modeling are relative and the actual values have no practical significance beyond highlighting relative changes in predicted garnet abundance.

Sensitivity Testing
The shape of an intrusion can influence the results of numerical simulations like those undertaken during this study.Here, we present the results of a series of sensitivity tests to verify the influence of intrusion shape on the results of our simulations.This involved the construction of two additional models to verify the influence caused by the shape of the apophyses of the intrusion as shown in Figure 13a,b.Our original simulation model was based on interpreted geophysical data and the resulting geological cross-sections.However, the size of the northwestern apophyses is somewhat uncertain and was originally modeled to have the same dimensions as the southeastern apophyses (Figure 8; the interpreted section is too larger to show the relatively small differences in the sizes of the two limbs).We tested the influence of this uncertainty and any possible variation in the size of this part of the intrusion by constructing two additional models that contained northwestern intrusion apophyses that varied in size, with Figure 13a showing the model with the wider northwestern apophyses, and Figure 13b showing the model with the narrower apophyses.Variations in garnet area ratios reflect the proportion of these regions that contains garnet at different relative concentrations.Note that the concentration values obtained during this modeling are relative and the actual values have no practical significance beyond highlighting relative changes in predicted garnet abundance.

Sensitivity Testing
The shape of an intrusion can influence the results of numerical simulations like those undertaken during this study.Here, we present the results of a series of sensitivity tests to verify the influence of intrusion shape on the results of our simulations.This involved the construction of two additional models to verify the influence caused by the shape of the apophyses of the intrusion as shown in Figure 13a-b.Our original simulation model was based on interpreted geophysical data and the resulting geological cross-sections.However, the size of the northwestern apophyses is somewhat uncertain and was originally modeled to have the same dimensions as the southeastern apophyses (Figure 8; the interpreted section is too larger to show the relatively small differences in the sizes of the two limbs).We tested the influence of this uncertainty and any possible variation in the size of this part of the intrusion by constructing two additional models that contained northwestern intrusion apophyses that varied in size, with Figure 13a showing the model with the wider northwestern apophyses, and Figure 13b showing the model with the narrower apophyses.The results of this sensitivity testing are shown in Figure 14, which shows the results of modeling with a wider northwestern apophysis (Figures 13a, 14a1-a4), the original modeling (Figures 10, 14b1-b5), and with a narrower northwestern apophysis (Figures 13b, 14c1-c4).The northwestern apophyses contain garnet independent of the size of this part of the intrusion, emphasizing that even with this uncertainty this region should be considered highly prospective for future exploration.This exploration potential was again quantified using garnet-bearing area ratios (Figure 15; Table 5).The sensitivity testing models essentially replicate the results of the original modeling, where more of the northwest limb contains garnet at low concentrations (>2 × 10 −21 mol/m 3 ), but the proportion of this area containing garnet drops rapidly to below the proportion of the southeast limb when considering areas with higher concentrations of garnet (>6 × 10 −21 mol/m 3 ), with these results consistent over all three models.Our sensitivity modeling therefore indicates that the broad conclusions from our original model are independent of variations in size of the northwestern apophysis, namely that a significant region of the northwest limb should The results of this sensitivity testing are shown in Figure 14, which shows the results of modeling with a wider northwestern apophysis (Figures 13a and 14a1-a4), the original modeling (Figure 10 and Figure 14b1-b5), and with a narrower northwestern apophysis (Figure 13b and Figure 14c1-c4).The northwestern apophyses contain garnet independent of the size of this part of the intrusion, emphasizing that even with this uncertainty this region should be considered highly prospective for future exploration.This exploration potential was again quantified using garnet-bearing area ratios (Figure 15; Table 5).The sensitivity testing models essentially replicate the results of the original modeling, where more of the northwest limb contains garnet at low concentrations (>2 × 10 −21 mol/m 3 ), but the proportion of this area containing garnet drops rapidly to below the proportion of the southeast limb when considering areas with higher concentrations of garnet (>6 × 10 −21 mol/m 3 ), with these results consistent over all three models.Our sensitivity modeling therefore indicates that the broad conclusions from our original model are independent of variations in size of the northwestern apophysis, namely that a significant region of the northwest limb should contain mineralization (over a slightly larger area than the mineralization within the southeast limb), but this mineralization may be less well developed (lower grade?) than that within the southeast limb.
contain mineralization (over a slightly larger area than the mineralization within the southeast limb), but this mineralization may be less well developed (lower grade?) than that within the southeast limb.contain mineralization (over a slightly larger area than the mineralization within the southeast limb), but this mineralization may be less well developed (lower grade?) than that within the southeast limb.Variations in garnet area ratios reflect the proportion of these regions that contains garnet at different relative concentrations.Note that the concentration values obtained during this modeling are relative and the actual values have no practical significance beyond highlighting relative changes in predicted garnet abundance.

Discussion
Mineral exploration in regions containing thick and unprospective cover sequences requiring the targeting of deep-seated 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 [46][47][48][49][50]. Recent advances in three dimensional prospectivity modeling for mineral exploration has highlighted the potential use of prospectivity modeling [12][13][14][15][16] and numerical modeling [21] in this type of exploration targeting.However, the simulation undertaken during this study provides another method of predicting the location of areas that are prospective for deep exploration, supplementing more traditional exploration approaches, and removing risks related to the conditional dependence associated with some other types of prospectivity modeling.
Our modeling has correctly predicted the location of areas of known mineralization despite being based on modeling rather than any training of our dataset using areas of previously known mineralization.This indicates that this simulation approach can provide an independent identification of highly prospective targets for future exploration as well as having potential use in exploration for deep-seated mineralization elsewhere.This is consistent with the results of [21], who used a numerical simulation within a three-dimensional prospectivity analysis of the skarn-dominated Yueshan orefield in China.This research indicated that the incorporation of numerical simulation enhanced the results of prospectivity modeling, indicating the usefulness of numerical simulation in exploration at depth.
This study used numerical simulation to identify areas for future exploration in areas surrounding regions of known skarn mineralization, but the effectiveness of this approach on larger scales (e.g., entire orefields) requires further assessment.The Magushan orefield is located within the southern part of the newly discovered Nanling-Xuancheng mining district, an area with high exploration potential after the discovery and exploration of the Chating porphyry [31] and Magushan skarn [22] deposits.However, the majority of the area is covered by Quaternary and thick Cretaceous sediments (Figure 3), meaning this region contains little exposure of older rocks that may host or be associated with mineralization.This presents a significant challenge to traditional mineral exploration approaches, meaning that the construction of a simplified orefield-scale numerical model and the associated numerical simulation of mineralized processes based on the Magushan deposit and existing geological and geophysical data for this area could potentially increase exploration success in this region.This approach is far more economical than the significant amount of drilling that is needed to test the entire area, which covers some 23 × 21 km.This situation 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 more important over time, and this type of numerical simulation should be more widely applied during future mineral exploration.
However, there are still shortcomings in our research.The numerical model is in a 2D environment, and the influence from the surrounding environment (both sides of the section) is ignored.In addition, the 3D shape of the contact zone between geological units (i.e., sedimentary units and the intrusion) may influence the distribution of the concentration of garnet and associated mineralization in this region.A more precisely interpreted and better constrained geophysical section would also provide more insight into the geology and mineralization in this area, enabling the construction of an improved and more accurate model and related multi-field coupling simulations that in turn more accurately predicts the location of highly prospective areas to be targeted during future mineral exploration.

Conclusions
This study numerically modeled the processes that formed the skarn-type Magushan copper-molybdenum deposit and producing the following key findings: (1) The simulated distribution of garnet matches areas of known mineralization identified during drilling and exploitation of the deposit to date, indicating that our modeling could potentially be used to predict areas of hitherto unknown mineralization as well as verifying the potential use of a coupled multi-field approach in economic geology research and mineral exploration.
(2) Our modeling also identified a prospective exploration target to the northwest of the present Magushan deposit that appears to be highly prospective when combining our modeling with the results of previous geological and geophysical research.This indicates that numerical simulation represents an effective approach to prospectivity modeling and exploration targeting in areas with thick and barren sedimentary (or other) cover sequences where traditional geological and geochemical techniques are less effective.

Figure 2 .
Figure 2. Map showing the geology and mineral deposits of the Nanling-Xuancheng mining camp.Modified from [31].

Figure 3 .
Figure 3. Map showing the geology and mineral deposits of the Magushan orefield.Modified from [31].

Figure 2 .
Figure 2. Map showing the geology and mineral deposits of the Nanling-Xuancheng mining camp.Modified from [31].

Figure 2 .
Figure 2. Map showing the geology and mineral deposits of the Nanling-Xuancheng mining camp.Modified from [31].

Figure 3 .
Figure 3. Map showing the geology and mineral deposits of the Magushan orefield.Modified from [31].

Figure 3 .
Figure 3. Map showing the geology and mineral deposits of the Magushan orefield.Modified from [31].

Figure 4 .
Figure 4. Interpreted geological section based on a constrained joint inversion of gravity and magnetic potential field data (Section A-A' in Figure 3).

Figure 4 .
Figure 4. Interpreted geological section based on a constrained joint inversion of gravity and magnetic potential field data (Section A-A' in Figure 3).

Figure 5 .
Figure 5. Map showing (a) the geology of the area around the skarn-type Magushan coppermolybdenum deposit (omitting Quaternary sediments for clarity) and (b) a cross-section along B-B'showing the geology of the main orebody within the deposit.Modified from[22].

Figure 5 .
Figure 5. Map showing (a) the geology of the area around the skarn-type Magushan copper-molybdenum deposit (omitting Quaternary sediments for clarity) and (b) a cross-section along B-B' showing the geology of the main orebody within the deposit.Modified from [22].

Figure 6 .
Figure 6.Workflow used in the numerical modeling of the Magushan ore deposit.

Figure 7 .
Figure 7. Simplified schematic illustration the formation of lobate skarn boundaries such as those within the Magushan deposit and corresponding fluid-flow patterns.Modified from [4].

Figure 6 .
Figure 6.Workflow used in the numerical modeling of the Magushan ore deposit.

Figure 6 .
Figure 6.Workflow used in the numerical modeling of the Magushan ore deposit.

Figure 7 .
Figure 7. Simplified schematic illustration the formation of lobate skarn boundaries such as those within the Magushan deposit and corresponding fluid-flow patterns.Modified from [4].

Figure 7 .
Figure 7. Simplified schematic illustration the formation of lobate skarn boundaries such as those within the Magushan deposit and corresponding fluid-flow patterns.Modified from [4].

Figure 8 .
Figure 8. Simulation model used during modeling constructed according to the conceptual model outlined in the text.Model state is after magma solidifies.

Figure 8 .
Figure 8. Simulation model used during modeling and constructed according to the conceptual model outlined in the text.Model state is after the solidification of the magma in this region.

Figure 9 .
Figure 9.The coupled relationships between fluid flow, mechanics, heat transfer, material migration, and chemical reactions during mineralization and mineral deposit formation, with colors indicating the different physical fields involved in these processes.

Figure 9 .
Figure 9.The coupled relationships between fluid flow, mechanics, heat transfer, material migration, and chemical reactions during mineralization and mineral deposit formation, with colors indicating the different physical fields involved in these processes.

Figure 11 .
Figure 11.Graph showing variations in garnet area ratios (i.e. the proportion of a region that contains garnet at different relative concentrations) within the northwest and southeast limbs of the anticline with changes in minimum garnet concentrations.Note that the concentration values obtained during this modeling are relative and the actual values have no practical significance beyond highlighting relative changes in predicted garnet abundance.

Figure 12 .
Figure 12.Map showing the location of the delineated target within the northwest limb of the syncline in relation to the rest of the Magushan orefield.

Figure 11 .
Figure 11.Graph showing variations in garnet area ratios (i.e., the proportion of a region that contains garnet at different relative concentrations) within the northwest and southeast limbs of the anticline with changes in minimum garnet concentrations.Note that the concentration values obtained during this modeling are relative and the actual values have no practical significance beyond highlighting relative changes in predicted garnet abundance.

Figure 11 .
Figure 11.Graph showing variations in garnet area ratios (i.e. the proportion of a region that contains garnet at different relative concentrations) within the northwest and southeast limbs of the anticline with changes in minimum garnet concentrations.Note that the concentration values obtained during this modeling are relative and the actual values have no practical significance beyond highlighting relative changes in predicted garnet abundance.

Figure 12 .
Figure 12.Map showing the location of the delineated target within the northwest limb of the syncline in relation to the rest of the Magushan orefield.

Figure 12 .
Figure 12.Map showing the location of the delineated target within the northwest limb of the syncline in relation to the rest of the Magushan orefield.

Figure 13 .
Figure 13.Models used in sensitivity testing to determine the effect of variations in intrusion size: model with (a) a wider northwestern apophysis limb and (b) a narrower northwestern apophysis.

Figure 13 .
Figure 13.Models used in sensitivity testing to determine the effect of variations in intrusion size: model with (a) a wider northwestern apophysis limb and (b) a narrower northwestern apophysis.

Figure 14 .
Figure 14.The results of sensitivity testing undertaken during this study, comparing the distribution of generated garnet in (a1-a4) the model with a wider northwestern apophysis, (b1-b4) the original model, and (c1-c4) a model with a narrower northwestern apophysis.

Figure 15 .
Figure 15.Graphs showing variations in garnet-bearing area ratio values within the northwest and southwest limbs of the syncline in the study area relative to variations in garnet abundances as a result of modeling with (a) a wider northwestern apophysis, (b) the original modeling setup, and (c) with a narrower northwestern apophysis.Note that the concentration values obtained during this modeling are relative and the actual values have no practical significance beyond highlighting relative changes in predicted garnet abundance

Figure 14 .
Figure 14.The results of sensitivity testing undertaken during this study, comparing the distribution of generated garnet in (a1-a4) the model with a wider northwestern apophysis, (b1-b4) the original model, and (c1-c4) a model with a narrower northwestern apophysis.

Figure 14 .
Figure 14.The results of sensitivity testing undertaken during this study, comparing the distribution of generated garnet in (a1-a4) the model with a wider northwestern apophysis, (b1-b4) the original model, and (c1-c4) a model with a narrower northwestern apophysis.

Figure 15 .
Figure 15.Graphs showing variations in garnet-bearing area ratio values within the northwest and southwest limbs of the syncline in the study area relative to variations in garnet abundances as a result of modeling with (a) a wider northwestern apophysis, (b) the original modeling setup, and (c) with a narrower northwestern apophysis.Note that the concentration values obtained during this modeling are relative and the actual values have no practical significance beyond highlighting relative changes in predicted garnet abundance

Figure 15 .
Figure 15.Graphs showing variations in garnet-bearing area ratio values within the northwest and southwest limbs of the syncline in the study area relative to variations in garnet abundances as a result of modeling with (a) a wider northwestern apophysis, (b) the original modeling setup, and (c) with a narrower northwestern apophysis.Note that the concentration values obtained during this modeling are relative and the actual values have no practical significance beyond highlighting relative changes in predicted garnet abundance.

Table 1 .
Stratigraphy of the study area; adapted from [23].

Table 1 .
Stratigraphy of the study area; adapted from [23].

Table 2 .
Rock material parameters used in simulation modeling; adapted from [23].

Table 3 .
[42]sure gradients within different layers in the lithosphere (where  is pressure increment, ℎ is depth increment,  is the gravity constant, and  is the density of the object).Adapted from[42].

Table 4 .
Garnet area ratio values for the northwest and southeast limb of the syncline.

Table 4 .
Garnet area ratio values for the northwest and southeast limb of the syncline.

Table 5 .
Garnet area ratio values for the northwest and southeast limb of the syncline.

Table 5 .
Garnet area ratio values for the northwest and southeast limb of the syncline.

Table 5 .
Garnet area ratio values for the northwest and southeast limb of the syncline.