Stochastic Modeling of Chemical Compounds in a Limestone Deposit by Unlocking the Complexity in Bivariate Relationships

: Modeling multivariate variables with complexity in a cross-correlation structure is always applicable to mineral resource evaluation and exploration in multi-element deposits. However, the geostatistical algorithm for such modeling is usually challenging. In this respect, projection pursuit multivariate transform (PPMT), which can successfully handle the complexity of interest in bivariate relationships, may be particularly useful. This work presents an algorithm for combining projection pursuit multivariate transform (PPMT) with a conventional (co)-simulation technique where spatial dependency among variables can be deﬁned by a linear model of co-regionalization (LMC). This algorithm is examined by one real case study in a limestone deposit in the south of Kazakhstan, in which four chemical compounds (CaO, Al 2 O 3 , Fe 2 O 3 , and SiO 2 ) with complexity in bivariate relationships are analyzed and 100 realizations are produced for each variable. To show the e ﬀ ectiveness of the proposed algorithm, the outputs (realizations) are statistically examined and the results show that this methodology is legitimate for reproduction of original mean, variance, and complex cross-correlation among the variables and can be employed for further processes. Then, the applicability of the concept is demonstrated on a workﬂow to classify this limestone deposit as measured, indicated, or inferred based on Joint Ore Reserves Committee (JORC) code. The categorization is carried out based on two zone deﬁnitions, geological, and mining units.


Introduction
Accurate evaluation of viable mineral resources is fundamentally important in optimal sustainable development and mine planning procedures since it has an enormous impact on the value of produced metals and formation of technical plans, from extraction to closure of the mine [1][2][3]. In the world of emerging new technologies and increasing complexity for modeling the grades and chemical compounds in deposits, precise estimation and classification of resources is becoming more crucial, as it can vary depending on the technological modifications that change the notion of mine planning [4]. Nowadays, the "best practice" in resource classification takes into account assumptions made in grades and tonnages of orebodies, as well as in the methods applied, future expenses, and commodity costs. One of the main ways to solve these issues is to consider uncertainty and risk assessment [5][6][7]. Consequently, integrating enhanced geostatistical techniques in order to quantify uncertainty in the process of evaluating resources allow practitioners to effectively unlock the complexity in the modeling of mineral deposits [8,9].
In every mining project, the discrepancy between the estimated resource model and actual production value is the main concern, which leads to low probability of meeting production targets anamorphosis, which transforms variables into standard Gaussian form [28], or a quantile-quantile based approach [29]. With two or more variables, multivariable geostatistical modelling is based on the multivariate Gaussianity assumption [16,30]. In this technique, the variables should be independently transferred to normal Gaussian distribution. Then, (co)-simulation can be performed over normal score data, taking into account the cross-dependency functions that are defined by, for instance, a linear model of co-regionalization. The obtained realizations can subsequently be back-transformed to original scale in order to approximately reproduce the original linear multivariate relationships among the variables. Despite the fact that independent normal score transform ensures that each item of normal score data is multi-Gaussian, the multivariate Gaussianity assumption is violated when there is complexity in bivariate relations between pairs of variables. For example, with features such as geological constraints, heteroscedasticity, and nonlinearity, as illustrated in Figure 1 [31], that employ the typical (co)-simulation algorithm, using the multi-Gaussianity assumption is problematic. In other words, in a typical (co)-simulation paradigm, in order to implement the algorithms effectively, the normal score transformed data should follow the elliptical (Figure 1a). However, in the case of nonlinearity, heteroscedasticity, or linearity constraints (Figure 1), the multi-Gaussianity assumption cannot be respected and the outputs of typical (co)-simulation algorithms are unsatisfactory. To cope with this difficulty, one idea is to use another transformation technique such as projection pursuit multivariate transformation (PPMT) [23,24] rather than the regular normal score transformation (e.g., quantile-quantile or Gaussian anamorphosis [28,29]) that is common for univariate transformation.

Projection Pursuit Multivariate Transform Steps
The steps of projection pursuit multivariate transform are composed of preprocessing and projection pursuit. In the preprocessing steps, variables are transformed to normal score values and linear dependency is removed. In projection pursuit, complexity dependency is removed.

