Digital twins with distributed particle simulation for mine-to-mill material tracking

Systems for transport and processing of granular media are challenging to analyse, operate and optimise. In the mining and mineral processing industries these systems are chains of processes with complex interplay between the equipment, control, and the processed material. The material properties have natural variations that are usually only known at certain locations. Therefore, we explore a material-oriented approach to digital twins with a particle representation of the granular media. In digital form, the material is treated as pseudo-particles, each representing a large collection of real particles of various sizes, shapes and, mineral properties. Movements and changes in the state of the material are determined by the combined data from control systems, sensors, vehicle telematics, and simulation models at locations where no real sensors can see. The particle-based representation enables material tracking along the chain of processes. Each digital particle can act as a carrier of observational data generated by the equipment as it interacts with the real material. This makes it possible to better learn material properties from process observations, and to predict the effect on downstream processes. We test the technique on a mining simulator and demonstrate analysis that can be performed using data from cross-system material tracking.


Introduction
It is estimated that mining accounts for 6% of the global energy consumption. Grinding is the most energy demanding process, accounting for 32% of the mine's energy consumption, followed by haulage (24%), underground mine ventilation (9%), and digging (8%) [Holmberg et al.(2017)Holmberg, Kivikytö-Reponen, Härkisaari, Valtonen, and Erdemir]. The grinding efficiency is many times as low as 1% and depends strongly on the particle size distribution and ore hardness [Radziszewski(2013)]. Therefore there is a strong motivation for determining and controlling the properties of the ore fed into the mill. Mine-to-mill refers to a holistic approach to mine optimisation, see [McKee(2013)] for a comprehensive description. Instead of optimising each process individually, one takes the ore body characteristic into account, and consider the throughput and energy consumption of the full chain of processes when deciding on target fragmentation and controlling the feed to the mill. The mine-to-mill approach has, however, proven difficult to implement in practice. One reason is that it requires an elaborate digital infrastructure for collecting, processing and communicating data between edge devices, distributed control systems and centralized data services. The technology, popularized as digital twins [Grieves(2014)] or Industrial Internet of Things (IIoT) [Jeschke et al.(2017)Jeschke, Figure 1: Illustration of a mine-to-mill process with particle-based material tracking. The chain of processes include: drilling; blasting; load, haul and dump; conveying; crushing; sizing; storage and reclaim; and grinding; which then lead to a concentrator. The processes are carried out by physical assets that perform unit operations on the material. Connected virtual assets mirror the operations by transforming a distributed digital representation of the material. The architecture of the system of assets and services for tracking, data storage and analytics is illustrated in the lower right.
A mine-to-mill process chain is illustrated in Fig. 1. The performance of each individual process depends on its design, control, and on the properties of the material feed. Each process also alters the state of the material. The complex interdependency between the process dynamics and variations in the material properties makes it hard to determine what is cause and effect. The mineral concentration, hardness, strength, and toughness, of the ore body is partially known in advance, from exploratory drilling, and stored in a block model [Rossi and Deutsch(2014)]. These are input parameters to many models for blasting [Ouchterlony and Sanchidrián(2019)], crushing [Evertsson(2000)] or grinding [Napier-Munn et al.(1996)Napier-Munn, Morrell, Morrison, and Kojovic], such as the Kuz-Ram blasting and JKMRC grinding models. Measurement-while-drilling (MWD) of blast holes produce data that contain information about the geo-mechanical properties at higher spatial resolution than the exploration data. This may be used for improving the blast plans and for optimizing the subsequent processes [Rai et al.(2016) Rai, Schunnesson, Lindqvist, and Kumar, Zhou et al.(2011)Zhou, Hatherly, Ramos, and Nettleton]. Drilling and blasting is conducted in cycles, each cycle producing a volume of fragmented rock that is loaded using bucket excavators or load-haul-dump trucks before the next blasting. The particle size and shape distribution is subsequently transformed by the crushers and grinders before entering the concentrator. Variations in the diggability of the blasted rock [Singh and Narendrula(2006), Khorzoughi and Hall(2016), Brunton et al.(2003)Brunton, Thornton, Hodson, and Sprott] is useful for improving blast plans or for predicting the throughput at the crusher and mill. Transport and storage systems alters the location of the material, but may also induce (unwanted) mixing and segregation that propagate to subsequent processes and affect the dynamics there. Mine-to-mill analysis and optimization require knowledge of these changes but on-line measurements of granular media properties are inherently difficult. The material movements during vehicle and belt conveyor transport are easily tracked and these systems may be instrumented with sensors for monitoring some of the material properties.
But in chutes, stockpiles and silos, little is known about the material beyond the knowledge of what is fed into them. If any sensors can be installed at all, the observations are either indirect or limited to the surface of the granular media, which constitutes a small and perhaps not representative fraction of the bulk. Conventional tracking systems, that rely on combinations of plug-flow models (or first-in/first-out) and perfect mixing models that are calibrated using tracer experiments [Kvarnström and Bergquist(2012)], fall short on transport and storage systems with complex surface and discharge flow. This creates blind spots in the tracking of material movements and knowing the properties at any given location.
To remedy this, we explore digital twins with distributed particle simulation, of systems doing transport and processing of granular media. The idea is to represent the material and its movements in a structured data format based on a particle representation. In addition to the material's identity and current position, the particle data structure supports the reading and writing of observations from sensors and equipment along the chain of processes. When connected equipment performs unit operations on the material, the digital copy is updated accordingly. When the material reaches blind spots in the system, for example a silo or a stockpile, the digital copy is driven by a simulation model fed with data from the control system and available sensors in real-time. For the sake of computational performance and memory, the digital twin uses pseudo-particles that represents a large collection of real mineral particles with a distribution of grain size and shape. For fast simulation we use the nonsmooth discrete element method, which allow large time-step integration [Servin et al.(2014)Servin, Wang, Lacoursiére, and Bodin]. To achieve real-time performance in large processes, a data-driven model order reduction technique [Wallin and Servin(2021)] is employed. The combination of the two models in an executable application, that can be part of a distributed system like a digital twin, we call a granular surrogate.
We test the feasibility of doing material tracking with this technique on an idealized mine simulator with both discrete and continuous transport processes from two different sources. The tests include analysis of how crusher power draw and milling performance depends on variations in ore hardness and fragmentation at the blast location. From crusher time-series it is possible to identify the relation between the ore hardness and measurement-while-drilling signal despite having multiple ore sources being mixed. The sensitivity to tracking errors is analysed by perturbing the pseudo-particle origins. The response in mill performance to a change in ore hardness at the source is tested with different state of an intermediate stockpile. The effect is significant but it is still possible to reconcile the ore's grindability, in terms of the Bond Ball Mill index, to the locations of origin in the ore body thanks to the particle-based representation.
Related work include the data-driven mine-to-mill framework proposed and tested in [Erkayaoglu(2015)] and [Erkayaoglu and Dessureault(2019)]. It was concluded that material tracking needs to be performed in a more comprehensive way to better connect mining and mineral processing data, and stockpile management, and the required infrastructure was found to be the main limitation for implementing the framework. A closed-loop framework for real-time reserve management was introduced in [Benndorf and Buxton(2016)]. It incorporate sensor-based material characterisation, geostatistical modeling under uncertainty, data assimilation for a sequential model updating, and mining system simulation and optimisation. Using a modified Kalman filter technique, online sensor and process information is back-propagated into the resource block model, which in turn is used for production planning and other operative decisions. Innes et al. [Innes et al.(2011)Innes, Nettleton, andMelkumyan, Innes(2015)] developed a method for representing, tracking and fusing information of excavated material. The bulk material is represented as lumped masses, although not as particles and no other properties than mass and location. Tracking along the transport and storage process is achieved using an Augmented State Kalman Filter with a mass preservation constraint. Reconciliation of the Bond Mill Work index to the ore block model using mine-to-mill tracking was demonstrated in [Wambeke et al.(2018) Wambeke, Elder, Miller, Benndorf, and Peattie], using a first-in-first-out model for the stockpile. A 3D real-time stockpile model and mapping technique of product quality was developed in [S.(2016)] for the purpose of planning and control in stacking, reclaiming and blending. The pile surface was reconstructed from the point-cloud of a mobile scanning device. A voxelized model is then associated with a reclaiming machine to achieve autonomous operation and predictive quality. The simulator in [Berton et al.(2013)Berton, Jubinville, Hodouin, Prévost, and Navarra] demonstrate the combination of multiple, relatively simple, models for capturing the ore flow through a complex storage facility. No previous work was found where material tracking is made using a particle representation.
2 Digital twin as a distributed particle simulation From the material point of view, the digital twin of a mine may be considered a system of distributed particle simulation. For computational efficiency, we use pseudo-particles, each representing a large collection of real rock particles of various size, shape and mineral composition. Labelled n, a pseudo-particle have position xn(t) ∈ R 3 , velocity vn(t) =ẋn, mass mn ∈ R 1 , size distribution function fn ∈ R N f , mineral concentration cn ∈ R Nc , and mechanical properties hn ∈ R N h that characterize the real particles it represents, each discretized in classes of dimensionality N f , Nc, and N h , respectively. The particle state vector is represented Xn = [mn, xn, vn, fn, cn, hn]. Each particle has a position of origin, xn(t0). The material properties may be known to some degree of certainty at that location from prior exploration and represented in a block model of grid cells indexed i, with centre point xi and volume Vi. Or, they are the subject for determination by measurements in the process downstream. The fragmentation is set by the blasting model and/or measurements after blasting.
The fixed and mobile mining equipment constitute a discrete set of assets that transform the material by a set of unit operations, the primary two being move and fragment. Each pseudoparticle is tracked, or simulated, by a virtual asset which is a realtime model that mirrors a physical asset and its unit operations on the material. We distinguish between high-frequency micro-tacking inside each asset and macro-tracking of special events, e.g., when particles enter or exit an asset, or pass a point of special interest. This is illustrated in Fig. 2.  Figure 2: Illustration of micro and macro-tracking. The sparse macro-tracking events include the appearance of fragmented material at the source (pink or blue) after blasting, haul truck asset dispatching and dumping (yellow), material entering later being discharged from a storage asset (green) and a conveyor asset (grey). The micro-tracking is the relatively high frequency tracking of the material moving with or flowing inside an asset.
For a bucket excavator or truck asset, the micro-tracking can be accomplished simply by inheriting the position of the bucket or the truck during the transport from the loading to the dumping area. Particles transported on a belt conveyor may similarly be tracked between the loading and discharge points by simply tracking the movement of the belt using the conveyor control system. In the transporter's virtual asset, the belt is discretized longitudinally in bins of known location and particle content. Chutes, silos and stockpiles are different. In these assets the particle flow is passive, driven by gravity and feeders at the outlet for a controlled discharge rate. Micro-tracking require a model, that may range from simple first-in-first-out or perfect mixing models, to a physics-based flow model.
Mathematically, we represent a unit operation by an asset a as an operatorûa(t, ∆t) that propagate a particle state Xn(t) forward in time into the new state Xn(t+∆t) =ûa(t, ∆t)Xn(t), with propagation time-step ∆t. If the particle n enter the asset at time tin in a state Xn(tin), the exit state at time tout is given by where J a n is the number of unit operations by the asset. Micro-tracking is to resolve the evolution of the particles this way. In general, the step size ∆t may vary and the unit operation is a function of both time and the state Xa of the particles in the asset , i.e.,ûa(t, Xa(t)). Note that this representation support both kinematic and dynamic models for the particle motion, for move operations. Discrete element simulations, taking realtime asset data as input, is an example of the latter. Letting the particles co-move with a vehicle on conveyor is an example of the former. The evolution of the particle size distribution may also be represented by this, for fragment operations.
Macro-tracking is represented as message passing, between virtual assets or a centralized tracking service, about events of special interest and the involved particles. The message should have a structured data format, like including an identifier for the event, its type, location and time, plus a set of particle identifiers and state tuples. The macro-tracking data may then be stored in a database for later analytics or other services based on tracking data. An illustration of the architecture of the system of assets and services is found in Fig. 1.
Most assets are equipped with sensors that produce observations while handling the material. An observation is either a direct measurement of an intrinsic attribute or an indirect measurement. The latter may be registered as an extrinsic attribute that is assigned to the pseudo-particles at the time and location for the observation. Examples of extrinsic attributes are measurement-while-drilling (MWD) data, diggability, crushability, and grindability. In complement to real sensors, it is also possible to produce virtual sensors that enable monitoring or controlling a process at higher precision. A virtual measurement signal of an attribute p is accomplished by extracting the pseudo-particles N (x, t) in a selected volume V (x, t), centered around x, at time t, and computing the mass-weighted attribute with a moving average . . . tw of time window tw for filtering the effects of coarse particles. Each pseudo-particle has an individual size distribution, (3) to f (x, t) and the mean particle diameter by d (x, t) = d · f , where d is the diameters of the of particle bins. The cumulative size distribution F is the mass fraction of particles equal or larger than each size class d k , such that F k = Ns k =k f k . The DF value is the diameter d for which the cumulative size distribution equals F . It is easily determined from the cumulative size distribution F by direct inspection or interpolation of the nearest discrete size classes.
When particle simulation is needed for micro-tracking we use the discrete element method (DEM) with spherical pseudo-particles. The simulation attributes include particle diameter dp, coefficient of restitution ep, friction coefficient µp-t, rolling resistance coefficient µp-r, and particle cohesion cp. The parameters are calibrated to values that approximate the mechanics and flow dynamics of the real material on bulk level, as measured by the bulk mass density ρ b , internal friction µ b , cohesion c b , and angle of repose. The particle diameter is chosen as large as possible, but small enough for resolving the important flow features. Even with pseudo-particles as coarse as a loader bucket, the number of particles in a silo or stockpile asset may exceed what can be simulated in realtime with present hardware. To guarantee real-time performance of the particle dynamics, we employ data-driven model order reduction [Wallin and Servin(2021)]. The idea is to run numerous DEM simulations in advance, covering the relevant state-space as well as possible using authentic CAD models and control signals for in-and outflow to the asset. From the simulation data, a low-dimensional model is trained to predict the flow field that occur in a stockpile or silo as a function of the control signal and mass distribution. The computationally intense process of computing all contact forces is substituted by a simple evaluation of the predicted velocity field at the positions of the particles, which are simply advected using this velocity. This is applied to the particles in the bulk. To accurately capture the dispersion of particles impacting the pile surfaces, resolved DEM is used. The combination of the two models in an executable application, part of a distributed system, is referred to as granular surrogates.

