In mineral exploration projects, drilling is the most important, yet expensive discovery tool. Drilling locations are usually selected based on surface anomalies. A surface anomaly may be achieved from a certain geological, geophysical or geochemical exploration dataset. Drilling location can also be identified from integrated analysis of various exploration layers. Apart from integration methods, modeling methods between surface data or subsurface data to identify borehole locations are also introduced and used in this research. This proposed method could increase the probability of interesting the offered boreholes with subsurface mineralization. Sari-Gunay epithermal gold exploration datasets are used to test the proposed model. Sari-Gunaiy deposit is located in the southwest of Kordestan province, NW Iran. There are two epithermal gold mineralization zones at Sari-Gunay and Agh-Dagh hills. Detailed exploration has been completed at the Sari-Gunay hill, but the exploration at the Agh-Dagh hill is still in its preliminary stage. Epithermal gold mineralization in these two hills is low sulfidation in character. The mineralization is mostly hosted by fragmental/brecciated sub-volcanic rocks of the Miocene age. The dominant lithological units of study area are dacite and andesite porphyries, diatreme tuff breccia and sub-alkaline lithic porphyry rocks. The Sari-Gunay hill with a resource of 120-130 Mt at 2 g/t Au is the largest ever discovered gold deposit in the Middle East. In this research, quantity and quality modelling methods are used to model subsurface gold mineralization. In the quantity modelling, three mineralization parameters such as Mean, Sum and PI were modeled continuously by multivariate regression (MR) and support vector regression (SVR) methods. The PI parameter was also modelled discontinuously by discriminant analysis (DA) and support vector machine (SVM) methods. The results obtained from continuous and discontinuous modelling identified four important gold mineralized zones at the Sari-Gunay hill, Agh-Dagh hill, a faulted zone between Sari-Gunay and Agh-Dagh hills and in the northeast of the Sari-Gunay hill. The RBF-nu-SVR function produced the best results for all the four identified mineralized zones in continuous modelling, whereas in discontinuous modelling, RBF-nu-SVM function produced the best results for Sari-Gunay and Agh-Dagh hills and the LDA function for the faulted zone. The achieved results indicated that the functions of the SVR and SVM methods, because of minimizing the structural risk, were produced better modelling results than the functions of the MR and DA methods with less empirical risk. The maps obtained by Mean parameter and Sum and PI parameters could be used to design the drilling area and the density of drilling grids, respectively. In quality modelling methods, the combination of geological data and surface and subsurface geochemical data were used to model the subsurface gold mineralization. Fractal, positive and negative geochemical anomalies and geochemical zoning methods were used for this purpose. The results of the fractal modelling methods showed that blocks with very high gold grades were located in the brecciaed andesitic/dacitic rocks, while most of the mineralization occurred in the dacite porphyry rocks. The dacitic tuffs were ranked third in terms of gold mineralization. Fractal modelling at the Sari-Gunay hill indicated that the economic gold mineralized zone is ellipsoidal in shape (1050m ×650m ×350m). The results of the positive and negative geochemical anomaly methods showed that chalcophile elements related to the mineralization have positive anomalies, whereas the lithophile and siderophile elements revealed negative anomalies in the study area. Hydrothermal solutions and the geochemical characteristics of the elements played the most important role in the formation of the positive anomalies. Also hydrothermal alteration along with weathering processes were more important in the formation of the negative anomalies. The portion of the geochemical holes of the elements with positive anomalies to the geochemical holes of the elements with negative anomalies were modeled for all the four the mineralized zones. The results of geochemical zoning method showed that in the series of the elements obtained by surface axis, Au was a middle element and Ba, Ce, Cu, Hg, La, Pb, Sc and Sb elements and B, Mn, Mo and Zn elements were supra-ore and sub-ore elements, respectively. The compositional holes index indicated that the erosion level at the Sari-Gunay hill was 50-60 meters higher than the Agh-Dagh hill. The geochemical resources, obtained from the geochemical zoning, at the Sari-Gunay hill was at least 22.5 times higher than the Agh-Dagh hill.