Preprocessing Steps
Normal Score Transformation Consider data matrix A, which has M variables and N observations such that A M×N is transformed through Equation (1) into normal score [29]: where E and C −1 are original data cumulative distribution function (CDF) standard and normal data CDF, respectively.

Data Sphering
One of the requirements of the projection pursuit algorithm is to center the data with an orthogonal covariance matrix and variance, which can be done by applying the last preprocessing step, data sphering, which uses Equation (2): where R and B are eigenvector and diagonal eigenvalue matrices, respectively, obtained from the A covariance matrix's spectral decomposition.

Projection Pursuit
After completing the preprocessing steps, the projection pursuit algorithm can be computed, taking into consideration that projection is Q = Aβ, where β is the unit length vector of h × 1 dimension related to projection Q. With multi-Gaussianity of A, every unit length vector β must yield a Q that is univariate Gaussian. In order to measure the univariate non-Gaussianity, the projection index T(β) test statistic is designated; when projection Q is finely Gaussian, T(β) is equal to zero. The projection pursuit algorithm uses optimized search, which is focused on identifying the β where the projection index T(β) is the highest. Finding the highest projection index will also find the maximum non-Gaussian projection Q. Consequently, by using the determined optimum unit length vector β, data A is transformed to A", in which the projection is standard Gaussian, as Q = A"β. In order to achieve this transformation, Equation (3) is used: where values of γ are unit vectors obtained by applying the Gram-Schmidt algorithm to compute the orthogonal matrix of Z [32]. The next step is the transformation, which can be achieved by multiplying A and Z by each other (Equation (4)): The next step considers the X transformation, which yields the projection of standard Gaussian Q while holding the same orthogonal matrix: The last step is multiplication of Z T by X(AZ), which results in back-transformation to the initial basis (Equation (6)): After back-transformation, optimized search can be used again to detect other complexities of A".

Stopping Criteria
To select the optimum value for projection index T(β), considering the optimum h dimensions and number of N observations is critical due to the fact that a high value of h dimensions results in trouble with identifying the complexity, and a small amount of N observations leads to low reliability in detection. To choose the target projection index, a special algorithm is applied [22].

Application
Overall, the implementation of projection pursuit multivariate transformation on multivariate data is based on the forward-and backward-transformation techniques. For the sake of simplicity, Figure 2 shows a schematic illustration to explain the implementation of the PPMT technique [22].

Proposed Algorithm
As already explained, the (co)-simulation approaches based on multi-Gaussianity assumptions (e.g., turning bands [33] and sequential Gaussian [34] co-stimulation) are not suitable for modelling purposes with such bivariate complexity. On the contrary, PPMT is able to model these characteristics based on factor transformation. In this technique, the dataset is first transferred to PPMT factors, in which the correlation among the factors become almost zero [22,24]. Since the cross-dependency among the factors is zero, the independent simulation can be applied over each factor and the outputs should be back-transformed to original scale. However, one of the main difficulties is related to the decorrelation step, where sometimes the forward transformation is not able to completely remove the inherent correlation between pairs of variables and a mild dependency may remain among the transformed factors that must not be negligible. One solution is to employ one more decorrelation steps, such as minimum/maximum autocorrelation factor (MAF) or principal component analysis (PCA), to make sure the correlation at lag 0 and other arbitrary lag is substantially removed. However, these methods do not guarantee that the correlation will be entirely eliminated through all other lags [35]. A way around this impediment is to co-simulate the PPMT factors even with small correlation that are left after forward transformation of the original variables to factors by inference of cross-dependency functions using the linear model of co-regionalization [25]. The output of the simulation is then back-transformed to original scale in the same way as in the independent simulation. In this paradigm, no MAF or PCA transformation is needed. The steps of the proposed algorithm are as follows:

1.
Exploratory data analysis of multivariate data 2.
Investigation of the level of complexity in bivariate relation analysis 3.
Examination of removing cross-correlations among variables by using cross-correlogram 5.
Inference of cross-dependency functions by linear model of co-regionalization (LMC) 6.
(Co)-simulation of PPMT transformed factors taking into account the fitted LMC 7.
PPMT backward transformation of simulated results (realizations) into original scale 8.
Validation of the output by statistical analysis tools In order to show the capability of the proposed model, a case study is presented and the outputs after validation are taken into account for mineral resource classification, a critical issue in international reporting of deposits.

