Assessing Joint Simulation and Estimation Approaches in Geometallurgical Modeling

Abstract

Given that energy costs are a significant component of overall processing costs in mineral plants, reducing these costs through process optimization or technology adoption enhances the technical and financial feasibility of a deposit. Geometallurgical modeling plays a key role in understanding the relationship between material characteristics, mine planning, and processing stages, ultimately contributing to more efficient resource management and cost reduction in mineral processing. This study aims to develop a block model for evaluating comminution energy consumption (CEC) and identifying blocks with the highest energy usage potential during the grinding process in a specified region. Therefore, by applying advanced geostatistical techniques, including joint estimation and simulation based on geometallurgical data from multiple mineral processing stages, we predict CEC across the study area. The dataset encompasses 2.754 drill samples and a block model with 4.680 blocks. In this effort, imulation techniques, such as Plurigaussian and Turning Bands, provided more realistic outcomes than cokriging, considering the unique characteristics of geometallurgical data and the limitations of kriging methods.

Share and Cite:

Khorram, F. , Pais, G. and Reyes, M. (2025) Assessing Joint Simulation and Estimation Approaches in Geometallurgical Modeling. Journal of Geoscience and Environment Protection, 13, 316-330. doi: 10.4236/gep.2025.131017.

1. Introduction

At the core of geometallurgy lies a fundamental inquiry: How does a block behave during processing? This discipline has emerged as an indispensable tool for refining production planning, enabling the scrutiny of numerous pathways for mine development and commercial production strategies (Dominy et al., 2018). Traditionally, field modeling has predominantly focused on grade information, with the delineation of grade distribution enshrined within the reserve block model. However, crucial geological parameters such as mineralogy, lithology, alteration, and processing complexities like recovery often elude comprehensive consideration (Zhang & Wang, 2024). The meticulously crafted stock block model lays the groundwork for deeper mining endeavors, encompassing mining plans, processing methodologies, and exhaustive technical and economic analyses (Afum & Ben-Awuah, 2021).

In recent years, the necessity for a holistic approach spanning the mineral value chain has driven the evolution of geometallurgy. This nascent discipline intertwines exploration data with processing attributes, thereby mitigating the risks inherent in mining ventures (George et al., 2022). Geometallurgy’s comprehensive perspective significantly reduces operational hazards across extraction, processing, and final product dissemination. The profound influence of geological data on mineral processing is evident; for instance, the combination of unaltered fine-grained volcanic rock with soft sulfide minerals poses a significant crushing challenge that demands meticulous attention (Morales et al., 2019).

Geometallurgical methodologies typically fall into three categories: the conventional method, which uses metric grade values to compute mineral processing variables within a plant; the representative method, which involves processing tests on a variety of samples to directly gauge material processing response, such as the Bond Work Index (BWI); and the mineralogical method, which gathers limited information about the deposit’s mineralogy and its relationship with material processing attributes (Dominy et al., 2018).

In this study, we adopt the representative method, capitalizing on data from the Abrasion Index (AI) and BWI and their correlation with comminution energy consumption (CEC) to model this variable across diverse blocks (Dominy et al., 2018). The BWI quantifies a mineral’s resistance to crushing, encapsulating the work required per kilowatt-hour per ton to crush feeds and produce an output product. Innate reservoir parameters such as rock density, type, and porosity assume pivotal roles in estimating this geometallurgical parameter (Lamberg, 2011). The objective of this paper is to develop a specific energy consumption (SEC) block model and ascertain the maximum energy consumption anticipated from the grinding process through cokriging and joint simulation methodologies (Lishchuk et al., 2020; Lamberg, 2011). In recent years, a variety of advanced geostatistical methods has emerged to enhance the accuracy and reliability of spatial modeling in various fields, including mineral resource estimation and geometallurgy. By leveraging these advanced techniques, we aim to improve the prediction and management of energy consumption in the grinding process, contributing to more efficient and sustainable mining operations (Lishchuk et al., 2020).

2. Methodology

Cokriging is a geostatistical interpolation technique widely used in spatial modeling and prediction, particularly in fields such as geology, hydrology, and environmental science. This method extends the traditional kriging approach by incorporating auxiliary variables, which are correlated with the target variable, to improve prediction accuracy. Cokriging exploits the spatial relationship between the target variable and one or more secondary variables to provide more reliable estimates, especially in areas with limited data availability or uneven spatial distribution. By considering the cross-covariance between the target and auxiliary variables, cokriging can better capture spatial variability and reduce estimation errors. This technique has been applied in various domains, including mineral resource estimation, groundwater modeling, and soil mapping, to enhance the spatial prediction of important parameters (Isaaks & Srivastava, 1989; Chiles & Delfiner, 2012).