Integration test
An integration test of using granular surrogates for micro-and macro-tracking was performed, see Fig. 3. A virtual stockpile asset was created using AGX Dynamics for the particle simulation. To accelerate the simulations to realtime performance, a reduced order model was trained from groundtruth simulation data, as described in [Wallin and Servin(2021)], with tracking accuracy of 10%. The storage asset was connected to an ABB 800xA process simulator [ABB()] using OPC [OPC()] for receiving a discharge control signal and for exchanging macro-tracking data. The ABB simulator represent the assets of the incoming and outgoing conveyors, discretized in bins moving with the controlled conveyor belt speed. Each bin on the incoming conveyor contain material with a specific mass, ore grade and a tracking-id. That information is exchanged with the stockpile asset when a bin reach the discharge point. The storage asset responds by emitting particles with the corresponding total mass and mineral concentration. The particles fall, impact, and disperse over the pile surface as simulated by the physics-based stockpile model. The storage asset receives a discharge rate control signal from the ABB simulator. This drives a flow field in the bulk and particles are propagated in realtime by the data-driven reduced order model. The particles that reach the outlet are eliminated and a macro-tracking message is passed so that the ABB simulator receive information of the material that enters the outgoing conveyer. That information include tracking-ids, mass and ore grade. The micro-tracking data resides in the storage asset. A monitoring service is created as a stand-alone web application. It receives virtual sensor information from the stockpile asset using a REST 1 interface. The monitoring information include the total mass in the storage, the average grade, and the spatial distributions of these quantities as well as the predicted velocity field in the bulk of the pile. A video that demonstrate the integration test is included as supplementary material S2. It should be pointed out that the ABB 800xA simulator can be switched to sharing realtime data from a mine with the same layout. Figure 3: The integration test system include two conveyor assets in an ABB 800xA simulator (left), a storage asset (centre) driven by a granular surrogate, and a monitoring service (right) displaying macro-tracking data and virtual sensor data from the storage. The communication interfaces rely on OPC and REST, respectively. In the granular surrogate, black particles are dynamic while the others are propagated using the reduced model and color coded by residence time. A video is included as supplementary material S1.

