Integrated geometallurgical modelling of heavy mineral sands
accounting for profitability, extractability and processability under
uncertainty
Tom Wambeke, J
ӧ
rg Benndorf
Delft university of technology
ABSTRACT :
The development of a mineral resource is inherently connected with a substantial amount of risk due to the large financial investment required at the time that the geological knowledge is rather limited and the comprehension of the so-called “modifying factors” is incomplete. Ongoing research has shown that a great deal of the risk can be mitigated by considering a probabilistic resource model instead of a single determinist one. Unfortunately, too often probabilistic resource models are only restricted to the characterization of uncertainty in metal grades.
This paper presents a coherent framework with geostatistical modeling techniques to assess the geological uncertainty and spatial variability of all the factors impacting a heavy mineral sands operation. The framework integrates profitability through grade characterization with extractability and processability through the modeling of slime content (fine clay particles), oversize (fraction greater than 2 mm) and in-situ density. Since the spatial distribution of these properties is governed by the geological facies, it is necessary to first address the problem of identifying of identifying the geological domains. The challenge lies in the simulation of the transitional behaviour of especially the slime content and oversize, but also grade content, across soft boundaries. A new simulation methodology was developed to accommodate these boundary effects. Moreover, the methodology is especially designed to facilitate the updating of the resource model later on when production data becomes available.
The integrated effort results in a comprehensive characterization of resource uncertainty. Eventually this approach will result in the only “real” risk-robust design option and operational strategy.
1
Introduction
A recent study has shown that the global heavy mineral sand industry is on the verge to experience an enormous growth due to the growing urbanisation trend, ceramic tiles consumption and demand from the aerospace industry (Daedal Research 2014). In order to facilitate this growth, new mines need to be developed and the production rates of existing ones needs to be ramped up. Beware that increasing the supply also entails a great deal of risk. Investments need to be made at times where the lack of data induces considerable uncertainties.
In order to properly characterize these uncertainties, the geological complexity of a heavy mineral sand deposit needs to be understood and attempts need to be made to integrate it in a realistic way into the modeling exercise. Too often, the focus solely lies on the spatial characterization of the grades. However, a resource efficient evaluation along the complete mining value chain requires the understanding of the properties impacting the value, extractability and processability of the material.
Table 1 illustrates how this concept of integrated geometallurgical modelling needs to be applied in the context of heavy mineral sand. Since the location and the interaction between the geological facies have an important impact on the spatial variability of all relevant properties, the facies need to be modelled first. Subsequently the total heavy mineral (HM) grade is modelled due to its direct impact on the value of the deposit.
The slime content and the oversize, on the other hand, both influence the extractability and processability of the sand. A high slime content increases the risk of the formation of clay balls in the processing plant and results in a long settling time of the slurry in the pond. Oversized particles are screened and do not even enter the processing plant.
Changes in the slime content and oversize are inextricably connected with different grain size distribution curves. These curves together with the heavy mineral content can be used to estimate the in-situ density and compaction which determine the required power and forces needed to excavate the material.
Tab.1 : Concept of integrated geometallurgical modelling in the context of a heavy mineral sand.
Subject Modelling component Relevance under subject
Geology - Domains / structures control spatial variability of parameters below
Economy - Grade value of product sold
Processing - Slime
- Oversize
design processing plant design dredge pond scheduling operation processing costs per block design processing plant scheduling operation processing costs per block Excavation - In-situ density &
Compaction
design / selection equipment excavation costs per block
The paper first presents a brief introduction into the geological genesis of heavy mineral sand deposits to better comprehend its complex nature and unique features. Thereafter, a coherent modelling approach is presented tailored to the geological complexity and the needs for characterizing parameters influencing extraction and processing. Subsequently a case study is presented. Finally some conclusion are formulated.
2
GEOLOGY OF HEAVY MINERAL SAND DEPOSITS
Mineral sands consists of a group of minerals which were concentrated in alluvial environments due to their difference in specific gravity compared to the secondary minerals (mainly quartz). The principle valuable minerals are rutile, ilmenite, leucoxene (sources of titanium) and zircon. Due to their high specific gravity, these minerals tend to lag or concentrate during storms when lighter components such as quartz are carried offshore or along shore.
The heavy mineral (HM) assemblage predominantly reflects the crystalline basement of the local provenance. During periods of fair weather, sediments are brought down by rivers and accumulate in beaches. The HM grains are scattered throughout the quartz sand with background concentrations of less than 0.5% HM (Figure 1, left). These beach accumulations provide the basis for the thicker strandlines mined today.
These strandlines were formed over centuries or millennia of repeated storm events causing large waves, massive erosion and substantial reworking. Enrichment occurs through winnowing out the lighter quartz minerals leaving the heavier mineral grains behind. During this process, large volumes of quartz are moved offshore (Figure 1, middle).
The preservation of the accumulation occurs over geologically longer periods when fair weather occurs in combination with a marine transgression. The changing sea level moves the quartz sand back onshore, rebuilding the beach and preserving the strandlines (Figure1, right) (Jones 2008).
Abb.1 : Conceptual overview of the accumulation (left), concentration (right) and preservation (right) of a marine beach placer (Jones 2008).
The ongoing winnowing and the subsequent gradual burial obscure the exact location of the geological boundaries. The properties of interest, especially the slime content and the oversize, exhibit a rather complex transitional behaviour across the so-called soft boundaries.
To further complicate the matter, the geological information collected during the exploration campaign is not detailed enough to give any indication on the location of different geological domains. Therefore, the definition of mineralised boundaries is often based on a straightforward grade cut-off.
3
INTEGRATED GEOMETALLURIGAL MODELLING
Current practice in the heavy mineral sand industry is to use some software to manually delineate the mineralized zones in cross sections based on a predefined cut-off grade. Subsequently a wireframe is constructed and inverse distance weighting or kriging techniques are applied to estimate the grades inside. Reconciliation exercises have shown that this approach generally works effectively (Jones 2008). Nevertheless, there is still an urgent need to further improve the current practices for the following three reason:
1. The current estimation methodology results in a single “best” guess and thus it neglects the inherent uncertainty of the spatial modelling exercise.
2. The definition of a cut-off grade, before any design or planning decisions are made, has an enormous impact on the later modelling steps and decisions to be made.
3. The application of geometallurgy could lead to considerable improvements along the complete mining value chain.
The issue of characterizing the uncertainty of geological facies already lingers for more than 20 years in the geostatistical community (Goovaerts 1997, Lantuéjoul 2002, Armstrong et al. 2004, Strebelle 2002). Several solutions have been developed, all aiming to assess the volumetric uncertainty through the generation of multiple realizations which reproduce the complex geological features found in mineral resources. In order to create a more complete picture of the overall uncertainty in the deposit, these domain realizations can be further combined with separate grade simulation inside the geological domains.
Although the previously presented methodology seems to work relatively well for modelling certain mineral resources (Jones et al. 2013), its use cannot be easily extend to the modelling of heavy mineral sands. The consecutive stochastic simulation of distinct geological facies and orebody properties would result in different zones separated by hard boundaries, each displaying a different spatial variability. In this way, one of the key characteristics of heavy mineral sand deposits, i.e. the transitional behaviour across the soft boundaries (complex vertical trend) is completely ignored. It
is clear that there is a need for a deeper integration between geological modelling and grade, slime or oversize simulations.
Considering the arguments elucidated before, it would not be desirable to split the entire sand body into distinct stationary zones which properties can be simulated separately. Neither the assumption of one single stationary zone for the entire deposit would be correct. Instead of first explicitly modelling geological facies and then superimposing a stationary mean (globally or locally constant), a more complex trend component is modeled that implicitly accounts for a transitionally changing mean (and thus also the soft boundaries between geological domains).
This means that the attribute under study Z(x) now has a mean or trend component m(x) that is modelled in function of its location. Furthermore a weakly stationary random function S(x) and a noise component R(x) are superimposed on the trend.
Z = S + m + R
Conventionally, the trend describes the mathematical expectation at a certain location and is considered to be deterministic but unknown. This component is generally modeled through a linear combination of constant, linear, bi-linear and quadratic functional terms depending on absolute or relative spatial coordinates. At each spatial location, the functional terms are evaluated and universal or dual kriging equations are solved to calculate the combinatorial weights. Adjusting the size of the search neighborhood allows for the estimation of local complex trends. These kind of trends can potentially capture the transitional behaviour across boundaries which are actually unknown.
Instead of describing the complex trend as a linear combination of some spatial functional terms, a more extended approach was developed where probability or distance fields add further information to the equations.
Although it would not be preferable to use mineral data and a cut-off grade to delineate different stationary domains, it could be used to generate soft information for conditioning the trend function. The grade information generally provides us enough information to get an indication on the location of the strandlines. Once the outline of the strandline is approximated, a distance field can be constructed. For each point in the field, the distance to the closest point on the centre-line of a strandline is calculated. This distance is a good indicator to condition the complex trends observed in heavy mineral sand deposits.
Furthermore, the distance field can be transformed into a probability field. The closer the distance to one of the centre-lines, the more likely a point belongs to a strandline. The advantage of such a probability field is that it is easier to combine it with other sources of information.
Once the deterministic trend component is obtained, the remaining uncertainty and spatial variability is captured by a spatial random function. Such a function defines a spatial law that characterizing the joint uncertainty about a set of values. In geostatistics, the random function model often reduces to a set of cumulative density functions that involve not more than two locations at a time. In most cases, the first and second order moments are considered which are the expectation and the covariance function. This spatial random function is subsequently used during simulation to generate different possible realizations.
Although presented as a clear two phase procedure, it does not come that easy. In order to solve the dual or universal kriging equations, the covariance function of the residuals needs to be fitted. In order to calculate their empirical covariance function, the trend needs to be removed from the data. This is the same trend which we initially are trying to model.
Conventionally, an initial guess results in a basic covariance model for the residuals. This model is used in the kriging equations to estimate the complex trend. The trend is removed from the data
with the empirical covariance function and adjusted accordingly. This procedure continues until the inserted covariance model matches the one calculated from the resulting residuals.
Please notice that changing the model parameters of the inserted covariance function also influence the course of the empirical one calculated afterwards in an unpredictable manner. Therefore this iterative adjustment procedure is very cumbersome and time consuming and in this way a point of satisfaction might never be reached.
Recently a new method was presented (Benndorf & Menz 2013) that can be used to evaluate the fitting of geostatistical model parameters in the presence of a trend. The method extends conventional cross-validation approaches in two ways. Firstly, the parameter evaluation is not only limited to data locations, but is also performed on all prediction points in the grid. Secondly and more important, the method replaces the error measures used in cross-validation, based on single point substitutions, with a theoretical and empirical error curve. These curves result from calculated errors related with different rings of influence. The separate rings correspond to different variogram lags and facilitate the assessment of the fit of the complete variogram structure and its related parameters. “If the model parameters, which are used for the calculation of both measures of errors in prediction, are neither fitting the data measured, nor the structural behaviour of the attribute under consideration, a discrepancy between the theoretical and empirical error will occur (Benndorf & Menz 2013).”
To achieve a good model fit, the more sensible theoretical error curve needs to matched to the empirical error curve via iterative systematic adjustments of the covariance model parameters. Horizontal deviations indicate that the range of the covariance model needs to be adjusted. Vertical deviations are controlled using the variances. This provides a mean for systematically improving the model fit, compared to the previously discussed “guess & check” method where the influence of different model parameters cannot be explicitly discriminated.
Simultaneously with the calculation of the error curves, a curve of the universal kriging variance as function of the ring of influence can be constructed. Such a plot can be very useful since it allows an objective comparison between different types of trend functions. The one with the lowest universal kriging variance provides the best type of trend function for the data.
4
CASE STUDY
At the time of writing, the proposed methodology is being extensively tested at a complex heavy mineral sand deposit in Western Australia. The area of investigation, the so-called Rhea target, is located about 160 km north-northwest of Perth, and 25 km east of the coastal community of Cervantes. A total of 401 aircore holes were drilled on a semi-regular 20m x 100m grid overlying an area of about 2.6 km2. The exploration data is made publically available by the Australian government (Image Resources 2011).
The Rhea target lies in the Swan Coastal Plain and contains several beach deposits which were formed during a Quaternary cyclic regression phase. These beach deposits overlie Mesozoic fluvial deposits which constitute the basement in many regions along this paleocoastline. The “geology” of the Rhea target consist of three obscure horizons which host mineralization (Figure 2):
1. Strandline mineralization occurs mainly in the west of the tenement, are located close to the surface, are narrow and contain medium to high grades.
2. Mid-level mineralization is interpreted as a broad sheet which lies directly on the unconformity (basement/channel) at a depth averaging 15-20m.
3. Channel mineralization occurs in alluvial sediments within valleys cut into the Mesozoic basement. These channels predate the Quaternary marine sediments.
Abb.2 : Interpretation of a cross section of rhea target showing the strandline, mid-level and channel mineralization (Image Resources 2011)
Several geostatistical (simulation) techniques can now be applied to model the different mineralization horizons. Since our methodology does not consider the modeled domains as hard boundaries but merely as soft additional information, there is no need to resort to complex paradigms such as for example multi-point statistics or pluri-Gaussian simulations. Instead a more straightforward approach was opted that builds up the three-dimensional geology based on two-dimensional layers. The geological modelling goes as follows (Figure 3):
1. In each borehole use the available information, mainly the grade profile, to identify the depth of the basement, the bottom and thickness of the middle zone and the top and thickness of the strandline. At some locations this interpretation can be rather subjective, which only further supports our decision to use the boundaries as soft information.
2. The geometry of the channel will directly result from the modelling of the basement depth. Since a clear trend is present , this needs to be modeled first. Subsequently the residuals at the sampling points are calculated and their unknown values in the remainder of the field are simulated. Finally the trend is added back to a set of realizations to obtain an indication of the basement uncertainty.
3. Based on the borehole interpretation, both the bottom depth of the middle sheet and its thickness are simulated. Combining both sets of realizations results in a three dimensional sheet volume.
4. Based on the borehole interpretation, both the top depth and the thickness of the strandlines are simulated. Combining both sets of realizations results in a three dimensional picture of the strandlines.
5. Combining the five obtained sets of realizations (basement depth, bottom sheet, thickness sheet, top strandline & thickness strandline) results in realizations reflecting different geological scenarios. These scenarios are subsequently used to condition the complex trends.
Figure 3 give a schematic overview of the geological modelling procedure. Once the different layers are obtained, one can calculate for each grid cell in a three dimensional space the distance to each horizon. The obtained distances are subsequently inserted into the design matrix used during the construction of the dual kriging equations. Eventually this results in a trend that is conditioned on the modeled geology.
Abb.3 : Methodology to model the geological horizons which will be used later on to condition the complex trend.
This procedure will be further illustrated on a synthetic, but realistic, two dimensional example. To start with, a short explanation on how the data was obtained is appropriate. First, a complex vertical trend was constructed, which indicates the presence of two different mineralized horizons, one with its centre about 20m below the surface, the other just above the basement (Figure 4, right). Subsequently, a spatial random function and some white noise were added. Finally ten boreholes were drilled in the resulting synthetic field. The resulting borehole profiles are displayed on the top of figure 4 and resemble the observed grade behaviour of the Rhea target.
A more detailed investigation of the synthetic data is already very instructive. Although the data was clearly constructed based upon two mineralized horizons with soft boundaries, the superposition of some spatial variation really clutters the interpretation and delineation.
Nevertheless based on two cut-off values, a marker for the top strandline and one for the sheet above the basement could be constructed. These markers were subsequently transformed to distance fields and inserted into the design matrix of the dual kriging equations. The residual covariance function was fitted using the above discussed approach of Benndorf & Menz (2013). Subsequently the residuals at the sample locations were calculated and a simulation was initiated. On completion, the simulated residuals were superimposed on the previously modeled trend. The result is shown at the bottom of figure 4.
The resulting realization exhibits the realistic behaviour of the spatial variation of the grades in a heavy mineral sand deposit. Clear mineralized horizons can be observed, although their exact delineations are difficult to determine due to the transitional behaviour across the soft boundaries. Similar simulations are performed to additionally characterize the spatial variability and geological uncertainty in the other key properties impacting the mining and processing operations (slime content & oversize).
Abb.4 : Synthetic example of the spatial grade distribution. Complex vertical trend mimicking the behaviour observed in the boreholes of the Rhea target (right). Ten boreholes drilled in the synthetic field, which was created by the
superposition of the complex trend, a spatial random function and white noise (top). A simulated realization of the spatial grade distribution based on the ten boreholes (bottom).
5
CONCLUSION
The paper presents an integrated geometallurgical modelling framework to assess the geological uncertainty and spatial variability of all the factors impacting a heavy mineral sands operation. The initial results of a novel research concept, that integrates geology deeper into the geostatistical modelling exercise, are presented.
Generally geological facies are not sampled on a detailed level, making it almost impossible to properly delineate different geological horizons or facies. Typically, horizons are defined based on a single grade cut-off. The outline of such horizons are subsequently considered as hard boundaries and the transitional behaviour across such a boundary that is actual soft, is completely ignored. Although this procedure might still be slightly acceptable for the modelling of the spatial grade distribution, it is definitely not appropriate for the modelling of the slime content and oversize.
In order to generate more realistic realizations, a new simulation methodology was developed to accommodate these boundary effects. The transitional behaviour across the soft boundaries is integrated as a complex trend which is conditioned on the geology. This conditioning is achieved through distance fields which enter the dual kriging equations via the design matrix. Once the complex trend is modeled, the remaining uncertainty is characterized via the simulation of residual fields. Afterwards the trend is added back to these residuals to obtain a more realistic image of the variability of the key property.
The integrated effort results in a comprehensive characterization of the resource uncertainty. Propagating this uncertainty further into the mining value chain will result in the only “real” risk-robust design option and operational strategy.
6
REFERENCES
Armstrong, M, Galli, A, Beucher, H, Le Loc’h, G, Renard, D, Doligez, B, Eschard, R, Geffroy, F 2011.
Plurigaussian simulations in geosciences, Springer, Dordrecht.
Benndorf, J, Menz, J, 2014. Improving the assessment of uncertainty and risk in the spatial prediction of
environmental impacts: a new approach for fitting geostatistical model parameters based on dual krigigin in the presence of a trend, Stoch, Stochastic Environmental Research and Risk Assessment,
28, 3, pp. 627-637.
Deadal Research 2014. Global mineral sand industry: trends and opportunities (2013-2018), Avaialvle from: http://www.daedal-research.com/global-mineral-sand-industry-trends-and-opportunities-2013-2018 [March 9, 2014]
Goovaerts, P 1997. Geostatistics for natural resource evaluation, Oxford University Press, New York. Image Resources 2011. Final report, EIS co-funded drilling Cooljarloo Project E70/2898 North Perth Basin
Western Australia. Available from: http://www.dmp.wa.gov.au/5136.aspx [November 11, 2013] Jones, G 2008. Mineral sands – an overview of the industry.
Available from: http://www.iluka.com/docs/company-presentations/mineral-sands---an-overview-of-the-industry-by-greg-jones-manager-development-geology March 9, 2014]
Jones, P, Douglas, I, Jewbali, A, 2013. Modeling combined geological and grade uncertainty: application of
Multi-point simulation at the Apensu gold deposit, Ghana, Mathematical Geolgoy, 45, 8, pp. 949-956.
Lantejoul, C 2002. Geostatistical simulation: models and algorithms, Springer, Dordrecht.
Strebelle, S, 2002, Conditional simulation of complex geological structures using multiple-point statistics. Mathematical Geology 34, 1, pp. 1-21.
