Integrating AI and Simulation for End-to-End Mine to Mill Optimisation: A Meta-Modelling Framework

1. Introduction

The mining industry comprises complex, interconnected operations that are crucial for delivering the raw materials used in a wide range of downstream products and industries. Every stage from exploration and extraction to processing and delivery has a direct impact on the efficiency, sustainability, and profitability of the whole value chain. However, traditional management approaches often rely on siloed systems and disconnected decision-making, which can result in inefficiencies and missed opportunities to improve performance. In response to these challenges, mine to mill optimisation has emerged as a holistic strategy that seeks to integrate mining and milling processes to maximise overall value. Valery et al. [1] discussed a comprehensive optimisation strategy for mining operations that integrates mine and mill processes to maximise profitability. They noted that the production, efficiency, and cost benefits of mine to mill optimisation are only going to become more critical in the future as the mining industry faces growing challenges with decreasing feed grades and increasing cost pressures. In 1998, Morrison [2] investigated the energy efficiency of different comminution circuits by implementing simulation techniques. The most efficient setups were identified in this study by comparing the performance of various circuit configurations. Scott et al. [3] focused their research on rock mass characteristics, blasting, excavation, crushing and milling performance, all of which are required to improve ‘mine to mill’ optimisation. The overall benefits are more likely to be observed as longer-term trends in performance rather than immediate responses to short-term changes. Using a combination of case studies, literature reviews, and empirical data, the study analysed the impacts of intense blasting on downstream comminution benefits. The results highlighted the potential economic benefits of optimising rock breakage processes, with substantial increases in mill throughput and reductions in milling energy consumption.

In 2008, Coetzee et al. [4] applied Model Predictive Control (MPC) to optimise a run-of-mine ore milling circuit. Their methodology was to design the MPC framework to achieve the objective of reducing the variability in the product particle size which leads to increased recovery of gold from downstream processes. Blast fragmentation models and their impacts on comminution circuits were reviewed by Esen et al. [5], with the emphasis on the importance of -10mm fines for mill throughput. The models were validated in their study with new data sets and applied in two mine to mill projects. Their case studies demonstrated that the throughput increased significantly due to improved fragmentation. Kanchibotla et al. [6] focused their study on optimising the drilling and blasting (D&B) process to enhance milling efficiency and determine the impact of blasting on the value chain. The study included baseline and validation trials, along with simulations to improve mill throughput and reduce energy consumption. The results showed an increase in SAG mill throughput and reductions in specific energy consumption. The impacts of optimised blasting on downstream processes shows that the increased throughput from baseline to validation was mainly due to reduced ore hardness, followed by finer ROM fragmentation from more intense blasting and tighter pebble crusher Closed Side Settings.

In 2015, Adebayo et al. [7] investigated the impact of blast-hole drilling deviation on muck-pile loading costs at the Hwange Colliery Mine in Zimbabwe. The excavation costs for different bucket capacities were analysed and the results showed a strong correlation between extra loading cost and blast-hole deviation, with significant variations in drilling and excavation costs. Xingwana [8] focused on tracking ore loss and dilution to enhance the integration of mining and milling processes in deep gold mines where ore movement through detailed tracking reports was monitored and survey data were analysed. The results provided an enhanced understanding of ore movement and helped to identify key factors that influence ore loss and dilution. In 2017, Beyglou et al. [9] focused their study on improving loading and crushing efficiencies at the Aitik Mine by optimising blast-induced fragmentation. In their study, field trials target specific fragmentation classes to find the most economical blast design. The results showed how to achieve higher loading and crushing efficiencies and more cost-effective milling by correlating target fragmentation with crusher and mill operations.

In 2018, Dotto et al. [10] analysed the effects of fragmentation size distributions on truck and shovel productivity. In their study, they noted that the fragmentation size distribution affects shovel production in terms of digability and fill factor. In addition, digability impacts the loading time and the energy used in loading and, subsequently, truck productivity is also affected by payload, loading time and waiting time at the loader. Esen et al. [11] conducted a drill to mill plant optimisation project at the Pustynnoye Gold Mine to increase mill throughput and improve overall mine performance. They used in-situ ore body characterisation. As a result, the project achieved a 16.2% increase in crusher throughput and an 11.6% increase in mill throughput, along with improvements in fragmentation, wall control, dig rates, and operational efficiency. Park et al. [12] used drilling performance to improve rock-breakage efficiencies. Based on the results of their experiments, models were developed to predict rock-breakage efficiencies during crushing and grinding. By using this technique, blast designs that include the comminution energy consumption and efficiency in downstream rock-breakage processes can be obtained. Rafeeian et al. [13] used the mine to mill process to identify and rank parameters that have an impact on blasting fragmentation. Statistical analysis and design of experiments (DOE) were used to help them determine the significance of parameters such as specific charge and hole diameter. The results highlighted the importance of specific charge in the crushing process, which was ranked the highest among the input parameters. Faramarzi et al. [14] implemented the Extended Drop Weight Test (ExDWT) to study the ore heterogeneity and its impacts on comminution circuits. They also used simulations to evaluate the impacts of ore competence variability on SAG mill performance. Their results showed significant throughput fluctuations due to ore variability and noted the importance of optimal sampling strategies to accurately represent ore domains. The Grade Engineering concept and the Integrated Extraction Simulator (IES) were implemented in Amini et al. [15] in which neural networks were employed to develop a fully integrated value chain model for economic evaluation and mine planning. The results showed significant insights into optimising the value chain and identifying low risk high-value scenarios. Ghahramanieisalou and Sattarvand [16] noted that the integration of Digital Twins (DT) in mining engineering can help optimise productivity and enhance efficiency. Their study emphasised that DT is an innovative technology that can be used for data storage, simulation and analysis, as well as a visualisation platform. The study discussed the benefits of using DT in several applications, including fleet management, mine to mill optimisation and geotechnical assessment. Additionally, Servin et al. [17] provided a digital material model for tracking ore movement through a stockpile to improve grinding efficiency in SAG mills. They used a pseudo-particle model, calibrated with real material test piles and physical tracer measurements. The implementation of the model improved material flow in the designed processes.

