## Abstract

The importance of modeling geological units and mineralization factor, as a key step in evaluating mineral deposits, is undeniable. Multivariate geostatistical methods are very useful tools in the case of geochemical modeling. In this study, the minimum/maximum autocorrelation factor (MAF), a multivariate geostatistical method, and the sequential indicator simulation (SIS) were used in order to model the Sury Gunay epithermal gold deposit, NW of Iran. The analyses were done using core samples obtained from drill holes. The MAF is a geospatially modified version of principal component analysis, which decorrelates factors for all lag distances. By applying the MAF on alr-transformed data, four MAF were obtained and the fourth factor represents the mineralization factor. Joint simulation of the mineralization elements was carried out by applying the sequential Gaussian simulation to the fourth factor scores. The main rock types and alterations of the deposit were also modeled using SIS and the probabilistic models of the rock types, and alterations were obtained using E-type maps resulting from 20 realizations. The results indicate that there are better relations with higher values of the mineralization factor and the existence of the volcanogenic breccia and dacite porphyry rock types. Furthermore, the simulated alterations demonstrate that the higher probability of silicification existence may have better correlation with higher concentrations of the mineralization elements.

Modeling of mineralization controlling factors and geological units such as rock types, alterations, and mineralogy is crucial for the understanding and evaluation of mineral deposits. This can be done by either geological or computational approaches. Computational modeling, which is the method for making acceptable predictions, has become common in geoscience over the past two decades (Hobbs *et al.* 2000; Zhao *et al.* 2008, 2009; Liu *et al.* 2012; Fouedjio & Séguret 2016; Khan & Deutsch 2016; Mooney & Boisvert 2016). Geostatistical simulation and estimation algorithms are the most important approaches for computational modeling. Over the last decade, in order to model the geological features of a deposit – such as rock types or facies in petroleum reservoir – geostatistical simulation algorithms have been widely used (Modis & Sideri 2013; Talebi *et al.* 2015; Fouedjio & Séguret 2016; Moniz Caixeta *et al.* 2016; Musafer *et al.* 2016). Different simulation methods exist for continuous and categorical variables. One of the most common methods of simulating the categorical variables is Sequential Indicator Simulation (SIS) (Alabert 1987; Journel & Gomez-Hernandez 1993). This method can be used in the absence of large-scale curvilinear features (Yin 2013; Nejadi *et al.* 2015) and is suitable for incorporating uncertainty that propagates through to the resulting numerical models (Deutsch 2006; Olea & Luppens 2012).

Furthermore, statistical analysis has a significant role in mineral exploration. The main aim of methods of multivariate statistical analysis – such as principal component analysis (PCA) and factor analysis (FA) – is to reduce the dimensionality of the data and to derive the related mineralization factors. The minimum/maximum autocorrelation factor (MAF) is a modified version of PCA, which was first developed by Switzer & Green (1984). Later, Desbarats & Dimitrakopoulos (2000) used this method in geostatistical fields. The advantage of the MAF over PCA is that it spatially orthogonalizes the variables into uncorrelated factors for all lag spacings (Bandarian *et al.* 2008). In this method, coregionalization matrices are produced by computing the two-structure linear model of coregionalization (Bandarian *et al.* 2008; Shakiba *et al.* 2015). The MAF can be employed for grade estimations by decorrelating the variables and estimating each one independently, thus obviating the need to use cokriging (Da Silva & Costa 2014). Univariate simulation may not be suitable because of its flaws in ignoring the related variables. So, the MAF can be applied for joint simulation of variables (Tajvidi *et al*. 2015). In other research, Sohrabian & Tercan (2014) used the MAF for joint simulation of geophysical attributes. Boucher & Dimitrakopoulos (2009) also applied the MAF to orthogonalize the factors before conditional block simulation of non-Guassian variables via Lu simulation.