Case Study: Aktas-South Deposit in Kazakhstan
This case study is the Aktas-South deposit, which is located in Reddipalayalam, in the state of Tamil Nadu, south of Kazakhstan. The deposit is limestone being used at a cement plant that belongs to the Aktas Group. For confidentiality reasons, the exact location of the deposit is not disclosed. The regional geology of this deposit falls in the Masanchi and Uzun formations in the group of Sabanbay rocks. The limestone in the deposit refers to Masanchi formation. The whole sedimentary formation of this zone is related to marine origin due to marine transgressions and regressions with corresponding fluctuations. As the limestone beds lack marl intrusions, they are composed of clusters of marine organisms. The limestone in the Aktas South area is divided into four groups according to their outward aspects, as shown in Table 1.

Exploratory Data Analysis in Limestone Deposit
The dataset consists of 4553 samples with homotopic sampling patterns. This configuration means that the data are available through all sample points and all variables have been assayed on equal sets of sample locations [36].
The dataset is composed of four chemical compounds: iron oxide (Fe 2 O 3 ), aluminum oxide (Al 2 O 3 ), silicon oxide (SiO 2 ), and calcium oxide (CaO), whose values are assayed in percentages. Similar to all geostatistical projects, the first step is exploratory data analysis to identify global statistical characteristics of the underlying variables. First, possible outliers and duplicated data were recognized. The presence of outliers in the dataset makes the inference of statistical parameters problematic and nonrepresentative [37,38]. These aberrant values intentionally influence the variance and result in sharp fluctuations in variogram analysis [39]. In addition, detection and removal of duplicate data is also important prior to any geostatistical analysis. One of the problems is that these repeated values generate singular matrices in kriging systems, leading to unestimated blocks surrounding the duplicate locations [40]. After removing duplicate locations in this study, the variables are plotted in a probability plot. This statistical tool helps to detect and fix extreme and innermost values [41]. Following this, some outlier values are detected for Fe 2 O 3 and CaO ( Figure 3). In the distribution of Fe 2 O 3 , values more than about 18% are considered as extreme values, and in the distribution of CaO, values less than about 19% are recognized as innermost values, while the other two distributions (SiO 2 and Al 2 O 3 ) sound reasonable given that no significant outliers were identified. Once these outlier values are detected, the capping approach [41] is applied to treat them accordingly in order to preserve them in the dataset after examining whether they are valid samples and not erroneous. In this technique, the values in the upper and lower tails of the distribution should be moved back to the previous maximum value and forward to the next minimum value, respectively. The next step in statistical analysis is related to the declustering process. This step was not necessary, as the sampling pattern is almost regular in the region ( Figure 4). The most critical chemical compounds for this deposit are related to the variability of CaO and SiO 2 . These two variables introduce high-quality limestone for high amounts of CaO (>10%) and low amounts of SiO 2 (<40%). In addition, Fe 2 O 3 plays an important role, yet does not have as significant an influence as SiO 2 in favor of CaO. As can be seen, the majority of drillholes illustrate a high concentration of CaO, which is distributed homogeneously in the region. Interestingly, SiO 2 reveals poor concentration in the same locations, meaning it satisfies the quality of limestone for the entire deposit, corroborating high CaO and low SiO 2 . This visual inspection also indicates a strong spatial dependency between CaO and SiO 2 , in which there is a negative impact on the local distribution of these two chemical compounds. This shows a high concentration of CaO versus a low concentration of SiO 2 , which is important for this deposit. This phenomenon motivates a further investigation into the cross-correlation structures among these chemical variables toward better decision making for the selection of efficient geostatistical algorithms for 3D modelling and mineral resource evaluation. Then, statistical parameters of the data were computed, as shown in Table 2. The coefficient of variation (COV) for all variables, particularly CaO and SiO 2 , is less than 2.0, which indicates that the distribution of data has no significant harsh variability and predictive models can be suitable and meaningful [41]. In this limestone deposit, as previously mentioned, CaO and SiO 2 are two critical variables for ore/waste selection. For this purpose, areas with more than 10% CaO and less than 40% SiO 2 define ore zones, and areas with less than 10% CaO and more than 40% SiO 2 introduce waste zones based on mining excavation destination. Before initiating the modelling process, it might be of interest to calculate two global recovery functions, fraction of recoverable ore above or below the cutoff and mean grade above or below the cutoff [28,42]. These two parameters are calculated by bivariate cumulative distribution functions computed over CaO and SiO 2 as follows: Fraction of total tonnage in the specified cutoff for ore: Mean grade in the specified cutoff for CaO in ore: These two parameters show that about 66.2% and 33.8% of the entire deposit may be ore and waste, respectively, which is an interesting economical characteristic of this mine deposit. In order to investigate the cross-dependency among these four chemical compounds, the global correlation coefficient was calculated, as presented in Table 3. Fairly good negative and positive correlations can be seen between Fe 2 O 3 and CaO (−0.64) and between Fe 2 O 3 and SiO 2 (approximately 0.53). The highest dependency is seen between CaO and SiO 2 , which has negative correlation with almost −0.94 correlation coefficient. This corroborates the visual interpretation already provided in the location map of sampling points ( Figure 4). Other correlations, such as those between Fe 2 O 3 and Al 2 O 3 , Al 2 O 3 and CaO, and Al 2 O 3 and SiO 2 , are somewhat low. These values only give a general perspective on the linear dependency that exists among the variables and may not be suitable for examining whether or not complex characteristics such as nonlinearity, heteroscedasticity, and geological constraints may exist. In order to examine the latter characteristics, the bivariate relation in Figure 5 is presented as scatter plot between pairs of the variables Al 2 O 3 , CaO, Fe 2 O 3 , and SiO 2 . This statistical diagram is suitable to explore relationships such as complexities and linearity features between pairs of variables. This is an interesting illustration of different complexities between co-variables, starting from heteroscedastic characteristics between Fe 2 O 3 and Al 2 O 3 , and CaO and SiO 2 ; nonlinearity observed between SiO 2 and CaO; and geologic constraints between CaO and Al 2 O 3 .