Opposed to multivariable estimation method, it can be adopted Plurigaussian simulation, pioneered by Journel and Alabert (1990), offers a powerful approach to capture the complex spatial heterogeneity present in geological formations by accounting for multiple facies or geological units. This method leverages conditional Gaussian simulation within each geological unit and subsequently combines these realizations to generate a composite model that reflects the spatial continuity and discontinuity across different facies. Plurigaussian simulation has been widely adopted in diverse applications, including groundwater modeling (Mariethoz et al., 2010) and petroleum reservoir characterization (Dubrule, 2003), showcasing its versatility and effectiveness in capturing complex geological structures.

Also, Turning Band Cosimulation stands out as a notable technique for integrating multiple sources of information and capturing the spatial complicated correlation between variables of interest. Initially proposed by Kitanidis (1995), this method facilitates the joint simulation of multiple correlated variables, enabling the incorporation of auxiliary data such as lithological information or geophysical data into the modeling process. By jointly simulating multiple variables, Turning Band Cosimulation provides a holistic representation of spatial variability, leading to more robust and accurate spatial models. This method has found widespread application in various fields, including hydrogeology (Christakos & Hristopulos, 1998) and environmental modeling (Kitanidis, 1995), underscoring its efficacy in addressing complex spatial problems through the integration of diverse datasets and variables. TBM delivers comparable accuracy to more computationally intensive multidimensional spectral methods, while being more economical. It also outperforms methods that rely on covariance matrix inversion, offering a more efficient approach. Prior to applying TBM, the distributions are assessed using the test of histograms and correlation clouds, which ensures that the two-variable distribution follows a Gaussian pattern and that its two-dimensional representation is elliptical. On the other hand, madogram and rhodogram test that is used in this paper, validates the relationship between madogram, rhodogram, and variogram (Lauzon & Marcotte, 2020).

The methodology of this paper adopts a comprehensive approach to addressing the diverse geological and metallurgical characteristics inherent in the study area, which includes seven distinct rock types and two RX types. To accurately assess these variations, the study necessitates the segregation of metallurgical fields and the evaluation of whether separate simulation and estimation are required for each. The predictive modeling process integrates all relevant variables to ensure the precise simulation of samples. Notably, Secondary Predictive Information (SPI) is not utilized, as it is mathematically correlated with CCE, acting merely as a secondary variable in its estimation without providing additional insights (Santamaría et al., 2022). Given the non-additive nature of geometallurgical variables like SPI and their unsuitability for linear methodologies such as cokriging at larger scales, the paper adopts conditional simulation as an alternative. This technique avoids the smoothing effects typically associated with kriging, thereby producing more authentic results characterized by individual fluctuations in each realization (Barnett & Deutsch, 2015).

The procedure follows several stages, starting with Exploratory Data Analysis (EDA) to gain insights through statistical tables and graphical representations categorized by rock type, lithology, and structural drivers. To address biases from irregular sampling, declustering is applied with cell sizes adjusted to the sample mesh. Cosimulation and co-estimation are central to predictive modeling, incorporating all geometallurgical variables. Key steps include transforming data into Gaussian form, applying anamorphosis modeling, validating multi-Gaussian hypotheses, conducting variographic analysis, simulating spatial patterns using spectral rotating bands, and performing cross-validation to assess model accuracy. Cosimulation is used to generate multiple scenarios, and realizations are processed for further analysis. Cross-validation ensures the model’s predictive accuracy and data integrity, enhancing reliability and robustness.

3. Variable of Interest

For this paper, we have two separate databases: one contains drill samples, and the other is a block model. The variables used are classified into categorical and continuous groups. Categorical variables include Rhyodacite, Andesite, Tuff, Sedimentary, Dioritic Dam, Rhyodacitic Dam, Others, in situ sulfide, waste, cavities, mine fills, and air. Continuous variables include the Semi-autogenous Rolling Power Index (SAG) (SPI, equivalent to SGI), as well as long-term and short-term SGI (associated with drillholes and blast holes, respectively). SPI represents the time required to grind a material to a P80 of 1700 μm in a semi-autogenous (SAG) grinding circuit. It has a direct relationship with the energy consumption (CCE) of the SAG grinding stage (Xu et al., 2018):

CCE[ kWh t ]=5.49 ( SPI T 80 ) 0.56 (1)

