INTRODUCTION
Mineral sand refers to a group of heavy minerals which have a density of more than 2.85 g/cm3. They are processed to obtain two core product streams: zircon and a collection of Ti-bearing minerals such as ilmenite, rutile and leucoxene. Mineral sand resources remain to be of great importance due to the continuing global urbanisation and the resulting increase in the demand of paint, ceramic tiles and other construction materials (McKinsey & Company, 2012; Daedal Research 2014).
A resource efficient evaluation calls for a proper characterisation of the uncertainty in the properties impacting the profitability, extractability and processability of the raw material in-place (Ferrara, Guarascio and Massacci, 1993; Vallee, 2000; Dominy, Noppé and Annels, 2002). The total heavy mineral grade is modelled due to its direct impact on the economic value of the deposit. The slime content and the
oversize both influence the extractability and processability of the sand. Table 1 summarises the critical geometallurgical parameters which are readily available for geostatistical modelling.
Modelling the relevant attributes requires a sound understanding of the geological processes involved in the formation of the deposit (Jones, 2009). During periods of fair weather, sediments are brought down by rivers and accumulated in beaches. The heavy mineral grains are scattered throughout the quartz beach sand with background concentrations of less than 0.5 per cent (Figure 1a). Enrichments occurs when large volumes of lighter components are carried offshore or along shore by storm-induced littoral drifts. The heavier mineral grains tend to lag and concentrate while the lighter quartz particles are winnowed out (Figure 1b). A marine transgression in combination with fair weather during
Optimal Inference of Modelling
Parameters to Simulate Complex
Trends across Soft Boundaries –
A Case Study in Heavy Mineral Sands
T Wambeke
1,2and J Benndorf
3ABSTRACT
A risk-robust development of a heavy mineral resource requires an assessment of the geological uncertainty and spatial variability of the key factors impacting the mining and processing operation. Attributes of interest are the total heavy mineral grade, the slime content and the amount of oversized particles. The modelling of a mineral sand deposit presents two main challenges. Firstly, the exact location of the geological boundaries are obscure due to the complex interaction of the deposit formation processes. Secondly, the properties of interest exhibit a rather complex transitional behaviour across the so-called soft boundaries. Traditionally, the modelling exercise is restricted to a consecutive stochastic simulation of distinct geological facies and orebody properties. This results in different zones separated by hard boundaries, ignoring the typical transitional behaviour.
This contribution presents an alternative approach which models the entire sand body as a single domain while accounting for the complex trends present. Error curves are used to infer covariance model parameters without the need of a prior removal of the trend. The principle is based on minimising the differences between empirical and theoretical errors resulting from two methods of prediction calculated for different rings of influence. The added value results from the unambiguous decomposition of the spatial variation into a deterministic trend and a stochastic residual. The trend model neither overfits nor underfits the data. The calculated residuals are thus unbiased and the uncertainty can be correctly assessed.
The effort results in a characterisation of the resource uncertainty which provides a more realistic image of the spatial variability of the key properties in and across the geological domains of an existing heavy mineral sand deposit.
1. PhD candidate, Delft University of Technology, Stevinweg 1, 2628CN Delft, The Netherlands. Email: t.wambeke@tudelft.nl 2. Research Engineer, IHC Merwede, Smitweg 6, 2961 AW Kinderdijk, The Netherlands.
a geologically longer period of time results in the preservation of the enriched zones. The changing sea level moves the quartz sand back onshore rebuilding the beach and preserving the strandline (Figure 1c).
The strandlines mined today were exposed to repeated cycles of deposition, substantial reworking and massive storm erosion. The ongoing winnowing and supply of new sediments obscures the exact location of the geological boundaries. Due to these complex interactions, the properties of interest exhibit a rather complex transitional behaviour across the so-called soft boundaries. Figure 2 shows a vertical profile of measured attributes in a drill hole from the Rhea
target (to be discussed later on in the case study). The complex trends in the vertical direction are immediately apparent.
Traditionally, the delineation of mineralisation domains is often based on a straightforward cut-off grade. The introduction of a cut-off grade at this stage implicitly constrains design or operational decisions still to be made. The complexity of the geological processes directly affects the behaviour of the attribute and renders this conventional ‘geological’ zonation highly unreliable and subjective. This contribution presents an alternative technique which models the entire sand body as a single domain while accounting for complex trends present.
In the subsequent sections a new integrated modelling approach will be outlined and demonstrated. The approach is based on the dual kriging paradigm to build a complex trend model which is superimposed on conditional simulations of the residual. The applicability of the new methodology resides in a straightforward and objective inference of the model parameters, which does not require a prior removal of the complex trend. Error curves are used to ensure a good fit of the complete covariance structure along the different lag distances. The principle is based on minimising the differences between theoretical and empirical errors resulting from two methods of prediction calculated for different rings of influence. The integrated modelling of the trend and covariance model parameters results in an unambiguous decomposition of the spatial variation into a deterministic trend and a stochastic residual. The trend model neither overfits nor underfits the data. The calculated residuals are thus unbiased and the uncertainty can be correctly assessed via simulations.
First, the paper presents a short review on the current state-of-the-art in the modelling of a non-stationary mean (trend component). Thereafter, the new approach for fitting the
Subject
Modelling component
Relevance under subject
Geology Domain/structures Controls the spatial behaviour of attributes below Economy Slime content (<63 mu) Design processing plant
Design dredge pond Scheduling operation Processing costs per block Oversize (>2 mm) Design processing plant
Scheduling operation Processing costs per block Excavation In situ density and compaction Design/selection equipment
Excavation cots per block
TABLE 1
Overview of the available geometallurgical parameters in a
heavy mineral sand operation. The impact of each parameter
on a later stage in the value chain is defined.
FIG 1 – Conceptual representation of the formation of a marine beach placer deposit. Three consecutive events
occur: (A) the accumulation; (B) concentration; and (C) preservation (Jones, 2009).
A
B
C
FIG 2 – Profile of the measured grade, slime and oversize content in borehole RH304 of the Rhea target, Western Australia. A
geostatistical model parameters in the presence of a trend is explained. A first case study demonstrates the practicability of the new approach in a known environment. The second case study illustrates that the overall simulation methodology is capable to model complex behaviour observed in a real heavy mineral sand deposit.
A SHORT REVIEW OF MODELLING A
NON-STATIONARY MEAN
In many geostatistical applications, the attribute under study
Z(x) is modelled as a combination of a weakly stationary
random function S(x) and a deterministic trend component
m(x) (Journel and Huijbregts, 1978; Goovaerts, 1997). Further
a noise component R(x) can be superimposed to account for measurement error or variability on a scale much smaller than the sampling interval. This results in:
Z(x) = m(x) + S(x) + R(x), (1)
with expectations
E[S(x)] = 0, E[R(x)] = 0 and E[Z(x)] = m(x). (2)
The first necessary step in modelling is to establish an attribute S(x) that is stationary over the domain considered. The trend m(x) needs to be removed prior to variogram modelling and geostatistical simulation (Gringarten and Deutsch, 2001). The properties of the stochastic residual directly result from the estimation of a trend model.
In practice, however, it is usually not possible to unambiguously identify and separate the smoothly varying trend from the more erratic residual using the available data (Cressie, 1991). Often, the analyst needs to settle for a plausible, but admittedly non-unique decomposition of the observed spatial behaviour (Zimmerman and Stein, 2010; Journel and Rossi, 1989). Leuanthong and Deutsch (2004) concluded that the inference of the trend model is critical due to its influence on the statistics of the residuals. A misfit of the trend model could result in a severe bias during the assessment of uncertainty.
A variety of trend modelling approaches are popular in practice, many mainly as a result of their ease of application. Data in geological sections can be manually or automatically contoured (Rossi and Deutsch, 2014). Trend surface analysis, a form of multiple regression, fits polynomials by least squares to the spatial coordinates (Knotters, Brus and Oude Voshaar, 1995; Webster and Oliver, 2007). In densely sampled areas, moving window averages can be calculated (Benndorf and Menz, 2014). A 3D trend model can be constructed by rescaling and combining a 2D areal and 1D vertical trend (Deutsch, 2002).
Since the trend and residual are inherently connected, one component should not be modelled without consideration of the other. Common robust estimation algorithms such as universal kriging (UK) or dual kriging can be applied to model both components simultaneously (Isaaks and Srivastava, 1989; Goovaerts, 1997; Webster and Oliver, 2007). However, these algorithms can be used only in cases where the underlying covariance function of the residuals is known
a priori. Reliable inference of these statistics is an issue.
Armstrong (1984) points out that this is merely a chicken-and-egg problem. First, the underlying covariance model of the residuals is needed to solve the kriging equations. Secondly, in order to calculate the residuals, a trend model should be available. Back to square one, since the trend model results from the solution of the kriging equations, which requires a covariance model (Figure 3).
One might think the solution lies in an iterative adjustment of the model parameters until the calculated empirical covariances of the residuals matches those of the model. However, research has shown that such an approach results in a biased characterisation of the spatial variability (Chauvet and Galli, 1982; Cressie, 1991; Lark, Cullis and Welham, 2006).
A first solution to the problem of calculating a reliable variogram consists of selecting data pairs that are unaffected by the trend. This assumes that the pattern of variation of the residuals can be modelled isotropically (Journel and Huijbregts, 1978; Goovaerts, 1997). A more complicated technique consists of filtering out the particular trend model assumed by linear combinations of available data. These higher-order differences are not readily available in the case of non-gridded data and typically result in very large nugget effects (Delfiner, 1976).
Following section presents an alternative to the otherwise ambiguous decomposition of the spatial behaviour into a deterministic trend and a stochastic residual. An objective measure is presented to evaluate the trend model, avoiding a misfit that could result in a severe bias during the assessment of uncertainty.
A NEW APPROACH FOR FITTING
GEOSTATISTICAL MODEL PARAMETERS IN
THE PRESENCE OF A TREND
This section reviews the new approach for fitting the model parameters (Benndorf and Menz, 2014). The principle is founded upon the comparison between the average empirical and theoretical differences of two methods of prediction. The followings statement expresses the basic idea behind the method:
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 theoretical and empirical error will occur.
For a good performance in fitting model parameters, it can be shown that the differences between the two methods of prediction should be potentially large (Menz, 2012). Using UK as prediction method one and the estimation of the trend component based on a leased-square-method (LSM) as second method of prediction will lead to sufficiently large differences. The derivation of the method is briefly outlined. The Interested reader is referred to Benndorf and Menz (2014).
Differences between two prediction methods
based on dual kriging
Having chosen UK as prediction method one and the trend estimation utilising the least-square method as prediction method two, this section evaluates the differences Δ Method 1 – Method 2:
FIG 3 – In order to solve the dual kriging equations a covariance model
is required. However, with the conventional techniques, the parameters
can only be inferred from the residuals which are unknown a priori.
( ) ( ) ( ) ( ) ( ) ( )
Z xUK mx LSM S x mx m x Sx
= - = + - =
D t t t t t t . (3)
To estimate the difference, a noise-eliminated signal needs to be derived. Dual kriging offers the solution (Matheron, 1970; Dubrule, 1983; Cressie, 1991). It represents an alternative formulation to UK and allows separate prediction of a noise-eliminated signal Ŝ(x) (without the noise component R(x)) and the trend component mt(x):
( ) ( ) ( )
Zt xDK=St x +mt x. (4)
The dual kriging equations can be derived from UK (eg Goovaerts, 1997) and are written as:
( ) ( ) ( ) ( )
Zt x0 DK=c k f x0Tt+ 0Tnt=St x0 +mt x0, (5) Where c0 is a covariance vector describing the relationship between the current point of interest x0 and the data points in the selected neighbourhood, kt represent a collection of dual kriging weights, f(x0) contains the evaluated trend functions at x0 and nt is a vector of trend parameters. Equation 5 is identical to the UK estimator and is composed of two separate parts, the signal and the trend components.
After some algebraic manipulations the trend parameters can be obtained via:
( )
TF C z xT 1
=
nt - , (6)
where:
C is a covariance matrix of the sampling locations
T represents the covariance matrix of the trend
parameters: T = FTC-1F-1
z(x) is a vector containing the sampling data
Once the trend parameters are calculated, the empirical difference Δi,empirical between the two methods of prediction can be assessed using dual kriging to estimate the signal Ŝ(x0) at the prediction locations:
st^x0h=c k c0Tt= 0TC-16z^xh- nFt@. (7) The matrix F is the so-called design matrix and contains evaluated trend functions at the sampling locations. The application of geostatistical prediction or simulation requires a regular grid containing m grid nodes which cover the study area. Estimating the signal at each location x0, covering the whole grid, m empirical differences Δi,empirical (prediction method 1 – prediction method 2) can be calculated. The average empirical error is obtained by:
fR, m1 im 1s xi2
emp = =
D^ h
/
t^ h . (8)The parameter R in Equation 8 describes the radius of the ring of influence according to Figure 4 and is part of the fitting process using the error curves. This parameter is user-defined and will be discussed in the following section.
Now a formula for the theoretical error of the differences between both prediction methods needs to be derived. The application of the law of error propagation to Equation 7 results in the definition of a theoretical variance for estimating the signal:
D2^ hx0 =cT0C-1c c0 0- TC FTFC-1 -1c0. (9)
The theoretical error of differences finally results from the squared average signal variance calculated at the m grid nodes. The formula for the theoretical error can be derived analogue to Equation 8:
fR, m1 im 1D x2 i
theo = =
D ^
^ h
/
h. (10)Before proceeding, a quick summary of the previous derived relations is presented.
The application of least-square collocation at each of the m grid nodes xi results in:
st^x0h=c k cT0t= 0TC-16z^xh- nFt@ signal and
D2^ hx0 =cT0C-1c c0 0- TC FTFC-1 -1c0 variance of the signal The average empirical and theoretical errors of differences are obtained by:
fR, m1 im 1s xi2
emp = =
D^ h
/
t^ h empirical error andfR, m1 im 1D x2 i
theo = =
D ^
^ h
/
h theoretical errorError-diagram of differences for fitting model
parameters
In order to fit the model parameters, theoretical and empirical errors of the differences between two methods of prediction have to be compared. To be able to separate the influence of the model parameters on the derived errors, the method can be combined with the use of error curves. For the prediction of attributes at the m grid nodes this approach considers exclusively data in a defined ring of influence centred around each grid node xi. Figure 4 illustrates the concept of the ring of influence.
The parameter R describes the radius of the ring of influence, which is used to determine the squared empirical ŝ(xi)2 and
theoretical D2(x
i) differences for a collection of grid nodes.
For an incrementally increased radius R the corresponding average empirical (represented by circles in Figure 5) and theoretical (represented by crosses in Figure 5) errors can be plotted as a function of the radius of the ring of influence and connected to error curves of differences.
FIG 4 – The concept of the ring of influence. The empirical and
theoretical differences are calculated with data only located in
a particular ring of influence (Benndorf and Menz, 2014).
The example in Figure 5 illustrates the error curves for a case of ‘good model fit’. To achieve a good model fit, the more sensible theoretical error curves have to be adjusted to match the empirical error curves using an iterative back-fitting approach. This is done by systematically changing the model parameters of the variogram. Horizontal deviations are mainly influenced by the range of the variogram; vertical deviations are controlled using the variances. This provides the mean for a systematic fit of model parameters, where the influence of different model parameters is explicitly discriminated.
The universal kriging variance for defining a
best trend model
Above discussed error curves of differences between two prediction methods allow fitting the best covariance model for a given type of trend function. The determination of the best of a few selected trend-cases can be performed by comparing
the UK variance related to several trend models, conditioned to the availability of a sufficient amount of data. Figure 6 plots the obtained standard deviation of the UK variance as function of the ring of influence. A comparison of these plots for different types of trend functions (constant, linear, bilinear or quadratic) will give insight in their performance. The one with the lowest theoretical error provides the best type of trend function for the data.
Summary of fitting geostatistical model
parameters based on dual kriging in the
presence of a trend
Figure 7 summarises the new proposed algorithm for fitting geostatistical model parameters based on dual kriging in the presence of a trend. The algorithm can be defined in seven steps:
1. Prepare the data set for geostatistical analysis. If data are non-Gaussian, perform an N-score transformation. This
FIG 5 – Assessment of theoretical and empirical error curves. The error results from the corresponding average squared differences
calculated for incrementally increasing rings of influence at all the grid nodes (Benndorf and Menz, 2014).
step is necessary since the method assumes a normal distribution of errors of the differences between two methods of prediction.
2. Define an initial covariance model inclusive model parameters. This can be achieved, for example, by calculating the experimental covariance function. Set the type of the trend function model initially to constant (moving window).
3. Calculate the empirical and theoretical error curves of differences between two methods of prediction based on dual kriging (Equations 8 and 10).
4. Evaluate the fit of the error curves. If curves do not fit, adjust the model parameters and go back to step 3. If curves show a nice fit, fix model parameters for the considered type of trend function and go to step 5. 5. Calculate the UK variance for the chosen model parameters
and the trend function. Choose another trend function and go back to step 3, use the initial covariance model parameters. (Optional: to evaluate the performance of a different type of covariance function, go back to step 2.) 6. After all types of trend functions (and optionally
covariance functions) are considered, choose the one with the lowest UK variance.
7. Apply the found trend function and model parameters to the subsequent geostatistical prediction or simulation. Note that the algorithm does not call for an initial trend reduction of the data values z(x). The fitting procedure of the empirical and theoretical error curves of the differences between two methods in prediction separates the influences of the trend and the residual.
Therefore the resulting model parameters represent the deterministic and stochastic model components capturing the properties of the considered attribute as found in the data. Based on these model components a realistic and practically relevant assessment of uncertainty of attributes impacting the profitability, extractability and processability of mineral sands is possible.
CASE STUDY 1 – INFERRING MODEL
PARAMETERS IN PRESENCE OF A TREND
The first case study illustrates the strength of the introduced procedure with respect to the inference of the covariance model parameters. A synthetic, but realistic, 2D example deposit was constructed to evaluate whether the proposed procedure is able to correctly split the variable into a trend and residual component.
The spatial variation of the artificial variable is described over an area of 500 m × 100 m on a grid containing 2000 5 × 5 m2 cells. A complex vertical trend was constructed, resembling the observed patterns in the boreholes at south side of the Rhea target (case study II). The trend suggests the presence of two mineralised horizons, one with its centre about 20 m below surface, the other just above the basement (Figure 8, right). Subsequently a spatial random function and some white noise were added. The spatial variation of the random function is described by an isotropic exponential covariance function with a range of 30 m and a sill of 1. The white noise is normally distributed with a variance of 0.3. The resulting synthetic field is displayed at the top of Figure 8. Upon completion, 19 boreholes were drilled in the synthetic field with a spacing of 25 m. In each borehole, 20 × 5 m long composites were logged (Figure 8, bottom).
Notice the unfavourable conditions of the sampling campaign. The drill hole spacing of 25 m prevents the calculation of a reliable variogram in a direction perpendicular to the trend. Due to the sampling design, the critical short scale variability cannot be characterised in the horizontal direction (Figure 9a). The complex trend along the direction of the borehole makes also the downhole semi-variogram uninterpretable. The negative correlation indicative of a trend clearly shows up as a variogram that increases beyond the sill variance (Figure 9b). Assuming the trend is unknown, the conventional approaches of fitting a variogram are not suited to unambiguously infer the variogram model parameters of the underlying residual random function.
The sampling design only allows a characterisation of the short scale variability in the vertical direction, ie the direction of the strong trend. Following results demonstrate that the new approach for fitting the model parameters based on dual kriging is not affected by this strong vertical trend. Remember that the new approach does not require a prior removal of the trend.
Accounting for the geometry of the available data, a set of ellipsoidal ‘rings’ of influence is created. The outer limit is an ellipsoid with a major axis parallel to the horizontal axis of 260 m and a minor axis of about 65 m. The set is obtained by subdividing this outer limit into 15 ellipsoidal rings. Subsequently, an initial educated guess is made regarding the covariance model parameters. The corresponding theoretical and empirical error curves are calculated and the match between both is evaluated. Figure 10a shows error curves calculated with an isotropic covariance model with a range of 45 m and a sill of 1. After fitting the error curves according to the rules previously explained, the model parameters are improved and the covariance structure of the underlying residual random function is revealed. Figure 10b displays matched error curves resulting from a nugget effect of 0.3 and an isotropic exponential model with a range of 30 m and a sill of 1. Approximately the parameters that were used to construct the artificial sand body.
Once the necessary geostatistical model parameters are inferred the trend component is constructed. An algorithm based on the dual kriging paradigm was developed to calculate the deterministic trend component at all the nodes of the specified grid. The same algorithm was used to extract the residual component from the trend at the data locations.
Subsequently a sequential Gaussian simulation algorithm is used to generate several realisations characterising the uncertainty. The deterministic trend is added onto the stochastic realisations to assess the uncertainty of the original variable.
Figure 11 displays some of the results. A virtual borehole was drilled at an unsampled location in the simulated fields. The resulting calculation data is compared with the known values that were used to create the synthetic field. The modelled trend matches well with the known trend function (Figure 11, left). The spatial variability of the simulated residuals is very similar to that of the combined random function and white noise (Figure 11, middle). Finally, adding the trend back onto the simulated residuals, results in a very realistic representation of the uncertainty on the original variable (Figure 11, right). The simulations follow the general transitional behaviour of the attribute, while realistically characterising the short-scale variability.
CASE STUDY II – GEOMETALLURGICAL
MODELLING OF HEAVY MINERAL SANDS
The second case study is merely indented to illustrate that the proposed methodology results in a realistic model. The combination of a complex trend and a residual random function is perfectly suited to address the problem of characterising the transitional behaviour across soft boundaries.The case study is based on publically available data set which was obtained during a co-funded drilling project of the Australian government under the exploration incentive scheme (Image Resources, 2011; Wamex, 2014). The Rhea target
FIG 8 – Synthetic example of the spatial distribution of the artificial attribute (top). Nineteen boreholes drilled in the synthetic field
(bottom). The complex vertical trend mimicking the behaviour observed in the boreholes of the Rhea target (right).
A B
FIG 9 – Empirical semi-variogram functions calculated in the horizontal (A) and vertical (B) directions from the exploration data.
has an area of 1.5 km2 and is located in Western Australia, about 160 km north-north-west of Perth and 25 km east of the coastal community of Cervantes. The project area lies within the Swan Coastal Plain and contains beach deposits which were formed during a Quaternary cyclic regression phase. The beach deposits overlie Mesozoic fluvial deposits which form the basement in many regions along this palaeocoastline.
The Quaternary sediments contain two episodes of mineralisation. The ‘upper-level’ mineralisation originates from a partially preserved strand zone and is located close to or at the surface. The ‘mid-level’ mineralisation refers to a well-preserved and more extensive sheet-like structure located just above the basement. Incised into the Mesozoic clay basement, a third mineralisation zone occurs in a channel. The channel predates the Quaternary marine sediments and contains a large volume of low-grade material.
A total of 401 aircore holes were drilled for a total depth of 11 × 155 m. The drill holes were positioned along parallel lines forming a 20 m × 100 m sampling grid. Drilling stopped when the Mesozoic basement was intersected. Therefore the depth of the holes varies from about 20 m to 100 m. The slime
content, the oversize and the total heavy mineral grade were recorded every 1 m.
The sand body is subdivided into 64 000 blocks of 20 m × 20 m × 1 m, considering the data density and the selectivity of wet mining equipment. Subsequently, the method of fitting the empirical and theoretical error curves is used to correctly split the observed behaviour into a deterministic trend component and a stochastic residual component. Once the model parameters are inferred the trend model is constructed and several realisations representing the spatial uncertainty are generated. The trend model is added back onto the realisations in order to complete the geostatistical modelling exercise. The results corresponding to one single realisation are displayed in Figure 12.
The cross-section through the west side of the realisation of the slime content clearly indicates an increase of slime content just below the surface and above the basement. Also interesting to notice is that simulations pick up the increase in slime content at the east bank in the south part of the channel (section B-B’, Figure 12a). The general trend of a larger amount of oversized particles between the middle sheet and
FIG 10 – Inference of the covariance model parameter. Error curves not matched in the horizontal directions due to an incorrect range
of 45 m and a missing nugget (A) good match, obtained after range correction to 30 m and addition of a nugget of 0.3 (B).
A
B
FIG 11 – Virtual borehole through simulations at unsampled location. Left: the trend model (dots) matches the known function (crosses). Middle:
the spatial variation of the simulated residuals (grey) is similar to that of the known residual random function (dots). Right: adding the trend onto
the simulated residuals results in a very realistic representation of the uncertainty (grey). The known variable is plotted as a reference (dots).
the upper level mineralisation is reproduced in the generated realisation (Figure 12b). Notice that the increase in oversize in section B-B’ corresponds with a decrease in slime content (Figures 12a and 12b). The mid-level mineralisation shows up in the realisation as a continuous well-preserved sheet, just as described by the geologist (Figure 12c). The realisation also displays the remnants of parallel elongated ellipsoidal structures (Figure 12c). The outline of these structures are however not clearly observed since they were only partially
preserved. Although the channel is generally considered to be of low-grade, the realisations reveal that locally high-grade zones can occur. In the west part of the deposit, the upper strand zone seems to merge with the middle sheet around the second quarter of section A-A’. This event more or less coincides with a decrease in oversize.
This second case study illustrates that the proposed methodology can be successfully used without the conventional, rather arbitrary, delineation of different
FIG 12 – Geometallurgical modelling of the Rhea target. Single realisation of the (A) slime content (per cent); (B) oversize (per cent); and (C)
total heavy mineral content (per cent). Block size is 20 m × 20 m × 1 m. All block models are displayed vertical exaggeration of 5.
A
B
geological zones. The complex features, observed in the boreholes of this particular mineral sand deposit, are reproduced throughout the entire sand body in a realistic manner.
CONCLUSIONS
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 facts’ is incomplete. Ongoing research has shown that a great deal of risk can be mitigated by considering a probabilistic resource model instead of a single deterministic one. Often probabilistic resource models are mainly focused on the characterisation of uncertainty in grades. This paper presents a coherent framework with geostatistical modelling techniques to assess the geological uncertainty and spatial variability of all factors impacting a heavy mineral sands operation. The framework integrates profitability through grade characterisation with extractability and processability through the modelling of slime content and oversize.
This geometallurgical modelling of a mineral sand deposit presents two main challenges. Firstly, the exact location of the geological boundaries are obscure due to the complex interaction of the deposit formation processes. Secondly, the properties of interest exhibit a rather complex transitional behaviour across the so-called soft boundaries. The complexity of the geological processes directly affects the behaviour of the attribute and renders the conventional zonation based on a single cut-off grade highly unreliable and subjective. In order to overcome the existing limits, this contribution presents an alternative simulation methodology which assesses the uncertainty in the entire sand body through a single domain while accounting for the complex trends present.
Due to the direct impact of the trend on the uncertainty assessment, the importance of a correct decomposition of the data into a deterministic trend and a stochastic residual cannot be overstated. A misfit of the trend model could result in a severe biased characterisation of uncertainty. Since both components are inherently connected, one should not be modelled without consideration of the other. This contribution presents a new procedure to unambiguously split the behaviour of the variable under study into a trend and residual component. The practicality of the proposed procedure resides in a straightforward and objective inference of the model parameters which does not require a prior removal of the complex trend. The method extends known cross-validation approaches by two aspects. Firstly, the quality of the fit can be evaluated at grid nodes as well. Secondly, the error curves are ideally suited to separately investigate the influence of different parameters. Additionally the approach ensures a good fit of the complete variogram structure along the different lag distances.
Both case studies demonstrate the applicability and validity of the presented simulation methodology. The added value and strength of the error curves approach is illustrated through an unfavourably sampled known field. Where the conventional semi-variogram yields uninterpretable results, the new fitting technique discovers the model parameters of the underlying residual random function. The second case study was intended to test if the combination of a complex trend and a simulated residual field would results in a geological realistic image of the deposit. The simulations are capable of reproducing complex features described by the field geologist.
The current research was performed under the assumption that only the mean of the spatial random variables is non-stationary and that the residuals are stochastically independent of the trend (homoscedasticity). Solutions exist to independently deal with non-stationary two-point statistics (Machuca-Mory and Deutsch, 2013) and effects of heteroscedasticity (Leuangthong and Deutsch, 2004). However these solutions already (implicitly) assume that the trend has been correctly and unambiguously removed from the data. Future research needs to focus on the integration of a correct decomposition into a trend and a stochastic residual, based on model parameters which are correctly inferred under conditions of non-stationary one-point and two-point statistics. Subsequently spatial variability can be characterised through non-stationary residual field simulations.
ACKNOWLEDGEMENTS
The findings of this investigation are consequent of a research and development project sponsored by IHC Merwede whose assistance and support is gratefully acknowledged.
REFERENCES
Armstrong, M, 1984. Problems with universal kriging, Mathematical Geology, 16(1):101–108.
Benndorf, J and 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 kriging in the presence of a trend, Stochastic Environmental Research and Risk Assessment, 28(3):627–637.
Chauvet, P and Galli, A, 1982. Universal kriging, CGMM internal report number C-96.
Cressie, N A C, 1991. Statistics of Spatial Data, 900 p (John Wiley and Sons Inc: New York).
Daedal Research, 2014. Global mineral sand industry: trends and opportunities (2013–2018) [online]. Available from: <http:// www.marketresearch.com/Daedal-Research-v3440/Global-Mineral-Sand-Trends-Opportunities-8031855/> [Accessed: 8 May 2014].
Delfiner, P, 1976. Linear estimation of nonstationary phenomena, in Advanced Geostatistics in the Mining Industry (eds: M Guarascio, M David and C J Huijbregts), pp 49–68 (Reidel Publishing: Dordrecht).
Deutsch, C V, 2002. Geostatistical Reservoir Modeling, 376 p (Oxford University Press: New York).
Dominy, S C, Noppé, M A and Annels, A E, 2002. Errors and uncertainty in mineral resource and ore reserve estimation: the importance of getting it right, Journal of Exploration and Mining Geology, 11(4):77–98.
Dubrule, O, 1983. Two methods with different objectives: splines and kriging, Mathematical Geology, 15(2):245–257.
Ferrara, G, Guarascio, M and Massacci, P, 1993. Mineral processing design and simulation using geostatistical modelling of ore characteristics, in Proceedings XVIII International Mineral Processing Congress 1993, pp 497–504 (International Mineral and Processing Council: Sydney).
Goovaerts, P, 1997. Geostatistics for Natural Resource Evaluation, 483 p (Oxford University Press: New York).
Gringarten, E and Deutsch, C V, 2001. Teacher’s aide: variogram interpretation and modeling, Mathematical Geology, 33(4):507–534. Image Resources, 2011. EIS co-funded drilling, Cooljarloo Project
E70/2898 North Perth Basin, Western Australia, final report, Department of Mines and Petroleum.
Isaaks, E H and Srivastava, R M, 1989. An Introduction to Applied Geostatistics, 561 p (Oxford University Press: New York).
Jones, G, 2009. Mineral sands: an overview of the industry [online]. Available from: <http://www.iluka.com/docs/company- presentations/mineral-sands---an-overview-of-the-industry-by-greg-jones-manager-development-geology> [Accessed: 29 May 2014].
Journel, A G and Huijbregts, C J, 1978. Mining Geostatistics, 600 p (The Blackburn Press: Cadwell).
Journel, A G and Rossi, M E, 1989. When do we need a trend model in kriging?, Mathematical Geology, 21(7):715–739.
Knotters, M, Brus, D J and Oude Voshaar, J H, 1995. A comparison of kriging, co-kriging and kriging combined with regression for spatial interpolation of horizon depth with censored observations, Geoderma, 67(3):227–246.
Lark, R M, Cullis, B R and Welham, S J, 2006. On spatial prediction of soil properties in the presence of a spatial trend: the empirical best linear unbiased predictor (E-BLUP) with REML, European Journal of Soil Science, 57(6)787–799.
Leuangthong, O and Deutsch, C V, 2004. Transformation of residuals to avoid artifacts in geostatistical modelling with a trend, Mathematical Geology, 36(3):287–305.
Machuca-Mory, D F and Deutsch, C V, 2013. Non-stationary geostatistical modeling based on distance weighted statistics and distributions, Mathematical Geoscience, 45(1):31–48.
Matheron, G, 1970. La théorie des variables régionalisées et ses applications, Cahiers du centre de morphologie mathématique de Fontainebleau Fasc no 5, Ecole des Mines des Paris. Translation (1971): The theory of regionalised variables and its applications.
McKinsey & Company, 2012. Urban world: cities and the rise of the consuming class [online]. Available from: <http://www. mckinsey.com/insights/urbanization/urban_world_cities_and_ the_rise_of_the_consuming_class> [Accessed: 29 May 2014]. Menz, J, 1999. Forschungsergebnisse zur Geomodellierung und ihre
Bedeutung, in Fortschritte der Geoinformatik. Tagungsheft, Mathematische Geologie 4 (ed: H Thiergärtner), pp 19–30, Geo-Berlin’98, Symposium Geoinformatik und Modellierung (TU Berlin: Dresden).
Menz, J, 2012. Signalgesteuerte Kollokation – ein Verfahren zur Ableitung der Modellparameter für die geostatistische Vorhersage und Simulation, in Proceedings Sammelband 13 Geokinematischer Tag, pp 333–342 (Instituts für Markscheidewesen und Geodäsie and der TU Bergakademie: Freiberg).
Rossi, M E and Deutsch, C V, 2014. Mineral Resource Estimation, 332 p (Springer: Dordrecht).
Vallee, M, 2000. Mineral resource + engineering, economic and legal feasibility = ore reserve, CIM Bulletin, 93:53–61.
Wamex, 2014. Western Australian mineral exploration index [online]. Available from: <http://www.dmp.wa.gov.au/5136.aspx> [Accessed: November 2013].
Webster, R and Oliver, M A, 2007. Geostatistics for Environmental Scientists, 315 p (John Wiley & Sons Ltd: Chichester).
Zimmerman, D L and Stein, M, 2010. Classical geostatistical methods, in Handbook of Spatial Statistics (eds: A E Gelfand, P J Diggle, M Fuentes and P Guttorp), pp 29–44 (CRC Press: New York).
