School of Mining and Geosciences MODELING COMPLEX RELATIONSHIPS IN GEOMETALLURGICAL VARIABLES: ENHANCING METHODS WITH ACCEPTANCE-REJECTION AND HIERARCHICAL GAUSSIAN CO-SIMULATION Thesis supervisor: Nasser Madani Student: Shingiskhan Kuanyshev Thesis submitted to the School of Mining and Geosciences of Nazarbayev University in Partial Fulfillment of the Requirements for the Degree of Master of Science in Mining Engineering Submission date: 08/04/2024 Astana, 2024 ii Declaration This is to certify that to the best of my knowledge; the content of this thesis is my own work. This thesis has not been submitted for any degree or other purposes. I certify that the intellectual content of this thesis is the product of my own work and that all the assistance received in preparing this thesis and sources have been acknowledged. Shingiskhan Kuanyshev iii Abstract Resource estimation is the basis of efficient and sustainable mining. Accurate mineral resource estimation is critical to optimizing mine planning, minimizing waste, and ensuring the economic viability of mining operations. Traditionally, resource estimation focused primarily on grade, the concentration of valuable minerals in ore bodies. However, the evolving complexity of modern mining requires a holistic approach that includes not only the grade but also the geometallurgical properties of the ore. Geometallurgical properties, which include attributes such as mineralogical composition, texture, and work index, play a key role in shaping mining operations. Understanding and modeling these properties is essential to unlocking the full potential of mineral deposits. In the context of mining, geological complexity often leads to complex non-linear bivariate relationships between different attributes of ore bodies. These relationships can pose significant challenges for resource modeling, as traditional methods such as Principal Component Analysis (PCA) and Minimum/Maximum Autocorrelation Factor (MAF) are ill-suited to handle such complexities. These methods are inherently linear and may fail to capture the nuanced interactions and dependencies within geologic datasets. This research paper presents a new approach to address the complexity of resource estimation in mining, particularly when dealing with non-linear bivariate relationships between geometallurgical properties. The proposed method combines the accept-reject method with hierarchical sequential Gaussian co-simulation. This approach enables the careful modeling of complex relationships (non-linearity) within geological data, leading to more accurate resource estimates and better-informed mining decisions. A case study is presented in which the acceptance-rejection method with hierarchical sequential Gaussian co-simulation is applied to the modeling of two geometallurgical properties, recovery and chalcopyrite. The study shows how this innovative approach increases resource estimation accuracy by capturing non-linear dependencies and spatial variability between these two variables. Due to the hierarchical nature of geological data, the method adapts to different scales of variability, resulting in more realistic and practical resource models. The findings of this research not only highlight the importance of integrating geometallurgical properties into resource estimation, but also provide a valuable solution for solving complex nonlinear bivariate relationships in geostatistical analysis of other regionalized variables. This approach has the potential to revolutionize resource modeling practices in the mining industry, leading to more sustainable, efficient, and economically viable mining iv operations. As mining continues to face evolving challenges and requirements, it is essential to incorporate advanced techniques such as the accept-reject method with hierarchical sequential Gaussian co-simulation to exploit the full potential of the Earth's mineral resources. Keywords: Acceptance-Rejection, Sequential Gaussian co-simulation, geometallurgy, non-linear bivariate relationship. v Acknowledgments My heartfelt thanks go to my supervisor, Professor Nasser Madani, for his unwavering support, invaluable advice, and professional guidance throughout my thesis journey. Professor Madani's assistance in defining the thesis topic has been crucial; his profound knowledge of the mining industry, especially in resource estimation and geostatistics, has significantly shaped the development of my research. His extensive experience and expertise have been essential in the completion of my work. I am also deeply grateful for the opportunity to serve as a research assistant under Professor Madani, whose professionalism in the mining field I greatly admire and respect. Furthermore, I wish to acknowledge Nazarbayev University and the School of Mining and Geosciences for providing me with the chance to engage with highly experienced professors who are always ready to offer their assistance and guidance. The research-oriented environment of NU and its supportive academic community have been pivotal in the successful completion of my thesis. vi Table of Contents Declaration ....................................................................................................................................................................ii Abstract ....................................................................................................................................................................... iii Acknowledgments ......................................................................................................................................................... v List of Figures ........................................................................................................................................................... viii List of Tables ................................................................................................................................................................ ix 1 INTRODUCTION ...................................................................................................................................................... 1 1.1 Background ............................................................................................................................................................. 1 1.2 Problem Statement .................................................................................................................................................. 3 1.3 Objective of the Research ........................................................................................................................................ 4 1.4 Hypothesis ............................................................................................................................................................... 5 1.5 Significance to the Industry ..................................................................................................................................... 5 1.6 Scope of Work ......................................................................................................................................................... 5 1.7 Thesis Timeline ....................................................................................................................................................... 6 2 LITERATURE REVIEW ........................................................................................................................................... 7 2.1 Kriging .................................................................................................................................................................... 7 2.2 Independent Simulation ........................................................................................................................................... 8 2.3 Co-Simulation ......................................................................................................................................................... 8 2.4 PCA ....................................................................................................................................................................... 10 2.5 MAF ...................................................................................................................................................................... 10 2.6 PPMT .................................................................................................................................................................... 11 2.7 Flow Anamorphosis .............................................................................................................................................. 12 3 METHODOLOGY ................................................................................................................................................... 14 3.1 Hierarchical Sequential Gaussian Co-Simulation ................................................................................................. 14 3.2 Acceptance-Rejection method ............................................................................................................................... 16 3.3 Proposed Methodology .......................................................................................................................................... 17 3.4 Data Requirements ................................................................................................................................................ 19 3.5 Expected Results ................................................................................................................................................... 20 4 RESULTS ................................................................................................................................................................. 22 4.1 Case study ............................................................................................................................................................. 22 4.2 Exploratory data analysis ...................................................................................................................................... 24 4.3 Variogram analysis ................................................................................................................................................ 26 4.4 The outcomes of Conventional and Proposed methods ......................................................................................... 30 vii 4.5 Geostatistical validation of both methods ............................................................................................................. 37 4.5.1 Reproduction of global statistical parameters ............................................................................................ 37 4.5.2 Reproduction of bivariate relationship ....................................................................................................... 40 5 CONCLUSIONS ...................................................................................................................................................... 43 REFERENCE LIST ..................................................................................................................................................... 45 viii List of Figures Figure 1. Gantt Chart of the work.................................................................................................................................. 6 Figure 2. Geostatistical Simulation via FA (Talebi et al., 2018) ................................................................................. 13 Figure 3. Example of the Scatter Plot for Acceptance-Rejection method (“Monotonic relationships,” n.d.) ............. 18 Figure 4. Location Map of Chalcocite ......................................................................................................................... 23 Figure 5. Location Map of Recovery ........................................................................................................................... 24 Figure 6. Scatter plot between Recovery and Chalcocite ............................................................................................ 26 Figure 7. Scatter Plot of Chalcocite and Recovery in Normal Scores with their Marginal Distribution ..................... 27 Figure 8. Lagged Scatter Plot of Normal Score of Recovery ...................................................................................... 28 Figure 9. Lagged Scatter Plot of Normal Score of Chalcocite .................................................................................... 28 Figure 10. Direct and cross-variograms of the Recovery and Chalcocite in normal scores in both the horizontal (black) and vertical (blue) directions; points: experimental variogram, and solid line: variogram model (derived linear model of coregionalization) ............................................................................................................................... 30 Figure 11. Random realization of Chalcocite and Recovery produced by the proposed algorithm based on hierarchical sequential Gaussian simulation using acceptance-rejection method ........................................................ 31 Figure 12. Random realization of Chalcocite and Recovery produced by the conventional hierarchical sequential Gaussian simulation algorithm .................................................................................................................................... 32 Figure 13. E-type realization of Proposed Approach of Recovery based on hierarchical sequential Gaussian simulation using acceptance-rejection method ............................................................................................................ 33 Figure 14. E-type realization of Conventional Approach of Recovery based on conventional hierarchical sequential Gaussian simulation ................................................................................................................................................... 34 Figure 15. Variance Map of Chalcocite obtained from conventional hierarchical Gaussian co-simulation ................ 35 Figure 16. Variance Map of Recovery obtained from conventional hierarchical Gaussian co-simulation ................. 35 Figure 17. Variance Map of Chalcocite obtained from proposed hierarchical Gaussian co-simulation using acceptance-rejection method ...................................................................................................................................... 36 Figure 18. Variance Map of Recovery obtained from proposed hierarchical Gaussian co-simulation using acceptance-rejection method ....................................................................................................................................... 36 Figure 19. Histogram of non-linear correlation coefficients between Chalcocite and Recovery from 100 realizations produced by conventional and proposed co-simulation algorithms; solid line: mean of correlation coefficient obtained from 100 realizations, dashed line: original correlation coefficient .............................................................. 38 Figure 20. Comparison of Chalcocite histogram reproduction by conventional and proposed methods for one realization .................................................................................................................................................................... 39 Figure 21. Comparison of Recovery histogram reproduction by conventional and proposed methods for one realization .................................................................................................................................................................... 40 Figure 22. Reproduction of the bivariate relationship between Chalcocite and Recovery with Conventional method ..................................................................................................................................................................................... 41 Figure 23. Reproduction of the bivariate relationship between Chalcocite and Recovery with Proposed method ..... 42 ix List of Tables Table 1. Linear Correlation between Geometallurgical parameters ............................................................................ 20 Table 2. Non-Linear Correlation between Geometallurgical parameters .................................................................... 21 Table 3. Statistical parameters of the dataset before and after declustering ................................................................ 25 Table 4. Reproduction of the average statistical parameters over 100 realizations ..................................................... 37 1 1 INTRODUCTION 1.1 Background Resource modeling plays a pivotal role in the field of mining engineering as it serves as the fundamental basis for efficient mine design and the promotion of sustainable resource utilization (Gu et al., 2020). This resource estimation process offers mining engineers and geologists a comprehensive understanding of the spatial distribution, dimensions, morphology, composition, and grade in an ore deposit. Resource modeling serves as the fundamental basis for mine planning, playing a crucial role in shaping decisions that have a significant impact on the overall success of the mining operation. Historically, the primary emphasis in mine planning has been placed on the quality of ore resources. However, the modern perspective acknowledges that this limited perspective is inadequate. In contemporary mining operations, it is imperative to take into account not solely the ore grade, but also the geometallurgical characteristics of the deposit (Dominy et al., 2018). The change in viewpoint is motivated by the recognition that effective mining extends beyond the mere extraction of rich minerals. It encompasses a comprehensive understanding of the complicated properties of the ore, with the aim of optimizing operational procedures, minimizing waste generation, and guaranteeing sustained profitability in the long run. Throughout history, the central emphasis of resource modeling and mine planning has predominantly revolved around the assessment of ore grade. The ore grade is a measure of the concentration of precious minerals present in an ore deposit, usually expressed as the amount of metal contained per unit mass, such as grams per ton, or percentage. High-grade ores are considered economically advantageous due to their higher concentration of important minerals, which results in reduced processing and extraction efforts needed to obtain a specific quantity of metal. The appeal of high-grade ores is readily apparent, as they provide the potential for greater immediate financial gains and shorter times of recouping initial investments. Mining operations that focus on extracting high-grade minerals are sometimes regarded as more financially lucrative endeavors. Nevertheless, the exclusive focus on grades overlooks crucial factors related to the variety of ore bodies and does not consider the intricate nature of contemporary mining practices (Deutsch et al., 2015). 2 The significance of ore grade in the mining industry cannot be overstated; nonetheless, it is important to acknowledge that mine planning looking only at the grade presents a number of constraints, according to Deutsch et al. (2015), such as: ● Disregarding the variability of ore: It is uncommon for ore deposits to be homogeneous. ● Differences in mineral composition, texture, and geological structure can result in notable discrepancies in metallurgical characteristics. Failure to consider these variables can lead to suboptimal mining operations and reduced resource extraction efficiency. ● Waste management challenges arise from the implementation of a grade-centric approach, which frequently results in the generation of excessive waste material. This technique commonly entails the extensive extraction of non-valuable rock in conjunction with the desired ore. This phenomenon not only leads to an escalation in operating expenditures but also gives rise to environmental apprehensions pertaining to the management of trash and the restoration of land. ● Process optimization is a crucial aspect of contemporary mining operations, encompassing a range of intricate procedures such as crushing, grinding, flotation, and smelting. Each of these processes exhibits sensitivity to the properties of the ore. Achieving optimal resource recovery and energy efficiency in these processes is a tough task without a thorough comprehension of ore variability. ● The mining of high-grade ores has the potential to deplete resources at an accelerated rate, hence resulting in enduring environmental deterioration. The implementation of sustainable mining techniques necessitates the adoption of a comprehensive approach that takes into account the entirety of the ore deposit and strives to reduce the development of waste. Therefore, in light of the constraints associated with grade-centric methodologies, the mining sector has progressively embraced geometallurgy as a strategy to improve resource modeling and project planning. Geometallurgy is an interdisciplinary methodology that integrates geological, mineralogical, metallurgical, and geostatistical information to yield a thorough comprehension of the behavior of ore bodies (Hunt et al., 2019). Geometallurgical modeling considers many ore body features, reported by Rincon et al. (2019), including: 3 ● Mineral Composition: Geometallurgy investigates ore body mineral kinds and amounts. Processing methods vary by mineral, and their interplay might affect metallurgical behavior. ● It considers particle size, shape, and mineral relationships of the ore. These parameters affect crushing, grinding, and separating. ● Geometallurgy plots ore properties spatially, highlighting high and low variability. This data helps mine planners reduce waste and maximize resource recovery. ● Comminution Properties: It assesses the ore's response to crushing and grinding, crucial mineral processing procedures. Comminution qualities aid processing circuit design. ● Geometallurgy measures mineral liberation, which indicates how easily valuable minerals may be removed from gangue. This informs process design and recovery forecasts. ● Geometallurgy uses metallurgical test data to mimic processing conditions. Tests reveal the ore's behavior under different settings, enabling process improvement. ● Geostatistics: Modeling ore property spatial distribution improves resource assessment and mine design. 1.2 Problem Statement The intricate attributes of ore bodies frequently present significant difficulties in the realm of material modeling and mining technologies. One of the challenges that can occur arises from the complex nonlinear relationships that may exist among the many components of the deposit. Conventional approaches, such as employing Gaussian simulations or employing factorization techniques like principal component analysis (PCA) and minimum/maximum autocorrelation factors (MAF), may prove inadequate in achieving precise outcomes when attempting to construct simulations or models aimed at comprehending and forecasting these complicated relationships (Battalgazy and Madani, 2019). The geographic distribution and variability of orebody parameters, including grade, mineral composition, and texture, are crucial considerations for resource experts in the mining business. These factors, sometimes, exhibit robust nonlinear bivariate correlations that necessitate careful analysis and understanding. The interdependence of these traits is rare, and it is common for them to have a non-linear correlation. For instance, the correlation between the concentration of one metal and another metal in a mineral body does not exhibit a linear relationship, indicating that changes in one metal's concentration do not directly correspond to proportional changes in the concentration of the other metal (Malliaris et al., 2005). Conversely, the dynamics of their 4 relationship may be subject to the influence of geological complexity and interactions, leading to non-linear patterns of behavior. There exist various challenges that pose a threat to the conventional approaches employed in the field of resource modeling. ● Traditional approaches, such as employing conventional sequential Gaussian co- simulations, have demonstrated their effectiveness in accurately representing linear correlations and relatively uncomplicated geological formations. Nevertheless, when confronted with complex nonlinear bivariate relationships, these methodologies might have difficulties in accurately representing the inherent characteristics of the deposit structure. The Gaussian simulations in this study are predicated on the assumptions of normality and stationarity. However, it is important to acknowledge that these assumptions may not be valid in situations involving highly variable and nonlinear ore bodies. ● The factorization methods, namely Principal Component Analysis (PCA) and Minimum/maximum Autocorrelation Factor (MAF), are well recognized and utilized techniques in the field of data analysis and feature extraction due to their effectiveness in reducing dimensionality. The utilization of these tools facilitates the identification of patterns and underlying correlations within multivariate data sets. Nevertheless, it should be noted that these models are not naturally suited to effectively handle nonlinear relationships. When utilized in the context of complicated nonlinear relationships, these techniques might not sufficiently capture the fundamental framework, hence potentially resulting in misinterpretations and poor sampling outcomes. 1.3 Objective of the Research The main objective of this study is to thoroughly assess and enhance the use of the acceptance-rejection approach coupled with hierarchical sequential Gaussian co-simulation for modeling complex and non-linear bivariate associations in synthetic geometallurgical drill hole datasets. The primary aims of this study involve doing a comprehensive evaluation of the performance of the acceptance-rejection approach in comparison to conventional factor-based methodologies. Additionally, the study aims to employ hierarchical sequential Gaussian co- simulation to conduct a detailed analysis of the geometallurgical variables of interest. The study will also investigate how these methods can be synergistically combined to enhance modeling accuracy. A key part of the research involves validating these models against the synthetic data, ensuring their relevance and utility in simulated geological and mining 5 scenarios. Also, this research aims to check the application of these advanced statistical methods in the realm of geometallurgy, particularly in the context of synthetic datasets, to aid in more precise and efficient resource estimation and analysis in simulated mining environments. 1.4 Hypothesis In situations where conventional simulation and factorization techniques are insufficient, different approaches are required to effectively reflect the complex structure of nonlinear bivariate relationships. One particular methodology is the acceptance- rejection method with further implementation of Hierarchical Sequential Gaussian co-Simulation. In this work, synthetic geometallurgical dataset will be used, which was created by the logic of the work of Garrido et al. (2020). Therefore, based on the theory, it can be possible to predict that the proposed approach of acceptance-rejection method with further hierarchical sequential Gaussian co-simulation will provide good results in modeling such complex nonlinear relationships between geometallurgical properties. 1.5 Significance to the Industry The thesis research, which focuses on the application of modern statistical modeling techniques to analyze geometallurgical variables, presents significant optimism in terms of its potential to revolutionize the mining sector. This research aims to improve resource evaluation and ore body characterization by considering a comprehensive overview of geometallurgical parameters, beyond the exclusive consideration of ore grade. The use of this comprehensive strategy has the potential to result in improved operational efficiencies, decreased costs, and increased environmental sustainability. Furthermore, the use of customized processing procedures derived from these models has the potential to raise the rates of recovery and metallurgical performance, hence boosting the economic value and financial attractiveness of mining ventures and result in better mine planning and reduce the risks. 1.6 Scope of Work In order to comprehend current theories and pinpoint knowledge gaps, this research project will involve a thorough literature evaluation. Next, exploratory data analysis (EDA) is being used to look for patterns and anomalies in the dataset. This is accomplished through a simulation process in which relevant scenarios are modeled and simulated to produce data for analysis. Subsequently, a thorough analysis of the data collected will be conducted through the outputs of 6 the simulation, interpreting the results in light of this research topic. Ultimately, I will compare and contrast the results with existing literature to assess consistencies and discrepancies. 1.7 Thesis Timeline Figure 1. Gantt Chart of the work It is important to mention, the preliminary due date of submission of Research Thesis II to the reviewers is the week of 22nd-26th of April. 7 2 LITERATURE REVIEW 2.1 Kriging Kriging is the approach to estimate the attribute values at an unsampled location. The underlying premise of kriging is to establish a model that captures the spatial dependence or correlation among data points within a certain region. The spatial relationship between data points is described by a variogram or semivariogram, which quantifies the decline in a similarity between data points as the distance between them rises (Ryu et al., 2002). The variogram plays a crucial role in capturing the inherent spatial patterns and provides valuable insights for the kriging procedure. According to Goovaerts (1997), there are two main types of kriging: simple and ordinary. OK makes the assumption of second-order stationarity, a particular type of stationarity. This indicates that the spatial dependency (correlation between values at different places) depends primarily on the distance between them, while the statistical features of the variable, such as the mean and variance, remain constant throughout the study region. However, it should be noted that the covariance function exhibits stationarity. On the other hand, simple kriging technique incorporates a location-dependent estimate of the mean as a replacement for the constant mean. In simple kriging, the assumption is made that the mean value of the stationary random variable remains constant and is known. Basically, kriging provides one with two final maps: estimation and variance of the estimation error, showing the final produced estimation and the uncertainty of this modeling. In multivariate geostatistics, cokriging is introduced (Wackernagel, 2003). Cokriging is an enhanced version of kriging that operates on the same fundamental premise as kriging. However, it requires extra inputs such as auxiliary data, high variance in the data (which measures the covariance between primary and secondary variables), and temporal information (in the case of present trend in the data). These additional components enable cokriging to effectively capture the covariance between many linked regionalized variables. The method can be understood as a multivariate adaptation of kriging, where the secondary variable is utilized to estimate the first variable (Zhao et al., 2014). The fundamental motivation for including simulation in data analysis is the capacity to include several secondary variables with the primary variable, resulting in a more accurate data population. This improved accuracy is attributed to the reduced variance of estimation error compared to standard kriging methods if the dataset is heterotopic (Ochie and Rotimi, 2018). A heterotopic dataset in geostatistics refers to a dataset where data (first and second 8 variable) is not sampled at the same locations, resulting in a dataset with spatially disparate information (Araújo et al., 2021). Also, there are some important limitations while using any kriging technique. The kriging interpolator is limited to providing a local best estimate rather than a global best estimate. The approach frequently exhibits a tendency to overestimate the smaller values and underestimate the higher values, a phenomenon commonly referred to as the smoothing effect. 2.2 Independent Simulation The solution to the above-mentioned issues with kriging including smoothing effect has been resolved by the introduction of another methodology - simulation. Geostatistical simulation has the capability to generate many realizations or outcomes for one or several variables. In order to generate a range of distinct scenarios, the simulation method incorporates random variability (Lin et al., 2001). In order to conduct a simulation, it is typically necessary to have the probability distribution of the variable(s) as well as the geographical model (i.e., variogram) for the variable(s) as inputs. While the various outcomes generated by the simulation exhibit differences, they nonetheless adhere closely to the input distribution and input variogram. Overall, there are many different types of simulation including Sequential Gaussian Simulation (Hengl et al., 2010), Turning Bands Simulation (TBM) (Tompson et al., 1989), and Truncated Gaussian Simulation (TGS) (Beucher and Renard, 2016) etc. The SGS method utilizes the sequential approach to simulate Gaussian random fields, which is based on conditioning the simulation on the data obtained from a basic kriging exercise. While geostatistical modeling is generally regarded as a dependable approach, it has significant challenges in determining the concept of a moving neighborhood and accurately replicating short- scale continuity (Madani and Carranza, 2019). The suitability of this geostatistical simulation algorithm for univariate spatial modeling is satisfactory. However, the complexity of the modeling task increases when several cross- correlated variables are involved. One explanation is that the univariate simulation methodologies fail to adequately consider the dependency features among variables, resulting in inadequate outcomes (Madani and Carranza, 2019). Thus, the possible solution to this issue is the co- simulation which is discussed further. 2.3 Co-Simulation When simulating multi-element deposits independently, one does not consider the intrinsic correlation coefficient that exists between the variables. The co-simulation method, on the other 9 hand, takes into account this parameter by computing the cross-variograms (Maleki and Madani, 2016). There are many types of Gaussian co-simulation, however, first it is important to understand the concept behind them. In order to conduct a stochastic simulation, the variables must first be transformed into normal score values, which have a mean of zero and a variance of one. In general, the transformation can be carried out via Gaussian anamorphosis, which converts variables into conventional Gaussian form, or a quantile–quantile-based approach (Deutsch and Journel, 1998). The multivariate Gaussianity assumption is the foundation upon which multivariable geostatistical modeling is constructed (Journel and Huijbregts, 1978). Multivariable geostatistical modeling might include two or several variables. Within the context of this method, the variables ought to be independently transformed into a normal Gaussian distribution. After that, a co-simulation can be implemented over the normal score data, taking into consideration the cross-dependency functions. After the realizations have been acquired, they can then be back-transformed to their original scale in order to approximately reconstruct the linear multivariate relationships that existed between the variables (Battalgazy and Madani, 2019). As it was mentioned before, there are many techniques of gaussian co-simulation, but let us discuss two methods: Sequential Gaussian co-simulation and Turning Bands co-simulation. Sequential Gaussian co-simulation (Hengl et al., 2010) is a sequential approach that can be utilized to simulate a vector Gaussian random field, based on the correlation of structure and linear coregionalization model (Paravarzar et al., 2015). In theory, this method is considered to be completely accurate. Nevertheless, in practical applications, certain simplifications become necessary because of the processing limitations of cokriging when dealing with a high number of data points. This phenomenon occurs when there is a big amount of pre-existing data or when the number of places selected for simulation is large (Paravarzar et al., 2015). Turning Bands co-simulation’s (Tompson et al., 1989) primary principle is predicated on two steps. Firstly, several lines are drawn with random orientations. Secondly, a simulation of a one-dimensional Gaussian random field is conducted along each of these lines. In essence, the turning bands algorithm aims to simplify the task of simulation in three-dimensional (R3) or two- dimensional (R2) spaces by transforming it into a problem in one-dimensional space (R) (Emery, 2008). Moreover, few research and case studies were done to compare these two co-simulation techniques. According to Paravazar et al. (2015), the Turning Bands co-simulation algorithm demonstrated superior performance compared to the sequential Gaussian co-simulation algorithm in terms of accurately replicating the cross-correlation as determined by spatial continuity and statistical characteristics. 10 2.4 Principal component analysis (PCA) The joint modeling of several linked components in spatial analysis presents challenges. It is imperative to adequately represent the intricate interconnections among the variables and accurately replicate the observed spatial correlation in the generated realizations. Due to this rationale, several methodologies have been devised to convert the variables into univariate Gaussian or Multivariate Gaussian distributions, hence enabling the application of Gaussian Simulation techniques (Bolgkoranou, 2019). There are several techniques used for the factorization of variables for the abovementioned application, but in this work two methods will be described: PCA (Principal Component Analysis) and MAF (Minimum-Maximum Autocorrelation Factors). Principal Component Analysis (PCA) is a widely employed multivariate technique that aims to provide compact and more comprehensible sets of transformed variables by mitigating the presence of noise. The method described involves the decomposition of the variance-covariance matrix of variables into a set of independent transformed variables, referred to as principal components. These components are characterized by their ability to explain decreasing proportions of the overall variance observed in the data. In general, the majority of the variation is represented by the initial main components, resulting in a more concise illustration of the data. Nevertheless, it is worth mentioning that doing Principal Component Analysis (PCA) directly on the variance-covariance matrix of regionalized data, without taking into account spatial correlation, can be considered overly simplistic (Desbarats & Dimitrakopoulos, 2000). The complete decorrelation of principal components across all lags is solely attained when multivariate spatial data conforms to an "Intrinsic '' coregionalization model (Journel and Huijbregts, 1978). While this approach may not always adhere to precise applicability, it is not frequently employed in practical settings. 2.5 Minimum-Maximum Autocorrelation Factors (MAF) In order to mitigate deficiencies in the processing of remote sensing data through the utilization of spatial filtering and standard principal components, Switzer and Green (1984) introduced the “Minimum-Maximum Autocorrelation Factors”. The approach utilized in this methodology entails the implementation of a linear transformation on the variables that exhibit correlation. This transformation results in the conversion of the initial data into a domain where the variables are no longer correlated. The converted data can be calculated separately using kriging, bypassing the need for the linear model 11 of coregionalization. After the completion of the kriging process, the MAF estimates are subjected to a back-transformation, which restores their correlation in the original data space (Da Silva and Costa, 2014). The MAF technique is capable of transforming a group of variables into uncorrelated factors across various lag spacings. This is achieved under the condition that the covariance function of the variables can be fully described by a two-structure linear model of coregionalization (Ferreira, 2010). The aforementioned assumption may impose limitations, since the utilization of a model comprising only a nugget plus one structure or two structures could not be sufficient to accurately represent all cross-semivariograms (Bandarian and Mueller, 2008). The primary distinction between MAF and PCA lies in their respective assurances of spatial orthogonality. MAF ensures spatial orthogonality at two distinct distances in space, denoted as 𝑥1 and 𝑥2. In contrast, PCA only guarantees spatial orthogonality at a single distance, denoted as 𝑥, which encompasses all possible distances, including zero. The MAF technique involves the individual simulation and subsequent back-transformation of independent factors. This process generates conditional simulations of correlated qualities, which accurately replicate the cross- correlations observed in the original variables (Boluwade and Madramootoo, 2015). 2.6 Projection pursuit multivariate transform (PPMT) The aforementioned procedures, such as simulation via factorization (MAF and PCA) and standard co-simulation approaches, are suitable for situations where the bivariate relationship is linear. However, in circumstances of non-linearity, alternative methodologies, such as Projection Pursuit Multivariate Transform (PPMT) needs to be employed (Faria et al., 2021). Moreover, it is important to mention that in geometallurgy problems, such a linear bivariate relationship is rare and thus, MAF and PCA methods are not the fittest in such types of modeling. PPMT was introduced by Barnett et al. (2014) in order to come up with the solution of complex multivariate correlations for variables. The majority of geostatistical simulation techniques are dependent on the assumption of the multivariate Gaussian distribution. The PPMT aims to transform the data into a set of variables that exhibit no correlation and adhere to a multivariate Gaussian distribution. Subsequently, the application of geostatistical simulation techniques becomes feasible for the individual transformed variables, eliminating the necessity for cokriging or the linear model of coregionalization (Journel and Huijbregts, 1978) in order to replicate potential cross-correlations among them. According to Bassani et al. (2018), this is the procedure of applying PPMT: 1. The variables should be switched to a normal score. 12 2. The variables are subjected to a sphere transformation. 3. Identify a compelling visualization of the dataset. 4. The data should be normalized along the projection identified in step 3. 5. Continue iterating through steps (3) and (4) until the dataset exhibits a multivariate Gaussian distribution with an identical covariance matrix. However, one of the issues of PPMT, as well as of PCA and MAF is that they require the dataset to be homotopic, which is the rare case in real life modeling (Faria, 2021). Moreover, the second limitation of this technique is that the PPMT algorithm does not explicitly take into account limitations related to the summation or percentage of variables (Bassani et al., 2018). 2.7 Flow Anamorphosis Introduced by van den Boogaart and Tolosana-Delgado (2018), the flow anamorphosis technique is a version of Gaussian anamorphosis that may effectively convert original multivariate data into a multivariate normal distribution. Notably, this method remains invariant regardless of the log-ratio transform chosen (Talebi et al., 2018). The proposed approach involves transforming a kernel density estimate derived from the observed multivariate density into the density of a conventional multivariate normal distribution (Talebi et al, 2018). Flow anamorphosis has the capability to replicate intricate patterns observed in input data, such as the existence of outliers, multiple populations, nonlinearity, and heteroscedasticity (Talebi et al., 2018). Additionally, it possesses the property of invariance under the selection of log-ratio transformation, as well as the ability to generate spatially orthogonal factors, which simplifies the process of geostatistical simulation. Various statistical tests for assessing multivariate normality can be utilized to evaluate the conformity of the transformed data to a normal distribution (Korkmaz et al., 2014). The transformation is governed by two parameters: 1 and 2, which represent the starting and final spreads of the smoothing kernels used in the kernel density estimates, respectively. These values require tuning. The selection of an appropriate value for a parameter is contingent upon several factors, including the number of variables, sample size, and complexity of the input data. On the other hand, the parameter governs the extent of the altered distributions. The simulated outcomes are then converted back to the original domain using the inverse of the transformation function Flow Anamorphosis (Talebi et al., 2018). Next, according to Talebi et al. (2018), the figure illustrating the procedure of the flow anamorphosis technique is presented with the description of each step: 13 a. The scatterplot and marginal kernel smoothing density estimations are utilized to analyze the input data consisting of two simulated variables that exhibit a complex relationship. The properties of the data encompass the existence of outliers, non-linear associations, several populations, and heteroscedasticity. b. Illustrate the concurrent deformation of the foundational space. c. Show the ultimate arrangement of the modified data, correspondingly. d. Transformation of data (Various geostatistical algorithms can be employed to simulate the modified data in a multivariate normal space). e. Demonstrate the process of co-deformation, wherein data is transformed from a simulated multivariate normal space back to its original space. f. The final allocation of the simulated data. Figure 2. Geostatistical Simulation via FA (Talebi et al., 2018) However, same as PPMT, flow anamorphosis is also dependent on the presence of homotopic dataset and, is restricted to specific marginal distributions. 14 3 METHODOLOGY In this section, the theory of Hierarchical Sequential Gaussian co-simulation and “Acceptance-Rejection” methodology will be illustrated. 3.1 Hierarchical Sequential Gaussian Co-Simulation First, to understand the idea behind the “Hierarchical Sequential Gaussian co-simulation”, it is vital to figure out the working concept of Sequential Gaussian co-simulation. Thus, according to Journel and Huijbregts (1978), geostatistics involves the interpretation of 𝐾 variables [𝑍𝑘(𝑥), 𝑘 = 1 ] as realizations of K random variables, which are represented as 𝑍𝑘(𝑥), 𝑘 = 1 = 1 to K]. The collection of random variables in a given space can be represented as K random functions, which are sometimes denoted as [𝑍𝑘(𝑥), 𝑘 = 1 𝑡𝑜 𝐾], where k ranges from 1 to K]. Within the given theoretical framework, the task involves creating simultaneous instances of K random functions while ensuring the following conditions: 1) accurately replicating the distribution and variogram of each random function; 2) accurately replicating the cross-variograms between the random functions; and 3) appropriately honoring the sample values. The initial stage involves making the assumption that every random function 𝑍𝑘(𝑥) is strictly stationary and that there exists a random function 𝑌𝑘(𝑥), that satisfies this condition: 𝑌𝑘(𝑥) = 𝛷𝑘[𝑍𝑘(𝑥)], and 𝑌𝑘(𝑥)~𝑁(0,1), Where 𝛷𝑘 is a gaussian transform which is a one-to-one function. The variables that have been converted are commonly known as normal scores. The random functions [𝑌𝑘(𝑥), 𝑘 = 1 𝑡𝑜 𝐾] are inherently designed to possess the properties of rigorous stationarity and standard normal distribution. Furthermore, it is commonly assumed that these variables exhibit a joint multivariate normal distribution. The whole joint probability distribution of the set [𝑌𝑘(𝑥), 𝑘 = 1 𝑡𝑜 𝐾] can be determined when the mean, auto-covariance, and cross-covariance functions are known. The means are zero by virtue of their construction. The auto- and cross-covariance functions are connected to the auto- and cross-variograms through a mathematical relationship: 𝐶𝑘𝑘′(ℎ) = 𝐶𝑘𝑘′(0) − 𝛶𝑘𝑘′(ℎ), 𝑘 = 1 𝑡𝑜 𝐾, 𝑘′ = 1 𝑡𝑜 𝐾 In the above equation, 𝐶𝑘𝑘′(ℎ) 𝑎𝑛𝑑 𝛶𝑘𝑘′(ℎ) represent the covariance and cross- variogram, respectively, between the variables 𝑧𝑘 and 𝑧𝑘′ at a distance h. The further step to proceed the sequential gaussian simulation itself. The proposed algorithm is given by Verly (1993). 15 ● The subsequent stage involves the implementation of sequential simulation. The procedure is executed in the following manner: ● Randomly identify a place, denoted as 𝑥0, that has not previously undergone simulation. ● Utilizing the available normal scores from both observed and previously simulated data, calculate the conditional distribution of the random vector at the specified position, xo. The random vector is comprised of the random variables [𝑌𝑘(𝑥0), 𝑘 = 1 𝑡𝑜 𝐾], with the exception of those that have been sampled. ● Use the conditional distribution to make a realization. This realization is made up of simulated normal score vectors. ● Continue the procedure of simulation until each location have been covered. Now, let us look at the theory behind the hierarchy of variables and hierarchical sequential gaussian co-simulation itself. According to Almeida and Journel (1994), this technique incorporates spatial correlation by preserving the collocated values of previously simulated primary variables of a different kind while simulating a specific primary variable. Hence, it is imperative to establish a hierarchical structure of primary variables, commencing with the variable 𝑌1(𝑢) that holds the utmost significance. This hierarchy is based on the autocovariance 𝜌𝑌1 (ℎ) as well as the cross-covariances 𝜌𝑌𝑖𝑌𝑘 (ℎ)and 𝜌𝑌𝑖 𝐵𝑖 (ℎ), which play a crucial role in achieving accurate replication. The estimator of collocated cokriging can be interpreted as a comprehensive kriging method that incorporates all adjacent data points of the same type as the variable being estimated, along with a regression analysis that only utilizes collocated data points of different types, regardless of whether they are primary or secondary (Almeida and Journel, 1994). The execution of the simulation algorithm encompasses the subsequent stages: ● In order to establish a hierarchical structure of primary variables, we begin with the utmost significant variable, denoted as 𝑍𝑙(𝑢), and go downwards to the comparatively less significant variable, denoted as 𝑍𝑘(𝑢) ● For the sake of standardizing the variables, we will convert all major variables to their normal scores denoted as 𝑌𝑘(𝑢), and all secondary variables to their normal scores denoted as 𝐵𝑙(𝑢). ● In cases when secondary data (referred to as L) is not accessible at each node u where primary variables are to be simulated, it is necessary to pre-simulate them individually, taking into account the matching available data. The pre-simulation method does not take into consideration any interdependence among secondary variables. Nevertheless, this 16 approximation is deemed acceptable due to two reasons. Firstly, secondary variables are associated with subjective information. Secondly, they are typically sampled at a high frequency. ● In the context of a simulation grid, a random path can be defined as a traversal that visits each node u, u', u", and so on, only once. ● At node u, the objective is to calculate the Gaussian conditional distribution (ccdf) of 𝑌𝑙(𝑢) (u) given the n adjacent primary data of the same type 𝑦𝑖(𝑢𝛼), where α ranges from 1 to n. Additionally, only the collocated secondary data 𝑏𝑙(𝑢), where l ranges from 1 to L, are considered. ● At the same node u, we now aim to calculate the ccdf of the Gaussian distribution for the random variable 𝑌2(𝑢) (u)This calculation is based on the following information: the n data points of the same kind, denoted as 𝑦2(𝑢𝛼),, where α ranges from 1 to n; the previously simulated collocated value 𝑌𝑙(𝑢); and the collocated secondary data, denoted as 𝑏𝑙(𝑢),, where l ranges from 1 to L. ● Next, by the same logic, if K is > 2, we need to estimate ccdf of the Gaussian distribution for 𝑦3(𝑢). The primary data consists of n observations of the same type, denoted as 𝑦3(𝑢𝛼), where α ranges from 1 to n. Additionally, the previously simulated values 𝑌1(𝑢) and 𝑌2(𝑢) are also included as collocated simulated values. The collocated secondary data consists of 𝑏𝑙(𝑢), where l ranges from 1 to L. ● Proceed to the subsequent node, denoted as u', along the stochastic trajectory and iterate the preceding 3 instructions. ● And finally, Iterate through each node (u, u', u") within the simulation grid. 3.2 Acceptance-Rejection method Next, moving on to the acceptance-rejection methodology, it is a statistical technique utilized for the purpose of generating random variates from a given probability distribution. The process entails the acceptance or rejection of candidate values through a comparative analysis with a distribution that is relevant to the context (Zheng, 2001). The methodology described has been applied in many domains, such as geostatistical modeling, to effectively handle spatial uncertainty and examine intricate inequality constraint (Madani and Abulkhair, 2020). One instance of the acceptance-rejection methodology in geostatistical modeling pertains to its utilization in the examination of geographical heterogeneity in radon concentrations. Through the application of this approach, scholars have constructed a model that accurately 17 represents the arrangement of radon levels in space. This research resulted in significant knowledge pertaining to environmental and public health matters (Fanshawe and Diggle, 2011). One of the notable benefits of employing the acceptance-rejection methodology lies in its inherent adaptability when confronted with complex distributions and relationships. According to Zheng (2001), the utilization of this tool enables the development of random variates from distributions that differ from the standard, hence increasing its versatility in the context of geostatistical modeling. Furthermore, the methodology offers a structured framework for effectively managing geographical uncertainty, empowering researchers to make well-informed decisions through the study of spatial data. Nevertheless, the acceptance-rejection process exhibits certain drawbacks, particularly in situations where the distribution of the candidate significantly deviates from the target distribution. In instances of this nature, the efficacy of the approach may diminish, resulting in extended computing durations and diminished efficacy in producing random variables. Moreover, the methodology might require meticulous selection of relevant distributions and improving of parameters to attain optimal outcomes, hence presenting difficulties in real-world implementations. In brief, the acceptance-rejection methodology is a helpful technique in the field of geostatistical modeling. It provides a systematic strategy for addressing geographical uncertainty and analyzing intricate non-linear bivariate interactions. Researchers should exercise caution and acknowledge the limitations of this approach, especially when dealing with situations characterized by substantial disparities in distribution. 3.3 Proposed Methodology The Acceptance-Rejection approach is applied first in the methodology, and it is specifically designed to handle complicated, non-linear bivariate interactions. With this approach, candidate samples are generated, and their acceptance or rejection is determined by a criterion that is consistent with the intended bivariate relation, known as scattergram, guaranteeing that the underlying non-linear connections are followed. As an illustration of the working logic of this Acceptance-Rejection method, let us take an example of a random variable 1 (X-axis) and variable 2 (Y-axis) with them having a pretty good non-linear correlation of ~0.8 (Figure 3). 18 Then, after producing the scatterplot, there will be a clear visual non-linear relationship of variables. The co-simulation by the proposed method takes the simulated value of variable 1 at target grid node, for example 5%, and proceeds to simulate the variable 2 at the y-axis conditioned to the simulated value of variable 1 (5%). When the variable 2 is simulated, it is required for it to be in the interval of y-axis that it within the non-linear condition (from 25 to 50 % by y-axis) following the bivariate relationship manifesting in the scatterplot between these two variables (Figure 3). If this condition is held, the simulated value of variable 2 is “accepted”. On the other hand, in the case of not trapping into the interval of interest, the value is “rejected” and the value is re-simulated until the condition is met and accommodate in the interval, and therefore the simulated value of variable 2 is “accepted”. Figure 3. Example of the Scatter Plot for Acceptance-Rejection method 19 The next step is to use Exploratory Data Analysis (EDA). This is a critical stage in identifying trends, abnormalities, and the structure of the dataset. An appropriate simulation approach will be chosen based on the knowledge gathered from the EDA. The selected method will next be used in the simulation phase, where the goal will be to replicate the dynamics and patterns that were seen and recognized in the preceding phases. This will enable a thorough comprehension and modeling of the behaviors included in the dataset. 3.4 Data Requirements The data used in this research is a synthetic geometallurgical drill hole dataset presented in the work by Garrido et al. (2020). In this section, the main concept of authors will be explained and thus, the formation of the dataset as well. This article is about the simulation of synthetic exploration and geometallurgical database of porphyry copper deposits for educational purposes. It proposes a methodology to simulate geometallurgical data with geostatistical tools, preserving the coherent relationship among primary attributes, mineralogy, and response attributes. The simulated geometallurgical drill holes generated are realistic and consistent with input statistics, geology, and mineralogy, producing realistic processing metallurgical performance responses. The authors state that the access to real geometallurgical data is very limited in practice, making it difficult for practitioners, researchers, and students to test methods, models, and reproduce results in the field of geometallurgy. Therefore, the author proposes a methodology to simulate geometallurgical data with geostatistical tools preserving the coherent relationship among primary attributes, such as grades and geological attributes, with mineralogy and some response attributes, for example, grindability, throughput, kinetic flotation performance, and recovery. The proposed methodology by the authors is based on three main components: 1. Definition of spatial relationship between geometallurgical units: This involves defining the spatial relationship between geometallurgical units, such as ore types, lithologies, alteration types, and mineralization styles, using geological and geostatistical tools. 2. Co-simulation of regionalized variables with geometallurgical coherence: This involves cosimulating regionalized variables, such as grades, geological attributes, and mineralogy, with geometallurgical coherence using geostatistical tools. 3. Simulation of georeferenced drill holes based on geometrical and operational constraints: This involves simulating georeferenced drill holes based on geometrical and 20 operational constraints, such as drill hole spacing, orientation, and length, using geostatistical tools. By following these three main components, the authors propose a reproducible methodology for the simulation of a synthetic geometallurgical drill holes dataset, with special interest in preserving the coherence between geology, mineralogy, and grades. Response attributes were included in the drill hole database, comminution process, and flotation performance. 3.5 Expected Results As it was mentioned earlier, the main idea of this research is to deal with complex non- linear bivariate relationships between any geometallurgical properties such as minerals (bornite, pyrite, chalcocite, chalcopyrite, molybdenite, tennantite, clays), work index (the amount of energy required to crush the material) or recovery. Thus, some preliminary exploratory data analysis was conducted. As can be seen from Table 1, it is difficult to say there is a strong linear relationship between any of the parameters (excluding recovery and chalcocite or molybdenite and bornite). Table 1. Linear Correlation between Geometallurgical parameters clays chalcoci te bornite chalcop yrite tennant ite molybd enite pyrite bwi recovery clays 1 -0.374 0.4458 -0.0376 -0.0501 -0.0182 -0.2676 -0.0157 0.304 chalcocite -0.37 1 0.1775 -0.0180 0.2434 0.1872 0.1731 0.085 -0.5859 bornite 0.445 0.1775 1 0.3364 0.3961 0.6272 -0.6038 0.066 0.4941 chalcopyrite -0.0376 -0.0180 0.3364 1 0.5073 0.6951 -0.0638 -0.0416 0.3195 tennantite -0.0501 0.2434 0.3961 0.5073 1 0.5511 -0.0047 0.1556 0.01949 molybdenite -0.0182 0.1872 0.6272 0.6951 0.5511 1 -0.335 0.02017 0.3167 pyrite -0.2676 0.1731 -0.6038 -0.0638 -0.0047 -0.335 1 -0.0783 -0.5244 bwi -0.0157 0.085 0.066 -0.0416 0.1556 0.02017 -0.0783 1 -0.0505 recovery 0.304 -0.5859 0.4941 0.3195 0.01949 0.3167 -0.5244 -0.0505 1 However, if we look at table 2, which illustrates the spearman correlation (non-linear) between the variables, it can be observed that there is pretty good overall non-linear correlation. 21 Table 2. Non-Linear Correlation between Geometallurgical parameters clays chalcoci te bornite chalcop yrite tennant ite molybd enite pyrite bwi recovery clays 1 -0.18 0.02 -0.02 0.03 0.01 -0.18 -0.01 0.2 chalcocite -0.18 1 0.07 -0.13 0.03 -0.01 0.09 0.11 -0.74 bornite 0.02 0.07 1 0.23 0.16 0.6 -0.25 0.05 0.18 chalcopyrite -0.02 -0.13 0.23 1 0.19 0.48 -0.07 -0.05 0.26 tennantite 0.03 0.03 0.16 0.19 1 0.27 -0.03 0.05 0.01 molybdenite 0.01 -0.01 0.6 0.48 0.27 1 -0.21 -0.02 0.16 pyrite -0.18 0.09 -0.25 -0.07 -0.03 -0.21 1 -0.07 -0.27 bwi -0.01 0.11 0.05 -0.05 0.05 -0.02 -0.07 1 -0.08 recovery 0.2 -0.74 0.18 0.26 0.01 0.16 -0.27 -0.08 1 Thus, as it will be discussed, some factorization techniques like Principal Component Analysis (PCA) or Minimum/Maximum autocorrelation factors (MAF) may be suboptimal for such non-linear relationships between these geometallurgical properties because of their constraints. Therefore, it is expected that with the help of acceptance-rejection technique and further implementation of sequential Gaussian co-simulation, this case’s complex bivariate relationship will be modeled, and the outcome of the simulation will be following the assigned non-linearity. 22 4 RESULTS 4.1 Case study The dataset in this research is synthetic, which was introduced in Garrido et al (2020). The approach employed in this study involves the representation of geometallurgical data through the identification of variable types, the establishment of a uniform database, and the simulation of geological variables and their corresponding responses for comminution and flotation processes. This methodology guarantees a cohesive correlation between primary characteristics, mineral composition, and response parameters. The dataset is homotopic including many categorical and continuous geometallurgical variables. In this work, both conventional hierarchical sequential Gaussian co-simulation and proposed method, hierarchical sequential Gaussian co-simulation using acceptance-rejection technique will be applied to simulate two variables: chalcocite (𝐶𝑢2𝑆) and recovery (Rec). As can be seen from the location maps below (Figures 4 and 5), there is quite good negative correlation between these two variables, as in the central left part of the dataset, high values of Chalcocite indicate low values of recovery in that region. The same can be observed in the central right part of the deposit. 23 Figure 4. Location Map of Chalcocite 24 Figure 5. Location Map of Recovery 4.2 Exploratory data analysis The sampling pattern within the deposit was irregular in specific areas. This irregularity is apparent in the location map in Figures 4 and 5, showing a concentration of samples in the deposit’s middle part. To mitigate the sampling irregularity, cell-declustering was employed based on methods from David (1977) and Deutsch (1989), using a cell size of 30m x 30m x 15m. This approach was chosen after examining the impact of various cell sizes on the overall statistical parameters of the dataset. The adjustments made before and after cell-declustering are documented in Table 3, highlighting a reduction in the average of Chalcocite grade and a slight increase in the average Recovery grade, due to the higher weights assigned to samples from the deposit’s periphery. 25 Table 3. Statistical parameters of the dataset before and after declustering Statistical parameter Before declustering After declustering Chalcocite (%) Recovery (%) Chalcocite (%) Recovery (%) Mean 0.064 88.08 0.06 88.09 Standard deviation 0.17 5.3 0.17 5.30 Coefficient of variation 2.66 0.06 2.83 0.06 Linear Correlation coefficient -0.59 -0.59 Non-Linear Correlation coefficient -0.74 -0.74 Minimum 0 60.335 0 60.336 Maximum 2.32 94.147 2.705 94.147 Moreover, as it was mentioned earlier, the synthetic dataset includes many geometallurgical parameters, however, due to a proper negative non-linear correlation and acceptable negative linear correlation, Chalcocite and Recovery were chosen for this research (Table 3) to examine the proposed co-simulation algorithm. The scatter plot between Chalcocite and Recovery reveals a pronounced and clearly seen non-linear relationship in their bivariate distribution, as illustrated in Figure 6. Chalcocite predominantly exhibits low content across the samples, in contrast to the more evenly distributed Recovery. 26 Figure 6. Scatter plot between Recovery and Chalcocite 4.3 Variogram analysis The process of variogram analysis for co-simulation mandates converting all variables into their normal scores, considering declustering weights as obtained beforehand. According to Emery & Peláez (2011), data transformation to Gaussian space entails turning the initial data into a format that conforms to a Gaussian distribution. Ensuring the validity of the findings generated via sequential Gaussian co-simulation is dependent on this transformation. Furthermore, after the normal score transformation, the suitability for co-simulation is assessed by evaluating the collocated scatter plot of normal scores for Chalcocite and Recovery, depicted in Figure 7. The correlation between Gaussian Chalcocite and Gaussian Recovery distributions is quite complex, 27 evidenced by a non-linear correlation coefficient of -0.61, affirming their bivariate Gaussian compatibility. Figure 7. Scatter Plot of Chalcocite and Recovery in Normal Scores with their Marginal Distribution Bivariate Gaussianity is evaluated using lagged scatter plots for each variable’s normal scores. According to Krupskii and Genton (2017), the lagged scatter plot of normal scores is a visual tool used in several disciplines, including statistics and geology, to evaluate the multi- Gaussianity. When analyzing the scatter plot of normal scores, it is expected to display an elliptical form if the assumption of joint normality is true. The lagged scatter plot of both Gaussian Recovery and Gaussian Chalcocite is shown in Figures 8 and 9 respectively. In these graphs, X-axis stand for Gaussian values of the underlying variable at the distance x, while Y-axis is illustrating Gaussian values of the underlying variable at the distance 𝑥 + ℎ, where h is a lag distance (distance between pairs). As can be seen from these figures, both distributions are presenting an elliptical shape. This evidence suggests that bivariate Gaussian assumptions is held, allowing for sequential Gaussian co-simulation in this research. 28 Figure 8. Lagged Scatter Plot of Normal Score of Recovery Figure 9. Lagged Scatter Plot of Normal Score of Chalcocite 29 Next, anisotropy analysis through sample variograms in various directions was considered. According to Cayon (2010), deriving sample variograms in different directions is a fundamental approach to assess anisotropy in geostatistics and various geological studies. This method allows researchers to quantify the directional dependence of spatial structures, characterize isotropy, and determine the extent of anisotropy in diverse environmental and geological contexts. In our case, no horizontal anisotropy is detected, indicating a horizontally isotropic variogram. However, vertical anisotropy is identified, leading to the calculation of experimental direct and cross variograms in both vertical and horizontal orientations over the normal score data. A two- structured linear model of coregionalization is then fitted semi-automatically, as illustrated in Figure 10. Below is the formula of variograms: ( 𝛾 𝑅𝑒𝑐 𝛾 𝑅𝑒𝑐/𝐶𝑢 2 𝑆 𝛾 𝐶𝑢2𝑆/𝑅𝑒𝑐 𝛾 𝐶𝑢2𝑆 )=( 0.02 −0.02 −0.02 0.14 ) 𝑁𝑢𝑔𝑔𝑒𝑡 + ( 0.09 −0.09 −0.09 0.35 ) 𝑆𝑝ℎ(78.3𝑚, 78.3𝑚, 111.7𝑚) + ( 0.9 −0.49 −0.49 0.51 ) 𝐸𝑥𝑝(225.8𝑚, 225.8𝑚, 152.8𝑚) (1) 30 Figure 10. Direct and cross-variograms of the Recovery and Chalcocite in normal scores in both the horizontal (black) and vertical (blue) directions; points: experimental variogram, and solid line: variogram model (derived linear model of coregionalization) 4.4 The outcomes of Conventional and Proposed methods The modeling of Chalcocite and Recovery grades employs the sequential Gaussian co- simulation algorithm, as well as proposed acceptance – rejection method, both of which were complemented by inverse transform sampling. To obtain the realizations by simulation, the grid was created taking into consideration the dimensions of the sampling pattern. Figures 11 and 12 display random realizations of the proposed and conventional algorithm, respectively. The figures highlight high chalcocite concentrations in the deposit’s central-left area, with relatively low Recovery grades in the same region. Vice versa, in those north, west and east parts, elevated grades of recovery signal about absence of chalcocite, showing that simulation was able to successfully reproduce negative non-linear correlation between the variables. Also, these findings align with the location maps in Figures 4 and 5. Moreover, the causes of nugget effect can be seen as some disturbance is present all over the realizations, however, disregarding the nugget effect can lead to significant biases in the assessment and prediction of geometric phenomena (Kaufman & Shaby, 2013). 31 Figure 11. Random realization of Chalcocite and Recovery produced by the proposed algorithm based on hierarchical sequential Gaussian simulation using acceptance-rejection method 32 Figure 12. Random realization of Chalcocite and Recovery produced by the conventional hierarchical sequential Gaussian simulation algorithm The effectiveness of the proposed algorithm is further validated by comparing its realizations with those from conventional sequential Gaussian co-simulation using identical 33 parameters. Since the conventional approach cannot replicate the non-linear relationship of the original dataset, resulting in the overestimation of secondary variables, the proposed algorithm adjusts for this by re-simulating values outside the defined constraint (Madani and Abulkhair, 2020), which in our case is non-linearity. Consequently, the E-type map produced by the proposed and conventional methods show slightly lower Recovery grades across the region, except for the central area, as demonstrated in Figures 13 and 14. Figure 13. E-type realization of Proposed Approach of Recovery based on hierarchical sequential Gaussian simulation using acceptance-rejection method 34 Figure 14. E-type realization of Conventional Approach of Recovery based on conventional hierarchical sequential Gaussian simulation In addition, variance maps of both variables were obtained after using both approaches. As can be seen in the figures below, for Chalcocite (figures 15 and 17) the variance is close to 0 in almost all regions except for down-left part of the deposit. As well as in Recovery case (figures 16 and 18), the same region has very high variance. The reason behind it is the fact that this region was left unsampled (figures 4 and 5). 35 Figure 15. Variance Map of Chalcocite obtained from conventional hierarchical Gaussian co-simulation Figure 16. Variance Map of Recovery obtained from conventional hierarchical Gaussian co-simulation 36 Figure 17. Variance Map of Chalcocite obtained from proposed hierarchical Gaussian co-simulation using acceptance-rejection method Figure 18. Variance Map of Recovery obtained from proposed hierarchical Gaussian co-simulation using acceptance-rejection method 37 4.5 Geostatistical validation of both methods The geostatistical comparison of algorithms necessitates validation through both histogram and variogram assessments, encompassing first- and second-order statistical analyses. The objective for the proposed algorithm in this research is to accurately replicate the non-linearity between Chalcocite and Recovery and reproduce the better global statistical parameters and bivariate relationship of the original dataset compared to the conventional one. Validation therefore involves an examination of how well each algorithm reproduces global statistics and bivariate relationships. 4.5.1 Reproduction of global statistical parameters Now, the first step of validation of the results is assessing the reproduction of global statistical parameters, for which the primary confirmation step is to compare the mean, standard deviation, and correlation coefficients of the realizations of two methods against those of the original dataset. This comparison is quantified by calculating the average mean grade, standard deviation, coefficient of variation for Chalcocite and Recovery, along with their average correlation coefficients from the realizations generated by both conventional and proposed co- simulation algorithms. Table below showcases these statistical parameters against the statistics of the original dataset (after declustering). Table 4. Reproduction of the average statistical parameters over 100 realizations Statistical parameter Original dataset Conventional algorithm Proposed algorithm Chalcocite (%) Recovery (%) Chalcocite (%) Recovery (%) Chalcocite (%) Recovery (%) Mean 0.06 88.09 0.06 87.5 0.06 87.51 Standard deviation 0.17 5.30 0.17 5.95 0.17 5.93 Coefficient of variation 2.83 0.06 2.83 0.068 2.83 0.068 Non-linear Correlation coefficient -0.74 -0.68 -0.73 Linear Correlation coefficient -0.59 -0.51 -0.58 As can be seen from table 4, the results of mean, standard deviation and coefficient of variation are almost identical between realizations of Chalcocite and Recovery, which are 38 produced by two methods. However, there is a considerable difference in linear and non-linear correlation. Figure 19 illustrates the difference in non-linear correlation between original dataset and outputs of Proposed and Conventional algorithms. Ultimately, the proposed algorithm achieves superior reproduction of original correlation coefficient, evidencing its effective reproduction of the non-linearity between Chalcocite and Recovery as visualized in Figure 23. A more detailed exploration of this non-linearity replication is discussed in further section. Figure 19. Histogram of non-linear correlation coefficients between Chalcocite and Recovery from 100 realizations produced by conventional and proposed co-simulation algorithms; solid line: mean of correlation coefficient obtained from 100 realizations, dashed line: original correlation coefficient Moreover, according to Madani and Abiulkhair (2020), a key concern about the proposed algorithm centers on its capability to accurately reproduce the marginal distribution of the secondary variable, particularly due to the re-simulation of values that breach a non-linear relationship. Utilizing a histogram plot of one realization per variable offers a means to verify the fidelity of marginal distributions in comparison to the original histograms. For both variables, Chalcocite and Recovery, the outcome from both algorithms is almost identical, bolstering the validity of the comparison as depicted in the figure below. The histogram derived from a single 39 Chalcocite and Recovery realization mirrors the original marginal distribution, validating the efficacy of co-simulation methods (Abulkhair, 2021). However, validating histograms based on a single realization may not suffice. Comprehensive histogram validation is necessary to confirm that the proposed algorithm can replicate the secondary variable’s marginal distribution, even when employing inverse transform sampling for re-simulation (Abulkhair, 2021). Figure 20. Comparison of Chalcocite histogram reproduction by conventional and proposed methods for one realization 40 Figure 21. Comparison of Recovery histogram reproduction by conventional and proposed methods for one realization 4.5.2 Reproduction of bivariate relationship The methodology presented in this study is vital when modeling variables with a non-linear relationship in their bivariate relationship, especially significant for geometallurgical variables as there are many of them and frequently, these variables depend on each other in such non-linear manner. As discussed in the previous part, statistically the correlation of simulated variables of two methods differs massively, also the scatter plot reveals a significant difference (as depicted in Figures 22 and 23). The conventional approach exhibits a poor replication of the bivariate relationship between Chalcocite and Recovery (Figure 22), with some simulated values (blue) exceeding the non-linearity defined by the original values (red). The reason behind it is explained in theory in the Literature review part. Basically, the proposed algorithm incorporates an additional mechanism through inverse transform sampling, which mandates the re-simulation of all points that surpass a constrain, a non-linearity in our case, ensuring that the proposed co-simulation method accurately reflects the non-linearity between 41 variables, as evidenced in Figure 23, where scatter plots of both simulated (blue) and original (red) values are juxtaposed. Figure 22. Reproduction of the bivariate relationship between Chalcocite and Recovery with Conventional method 42 Figure 23. Reproduction of the bivariate relationship between Chalcocite and Recovery with Proposed method However, even though the results obtained from proposed methodology clearly outperform the realization from conventional co-simulation, there are still some points that are out of the non- linearity bivariate conditions. As can be seen in the figure 23, in the region of high Chalcocite, the simulated data (red) is clearly going out of the non-linearity conditions. The reason behind it is lack of dataset in this region of scatterplot, so there is a very narrow Recovery condition for Chalcocite, thus convergence may not happen properly. 43 5 CONCLUSIONS In this thesis, the focus was on refining resource modeling in mining engineering by addressing the complex nonlinear relationships prevalent in ore deposits. The study critiques the traditional emphasis on ore grade, advocating for a broader perspective that encompasses geometallurgical characteristics for sustainable and efficient mining operations. Moreover, in this work, the inadequacy of conventional methods like Gaussian co-simulations, PCA, and MAF in capturing the intricate nonlinear interactions within ore bodies were identified. The research proposes an advanced methodology, integrating hierarchical sequential Gaussian co-simulation with the acceptance-rejection method, tailored for complex bivariate relationships in synthetic geometallurgical drill hole datasets. The approach aims to outperform conventional simulation techniques by accurately representing the spatial and statistical complexities of ore deposits, especially in non-linear contexts. As it was illustrated in Results part, the proposed approach outperformed conventional one in reproducing the global statistics and bivariate relationship of original dataset. The thesis emphasizes the shift in mining practices from grade-centric to geometallurgically informed decisions, aligning with modern industry requirements for detailed understanding of ore properties. The acceptance-rejection method is highlighted for its ability to handle non-linear and complex geometallurgical relationships, improving the precision of resource estimation and mine planning. A synthetic dataset mimicking real-world geometallurgical conditions serves as the testbed for evaluating the proposed methodology in this work. The study provides a comprehensive analysis of this dataset, mainly Chalcocite and Recovery, demonstrating the method's effectiveness in accurately modeling the intricate relationships between geometallurgical parameters. Therefore, the findings underscore the limitations of traditional resource modeling approaches and showcase the potential of combining hierarchical sequential Gaussian co- simulation with acceptance-rejection sampling to enhance the modeling accuracy of complex geological phenomena. This thesis contributes significantly to the field of mining engineering by offering a sophisticated tool for geostatistical modeling, emphasizing the need for a comprehensive approach to understand and simulate the multifaceted nature of ore deposits. It not only addresses existing challenges in resource modeling but also lays the groundwork for future advancements in the domain, promoting more effective and sustainable mining practices. 44 For example, in the case of this deposit, one should analyze the negative correlation between Chalcocite and Recovery and thus, should adapt to different Recoveries throughout the mine, and possibly change the mineral processing route depending on the presence and absence of Chalcocite (low or high Recovery). However, some further research on this topic and the development of the Acceptance- Rejection approach still should take place. For example, this thesis faced some inaccuracy in output of the realizations of this method due to absence of properly sampled data in the tail of the scatterplot in the figure 23. Thus, the method still can be upgraded with solution of such issues. 45 REFERENCE LIST Abulkhair, S., 2021. A NOVEL APPROACH FOR GEOSTATISTICAL MODELING AND PRODUCTION SCHEDULING OF IRON DEPOSITS WITH INEQUALITY CONSTRAINTS. Nazarbayev University. https://doi.org/10.13140/RG.2.2.33215.64165/1 Almeida, A. S., & Journel, A. G. (1994, July). Joint simulation of multiple variables with a Markov-type coregionalization model. Mathematical Geology, 26(5), 565–588. https://doi.org/10.1007/bf02089242 Araújo, C. da P., Bassani, M.A.A., Koppe, V.C., Costa, J.F.C.L., Soares, A. de O., 2021. Geostatistical simulations with heterotopic hard and soft data without modeling the linear model of coregionalization. REM - International Engineering Journal 74, 269– 278. https://doi.org/10.1590/0370-44672020740075 Asa, E., Saafi, M., Membah, J., & Billa, A. (2012, January). Comparison of Linear and Nonlinear Kriging Methods for Characterization and Interpolation of Soil Data. Journal of Computing in Civil Engineering, 26(1), 11–18. https://doi.org/10.1061/(asce)cp.1943-5487.0000118 Bandarian, E., & Mueller. (2008). Reformulation of Direct Minimum/Maximum Autocorrelation Factors as a Generalised Eigenvalue Problem. Research Online. Retrieved November 20, 2023, from https://ro.ecu.edu.au/ecuworks/1086/. Barnett, R. M., Manchuk, J. G., & Deutsch, C. V. (2013, November 13). Projection Pursuit Multivariate Transform. Mathematical Geosciences, 46(3), 337–359. https://doi.org/10.1007/s11004-013-9497-7 Barnett, R. M., Manchuk, J. G., & Deutsch, C. V. (2016, October 20). The Projection-Pursuit Multivariate Transform for Improved Continuous Variable Modeling. SPE Journal, 21(06), 2010–2026. https://doi.org/10.2118/184388-pa Bassani, M. A. A., Coimbra Leite Costa, J. F., & Deutsch, C. V. (2018, June 19). Multivariate geostatistical simulation with sum and fraction constraints. Applied Earth Science, 127(3), 83–93. https://doi.org/10.1080/25726838.2018.1468145 Battalgazy, & Madani. (2019, November 4). Stochastic Modeling of Chemical Compounds in a Limestone Deposit by Unlocking the Complexity in Bivariate Relationships. Minerals, 9(11), 683. https://doi.org/10.3390/min9110683 https://doi.org/10.1007/bf02089242 https://doi.org/10.3390/min9110683 46 Beucher, H., Renard, D., 2016. Truncated Gaussian and derived methods. Comptes Rendus. Géoscience. https://doi.org/10.1016/j.crte.2015.10.004 Boluwade, A., & Madramootoo, C. A. (2015, December). Geostatistical independent simulation of spatially correlated soil variables. Computers & Geosciences, 85, 3–15. https://doi.org/10.1016/j.cageo.2015.09.002 David, M, 1977. Geostatistical ore reserve estimation, 384 p (Elsevier: New York). da Silva, C. Z., & Costa, J. F. (2014, April 8). Minimum/maximum autocorrelation factors applied to grade estimation. Stochastic Environmental Research and Risk Assessment, 28(8), 1929–1938. https://doi.org/10.1007/s00477-014-0879-2 Deutsch, J. L., Palmer, K., Deutsch, C. V., Szymanski, J., & Etsell, T. H. (2015, July 4). Spatial Modeling of Geometallurgical Properties: Techniques and a Case Study. Natural Resources Research, 25(2), 161–181. https://doi.org/10.1007/s11053-015-9276-x Deutsch, C V, 1989. DECLUS: a fortran 77 program for determining optimum spatial declustering weights. Computers & Geosciences, 15(3):325-332. Ding, X. L., Li, Z. W., Zhu, J. J., Feng, G. C., & Long, J. P. (2008, September 3). Atmospheric Effects on InSAR Measurements and Their Mitigation. Sensors, 8(9), 5426–5448. https://doi.org/10.3390/s8095426 Dominy, S. C., O’Connor, L., Parbhakar-Fox, A., Glass, H. J., & Purevgerel, S. (2018, December 1). Geometallurgy—A Route to More Resilient Mine Operations. Minerals. https://doi.org/10.3390/min8120560 Emery, X., & Peláez, M. (2011, April 28). Assessing the accuracy of sequential Gaussian simulation and cosimulation. Computational Geosciences, 15(4), 673–689. https://doi.org/10.1007/s10596-011-9235-5 Fanshawe, T. R., & Diggle, P. J. (2011, August 20). Bivariate geostatistical modelling: a review and an application to spatial variation in radon concentrations. Environmental and Ecological Statistics, 19(2), 139–160. https://doi.org/10.1007/s10651-011-0179-7 Faria, P. H., Coimbra Leite Costa, J. F., & Arcari Bassani, M. A. (2021, March 12). Multivariate geostatistical simulation with PPMT: an application for uncertainty measurement. Applied Earth Science, 130(3), 174–184. https://doi.org/10.1080/25726838.2021.1892364 https://doi.org/10.1016/j.cageo.2015.09.002 https://doi.org/10.1007/s11053-015-9276-x https://doi.org/10.3390/min8120560 47 Ferreira, J. (2010). Comparison and application of three decorrelation methods PCA, MAF and ACDC. Research Online. Retrieved November 20, 2023, from https://ro.ecu.edu.au/theses_hons/1217/. Garrido, M., Sepúlveda, E., Ortiz, J., & Townley, B. (2020, May 14). Simulation of Synthetic Exploration and Geometallurgical Database of Porphyry Copper Deposits for Educational Purposes. Natural Resources Research, 29(6), 3527–3545. https://doi.org/10.1007/s11053-020-09692-6 Goovaerts, P. (1997). Geostatistics for Natural Resources Evaluation. Oxford University Press. Gu, X., Wang, X., Liu, Z., Zha, W., Xu, X., Zheng, M., 2020. A Multi-Objective Optimization Model Using Improved NSGA-II for Optimizing Metal Mines Production Process. IEEE Access 8, 28847–28858. https://doi.org/10.1109/access.2020.2972018 Hengl, T., Heuvelink, G.B.M., van Loon, E.E., 2010. On the uncertainty of stream networks derived from elevation data: the error propagation approach. Hydrology and Earth System Sciences 14, 1153–1165. https://doi.org/10.5194/hess-14-1153-2010 Hosseini, S. A., & Asghari, O. (2018, April 26). Multivariate Geostatistical Simulation on Block- Support in the Presence of Complex Multivariate Relationships: Iron Ore Deposit Case Study. Natural Resources Research, 28(1), 125–144. https://doi.org/10.1007/s11053-018-9379-2 Hunt, J., Berry, R., Becker, M., & Baumgartner, R. (2019, December 1). A Special Issue Dedicated to Geometallurgy: Preface. Economic Geology, 114(8), 1473–1479. https://doi.org/10.5382/econgeo.4688 Isaaks, E.H. and Srivastava, R.M. (1989) An Introduction to Applied Geostatistics. Oxford University Press, New York, 413. - References - Scientific Research Publishing. (n.d.). https://www.scirp.org/(S(czeh2tfqw2orz553k1w0r45))/reference/referencespapers.asp x?referenceid=1876009 Journel, A. G., & Huijbregts, Ch. J. (1978). Mining Geostatistics. Blackburn Press. Kaufman, C. G., & Shaby, B. A. (2013, February 7). The role of the range parameter for estimation and prediction in geostatistics. Biometrika, 100(2), 473–484. https://doi.org/10.1093/biomet/ass079 48 Korkmaz, S., Goksuluk, D., & Zararsiz, G. (2014). MVN: An R Package for Assessing Multivariate Normality. The R Journal, 6(2), 151. https://doi.org/10.32614/rj-2014- 031 Krupskii, P., Genton, M.G., 2017. Factor copula models for data with spatio-temporal dependence. Spatial Statistics 22, 180–195. https://doi.org/10.1016/j.spasta.2017.10.001Lin, Y. P., Chang, T. K., & Teng, T. P. (2001, July 14). Characterization of soil lead by comparing sequential Gaussian simulation, simulated annealing simulation and kriging methods. Environmental Geology, 41(1–2), 189–199. https://doi.org/10.1007/s002540100382 Madani, N. (2023). Geostatistical Analysis in Geometallurgy. [PowerPoint slides]. Madani, N., & Abulkhair, S. (2020, July 25). A hierarchical cosimulation algorithm integrated with an acceptance–rejection method for the geostatistical modeling of variables with inequality constraints. Stochastic Environmental Research and Risk Assessment, 34(10), 1559–1589. https://doi.org/10.1007/s00477-020-01838-5 Madani, N., & Carranza, E. J. M. (2019, August 24). Co-simulated Size Number: An Elegant Novel Algorithm for Identification of Multivariate Geochemical Anomalies. Natural Resources Research, 29(1), 13–40. https://doi.org/10.1007/s11053-019-09547-9 Maleki, & Madani. (2016). Multivariate geostatistical analysis: an application to ore body evaluation. Iranian Journal of Earth Sciences, 173–184(8). Malliaris, A. G. T., Corazza, M., & Scalco, E. (2006). Nonlinear Bivariate Comovements of Asset Prices: Theory and Tests. SSRN Electronic Journal. https://doi.org/10.2139/ssrn.1021802 Ochie, K. I., & Rotimi, O. J. (2018, August 6). Geostatistics – Kriging and Co-Kriging Methods in Reservoir Characterization of Hydrocarbon Rock Deposits. All Days. https://doi.org/10.2118/193483-ms Rincon, J., Gaydardzhiev, S., & Stamenov, L. (2019, January). Coupling comminution indices and mineralogical features as an approach to a geometallurgical characterization of a copper ore. Minerals Engineering, 130, 57–66. https://doi.org/10.1016/j.mineng.2018.10.007 https://doi.org/10.32614/rj-2014-031 https://doi.org/10.32614/rj-2014-031 https://doi.org/10.2139/ssrn.1021802 49 Ryu, J. S., Kim, M. S., Cha, K. J., Lee, T. H., & Choi, D. H. (2002, May). Kriging interpolation methods in geostatistics and DACE model. KSME International Journal, 16(5), 619– 632. https://doi.org/10.1007/bf03184811 Switzer, P., & Green, A. (1984, January 1). Min/Max Autocorrelation factors for multivariate spatial imagery: Technical Report 6. ResearchGate. https://www.researchgate.net/publication/239665649_MinMax_Autocorrelation_facto rs_for_multivariate_spatial_imagery_Technical_Report_6 Talebi, H., Mueller, U., Tolosana-Delgado, R., & van den Boogaart, K. G. (2018, August 31). Geostatistical Simulation of Geochemical Compositions in the Presence of Multiple Geological Units: Application to Mineral Resource Evaluation. Mathematical Geosciences, 51(2), 129–153. https://doi.org/10.1007/s11004-018-9763-9 Tompson, A.F.B., Ababou, R., Gelhar, L.W., 1989. Implementation of the three‐dimensional turning bands random field generator. Water Resources Research 25, 2227–2243. https://doi.org/10.1029/wr025i010p02227 Van den Boogaart, K. G., & Tolosana-Delgado, R. (2018). Predictive Geometallurgy: An Interdisciplinary Key Challenge for Mathematical Geosciences. Handbook of Mathematical Geosciences, 673–686. https://doi.org/10.1007/978-3-319-78999-6_33 Verly, G. W. (1993). Sequential Gaussian Cosimulation: A Simulation Method Integrating Several Types of Information. Quantitative Geology and Geostatistics, 543–554. https://doi.org/10.1007/978-94-011-1739-5_42 Zhao, S., Zhou, Y., Wang, M., Xin, X., & Chen, F. (2014). Thickness, porosity, and permeability prediction: comparative studies and application of the geostatistical modeling in an Oil field. Environmental Systems Research, 3(1), 7. https://doi.org/10.1186/2193-2697-3-7 Zheng, L. (2001, February). Geostatistics: Modeling Spatial Uncertainty. Computers & Geosciences, 27(1), 121–123. https://doi.org/10.1016/s0098-3004(00)00063-7 https://doi.org/10.1007/s11004-018-9763-9