However, applying the MAF is problematic due to the compositional nature of the geochemical data. Spurious correlations between variables are a major problem existing in raw geochemical data because of the constant sum constraint. In other words, this restriction does not make the variables independent such that increasing the proportion of one variable causes the proportions of other variables to decrease. Therefore, correlations are not free to vary from −1 to +1 (Pawlowsky-Glahn & Egozcue 2006). So, applying some multivariate statistical methods, which are based on the correlation coefficients among variables, to untransformed data may yield unreliable results. Therefore, such closed system data should be opened by logratio transformations prior to further analysis (Aitchison 1983, 1986, 1999; Aitchison *et al.* 2000; Buccianti & Pawlowsky-Glahn 2005; Buccianti & Grunsky 2014). To achieve this purpose, three types of logratio transformations – namely additive logratio (**alr**) transformation (Aitchison 1983), centered logratio (**clr**) transformation (Aitchison 1983), and isometric logratio (**ilr**) transformation (Egozcue *et al.* 2003) – have been introduced. Many research papers in geochemical exploration have considered the compositional nature of geochemical data. For example, three types of logratio transformation have been applied to stream sediment geochemical data by Carranza (2011) for mapping and analyzing geochemical anomalies. To identify pathfinder elements anomalies using multivariate analysis, de Caritat & Grunsky (2013) have performed both alr and clr transformations to soil geochemical data in Australia. Moreover, Savu-Krohn *et al.* (2011) used ilr transformed data to model the geochemical fingerprint of ores by machine learning. Mueller & Grunsky (2016) applied the MAF to clr-transformed data for better understanding of the underlying geology through the recognition of factors representing lithological features and glacial transport processes. Another example application of the MAF to compositional data has been done by Mueller *et al.* (2014) whereby they performed MAF to alr-transformed data and concluded that alr-normal score transform-MAF yielded better results against a full Guassian co-simulation of the raw data.

The main aim of this study is to model rock type and alteration as key mineralization controller factors and to find relations between the mineralization and the rock types and alterations. In this paper, the MAF method was performed on alr-transformed data to produce and model the related mineralization factor for gold and relevant elements. Then, SIS was used for modeling rock types and alterations. Finally, the results were linked with each other, which define the geological modeling of the deposit.

## Case study

### Geological setting

The Sury Gunay epithermal gold deposit is located 60 km NW of Hamedan, in the SE corner of Kurdestan province in NW Iran (Fig. 1). Initial investigation has been conducted in May and August 1999 by RioTinto Mining and Exploration Limited (Wilkinson 2005). The estimated resource of the deposit was about 52 Mt oxide resource at an average grade of 1.77 g/t gold and it could be assumed as the biggest gold deposit in Iran (Kouhestani *et al.* 2012; Ghorbani 2013). Recently, Geranian *et al.* (2016) identified new potential areas for further drilling in the Sury Gunay deposit, and so its resources are probably more than what is reported here.

The Sury Gunay epithermal deposit is situated at the Sanandaj-Sirgan metamorphic belt of Jurassic to upper Cretaceous age. There are two main parallel trending zones in western Iran namely Sanandaj-Sirjan belt and Uromieh-Dokhtar volcanic belt with Palaeocene and Pleistocene age. The Sury Gunay deposit is situated at the short volcanic belt between and sub-parallel to the Sanandaj-Sirjan and Uromieh-Dokhtar belts, named Ghorveh-Bijar belt. The main properties of this belt are discrete volcanic centers and dacites and Pleistocene basalts with a set of large NW–SE-trending structures (Wilkinson 2005).

The geology of the Sury Gunay deposit is related to a high-K sub alkaline volcanic complex of middle Miocene age. There is a large area of phyllic and argillic alteration in a 16 km^{2} area around the vein system, which includes considerable amounts of gold mineralization in Sury Gunay. The deposit is classified as a low sulfidation epithermal deposit, which is thought to represent the eroded roots of a large central volcano (Richards *et al.* 2006). The major host rock of the Sury Gunay is a series of dacite porphyries including quartz dacite porphyry, dacite lithic tuff – as a part of the dome complex which is transected by vents of tuff – and breccia (Fig. 1). The main mineralization of the Sury Gunay consists of Au, Sb, Hg mineralization surrounded by As, Pb, Zn haloes. The large hydrothermal alteration system that exists in the area mainly consists of sericitic and silicification alterations (Wilkinson 2005; Richards *et al.* 2006).

### Dataset

In this study, a dataset containing 26 347 samples collected from 86 drill holes was used (Fig. 2). The drillhole samples collected at intervals of 1 meter have been analyzed for 46 elements using inductively coupled plasma-optical emission spectrometry analytical method and for gold by fire assay method in the OMAC laboratory in Ireland. In addition, the drill holes have been logged with information on mineralogy, rock types, and alterations (Jackson *et al.* 2006). Because of the high variety of rock types and alterations in the area, some of the rock types and the alteration that are more abundant and important have been selected for this study. These are, respectively, dacite porphyry (DP), quartz dacite porphyry (QDP), dacite lithic tuff (DLT), and volcanogenic breccia (VBX) as rock types and silicification, sericitization, and tourmaline as alteration signatures. The descriptive statistics of the Au per rock types and alterations show that the total mean and the variance of the Au data are, respectively, 1.35 ppm and 9.66 ppm. The maximum values of Au are related to VBX and silicification among the rock types and alterations, which are, respectively, 1.67 ppm and 1.59 ppm (Table 1). Besides, it should be mentioned that the outlier values are replaced using Q-Q plot method.