Simulator
Before implementation in a real mine, particle-based material tracking is tested in a simulator of a simplified 2D model of a mine shown in Fig. 4. A video of the simulator is included as supplementary material S2. Each asset can easily be replaced with more elaborate 3D models or virtual assets mirroring real physical assets. There are two sources of ore from which haul trucks transport blasted material to a single crusher. The crushed material is fed to a stockpile for intermediate storage, from which material is drawn and fed to a mill. Belt conveyors transport the material between the crusher, storage, and the mill. The simulator is used for testing the analysis that is enabled by the particle-based tracking and for understanding the sensitivity to tracking errors. Figure 4: Image from the simulator with material hauled with trucks from two sources, dumped at a jaw crusher, conveyed to and from a storage, and finally fed to a grinder. Green particles come from the left source and blue particles from the right source. A video is available as supplementary material S2.
Macro-tracking events include the haul trucks dispatching from the loading point, dumping at the crusher station, and material passing the crusher zone, the inlets and outlets to the stockpile, and the inlet to the grinder. Model order reduction is applied on the belt conveyors. The simulator is implemented using the physics engine AGX Dynamics [AGX()], using a nonsmooth (or time-implicit) DEM Servin, Wang, Lacoursière, and Bodin] for the particles, and rigid multibody dynamics for the asset's moving parts using feedforward controllers. The pseudo-particle size range between 0.2 and 0.48 m, which corresponds to a mass range between 12 to 160 kg. The distance between the two loading points is roughly 150 m and the stockpile is 60 m tall. The simulations were run with time-step 11 ms and 100 projected Gauss-Seidel solver iterations, which corresponds to an error tolerance around 5 % for the stockpile bulk behavior Servin, Wang, Lacoursière, and Bodin].
At the sources, each particle is assigned an identity and a location of origin. The left source is assigned a mineral concentration of 0.2, while the right source is assigned the value 0.8. The right source is given a mineral hardness of 1.0, while the hardness at the left source varies over the ore body according to some function h(x). It is tested how a step in the hardness distribution function, from value 1.0 to 2.0, in the left source propagates through the system. It is also investigated how well a hardness profile in the ore body can be backtracked from observations at the crusher and the grinder. The (true) particle size distribution is discretized in N d = 16 size classes with a logarithmic binning. The initial size distributions is shown in Fig. 5, including also the average size distribution after the crusher. Do note that the pseudoparticle diameter, used for visualization and contact dynamics, have no correspondence to its internal particle size distribution. Each haul truck receive a load of roughly 15 ton (roughly 150 pseudo-particles) at the loading point. For simplicity, the loading equipment is not modeled explicitly. The trucks dump the load on the crusher inlet and there occur mixing of the material sources from here on. The crusher is modeled as a jaw crusher, drawing material from above and crushing it with a jaw driven in an oscillatory motion. The pseudo-particles retain their original diameter, but when they pass the jaw, the internal particle size distribution is transformed. The fractions larger than the crusher gap width, dcr = 0.25 m, break into a uniform distribution of the smaller size classes. This unit operation is represented by fn := fn − ∆fn− + ∆fn+, with the change vector for crushed (oversize) fractions, ∆f k n− = f k n for d k ≥ dcr (otherwise zero), and the change vector for the n d undersized fractions ∆f k n+ = 1 n d k ∆f k n− ( otherwise zero). The resulting size distribution is shown in Fig. 5. We assume that the torque for driving the jaw crusher stands in proportion to the amount of mass in the crusher zone Vcr(t), its hardness and size distribution, according to the following model where the crushing coefficient k n cr = 1.0 for the oversize and k n cr = 0.1 for the undersize. The momentaneous power draw of the crusher is computed by Pcr(t) = τcrωcr, where ωcr = 10 rad/s is the crusher speed.
The stockpile holds up to 2, 000 tons (roughly 20, 000 k pseudo-particles). There are two entry discharge points and four exit discharge points, but only one is used for the presented tests. When material is discharged from the bottom outlet, at a controlled speed ranging up to 2 m/s (or 2 ton/s), a funnel flow field is activated in the pile, reaching the surface where a depletion is gradually formed, unless new material is fed into the storage at the same rate or faster. Two time instances of pile states are shown in Fig. 6. For the pointy pile state one can expect more mixing of the incoming material and longer residence time, contrary to the depleted pile state where material flow more directly through the storage and with less mixing. An idealized mill model is assumed, guided by a model that has been used for mill control at Boliden Aitik [Araker and Bostrom(2020)]. The idealized model consists of a single grinder, although most real mills consists of a circuit of grinders and classifiers. When a pseudo-particle n reach the the grinder, a size reduction model begins to operate on its size distribution. Each size fraction k is assigned a probability per unit time, P k n , for breaking into smaller fractions and a probability per unit time, Q k n , for leaving the mill. At each time-step the mass fractions carried by each pseudo-particle are redistributed according to fn := fn − ∆fn− + ∆fn+, where the broken and passing fractions are and the reduced fractions are Modeling a size-classifier, the probability per unit time for the finest fractions, smaller than d cl , to leave the mill is Q k n = kesc/(1 + exp[k cl (d k /d cl − 1)]). We use d cl = 0.025 mm, k cl = 10, and kesc = 10 3 . The breakage probability depends on the (true) particle size, impact energy, hardness, and the number of impacts [Bruchmüller et al.(2011)Bruchmüller, van Wachem, Gu, and Luo], which in turn depends on the mill's angular velocity, ωgr, and the active torque τgr-a. We assume the following model for breakage probability per unit time where is an impact energy distribution function, E k n,brk = E 0 brk hn/[1 + (d k /dgr) 1/4 ] is a specific breakage energy distribution, D50 is the (timedependent) mean diameter of the material currently in the grinder, dgr = 1.0 mm is a grinding length-scale parameter, and kgr = 1.0 a dimensionless breakage coefficient. The internal drive friction, τ 0 gr-dr , of an empty mill and the critical mill speed, ωcrit = 1.5 rad/s, are used for normalization. The mill breakage model can be seen in Fig. 7, withP(d k ) ≡ P(d k )/ τgr-aωgr τ mill ω crit . The exit probabilityQ(d k ) ≡ Q(d k )/kesc is also included.  Figure 7: The mill breakage model as function of particle size and mineral hardness illustrated in terms of the impact energy distribution (a), specific breakage energy (b), and the breakage probability rate (c). The probability rate for exiting the mill is also included.
The total grinding torque is the sum of three contributions, τgr = τ gr-dr + τacc + τgr-a. The total mass, mtot(t) = m0 + mgr(t), is the mass of the moving part of the mill, m0 = 10 3 kg, plus the mass of the material in the grinder, mgr(t) = n mn. The internal drive friction is modelled τ gr-dr = µ dr mtotgr mill and the torque required for accelerating the grinder is τacc = mtotr 2 millωgr , where we denote the internal friction coefficient µ dr = 0.02, gravity acceleration g = 9.8 m/s 2 , and mill radius r mill = 5 m. The active torque depends on the charge's mass, mgr = n mn, size distribution, and dynamic angle, which depends on the mill angular velocity. We assume, again guided by [Araker and Bostrom(2020)], that the active torque can be modeled τgr-a = µgr-amgrgr mill with effective internal friction of the charge µgr-a = kgr-aΘ(ωgr/ωcrit) mgr-aωgr mgrω crit . Here, mgr-a = n,k≥Ngr-a m k n is the mass of the charge fractions larger than an active grinding size dgr-a, with index Ngr-a. We use dgr-a = 1.0 mm, which is Ngr-a = 7 in the 16 class size representation, and kgr-a = 1. With Θ(x) = Re[(1 − x 4 ) 1/2 ] the effective internal charge friction and active torque drops to zero at and above the critical velocity. The power draw of the grinder is computed Pgr(t) = τgrωgr. When nothing else is stated, the mill is run at constant nominal speed ωgr = 1.2 rad/s.