Hosseini et al. [18] integrated strategic and tactical mine planning to improve the overall value chain performance using a multi-objective grey wolf optimisation (MOGWO) approach and data-driven models. The application of MOGWO demonstrated significant efficiency and profitability improvements. Navarro et al. [19] obtained key parameters for the production process of structural canga (CE), compact jaspillite (JC), highly compact jaspillite (JAC) and extremely compact jaspillite (JEC) and presented an economic analysis of a mine to crusher model for an open-pit truck-less iron ore mine in northern Brazil, where operating costs for drilling, blasting, loading and primary crushing were modelled to correlate blasting design data to particle size distributions. Finally, in Hanhiniemi [20], the relationship between ball mill throughput and classification efficiency was examined in detail.

While traditional approaches have significantly improved various aspects of mining operations, recent research and industry practices point to the application of artificial intelligence as a powerful next step.

With AI-powered systems, mining companies can optimise their operations by analysing and integrating data from different sources. Table 1 summarises the applications of artificial intelligence (AI) and its subset, machine learning (ML), in various areas of mining operations.

In particular, optimisation at the integrated system-level of the entire mining process (from mining to comminution) is still not commonly applied in the industry. Past optimisations were often siloed (focusing on one stage or one process at a time) or lacked advanced AI techniques for analysing complex interdependencies. Based on our comprehensive literature review, the research gaps that were identified in this area include: (1) the lack of a systematic approach to mine to mill system optimisation that fully integrates all stages, (2) limited applications of modern AI and machine learning techniques for predictive modelling and optimisation, and (3) a disconnection between orebody characterisations and operational parameters, leading to reactive rather than proactive sub-optimal operational settings. To address these challenges, a systematic, data-driven approach that leverages advanced simulation tools and machine learning techniques is adopted in this study to address these gaps by developing an integrated AI-driven “meta-model” of the mining value chain. The meta-model approach uses a simulator to generate vast amounts of scenario-based data and then trains ML models on this data to capture key relationships and sensitivities. The meta-model can then be used to provide near real-time predictions and optimisation insights without having to always seek solutions from computationally expensive simulations. Consequently, the entire mining operation can remain in its optimal state and become more effective and efficient by following a “knowledge-based” mine to mill intelligent optimisation system. This study advances mine to mill optimisation in several keyways. While earlier mine to mill studies have typically focused on isolated stages, such as blasting or milling, the framework developed in this study integrates the entire chain from blasting to screening, crushing, stockpiling and SAG milling within a single simulation environment. In contrast to prior studies, which often rely on limited site trials or simplified models, the Integrated Extraction Simulator (IES) from Orica was employed to generate a high-fidelity digital representation of the Cortez Mine Complex flowsheet. In addition to this, meta-models trained on more than three million simulations enable robust prediction of downstream performance under a wide range of operational settings. SHapley Additive exPlanations (SHAP) analysis was applied to provide transparent and quantitative insights into the relative importance of parameters across multiple mining stages. Together, these innovations provide a step change from previous approaches by combining physics-based simulation, large-scale scenario generation, machine learning and interpretable AI into a unified optimisation framework.

2. Methodology

The methodology for this research involves an integrated, end-to-end approach that combines real operational data, high-fidelity simulation, and machine learning to develop a predictive meta-model of the mine to mill system. This meta-model is trained on millions of simulated scenarios and enables rapid, real-time predictions without the need for repeated full simulations. To ensure transparency and interpretability, SHapley Additive exPlanations (SHAP) analysis is applied to quantify the influence of input parameters and identify the key drivers of system performance. The key steps in this approach are: (1) data acquisition and case study definition, (2) building and validating a simulation model for the operation, (3) conducting large-scale scenario simulations (sensitivity analysis) to generate training data, (4) training and evaluating machine learning models (meta-models) based on the simulation data, and (5) interpreting the meta-models to identify important drivers and gain insights for operation optimisations. These steps are described in detail below.

1.1. Case Study and Data Acquisition

To ensure the practical significance of the research, the Cortez Gold Complex (specifically the Wenban pit at Barrick’s Cortez operation in Nevada, USA) was chosen as the case study. This site was selected because of the availability of detailed published data and reports on its geology, mine plan, and processing operations. In particular, the technical reports of Fiddes et al. [55] for the Cortez Mine and the PhD thesis of Hanhiniemi [56] were referenced extensively and include the process flow diagrams, production rates, equipment specifications, and rock characteristics. Key aspects of the Cortez/Wenban orebody and operations include gold ore with moderate hardness, drilling-and-blasting in the open-pit, feed to the primary crusher, followed by SAG (semi-autogenous grinding) milling and downstream recovery processes. The process flow diagram for the ore handling and grinding circuit was extracted from the report and used as the blueprint for our simulation model.

The process flow diagram of the Cortez gold mine begins with the blasting stage. The extracted rocks are then transported to a primary crusher, and the crushed materials are fed directly into a SAG mill, which operates in a closed loop with a SAG screen and a cone crusher. The next stage involves a sump where the output from the SAG mill is mixed with water. This slurry is subsequently pumped to hydro cyclones and a ball mill, for further grinding and size reduction.

2.1. Simulations

2.1.1. Simulation Background

Modelling of mineral processing plants was established by the JKMRC as one of the pioneer developers in the early 1970s [57]. Since then, models were regularly developed and improved for other processes in the mining value chain, such as flotation and blast fragmentation [58]. Process modelling methodologies have also been developed by other research and commercial entities, in the form of equations or specialised commercial software such as METSIM [59], MODSIM [60], SysCAD [61], CEET [62], FLEET [63], JKSimMet [64], and JKSimFloat [65].