PPMT Forward Transformation
As mentioned before, one of the objectives of this paper is to jointly model variables by considering original correlations among variables to reconstruct the shape of complex bivariate relations. It is explained that correlation shapes between pairs of variables are complex ( Figure 5). These types of features motivate the use of factorization approaches such as projection pursuit multivariate transform (PPMT) to jointly model the underlying four variables (Al 2 O 3 , CaO, Fe 2 O 3 , and SiO 2 ). The result of this modelling approach is then applied for mineral resource classification. However, prior to any geostatistical modelling, whether it is independent simulation or co-simulation after forward PPMT transformation, the removal of correlations after this first transformation should be assessed. Therefore, in this study, as a common practice in forward PPMT transformation, the variables are first transformed to PPMT factors and then undergo one further transformation by MAF, implemented to completely remove cross-correlation. This can be evaluated by a cross-correlogram ( Figure 6).
The results of the cross-correlogram between MAF factors show that a small amount of correlation is still resistant in some lags. For instance, one can recognize that the correlation in the first lag (~100 m) is around 25% between SiO 2 and Al 2 O 3 even after this MAF transformation ( Figure 6). This signifies that MAF is not able to decorrelate the variables over all the lag separations. In this respect, it is not advocated to use the independent simulation due to the remaining small correlations among factors. In order to cope with the proposed algorithm in this case, it was decided to employ co-simulation right after the initial PPMT forward transformation of underlying variables irrespective of considering any further MAF transformation. For this, once again, the cross-correlation among the transformed PPMT factors is examined through the cross-correlogram and before the MAF step. The results show that the correlation manifests itself through some lags (e.g.,~12% between SiO 2 and CaO at a lag separation of 150 m and~11% between SiO 2 and Al 2 O 3 at a lag separation of 400 m; Figure 7), although at a lag separation of 0, the correlation is significantly removed (Figure 8). Even these small amounts of correlation among the PPMT factors are important and provoke applying the co-simulation algorithms via PPMT transformed factors and considering the linear model of co-regionalization.

