3D Mineral Prospectivity Modeling for the Low-Sulfidation Epithermal Gold Deposit: A Case Study of the Axi Gold Deposit, Western Tianshan, NW China

The Axi low-sulfidation (LS) epithermal deposit in northwestern China is the result of geological controls on hydrothermal fluid flow through strike-slip faults. Such controls occur commonly in LS epithermal deposits worldwide, but unfortunately, these have not been quantitatively analyzed to determine their spatial relationships with gold distribution and further guide mineral prospecting. In this study, we conduct a 3D mineral prospectivity modeling approach for the Axi deposit involving 3D geological modeling, 3D spatial analysis, and prospectivity modeling. The spatial analysis of geometric features revealed the gold mineralization trends in convex segments (0–20 m) with a specific distance from fault 2, the lower interface of late volcanic phase, and the upper interface of phyllic alteration with steep slopes (>65◦), implying that gold deposition was significantly controlled by the morphological characteristics and distance fields of geologic features. The present alteration–mineralization zone at Axi has a larger width in bending sites (sections No. 35–15 and No. 40–56) than elsewhere, indicating the location of two fluid conduits extending to depth. The prediction-area plots and receiver operating characteristic curves demonstrated that (genetic algorithm optimized support vector regression (GA-SVR)) outperformed multiple nonlinear regression and fuzzy weights-of-evidence, which was proposed as a robust method to solve complicated nonlinear and high-dimensional issues in prospectivity modeling. Our study manifests spatial controls of structure, host rock, and alteration on LS epithermal gold deposition, and highlights the capability of GA-SVR for identifying deposit-scale potential epithermal gold mineralization.