Step response
The simulator is initialized with ore of constant hardness h = 1 transported from both sources. A step change to hardness h = 2 is introduced at the left source and the response in the crusher may be seen in Fig. 8. The haul truck reach this harder material at time t = 350 s. About 60 s later (five seconds more than one truck cycle) the crusher power draw has increased by 40%, from 1.5 kW to 2.1 kW. A virtual sensor in the crusher reveal that the mineral hardness varies between 1 and 2, with a time-averaged hardness around 1.4. The large oscillations in hardness and power draw, with 55 s time period, is because of the trucks dumping intechangely from the left and right source with hardness 1 and 2, respectively. For the step response in the mill we analyse the effect of the ore passing through the stockpile in Fig. 6. The discharge of the piled material start at 500 and the feed reach the empty mill shortly thereafter. The material in the mill remain of hardness h = 1 for 90 more seconds until it starts increasing at time 590 s up to h = 1.5 with only small oscillations (well mixed). At time 875 s, the storage level is down at zero (apart from the dead zones) and the rate of mass flowing into the mill drops quickly to the rate of mass leaving the crusher and flowing into the storage. Therefore, after 875 s the oscillations in hardness and power are again reflecting the truck cycles. For reference we run the identical simulation with constant hardness h = 1 at both sources. As seen in Fig. 9, the difference in the two power draw signals follow the increase in material hardness inside the mill. The increase is about 20 % on average. Note that the mechanism for the higher power draw is that the harder ore takes longer time to grind and the mill therefore accumulate more mass over time.  Figure 9: The response in the mill power draw (a) from a step change in the ore hardness (b) after passing through the storage in Fig. 6. The ore mass in the mill is also included and a 30 s moving average of the power draw and hardness. . For reference, the case of constant hardness (h = 1) is also included.