The Drop Weight test seeks to determine the energy required to fracture a mineral sample with specific characteristics, such as size and density. The rupture mechanisms of the test are measured through the impact rupture parameter (A×B). High values of this parameter indicate low specific energy consumption (Bueno et al., 2015). The Bond Work Index (BWI) expresses the resistance of a material to fracture and/or abrasion and measures the energy required to grind a material in a ball mill. It is measured in [kWh/t]. Crushing (grinding) consumes up to 4% of electrical energy globally, and approximately 50% of energy consumption at the mine site is used in crushing. The unit for BWI is Joules per kilogram and is related to the hardness of the samples; thus, the geological unit and rock type labels of each sample play an important role in its prediction (Xu et al., 2018; Bueno et al., 2015).

The blocked area for the 6.25 m × 6.25 m × 15 m mesh block model includes the following variables: rock density, as a continuous variable, which may have a correlation with crushing energy consumption. An indicator associated with the main failure is used as a categorical variable. Each part of this area can have its own separate stationarity, which can be checked using the variogram of the continuous variables in each section.

4. Discussion and Results

To address data issues in statistical analysis, the dataset was carefully processed by removing outliers through box plot analysis, imputing missing values, and merging duplicate entries. Records with undefined lithology or rock type codes were categorized as “other” to maintain data integrity. The dataset, marked by spatial heterogeneity, was divided into groups based on rock type, lithology, and structural factors, and statistical analyses, including histograms and multivariate analysis, were performed. Weak correlations between some variables led to the decision to treat them as auxiliary variables for comminution energy consumption (CCE) modeling.

Due to the non-additive nature of these variables, traditional kriging methods were deemed unreliable, prompting the use of simulation for more accurate modeling. The data revealed two spatial clusters, with the northern cluster showing higher lithological diversity and generally higher values for A×B, Bond Work Index (BWI), and CCE, reinforcing the importance of spatial variability. These results highlight the need for advanced simulation techniques to accurately capture the complexities in modeling CCE. Figure 1 illustrates a visualization of drilling and blast hole dat.

Figure 2 presents a three-dimensional visualization of the block model, highlighting the values of the variables under study. In the initial phase of the investigation, the focus was placed on the fundamental data, and as a result, only the spatial distribution of the block data was analyzed and visualized. This approach allowed for a clear representation of the data’s spatial arrangement, providing a foundational understanding before delving into more complex analyses.

Figure 1. Three-dimensional visualization of drilling and blast hole data and values of studied variables.

Figure 2. Three-dimensional visualization of block model and values of studied variables.

In the first part of the paper, we demonstrate that while variables do not change drastically at geological boundaries, the nature of failure due to structural events varies. Fault boundaries divide the area into two homogeneous sections, prompting the division of drilling data and blocks into two groups based on their proximity to the fault. For the block model, rock density is modeled to identify blocks linked to high comminution energy consumption (CCE), considering the correlation between CCE, rock density, and hardness. Due to the non-additive nature of these variables, simulation methods are used instead of kriging or cokriging for more accurate CCE modeling. Figure 3 demonstrates the variogram of SPI in different lithological domains.

The case study area is preprocessed for geostatistical estimation, ensuring compliance with stationarity assumptions. The block model includes relationship data between lithologies, which is used for contact analysis. Stationarity analysis is performed on lithological units, and variograms are calculated with specific parameters. Contact analysis between lithologies, especially in fault zones, is crucial for the final modeling process. The study consolidates lithological zones and rock types, treating fault zones separately due to their distinct properties, which could cause significant statistical changes.

Figure 3. Variogram of SPI in different lithological domains.

In contact analysis, distinguishing between hard and soft boundaries is essential. Hard boundaries feature clear discontinuities, while soft boundaries show continuity. The presence of edge effects is evaluated by analyzing variations in mean values near the boundary (Dumakor & Arya, 2021). To conduct this analysis, methods such as correlograms, mean value analysis, and delayed dispersion plots are used to explore spatial relationships. Separated areas greater than 50 meters are deemed ineffective for contact analysis, while smaller separations indicate smoother contacts (Maleki & Emery, 2020). Figure 4 shows the correlogram contact analysis

Figure 4. Correlogram contact analysis for CCE and density.

for CCE and density. However, an interface between certain lithologies showed less continuity and was excluded from analysis due to undefined classifications. The “others” geological unit was removed from further analysis to ensure consistency, and results were generalized to broader rock type ranges.

For variogram analysis, experimental variograms (both direct and crossed) are computed for the transformed Gaussian variables CCE, A × B, and BWI in both

Figure 5. Demonstration of directional variograms of variables, failure envelope.