## Methodology

### Logratio transformation

The constant sum constraint is the main problem existing in compositional data (e.g. raw geochemical data). These data are defined as:
(1)where the variables are positive. Due to this constraint, an increase in the values of an element is neutralized by a decrease in the values of the other elements. Therefore, the correlation coefficients between the elements are not reliable, and are negatively biased (Egozcue & Pawlowsky-Glahn 2005; Thomas & Aitchison 2006). Therefore, applying statistical methods such as PCA and MAF to raw compositional data may yield unreliable results (Aitchison 1983). Thus, logratio transformations should be applied to open the closed system data (Carranza 2011; Zuo *et al.* 2013; Wang *et al.* 2014; Zuo 2014). Among the methods of logratio transformations, alr transformation was used in this study before applying the MAF. The alr transformation can be shown as (Aitchison 1982):
(2)where *x*_{D} is denominator with important features. Firstly, the denominator must be the same for all variables and it should be strictly positive for all components (Job 2012). The denominator should be chosen among the variables which do not have important roles in the study. Furthermore, after applying the alr transformation, the denominator used is excluded from further analysis, and only the transformed components are used (Job 2012). Moreover, in some cases, analysis is done using a few numbers of the components not all the composition. These components construct a subcomposition which is a subset extracted from a full composition (Job 2012). Therefore, the total sum of the subcomposition components does not satisfy the constant sum. Thus, a filler component is defined for a subcomposition which is K-sum of all components.

### Minimum/Maximum Autocorrelation Factors (MAF)

As mentioned before, the MAF is a modified version of PCA whereby the produced factors are decorrelated not only for lag zero (*h *= 0) but also for all lag distances. In other words, MAF factors are decorrelated spatially (in whole sample spaces). Spatial decorrelation of the MAF factors are done using two-structure linear model of coregionalization (LMC). The first step of the MAF method is to transform the variable distribution into normal scores (zero mean and unit variance). Then, the important phase is to decorrelate the factors, which is done by three key steps. The first step is decorrelation using PCA, which can be done by using spectral decomposition of the variance-covariance matrix *B* (the overall sill values of the variograms and equivalent to correlation coefficients matrix of the variables) of the variables as : where is a diagonal matrix of the eigenvalues and *Q* is orthogonal matrix of the corresponding eigenvectors. In this step, PC's can be calculated using matrix , matrix *Q* and the normal scores of the variables (*Z*_{(u)}). However, they are not decorrelated in space yet. In the next step, from the two-structure LMC (Equation 3), matrix *B _{1}* (the sill values of the first structures) is achievable (Rondon 2012):
(3)where the functions and are the nested variogram models. Then, in order to produce the spatially decorrelated factors, matrix

*V*can be defined as , where . Matrix

*V*can be orthogonized as , with as eigenvectors. Consequently, the MAF can be derived using eigenvectors of

*V*by defining (Rondon 2012): (4)where

*Z*(

*u*) is an N-dimensional random field with Gaussian normal distribution. After obtaining the MAF, each factor can be simulated independently using univariate variogram of the MAF scores and simulation methods.

Let represent the simulated MAF scores, then the back-transformation of the simulated scores to the simulated attributes are computed as (Rondon 2012):
(5)The key step for interpreting the MAF results is analyzing the *M*^{−1} transformation matrix. In this matrix, for every variable, each factor has its own weights. These weights demonstrate the impact of each factor towards the corresponding variables. Those factors having high absolute values will have the strongest impact on the variables. The strength of the MAF is reproducing well the statistical features such as univariate statistics, scatter plot or spatial continuity (Boucher 2003; Dimitrakopoulos & Fonseca 2003; Boucher *et al.* 2005; Boucher & Dimitrakopoulos 2009).

## Results and discussion