Identification of hardness from measurement-while-drilling
To test the feasibility of identifying the rock hardness from the MWD signal we assume the quadratic hardness profile h(x) in Fig. 10 at the left source, ranging continuously between 1 and 2. While drilling, a measurement profile s(x) is measured. For the test we assume s = h 0.4 and introduce a measurement location error such that sn = h(xn + R(−∆x, ∆x)) 0.4 , where R(−∆x, ∆x) is a uniform disturbance. We use ∆x = 12.5 m, which corresponds to 6 haul truck cycles back and fort to the crusher. The right source has constant hardness h = 1. Ore from the two sources enter the crusher, partially mixed, and give rise to the torque signal in Fig. 11.   Figure 11: The observed crusher torque compared to the expected torque from different model ansatz given the MWD signal and hardness. Four truck load reach the crusher during the shown time period, two from the left source with harder ore and two from the right source.
Assuming the crusher torque is related to rock hardness according to Eq. (4) any ansatz s(h, β) can be explored and optimized with respect to model parameter β. Different outcomes for the ansatz s = h β are shown in Fig. 11. The best fit was found at β = 0.402, which deviate by 0.5 % from the ground truth because of the tracking location error that was introduced. Once a relation between the MWD signal and hardness has been established it is possible to predict the hardness of the material that is about to enter the crusher given the MWD signal measured at its location of origin. This is illustrated in Fig. 12 over a sampling duration of 2, 500 s, or roughly 40 haul trucks from the left source.