2.1.2. Process Simulation with IES

Given the objectives of this study, a simulation platform capable of modelling the entire mine to mill chain and quickly evaluating a large number of scenarios is required. The Integrated Extraction Simulator (IES) is chosen as the core simulation tool, as it is a cutting-edge platform designed for end-to-end mining process simulations.

The IES allows the user to construct a customised flowsheet of unit process models. For this case study, the flowsheet in IES includes the following linked units: Blasting, Screening, Primary Crusher (Jaw), Stockpile and SAG mill. Each unit in IES is underpinned by a model that can be either physics-based or data-driven as described in the following paragraphs:

(1) Blasting module: The Kuz-Ram fragmentation model was parameterised using site-specific burden, spacing, hole diameter, explosive density, and velocity of detonation. Fragmentation curves were calibrated against Split-Desktop photo-analysis of muck piles to ensure realistic feed distributions.

(2) Screening stage: A single-deck screen split blast-induced materials into over-size and under-size material fractions according to the cut-size (D50) and efficiency parameter (Alpha). Screen parameters were gathered from the mining site based on real screen specifications.

(3) Primary crusher: A jaw crusher model received oversize materials from the screening stage. Its breakage function and throughput capacity were set based on equipment specifications reported in the Cortez technical report, ensuring the matching of throughput with operational ranges.

(4) Stockpile: The stockpile was simulated as a buffer and blending unit, capturing residence time effects and smoothing feed fluctuations. Material flow and particle size distributions in the stockpile were validated against reported site data.

(5) SAG mill circuit: The SAG mill received materials from the stockpile via a sump with added water. Mill behaviour was controlled by critical speed, load percentage, and water content. The mill operated in a closed circuit with a discharge screen and a pebble crusher, replicating the actual comminution loop. Outputs were compared against measured SAG mill product size distributions for calibration.

The fully connected flowsheet enabled the propagation of material characteristics from the blast face to the SAG mill, ensuring that upstream design changes directly influenced downstream performance. The model was iteratively calibrated at each stage until simulated size distributions and mass balances aligned with the available plant data. Figure 1 shows the schematic diagram of this flowsheet implemented in the IES platform.

Before proceeding to large-scale simulations, a calibration and validation of the implemented IES model against actual operating data was performed. This step is crucial to ensure the accuracy and credibility of the simulation model. Several adjustments were made iteratively: the output from the blast fragmentation model was compared to post-blast muck pile size measurements obtained via the image analysis of actual muck pile photos using the Split-desktop software [56]. Once validated, the IES model is used as a “safe sandbox” to conduct extensive numerical experiments that would be impractical to conduct in an actual operation.

2.1.3. Blasting Simulation

A comprehensive set of photographs were collected from post blasting muck piles, which were then analysed to derive the particle size distribution as mentioned above [56].

These data were then fed into the IES platform to enable the simulation of the blasting operation, resulting in a detailed and realistic model of the process. A comparison of predicted outcomes from simulation and real data is shown in Figure 2.

2.1.4. Stockpile Simulation

In the IES flowsheet, the stockpile product was estimated as a blend of two incoming streams: (i) the crusher product, and (ii) the undersize material from the screen. The simulator combined the particle size distributions (PSDs) of these two sources according to their relative mass flow rates, generating a homogenised PSD that represents the stockpile feed. This blending approach ensured that the stockpile did not act as a simple pass-through, but rather as a buffer that mixed materials from different time steps and smoothed short-term fluctuations in feed quality. As a result, the PSD of the reclaimed stockpile materials reflected both the fines passing through the screen and the coarser fractions retained and broken in the crusher. The simulated stockpile PSD was validated against site-reported measurements at the reclaim point. Figure 3 shows a comparison between actual and predicted particle size distributions within the stockpile. This curve highlights the effectiveness of the screening and crushing process in achieving the desired particle sizes. Throughout this paper, On-Screen (OS) denotes the screen oversize stream (material retained on the screen), Under-Screen (US) denotes the screen undersize stream (material passing through the apertures), and Crusher outputs refer to the crusher product stream discharged after comminution. These definitions are applied consistently across all figures and tables to avoid ambiguity.

2.1.5. SAG Mill Simulation

At the Cortez Gold Mine, the ore processing circuit included a Semi-Autogenous Grinding (SAG) mill as the key part of the comminution process. The SAG mill received materials from a sump, which acted as a buffer to maintain a steady feed rate. After grinding, the resulting slurry was sent to a screen for further processing. As real SAG mill data were not available, the SAG mill prediction was compared with the output data of the SAG mill screen (On-screen and Under-screen). The performance of this circuit is shown in Figure 4 and Figure 5.

3. Interpretation and Validation of Scenarios

Once the simulation models are validated, a sensitivity analysis is conducted by systematically varying the variables at each stage to assess their individual and collective impacts on the overall system performance. This involves running simulations of a vast number of scenarios to explore the effects of changes in parameters such as mineral characteristics, blast design, and grinding circuit variables. For our simulation platform, the flowchart of the modelled mining system showing all the inputs and outputs for each stage are illustrated in Figure 6.

A total of three million simulation scenarios were generated in IES to capture the variability and sensitivity of the mine to mill system. Parameter ranges were defined based on published technical data for the Cortez mine and typical operating windows reported in the literature. For example, blasting variables such as burden, spacing, hole diameter, explosive density, and velocity of detonation were varied within practical site limits. Rock properties, including unconfined compressive strength (UCS) and Young’s modulus, were also incorporated. Downstream parameters included screen cut-size (D50) and sharpness (Alpha), crusher close-side setting, SAG mill load, speed, and water content. Sampling was performed using a randomised grid approach that combined elements of full-factorial exploration with random sampling to ensure broad coverage of the parameter space without the prohibitive cost of a strict full factorial design. This approach generated a diverse and balanced dataset suitable for training machine learning models.