The main purpose of performing the MAF in this study is to model the gold and related elements as the mineralization factor, simultaneously. Based on the geological report of Sury Gunay deposit by Wilkinson (2005), geochemical studies of soil, rock chip and drill cores revealed that the main mineralization system of Sury Gunay is Au-As-Sb-Hg-Tl ± Ag zoned epithermal system. Among the mentioned elements, Tl and Ag are omitted from the dataset because of the high percentage of censored data; Tl has 72% and Ag has 63% of censored data. Therefore, four elements (Au, As, Hg, Sb) are selected among all measured elements for applying the MAF. In the first step, all censored data existing in the dataset should be imputed prior to further analysis. For this purpose, a method called alr-EM (Palarea-Albaladejo & Martín-Fernández 2008) was utilized. The advantage of this method is that it takes into account the compositional nature of geochemical data. Next, alr transformation was applied to the imputed dataset using a filler component (10^{6}-sum) as the denominator to eliminate the constant sum constraint in the data and transform the variables to an open system. Therefore, a classical multivariate geostatistical frameworks could be applied to compositions following the principle of working on alr-transformed scores (Pawlowsky-Glahn & Olea 2004). By using the alr method, one could guarantee the unbiasedness of the transformed data (Ward & Mueller 2012).

### Application of the MAF

In this study the MAF was used for multivariate mapping of the variables. It should be applied on standard normal distribution data. Thus, normal score transformation was applied to alr-transformed data, at the first step. However, this method transforms the variables separately and the resulting normal scores are not multivariate normal (van den Boogaart *et al.* 2016). Recently, van den Boogaart *et al.* (2016) proposed an affine-equivariant multivariate normal score transformation method, called flow anamorphosis, for compositional data that offers a solution for the problem. Because of the mathematical complexity of this method, the authors ignored this method and simply have used normal score transformation. For the next step, omnidirectional two-structured LMC were computed for all pairs. Both structures are spherical with the range equal to 5 m for the first structure and the range equal to 350 m for the second structure. These ranges were determined based on the distances between the drillholes and best fit for all pairs. The LMC model is presented as:
(6)Matrix *B _{1}* and

*B*correspond to the following fitted LMCs: where

_{2}*B*

_{1}and

*B*

_{2}are matrices of sill values of the first and second structures of the variograms, respectively. Note that the columns and rows correspond respectively to As, Au, Hg and Sb. To produce the MAF, PCA should be performed two times. It means that the MAF are created by further rotating of the PCA factors (Shakiba

*et al.*2015). The important point after obtaining the MAF is to check if the factors are decorrelated for all lag spacings or not. To explore this issue, correlograms between the factors obtained by the MAF approach were calculated (Fig. 3). As shown in Figure 3, all factors are well decorrelated spatially so that the MAF results are reliable and further analysis can be performed. Also, the histograms of all factors should have normal distribution after applying the MAF. The cumulative distribution functions of the four factors have normal distributions with mean 0 and variance 1 (Fig. 4). The eigenvalues corresponding to the factors are shown in Table 2.

Once the decorrelation between the MAF is guaranteed, the factors can be interpreted using matrix *M*^{−1}. Matrix *M*^{−1} is as follows:

As shown in matrix *M*^{−1}, factor 4 demonstrates strong correlations of the same signs among all elements, which may be interpreted as the mineralization factor. As a result, mapping the factor 4 scores can show multivariate anomaly of the mineralization. Mapping of the mineralization factor was done, in this study, using a geostatistical method named sequential Gaussian simulation (SGS). This method was selected because it provides more accurate results compared to kriging method. Kriging is a linear geostatistical interpolation method that has the disadvantage of smoothing the effects. This method is not able to reproduce the original distribution of the data, especially when the data are highly skewed. In contrast, SGS yields more precise results (Deutsch & Journel 1998; Sadeghi *et al.* 2015; Deutsch *et al.* 2016; Khan & Deutsch 2016; Soleimani *et al.* 2016). For applying the SGS, a block model of the deposit was created with the grid size of 10 × 10 × 5 m. All simulations of the mineralization factor, rock types, and alterations were done in this block model. Moreover, variograms of the factor 4 were calculated considering the anisotropy of the mineralization factor. The major variogram was calculated for azimuth 30 with dip of 50 and the range of 350 m and the minor variogram for azimuth 210, dip 40 and range of 220 m. Simulation was performed using the mentioned variograms; and the multivariate anomaly map of the deposit is plotted (Fig. 5). The simulated values of the mineralization factor are multiplied by −1 in order to better demonstrate the mineralization factor anomalies. The relations between rock types, alterations and mineralization factor can be investigated from the multivariate anomaly map.

### Modeling of rock types and alterations