Figure 6. Demonstration of directional variograms of variables, under fault.

horizontal and vertical directions. These are then fitted to a linear coregionalization model. Figure 5 illustrates the variogram models fitted to the data variograms above and below the fault, with each section modeled independently. Notably, CCE, BWI, and A × B exhibit negative spatial correlation in the upper section of the fault, while a slight positive correlation is observed in the lower section. Directional variogram of different variables are shown in Figure 5 and Figure 6.

Simulation and Post Processing

The variables CCE, A×B, and BWI are co-simulated within each block of the block model, conditioned on the available sampling data. The subsequent figures present maps of a realization (depicting selected sections) along with block-by-block statistical calculations, including mean, extremes, quartiles, and standard deviation, derived from the set of realizations.

Using the simulation methodology, as detailed in Table 1 and Table 2, fifty realizations were generated. A graphical representation of ten of these realizations

Figure 7. TBCOSIM block slices simulating variables.

Figure 8. Some implementations of TBCOSIM CCE simulation.

Table 1. Tbcosim parameters.

Search neighborhood (meter in each axis x,y,z)

800

800

800

Angles for search ellipsoid (degree in each axis x,y,z)

90

90

90

Number of realizations

50

Number of turning lines

500

Maximum number of locations to simulate simultaneously

5000

Seed for random number generation

9,784,498

Table 2. Plurisim parameters.

Find ellipsoid: radius in rotated system (meter in each axis x,y,z)

800

800

800

Angles that define the rotated system (degree in each axis x,y,z)

90

90

90

Number of realizations

50

Number of lines for spin band simulation

500

Number of target locations to process simultaneously, 30 iterations

5000

Seed for random number generation

9,784,498

is provided in Figure 7, Figure 8. As anticipated, energy consumption is found to be higher in the northern blocks of the study area. Figures 9-11 presents the statistical characteristics of 100 realizations for each simulated block, further corroborating that the northern blocks exhibit a greater potential for reducing energy consumption.

The primary objective of this study is to identify the blocks and regions with the maximum energy consumption. Given that CCE is expressed in kWh/ton, and considering the varying densities of different lithologies, blocks are categorized accordingly. Andesite and rhyodacite, sharing similar densities, are grouped together and marked with a reddish color. Tuff and sedimentary units, which

Figure 9. Some realizations of PLURICOSIM Simulation Lithologies.

Figure 10. Some realizations of PLURICOSIM Simulation Lithologies.

Figure 11. Some realizations of TBCOSIN simulation Consuming energy.

exhibit different densities, form a second group, represented in brown. The remaining lithologies, with distinct density characteristics, are classified in a third group and marked in green.

To simulate these lithologies, the Gauss-Pleuri method (detailed in Table 2) was applied within the same blocks, attributing CCE and the respective lithology densities to each block. Energy consumption is then calculated based on the volume of each block (40 m × 40 m × 40 m).

As it could be seen in Figure 12 and Figure 13, the cross-validation of the simulations for A × B, BWI, and CCE is illustrated through uncertainty and correlation plots. The uncertainty plot (accuracy plot) evaluates whether the simulations accurately reflect the uncertainty in values, with the plot’s alignment

Figure 12. Simulation cross-validations, failure envelope.

Figure 13. Simulation cross-validations, low failure.

along the diagonal indicating correct reproduction. The correlation plot compares real data values to the average of the simulated values, verifying that the simulations predict the variable of interest without significant global or conditional bias.

5. Conclusion