Introduction
Low-sulfidation (LS) epithermal deposits with major Ag and Au production are formed in shallow parts (<1 km) of magmatic-hydrothermal systems at moderate-low temperature conditions of 100 • C to 300 • C [1][2][3][4][5]. The ore-forming fluids of LS epithermal deposits have been fully investigated and identified as dominantly meteoric water with minor but crucial additions of magmatic water, and their mixing, boiling, and interaction with country rocks are proposed as principal mineralization Minerals 2020, 10 mechanisms [5][6][7][8]. In spite of the consensus about metallogenic aspects, the spatial pattern of LS epithermal mineralization at the deposit scale remains poorly understood. One fundamental understanding of the localization of LS epithermal mineralization is the control on metal transport and deposition during epizonal fluid flow [1,[9][10][11][12][13]. It, to some extent, requires a special tool to quantitatively analyze the spatial association between the metal distribution and geologic features plausibly related to mineralization [14][15][16]. However, research regarding geological controls on LS epithermal mineralization at the deposit-scale is rarely conducted, hindering the further understanding of the LS epithermal ore-forming processes and mineral exploration.
3D spatial analysis has been demonstrated as an effective approach for identifying the spatial features of geological objects in 3D space and extracting diverse potential ore-controlling features for 3D prospectivity modeling [17][18][19][20][21][22]. However, the spatial association between ore and ore-controlling features would be significantly more complicated with the expansion of geological features in space [23][24][25]. In terms of the non-linearity and high-dimensionality of spatial data, traditional 3D prospectivity methods remain inadequate, thereby resulting in larger error and bias in prospectivity modeling [19,20]. Recent studies employ machine learning algorithms in prospectivity modeling, for they have good spatial dataset processing performance unrelated to statistical distribution and automatically generate high-level representations from complex data [19,[23][24][25][26][27][28]. Genetic algorithm optimized support vector regression (GA-SVR), a supervised learning regression method extended from support vector machine (SVM) with parameters optimized by a genetic algorithm, is a machine learning algorithms with better stability and generalization abilities in trouble-shooting complicated nonlinear and high dimensional issue and may be a promising tool in 3D prospectivity modeling [29][30][31][32].
The Axi gold deposit, as the largest LS epithermal deposit in the West Tianshan, has been the subject of extensive geological studies and drill exploration, that provide important insights into ore genesis and abundant spatial information of geological bodies. In this study, the Axi deposit was chosen as a representative research case of LS epithermal deposits to implement improved 3D prospectivity modeling, including 3D geological modeling, multiple spatial analysis of potential ore-controlling geological features, and GA-SVR prospectivity modeling. The work presented herein focuses on the spatial pattern of LS epithermal mineralization and its controls from geological features, as well as reducing bias and uncertainty in prospectivity modeling by a machine learning method. Its implementation would help us (1) to better understand the gold spatial deposition during fluid flow, (2) to guide deep mineral exploration at Axi, (3) to verify the availability of GA-SVR in prospectivity modeling, and finally, (4) to provide the reference for mineral prospecting in other LS epithermal gold deposits.

Regional and Deposit Geology
The Tulasu Basin in Western Tianshan, NW China, hosting several Late Paleozoic porphyry, skarn, and epithermal deposits, is sandwiched between the NWW-striking North Yili Basin Fault and South Keguqinshan Mountain Fault ( Figure 1) [33][34][35]. The basement of the Tulasu basin is mainly subdivided into the lower part with mesoproterozoic and neoproterozoic sedimentary clastic rocks and the upper part with ordovician shallow marine clastic rocks with interlayered carbonates [36,37]. It is unconformably covered by paleozoic strata of the lower carboniferous Dahalajunshan formation and overlying the lower carboniferous aqialehe formation. The Dahalajunshan formation, generated in a continental arc setting related to the subduction of the North Tianshan Ocean [33,38,39], consists of tuff, basaltic andesite, andesite, dacite, and rhyolite, with minor volcanoclastic breccia. The aqialehe formation is mainly comprised of conglomerate, sandstone, and limestone [38]. The granitic intrusions (e.g., diorite, granodiorite, and granitic porphyry) exposed in the Tulasu basin commonly intrude into the Dahalajunshan andesite, with the early carboniferous age of ca. 357-348 Ma (Figure 1a) [40]. slip extension, and post-ore compression and overprinting by secondary reverse faults (Figure 1b) [41]. The well-developed alteration-mineralization zoning pattern in the Axi epithermal deposit, from the core outward, mainly consists of silicification, quartz-carbonate, pyrite-sericite-quartz, and propylitic alteration ( Figure 2) [33][34][35][36]. In general, pyrite-sericite-quartz and silicification close to gold mineralization occur in the hanging wall of fault F2, and the pervasive propylitic alteration affects all lithologies in the Axi district.  The Axi gold deposit, located in the central part of the Tulasu Basin, is hosted in andesitic and dacitic volcanic rocks and breccia of the uppermost member of the Dahalajunshan formation ( Figure 1b) [34,35,37]. The strike-slip fault (F 2 ), which is N-trending and east-dipping, controls nearly all economic orebodies in the Axi gold deposit. The radial fault series developed in the western segment rarely exhibit alteration and mineralization, rather, these likely reflect a structural overprinting on the late evolution of a diatreme at Axi. The fault F 2 displays variable strikes from south to north due to multiple phases of deformation associated with pre-ore strike-slip extension, syn-ore normal-strike slip extension, and post-ore compression and overprinting by secondary reverse faults (Figure 1b) [41]. The well-developed alteration-mineralization zoning pattern in the Axi epithermal deposit, from the core outward, mainly consists of silicification, quartz-carbonate, pyrite-sericite-quartz, and propylitic alteration ( Figure 2) [33][34][35][36]. In general, pyrite-sericite-quartz and silicification close to gold mineralization occur in the hanging wall of fault F 2 , and the pervasive propylitic alteration affects all lithologies in the Axi district. Gold orebodies in the controlled exploration engineering area extend approximately 1300 m from north to south (Figure 1b), with dipping angles of 45° to 80°, an average thickness of 20 m and the maximum depth of about 425 m. Nine orebodies that generally parallel fault F2 in the hanging wall have been identified. These gold orebodies exhibit a similar spatial variation to fault F2, both in vertical and horizontal directions ( Figure 2). Based on the ore textures, compositions, and crosscutting relationships, two subtypes of primary gold ore were identified: disseminated ore (20% of documented gold resources) and quartz-vein ore (80% of documented gold resources) [37]. The disseminated ore is hosted in rocks with strong phyllic alteration that typically has a gray to yellowish-brown color and a relatively lower gold grade compared to vein ore (2 g/t vs. 5.6 g/t). Quartz-vein ore with a width of 0.1 m to 20 m consists mostly of gray, smoky gray quartz-sulfide veins, stockworks, and breccias. The major orebodies commonly contain phyllic altered fragments that were replaced and cemented by quartz-sulfide veins. The absolute ages via auriferous pyrite Re-Os dating for the two types of ore are ca. 356-348 Ma (disseminated) and ca. 332 Ma (veins) [37]. Gold-bearing minerals in disseminated ore are described as coarse-grained euhedral pyrite, whereas the gold in quartz-vein ore is typically related to electrum, native gold, zoned pyrite, arsenopyrite, and marcasite [33][34][35][36]. Gangue minerals in both types of ore are mainly quartz, chalcedony, illite, chlorite, adularia, calcite, and barite.

Metallogenic Model
The formation of the Axi epithermal gold deposit involved two major stages; the disseminated ore was mainly formed in the early ore-stage whereas the vein-type ore was generated in the main ore-stage [35,37]. In the early ore-stage, ore-forming fluids were predominantly magmatic water and had high-temperature and a high concentration of Cu, Ni, Co, Mn, HS − , Cl − [15]. Relatively hot hydrothermal fluids ascended from depth along faults or cracks to shallower levels. These structures, as fluid conduits, were beneficial for gold deposition and hydrothermal alteration when fluid infiltrated and diffused along fracture zones and interplayed with country rocks [35]. During the later ore-stage, ore-forming fluids were progressively diluted and cooled by significant inputs of meteoric water [33,35,36]. Gold was transported as Au-bearing chloride-and sulfur-complexes in ore-forming  Figure 2). Based on the ore textures, compositions, and crosscutting relationships, two subtypes of primary gold ore were identified: disseminated ore (20% of documented gold resources) and quartz-vein ore (80% of documented gold resources) [37]. The disseminated ore is hosted in rocks with strong phyllic alteration that typically has a gray to yellowish-brown color and a relatively lower gold grade compared to vein ore (2 g/t vs. 5.6 g/t). Quartz-vein ore with a width of 0.1 m to 20 m consists mostly of gray, smoky gray quartz-sulfide veins, stockworks, and breccias. The major orebodies commonly contain phyllic altered fragments that were replaced and cemented by quartz-sulfide veins. The absolute ages via auriferous pyrite Re-Os dating for the two types of ore are ca. 356-348 Ma (disseminated) and ca. 332 Ma (veins) [37]. Gold-bearing minerals in disseminated ore are described as coarse-grained euhedral pyrite, whereas the gold in quartz-vein ore is typically related to electrum, native gold, zoned pyrite, arsenopyrite, and marcasite [33][34][35][36]. Gangue minerals in both types of ore are mainly quartz, chalcedony, illite, chlorite, adularia, calcite, and barite.

Metallogenic Model
The formation of the Axi epithermal gold deposit involved two major stages; the disseminated ore was mainly formed in the early ore-stage whereas the vein-type ore was generated in the main ore-stage [35,37]. In the early ore-stage, ore-forming fluids were predominantly magmatic water and had high-temperature and a high concentration of Cu, Ni, Co, Mn, HS − , Cl − [15]. Relatively hot hydrothermal fluids ascended from depth along faults or cracks to shallower levels. These structures, as fluid conduits, were beneficial for gold deposition and hydrothermal alteration when fluid infiltrated and diffused along fracture zones and interplayed with country rocks [35]. During the later ore-stage, ore-forming fluids were progressively diluted and cooled by significant inputs of meteoric water [33,35,36]. Gold was transported as Au-bearing chloride-and sulfur-complexes in ore-forming fluids and the mechanism of gold deposition was interpreted as fluid boiling caused by hydraulic fracturing and/or seismic activity in fluid conduits [33,35,36].
The following empirical exploration markers are recognized in the Axi deposit: (1) Gold occurrence is structurally controlled by the development of fault F 2 , generally occurring in the hanging wall of F 2 . The spatial distribution of gold mineralization is closely associated with the geometry of F 2 . The dilational sites along F 2 with changing dip angle always hold high-grade orebodies.
(2) The late volcanic phase rock (LVP) is the major host rock for gold orebodies. Ore-forming fluids preferentially interacted with LVP that had undergone brecciation and structural deformation. The gold deposition is possibly related to the location and shape of LVP.
(3) Orebodies are commonly within regions that have undergone phyllic alteration and silicification, and high-economic lodes are located in the center of alteration zonation, i.e., silicification. Typically, the scale of the gold mineralization zone is positively correlated to the thickness of alteration zone.

Datasets
In this study, data employed to conduct 3D prospectivity modeling of the Axi gold deposit covered several surface geological maps, 29 cross-sections with 40 m intervals, 211 drill holes with a maximum depth of 653.7 m, as well as the gold assay data of 13,628 samples. These datasets were provided by Western Region Gold Co., Ltd. (Urumqi, Xinjiang Uygur Autonomous Region, China). Databases may be requested by contacting the primary author.

3D Geological Modeling
It is notable that 3D geological modeling was the foundation for spatial analysis and prospecting models, and its quality is closely related to the complexity of geological features, geological interpretation, and exploration data density [20,42]. Two categories of 3D geological modeling methods exist, explicit and implicit. Explicit modeling is more applicable to 3D reconstruction of deposits with abundant geological data with low-uncertainty and can be made more accurate by skilled geological interpretation [20,43,44]. Conversely, implicit modeling is suitable for 3D geological modeling at the district-scale, where geological data is sparse and the geology is poorly understood [43,44].
Since an extensive geological database with comprehensive geological data and many scholarly geological interpretations exist for the Axi deposit, we preferred using a data-and knowledge-driven explicit modeling method. We constructed a series of 3D geological models with the SKUA-GOCAD TM software program (version 2018, Emerson Paradigm Holding LLC, MO, USA), consisting of orebodies, fault F 2 , late volcanic phase (LVP), and phyllic alteration zone ( Figure 3). Minerals 2020, 10, x FOR PEER REVIEW 6 of 21

Spatial Analysis of Geological Features
Many highly regarded studies have determined that the spatial distribution and morphology of geological features control orebody occurrence in epithermal deposits [3,43,45]. Thus, the spatial analysis of geological features is crucial for providing effective mineral exploration criteria [17,18,[46][47][48]. In this study, multi-spatial analysis methods were integrated to quantify their relative spatial control on gold mineralization and further generate a series of predictive maps by spatial distance analysis, slope analysis, and morphological undulation analysis.
The location and enrichment of gold orebodies in the Axi deposit are spatially associated with geological features of F2, LVP, and phyllic alteration zone. The spatial distance analysis method was performed to quantitatively represent the influence of ore-controlling factors [19,20]. We defined distance value (d) as the minimum Euclidean distance between each voxel in 3D prospecting space and the ore-controlling factors and used a sign on the Euclidean distance to depict orientation relationships. For a given voxel, d < 0 indicates the voxel is located in the hanging wall of the geological body, inversely d > 0 illustrated it situated in the footwall.
Gold deposition from ore fluid linked to the destabilization of Au-bearing complexes was affected by geometric features [1,3,47,48]. Slope, as a significant indicator for quantitatively representing geometry, was computed. In this study, 3D geological models of F2, LVP and the phyllic alteration zone were expressed as TIN models with enormous triangular patches, and the plane equation for each patch can be formulated as z = ax + by + c. The slope was defined as the angle between the triangular patch and horizontal surface, calculated as g = arccos 1 + +1 . Changes in the slope of geological surfaces, such as the changes in the slope of a fault plane, commonly determine the formation and location of orebodies [18][19][20]. We conducted a morphological undulation analysis method to compute undulation based on trend analysis [18,20,48].

Spatial Analysis of Geological Features
Many highly regarded studies have determined that the spatial distribution and morphology of geological features control orebody occurrence in epithermal deposits [3,43,45]. Thus, the spatial analysis of geological features is crucial for providing effective mineral exploration criteria [17,18,[46][47][48]. In this study, multi-spatial analysis methods were integrated to quantify their relative spatial control on gold mineralization and further generate a series of predictive maps by spatial distance analysis, slope analysis, and morphological undulation analysis.
The location and enrichment of gold orebodies in the Axi deposit are spatially associated with geological features of F 2 , LVP, and phyllic alteration zone. The spatial distance analysis method was performed to quantitatively represent the influence of ore-controlling factors [19,20]. We defined distance value (d) as the minimum Euclidean distance between each voxel in 3D prospecting space and the ore-controlling factors and used a sign on the Euclidean distance to depict orientation relationships. For a given voxel, d < 0 indicates the voxel is located in the hanging wall of the geological body, inversely d > 0 illustrated it situated in the footwall.
Gold deposition from ore fluid linked to the destabilization of Au-bearing complexes was affected by geometric features [1,3,47,48]. Slope, as a significant indicator for quantitatively representing geometry, was computed. In this study, 3D geological models of F 2 , LVP and the phyllic alteration zone were expressed as TIN models with enormous triangular patches, and the plane equation for each patch can be formulated as z = ax + by + c. The slope was defined as the angle between the triangular patch and horizontal surface, calculated as g = arccos . Minerals 2020, 10, 233 Changes in the slope of geological surfaces, such as the changes in the slope of a fault plane, commonly determine the formation and location of orebodies [18][19][20]. We conducted a morphological undulation analysis method to compute undulation based on trend analysis [18,20,48]. It was subdivided into two steps: (1) extracting trend surface T(x, y, z) from discrete smooth interpolation (DSI) to original geological interface Z(x, y, z); (2) for each point in Z(x, y, z), calculating residual value R(x, y, z), defined as the minimum Euclidean distance from T(x, y, z) to Z(x, y, z) in elevation. R(x, y, z) > 0 indicated a uplift shape, oppositely, R(x, y, z) < 0 illustrated a depression shape.

Prospectivity Modeling and Assessing
Prospectivity modeling entails the analysis and synthesis of various predictive maps derived from the spatial analysis of geological models. Diverse predictive maps can be treated as n-Dimension feature vectors with a nonlinear relationship, hence prospectivity modeling actually is a solution for the complex nonlinear and high dimensional problem [24,49]. This problem is always treated as a classification process that divides the study area into favorable and unfavorable parts [50][51][52], but essentially, it is a multivariate nonlinear regression problem to quantify the spatial association between predictive maps and gold occurrence [24,53]. Support vector regression (SVR), as a supervised learning regression algorithm branching from SVM, has been widely applied for solving the nonlinear regressing estimation and time series forecasting issues [29,31,54,55] and has also been used in geological engineering and the mineral industry [56,57]. Additionally, it can also be a preferable solution for performing complicated nonlinear prospectivity modeling. In this study, a hybrid GA-SVR (genetic algorithm optimized support vector regression) was adopted for resource evaluation and prediction analysis. Simultaneously, both multiple nonlinear regression (MNR) and the fuzzy weights-of-evidence (fuzzy WofE) method were chosen for comparative assessment through receiver operating characteristic (ROC) and prediction-area (PA) analysis [58,59].

Support Vector Regression Model
SVM, proposed on the basis of statistical learning theory, is a supervised learning algorithm for pattern classification [50,51,60]. It implements the global optimum based on the structural risk minimization principle [60]. Extended from SVM for trouble-shooting of nonlinear regression issues, SVR attempts to fit an optimal approximating hyperplane to the training set, through mapping nonlinear input datasets into a multi-dimensional or hyper-dimensional feature space [56,60,61]. Thus, the optimal hyperplane quantitatively describes the relationship between features and target values in the training dataset, to predict the future target values with input features. This method has been shown to be superior in minimizing the expected error and reducing the problem of over-fitting [56,60,62]. ε-SVR, as a popular SVR model, locates the hyperplane with an epsilon-insensitive (ε-insensitive) loss function [60,61]. The ε-SVR function can be obtained as formula (1), where f(x) denotes regression function, x is the training sample, a i andâ i is Lagrange multiplier, k(x, x i ) is the kernel function, x i is the i-th eigenvector of x (the vector of the i-th predictive map), and b is the offset. The kernel function k(x, x i ) is the key for mapping nonlinear input data into a high dimensional space. Diverse kernel functions, such as the polynomial kernel function, the Sigmoid kernel function and the Gaussian radial basis function (RBF), have been used in previous research. Among them, RBF is the most conventionally devoted to ε-SVR modeling, for its high efficiency and lower deviation in sample training [63]. It is defined as follows: where σ stands for the width of the RBF, controlling the solution complexity.
In the aforementioned ε-SVR model with RBF, three parameters need to be optimized: the penalty parameter (c), the insensitive loss function parameter (ε), and the width parameter of the RBF (σ), respectively representing the trade-off between training deviation and model complexity, the sparseness of solution, and the complexity of solution [61,62].
The traditional parameter optimization methods (e.g., grid search, gradient descend, and cross-validation) are insufficient in search efficiency and accuracy [31,54]. Genetic algorithm (GA), an adaptive heuristic algorithm mimicking biological evolution, is one type of global optimization algorithm that has proved to be superior in parameter selection, because of its simple structure, good robustness and intelligence searching [30][31][32]64]. GA treats parameters to be optimized as chromosomes, and determines the optimal model via a computational iterative process based on the survival of the fittest in a population, realizing intelligent exploitation of a random search in a predefined searching space [30,32]. With a good performance in finding the global solution for the optimization problem, GA has been gradually applied to synchronously optimize parameters of Machine Learning, and then hybrid algorithms (e.g., GA-SVR and genetic regulatory networks) have been proposed for predication [64][65][66][67]. Particularly, plenty of studies demonstrate that GA-SVR improves the generalization ability and forecasting accuracy, which has been extensively employed in various research fields [64][65][66].
The hybrid GA-SVR (genetic algorithm optimized support vector regression) that was adopted for prospectivity modeling involved several processes ( Figure 4): (1) normalizing all predictive maps proposed in Section 3.3 into values with zero mean, to eliminate the differences of units and magnitudes between these selected parameters; (2) dividing known data into two datasets, namely a training dataset and a test dataset. The former is applied to train the prospectivity model, whereas the latter is used for model assessment; (3) encoding parameters (c, ε, σ) to be optimized, and determining initial population and maximum iterations; (4) training initial population with SVR, and then calculating the fitness function; (5) if stopping criteria have been met, parameters optimization would be finished with maximum fitness function defined as preferential parameters. Otherwise, go to the next step; (6) generating a new population of chromosomes by genetic operations, including selection, crossover, and mutation, and then return steps (4) and (5).
The GA was implemented by using GATBX toolbox of MATLAB (developed by The University of Sheffield, https://www.shef.ac.uk/polopoly_fs/1.60188!/file/manual.pdf). We adopted a binary-coded genetic algorithm with stochastic universal sampling selection, single point crossover, and discrete mutation. Selective probability, crossover probability, and mutation probability were 0.9, 0.7 and 0.035 respectively. Initial population size was chosen as 20, and ranges of (c, ε, σ) were defined as [0, 100], [0.01, 1] and [0, 100] respectively. Ultimately, c = 1.17, ε = 0.01, and σ = 16.86 were fixed as optimum parameters for the construction of SVR. SVR was performed by using LIBSVM toolbox of MATLAB [68]. Furthermore, we compared GA-SVR and GS-SVR (grid search optimized SVR) via RMSE (root mean squared error), MAE (mean absolute error), R 2 (squared correlation coefficient) with known data ( Table 1). The results of this comparison showed a more robust ability of GA-SVR relative to GS-SVR for the mineral prospectivity data.

Comparison and Assessment
Known mineralization samples in the prospectivity space are thought to be an effective verification for model validity assessment [19,20,59]. To test the performance of SVR, another two typical methods, MNR and fuzzy WofE [69,70], a statistics evaluation model and a Bayesian probabilistic model were chosen as comparison methods.
The aforementioned methods were compared in terms of PA plot and ROC curve. In the P-A plot, two curves, namely the prediction-rate curve and the occupied area curve, are mapped for estimating the ability of each to predict mineralization, mineral occurrences and the sizes of predicted target areas [58]. Additionally, the point where curves intersect straightway reflects the probability of mineral deposit occurrence. When the intersection point shows a higher prediction rate and smaller occupied area, it indicates high mineralization probability. ROC curve, a validation technique widely used in machine learning, has been introduced into prospectivity models assessment [24,59,71]. It takes both the true positive rate and false positive rate into consideration. Additionally, the area under the ROC curve (AUC), with values varying from 0 to 1, can be used to measure the performance of prospectivity models. An AUC value of 1 indicates the result has perfect accuracy, whereas an AUC value below 0.5 determines the result is a random model [59,71].

Gold Distribution
Based on the gold assay data and the 3D model of orebodies, gold grade (Au) was assigned to every voxel of orebodies via the 3D Kriging interpolation method [20,72]. The 3D assay model displays a discontinuous gold spatial distribution and an obvious gold-enrichment in the footwall of

Comparison and Assessment
Known mineralization samples in the prospectivity space are thought to be an effective verification for model validity assessment [19,20,59]. To test the performance of SVR, another two typical methods, MNR and fuzzy WofE [69,70], a statistics evaluation model and a Bayesian probabilistic model were chosen as comparison methods.
The aforementioned methods were compared in terms of PA plot and ROC curve. In the P-A plot, two curves, namely the prediction-rate curve and the occupied area curve, are mapped for estimating the ability of each to predict mineralization, mineral occurrences and the sizes of predicted target areas [58]. Additionally, the point where curves intersect straightway reflects the probability of mineral deposit occurrence. When the intersection point shows a higher prediction rate and smaller occupied area, it indicates high mineralization probability. ROC curve, a validation technique widely used in machine learning, has been introduced into prospectivity models assessment [24,59,71]. It takes both the true positive rate and false positive rate into consideration. Additionally, the area under the ROC curve (AUC), with values varying from 0 to 1, can be used to measure the performance of prospectivity models. An AUC value of 1 indicates the result has perfect accuracy, whereas an AUC value below 0.5 determines the result is a random model [59,71].

Gold Distribution
Based on the gold assay data and the 3D model of orebodies, gold grade (Au) was assigned to every voxel of orebodies via the 3D Kriging interpolation method [20,72]. The 3D assay model displays a discontinuous gold spatial distribution and an obvious gold-enrichment in the footwall of the orebody in the north domain ( Figure 5). Gold mineralization with Au > 3 g/t in the north domain is mainly confined into sections of No. 48 and No. 12, where mineralization extends about 300 m from south to north. This area is a gently undulating (overall trend NNE) with a maximum thickness of orebody (about 40 m). In the northernmost orebody changes to NNW-striking, roughly paralleling with F 2 (Figure 2) and maintains a relatively high grade proximal to the fault. Industrial-grade mineralization occurs in the domain between sections No. 35

Spatial Analysis of Geological Features
The spatial distance analysis results demonstrate that the higher values (yellow to red) of distance to F2, LVP, and the upper interface of the phyllic alteration zone are similarly confined in the areas between sections No. 56 (Figure 6a-c). Gold occurrences are principally located in the hanging wall of F2 and LVP lower interface with a distance of 0 to 50 m (Figure 6a,b), and in the footwall of phyllic alteration upper interface zone with a distance of 0 to −80 m (Figure 6c). Distinct bimodal distribution has been found between the gold grade and distances (dF, dGV, and dAlt) (Figure 7(a1-c1)), showing multiple-peak relationships with gold grade. For instance, the Au-dAlt has two peaks, the flat peak is in the range from 40 to 60 with a grade peak of 3.4 g/t, whereas the sharp peak centralizes closely to 95 with a grade peak of 3 g/t (Figure 7(c1)). However, the relationship between AuMet and distances (dF, dGV, and dAlt) is similar to left-skewed Gaussian distribution with peaks respectively in 30, 40 and 50 (Figure 7(a2-c2)).
The slope extraction results show that the highest values of gF are concentrated near-surface between sections No. 36 and No. 4, but at different depths between sections No. 11 and No. 27 ( Figure  6d). The reverse is true for gAlt, which is lower at shallow depth (Figure 6f). Gold is concentrated in relatively sloped areas (>65°; Figure 7(d1-f1)). Analogously, there is a bimodal distribution between the slopes (gF, gGV, and gAlt) and grade (Figure 7(d1-f1)), such as Au vs. gF which has two peaks,

Spatial Analysis of Geological Features
The spatial distance analysis results demonstrate that the higher values (yellow to red) of distance to F 2 , LVP, and the upper interface of the phyllic alteration zone are similarly confined in the areas between sections No. 56 and No. 16,and No. 35 and No. 15. The dF, dGV and dAlt show overall lower values in the south and higher values in the north (Figure 6a-c). Gold occurrences are principally located in the hanging wall of F 2 and LVP lower interface with a distance of 0 to 50 m (Figure 6a,b), and in the footwall of phyllic alteration upper interface zone with a distance of 0 to −80 m (Figure 6c). Distinct bimodal distribution has been found between the gold grade and distances (dF, dGV, and dAlt) (Figure 7(a1-c1)), showing multiple-peak relationships with gold grade. For instance, the Au-dAlt has two peaks, the flat peak is in the range from 40 to 60 with a grade peak of 3.4 g/t, whereas the sharp peak centralizes closely to 95 with a grade peak of 3 g/t (Figure 7(c1)). However, the relationship between AuMet and distances (dF, dGV, and dAlt) is similar to left-skewed Gaussian distribution with peaks respectively in 30, 40 and 50 (Figure 7(a2-c2)).
The slope extraction results show that the highest values of gF are concentrated near-surface between sections No. 36 and No. 4, but at different depths between sections No. 11 and No. 27 (Figure 6d). The reverse is true for gAlt, which is lower at shallow depth (Figure 6f). Gold is concentrated in relatively sloped areas (>65 • ; Figure 7(d1-f1)). Analogously, there is a bimodal distribution between the slopes (gF, gGV, and gAlt) and grade (Figure 7(d1-f1)), such as Au vs. gF which has two peaks, 65 • and 85 • (Figure 7(d1)). A right-skewed distribution between the slope and AuMet demonstrates the segments with steeper slopes host a higher percentage of gold than other intervals (Figure 7(d2-f2)).
The morphological undulation analysis displays there are several distinct high-value centers of wrF, wrGV, and wrAlt in the areas between sections No. 48 to 32, No. 12 to 3, and No. 19 to 43, respectively (Figure 6g-i). Concentrated gold distribution areas of undulation (wrF, wrGV, and wrAlt) are in the interval of 0 to 20 m, −10 to 10 m, and 0 to 20 m, respectively (Figure 7(g1-i1)). The relationships between the undulation (wrF, wrGV, and wrAlt) and grade show a multimodal distribution (Figure 7(g1-i1)). Differing from the normal distribution of the other AuMet undulations, the relationship between AuMet and wrF is a bimodal distribution with two peaks of −2 and 16 (Figure 7(i2)).
Minerals 2020, 10, x FOR PEER REVIEW 11 of 21 65° and 85° (Figure 7(d1)). A right-skewed distribution between the slope and AuMet demonstrates the segments with steeper slopes host a higher percentage of gold than other intervals (Figure 7(d2-f2)). The morphological undulation analysis displays there are several distinct high-value centers of wrF, wrGV, and wrAlt in the areas between sections No. 48 to 32, No. 12 to 3, and No. 19 to 43, respectively (Figure 6g-i). Concentrated gold distribution areas of undulation (wrF, wrGV, and wrAlt) are in the interval of 0 to 20 m, −10 to 10 m, and 0 to 20 m, respectively (Figure 7(g1-i1)). The relationships between the undulation (wrF, wrGV, and wrAlt) and grade show a multimodal distribution (Figure 7(g1-i1)). Differing from the normal distribution of the other AuMet undulations, the relationship between AuMet and wrF is a bimodal distribution with two peaks of −2 and 16 (Figure 7(i2)).

The Efficiency of Prospectivity Models
The P-A plots of three prospectivity models show that GA-SVR had a prediction rate of 74% in 26% of the study area; 62% of the known ore voxels were predicted in 38% of the study area by MNR, and only 55% of the known occurrences were forecasted in 45% of the study area via fuzzy WofE (Figure 8a-c). In ROC curves, GA-SVR models with the AUC value of 0.95 exhibits a better performance of prospectivity, whereas MNR and fuzzy WofE are less accurate, with the same AUC value of 0.73 (Figure 8d). AuMet; (f) gAlt v.s. Au and AuMet; (g) wrF v.s. Au and AuMet; (h) wrGV v.s. Au and AuMet; (i) wrAlt v.s. Au and AuMet. Au stands for the gold grade, and AuMet is the gold resource summation of all ore voxels in the same factor value.

The Efficiency of Prospectivity Models
The P-A plots of three prospectivity models show that GA-SVR had a prediction rate of 74% in 26% of the study area; 62% of the known ore voxels were predicted in 38% of the study area by MNR, and only 55% of the known occurrences were forecasted in 45% of the study area via fuzzy WofE (Figure 8a-c). In ROC curves, GA-SVR models with the AUC value of 0.95 exhibits a better performance of prospectivity, whereas MNR and fuzzy WofE are less accurate, with the same AUC value of 0.73 (Figure 8d).

Target Appraisal
Based on the criterion of the highest predictive score in the smallest occupied area, three exploration targets (targets I to III) with high-value mineralization probability have been inferred in the study area using the GA-SVR prospectivity model (Figure 9). Target I is located at a depth beneath the known orebody in the north segment and in the hanging wall of the fault F2, with the burial depth of 440 m to 750 m. Close to the largest and economically productive orebody group, target I is most likely a group of deep, ore-forming fluid conduits stretching from the deep magmatic intrusion of the Axi deposit (see below discussion) and is considered as the highest priority target in the study area. Target II in the south part of the deposit occurs below the known orebodies and in the hanging wall of the fault F2 as well, at a depth of >800 m. It is situated in the transitional zone of the fault, orebodies, and alteration zone. Target III is located to north stretching farther along the N-trend of known orebodies, with a depth of >1000 m, probably demonstrating that it is the north branch of the proven orebody.

Target Appraisal
Based on the criterion of the highest predictive score in the smallest occupied area, three exploration targets (targets I to III) with high-value mineralization probability have been inferred in the study area using the GA-SVR prospectivity model (Figure 9). Target I is located at a depth beneath the known orebody in the north segment and in the hanging wall of the fault F 2 , with the burial depth of 440 m to 750 m. Close to the largest and economically productive orebody group, target I is most likely a group of deep, ore-forming fluid conduits stretching from the deep magmatic intrusion of the Axi deposit (see below discussion) and is considered as the highest priority target in the study area. Target II in the south part of the deposit occurs below the known orebodies and in the hanging wall of the fault F 2 as well, at a depth of >800 m. It is situated in the transitional zone of the fault, orebodies, and alteration zone. Target III is located to north stretching farther along the N-trend of known orebodies, with a depth of >1000 m, probably demonstrating that it is the north branch of the proven orebody.

Insights from 3D Models and Spatial Analysis
Our spatial analysis shows that there is a bimodal distribution of gold grade (Au) and a single peak of gold resource (AuMet) in the pyrite-sericite-quartz alteration buffer field (Figure 7(c1)). Here, dAlt (40,60) may represent the position of main orebodies, because of the higher grades (avg. up to 3.5 g/t) and approximate distance to width (about 20 m). The intervals of dAlt (0, 20) and dAlt (20,40) most likely reflect the location of disseminated ore and the mixing zone of disseminated ore and quartz-sulfide stockwork, respectively, taking account of the value and increasing rate of Au and AuMet. Such vein-alteration distribution is well supported by field investigations (Figure 2). It is very interesting that almost all siliceous veins with gold mineralization are hosted in the phyllic alteration zone, with major vein-type orebodies located in their centers. Apparently, phyllic alteration, to some extent, controls the occurrence of late quartz veins that are products of distributed hydrothermal fluid infiltrating a brecciated zone [33]. However, the crosscutting relationship and 20 M.y. time gap between phyllic alteration and auriferous quartz veins preclude the coeval formation in one hydrothermal event [33,35,37,39]. Previous studies have proposed that the interactions between hydrothermal fluids and rocks could enhance the rock permeability by mass and volume changes linked to mineral dissolution and precipitation, and further control or facilitate the brecciation and formation of new fractures [73][74][75][76]. Since the fault F2 underwent at least three-phases of deformation during strike-slip movement during ore-forming processes [41], the occurrence of quartz veins was mostly likely promoted by late structural deformation on the earlier altered rock.
Two major types of economic ore (i.e., disseminated and vein-type) at Axi are genetically related to early fluid-rock reactions during phyllic alteration and later fluid mixing and boiling in the fault damage zone [33,35,37,39]. The above mineralization mechanism is ubiquitous in LS epithermal gold deposits and is essentially attributed to fluid flow and convergence in the epithermal environment [2,3,6]. The 3D models and spatial analysis indicate that the morphological characteristics and distance field of geological features significantly control the distribution of gold resources in the Axi gold deposit. Specifically, the AuMet is commonly confined to the zones distance to F2 (0-50 m), the LVP lower interface (10-60 m), and the phyllic upper interface (20-70 m), with high dipping angles of 65-75° and gentle undulation, but with increased frequency on convex areas (0-20 m; Figure 7). These characteristics underline the critical roles of these geological interfaces for driving and gathering ore-forming fluids at Axi. It is notable that all convex areas have high grade and resources (Figure 7(g1,g2)). This implies strain localization (or release) with orientation similar to global strain,

Insights from 3D Models and Spatial Analysis
Our spatial analysis shows that there is a bimodal distribution of gold grade (Au) and a single peak of gold resource (AuMet) in the pyrite-sericite-quartz alteration buffer field (Figure 7(c1)). Here, dAlt (40, 60) may represent the position of main orebodies, because of the higher grades (avg. up to 3.5 g/t) and approximate distance to width (about 20 m). The intervals of dAlt (0, 20) and dAlt (20,40) most likely reflect the location of disseminated ore and the mixing zone of disseminated ore and quartz-sulfide stockwork, respectively, taking account of the value and increasing rate of Au and AuMet. Such vein-alteration distribution is well supported by field investigations (Figure 2). It is very interesting that almost all siliceous veins with gold mineralization are hosted in the phyllic alteration zone, with major vein-type orebodies located in their centers. Apparently, phyllic alteration, to some extent, controls the occurrence of late quartz veins that are products of distributed hydrothermal fluid infiltrating a brecciated zone [33]. However, the crosscutting relationship and 20 M.y. time gap between phyllic alteration and auriferous quartz veins preclude the coeval formation in one hydrothermal event [33,35,37,39]. Previous studies have proposed that the interactions between hydrothermal fluids and rocks could enhance the rock permeability by mass and volume changes linked to mineral dissolution and precipitation, and further control or facilitate the brecciation and formation of new fractures [73][74][75][76]. Since the fault F 2 underwent at least three-phases of deformation during strike-slip movement during ore-forming processes [41], the occurrence of quartz veins was mostly likely promoted by late structural deformation on the earlier altered rock.
Two major types of economic ore (i.e., disseminated and vein-type) at Axi are genetically related to early fluid-rock reactions during phyllic alteration and later fluid mixing and boiling in the fault damage zone [33,35,37,39]. The above mineralization mechanism is ubiquitous in LS epithermal gold deposits and is essentially attributed to fluid flow and convergence in the epithermal environment [2,3,6]. The 3D models and spatial analysis indicate that the morphological characteristics and distance field of geological features significantly control the distribution of gold resources in the Axi gold deposit. Specifically, the AuMet is commonly confined to the zones distance to F 2 (0-50 m), the LVP lower interface (10-60 m), and the phyllic upper interface (20-70 m), with high dipping angles of 65-75 • and gentle undulation, but with increased frequency on convex areas (0-20 m; Figure 7). These characteristics underline the critical roles of these geological interfaces for driving and gathering ore-forming fluids at Axi. It is notable that all convex areas have high grade and resources (Figure 7(g1,g2)). This implies strain localization (or release) with orientation similar to global strain, accompanied gold deposition, and likely led to the brecciation or caused fluid preferential confluence in more outward bending areas at Axi.
The 3D models and distance field analysis reveal a rule that the alteration-mineralization zone in bending sites is clearly wider than the overall distribution (Figure 3), and notably, also have relatively high gold grades (Figure 7(c1)). Interestingly, these are locations where the main orebodies are comprised of quartz vein-type ores, suggesting at least a rough spatial or genetic correlation between vein-type orebodies and normal-strike slip extension deformation. The fragmentized contact zone ascribed to the overprinting structural deformation on the early altered zone may be the pathway and dilation space for hydrothermal fluid transmission [3,75,76]. In addition, the altered rocks near permeable channels commonly have different mechanical properties from nearby wall rocks, thereby showing distinct deformation during fault propagation, as discussed above. Thus, the quartz vein orebodies in the transition zone reflect the original location of fluid conduits, which is consistent with the empirical mineral assemblage of quartz/chalcedony-carbonates-adularia-sericite as the indicator for fluid conduits in the epithermal system [3,8,77,78]. We, therefore, infer that the major fluid conduits of the Axi epithermal gold deposit are located in the bending areas (No. .

Application of Support Vector Regression Model
The extracted metallogenic information in 3D prospectivity modeling of the Axi epithermal gold deposit has the following features: (1) non-linearity in that gold mineralization is spatially controlled by a series of geological factors with non-linear relationships [14,49,59]; (2) distribution diversification in that the ore-controlling indicators datasets exhibit different distribution statistics [23][24][25]; (3) high-dimensionality in that 9 predictor maps construct a nine-dimensional feature space [23][24][25].
In this study, the MNR method for integrating ore-forming information from the Axi epithermal gold deposit shows weak prospective ability. MNR has been proved fit for prospectivity modeling of deposits with a relatively mild grade and amount of metal variation, distributed normally [23,26]. However, the gold grade of Axi gold deposit changes sharply and yields a log-normal distribution ( Figure 6), hence, MNR is a poor fit for prospectivity modeling of this LS epithermal gold deposit, and especially performs worse in fitting samples with erratic high-grade. Although the fuzzy WofE identified more known ore voxels in a smaller occupied space when compared to MNR, it did not produce a robust prospectivity model. In fuzzy WofE, spatially continuous predictor map values are first simplified and discretized into a series of arbitrary classes, and the same weight is assigned to all values in each class. This process is excessively dependent on expert knowledge, resulting in some geological information loss or inconsistencies, which would generate systemic bias in the weighting of spatial evidential layers in the prospectivity model [20,79]. Furthermore, multiple layer integration based on fuzzy WofE is analogous to a nonlinear multiplicative cascade process. Thus, the accuracy and reliability of predictive results may decline, if there is dependence or weak conditional independence between diverse predictive maps [57].
Both the PA plot and ROC curve indicated GA-SVR was superior to MNR and fuzzy WofE in delimiting exploration targets in the Axi gold deposit. These results could be interpreted in several ways: (1) As previously mentioned, the known Au-grade in Axi gold deposit changes sharply and unevenly distributes in space, resulting in the poor performance of traditional prospectivity methods (e.g., MNR and fuzzy WofE). Nevertheless, GA-SVR is able to overcome the aforesaid disadvantages, because of it is irrespective of data distribution and outliers, making better generalization performance [30,32,65,66]. (2) The spatial analysis results indicate that there is a complicated non-linear spatial association between gold occurrence and geological features, due to the superimposition of various geological events in the ore-forming processes that are difficult to quantify [57]. SVR, characterized by non-linearity, self-learning, and robustness, constructs a higher dimensional feature space from the lower dimensional input space via non-linear mapping [20,50,51]. Then, the complicated non-linear spatial association in the original low dimensional input space is treated as an oversimplified linear regression solution in high dimensional feature space, and therefore it can be commendably quantified and the prediction ability is improved [54,62]. (3) There are many parameters that have to be determined using expert opinion in traditional prospectivity models (e.g., weight definition in fuzzy WofE), resulting in systemic bias and error [20,79]. By contrast, GA-SVR obtains optimal parameters that could produce the best prediction via GA with less subjectivity [30,32,61].
Consequently, the hybrid model of GA-SVR based on the combination of SVR and GA takes advantage both of GA and SVR in the prospecting modeling of the Axi gold deposit, and therefore is a promising method in 3D mineral prospectivity modeling.

Exploration Significanture for LS Epithermal Deposit
LS epithermal mineralization is typically developed at shallow depths in terms of the high-permeability networks favoring magmatic water ascent, external fluid influx, and the interplay of fluid and rock, thereby resulting in the extensional vein systems and intensive water-rock interactions [1][2][3][4][5][6][7]. These characteristics are clear in our 3D models and spatial analysis results and allow understanding of the metallogenic implications in the Axi LS epithermal gold deposit. In addition, the integration of obtained spatial data about F 2 , LVP and the phyllic alteration zone controls on mineralization was used to identify the fertility of unknown areas at Axi. One of the targets in the northern Axi district is buried in the Aqialehe sedimentary formation and indicates ore-forming fluid migration along F 2 from conduits to the north in the shallow areas. This suggests a major horizontal fluid pathway controlled by structures in the epizonal segments of this LS epithermal deposit. Thus, in the exploration of epithermal mineralization in a mine, more attention should be given to the (sub-) horizontal extending directions from present orebodies with particular concern paid to the local structural deformation.
Although deep segments of ore fluid conduits host less gold in epithermal deposits, the areas could offer clues to find a new type of mineralization, such as base metal mineralization and porphyry [80,81]. Hence, determining as much as possible about the LS epithermal deposit fluid conduits, even at depth, has the potential to be quite valuable in guiding further mineral exploration. In this study, a comprehensive analysis of geological models and spatial analysis obtained important information about the ore fluid conduits in the Axi epithermal system that include the relative specific spatial location, meaningful geological indicators, and overprinting structural influences on up-flow zones. The above findings highlight the great potential of 3D geological modeling and spatial analysis in understanding the overprinting structural-hydrothermal processes and determining the location of fluid conduits.
In summary, the framework of 3D prospectivity modeling presented here, including 3D geological modeling, spatial analysis, and prospectivity modeling, can improve the understanding of ore genesis and allow for the improvement of exploration strategies in LS epithermal deposits. We encourage all LS mine operators and prospectors to develop similar models.

Conclusions
(1) The 3D prospectivity modeling for the Axi gold deposit quantitatively determines the spatial association between orebodies and geological features, and further improves evaluation of the mineralization potential in 3D space at the deposit-scale, and can aid in developing and enhancing mineral exploration strategies for the LS epithermal gold deposit.
(2) The hybrid model GA-SVR is proven to be better at quantifying the complicated non-linear spatial association between gold occurrence and geological features, with less subjectivity when compared to MNR and fuzzy WofE, demonstrating a better prediction capability in prospectivity mapping.