Before running large-scale simulations, the IES flowsheet was carefully calibrated and validated against actual operating data from the Cortez mine. Blasting fragmentation models were benchmarked using Split-Desktop image analysis of muck piles, crusher and screen outputs were compared with site-reported size distributions, and SAG mill discharge curves were cross-checked with published plant data. Calibration was iterative: model parameters were adjusted until the simulated particle size distributions and mass flows closely matched observed values. Only after this validation step was complete was the simulator used as a “digital twin” environment for extensive scenario generation. This ensured that the simulated outputs were realistic and credible representations of plant behaviour.

3.1. Algorithms

Machine learning (ML), a branch of artificial intelligence (AI), enables systems to learn patterns and make predictions from data without being explicitly programmed. Unlike traditional statistical methods, ML algorithms can handle complex, multi-dimensional problems and uncover hidden relationships between input (independent) and output (dependent) variables [66]. As mentioned earlier, in addition to the simple linear regression analysis, three ML-based algorithms, namely Random Forest (RF), Decision Tree (DT), and eXtreme Gradient Boosting (XGBoost) were used in this study to evaluate their performance in predicting the outputs of our mine to mill system. These algorithms are briefly explained below.

3.1.1. Decision Tree (DT) Algorithm

The decision tree algorithm is based on heuristics, such as the greedy approach, where the tree is constructed gradually, and locally optimal decisions are made at the construction of each node [68]. A decision tree is a tree-like model where internal nodes represent features of the dataset, branches indicate decision rules, and leaf nodes show the final outcomes or classifications. While leaf nodes signify the end results and do not branch further, decision nodes are used to split the data based on different criteria and can lead to multiple branches.

A decision tree is built by recursively splitting the dataset based on the best variable selected using a specific criterion such as variance reduction, mean squared error, mean absolute deviation, or information gain chosen according to the nature of the data. The data are divided into two subsets at each split, forming child nodes, and this process continues until a predefined stopping condition is met, such as reaching a maximum tree depth or a minimum number of samples. Once the tree is fully constructed, it undergoes pruning to remove nodes that do not significantly improve prediction performance, helping to reduce overfitting and enhance generalisation on unseen data.

3.1.2. Random Forest (RF) Algorithm

Random Forest (RF) is a popular supervised machine learning algorithm used for both classification and regression tasks. It operates by building multiple decision trees on different subsets of the dataset and then combining their predictions to produce the final result. This ensemble approach helps reduce the risk of overfitting and enhances the overall prediction accuracy of the model.

The algorithm builds multiple decision trees using a technique termed bootstrap aggregation, or bagging, where each tree is trained on a randomly sampled subset of the data. At each split within a tree, only a random subset of features is considered, which promotes diversity among the trees and reduces their correlation. The trees grow recursively until a specified stopping condition such as maximum depth or minimum samples per leaf is met. For prediction, a new observation is passed through all trees, and the final output is the average of their predictions in regression tasks. To achieve optimal performance, several hyperparameters must be tuned, including the number of trees, tree depth, minimum samples for splitting, minimum leaf size, and the number of features considered at each split [67,68].

3.1.3. Linear Regression

Linear regression is one of the most basic and widely used machine learning techniques for predictive analysis. It is a statistical method that models the relationship between a dependent (or output) variable and one or more independent (or input) variables using a linear approach. In LR, the relationships are modelled using linear predictor functions, with unknown model parameters estimated from the data [69,70,71]. The general structure of a multivariate linear regression is as follows:

$$Y=b_0+b_1{X}_1+b_2{X}_2+\dots +b_p{X}_p+e$$

(1)

where, $Y$ represents the dependent variable, and $X_1$, $X_2$, …, $X_p$ are the independent variables. The terms $b_0$, $b_1$, $b_2$, …, $b_p$ are the regression coefficients, and e stands for the error term. The model is typically trained using a method such as Ordinary Least Squares (OLS), which works by minimising the total squared difference between the predicted and actual values of the dependent variable. Each coefficient indicates how much $Y$ is expected to change with a one-unit increment in the corresponding $X$ value, assuming all other variables remain the same [72,73].

3.1.4. eXtreme Gradient Boosting (XGBoost) Algorithm

eXtreme Gradient Boosting (XGBoost) is a well-known and powerful machine learning (ensemble learning) algorithm that combines multiple weak decision trees to create a robust predictive model, which has been used for a variety of purposes in engineering [74,75,76]. It is a sequential learning technique where each new tree is built to correct the errors made by the previous one. Unlike the RF algorithm, which selects splits randomly, XGBoost focuses on improving prediction accuracy by adjusting to residual errors. At each iteration, the model updates weights, giving more importance to incorrectly predicted instances and less to correctly predicted ones. This process continues through multiple iterations, with the performance of the model better refined at each step [77,78].

The XGBoost method is based on the Gradient Boosting Decision Tree (GBDT) structure. However, in contrast to GBDT, the XGBoost objective function has a regularisation term that protects against overfitting.

Additionally, in contrast to Gradient Boosted Decision Trees (GBDT), which uses only the first-order derivative of the loss function, XGBoost incorporates both the first and the second-order derivatives by applying a second-order Taylor expansion of the objective function. Due to its built-in regularisation, XGBoost is more effective in preventing overfitting compared to other models.

The main hyperparameters used for the machine learning algorithms implemented in this study are shown in Table 3, and include Linear Regression, Decision Tree, Random Forest, and XGBoost. These parameters were selected based on model suitability for the specific prediction tasks within each mining stage. Key configuration options such as the number of estimators, tree depth, learning rate, regularisation terms, and data sampling strategies were tailored to balance model complexity and predictive accuracy. Highlighting these hyperparameters allows for reproducibility and provides insight into the relative flexibility and control offered by each method.

3.2. Model Validation

In this study, the Coefficient of Determination (R²) is used to assess the prediction accuracy of the proposed models. The mathematical expressions for the performance indicators are provided below:

$${R^2}_{(ytrue,ypred)}=1-{\displaystyle \frac{\sum {(\mathrm{y}\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}-\mathrm{y}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d})}^2}{\sum {(\mathrm{y}\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}-\mathrm{ȳ})}^2}}={\displaystyle \frac{RSS}{TSS}}$$

(8)

where $RSS$, $TSS$, $y_{true}$, and $y_{pred}$are, respectively, the sum of the squares of the residuals, the total sum of squares of the deviations from the average$\overline{y}$of true values, the measured/real output, and the predicted output.

Using this performance indicator, the accuracy of the developed machine learning models can be evaluated.

RMSE and MAE are other accuracy metrics that are used in this study to provide a robust comparison among different algorithms. The accuracy metrics of the different algorithms for all stages are shown in Table 4.

Figure 7 presents a comparative analysis of the prediction performance (R² values) of four machine learning algorithms - Linear Regression, Decision Tree, Random Forest, and Gradient Boosting - across multiple stages of the mining processes modelled. For each output metric (P20, P50 and P80), the R² values indicate the goodness-of-fit of the model. All ML models show excellent accuracy of more than 90 percent; however, Gradient Boosting consistently delivers the highest R² scores in the majority of cases, particularly in blasting and SAG mill predictions, demonstrating superior prediction accuracy. Random Forest also performs well, especially in the crushing and screening stages. In contrast, Decision Tree and Linear Regression show relatively poorer performance across most processes, with Linear Regression being the least accurate in modelling the complex, non-linear relationships in this dataset. Overall, the models demonstrated strong prediction performance, with accuracy exceeding 90% for outputs such as P20, P50, P80 and mass flow, indicating that the simulation framework and scenario design were both robust and reliable.

The performance of all machine learning models was evaluated using a combination of training/testing splits and k-fold cross-validation to ensure both robustness and generalisability. The dataset of simulated scenarios was randomly divided into 80% training and 20% testing subsets. Models were trained on the training set, while predictive performance was assessed on the held-out test set. In addition, a 5-fold cross-validation was applied during hyperparameter tuning to avoid overfitting and to confirm that performance was consistent across different partitions of the data. For all models, the average R² values for both training and testing datasets exceeded 90% percent, indicating strong agreement between predicted and simulated outputs across all stages, particularly for outputs with complex non-linear dependencies. This performance is attributed to their built-in regularisation, and resistance to overfitting, making them suitable for high-dimensional data.

4. Sensitivity Analysis and SHAP Discussion

In this study, the four ML models were selected for their ability to handle complex relationships between mining variables, their resilience to overfitting, and their demonstrated effectiveness in industrial process optimisations. To enhance model interpretability, SHapley Additive exPlanations (SHAP) analysis was applied. For clarification, it is necessary to explain briefly the interpretation of SHAP (SHapley Additive exPlanations) plots. Each plot ranks the model input variables (features) according to their relative influence on the predicted output. The order of the features on the vertical axis reflects their importance, with the most influential parameters at the top. Each point on the plot corresponds to one simulated scenario, with the colour scale (red to blue) indicating whether the feature value is high (red) or low (blue). The horizontal position of each point represents the SHAP value, which quantifies the impact of that feature on the model prediction: positive SHAP values shift the output upward (e.g., towards coarser PSDs), while negative values shift it downward (e.g., towards finer PSDs). The spread of values along the x-axis shows the variability of the impact of the corresponding feature, with wider spreads indicating stronger or more variable influence on the outcome.

In simple terms, the plots show not just which parameters are the most critical, but also how changes in those parameters (low vs. high values) influence the outputs (P20, P50, P80, or mass flow).

The SHAP summary plots (see Figure 8a to 8c) for Blasting P20, P50 and P80 illustrate the influences of different input parameters on the prediction of blasting fragmentation sizes. Among the most influential parameters, the blast hole diameter shows a strong impact on P20, meaning that larger blast holes produce finer materials. This is because larger holes hold more explosives, increasing the energy intensity per unit volume of rocks and resulting in more effective breakage near the blasthole. In addition, both velocity of detonation (VOD) and explosives density exhibit a strong negative correlation with P20. Physically, higher VOD explosives create faster and more intense shockwaves that propagate more effectively through the rock, enhancing fracture development and leading to finer fragmentation. Similarly, higher explosive densities (higher powder factor) provide more energy per unit volume of rocks, increasing the intensity of rock breakage and reducing the resulting fragment sizes. The influences of spacing and burden are also significant, but in contrast to the aforementioned parameters, larger values lead to higher P20. This is expected, as when the spacing or burden between blast holes increases, the energy from each explosive charge becomes more isolated and the local explosive charge density decreases, which reduces the overall interaction between shockwaves and results in coarser fragmentations. The unconfined compressive strength (UCS) of the rock has a direct physical impact on rock fragmentation. Stronger rocks with higher UCS require more energy to break and are therefore more resistant to fragmentation, typically resulting in larger fragments and a higher P20 value. Young’s Modulus, although less dominant, shows a subtle trend where rocks with higher stiffness tend to increase P20 slightly, likely because stiffer materials absorb and transmit energy more efficiently, reducing the intensity of localised fractures.

For the Blasting P50, it is influenced primarily by energy concentration and geometric confinement. The SHAP plot shows that larger hole diameters push the SHAP values strongly to the left, signalling a fall in P50 and therefore a finer median fragmentation size. Physically, the explosive energy inside a larger blasthole is spread over more rocks, increasing local pressure and expanding finer breakage. Higher explosive density drives the SHAP values leftward, lowering P50, which is because denser explosives pack more energy per unit length, producing stronger shock waves and finer fragments. Geometric controls follow this regime in that greater burden and spacing both add positive SHAP contributions because there are more unbroken rocks surrounding each charge, leading to decreased energy overlap hence larger fragments. Rock properties appear further down the list as stiffer rocks with higher Young’s Modulus absorb more of the shock energy and push P50 up slightly, while VOD and UCS appear to have only a modest leverage, indicating that once the 50 % passing threshold is reached, further energy increases, or rock strength variations have secondary effects compared with the charge geometry and explosive density.