Variogram Inference
As mentioned in a previous section, the transformed factors retain the correlation at other lag distances except zero, and a further decorrelation transformation technique such as MAF cannot completely remove the inherent correlations. Following the proposed algorithm in this research, the next step is variogram analysis over the transformed variables. This step is needed to establish the co-kriging system in co-simulation algorithms, taking into account the linear model of co-regionalization. This latter introduces even small correlations in co-simulation algorithms. In this respect, inference of direct and cross-variograms for all four PPMT transformed variables is implemented. The sill of experimental cross-variogram as a measure of joint variability to some extent reflects the magnitude of the correlation between the variables [43]. In the case of standardized variables (variance equal to 1), the sill of cross-variograms must be the cross-correlation between collocated values of those variables [44].
It is also worth mentioning that anisotropy was examined through the original dataset (not transformed), and it was seen that the existence of anisotropy in the region was improbable. For fitting of theoretical variograms over experimental variograms, the linear model of co-regionalization [45] with semi-automatic technique is chosen, in which the semi-definiteness condition [36,46] is respected through the process of fitting. In this model, direct and cross-covariances are defined as the sum of basic covariances. Therefore, by this technique and after computation of experimental variograms, two nested structures (spherical, with the first scale factor as 54 m and the second scale factor as 216 m) are fitted accordingly, with no nugget effect, as can be seen in Figure 9 (Equation (7)).

Stochastic Modeling in Limestone Deposit
After variogram inference over PPMT transformed covariates and following the proposed algorithm in this study, the step for implementing co-simulation takes into account grid nodes with dimensions of 48 m × 48 m × 11 m for each block to jointly model all transformed factors. Turning bands co-simulation (TBCOSIM) is selected because of its versatility and reliability in reproducing global and statistical parameters when compared to other Gaussian co-simulation algorithms [47]. In this method, each variable is first simulated nonconditionally and then through a post-processing step by co-kriging conditioned on the available information and borehole records [33]. In this regard, ordinary co-kriging is used with moving neighborhood ranges equal to the range of variograms (200 m) attending up to 50 nearest surrounding sample points in the process of conditioning. Multiple-search strategy is also selected for this purpose, since it shows better results compared to single-search strategy [48]. One of the main criticisms against TBCOSIM, however, is related to producing artifacts because of turning lines. To deal with this difficulty, the number of lines should be reasonably increased [33,49]. Therefore, 1000 lines was set for the turning bands simulation method to eliminate possible stripping effects. The number of realizations was considered to be 100. The realizations that interpret the spatial variability of each factor were then back-transferred to original scale though backward PPMT transformation. Next, E-type maps were generated by averaging through 100 realizations within each block in co-simulated element over the variability in original scale. E-type maps of CaO, Al 2 O 3 , Fe 2 O 3 , and SiO 2 were produced through 100 realizations, as shown in Figure 10. Before further analysis over the simulated results, it is necessary to check whether the outputs of this proposed algorithm are statistically valid.