Predicting the properties of the mill feed
With particle-based tracking and Eq. (3) is straight-forward to monitor the properties of the ore kept in the storage as function of time and space. Given a storage flow model it is also possible to estimate which particles in the funnel zone that will leave the storage in the near future, say the next 30 s. This way the properties of the ore feed that will reach the mill can be predicted in advance, even before exiting the storage. Sample spatial distributions are shown in Fig. 13. The time evolution of the properties in the pile and of the future feed are shown in Fig. 14.

Backtracking the grindability
The grindability of an ore may be measured by the amount of energy, per unit mass, that is required for reducing the fragmentation down to a specific particle size. For the sake of Figure 13: The spatial distribution of velocity (left), concentration (centre) and hardness (right) in the pile at a time instant. The material that will be discharged in near time may be estimated from the funnel flow zone (left).  Figure 14: Evolution of the storage average properties over time (left) and the predicted feed for the next 30 s (right) given the funnel flow zone in Fig. 13. A 2 s moving average is included (orange).
production planning the grindability is often estimated in advance using the relatively sparse exploration data of the ore body. The Bond ball mill work index (BWI) is defined where the work per unit mass (W/kg) is computed W (t) = Pgr / ṁgr-out using a moving average with time-window tav. D80 is the 80% passing diameter (m), computed from the cumulative size distribution function for mass that enter and leave the mill. The Bond grindability is then defined as the inverse G ≡ BWI −1 . We demonstrate the grindability measured in the mill can be backtracked to ore's location of origin for the purpose of improving future estimates for grindability. Assume there is an a priori, G 0 i , estimate of the grindability of the ore in block i. At time t k a grindability G(t k ) is observed at the mill. The average mass in the mill originating from the block i is m i gr (t k ). We reconciliate the grindability at block i and time t k by the weighted average where the weight factor is given by the ratio of mass with origin i to the mass in the mill and relative to the block mass αi(t k ) = m i gr 2 /m 0 i mgr . Here, mgr (t k ) is the amount of mass in the mill and m 0 i is the total mass of the ore block i. Note that the weight factor is small (or zero) when there is little (or no) mass from the considered block. The weight factor for the original grindability reflect the confidence of that estimate, α 0 i 1 if the confidence is low and α 0 i 1 if the confidence is high. We test backtracking the grindability of the ore originating from the left source for the case of a quadratic hardness profile shown in Fig. 10, but mixed in the mill with ore from the right source. We use a block size of 10 m and α 0 = 0.1. The BWI signals are shown in Fig. 15 and the grindability backtracked to the source using Eq. 9 is found in Fig. 16.