For the coarser tail represented by Blasting P80, the same physical principles can be observed. The hole diameter variable again dominates, and wider holes exhibit stronger negative SHAP swing, confirming that fine fragments are produced when energy is increased inside the blasthole. Burden variable overtakes blasthole spacing as the second driver, because the distance to the free face primarily controls the blasting confinement and rock movement. A larger burden creates higher confinement and favours the formation of larger fragments, hence a higher P80. The explosive density is influential and its negative SHAP values show that high-density charges help to produce finer materials.

Figure 9a to 9c show the impact of blast induced materials on the On-Screen particle size distributions. Figure 9a clearly shows that Blasting P20 has the strongest influence on On-Screen P20 where the finest 20% of the blasted materials can contribute predominantly to the lower end of the On-Screen particle size distribution. Since Blasting P80 represents coarser fragments that almost always remain on the screen regardless of changes, its impact on On-Screen P20 is minimal and mostly constant. However, variables such as D50 and Blasting P50 still exert some influence. D50 determines the screen's cut-off size, so changes in D50 can slightly shift the boundary at which particles transition from passing through to being retained. Because Blasting P50 often aligns closely with this cut-off, its influence appears as a more gradual shift in the mid-range of the distribution, rather than a dominant control on the fine tail unlike Blasting P20.

The SHAP plots for On-Screen P50 and P80 clearly show that Blasting P80 is the most influential variable, and this aligns with the physical logic of screening. Since On-Screen P50 and On-Screen P80 represent the middle and upper tails of the On-Screen size distribution, they are most affected by the coarse fraction of the incoming feed. Blasting P80 defines the largest fragments generated during blasting, which are more likely to remain on the screen. As a result, higher values of Blasting P80 shift the retained size distribution toward coarser values, increasing both On-Screen P50 and On-Screen P80. In contrast, variables such as D50, Alpha, and Blasting P50 have only minor effects because they primarily influence finer to mid-sized materials, which does not strongly control the upper tail of the size distribution of the retained stream. Thus, the model correctly identifies Blasting P80 as the dominant factor shaping the coarser ends of the On-Screen size spectrum.

Figure 10a to 10c indicate a clear hierarchy of drivers for the Under-Screen size distribution. The screen aperture size (D50) exerts the greatest influence because it defines the physical aperture through which all materials must pass. Increasing D50 allows progressively coarser fragments to pass to the Under-Screen stream, thereby shifting Under-Screen P20, Under-Screen P50 and Under-Screen P80 to larger sizes, whereas reducing D50 has the opposite effect. Once the screen cut-off size (D50) is fixed, the availability of fine particles passing through the screen is primarily governed by Blasting P20. This is because particles smaller than the blasting P20 size make up the majority of the sub-D50 fraction. Therefore, higher Blasting P20 values result in a greater quantity of fines, which in turn raises the Under-Screen size percentiles. Blasting P50 impacts the Under-Screen size distribution when its corresponding particle size is close to the screen aperture size. In this case, an abundant mid-size materials increase the likelihood of borderline fragments passing the aperture especially when D50 is relatively large, therefore modestly raising Under-Screen P50 and, to a lesser extent, Under-Screen P80.

The SHAP plots for the crusher (see Figure 11a to 11c) show that the size of the oversize materials coming from the screen is the most critical parameter for evaluating the crusher output particle size. When the incoming median size (On-Screen P50) is large, the larger rocks within the feed tend to grind against each other more effectively, generating additional fines during crushing. As a result, the particle size corresponding to Crusher P20 becomes smaller, reflecting the increase in fine output. However, if the fine tail of the feed (On-Screen P20) is already coarse, the crusher receives fewer fines to begin with, and the Crusher P20 increases indicating a coarser product in the fine range. A high On-Screen P20 tends to make larger size materials entering the chamber, so the product P50 and P80 shifts to the right even considering the breakage. When the On-Screen P50 is large, the crusher gets a greater proportion of “medium-large” rocks that all must be crushed down to pass through the crusher. These rocks break down into smaller fragments, so the new midpoint of the product size distribution (Crusher P50 and P80) shifts downward and becomes finer. For On-Screen P80, when it increases, more oversized materials enter the crusher. These large size rocks take up space in the chamber and may take several compressions to break, so mid-size fragments receive fewer crushing impacts and pass through the crusher more easily.

For the SAG mill (Figure 12a to 12c), the SHAP analysis reveals that internal operational settings have significantly greater influences than feed characteristics on the SAG mill performance. For SAG mill P20, mill speed is the most important factor, followed by the amounts of load and water. When the mill rotates at a very high speed the rocks spend less interacting with each other, so the finest 20 % of the product becomes coarser and the value of P20 increases. A fuller mill or thinner slurry has the opposite effect and produces smaller P20 values, because they allow longer interaction time between materials inside the mill. For SAG mill P50, mill load (percentage) is the dominant factor. A higher charge increases breakage energy and produces finer particles (lower P50). Rotation speed has a dual effect: moderate speeds grind more efficiently (reducing P50), but excessive speed causes premature discharge, increasing P50. Water also contributes modestly increased water lowers P50 slightly by improving material flow and reducing viscosity, which enhances grinding. For SAG mill P80, mill speed dominates by a wide margin. Higher speeds lift coarse particles higher up the shell, leading to stronger impacts and more breakage. This reduces the coarse tail, lowering P80. The higher mill power increases the breakage energy, allowing rocks to grind to a finer size. More water lowers slurry viscosity, allowing particles and media to move more freely and collide with greater energy. This increases fines production and reduces P20, P50, and P80, provided the water addition doesn't shorten residence time excessively.