Validation
Validation analysis in this section is concerned with the reproduction of original statistical characteristics such as mean, variance, correlation coefficient, and shape of bivariate relation. This type of validation process is required to give practitioners insight with regard to the level of reliability of the aforementioned approach in the section on mineral resource estimation and classification. A comparison between the means of original variables (Fe 2 O 3 , Al 2 O 3 , SiO 2 , and CaO) and the means of simulated and back-transformed variables calculated from PPMT through 100 realizations is shown in Figure 11. As can be seen, PPMT is able to reproduce the original mean values of each variable.  (Figure 12), it is intuitively determined that PPMT produces satisfactory outputs in terms of reproduction of variance. However, minor deviations, such as for Fe 2 O 3 , as can be seen from this figure, are referred to the influence of conditioning data [24,50,51]. However, this tiny departure of average of variances for 100 realizations from original variance is not remarkably significant.
Finally, yet importantly, comparisons of correlation coefficients provided by PPMT and original correlation coefficients are examined among all variables (Fe 2 O 3 , Al 2 O 3 , SiO 2 , and CaO) through 100 realizations. As it can be seen from Figure 13, co-simulation methodology shows satisfactory results in the reproduction of original correlation coefficients among co-variables. This good performance among variables can be explained by the fact that co-simulation considers the intrinsic correlation between variables by incorporating the linear model of co-regionalization [51][52][53].   Another part of the comparison is related to examining the ability of the proposed method to reconstitute the shape of the original bivariate relations between pairs of chemical compounds. By visual inspection of reproductions of the shape of correlation for the underlying pairs of variables, one can see the difference between original and reproduced shapes in Figure 14. It should be mentioned that only one realization was taken to illustrate the reproduction of correlation shape, and it was randomly chosen as a first realization for all cases for the sake of fairness. As can be seen from the scatter plots of Fe 2 O 3 and Al 2 O 3 , PPMT co-simulation successfully reproduces the shape of original correlation. The same features are evidently demonstrated in the reproduction of the shape of original correlation between Fe 2 O 3 and CaO, Fe 2 O 3 and SiO 2 , and SiO 2 and CaO. However, inadequate results in the reproduction of the shape of correlation between SiO 2 and Al 2 O 3 , and CaO and Al 2 O 3 can be seen. Overall, the proposed approach based on the combination of PPMT and co-simulation is capable of reconstructing the original shape of correlation.

Mineral Resource Classification
Based on the JORC code definition (www.JORC.org), mineral resources can be classified into measured, indicated, or inferred based on the level of confidence. In fact, the guideline in this code was developed to assure transparency for investors in the declaration of mineral resources in order to prevent fraud [41]. A measured mineral resource is part of a deposit that presents a high level of confidence in the estimation of recovery functions such as tonnage, mean grade, and metal quantity above cutoff. Indicated mineral resources can be assigned to those areas that demonstrate a reasonable level of confidence in the estimation of recovery functions such as tonnage, mean grade, and metal quantity above the cutoff. Inferred mineral resources represent locations for which the recovery functions are evaluated at a low level of confidence. In this regard, there are mainly two algorithms for mineral resource classification, connected with either deterministic or stochastic paradigms. In the latter, one needs to employ, for instance, geostatistical simulation techniques to produce different scenarios of a mine deposit (i.e., realizations) with the aim of quantifying the uncertainty. This leads to probabilistic reporting of the mineral resource measures.
Mineral resource classification in this limestone deposit is done on the basis of two main targets, following the stochastic approach. One is related to classification based on ore zone definition (geological), in which mostly CaO and Fe 2 O 3 are two critical components for the underlying zone. The second is mainly concerned with mining excavation units where CaO, SiO 2 , and Al 2 O 3 are vital variables for this zone definition. In the following, we present the results for both types of mineral resource classification that define the different categories based on ore zone definitions and mining units.

Ore Zone Definitions
One of the main objectives of this study is to employ the proposed algorithm based on a combination of PPMT transformation and co-simulation algorithm, and to show the results of mineral resource classification by demonstrating their distinctions. As shown previously in the section on validating simulation results with regard to reproducing original mean, variance, correlation coefficient, and shape of bivariate relations between pairs of variables, the PPMT and co-simulation methodology corroborates that the outputs of realizations are statistically sound and can be applied for further processing, such as mineral resource classification. Aktas-South deposit is grouped into three main ore zones based on the definitions of chemical cutoffs as indicated in Table 4. CaO in this limestone deposit is the main chemical compound. Resource estimation of this variable based on JORC code is shown in Table 5 for each ore zone. In order to show the distinctions between outputs of resource classification of CaO (measured, indicated, inferred, and total tonnage), the results are constructed in one unique graph, shown in Figure 15. As previously explained, the classification of resources is modelled according to uncertainty, in which measured resources are quantified 90% of the time within ±15%, indicated resources within ±15% and ±30%, and inferred resources within ±30% and ±100%, while other materials within more than ±100% are not distinguished as resources [41].  As can be observed from Figure 15, marl-chert has the highest total tonnage of CaO, 208 Mt, according to the PPMT co-simulation method (proposed algorithm), followed by pale yellow limestone and brown limestone, with 181 Mt and 92 Mt, respectively. In the measured category, marl-chert ore and brown limestone ore tonnage are almost the same, 53 Mt, while pale yellow limestone has higher tonnage, 102 Mt. In the indicated and inferred categories, marl-chert has higher tonnages, 30 Mt and 125 Mt, respectively, followed by pale yellow limestone, 11 and 68 Mt, respectively, and brown limestone has the lowest tonnages, 5.5 Mt and 34 Mt, respectively. It should be mentioned that ore zone definition in this step has two restrictions (chemical cut-offs), and because of that results of resource classification differ.