Discussion and conclusion
A digital twin of a mine may be considered as a system of particle simulations, distributed over virtual assets that evolve according to their physical counterpart using realtime particle simulation where no real sensors or control signals are available. The simulator tests demonstrate how data analysis and state predictions are made possible by the particle-based material tracking, i.e., in situations where material from different sources are mixed and alter the comminution processes by their natural variations. Knowing better what material is currently in storage and about to exit enable more efficient milling, by more precise control of the mill speed or by blending ores from different sources. Understanding how fragmentation and hardness affect the entire chain of processes from loading, crushing to milling, provides the data for deciding what is the most optimal fragmentation from a mine-to-mill perspective. The realization of particle-based material tracking in a real mine involve creating virtual assets for each of the physical assets, a central historian and tracking service, and interconnecting these in a distributed system. In a modern mine, most of the components and infrastructure for this is already in place for the purpose of local production planning and optimization. It is mainly for the storage systems there is a lack of models that provide sufficient resolution and performance for driving the virtual assets. Particle-based simulation, such as DEM, accelerated using reduced order modeling techniques may fill this need. The main uncertainties is how sensitive the result is to the resolution of particle size and shape, and to the DEM model parameters including particle mass density, friction etc. As demonstrated here, this can be investigated in advance using a simulator with material tracking. Future work should investigate the required level of fidelity of the simulation modes to sufficiently capture how the material disperse over pile surfaces and the velocity field in funnel discharge flow. Other unknown factors to investigate are the effect of moisture and temperature on the flow properties. Calibration and validation tests involving physical tracers is an important tool here.