Several machine learning approaches have been developed for the identification of geochemical populations. In these approaches, the geochemical elements are usually the sole quantitative variables used as inputs for geochemical population recognition. This means that the presence of other qualitative variables, such as geological information, is overlooked in the analysis. Hierarchical clustering, as an unsupervised machine learning method, is a common approach for dimensional reduction in the analysis of geochemical data. In this study, an alternative to this technique, known as geostatistical hierarchical clustering (GHC), is applied to identify geochemical populations in 3D in the Bondar Hanza copper porphyry deposit, Iran. In this paradigm, the qualitative geological variables can also be incorporated for geochemical population identification, in addition to qualitative geochemical elements. In this study, an innovative solution is presented to tune the weighting parameters of each variable in GHC, based on the associations that the clusters (i.e., geochemical populations) should have with the geological information. The results are compared with k-means and number-size fractal/ multifractal (N-S) methods. As a result, GHC showed better agreement with alterations, rock types, and mineralization zones in this deposit. Finally, some important instructions are provided for further mineral exploration.