The primary aim of geometallurgy is to optimize the value of extracted minerals by integrating geological, mining, metallurgical, environmental, and social factors. Through SPI testing, standard size reduction is assessed and expressed as a hardness index in the SAG mill, while geometallurgical tests, such as the Band Work Index, enable the evaluation of commercial operations’ efficiency with respect to energy consumption during rock breakage or crushing. The determination of energy consumption is critical for optimizing the grinding and milling stages, and accurately predicting this parameter enhances equipment efficiency. One of the significant challenges in predicting geometallurgical variables lies in the inability to scale certain parameters effectively. In this context, this paper proposes the use of multi-Gaussian conditional simulation as a methodological approach to develop a block model of specific energy consumption. This approach provides a robust framework for addressing these challenges and offers a practical solution for improving energy efficiency in mining operations.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1] Afum, B. O., & Ben-Awuah, E. (2021). A Review of Models and Algorithms for Surface-Underground Mining Options and Transitions Optimization: Some Lessons Learnt and the Way Forward. Mining, 1, 112-134. [Google Scholar] [CrossRef]
[2] Barnett, R. M., & Deutsch, C. V. (2015). Multivariate Imputation of Unequally Sampled Geological Variables. Mathematical Geosciences, 47, 791-817. [Google Scholar] [CrossRef]
[3] Bueno, M., Foggiatto, B., & Lane, G. (2015). Geometallurgy Applied in Comminution to Minimize Design Risks. In 6th International Conference on Semi-Autogenous Grinding Technology (Vol. 11, pp. 1-19). High Press.
[4] Chilès, J., & Delfiner, P. (2012). Geostatistics: Modeling Spatial Uncertainty. Wiley. [Google Scholar] [CrossRef]
[5] Christakos, G., & Hristopulos, D. T. (1998). Modern Spatiotemporal Geostatistics. Oxford University Press.
[6] Dominy, S. C., O’Connor, L., Parbhakar-Fox, A., Glass, H. J., & Purevgerel, S. (2018). Geometallurgy—A Route to More Resilient Mine Operations. Minerals, 8, Article No. 560. [Google Scholar] [CrossRef]
[7] Dubrule, O. (2003). Geostatistics for Seismic Data Integration in Earth Models (pp. 1-30). Geostatistics Banff.
[8] Dumakor-Dupey, N. K., & Arya, S. (2021). Machine Learning—A Review of Applications in Mineral Resource Estimation. Energies, 14, Article No. 4079. [Google Scholar] [CrossRef]
[9] George, J., Whitten, D., Metters, R., & Abbey, J. (2022). Emancipatory Technology and Developing‐World Supply Chains: A Case Study of African Women Gemstone Miners. Journal of Operations Management, 68, 619-648. [Google Scholar] [CrossRef]
[10] Isaaks, E. H., & Srivastava, M. (1989). Applied Geostatistics. Oxford University Press.
[11] Journel, A. G., & Alabert, F. G. (1990). New Method for Reservoir Mapping. Journal of Petroleum Technology, 42, 212-218. [Google Scholar] [CrossRef]
[12] Kitanidis, P. K. (1995). Introduction to Geostatistics: Applications in Hydrogeology. Cam-bridge University Press.
[13] Lamberg, P. (2011). Particles—The Bridge between Geology and Metallurgy. In Konferens i Mineralteknik (p. 241). Luleå Tekniska Universitet.
[14] Lauzon, D., & Marcotte, D. (2020). The Sequential Spectral Turning Band Simulator as an Alternative to Gibbs Sampler in Large Truncated-or Pluri-Gaussian Simulations. Stochastic Environmental Research and Risk Assessment, 34, 1939-1951. [Google Scholar] [CrossRef]
[15] Lishchuk, V., Koch, P., Ghorbani, Y., & Butcher, A. R. (2020). Towards Integrated Geometallurgical Approach: Critical Review of Current Practices and Future Trends. Minerals Engineering, 145, Article ID: 106072. [Google Scholar] [CrossRef]
[16] Maleki, M., & Emery, X. (2020). Geostatistics in the Presence of Geological Boundaries: Exploratory Tools for Contact Analysis. Ore Geology Reviews, 120, Article ID: 103397. [Google Scholar] [CrossRef]
[17] Mariethoz, G., Renard, P., & Straubhaar, J. (2010). The Direct Sampling Method to Perform Multiple‐Point Geostatistical Simulations. Water Resources Research, 46, W11536. [Google Scholar] [CrossRef]
[18] Morales, N., Seguel, S., Cáceres, A., Jélvez, E., & Alarcón, M. (2019). Incorporation of Geometallurgical Attributes and Geological Uncertainty into Long-Term Open-Pit Mine Planning. Minerals, 9, Article No. 108. [Google Scholar] [CrossRef]
[19] Santamaría-Bonfil, G., Santoyo, E., Díaz-González, L., & Arroyo-Figueroa, G. (2022). Equivalent Imputation Methodology for Handling Missing Data in Compositional Geochemical Databases of Geothermal Fluids. Geothermics, 104, Article ID: 102440. [Google Scholar] [CrossRef]
[20] Xu, L., Luo, K., & Zhao, Y. (2018). Numerical Prediction of Wear in SAG Mills Based on DEM Simulations. Powder Technology, 329, 353-363. [Google Scholar] [CrossRef]
[21] Zhang, Q., & Wang, T. (2024). Deep Learning for Exploring Landslides with Remote Sensing and Geo-Environmental Data: Frameworks, Progress, Challenges, and Opportunities. Remote Sensing, 16, Article No. 1344. [Google Scholar] [CrossRef]

Copyright © 2026 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.