The SIS was used in this study for modeling the categorical variables like rock types and alterations. As mentioned earlier, some of the more abundant and more important rock types and alterations were chosen to be modeled among all types. The overall percentage of each rock type is thus: dacite porphyry (DP) 32%, quartz dacite porphyry (QDP) 24%, dacite lithic tuff (DLT) 14% and volcanogenic breccia (VBX) 5%. Also, silicification, sericitization, and tourmaline were selected as mineralization bearing alterations that are very common in the study area and have been reported in most drillhole samples. Thus, in order to investigate a reliable relation between the mentioned alterations and the mineralization factor, only samples with strong to intense alterations were used for modeling. In the case of tourmaline alteration, samples with medium to intense alteration were used. The total percentage of each alteration is thus: sericitic 74%, silicification 43% and tourmaline 32%. This means that more than one alteration was detected in most samples. The SIS was performed 20 times for each variable and the E-type of the 20 realizations was calculated in order to produce the probabilistic model of the rock types and alterations based on the most frequent occurrences (Figs 6 and 7).

The model of the rock types reveals that the most dominant rock type of the deposit is dacite porphyry. Dacite lithic tuff occurs mostly in southern parts of the deposit, at intermediate and high depths. Quartz dacite porphyry mostly occurs in the south parts of the deposit from the lowest to the highest depth. Volcanogenic breccia constitutes a minority of the rock types and mainly occurs in the SW part of the deposit (Fig. 6). Visual comparison of the rock types and mineralization factor models suggested that where the rock type is VBX, it seems that the mineralization factor shows high values (Figs 5 and 6). Also, it is worth mentioning that to those blocks with the same E-type values, no rock type was assigned.

The alteration maps reveal that the most dominant alteration of the deposit is sericitic, which exists in most parts of the deposit. Silicification is the second most dominant alteration of the deposit, which is mostly dispersed in the south parts, whereas tourmaline represents minor alteration, which is mainly concentrated in SW part of the deposit (Fig. 7). In the case of alterations, it is not possible to present an overall alteration model. This is because many alterations can exist in a single block.

Furthermore, to investigate the relations between the rock types, alterations and the mineralization factor in detail, the averages of the mineralization factor were calculated separately for each rock type (Table 3). Moreover, in the case of alterations, an E-type threshold value was chosen for different alterations, in which the blocks with E-type values exceeding the threshold (specified for each alteration) were considered as the corresponding alterations. The thresholds for sericitic, silicification and tourmaline alterations are, respectively, 0.7, 0.55, and 0.5, which were chosen based on the validation of the alterations models by different thresholds of drillhole data. Then, the mean values of the mineralization factor were computed for each alteration (Table 3).

Table 3 shows that the highest mean values of the mineralization factor among the rock types correspond to VBX, DP, DLT, and QDP. Consequently, in areas where VBX and DP exist, concentrations of the mineralization elements are expected to be higher. Also, the highest mean values of the mineralization factor among the alterations correspond to the sericitic, silicification and tourmaline. It seems that the very high frequency of sericitic alteration in the area caused the mean value of the mineralization factor to become higher compared to other alterations. In other words, sericitic alteration exists in both rich and depleted areas of the deposit so that the exact relation between sericitic alteration and the mineralization factor cannot be investigated in the model. The mean value of the mineralization factor in silicification is slightly less than the mean value of sericitic alteration (Table 3), but silicification has the highest Au mean value in the drillhole data. Thus, it can be concluded that silicification is a good sign of mineralization in the deposit.

## Conclusion

Geostatistical modeling of the mineralization controller factors is a useful tool for evaluating epithermal gold deposits. In this study, simulation methods were applied for modeling the categorical (rock types and alterations) and continuous variables (grade). Also, the MAF was used for joint modeling of the mineralization elements. Application of the MAF method to alr-transformed data yielded four factors, of which factor 4 was interpreted as the mineralization factor due to the strong correlation coefficients between all elements. The mineralization factor was modeled using the SGS method to evaluate the relations between rock types, alterations and the factor. Besides, rock types and alterations were modeled by the SIS method for 20 realizations in order to produce a probabilistic model for rock types and alterations. The interpretation was conducted via E-type maps of the rock types and alterations. The results demonstrate that higher concentrations of the mineralization elements are related to the existence probability of VBX and DP rock types rather than DLT and QDP. Furthermore, the frequent occurrences of sericitic alteration in the deposit did not allow us to have accurate interpretation about the relations between the mineralization factor and alterations. Finally, silicification may have a positive relation with the mineralization according to the mean value of the Au in silicified samples in drillhole data and the mean value of the mineralization factor in presented alteration models.

## Acknowledgements

The authors would like to thank CESCO mining company for providing the geochemical data and also Mr. Mokammad Reza Fani Pakdel for his kind cooperation.

- © 2017 The Author(s)