5. Analyses of Key Parameters Influencing Particle Size Distributions at Different Mining Stages

Understanding the influences of key input parameters on particle size distributions at different mining stages is crucial for the optimisation of operations from mine to SAG mill. The SHAP plots and analyses presented above provide valuable insights into the relative importance of different parameters, but they do not actually reveal the complex nonlinear relationships between input parameters and output size distributions. In this section, surface plots are used to visualize some of the key relationships. These plots help uncover intricate relationships and operational boundaries, offering a more comprehensive understanding of the combined impacts of certain parameters on the particle size distribution and throughput. This level of detail supports decision-making for parameter setting in operations and scenario optimisations.

5.1. Size Distributions of Material at Blasting Stage

Figure 13 illustrates the influence of explosive density on rock fragmentation across varying burdens and blast hole spacings. The three subplots respectively present the predicted P80, P50 and P20 values under explosive densities of 0.9 (green), 1.0 (red) and 1.1 (grey) g/cm³. Across all fragmentation metrics, a decrease in explosive density leads to a noticeable increase in fragment size because lower-density explosives generally produce lower detonation energy and gas pressure, reducing their ability to break the rock effectively. In addition, the plots demonstrate that both burden and spacing significantly affect fragmentation outcomes, with larger burdens and blast hole spacings generally producing coarser fragments due to the increase in gaps in energy overlap and the reduced interaction of stress waves between adjacent holes. The combined effects of burden, spacing, and explosive density underscores the importance of optimising explosive parameters to achieve desirable fragmentation for downstream processes.

The impacts of different burdens, spacings and UCSs on the blasting fragmentations are shown in Figure 14. However, based on the SHAP analysis in the previous section, the influence of UCS for P80 and P50 is negligible, but P20 is sensitive to UCS (30, 50, and 70 MPa). Higher UCS means the rock is harder to break and so the same explosive energy produces less fines and larger fragments. Other effects of burden and blast hole spacing on rock fragmentations are covered in the previous section.

As shown in Figure 15, the impacts of the velocity of detonation (VOD) on P20 is evident because a higher VOD means the explosive reaction occurs more rapidly, producing a stronger and more abrupt shock waves that initiates cracks more efficiently in the rock. This high-energy impact leads to more extensive micro-fracturing and better fragmentation, especially in the near-field zone around the blast hole. As a result, rocks are broken into smaller fragments. However, the only parameters that have effects on P50 and P80 are explosive density and blast hole diameter in the way that increasing the density of the charge and the blast hole diameter leads to a reduction in the size of the blasted materials. Hence a higher explosive charge density means more explosive per unit length of the hole, leading to stronger shock waves and higher gas pressure, hence enhanced rock fracturing. Similarly, a larger blast hole diameter allows for a greater total explosive charge (increased local explosive density), increasing the energy locally available to break the rock. This increased energy improves the intensity and the reach of crack formation and propagation, resulting in finer fragmentations.

5.2. Size Distributions of Materials after Screening

The size distributions of the On-Screen and Under-Screen outputs after screening are analysed in this section. Figure 16 illustrates the impact of key variables (based on SHAP analysis) on the On-Screen size distributions. As the physical meaning of different characteristics of blast-induced materials and their impact on screened product was discussed in the previous section, only the surface plots of the On-Screen and Under-Screen product sizes in relation to combined variables are shown here.

The impacts of key variables specifically D50, Blasting P20, and Blasting P50 on the Under-Screen particle size distributions, as revealed by SHAP analysis, are illustrated in Figure 17. Among these, Blasting P20 shows a strong and consistent influence, indicating that the fine fraction of the blasted material plays a critical role in shaping the Under-Screen product. In contrast, Blasting P50 exhibits minimal impact, suggesting that the median size of the blast feed does not significantly affect the screen's output distribution. The effect of D50 remains dominant, as it defines the screen's cut-off size and thus directly controls which particles pass through. Together, these results highlight the importance of fine generation during blasting and its downstream influence on screening performance.

5.3. Size Distributions of Crusher Products

The effects of the On-Screen particle size distribution on the crusher outputs are shown in Figure 18, where it can be seen that the impacts of On-Screen P50 are significantly more compared with the other two parameters.

5.4. Size Distributions of SAG Mill Outputs

Figure 19 shows the predicted particle size distribution outputs of P20, P50 and P80 as functions of mill critical speed and mill load at three different water content levels of 65%, 70% and 75%. It can be seen that increasing the mill critical speed generally results in coarser product sizes for P20 and P50 in contrast to P80. Mill load has a moderate influence, with higher loads tending to slightly decrease particle sizes. Water content notably affects the outcomes. Higher water content consistently leads to finer particle sizes, in which 70% and 75% water content produces lower particle size distributions than 65% water content. This effect is most evident in the P20 plot, where the differences are more pronounced, albeit less significantly, in the P80 plot.

Figure 20 illustrates the impact of particle size distributions of the stockpiles (On-Screen and Under-Screen products of the screening process) on the size distributions of the SAG mill outputs. It is clearly evident that, as the particle size of the stockpile materials increases, the particle size of the final product from the SAG mill also increases, indicating a direct relationship between the particle sizes of the inputs and SAG mill output products.

6. Conclusions

The meta-model developed in this study demonstrated exceptional predictive performance, achieving R² coefficient values exceeding 90% across all key process outputs, including fragmentation percentiles and mass flow. This high level of accuracy underscores the robustness and generalisability of the machine learning framework trained on millions of simulation scenarios. Importantly, the meta-model significantly reduces the need for repeated high-fidelity simulations, enabling near real-time prediction of system responses under various input conditions. This shift from simulation-based to prediction-based optimisation drastically lowers computational requirements while preserving decision quality, making the framework highly suitable for integration into live operational environments or digital twin systems.