Mining Units
In this deposit, there is another classification based on chemical cutoff definitions, but this time from the mining destination perspective, whereas the previous tonnage classification was based on the definitions of ore zones in favor of geological interpretations. The definitions of these mining units are presented in Table 6. In this paradigm, the units are mainly established in accordance with three main chemical compounds, CaO, SiO 2 , and Fe 2 O 3 . As also shown in Table 3, there is strong correlation among these three variables. The results of the proposed approach, in fact, show that the correlation coefficients are reproduced at a satisfactory level of confidence. The importance of this is shown in this section; for instance, these three variables are the key factors in favor of grouping the Aktas-South limestone deposit for mining excavation, which consequently impacts the rigorous classification of mineral resources. The list of mining units based on chemical cutoffs is shown in Table 6. It should be mentioned that in this classification, the number of restrictions (chemical cutoffs) reaches four, which means that the results of classification will be according to seven restrictions. As can be seen, CaO and SiO 2 are two chemical compounds that define the ore/waste classification technique. The tonnage for each classification based on JORC code in this deposit is summarized in Table 7, and all results of the mentioned method with resource classification are summarized in single graphs for highly green limestone (HGLS), brown limestone (BROWNLS), ferrous limestone (FEROLS), cherty limestone (CHERTYLS), cherty limestone 2 (CHERTYLS2), MARL, and WASTE in Figure 16, following the same method as explained in the previous section [41].   Depending on the chemical restrictions of each lithology, resource tonnages have different trends of results by category. The highest total tonnages are shown in HGLS and MARL, followed by almost the same tonnages in BROWNS, FEROLS, and CHERTYLS (see Figure 16). The lowest total tonnages are computed in CHERTYLS2 and WASTE. Turning to the measured category, the highest tonnage is seen in HGLS, followed by MARL, BROWNS, FEROLS, and CHERTYLS, and the lowest tonnages are reported in CHERTYLS and WASTE. With an increased number of geological restrictions, the results of resource classification change more chaotically. However, it is clearly shown that CHERTYLS2 and WASTE, which have narrow restrictions in CaO, have the lowest tonnages in every category.

Conclusions
This contribution provides a technique for the multivariate geostatistical analysis of co-regionalized variables, particularly for datasets where complexity exists between bivariate relationships. The proposed algorithm in this study is based on the combination of a conventional stochastic simulation approach (Gaussian (co)-simulation) and a recently developed factorization technique (projection pursuit multivariate transformation (PPMT)). A main difference from other multivariate geostatistical approaches is that in PPMT, there is no need to further transform the factors to remove the correlation. Instead, a linear model of co-regionalization should be defined to introduce the small cross-spatial dependency among the transformed variables, in which it can be inferred from the factors right after PPMT forward transformation.
In addition, the real case study, limestone deposit, illustrates the effectiveness of the proposed algorithm for mineral resource classification based on the two-zone separation following the JORC code definition. To do so, the outputs of the realizations are first examined to determine whether they are statistically valid, then taken into account to classify resources as measured, inferred, or indicated. Through the validation step, reproduction of mean, variance, and global (cross)-correlation is examined, illustrating that all simulation results, on average, converge to the original statistical parameter, indicating the robustness of the proposed algorithm. Behind the seeming simplicity and great number of practical investigations, much effort is still needed to deal with more complex sampling configuration, particularly in the presence of partially or heterotopic sampling patterns. In such cases, using PPMT is restricted and an imputation technique may be an alternative to combine PPMT with (co)-simulation.