In this study, an integrated, AI-driven meta-modelling framework to simulate and optimise the mine to mill processes from blasting to SAG mill in terms of particle size distributions was developed and validated. Leveraging real-world data from the Cortez Gold Mine and using the Integrated Extraction Simulator (IES) platform, over three million scenarios were simulated to identify key operational parameters. SHapley Additive exPlanations (SHAP) were applied to quantify the influence of each input parameter, offering transparent insights detailing the characteristics of the system. To further interpret model results, key complex relationships between input variables and the particle size distribution of process outputs were studied in significant detail. These relationships were analysed and discussed to guide process optimisation and deepen the understanding of the mining system dynamics. In the blasting stage, SHAP analysis confirmed that geometric and energetic factors dominate fragmentation outcomes. Hole diameter was consistently the most influential factor across all three output percentiles (P20, P50, P80), with larger diameters producing finer fragmentation due to higher explosive loading and energy intensity per unit rock volume. Explosive density and velocity of detonation (VOD) showed strong negative correlations with particle sizes, particularly P20 and P50, as denser and faster explosives generate more effective shock waves and enhance fracturing. Conversely, larger burden and blast hole spacing values were associated with coarser fragmentation. For rock mechanical properties, higher UCS values led to coarser fragments due to the increased resistance to breakage, while Young’s modulus played a lesser but consistent role in moderating fragmentation through energy absorption and transmission characteristics. For the screening process, the model highlighted that both feed characteristics and screen design parameters jointly determine the material separation. The aperture size of the screen (D50) was the most dominant driver of the Under-Screen particle size distribution. As D50 increased, progressively coarser particles were able to pass the screen, raising Under-Screen P20, Under-Screen P50, and Under-Screen P80 values. Blasting P20 followed closely, where the upstream with finer fragmentation led to more fines entering the Under-Screen stream, especially at P20. Blasting P50 had a noticeable influence when its value was near the screen aperture size, marginally increasing the fraction of mid-sized particles that could pass through. For the On-Screen stream, Blasting P80 was the key driver, particularly for On-Screen P50 and On-Screen P80. This makes physical sense, as coarse fragments are retained on the screen, and increases in P80 directly shift the retained size distribution to the right. In contrast, Blasting P20 had a strong influence primarily on On-Screen P20, as finer fragments from blasting were retained in the lower end of the On-Screen particle size distribution. Parameters such as screen sharpness (Alpha) and D50 provided minor tuning effects but were secondary to the feed characteristics. The crusher performance was governed by the size distribution of the On-Screen feed stream. Specifically, larger On-Screen P50 inputs led to increased breakage and finer product sizes (lower Crusher P50 and P80). However, coarser On-Screen P20 or On-Screen P80 inputs created more resistance to fragmentation, reducing efficiency and resulting in a coarser output, especially at P80. Internal operational parameters had the most significant impact in the SAG mill stage. Mill rotation speed dominated across all particle size percentiles, especially P80, by altering the dynamics of rock-lift and impact energy. The amounts of load and water content were also critical; higher load improved breakage by increasing charge mass and interaction frequency, while optimal water levels facilitated most effective slurry movement and grinding efficiency. Feed characteristics were comparatively less influential at this stage.

While this study focuses on the Cortez Mine, the framework is designed to be general and transferable to other mining operations. The core workflow combining calibrated IES simulations, large-scale scenario generation, and machine learning meta-modelling can be applied to different ores, equipment, and circuit configurations. To adapt the framework for a new site, the following steps would be required: (1). defining appropriate parameter ranges that reflect the geology and equipment of the operation; (2). calibrating the IES flowsheet with site-specific data (e.g., blast fragmentation, crusher performance, SAG mill throughput); and (3). regenerating a sufficiently large dataset of scenarios for training. In cases where ore hardness, equipment type, or circuit layout differ significantly, the meta-model would need to be retrained on these new simulation outputs. Thus, while the trained models in this study are site-specific, the methodology itself is broadly applicable across different mines and processing circuits, offering a systematic pathway to build digital twins and decision-support tools tailored to local conditions.

The proposed meta-model offers a practical pathway for supporting real-time decision-making across the mining value chain. By enabling near-instant predictions of downstream responses, the framework allows engineers to evaluate “what-if” scenarios without relying on time-consuming simulations. For example, blast design parameters such as burden, spacing, and hole diameter can be adjusted to deliver finer fragmentation that improves mill throughput, while operators can simultaneously assess crusher close-side settings or SAG mill speed and water addition to maintain grind size at minimum energy cost. This capacity to rapidly explore multiple operating alternatives ensures that decisions can be made proactively within daily or shift-based planning horizons, rather than retrospectively. As a result, operations are able to balance competing objectives such as fragmentation quality, throughput, energy efficiency, and cost with far greater agility than traditional trial-and-error approaches.

From a broader perspective, the integration of the meta-model into a control room or digital twin environment has the potential to deliver measurable sustainability and economic gains. Mine to mill studies have consistently shown that improved fragmentation and optimised milling conditions can increase throughput while lowering specific energy consumption, leading to reductions in both operating costs and greenhouse gas emissions. When linked to live data streams, the meta-model can serve as an intelligent decision-support system that connects orebody characteristics to final product outcomes, thereby maximising resource value while minimising waste. Although periodic calibration is required to account for ore variability and changing equipment performance, the framework provides a robust foundation for knowledge-driven, sustainable management of complex mining systems.

In future works, the meta-modelling framework will be extended to incorporate additional downstream processes, including the SAG mill discharge screens, hydro cyclones, and ball mills, providing a more comprehensive digital representation of the comminution circuit. Furthermore, detailed resource data, such as ore block characteristics and spatial heterogeneity, will be integrated into the modelling pipeline. This advance will enable the development of a fully connected digital twin of the entire mining system, linking orebody properties to final processing outcomes and facilitating more intelligent, resource-informed operational strategies.