MODERN MINING GEOSTATISTICS -
CONVENTION WITH CODELCO YEAR 10 (2015) -
ANNIVERSARY REPORT

Serge Antoine Séguret, Francky Fouedjio, Ramon Freire, Sebastian de La Fuente, Cristian Guajardo, Claudio Rojas

<figure><img src="image_1.png" /></figure>

To cite this version:

Serge Antoine Séguret, Francky Fouedjio, Ramon Freire, Sebastian de La Fuente, Cristian Guajardo, et al.. MODERN MINING GEOSTATISTICS - CONVENTION WITH CODELCO YEAR 10 (2015) - ANNIVERSARY REPORT. [Research Report] R151215SSEG, Mines ParisTech; Center for geo-sciences. 2016, 117 p. hal-01244786v2

HAL Id: hal-01244786

https://minesparis-psl.hal.science/hal-01244786v2

Submitted on 22 Jun 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
<figure><img src="image_2.png" /></figure>

MODERN MINING GEOSTATISTICS

CONVENTION WITH CODELCO YEAR 10 (2015)

ANNIVERSARY REPORT

(June 2016 version)

For the Ecole des Mines :

Serge Antoine Séguret, Francky Fouedjio

For CODELCO :

Ramon Freire, Sebastian de la Fuente, Cristian Guajardo, Claudio Rojas

December 2015

N° R151215SSEG

Mines ParisTech

Centre of Geosciences and Geoengineering

35 rue Saint Honoré

77305 Fontainebleau - France
Serge Antoine Séguret

Global report of year 10 (2015). Contrat ARMINES/CODELCO Chile no. 10701

"Geostatistical consulting, tutoring, development and implementation contract".

Décembre 2015. 117 pages.

<table><tr><td>Equipe</td><td>Géostatistique</td></tr><tr><td>Visa</td><td>Jacques Rivoirard</td></tr></table>
Introduction

During the 10th year of the convention between CODELCO and the Ecole des Mines of Paris, effort has been made to promote the work done previously by attending four congresses in three countries (Chile, France, Germany), two presentations being made by engineers of Codelco, namely Cristian Guajardo (Chuquicamata) and Sebastian de la Fuente (Radomiro Tomic).

In comparison to previous years, this anniversary report is presented as a catalog containing
the papers and the relevant oral presentations with comments. For Mine-Planning and
Geomin, no paper was published.

This report is not confidential and can be distributed in the good will of Codelco managers
and engineers.

Content

A Geostatistical comparison between blast and drill holes in a porphyry copper deposit

→page 5

B The oral presentation with comments of “Geostatistical comparison between blast and drill holes in a porphyry copper deposit”

→page 15

C The oral presentation with comments of “Diamond Drill Holes, Blast Holes & Cokriging”

→page 35

D The oral presentation with comments of “Anisotropy of the Rock Quality Designation (RQD) & its Geostatistical Evaluation”

→page 55

E Geostatistical Evaluation of Rock-Quality Designation and its link with Linear Fracture Frequency

→page 77

F The oral presentation with comments of “Geostatistical Evaluation of Rock-Quality Designation and its link with Linear Fracture Frequency”

→page 87

G Breccia Pipe Prediction: a new approach using non-stationary covariance

→page 107

Congresses involved

WCSB 7

7th World Conference on
Sampling and Blending

10-12 june 2015
Bordeaux, France

Centre de Congrès Cité Mondiale
18 parvis des Chartons - 33080 Bordeaux

GEOMIN 2015

6th International Seminar on
Geology for the Mining Industry

<figure><img src="image_3.png" /></figure>

Universitat Politècnica de Catalunya

GECAMIN
Consejo Nacional de Investigación Científica

MINEPLANNING
2015

4th International Seminar on
Mine Planning

July 8 - 10, 2015

Enjoy Hotel, Antofagasta, Chile

<figure><img src="image_4.png" /></figure>

IAMG 2015

The 17th annual conference of the International Association for
Mathematical Geosciences,
September 5-13, 2015, Freiberg (Saxony) Germany

<figure><img src="image_5.png" /></figure>

MINES PARIS

Convenio

<figure><img src="image_6.png" /></figure>

CODELCO
Produciendo Futuro
# Chapter A

## Geostatistical comparison between blast and drill holes in a porphyry copper deposit

(June 2016 version)

Serge A. Séguret (Mines ParisTech, France)

A paper presented at WCSB7, 7th world conference on sampling and blending, 10-12 June, Bordeaux, France

### Abstract

Diamond drill-hole grades are known to be of better quality than those of blast holes; is this true? We present a formal study of a porphyry copper deposit in Chile where the variogram of 3 meter long drill hole samples is compared to 15 meter long blast hole ones and we show that the blast holes can be assumed to regularizing the point information deduced from the drill holes, except for a nugget effect specific to the blast samples. Complementary analyses based on migrated data show that the drill holes also have their own errors.

After a brief description of the first steps in the blast sampling protocol, we show, by using extension variance concepts, that the blast error is not due to the arbitrary removal of material from the sampling cone produced by drilling.

The present study establishes a formal link between blast and drill holes which leads to linear systems:

* Removal by kriging of the blast (or the drill) error;

* Deconvolution of the blast measurements to transform them into point ones;

* Block modeling where drill and blast holes are used together.

In the paper, we thought it useful to detail some calculations and give some key formulas so that the reader can eventually adapt to other comparisons such as diamond drill holes compared to reverse circulation drill holes. Overall, this study shows how to combine measurements known on two different supports, a very complex challenge.
# Geostatistical comparison between blast and drill holes in a porphyry copper deposit

(June 2016 version)

Serge Antoine Séguret

MINES ParisTech, Center for Geosciences/Geostatistical team, 35 rue Saint Honoré, 7730 Fontainebleau, France. E-mail:

serge.seguret@mines-paristech.fr

Serge has worked for more than thirty years in the famous laboratory founded by Georges Matheron in 1968 at Fontainebleau, where the major developments of Geostatistics have been achieved. For the last ten years, he has been heavily involved in mining operations via a long term collaboration with Codelco (Chile) and Vale (Brazil).

## Abstract

Diamond drill-hole grades are known to be of better quality than those of blast holes; is this true? We present a formal study of a porphyry copper deposit in Chile where the variogram of 3 meter long drill hole samples is compared to 15 meter long blast hole ones and we show that the blast holes can be assumed to regularizing the point information deduced from the drill holes, except for a nugget effect specific to the blast samples. Complementary analyses based on migrated data show that the drill holes also have their own errors.

After a brief description of the first steps in the blast sampling protocol, we show, by using extension variance concepts, that the blast error is not due to the arbitrary removal of material from the sampling cone produced by drilling.

The present study establishes a formal link between blast and drill holes which leads to linear systems:

* Removal by kriging of the blast (or the drill) error;

* Deconvolution of the blast measurements to transform them into point ones;

* Block modeling where drill and blast holes are used together.

In the following, we thought it useful to detail some calculations and give some key formulas so that the reader can eventually adapt to other comparisons such as diamond drill holes compared to reverse circulation drill holes. Overall, this study shows how to combine measurements known on two different supports, a very complex challenge.

## Data

The data comes from an open-pit copper mine in Northern Chile of which a 600x400x125 m³ sub domain is analysed (Figure 1) as it is almost homogeneously covered by around 3,000 drill-hole samples (3m long) and 13,000 blast-hole samples (15m long).

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 2.** (a) Drill hole copper grade variogram; (b) Blast copper grade variogram. Three directions are represented, 45° North (N45), 135° North (N135), and vertical (D-90). Black continuous line is the isotropic variogram</figcaption></figure>

## Methodology

The geostatistical comparison between the two types of measurements is decomposed into two steps:

### 1. Deconvolution & Convolution:

* Starting from the drill variogram, identifying the basic structures that model its behavior and deducing the underlying "point" variogram by deconvolution;

* Making the theoretical convolution of the point variogram on 15-meter long supports and checking that it correctly fits the vertical and horizontal blast variograms, except for an additional nugget effect of 0.2.

### 2. Migration & Cross variogram

* As there is no point where both drill and blast measurements are known, we make some blast holes migrate to drill hole locations and calculate the cross variogram;

* The objective is to measure the nugget effect shared by the two types of measurements.

There are not enough drill samples to distinguish between horizontal and vertical drill variograms (they are drilled along many different directions). This is the first reason why an omnidirectional variogram will be considered for the drill samples, the second one is that all the formulas at our disposal require isotropy.

Consequently, we make two comparisons between:

* An omnidirectional drill variogram and a vertical blast one;

* An omnidirectional drill variogram and a horizontal blast one.

The distinction is important because the formulas differ between the two cases.

## General formulas

All the formulas have been known for a long time in the literature, but in different places, and some are not even published. For the convolution charts, the most useful reference is probably ¹; for the complete fundamental formulas, refer to ². Concerning the extension formulas, refer to ³.

In the following we apply a procedure illustrated in ⁴ where we use the following approximation of a variogram regularized over a
support “l” (the distance “h” being large in comparison with the dimension of the support):

$$ \gamma_l(h) = \gamma(h) - \bar{\gamma}(l, l) \quad (1) $$

with

$$ \bar{\gamma}(l, l) = \frac{1}{l^2} \int_0^l \int_0^l \gamma(u - v) du dv \quad (2) $$

$\bar{\gamma}(l, l)$ is the average of the point variogram when both extremities of vector h describe the support independently. In (2), 1D integrals are used because the core diameters are small compared to the lengths.

The way this formula is applied depends on the structure of the point variogram (spherical, exponential, linear, etc...) but also on the calculation direction compared with the regularization direction. In the following, we consider two situations:

* The calculation direction is parallel to the regularization direction, notation $\gamma_l''(h)$;

* The calculation direction is perpendicular to the regularization direction, notation $\gamma_l^\perp(h)$.

For the structures with a range, whether asymptotically (Exponential, Gaussian) or real (Spherical), we have:

$$ \text{range of } \gamma_l(h) = \text{range of } \gamma(h) + l \quad (3) $$

Note that (3) is not compatible with approximation (1) which amounts to assigning to the regularized model the same range as that of the point model. So (1) is essentially useful for <u>comparing the sills</u> of regularized structures.

## Step 1: deconvolution & convolution

### Fitting the drill-hole variogram

Three basic structures are necessary: nugget, exponential, linear:

$$ \gamma_{\text{drill}}(h) = \varepsilon_{\text{drill}}^2 + C_{\text{drill}}(1 - e^{-a_{\text{drill}}}) + b_{\text{drill}} |h| \quad (4) $$

with:

$$ \varepsilon_{\text{drill}}^2 = 0.05, C_{\text{drill}} = 0.05, 3a_{\text{drill}} = 35m, b_{\text{drill}} = \frac{0.015}{100} $$

<figure><img src="image_9.png" /><figcaption>**Figure 3.** Drill hole variogram fitting. Dotted line, the experimental curve; continuous line, the model</figcaption></figure>
## Nugget effect (or small-range structure) deconvolution & convolution

The attenuation of the nugget effect, whether “pure” or associated with a microstructure which reaches its sill long before the first variogram lags, is proportional to the ratio of the supports. In the present case study, the diameters of drill holes and blast holes are considered to be equal and we ratio the lengths but generally speaking, one has to consider the ratio of the volumes:

$$ \varepsilon_{\text{blast}}^2 = \frac{l_{\text{drill}}}{l_{\text{blast}}} \varepsilon_{\text{drill}}^2 \qquad (5) $$

With $l_{drill}=3$, $l_{blast}=15$, the nugget effect of the blasts must be five times smaller than that of the drills. For $\varepsilon_{drill}^2 = 0.05$ (Figure 2a) we obtain $\varepsilon_{blast}^2 = 0.01$, a value three times smaller than the 0.03 value deduced from the blast variogram (Figure 2b). If one takes the blast nugget effect as a reference, the drill nugget effect should be 0.15, a quantity above the local sill of the variogram and not realistic.

Conclusion: the support cannot explain the differences between the nuggets of the blasts and of the drills. The blast nugget is too large.

## Vertical variograms – Deconvolution & convolution

The calculation direction is parallel to the blast regularization direction (i.e. vertical).

### Exponential structure

If the practical drill range is 35m, the parameter associated with the underlying point exponential structure is expressed by (3):

$$ 3a_0 = 35 - 3 \rightarrow a_0 = 10.7 $$

If $\gamma(\cdot)$ denotes a variogram normalized by its sill, the underlying point sill $C_0$ of the exponential structure is produced by (1):

$$ C_{\text{drill}} = C_0(1 - \bar{\gamma}(3,3)) $$

For the exponential structure, the charts in ¹ yield:

$$ \bar{\gamma}(3,3) = 0.087 \rightarrow C_0 = 0.055 $$

For $l=15m$, we deduce:

$$ C_{15} = C_0(1 - \bar{\gamma}(15,15)) $$

and we obtain:

$$ \bar{\gamma}(15,15) = 0.34 \rightarrow C_{15} = 0.036 $$

We will see later if these results correspond to the experimental blast variogram, but we must first look at the linear structure which completes the model (4).

### Linear structure

For $h>l$ we have, where $b$ is the slope of the structure ¹:

$$ \bar{\gamma}(l, l) = b \frac{l}{3} \qquad (6) $$

The slope $b$, which does not change with the support, is given by the drill samples and the difference between two supports $l$ and $l'$ equals:

$$ \bar{\gamma}(l, l) - \bar{\gamma}(l', l') = b \frac{l-l'}{3} \qquad (7) $$
When $l=0$m and $l'=3$m, and with $b_{drill} = \frac{0.015}{100}$ obtained by (4), the attenuation is 0.00015, a negligible quantity. When $l=3$m and $l'=15$m, the attenuation is 0.0006, still negligible. In any case, the effect of the regularization on the linear structure is negligible. This is due to the weak slope of the linear structure.

The combination of all the regularizations is shown in figure 4a where the dotted line represents the actual model and the red line the model we should obtain with a more realistic nugget effect. One can see that apart from the problem of the nugget effect, the variation range is acceptable, even if the linear part of the theoretical structure does not appear in the vertical experimental blast variogram.

<figure><img src="image_10.png" /><figcaption>**Figure 4.** In blue, the points of the experimental blast variogram; dotted line, the theoretical model for the blasts deduced from the drills; in red, the theoretical model with a more realistic nugget effect (a) Theoretical regularization parallel to the vertical blast variogram (b) Theoretical regularization perpendicular to the horizontal blast variogram.</figcaption></figure>

## Horizontal variograms – Deconvolution & convolution

The calculation direction is perpendicular to the blast regularization direction (i.e. horizontal).

The same procedure is followed, the only difference is that approximation (1) is not acceptable and we have to use charts that produce the exact calculation (see ¹, chart number 11).

We obtain figure 4b where the dotted line represents the actual model and the red line the model we should obtain with a more realistic nugget effect. The fit is good.

## First conclusions

If we omit the problem of the nugget effect, we see that both blast and drill holes can be considered as a regularization of the same reality according to their respective supports. This result, which we did not dare to hope, surprised us pleasantly and shows that the measurements from the blast holes are not as bad as people often think, anyway the case for this company. But the approach followed up to now suffers from two uncertainties:

* The analyses are done independently. Imagine that all the blast locations have been shifted from a constant equal to the range (around 100m). In that case, the correlation between blast and drills will be zero while the same coherence properties are maintained when making individual regularizations as previously;

* The analyses refer to the drill nugget assumed to be a "natural" micro structure; is this true?

To answer these questions, cross variograms must be calculated but we do not have any location with both measurements, so a migration is necessary.

## Step 2: migration & cross variogram
## Migration

In order to obtain a significant number of measurements at the same location, around 1,000 blasts samples were migrated to drill locations when the migration distance did not exceed 10 meters. Figure 5a presents the scatter diagram between the migrated values and the drill ones. The correlation coefficient is low (0.4) because the nugget effects are large.

<figure><img src="image_11.png" /><figcaption>**Figure 5.** (a) Scatter diagram between migrated blasts and drills; (b) Direct variogram of migrated blasts (black triangles), corresponding drills (red points) and cross-variogram of both (blue stars). The cross variogram reveals a tiny negative nugget effect with no comparison with the drill or blast ones</figcaption></figure>

On Figure 5b, points (resp. triangles) present the migrated blast (resp. drill) variograms. They differ slightly from the previous ones because the number of samples is smaller and the migration affects the results. In the same figure, the stars represent the cross variogram which does not show a significant nugget effect, possibly a small negative one without any magnitude in common with the effects encountered on the individual variograms.

## Conclusions

It seems that the drill-holes have their own errors too, independent of the blast ones, and the two measurements share only the structured parts of the variogram: the exponential and linear structures.

## Analyse of the blast error

### Description of the blast sampling

Up to now the theoretical blast support has been set to 15m but in fact the blast drilling length is approximately 17m, producing a large cone from the floor of which around 5cm of material is removed by hand across the entire surface, the idea being to restore an overall volume of 15m. Without any consideration of the numerous sampling procedures, we stay at this stage and ask the question: could the error specific to the blasts be due to the arbitrary removal of material and the blast length variability?

### Randomization of the blast support

Let $l$ and $l'$ be two different supports. One finds in 3 the formula which expresses the variance of the difference between the two grades Y over $l$ and $l'$, called "extension variance from $l$ to $l'$, also equal to twice the variogram between the grades averaged over the two supports:

$$D^2(Y_l(x) - Y_{l'}(x+h)) = E[(Y_l(x) - Y_{l'}(x+h))^2] = 2\gamma_{ll'}(h) = 2\bar{\gamma}(l, l') - \bar{\gamma}(l, l) - \bar{\gamma}(l', l') \quad (8)$$
In (8), $l'_h$ represents the translation of the support $l'$ by a vector $h$. $\bar{\gamma}(l, l)$ and $\bar{\gamma}(l', l')$ represent the averaged variogram when two points move independently along both the supports involved.

Suppose that $l$ and $l'$ are randomly and independently selected uniformly in an interval, for example equal to [12.5m, 17.5m]. Then one has to calculate the mathematical expectation of (8) to obtain the resulting variogram. We have:

$$E[\bar{\gamma}(l, l)] = E[\bar{\gamma}(l', l')] \quad (9)$$

$$E[\gamma_{ll'}(h)] = E[\bar{\gamma}(l, l')] - E[\bar{\gamma}(l', l')] \quad (10)$$

(10) is the theoretical variogram that we want to compare to the actual experimental variogram in order to verify if the blast nugget could be associated with some support-length uncertainty.

$E[\bar{\gamma}(l, l')]$ is a continuous function, complex to calculate as it depends on the mutual configuration of $l$ and $l'$, but about which we know that for $h$ greater than the range plus $l$, it reaches and stays at the sill of the underlying point variogram. In practice, the only structure that we consider is the exponential; its point sill is 0.055. For the interval of support-length uncertainty [12.5m, 17.5m], we deduce from (10) that the sill is reduced by a quantity obtained by:

$$E[\bar{\gamma}(l, l)] = \frac{1}{l} \int_{12.5}^{17.5} \bar{\gamma}(l, l) dl \quad (11)$$

To evaluate the range of variations, the integral (11) is approximated by a finite sum:

$$E[\bar{\gamma}(l, l)] \approx \frac{1}{5} \sum_{l=13,14,15,16,17} \bar{\gamma}(l, l) \quad (12)$$

We use the same charts as previously to calculate the values of $\bar{\gamma}(l, l)$ involved and finally (12) yields:

$$E[\bar{\gamma}(l, l)] \approx 0.055 \frac{1}{5}(0.295 + 0.305 + 0.325 + 0.337 + 0.352) = 0.055 * 0.323 \quad (13)$$

Notice that even if we randomize the blast support over a larger interval still centered around 15m, the variance reduction does not change and stays approximately equal to the sill multiplied by 0.325. If we suppose that the support fluctuation is not symmetric around 15m, but around 13m for example, the multiplicative factor for the sill reduction decreases to 0.295. In any case, we conclude that:

* The uncertainty on the support length does not produce a nugget effect but a variance reduction;

* This variance reduction represents approximately 30% of the underlying variogram sill;

* The arbitrary removal of the material, as well as the uncertainty on the blast length, cannot explain an error specific to the blasts and necessarily linked to the subsequent sampling procedures.

**Summary: a formal link between blast and drill holes**

## Formal link

Finally, we have:

$$Y_{\text{blast}}(x, y, z) = Y(x, y, z) * p_{15m}(z) + R(x, y, z) \quad (14)$$

with

$Y(x,y,z)$, the point grade assumed to be isotropic and devoid of any measurement error; "*" denotes a convolution product;

$$Y(x, y, z) * p_{15m}(z) = \int_{-\infty}^{+\infty} Y(x, y, u) p_{15m}(z - u) du;$$
$$p_{15m}(z) = \frac{1}{15} \mathbf{1}_{[0, \frac{15}{2}]}(|z|);$$

$\mathbf{1}_{[0, \frac{15}{2}]}(|z|)$ the indicator function equal to 0 outside the interval and 1 inside it;

R(x,y,z), a “white noise” residual statistically and spatially independent from Y(x,y,z) and representing the blast error

The variogram of Y_blast(x,y,z) becomes:

$$\gamma_{\text{blast}}(h) = \gamma_{15m}(h) + \gamma_R(h) \qquad (15)$$

with

$\gamma_R(h)$, the nugget effect due to the blast error, having the variance $\sigma_R^2$;

$$\gamma_{15m}(h) = (\gamma * P_{15m})(h) - (\gamma * P_{15m})(0);$$

$\gamma(h)$, the point variogram, assumed to be isotropic;

$$P_{15m}(h) = p * \dot{p}(h) = \frac{1}{15^2}(-|h| + 15)\mathbf{1}_{[0,15]}(|h|)$$

The model supposes that the blasts and the drills have the same average because the independent residuals are of zero mean. It must be verified when using this model. It is approximately the case here (0.63 for the drills, 0.69 for the blasts).

## Removing the blast error by kriging

Model (14) can be used to remove the blast error by “Factorial Kriging” estimation⁵. One can easily build a linear system applicable to each blast measurement, choosing a local neighborhood of surrounding blast samples. The system is presented symbolically by using matrix formalism:

$$\begin{pmatrix} \gamma_{15m} + \gamma_R & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} \lambda \\ \mu \end{pmatrix} = \begin{pmatrix} \gamma_{15m} + \sigma_R^2 \end{pmatrix}$$

In this system, $\gamma_R$ disappears from the second member of the linear system whereby we remove, from the estimation, the part associated with the measurement error. It does not mean that in the remaining part $\gamma * P_{15m}$ there is no nugget effect; it means that only the “natural” part remains. In our case, the complete nugget effect has to be removed because blasts and drills do not share any micro-structure.

## Deconvolution by kriging

It may be interesting to remove the effect of regularization on the blast using a kriging system which estimates, for each blast measurement, a “point” value while simultaneously removing the part of the nugget effect associated with blast errors:

$$\begin{pmatrix} \gamma_{15m} + \gamma_R & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} \lambda \\ \mu \end{pmatrix} = \begin{pmatrix} \gamma * p_{15m} - (\gamma * p_{15m})(0) + \sigma_R^2 \end{pmatrix}$$

The difference with the previous system is that in the second member, $\gamma * P_{15m}$ (capital “P_{15m}”) is replaced by:

$$\gamma * p_{15m}(h_x, h_y, h_z) = \int_{-\infty}^{+\infty} \gamma(h_x, h_y, u) p_{15m}(h_z - u) du \quad (\text{small "p"}).$$

## Block estimate by cokriging drill and blast measures

Finally, one can imagine locally renewing the mine planning block model by using blasts and drills together through a cokriging system with a linked mean (same average for both measurements):

$$\begin{pmatrix} \gamma_{3m} & \gamma_{3m,15m} & 1 \\ \gamma_{3m,15m} & \gamma_{15m} + \gamma_R & 1 \\ 1 & 1 & 0 \end{pmatrix} \begin{pmatrix} \lambda \\ \lambda' \\ \mu \end{pmatrix} = \begin{pmatrix} \gamma * p_{3m} * p_v - (\gamma * p_{3m} * p_v)(0) \\ \gamma * p_{15m} * p_v - (\gamma * p_{15m} * p_v)(0) + \sigma_R^2 \end{pmatrix}$$
These systems were tested on a realistic simulation where the truth is known; they produce good results which will be published in the near future.

## Conclusion

In this deposit – and more generally, in this company (other test have been done), diamond drill hole grades and blast hole grades are consistent in the sense that, apart from the nugget effect, the structured part of their respective variograms follow the theoretical laws of regularization.

Concerning the nugget effects, we discover, by cross-analyses, that there is no natural micro-structure in the underlying point grade and the large nugget effects encountered on the variograms (approximately 50% of the variance for blasts and drills) are due to blast and drill measurement errors, independent of either measurement type.

The analysis of the blast error leads to the conclusion that the error is not due to the first step of the sampling procedure, it has to be found later in the process.

As a conclusion, some linear systems are proposed for removing the nugget effects from the data, reducing the effect of convolution and, more importantly, using blasts and drills together for the short-term mine planning. These systems, among numerous different potential ones, easy to demonstrate, result directly from the formal link established here between blast and drill holes. Before using these systems, the link must be verified by adhering to the methodology presented here.

## Acknowledgements

The author would like to thank CODELCO Chile, for its strong support in the implementation of good geostatistical practices along the copper business value chain. I warmly thank the Editor for his excellent recommendations which improved the manuscript. Many thanks to Sebastian de La Fuentes and the complete team of R&T. Without them, nothing would have been possible.

## References

1. A.G. Journel, Ch. J. Huijbregts, "Mining Geostatistics", 5th Edn. The Blackburn Press (1991).

2. J. Serra, « *Echantillonnage et estimation locale des phénomènes de transition miniers* ». Doctoral thesis, Faculté des Sciences de Nancy (1967).

3. G. Matheron, « *Les Variables Régionalisées et leur estimation* ». Masson et Cie Editeurs (1965).

4. S. A. Séguret, "Spatial sampling effect of laboratory practices in a porphyry copper deposit". *Proceedings of the 5th world conference on Sampling and Blending*, Santiago, Chile (October 2011).

5. J. P. Chilès, P. Delfiner, "Geostatistics", 2nd edition, Willey (2012).
Chapter B

The oral presentation with comments of
“Geostatistical comparison between blast and
drill holes in a porphyry copper deposit”
(June 2016 version)

Serge A. Séguret (Mines ParisTech, France)

Presented at WCSB7, 7th world conference on sampling and blending, 10-12 June, Bordeaux, France
# GEOSTATISTICAL COMPARISON BETWEEN BLAST AND DRILL HOLES IN A PORPHYRY COPPER DEPOSIT

(June 2016 version)

Serge Antoine Séguret – Mines-ParisTech, France

<figure><img src="image_12.png" /></figure>

CODELCO
Produciendo Futuro
Codelco Chili

WCSB7 7th World Conference on
Sampling and Blending

10-12 june 2015
Bordeaux, France
Centre de Congrès Cité Mondial
18 parvis des Chartrons - 33080 Bordeaux

<figure><img src="image_13.png" /></figure>

Mines-ParisTech

Good morning everybody.

I am pleased to present you with the work entitled

"GEOSTATISTICAL COMPARISON BETWEEN BLAST AND DRILL HOLES IN A PORPHYRY COPPER DEPOSIT".

This work is supported by the Chilean company Codelco, which produces copper, and the Paris School of Mines where I have worked for over thirty years in the Geostatistical laboratory founded by Georges Matheron at Fontainebleau.
Drill holes, blast holes

<figure><img src="image_14.png" /></figure>

Drill holes

* Few

* Long-term planning (month, year, decade)

* Large-scale strategy

* Good quality

<figure><img src="image_15.png" /></figure>

Blast holes

* Many

* Short-term planning (day, week,)

* Small-scale selectivity

* Bad quality

2

Typically in open pit mines, geologists, mining engineers, metallurgists, have at their disposal two types of measurements for the grades:

A first type, from drill holes – “diamond drill holes” in our case.

A second type, from the blast holes.
<figure><img src="image_16.png" /><figcaption>Drill holes</figcaption></figure>

## Drill holes, blast holes

<figure><img src="image_17.png" /></figure>

* Few

* Long-term planning (month, year, decade)

* Large-scale strategy

* Good quality

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>Blast holes</figcaption></figure>

<figure><img src="image_20.png" /></figure>

* Many

* Short-term planning (day, week,)

* Small-scale selectivity

* Bad quality

3

Because they are much more expensive, the diamond drill holes are less numerous than the blast holes, and it is usual to encounter sampling rates ranging from one over three to one over ten or worse.

Not only the sampling density is involved but also how samples are distributed in space too, as shown on the slide:

* The circles, representing the drill holes, are widely spaced, see the left horizontal cross-section

* In the same figure, the green crosses, representing the blast holes, are more densely spaced

* Vertically, on the right-hand figure, it is the reverse, with almost continuous drill hole information while the blasts are more widely spaced

These differences make it even more difficult to compare the statistical properties of the two types of measurements, including when calculating directional variograms because statistical inference conditions are not the same
<figure><img src="image_21.png" /></figure>

Another difference concerns the way the measurements are used.

The long road that leads to the opening of the mine is marked by drilling campaigns, to achieve the block model that will condition the exploitation at large scale as well as for medium- and long-term planning. Typically, kriging and Geostatistics are used to build the model at this stage.

In addition, the blast holes are used for short term planning with no need of Geostatistics, a simple moving average is often used to estimate the block quantity of metal

These separate uses of two types of measurements that are supposed to represent the same thing raise questions about their relationship. In particular, would it not be possible to enrich the short-term estimate, now based only on blast holes, by adding the drill hole measurements?
Drill holes, blast holes

<figure><img src="image_22.png" /></figure>

Drill holes

* Few

* Long-term planning (month, year, decade)

* Large scale strategy

* Good quality

?

<figure><img src="image_23.png" /></figure>

Blast holes

* Many

* Short-term planning (day, week,)

* Small-scale selectivity

* Bad quality

5

Finally, we often hear, without real justification, that the diamond drill holes are much better than the blast ones.
We ask the questions:

* Better how?

* Better for what?

* Is it true?
<figure><img src="image_24.png" /></figure>

These are the reasons why this study is divided into three stages:

1. Comparing the measurement qualities by comparing their variograms

2. Establishing a formal link between the two measurements

3. Deducing linear systems enabling us to use the two types of measurements together
Data

<figure><img src="image_25.png" /></figure>

7

The data are from an open-pit copper mine in Northern Chile of where sub domain is analysed because it is almost homogeneously covered by around 3,000 drill-hole samples (3m long) and 13,000 blast-hole samples (15m long)
Drill holes, blast holes

Drill holes

<figure><img src="image_26.png" /></figure>

Blast holes

<figure><img src="image_27.png" /></figure>

Blasts and drills differs by their sampling density, their support size and their orientation

Horizontally we have approximately one blast every 10 meters and vertically every 15 meters

For the drills, it is a bit more complex because they are not all vertical and the grid is not regular
but we have approximately one drill hole every 50 meters horizontally, and every 3 meters along
the drill hole.

3m is assumed to be the support of the drill samples; 15 meters is assumed to be the support of the blast samples.

The global sampling ratio is approximately one drill sample to four blast samples
<figure><img src="image_28.png" /></figure>

The geostatistical comparison between the two types of measurements is divided into two steps:

Starting from the drill variogram, identifying the basic structures that model its behavior and deducing the underlying point-support variogram

Making the theoretical convolution of the point variogram on 15-meter long supports and comparing it to the blast variogram

It is important to distinguish two situations: variogram calculation parallel or perpendicular to the regularization direction because the formulae are not the same
<figure><img src="image_29.png" /></figure>

① Quality comparison (2/5)
Results

<figure><img src="image_30.png" /></figure>

These are the variograms.

On the left, the drill variogram, on the right, the blast one

We notice that the behaviors are similar and that they both contain a high percentage of nugget effect
① Quality comparison (3/5)
Charts

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /><figcaption>CHART NO. 11. Exponential model. Grading over l: γ<sub>o</sub>(L).</figcaption></figure>

We do not detail here the calculations here, they are in the paper, we just show the charts to be used.

They were done by hand forty-five years ago by a man called Jacky Laurent, who retired last year...
① Quality comparison (4/5)
Results

<figure><img src="image_33.png" /></figure>

3 structures were identified on the drill variograms, nugget effect, exponential and linear with a weak slope
Then, an underlying point-support model was deduced
This model was theoretically regularized over 15 meters

The bottom left-hand figure shows the vertical comparison. The dotted blue line represents the experimental vertical blast variogram, the dotted black line represents the present model and the red line the model we would obtain with a more realistic nugget effect.

One can see that apart from the problem of the nugget effect, the variation range is acceptable, even if the linear part of the theoretical structure does not appear in the vertical experimental blast variogram

The bottom right-hand figure shows the horizontal comparison. Again, apart from the nugget effect, the fitting is good.

The first conclusion is that if we omit the problem of the nugget effect, we see that both blast and drill holes can be considered a regularization of the same phenomena in accordance with their respective supports.

But the approach followed up to now suffers from two uncertainties:

* The analyses are done independently.

* The analyses refer to the drill-hole nugget which we assumed to be a "natural" micro structure; is this true?

To answer these questions, cross variograms must be calculated but we do not have any location with both measurements, so a migration is necessary.
① Quality comparison (5/5)

Cross analyses & conclusion

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

In order to obtain a significant number of measurements at the same location, around 1,000 blast samples were migrated to drill locations when the migration distance did not exceed 10 meters.

The right-hand figure presents the variograms.

Red points indicate the migrated blast variogram, black triangles show the drill variogram and the stars represent their cross variogram which does not show a significant nugget effect, possibly a small negative one without anything like the effects encountered on the individual variograms.

The conclusion is that the drill-holes have their own errors, independent of the blast ones, and the two measurements share only the structured parts of the variogram: the exponential and linear structures.
## ② Formal link

$$Y_{blast}(x,y,z) = Y(x,y,z) * p_{15m}(z) + R_{blast}(x,y,z)$$

<figure><img src="image_36.png" /></figure>

$$\gamma_{blast}(h) = \gamma_{15m}(h) + \gamma_{R_{blast}}(h)$$

<figure><img src="image_37.png" /></figure>

Blast sampling
error

<figure><img src="image_38.png" /></figure>

$$Y_{drill}(x,y,z) = Y(x,y,z) * p_{3m}(z) + R_{drill}(x,y,z)$$

<figure><img src="image_39.png" /></figure>

$$\gamma_{drill}(h) = \gamma_{3m}(h) + \gamma_{R_{drill}}(h)$$

<figure><img src="image_40.png" /></figure>

Drill sampling
error

<figure><img src="image_41.png" /></figure>

Finally we are able to establish a formal link between blast and drill measurements.
### ③ Consequence
Block modelling

$$ \begin{pmatrix} C * P_{3m} & C * P_{3m} * P_{15m} & 1 \\ C * P_{15m} * P_{3m} & C * P_{15m} + C_R \delta(0) & 1 \\ 1 & 1 & 0 \end{pmatrix} \begin{pmatrix} \lambda \\ \lambda' \\ \mu \end{pmatrix} = \begin{pmatrix} C * P_{3m} * P_v \\ C * P_{15m} * P_v \\ 1 \end{pmatrix} $$

<figure><img src="image_42.png" /></figure>

$Y_{blast}(u_j)$
$\lambda'_j$

Block kriging
using: • +

$Y_{drill}(u_i)$
$\lambda_i$

$$ Y_v(u) = \sum_i \lambda_i Y_{drill}(u_i) + \sum_j \lambda'_j Y_{blast}(u_j) $$

... and this link makes it possible to deduce linear systems like this one where the way to use blast and drill together in a single system is defined.

Thus, it becomes possible to improve the block model, based on drill holes, by integrating the blast holes as they are.

It also becomes possible to predict the production of the following days more accurately.
Filtering the blast error

$$ (C * P_{15m}) + \underbrace{C_R \delta(0)}_{1} \begin{pmatrix} 1 \\ 0 \end{pmatrix} \begin{pmatrix} \lambda \\ \mu \end{pmatrix} = (C * P_{15m}) $$

<figure><img src="image_43.png" /></figure>

$$ Y_{\text{blast}}^{*}(u) = \sum_{j} \lambda_{j} Y_{\text{blast}}(u_{j}) $$

One can also use this link to remove the blast error. The filter can be applied to each blast measurement, using a local neighborhood of surrounding blast samples. The system to be used is presented symbolically with a matrix formalism.

Such systems have been tested on a realistic simulation where everything is known (true point value, true block values, blast with or without errors).

The results will be shown at a congress next July in Northern Chile, but I can show you here one of the result concerning the filtering of the blast error by kriging.

For comparison, estimation is made by kriging with no filtering.
Filtering the blast error

<figure><img src="image_44.png" /></figure>

* Here are the results

* On these scatter diagrams, the horizontal axis represents the true blast without any sampling error

* On the left scatter diagram, the vertical axis is a usual kriging. The correlation with the truth is 0.65

* On the right, when the filtering is activated, the correlation increases to 0.9. Why?

* Because with filtering, the kriging neighborhood can incorporate the target point where the filter is applied. This point takes a high kriging weight (more than 65%). Although noisy, this point is closer to the truth than any average based on surrounding points which explains why the filter estimate is closer to the truth

* So finally, the advantage of this linear system is to enable the kriging neighborhood to incorporate the target point information
Conclusions

* Consistency

* Blast sampling errors, drill sampling errors

* (Co)kriging linear systems

In this deposit – and more generally, in this company, diamond drill hole grades and blast hole grades are consistent in the sense that, apart from the nugget effect, the structured part of their respective variograms follows the theoretical laws of regularization.

Concerning the nugget effects, we discover, by cross-analyses, that there is no natural micro-structure in the underlying point grade and the large nugget effects encountered on the variograms are due to blast and drill measurement errors.

In conclusion, some linear systems are proposed for removing the nugget effects from the data, and, more importantly, using blasts and drills together for short-term planning in mining.

These systems, among many other potential ones, easy to demonstrate, result directly from the formal link established here between blast and drill holes.

Before these systems are applied, the link must be verified according to the methodology presented here.
Acknowledgment

<figure><img src="image_45.png" /></figure>

CODELCO
Produciendo Futuro
Codelco Chili

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

MINES
ParisTech

• Codelco

• Chile

• France

• Paris School of Mines

Thank you for your attention

O yes I must apologize in advance: I have exactly 40 minutes left to reach the train station and I hope the taxi I ordered is waiting for me. Please, do not be surprised if I run out of the room when finished. I am not ill, I am not running away from the issues, I just have a train to catch....

Any question?
Chapter C

The oral presentation with comments of
“Diamond Drill Holes, Blast Holes
& Cokriging”

(June 2016 version)

Serge A. Séguret (Mines ParisTech, France)

Sebastian de La Fuente (Codelco, Chile)

A presentation made by S. DE LA FUENTE at Mine-Planning 2015, 4th
International Seminar on Mine Planning, July 8-10, 2015, Antofagasta, Chile.

(No opportunity for written papers at this congress)

Abstract

This talk is a geostatistical study of diamond drill holes and blast holes and their potential use
together in a cokriging system.

Both measurements of a real copper deposit are formally compared, leading to a model where
a blast hole can be considered as a regularization of the drill information up to a nugget effect
characteristic of the blasts.

This formal link makes it possible to build a cokriging system that accounts for the different
supports and leads to a block model based on blast and drill holes.

The model is tested on a realistic simulation where the true block grades, which are known,
are compared to their estimate obtained by:

- Kriging using only drill holes;

- Kriging using only blast holes;

- Cokriging using drill and blast holes together.

The first conclusion is that the best estimate is obtained when only blasts or blast and drill
holes are used together, there is no significant difference, there are too many blasts. This
justifies the usual practice consisting in basing the short-term planning on blasts only. But
another conclusion appears when kriging is compared to moving average (another common
practice), both based on blasts: depending on the number of data used in the neighborhood,
the moving average produces a strong conditional bias, a useful reminder of the reason why
Kriging was created more than fifty years ago.
DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING
(June 2016 version)

Serge Antoine Séguret, MINES ParisTech, Fontainebleau, France
Sebastian de La Fuente, Codelco, Radomiro Tomic, Chile

Presenting author: Serge Antoine Séguret

<figure><img src="image_48.png" /></figure>

CODELCO
Produciendo Futuro
Codelco Chili

Mine-Planning 2015

4th International Seminar on
Mine Planning

July 8-10, 2015

Antofagasta, Chile

<figure><img src="image_49.png" /></figure>

Good morning everyone

It is my pleasure to present the work entitled

"Diamond Drill-Holes, Blast-Holes and Cokriging"

... a work supported by Codelco and the Paris School of Mines
Drill holes, blast holes

Drill holes

<figure><img src="image_50.png" /></figure>

* "Few"

* Long term planning (month, year, decade)

* Large scale strategy

* "Good" quality

<figure><img src="image_51.png" /></figure>

Blast holes

<figure><img src="image_52.png" /></figure>

<figure><img src="image_53.png" /></figure>

Question: how
combining the two types
of measurements?

* Numerous

* Short term planning (day, week,)

* Small scale selectivity

* "Bad" quality

37

Typically in open pit mines, geologists, mining engineers, metallurgists, have at their disposal two types of measurements for the grades:

A first type, from drill holes – “diamond drill holes” in our case.

A second type, from the blast holes.

Because they are much more expensive, the diamond drill holes are less numerous than the blast holes, and it is usual to encounter sampling rates ranging from one over three to one over ten or worse.

Not only the sampling density is involved but also how samples are distributed in space too, as shown on the slide:

* The circles, representing the drill holes, are widely spaced, see the left horizontal cross-section

* In the same figure, the green crosses, representing the blast holes, are more densely spaced

* Vertically, on the right-hand figure, it is the reverse, with almost continuous drill hole information while the blasts are more widely spaced

These differences make it even more difficult to compare the statistical properties of the two types of measurements, including when calculating directional variograms because statistical inference conditions are not the same

Another difference concerns the way the measurements are used.

The long road that leads to the opening of the mine is marked by drilling campaigns, to achieve the block model that will condition the exploitation at large scale as well as for medium- and long-term planning. Typically, kriging and Geostatistics are used to build the model at this stage.

In addition, the blast holes are used for short term planning with no need of Geostatistics, a simple moving average is often used to estimate the block quantity of metal

These separate uses of two types of measurements that are supposed to represent the same thing raise questions about their relationship. In particular, would it not be possible to enrich the short-term estimate, now based only on blast holes, by adding the drill hole measurements?

Finally, we often hear, without real justification, that the diamond drill holes are much better than the blast ones.

We ask the questions:

* Better how?

* Better for what?

* Is it true?
Plan

Summary of a real study (porphyry copper)

<figure><img src="image_54.png" /></figure>

Formal link

<figure><img src="image_55.png" /></figure>

Linear systems:

<figure><img src="image_56.png" /></figure>

• Filtering the blast error

• Blast deconvolution

• Block modelling using blasts & drills

Building a realistic simulation

<figure><img src="image_57.png" /></figure>

Testing the systems

<figure><img src="image_58.png" /></figure>

To answer these questions, we propose the following steps:

➤ We begin by summarizing a formal study of blast and drill-holes which was presented some months ago at a sampling congress in Bordeaux, France

➤ This study makes possible some linear systems likes:

• Removing the blast error by kriging

• Estimating point-support values using blast measurements

• Bock modelling using blast and drill measurements together

We test the different systems on a simulation
<figure><img src="image_59.png" /></figure>

➤ The data are from an open-pit copper mine in Northern Chile where a sub domain was chosen for analysis because it is almost homogeneously covered by around 3,000 drill-hole samples (3m long) and 13,000 blast-hole samples (15m long)

➤ The geostatistical comparison between the two measurement types is divided into two steps:

• Starting from the drill variogram, identifying the basic structures that model its behavior and deducing the underlying point-support variogram

• Making the theoretical convolution of the point variogram on 15-meter long supports and comparing it to the blast variogram

➤ It is important to distinguish two situations: variogram calculation parallel or perpendicular to the regularization direction because the formulae are not the same

➤ Comparisons are completed by cross analyses based on migrated data.
Porphyry copper study
(summary 2/2)

<figure><img src="image_60.png" /></figure>

<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /></figure>

➤ 3 structures were identified on the drill variograms, nugget effect, exponential and linear with a weak slope

➤ The upper left-hand figure shows the vertical comparison. The dotted blue line represents the experimental vertical blast variogram, the dotted black line represents the present model and the red line the model we would obtain with a more realistic nugget effect. One can see that apart from the problem of the nugget effect, the variation range is acceptable, even if the linear part of the theoretical structure does not appear in the vertical experimental blast variogram

➤ The upper right-hand figure shows the horizontal comparison. Again, apart from the nugget effect, the fitting is good.

➤ If we omit the problem of the nugget effect, we see that both blast and drill holes can be considered a regularization of the same phenomenon in accordance with their respective supports.

➤ In order to obtain a significant number of measurements at the same location, around 1,000 blast samples were migrated to drill locations when the migration distance did not exceed 10 meters (see bottom left-hand figure).

➤ The bottom right-hand figure presents the variograms. Red points indicate the migrated blast variogram, black triangles show the drill variogram and the stars represent their cross variogram which does not show a significant nugget effect, possibly a small negative one without anything like the effects encountered on the individual variograms.
Model

$$
Y_{\text{blast}}(x,y,z) = Y(x,y,z) * p_{15m}(z) + R_{\text{blast}}(x,y,z)
$$

<figure><img src="image_63.png" /></figure>

$$
\gamma_{\text{blast}}(h) = \gamma_{15m}(h) + \gamma_{R_{\text{blast}}}(h)
$$

<figure><img src="image_64.png" /></figure>

Blast sampling
error

<figure><img src="image_65.png" /></figure>

$$
Y_{\text{drill}}(x,y,z) = Y(x,y,z) * p_{3m}(z) + R_{\text{drill}}(x,y,z)
$$

<figure><img src="image_66.png" /></figure>

$$
\gamma_{\text{drill}}(h) = \gamma_{3m}(h) + \gamma_{R_{\text{drill}}}(h)
$$

<figure><img src="image_67.png" /></figure>

Drill sampling
error

<figure><img src="image_68.png" /></figure>

➤ The conclusion is that the drill-holes have their own errors, independent of the blast ones, and the two measurements share only the structured parts of the variogram: the exponential and linear structures

➤ In this deposit – and more generally, in this company - diamond drill hole grades and blast hole grades are consistent in the sense that, apart from the nugget effect, the structured part of their respective variograms follows the theoretical laws of regularization. This result may surprise but one must admit that it is impossible to say that grade measurements from drill holes are better than those from blast one:

• Yes, there is a reduction of variance for the blast, but this is due to the
support which is larger

• In both cases, the ratio of nugget effect to the variance is around 40 percent

➤ This result justifies the present practices where selectivity is based only on blasts

➤ This formal link establishes that the blasts (respectively the drills) are considered as a regularization over 15 meters (respectively 3 meters) of the point-support copper grade “Y of x, y, z” plus a residual “R” proper to each type of measurement and independent of each other.
DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING

<figure><img src="image_69.png" /></figure>

* A refined simulation of point-support grades was done every meter horizontally and every 20 centimeters vertically

* Then, 100 vertical drill holes were created by averaging, every 50 meters horizontally, all the values along 3 meters vertically. In this way we obtain more than 4000 drill samples

* In the same way, more than 8000 15 meter-long vertical blast holes were produced with a horizontal spacing of 12 meters

* The true 15 meter cubic block value is found by averaging the more than 150 million point-support values contained in the block

* The sampling ratio is approximately one drill to three blasts and the domain covered by the simulation is 400 by 400 meters squared horizontally and 100 meters vertically
Geostatistical simulation
(2/2)

Points

<figure><img src="image_70.png" /></figure>

Drill
samples

<figure><img src="image_71.png" /></figure>

Blast
samples

<figure><img src="image_72.png" /></figure>

Blocks

<figure><img src="image_73.png" /></figure>

➤ The grades have an exponential variogram with a practical range of around 30 meters

➤ For simplification, the drills have no errors

➤ We add to the blasts a random noise with a nugget effect of 0.2, representing the blast sampling error

➤ The grades are realistic with 0 as minimum, 3.5 as maximum, an average equal to 0.63 as in the real deposit and the distribution has a correct right queued

➤ We verify that the sills of drill, blast and block variograms obey the laws of regularization
1 Filtering the blast error

$$ (C * P_{15m} + C_R \delta(0)) \begin{pmatrix} 1 \\ 0 \end{pmatrix} \begin{pmatrix} \lambda \\ \mu \end{pmatrix} = (C * P_{15m}) $$

<figure><img src="image_74.png" /></figure>

$$ Y_{\text{blast}}^{*}(u) = \sum_{j} \lambda_{j} Y_{\text{blast}}(u_{j}) $$

> As we are in a stationary case, matrices are presented using covariance C, it is simpler

> The first test that we propose is to remove the blast error. This filter can be applied to each blast measurement, using a local neighborhood of surrounding blast samples.

> The system is presented symbolically with a matrix formalism. In the system, "gamma R" disappears from the second member of the linear system. In this way we remove from the estimation, the part associated with the measurement error.

> Generally speaking, this does not mean that in the remaining part "gamma star P", there is no nugget effect, it means that only the "natural" part remains i.e. the part due only to the sampling error. In our case, we recall that there is no natural nugget effect

> The neighborhood must contain the sample from which the noise has been removed; otherwise, the filtering is not efficient.

> For comparison, the estimation is made by kriging with no filtering, using the same neighborhood (but without the target sample this time, otherwise kriging will give back the value of the data point).

> We select, among all the simulated blasts, a subset of around 1000 samples on which estimations will be conducted, using the additional samples (in case of ordinary kriging) and all the samples in case of nugget filtering.

> The reference is the "truth" i.e. the blast without errors which we know because we work on a simulation where everything is known
<figure><img src="image_75.png" /></figure>

* Here are the results

* On these scatter diagrams, the horizontal axis represents the true blast without any sampling error

* On the left scatter diagram, the vertical axis is a usual kriging. The correlation with the truth is 0.65

* On the right, when the filtering is activated, the correlation increases to 0.9. Why?

* Because with filtering, the kriging neighborhood can incorporate the target point where the filter is applied. This point takes a high kriging weight (more than 65%). Although noisy, this point is closer to the truth than any average based on surrounding points which explains why the filter estimate is closer to the truth

* So finally, the advantage of this linear system is to enable the kriging neighborhood to incorporate the target point information
2 Blast deconvolution

<figure><img src="image_76.png" /></figure>

<figure><img src="image_77.png" /></figure>

$Y_{blast}(u_j)$

$Y_{point}^*(u) = \sum_j \lambda_j Y_{blast}(u_j)$

➤ Now we propose a second test

➤ It can be interesting to remove the effect of regularization on the blast with a kriging system that estimates, for each blast measurement, a "point" value, while simultaneously removing the part of the nugget effect associated with blast errors.

➤ The difference with the previous system is that in the second member of the system, the capital letter for "P" - a two-time convolution - is replaced by a small "p" - a one-time convolution

➤ A comparison is made with the true point value, and with previous estimates (estimating a blast with or without nugget effect).
Blast deconvolution

<figure><img src="image_78.png" /></figure>

<figure><img src="image_79.png" /></figure>

<figure><img src="image_80.png" /></figure>

* In the three scatter diagrams, the horizontal axis represents the true point-support value

* The upper left scatter figure presents the result when a deconvolution is made, together with error filtering. The correlation with the true value is good, it equals 0.85

* The upper right-hand figure presents the result of the error removal with no deconvolution. It corresponds to the previously presented system but this time as compared with the point support value, the reason why the correlation equals 0.8 and not 0.9 when compared to the blast values. In comparison with the left-hand figure, the deconvolution increases significantly the accuracy of the estimation

* The bottom figure shows the results when no filtering and no deconvolution is done. The correlation is very low, it equals 0.55. The reason is the same as in the previous case. Filtering and/or deconvolution can use the target points where the filter is applied
### 3 Cokriging using blasts and drills

$$ \begin{pmatrix} C * P_{3m} \\ C * p_{15m} * p_{3m} \end{pmatrix} \begin{pmatrix} C * p_{3m} * p_{15m} & 1 \\ C * P_{15m} + C_R \delta(0) & 1 \end{pmatrix} \begin{pmatrix} \lambda \\ \lambda' \\ \mu \end{pmatrix} = \begin{pmatrix} C * p_{3m} * p_v \\ C * p_{15m} * p_v \end{pmatrix} $$

<figure><img src="image_81.png" /></figure>

* The third test is made to locally renew the mine planning block model by using blasts and drills together

* In the figure, red points represent the drill samples, green crosses the blast samples. The final estimation is a combination of the two measurements.

* This is a cokriging system with linked mean because drills and blasts have the same average, which is mandatory for carrying out all these calculations

* The objective is estimating the average grade at the block scale V

* We compare with the other systems:

* Block grade estimate by kriging using only drill holes

* Block grade estimate by kriging using only blast holes
Cokriging using blasts and drills

<figure><img src="image_82.png" /></figure>

<figure><img src="image_83.png" /></figure>

<figure><img src="image_84.png" /></figure>

➤ Here are the results. For the three scatter diagrams, the horizontal axis is the true block grade

➤ The upper left diagram is the result obtained by ordinary kriging using drills; upper right diagram the result when using blasts. The jump by the correlation coefficient from 0.4 (OK using drills) to 0.9 (OK using blasts) is impressive. Even if the blasts are regularized over 15m, the fact that they are more numerous and respect the variogram (up to a nugget effect) justifies their use when possible, instead of the drills, and the use of the blasts in selection for mining operations

➤ The bottom diagram concerns cokriging using blast and drill together. The performance is similar to Ordinary Kriging using only blasts. In our case cokriging is not useful because the blasts are so numerous and of such good quality that adding a drill contribution does not improve the results

➤ Notice that the nugget effect of the blasts has no impact on the results because with a block estimate, the average of the variograms is obtained by a random sampling of the block which neutralizes the nugget effect.
Moving average or kriging?

<figure><img src="image_85.png" /></figure>

<figure><img src="image_86.png" /></figure>

<figure><img src="image_87.png" /></figure>

Actually, the short term planning is based on average using the blast included in the block, so question is asked to see if kriging can give some improvement. We are comparing three experiments:

- Ordinary Kriging using 24 surrounding blasts measures (previous work);

- Moving average using the same 24 surrounding blast measures;

- Moving average using 4 blast measures of the same level.

Using the same data points, replacing Ordinary kriging by an average reduces the correlation with the truth from 0.884 to 0.741. This is a very important reduction which should indicate the practitioners to use kriging.

Now, if practitioners want to keep their habits, one can see that when using only four points, the result is better that when using 24 points because the smoothing is less important: the correlation with the truth goes from 0.741 to 0.839, a performance still under the performance when using kriging.
DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING

<figure><img src="image_88.png" /></figure>

But using so few points is risky for conditional bias reasons. To illustrate this concept, let us consider again the previous scatter diagram, but this time with the truth for the vertical axis and the curve of the mathematical expectation of the truth given different estimates. We focus on the most representative [0.3, 1] range of grades.

The red curves represent the discrete calculation of the conditional expectation of the truth, given the different estimates.

When doing kriging using 24 points, the conditional expectation curve is close to the first diagonal and when we select the block according to their estimates, we get in average what we expect, with perhaps a tiny tendency to underestimate the high grades.

When we replace kriging by a moving average using 24 points, we are still close to the diagonal, with a tiny tendency to overestimate the low grades and underestimate the high grades.

When we still do a moving average, but with only 4 points this time, the conditional bias appears clearly: in the range of the low grade, we systematically underestimate the average grade of the blocks and can decide to classify as waste blocks which are in practice richer than expected. Reversely, in the range of the high grades, this moving average with only 4 points systematically overestimates the average grade of the block so that we classify as rich blocks which must be considered as waste.

This is for all these reasons that one must use enough points in the kriging neighborhood (let us say at least 20), and reason why kriging and Geostatistics have been created, more than 50 years ago. It is perhaps useful to recall it here.
Conclusion

* Blasts, not so bad

* Drills, not so good

* Coherency between both

* Many linear methods (Error filtering, deconvolution...)

* Blasts sampling density → Cokriging not useful

The study of a porphyry copper deposit showed a formal link between blast and drill holes, leading to numerous linear systems able at least to:

- Filter the blast error

- Blast deconvolution

- Block modelling using blasts & drills

Tested on a simulation, these systems has proven their interest, as well as the danger to replace kriging by a moving average, especially when using few points.
Acknowledgment

<figure><img src="image_89.png" /></figure>

CODELCO
Produciendo Futura
Codelco Chili

<figure><img src="image_90.png" /></figure>

<figure><img src="image_91.png" /></figure>

MINES
ParisTech

• Codelco

• Chile

• France

• Paris School of Mines
<figure><img src="image_92.png" /></figure>
Chapter D

The oral presentation with comments of

“Anisotropy of the Rock Quality Designation (RQD) & its Geostatistical Evaluation

Serge A. Séguret (Mines ParisTech, France)

Cristian Guajardo (Codelco, Chile)

Claudio Rojas (Codelco, Chile)

A presentation made by C. Guajardo at Geomin 2015, 4th International Seminar on Geology for the Mining Industry, July 8-10, 2015, Antofagasta, Chile.

(No opportunity for written papers at this congress)

Abstract

Rock Quality Designation (RQD) measures the borehole core recovery percentage incorporating only pieces of solid core that are longer than 100 mm measured along the centerline of the core. It is an important attribute used in geotechnics via Rock Mass Rating (RMR) for example.

This presentation concerns the behavior of this attribute in a Chilean porphyry copper deposit by analyzing more than 60,000 samples 1.5-meter long.

The drill holes have different directions and the nature of RQD requires accounting for the sample direction if the fracture network is anisotropic, a concept different from the geostatistical anisotropy which measures the variability along a direction set by two samples. A directional analysis is conducted and shows different variograms associated with different sample direction classes; this leads to different maps obtained by kriging, calling into question the usual practices that do not account for the sample direction. In fact, the directionality of the measurements makes RQD non additive (and not krigeable) if the samples used in the kriging neighborhood belong to different directional classes.

A discussion follows; pointing out that RQD is subject to the same directional bias as the Fracture Frequency (FF) and should be corrected in the same way by a sinus of the angle between the sample and the fracture (Terzaghi correction). But such a correction is difficult, if not impossible; tests are presented.
ANISOTROPY OF THE ROCK QUALITY DESIGNATION
(RQD) & ITS GEOSTATISTICAL EVALUATION

Serge Antoine Séguret, MINES ParisTech, Fontainebleau, France
Cristian Guajardo, Codelco, Chuquicamata, Chile
Claudio Rojas, Codelco, Radomiro Tomic, Chile

<figure><img src="image_93.png" /></figure>

CODELCO
Produciendo Futuro
Codelco Chili

Geomin 2015

4th International Seminar on
Geology for the Mining Industry
July 8-10, 2015
Antofagasta, Chile

<figure><img src="image_94.png" /></figure>

GEOMIN 2015

International Conference on
Geology for the Mining Industry

<figure><img src="image_95.png" /></figure>

GECAMIN

It is my pleasure to present the work done jointly by my collaborator Doctor Serge Séguret from France, and myself, on data provided by my team and by my colleague Claudio Rojas

The presentation is entitled

"ANISOTROPY OF THE ROCK QUALITY DESIGNATION & ITS GEOSTATISTICAL EVALUATION".

... a work supported by the Chilean company Codelco and the Paris School of Mines where Serge has worked for over thirty years in the Geostatistical laboratory founded by Georges Matheron at Fontainebleau.
Plan

1. RQD

2. Additivity

3. Directionality

4. Case study

5. Conclusion

GEOMIN 2015

International Conference on
Geostatistical Designation of Rock Quality Designations

<figure><img src="image_96.png" /></figure>

GECAMIN
GEOMIN CONFERENCE FOR ROCK QUALITY DESIGNATION

Here is the program that we propose to follow.
<figure><img src="image_97.png" /></figure>

The aim of RQD is to measure the degree of jointing or fracturing in a rock mass. It is one of the main attributes incorporated in RMR which is a comprehensive index of rock-mass quality used for the design and construction of excavations in rock, such as tunnels, mines, slopes and foundations.
2. Additivity ?
(for a given sample direction)

<figure><img src="image_98.png" /></figure>

$$RQD_{L_1 \cup L_2} = 100 \frac{L_{>10\text{cm}, 1} + L_{>10\text{cm}, 2}}{L_1 + L_2} \iff RQD_{L_1} & RQD_{L_2} ?$$

GEOMIN 2015

International Geostatistical Conference
Technology for the Mining Industry

<figure><img src="image_99.png" /></figure>

GECAMIN
GEOSTATISTICAL EVALUATION

<figure><img src="image_100.png" /></figure>

> The first question concerns the additivity of RQD i.e. the ability to estimate it at any location by a linear combination of measurements like kriging in geostatistics

> Suppose that all the samples have been drilled along the same direction in space and take two values collected along two different supports L₁ and L₂, each with its own sum of core pieces longer than 10 cm, respectively equal to L1>10 and L2>10...
Additivity
(for a given sample direction)

$$RQD_{L_1 \cup L_2} = 100 \frac{L_{>10\text{cm}, 1} + L_{>10\text{cm}, 2}}{L_1 + L_2}$$

$$= \frac{L_1 100 \frac{L_{>10\text{cm}, 1}}{L_1} + L_2 100 \frac{L_{>10\text{cm}, 2}}{L_2}}{L_1 + L_2}$$

$$= \frac{L_1 RQD_{L_1} + L_2 RQD_{L_2}}{L_1 + L_2}$$

→ Kriging possible
(same support & sample direction)

GEOMIN 2015

International Conference on
Geostatistics for the Mining Industry

<figure><img src="image_101.png" /></figure>

GECAMIN
GEOMIN CONFERENCE FOR THE INDUSTRIAL MINING INDUSTRY

➤ The value of RQD over the support L1 and L2 is...

➤ And when we develop, we obtain...

➤ This is by definition the way to combine additive quantities, equal to the average when L1=L2

➤ So kriging RQD is authorized, at punctual or block scale, using classical geostatistical tools when RQD is order-two stationary

➤ Although the calculation is easy, the interpretation is not: as RQD is a 1D measurement, estimating it at block scale using 1D samples just gives the average behavior of 1D samples over a block, and this is not a 3D property
### 3. Anisotropy
("Directionality")

Geotechnical anisotropy ≠ Geostatistical anisotropy

<figure><img src="image_102.png" /></figure>

$M(x)_{\overline{Dir1}}$

$\neq M(x)_{\overline{Dir2}}$

$\neq M(x)_{\overline{Dir2}}$

<figure><img src="image_103.png" /></figure>

$\neq$ $VAR[M(x) - M(x + \bar{h})]_{\overline{Dir2}}$

<figure><img src="image_104.png" /></figure>

* Now, we have to face this problem: the samples have different drilling directions

* This is a very general problem that occurs when geostatistical tools are applied to tensors: the quantity that we measure may depend on the direction of the measurement

* It is important, at this stage, to distinguish between the interpretation of the word "anisotropy" when used by a geostatistician, and the same word, when used by a geotechnician or a hydrogeologist.

* In geostatistics, this word means that when we calculate a variogram, based on differences of the variable of interest ("M" on the slide), the variance of such increments depends on the direction and this leads to concepts like "zonal" or "geometric" anisotropies

* In Geotechnics and in Hydrogeology, when one says that the phenomenon is anisotropic, it means that the measure itself depends on the direction, not only the increment

* The two anisotropies are obviously linked, but the underlying concepts differ fundamentally and this is the reason why we employ the word "directionality" as a synonym for anisotropy as used by the geotechnicians and hydrogeologists
<figure><img src="image_105.png" /></figure>

If, for a given estimation, based on a given amount of samples – typically the measurement used in the kriging neighborhood – we mix samples with different directions, we obtain a result that is not useful, and even incorrect, because we mix quantities that are not comparable, due to their directionality

If the different directions are represented in the same way, one could say that the obtained result represents the average behavior of the quantity along all the possible directions...

... but in practice, locally, the sampling of the directions is never homogeneous and the result just reflects the anisotropy of the sampling

This is the reason why we recommend classing the samples according to their direction, and calculating direction by direction
RQD Directional?

Yes, 3 reasons

Reason 1

<figure><img src="image_106.png" /></figure>

Rare (but existing) event

* Now let us analyze the directionality of RQD

* The primary reason is that a core with a joint running along it, is considered solid (see figure, in red circles). If the sample direction is slightly changed, the joint is no longer aligned with the sample and a core, previously considered to be longer than 10 cm and incorporated into RQD, can be rejected, which changes the final value of the attribute. But this is a rare event that can be neglected in most cases.
RQD Directional?
Reason 2

<figure><img src="image_107.png" /></figure>

GEOMIN 2015

International Conference on
Geometrical Design and
Technology for Mining Systems

<figure><img src="image_108.png" /></figure>

GECAMIN
International Conference on
Geometrical Design and
Technology for Mining Systems

➤ The second reason is that the spacing between the fractures is subject to directional bias.

➤ Suppose that all the fractures are parallel planes aligned along a single direction (left figure) and consider two consecutive fractures (middle figure)

➤ The spacing measured along direction n°1 is $l_1$ while it should be “11Sinus thetha” if we want to make the measurement independent of the sample direction by referring to direction n°2 perpendicular to the fracture direction.

➤ This correction must be applied to each segment (right-hand figure) but it cannot be applied directly to its sum with the common factor “sinus thetha” because after correction, some segments may be shorter than 10 cm (case of segment $l_2$ in the right-hand figure) and consequently, they will not be included in the summation

➤ The problem becomes more complex when we consider that a direction in three dimensions requires two angles (azimuth and dip).
RQD Directional?
Reason 2

Same problem for Fracture Frequency

<figure><img src="image_109.png" /></figure>

$$L' = \frac{L}{\sin \theta}$$

$$FF = \frac{N_{tot}(x)}{L} = \frac{N_{tot\ corrected}(x)}{L'}$$

$$\rightarrow N_{tot\ corrected}(x) = \frac{N_{tot}(x)}{\sin(\theta)}$$

(Terzaghi, 1965)

GEOMIN 2015

8th International Geostatistical Conference
on Geostatistical Methodology for the Mining Industry

<figure><img src="image_110.png" /></figure>

GECAMIN
GEOSTATISTICAL EVALUATION

➤ Like RQD, FF is subject to a bias which depends on the angle θ between the fracture and the sample direction

➤ While this angle tends to zero (sample axis parallel to the fractures), the counting tends to under-evaluate the number of fractures

➤ Terzaghi (1965) proposed to correct this bias by multiplying the number of fractures by « 1/sinus thetha ». As opposed to RQD, this correction is possible because the fractures are classed.
<figure><img src="image_111.png" /></figure>

► If we assume that all the samples have the same direction, and the fracture network is isotropic (left-hand figure), the directionality of RQD can be neglected, whereas if there are local anisotropies in the fracture network, the local value of RQD at location “x” will depend on their angles and will not be comparable to another value of RQD at an another location “x”' even close to “x” (right-hand figure)

► In the first case RQD measures the degree of fracturation in an objective way

► In the second case it also measures the anisotropy of the fracture network and this is not really its objective
<figure><img src="image_112.png" /></figure>

For all these reasons we recommend the following procedure:

the samples must first be classed according to their directions, and structural analyses must be conducted by direction

If the variograms by sample directions are close to each other, and close to the variogram obtained with all the samples, there is no reason to distinguish between the directions because the fracture network looks isotropic

If this is not the case, kriging must be conducted only with the samples of a given directional sampling class, yielding as many directional results as there are classes

If, for each location x, where the estimation is conducted, the results are always the same, it means that the directionality of the measurement can be neglected and kriging can be conducted with all the samples without distinction between the classes

If the differences are large, the result is directional, which is important for the geotechnicians because it means that, when looking at the rock-strength prior to drilling a tunnel, for example, they must account for the drilling direction

These ideas are illustrated below by a case study
<figure><img src="image_113.png" /></figure>

* We are at Antofagasta and everybody knows the map of Chile but this presentation will also be done in France and Germany, the reason why this map is here

* The data are from a copper mine in northern Chile. We have used more than 60,000 samples with RQD information

* All the samples are 1.5 m long

* They cover a 7 km by 3 km by 600 m domain
Case study
Sample direction classification

<figure><img src="image_114.png" /></figure>

<figure><img src="image_115.png" /></figure>

<figure><img src="image_116.png" /></figure>

<figure><img src="image_117.png" /></figure>

<figure><img src="image_118.png" /></figure>

<figure><img src="image_119.png" /></figure>

* All the drill holes belong to vertical West-East oriented planes and the sample directions can be summarized by six angular classes
<figure><img src="image_120.png" /></figure>

* For each angular class a RQD variogram is calculated

* As no geostatistical anisotropy is noticed, we use omnidirectional variograms

* One notices two variogram groups: sampling directions 1 and 6 on the one hand, and sampling directions 2, 3, 4, 5 on the other

* The variance increases with the dip

* The first five directions share approximately the same nugget effect
Case study
Models by sample direction

<figure><img src="image_121.png" /></figure>

<figure><img src="image_122.png" /></figure>

<figure><img src="image_123.png" /></figure>

<figure><img src="image_124.png" /></figure>

<figure><img src="image_125.png" /></figure>

<figure><img src="image_126.png" /></figure>

<figure><img src="image_127.png" /></figure>

➤ We fit a model for each one of the six sampling directions, and a seventh one taking into account all the samples

➤ Estimations are conducted by sampling direction, each time using only the samples associated with the direction
<figure><img src="image_128.png" /></figure>

* Here we present a typical horizontal cross-section

* As for the variograms, the directional RQD maps can be regrouped into two sets: {1, 6} (i.e. close to a horizontal sampling direction) and {2, 3, 4, 5} (i.e. almost vertical)

* The greatest differences concern the northern part which contains low RQD for the second set of directions

* In this region, the contrasts are extreme: for an almost horizontal direction (and along West-East), RQD is equal to 70% on average while it becomes lower than 20% when sampling is perpendicular

* This domain is probably composed of a stack of horizontal planes separated by fractures
Case study
No direction selection

<figure><img src="image_129.png" /></figure>

• Ordinary Kriging

• Moving neighborhood
≈ 30 points ANY Dir.

<figure><img src="image_130.png" /></figure>

<figure><img src="image_131.png" /></figure>

➤ The left-hand map is obtained when no distinction of the sample direction is considered

➤ Such a map is of no interest because it just shows the dominant directions of the sampling

➤ The right-hand maps differ the most
<figure><img src="image_132.png" /></figure>

* What is the truth?

* There is no truth

* RQD is directional and subjected to a bias linked to the angle between the sample and the fracture

* Previous results showed that direction 2 differs strongly from direction 6. What does that mean?

* We recall that the samples all belong to vertical East-West oriented planes. When we calculate RQD, we measure the intersection between these measurement planes and the planes defined by the fracture. For the horizontal cross-sections of the figure, RQD along direction 2 (105°) is close to 20% in the upper western part while along direction 6, RQD becomes close to 70%. One can deduce that in this part of the deposit, the intersection of the measurement planes with the fractures tends to be parallel to 50°, with a high density

* Such analyses, made plane by plane and by direction may help to detect particular domains that go unnoticed if the sample direction is not taken into account.

* Finally, while the RQD aim is to measure the degree of jointing or fracturing in a rock mass, it is clear that it measures mainly the anisotropy of the fracture network and it seems very important to account for the sample direction to ensure that the results are useful

* The conclusion is that the directionality of RQD, as well as of FF, must be accepted and used and it means that attributes like RMR become directional too, which may disturb the practitioner's habits. But does it matter so much? After all, when drilling a tunnel, it seems natural to consider that the rock strength depends on the drilling direction.
Acknowledgment

<figure><img src="image_133.png" /></figure>

CODELCO
Produciendo Futura
Codelco Chili

<figure><img src="image_134.png" /></figure>

<figure><img src="image_135.png" /></figure>

MINES
ParisTech

• Codelco

• Chile

• France

• Paris School of Mines
<figure><img src="image_136.png" /></figure>
Chapter E

Geostatistical Evaluation of Rock-Quality
Designation and its link with Linear Fracture
Frequency

Serge A. Séguret (Mines ParisTech, France)

Cristian Guajardo (Codelco, Chile)

A paper presented at IAMG 2015, 17th annual conference of the International Association for Mathematical Geosciences, September 5-13, Freiberg, Germany

Abstract

Rock Quality Designation (RQD) is an important attribute used in geotechnics for quantifying the rock quality. It measures the borehole core recovery percentage incorporating only pieces of solid core that are longer than 100 mm measured along the centerline of the core.

The presentation examines the behavior of this attribute in a Chilean porphyry copper deposit
by analyzing more than 60,000 1.5-meter long samples.

The drill holes have different directions and the nature of RQD requires accounting for the
sample direction if the fracture network is anisotropic, a concept different from the
geostatistical anisotropy which measures the variability along a direction set by two samples.

A directional analysis is conducted and shows different variograms associated with different sample direction classes. This leads to different maps, calling into question the usual practices which do not account for the sample direction.

The second part of the presentation concerns the link with the linear Fracture Frequency (FF), another important attribute which measures the number of discontinuities per meter. Under the assumption that the discontinuities along a line follow a Poisson process, Priest & Hudson established in 1976 a formula which expresses RQD as a function of FF. This formula is compared to E[RQF|FF], the mathematical expectation of RQD given FF, deduced from the data. The result of the comparison depends on whether FF is or is not corrected by the sinus of the angle between the sample direction and the fracture, as recommended by Terzaghi in 1965. When applied to FF, this correction breaks a natural correlation with RQD which appears when no correlation is applied. In the latter case, the Priest & Hudson formula is acceptable, in the first case it is not. So there is a dilemma: on the one hand, Terzaghi looks necessary to correctly calculate FF, on the other, the correction systematically increases FF and reduces the relative influence of RQD when both attributes are incorporated into an overall rating like the Rock Mass Rating (RMR), for example.

A discussion follows, pointing out that RQD is subject to the same directional bias as FF and
should be corrected in the same way by a sinus of the angle between the sample and the
fracture, but such a correction is difficult, if not impossible; tests are presented. Finally, a
correction of RQD is proposed, based on the Priest & Hudson formula.
<figure><img src="image_137.png" /></figure>

# Geostatistical Evaluation of Rock-Quality Designation and its link with Linear Fracture Frequency

S. A. SEGURET<fnref n="1" />* and C. GUAJARDO<fnref n="2" />

## 1 Introduction

For mining operations in the developed world, staff safety has become a major concern where knowledge, at different scales, of the rock-mass strength is required in order to adequately size the retaining walls and other underground infrastructures, or simply the benches and face slopes in open-pit mining.

One of the identified causes of fragility is fracturation. By fracturation we mean discontinuities visible on core samples, including “open” fractures or just “joints”.

Fracturation can be tested in the laboratory on samples, providing attributes such as the Point Load Test (PLT, Bieniawski 1975); or measured on core samples, such as the two attributes of interest in this paper: the Rock-Quality Designation (Deere & Al., 1967; Stagg & Zienkiesicz, 1968) and the (Linear) Fracture Frequency (Jaeger & Cook, 1969).

RQD relates to measurements of intact lengths, FF to the counting of discontinuities, and these two quantities are related: when FF increases, RQD tends to decrease, but in a non-deterministic manner as shown by many studies (Priest & Hudson 1981, Sen & Kazi 1984).

The objective of this paper is to study this relation, together with one of the features common to both measurements: their directionality. By directionality, we simply mean that the measurement depends on the direction of the survey i.e. that these quantities are the components of a tensor, a feature rarely, if ever, considered in practice.

If the directionality of FF can be forgotten by correcting this measurement for an angular bias (Terzaghi, 1965), such an operation is not possible with RQD, and we show that correcting FF but not RQD involves a distortion of their correlations, leading to increased influence of FF to the detriment of RQD in the Rock-Mass Rating, which is a comprehensive index of rock-mass quality used for the design and construction of rock excavations (Barton & Al., 1974). This creates a dilemma: either FF is corrected, inducing the described problems, or FF is not corrected and one must account for the sample direction during the estimation of FF at a grid nodes, the ultimate goal of this work. A third way is proposed: referring to a theoretical model linking FF and RQD, and using this reference to correct RQD. The chosen model is due to Priest & Hudson (1976). Although it is based on a questionable assumption, it has been widely studied in the literature and proved to correctly reflect the link between the two attributes (Wallis & King 1980; Goodman & Smith 1980).

In the following, we do not address the thorny problem of making measurements made along lines representative of a volume.

## 2 Rock-Quality Designation (RQD)

### 2.1 Definition

The aim of RQD is to measure the degree of jointing or fracturing in a rock mass. It measures the borehole core recovery percentage incorporating only pieces of solid core that are longer than 100 mm measured along the centerline of the core. Figure 1 shows an example. Notice that joints along the core axis are treated as solid core. For mapping RQD, one has to consider samples of the same length (the case in the following application).

<footnote n="1">¹Geosciences/Geostatistics – Mines ParisTech, France, [serge.seguret@mines-ParisTech.fr](mailto:serge.seguret@mines-ParisTech.fr)</footnote>

<footnote n="2">²Geotechnique - Codelco, Chile</footnote>

<footnote>*presenting author</footnote>
<figure><img src="image_138.png" /><figcaption>Figure 1. Detail of RQD calculation. Property highlighted in red is a cause of directionality</figcaption></figure>

## 2.2 Directionality

RQD depends on the angle between the sample direction and the fracture; the primary reason is that a core with a joint running along it is considered solid (Figure 1, in red circles). If the sample direction is slightly changed, the joint is no longer aligned with the sample and a core, previously considered as longer than 10 cm and incorporated into RQD, can be rejected, which changes the final value of the attribute.

The second reason is that the spacings between the fractures are subject to directional bias. Suppose that all the fractures are parallel planes aligned along a single direction (Figure 2a) and consider two consecutive fractures (Figure 2b).

<figure><img src="image_139.png" /><figcaption>Figure 2. RQD directionality, detail of the calculations</figcaption></figure>

The spacing measured along direction n°1 is $l_1$ while it should be $l_1 \sin \theta$ if we want to make the measurement independent of the sample direction by referring to direction n°2 perpendicular to the fracture direction (similar to the Terzaghi correction presented later). This correction must be applied to each segment (Figure 2c) but it cannot be applied directly to their sum with the common factor $\sin(\theta)$ because after correction some segments may be shorter than 10 cm (case of segment $l_2$ in Figure 2c) and consequently, they will not be included in the summation. The problem becomes more complex when we consider that a direction in three dimensions requires two angles (azimuth and dip) while, in practice, only the angle between the fracture and the sample is roughly evaluated on the samples.

If we assume that all the samples have the same direction, and the fracture network is isotropic (Figure 3a), the directionality of RQD can be neglected, whereas if there are local anisotropies in the fracture network, the local value of RQD at location “x” will depend on their angles and will not be comparable to another value of RQD at an another location “x” even close to “x” (Figure 3b). In the first case RQD measures the degree of fracturation in an objective way but, in the second case it also measures the anisotropy of the fracture network and this is not really its objective.
<figure><img src="image_140.png" /><figcaption>Figure 3. a) Not considering the RQD directionality implicitly supposes that the fracture network is isotropic. b) In case of anisotropy, RQD depends on the sample direction.</figcaption></figure>

## 2.3 Additivity & Estimation procedure

It results from the previous discussion that to map RQD, one must account for the sample direction. Now if we assume all the sample directions to be identical, this attribute is additive. To prove this, take two values known along two different supports $L_1$ and $L_2$, each with their own sum of core pieces longer than 10 cm, respectively equal to $L_{>10cm, 1}$ and $L_{>10cm, 2}$.

The value of RQD associated with the support $L_1 \cup L_2$ is:

$$RQD_{L_1 \cup L_2} = 100 \frac{L_{>10cm, 1} + L_{>10cm, 2}}{L_1 + L_2} = \frac{L_1 100 \frac{L_{>10cm, 1}}{L_1} + L_2 100 \frac{L_{>10cm, 2}}{L_2}}{L_1 + L_2} = \frac{L_1 RQD_{L_1} + L_2 RQD_{L_2}}{L_1 + L_2}$$

This is, by definition, the way to combine additive quantities, equal to the average when $L_1=L_2$. So kriging RQD is authorized, at punctual or block scale, using classical geostatistical tools when RQD is order-two stationary.

If the calculation is easy, the interpretation is not: as RQD is a 1D measurement, estimating it at block scale using 1D samples just gives the average behavior of 1D samples over a block, and this is not a 3D property.

The samples have different directions and to make possible an interpretation of the estimation result, the samples must first be classed according to their directions, and structural analyses must be conducted by direction. If the variograms by sample directions are close to each other, and close to the variogram obtained with all the samples, there is no reason to distinguish between the directions because the fracture network looks isotropic. If this is not the case, kriging must be conducted only with the samples of a given directional sampling class, giving as many directional results as there are classes. If, for each location x, where the estimation is conducted, the results are always similar, it means that the directionality of the measurement can be neglected and kriging can be conducted with all the samples without distinction between the classes. If the differences are large, the result is directional, which is important for the geotechnicians because it means that, when looking at the rock-strength to drill a tunnel, for example, they must account for the drilling direction. This situation is discussed later in the application section.

# 3 Linear Fracture Frequency (FF)

## 3.1 Definition

The (Linear) Fracture Frequency is basically the ratio of a number of fractures, counted by the geologist, divided by the sample length. But the calculation is not simple because a significant part of the sample may be crushed, making the fracture counting possible only on the non crushed part. FF then becomes the ratio of two quantities both of which change from one location to another in the ore
deposit. It is not additive and needs some care when estimated by kriging (Seguret & al., 2014). In the following and for the application, the samples used have no crushed parts.

For each fracture, the geologist evaluates its angle with the sample axis without considering the fracture plane, and the fracture is classed into one of the three angular classes: [0°, 30°], [30°, 60°], [60°, 90°]. Normally, a direction in 3D is defined by the vector of the azimuth and the dip, but the only measurement at our disposal is θ, the angle between the sample axis and the fracture. In Fig. 4a, both fractures belong to the same [0, 30°] set because no orientation is given to the sample axis and the smaller angle is taken.

<figure><img src="image_141.png" /><figcaption>Figure 4. The angular classification is crude. a) Both fractures belong to the same set. c) In the middle class (grey), two fractures can be perpendicular</figcaption></figure>

The possible directions included in each class are cones represented by a cross-section in Fig. 4b. The classes are not equivalent: [0°, 30°] and [60°, 90°] classes may contain fractures that have an angle of up to 60° maximum while class [30°, 60°] may contain perpendicular fractures. This affects the quality of the Terzaghi correction presented below and makes it controversial.

## 3.2 Terzaghi correction

Like RQD, FF is subject to a bias which depends on the angle θ between the fracture and the sample direction (Figure 2a). While this angle tends to zero (sample axis parallel to the fractures), the counting tends to under-evaluate the number of fractures. Terzaghi (1965) proposed to correct this bias by multiplying the number of fractures by $\frac{1}{\sin\theta}$. As opposed to RQD, this correction is possible

because the fractures are classed. In this paper, the coefficients 3.86, 1.414 and 1.035 are the multiplicative factors used for respectively [0°30°], [30°,60°], [60°,90°] classes.

## 3.3 RQD & FF

Sometimes only FF is measured while RQD may also be required, which is the reason why trials for a formal link between the two variables have been established since the seventies, a famous one being due to Priest & Hudson in 1976.

Under the assumption that the distribution of spacings between discontinuities along a line follows a negative exponential one, the authors deduce that:

$$RQD = 100e^{-0.1FF}(0.1FF + 1) \qquad (1)$$

In the theory developed by these authors, a distribution of spacing is associated with each value of FF, as follows:

$$f(x) = FFe^{-FFx} \quad \text{with } f(x), \text{ frequency of the spacing } x \qquad (2)$$

The quality of this assumption has been discussed in detail by de Marsily & Chilès in Bear & al. (1993) and it is not the aim of the present paper to continue this discussion. We use (1) to better illustrate the distortions of the correlation between FF and RQD, providing a kind of visual reference.
# 4 Application

## 4.1 Data

The data come from a northern Chilean copper mine. More than 50,000 samples are used, with FF (whether Terzaghi-corrected or not) and RQD information. All the samples are 1.5 m long. They cover a 7 km by 3 km by 600 m domain (Figure 5).

<figure><img src="image_142.png" /><figcaption>Figure 5. Base maps of the samples.</figcaption></figure>

All the drill holes belong to vertical West-East oriented planes and the sample directions can be summarized by six angular classes (Table 1):

<table><caption>Table 1. Sampling direction classes</caption><thead><tr><td>Angular class</td><td>Azimuth</td><td>Dip</td></tr></thead><tbody><tr><td>1</td><td>[245°, 295°]</td><td>[-57°, -42°]</td></tr><tr><td>2</td><td>[245°, 295°]</td><td>[-69°, -58°]</td></tr><tr><td>3</td><td>[245°, 295°]</td><td>[-80°, -70°]</td></tr><tr><td>4</td><td>[66°, 115°]</td><td>[-80°, -70°]</td></tr><tr><td>5</td><td>[66°, 115°]</td><td>[-67°, -57°]</td></tr><tr><td>6</td><td>[66°, 115°]</td><td>[-56°, -42°]</td></tr></tbody></table>

## 4.2 RQD

For each angular class a RQD variogram is calculated. As no geostatistical anisotropy is noticed, Figure 6 presents omnidirectional variograms.

<figure><img src="image_143.png" /><figcaption>Figure 6. RQD variograms for different sample directions 1 to 6. Directions are defined by Table 1. All drill holes belong to West-East vertical plans.</figcaption></figure>

One notices two variogram groups: sampling directions 1 and 6 on the one hand and sampling directions 2, 3, 4, 5 on the other. The variance increases with the dip. Notice that the first five directions share approximately the same nugget effect.
We fit a model for each one of the six sampling directions, and a seventh one taking all the samples. Estimations are conducted by sampling direction, each time using only the samples associated with the direction. Figures 7 present a horizontal cross-section.

<figure><img src="image_144.png" /><figcaption>Figure 7. Horizontal cross-section of RQD estimation. From a) to f), RQD maps by sampling directions. g) RQD map using all the samples together</figcaption></figure>

### 4.2.4 Discussion

As for the variograms, the directional RQD maps can be regrouped into two sets: {1, 6} (i.e. close to a horizontal sampling direction) and {2, 3, 4, 5} (i.e. almost vertical). The greatest differences concern the northern part which contains low RQD for the second set of directions. In this region, the contrasts are extreme: for an almost horizontal direction (and along West-East), RQD is equal to 70% on average while it becomes lower than 20% when sampling is perpendicular. This domain is probably composed of a stack of horizontal planes separated by fractures and it is meaningless to consider averaged RQD (i.e. averaged along the different directions, Figure 7g) because such a map just shows the dominant directions of the sampling. If RQD is not corrected in a way comparable to that of FF, it must be calculated direction by direction. This is the first important conclusion of this work.

## 4.3 RQD & FF

What about RQD compared to FF? Figure 8a shows a scatter diagram between FF (when Terzaghi-corrected) and RQD. The black curve represents the conditional expectation which is, for a given value of FF, the average of the different values of RQD. The behavior predicted by Priest & Hudson's formula is plotted in red. The least one can say is that this formula does not reflect the experimental average behavior. Does this mean that Priest & Hudson's hypotheses are not acceptable? Not necessarily because when RQD is compared to FF without the Terzaghi-correction (Figure 8b),
one can see that the formula slightly over-estimates RQD, on average in the range [0, 28], but the
differences are not comparable to the previous differences when Terzaghi was applied.

<figure><img src="image_145.png" /><figcaption>Figure 8. RQD versus FF. Red curve: Priest&Hudson model; black curve: conditional expection of RQD given FF. a) FF is Terzaghi-corrected. b) FF is not Terzaghi-corrected</figcaption></figure>

What happens? Terzaghi correction is necessary to make the discontinuity counting comparable from
one direction to another and finally transform the directional number into a counting systematically
perpendicular to the fractures. If we do not correct the measurements in this way, we add numbers
which are not comparable- the directionality of the measurement makes it non additive when the
sampling directions are mixed. For RQD, finally, the problem is the same. The distance between two
discontinuities should be corrected according to the angle of the sample with the fractures, so that
when we summarize the lengths over 10cm, we summarize quantities that are comparable. This
correction was not done and could not be done because the two fractures that delimitate the core
length do not necessarily have the same direction; they cannot be classed as the fractures are.
Consequently when we compare an uncorrected RQD with a Terzaghi-corrected FF, we compare two
populations that are not comparable because of the correction. When we compare RQD to FF without
correction, we compare coherent populations because the bias in the counting is fully compatible with
the bias in the associated length and then Priest & Hudson's formula is verified in an acceptable way.

**4.4 Correction trials**

From the total amount of data, a sub-set of around 4,500 samples was selected because they have at
least 3 fractures with the same direction (or more precisely, belonging to the same angular class). It
then becomes possible to correct RQD in the same way as FF, but this time multiplying by the reverse
coefficient: 0.26, 0.71 and 0.96 for the classes [0°, 30°], [30°, 60°] and [60°, 90°] respectively. As
previously mentioned, and even if we neglect the tolerance linked to the angular class, it is a crude
calculation which over-estimates RQD because of the pieces that might become smaller than 10cm
when multiplied by the sinus factor. Figures 9 present the results.

<figure><img src="image_146.png" /><figcaption>Figure 9. RQD versus FF. Red curve: Priest&Hudson model; black curve: conditional expection of RQD given FF. a) FF & RQD are not corrected b) only FF is Terzaghi corrected c) RQD & FF are corrected</figcaption></figure>
We focus on [0, 20] FF interval and Figures 9a and 9b show a zoom of 8b and 8a respectively. When a correction is applied to RQD, it induces so strong a reduction of variability that this time, Priest & Hudson over-estimates RQD given FF when it is Terzaghi-corrected (Figure 9c). So the correction is not efficient, if we refer to Priest & Hudson, and the best solution is certainly not to apply any correction, neither to RQD nor to FF (Figure 9a). It that case, the estimations must be directional as shown previously.

## 5 Conclusion

The first conclusion is that RQD, as well as FF, are directional measurements that must be estimated after classing the samples according to their directions. Otherwise we mix quantities that are not comparable, producing a not interpretable result. But if we accept this directionality, it means that attributes like RMR become directional too, something which may disturb the practitioner's habits. If we assume that the directionality of the results is accepted and incorporated in the workflow of dimensioning mining sites, the procedure proposed here, i.e. classing the samples according to their directions is somewhat rudimentary as it proposes estimations only on the sampled directions and no interpolation between these directions. This is a very general problem that occurs when geostatistical tools are applied to tensors, and the solution, in the future, is probably to consider the real regionalization space involved by the problem, a 5D space, 3 for the usual coordinates and 2 for the dip crossed by the azimuth, a concept that requires considerable theoretical development because we do not have models for distances crossed by angles.

Now if we want a non directional result, and if we decide to consider Priest & Hudson's formula as a reference, one possibility is to correct the conditional bias so that Priest & Hudson's curve and conditional expectation are identical. Consider again the scatter diagram between RQD as it is and Terzaghi-corrected FF.

<figure><img src="image_147.png" /><figcaption>Figure 10. Proposal for a RQD correction. Red curve: Priest&Hudson model; black curve: conditional expectation of RQD given FF. A possible correction is to change RQD so that the red curve joins the black one. User can choose another model than Priest&Hudson's</figcaption></figure>

For each FF, one could remove from RQD the differential between the conditional average and the red curve so that both averages are the same. In other words, the conditional distribution is shifted by a quantity $\Delta(FF)$. This operation requires referring to the Priest & Hudson's formula but one can choose other references. The idea is to preserve the natural correlation between FF and RQD.

## Acknowledgment

The authors would like to acknowledge Codelco, Chile, for their strong support in the implementation of good geostatistical practices along the copper-business value chain, as well as the French government for supporting half of this work, and anonymous reviewers who greatly contributed to improving the quality of the manuscript. Many thanks to Claudio Rojas and his team of R&T, as well as to the geotechnicians of Chuqui. Without them, nothing would have been possible.
## References

Barton N, Lien R, Lunde J (1974) Engineering classification of rock mass for the design of tunnel support. NGI publication 106, Oslo, Rock Mechanics, vol.6 issue 4 pp 189-236

Bear J, Tsang C F, De Marsilly G (1993) Flow and contaminant transport in fractured rock, chapter 4 "Stochastic models of fracture systems and their use in flow and transport modeling", pp. 169-235, Academic Press, Inc., USA

Bieniawski Z T (1975) The Point Load Test in Geotechnical Practice. Eng. Geol., Sept., pp. 1-11

Deere D U, Hendron A J, Patton JF D, Cording E J (1967) Design of surface and near-surface construction in rock, Failure and Breakage of Rock, ed. C. Fairhurst, Soc. Of Min. Eng., AIME, N. Y., pp. 237-302.

Goodman R E, Smith H R (1980) RQD and fracture spacing, Journal of the Geotechnical Engineering Diviosn, ASCE 106: 191-193

Jaeger J C, Cook N G W (1969) Fundamentals of Rock Mechanics, Methuen, London. Palmstrom A (1982) T

Priest S D, Hudson J A (1976) Discontinuity Spacings in Rock. J. Rock Mech. Min. Sci. & Geomech, vol. 13, pp. 135-148, Pergamon Press, Great Britain.

Priest S D, Hudson J A (1981) Estimation of Discontinuity Spacing and Trace Length Using Scanline Survey. International Journal of Rock Mechanics and Mining Sciences, vol. 18, pp. 183-197.

Seguret S A, Guajardo C, Freire R (2014) Geostatistical Evaluation of Fracture Frequency and Crushing. In proceedings of Caving 2014, Third International Symposium on Block and Sublevel Caving, Santiago, Chile.

Sen Z, Kazi A (1984) Discontinuity Spacing and RQD Estimates from Finite Length Scanlines, International Journal of Rock Mechanics and Mining Sciences, vol; 21, PP. 203-212.

Stagg K G, Zienkiewicz O C (1968) Rock Mechanics in Engineering Practice, Wiley, N. Y., 442 pp.

Terzaghi RD (1965) Sources of error in joint surveys. Geotechnique, vol. 15 issue 3 pp 287-304

Wallis P F, King M S (1980) Discontinuity spacings in a crystalline rock. International Journal of Rock Mechanics and Mining Sciences, vol. 17, pp. 63-66.
Chapter F

The oral presentation with comments of
“Geostatistical Evaluation of Rock-Quality
Designation and its link with Linear Fracture
Frequency”

Serge A. Séguret (Mines ParisTech, France)

Cristian Guajardo (Codelco, Chile)

Presented at IAMG 2015, 17th annual conference of the International Association for Mathematical Geosciences, September 5-13, Freiberg, Germany
GEOSTATISTICAL EVALUATION OF ROCK QUALITY
DESIGNATION & ITS LINK WITH LINEAR FRACTURE
FREQUENCY

Serge Antoine Séguret – Mines-ParisTech, France
Cristian Guajardo, Codelco, Chile

<figure><img src="image_148.png" /></figure>

IAMG 2015

17th conference of the
International Association for
Mathematical Geosciences

September 5-13, 2015

Freiburg, Germany
e

<figure><img src="image_149.png" /></figure>
Rock Quality Designation (RQD)

<figure><img src="image_150.png" /></figure>

89

➤ Rock Quality Designation aims at quantifying the degree of jointing or fracturing in a rock mass

➤ It measures the borehole core recovery percentage incorporating only pieces of solid core that are longer than 10 cm in length measured along the centerline of the core

➤ It is one of the main attribute incorporated in the Rock Mass Rating which is a comprehensive index of rock mass quality used for the design and construction of excavations in rock, such as tunnels, mines, slopes and foundations
Additivity ?
(for a given sample direction)

<figure><img src="image_151.png" /></figure>

$$RQD_{L_1 \cup L_2} = 100 \frac{L_{>10\text{cm}, 1} + L_{>10\text{cm}, 2}}{L_1 + L_2}$$

$$= \frac{L_1 100 \frac{L_{>10\text{cm}, 1}}{L_1} + L_2 100 \frac{L_{>10\text{cm}, 2}}{L_2}}{L_1 + L_2}$$

$$= \frac{L_1 RQD_{L_1} + L_2 RQD_{L_2}}{L_1 + L_2}$$

→ Kriging possible
(same support &
sample direction)

90

► The first question concerns the ability to estimate RQD at any location by a linear combination of measurements like kriging in geostatistics

► Suppose that all the samples have been drilled along the same direction and take two values collected along two different supports L₁ and L₂, each with its own sum of core pieces longer than 10 cm

► One can see here that the value of RQD over the two supports is a linear combination of the two quantities, equal to the average when the supports are similar

► So in this case kriging RQD is authorized, at punctual or block scale, using classical geostatistical tools when RQD is order-two stationary

► Although the estimation is easy, the interpretation is not: as RQD is a 1D measurement, estimating it at block scale using 1D samples just gives the average behavior of 1D samples over a block, and this is not a 3D property
Anisotropy
("Directionality")

Geotechnical anisotropy ≠ Geostatistical anisotropy

<figure><img src="image_152.png" /></figure>

<figure><img src="image_153.png" /></figure>

➢ Now, we have to face this problem: the samples have different drilling directions

> It is important, at this stage, to distinguish between the interpretation of the word "anisotropy" when used by a geostatistician, and the same word, when used by a geotechnician or a hydrogeologist

➢ In geostatistics, this word means that when we calculate a variogram, based on differences of the variable of interest, the variance of such increments depends on the direction and this leads to concepts like “zonal” or “geometric” anisotropies

➢ In Geotechnics and in Hydrogeology, when one says that the phenomenon is anisotropic, it means that the measure itself depends on the direction, not only the increment

➢ The two anisotropies are linked, but the underlying concepts differ fundamentally and this is the reason why we employ the word “directionality” as a synonym for anisotropy as used by the geotechnicians and hydrogeologists
RQD Directional

<figure><img src="image_154.png" /></figure>

92

➤ The main reason of RQD directionality is that the spacing between the fractures is subject to bias

➤ Suppose that all the fractures are parallel planes aligned along a single direction and consider two consecutive fractures

➤ The spacing measured along direction n°1 is $l_1$ while it should be “11Sinus thetha” if we want to make the measurement independent of the sample direction by referring to direction n°2 perpendicular to the fracture direction

➤ This correction must be applied to each segment but it cannot be applied directly to its sum with the common factor “sinus thetha” because after correction, some segments may be shorter than 10 cm. This is the case of segment $l_2$ in the right-hand figure which will not be included in the summation

➤ The problem becomes more complex when we consider that a direction in three dimensions requires two angles (azimuth and dip)
RQD Directional

Suppose all samples same direction

→ Network (an)isotropy

<figure><img src="image_155.png" /></figure>

RQD →
* Degree of fracturation
* Not directional

<figure><img src="image_156.png" /></figure>

RQD →
* Degree of fracturation
* Ratio of anisotropy
* Directional

93

➤ Samples are directional, but fractures too

➤ If we assume that all the samples have the same direction, and the fracture network is isotropic, the directionality of RQD can be neglected, whereas if there are local anisotropies in the fracture network, the local value of RQD at a given location will depend on their angles and will not be comparable to another value of RQD at another location even close

➤ In the first case RQD measures the degree of fracturation in an objective way

➤ In the second case it also measures the anisotropy of the fracture network and this is not really its objective
<figure><img src="image_157.png" /></figure>

For all these reasons we recommend the following procedure:

the samples must first be classed according to their directions, and structural
analyses must be conducted by direction

If the variograms by sample directions are close to each other, and close to the variogram obtained with all the samples, there is no reason to distinguish between the directions because the fracture network looks isotropic

If this is not the case, kriging must be conducted only with the samples of a given directional sampling class, yielding as many directional results as there are classes

If, for each location where the estimation is conducted, the results are always the same, it means that the directionality of the measurement can be neglected and kriging can be conducted with all the samples without distinction between the classes

If the differences are large, the result is directional, which is important for the geotechnicians because it means that, when looking at the rock-strength prior to drilling a tunnel, for example, they must account for the drilling direction

These ideas are illustrated below by a case study
<figure><img src="image_158.png" /></figure>

* The data are from a copper mine in northern Chile. We have used more than 60,000 samples with RQD information

* All the samples are 1.5 m long

* We have at our disposal six sampling directions all belonging to East-West vertical plans and a dip is enough to characterize the sampling directions
Case study
Models by sample direction

<figure><img src="image_159.png" /></figure>

<figure><img src="image_160.png" /></figure>

<figure><img src="image_161.png" /></figure>

<figure><img src="image_162.png" /></figure>

<figure><img src="image_163.png" /></figure>

<figure><img src="image_164.png" /></figure>

<figure><img src="image_165.png" /></figure>

➢ We fit a model for each one of the six directions, and a seventh one taking into account all the samples

➢ Estimations are conducted by sampling direction, each time using only the samples associated with the direction
<figure><img src="image_166.png" /></figure>

➢ We present here a typical horizontal cross-section of the results

➢ The directional RQD maps can be regrouped into two sets: {1, 6} and {2, 3, 4, 5}

➢ The greatest differences concern the northern part which contains low RQD for the second set of directions
Case study
No direction selection

• Ordinary Kriging

• Moving neighborhood
≈ 30 points ANY Dir.

<figure><img src="image_167.png" /></figure>

➤ This map is obtained when no distinction of the sample direction is considered

➤ Such a map is of no interest because it just shows the dominant directions of the sampling
<figure><img src="image_168.png" /></figure>

* So direction 50° differs strongly from direction 105°

* Along 50° (and in a West-East plan), RQD is equal to 70% on average while it becomes lower than 20% when sampling is along 105°

* 70° correspond to 1 meter for the length of solid core; 20% corresponds to only 30 centimeters

* One can deduce that in this part of the deposit, the fractures tend to be parallel to 50°, with a high density

* Such analyses, made plane by plane and by direction may help to detect particular domains that go unnoticed if the sample direction is not taken into account

* Finally, while RQD aims at measuring the degree of jointing or fracturing in a rock mass, its spatial variations measure mainly the anisotropy of the fracture network

* The conclusion is that the directionality of RQD must be accepted and used and it means that attributes like Rock Mass Rating must become directional too
FF

<figure><img src="image_169.png" /></figure>

$$FF = \frac{\sum n(\theta)}{L}$$

<figure><img src="image_170.png" /></figure>

$$FF_{corrected} = \frac{\sum \frac{n(\theta)}{\sin(\theta)}}{L} \quad (\text{Terzaghi, 1965})$$

100

➤ Now let us recall the Fracture Frequency (FF), a number of fracture divided by a sample length L

➤ The fracture counting depends on the angle of the fracture and the sample direction, as shown by Terzaghi in 1965, and it must be corrected by a sinus

➤ This is the same problem as for RQD but here with a concrete solution as the fractures are usually classified according to their direction

➤ Question is: by applying such a correction on FF (and not RQD) do not we break some relationship and change results when integrating both measures in attributes like Rock Mass Rating ?
RQD & FF
(Priest & Hudson 1976)

$$RQD = 100e^{-0.1FF} (0.1FF + 1)$$

$$f(x) = FFe^{-FFx} \rightarrow \text{Frequency of the spacing x}$$

<figure><img src="image_171.png" /></figure>

In 1976, Priest and Hudson established a link between FF and RQD

This formula results from the assumption that the distribution of spacings between discontinuities along a line follows a negative exponential distribution

The aim here is not to discuss this assumption but to use this formula as a reference
RQD & FF
(Terzaghi or not Terzaghi)

<figure><img src="image_172.png" /></figure>

102

► What about RQD compared to FF? Left-hand figure shows a scatter diagram between FF (when Terzaghi-corrected) and RQD. The black curve represents the conditional expectation which is, for a given value of FF, the average of the different values of RQD. The behavior predicted by Priest & Hudson's formula is plotted in red. The least one can say is that this formula does not reflect the experimental average behavior. Does this mean that Priest & Hudson's hypotheses are not acceptable?

► Not necessarily. Right-hand figure presents RQD versus FF without the Terzaghi-correction. One can see that the formula slightly over-estimates RQD in the FF range [0, 28], but the differences are not comparable to the previous differences when Terzaghi was applied

► What happens? When we correct FF, we induce a distortion of the distribution, a kind of anamorphosis which makes the values no more comparable to RQD

► When FF is not corrected, we compare coherent populations because the bias in the counting is fully compatible with the bias in the associated length and then Priest & Hudson's formula is verified in an acceptable way

► So we have to face a dilemma: Terzaghi correction is necessary to make the discontinuity counting comparable from one direction to another and finally transform the directional number into a counting systematically perpendicular to the fractures. But if we correct FF, we must also correct RQD
RQD & FF
(Try a correction)

<figure><img src="image_173.png" /></figure>

103

➤ Let us try such a correction for RQD

➤ From the total amount of data, a sub-set of more than 4,000 samples was selected because they have at least 3 fractures with the same direction. It then becomes possible to correct RQD in the same way as FF, as previously explained a crude calculation which over-estimates RQD. We focus on [0, 20] FF interval

➤ On the left-hand figure, FF is corrected, not RQD. In the middle, no correction is applied. On the right-hand figure, RQD and FF are corrected

➤ When a correction is applied to RQD, it induces so strong a reduction of variability that this time, Priest & Hudson over-estimates RQD given FF when it is Terzaghi-corrected

➤ So the correction is not efficient, if we refer to Priest & Hudson, and the best solution is certainly to not apply any correction, neither to RQD, nor to FF
Conclusion

> RQD, FF, directional

→ RMR directional

> 5D Geostatistics?

104

> The first conclusion is that RQD, as well as FF, are directional measurements that must be estimated after classing the samples according to their directions. Otherwise we mix quantities that are not comparable, producing a not interpretable result

> It means that attributes like the Rock Mass Rating become directional too.

> But the procedure of classing the samples according to their directions is somewhat rudimentary as it proposes estimations only on the sampled directions and no interpolation between these directions

> This is a very general problem that occurs when geostatistical tools are applied to tensors, and the solution, in the future, is to consider the real regionalization space involved by the problem, a 5D space, 3 for the usual coordinates and 2 for the dip crossed by the azimuth, a concept that requires considerable theoretical development because we do not have models for distances crossed by angles.
Acknowledgment

<figure><img src="image_174.png" /></figure>

CODELCO
Produciendo Futuro
Codelco Chili

<figure><img src="image_175.png" /></figure>

<figure><img src="image_176.png" /></figure>

MINES
ParisTech

• Codelco

• Chile

• France

• Paris School of Mines

...und vergessen Sie
nicht:

<figure><img src="image_177.png" /></figure>
<figure><img src="image_178.png" /></figure>
# Chapter G

## Breccia Pipe Prediction: a new approach using non-stationary covariance

Francky Fouedjio (Mines ParisTech, France)

Serge A. Séguret (Mines ParisTech, France)

A paper presented at IAMG 2015, 17th annual conference of the International Association for Mathematical Geosciences, September 5-13, Freiberg, Germany

### Abstract

In this paper, we are interested in prediction of the elevation of the pipe surface named “Braden” of the El Teniente mine in Chile. This latter is one of the largest known porphyry-copper ore bodies. Previous approaches have been applied on this dataset by Séguret and Celhay (2013). Here the problem is tackled by using a geostatistical approach based on non-stationary covariances. The proposed methodology offers an integrated treatment of all aspects of non-stationarity: mean, variance and spatial continuity. The estimation of non-stationary parameters is based on the local stationarity assumption and is achieved using a three-step estimation scheme. First, a non-parametric kernel estimator of the local variogram is built. Then, it is used in a weighted local least squares procedure to estimate non-stationarity parameters at a representative set of points referred to as anchor points. Next, a kernel smoothing approach is used to interpolate the non-stationary parameters at any location of interest. Once the non-stationary parameters are estimated, they have integrated into akriging procedure. The proposed methodology has revealed an increased prediction accuracy when compared to standard stationary method, and demonstrated the ability to extract the underlying non-stationarity. Indeed, a comparison of predictions and prediction standard deviations maps indicates that the proposed non-stationary method captures some varying spatial features such as locally varying anisotropy in the data that are missing using a stationary method, the outcome appears much more realistic.
# Breccia Pipe Prediction: a new approach using non-stationary covariances

Francky Fouedjio Kameni ¹* and Serge Séguret²

¹Mineral Resources Flagship - CSIRO, Australia, [francky.fouedjiokameni@csiro.au](mailto:francky.fouedjiokameni@csiro.au)

²Center of Geosciences and Geoengineering - MINES ParisTech, France

* presenting author

## Abstract

In this paper, we are interested in the prediction of the breccia pipe elevation named "Braden" of the El Teniente mine in Chile. The problem is tackled by developing a new geostatistical approach based on non-stationary covariances. The proposed method offers an integrated treatment of all aspects of non-stationarity: mean, variance and spatial continuity. The estimation of non-stationary parameters is free distribution and carried out under the local stationarity assumption. The resulting estimated non-stationary parameters are naturally integrated into a kriging procedure or conditional simulations. The proposed approach has revealed an increased prediction accuracy when compared to the stationary one and demonstrated the ability to extract the underlying non-stationarity.

## 1 Introduction

A canonical problem in the geosciences is the prediction of a physical quantity over the whole region of interest from a finite set of irregular spaced data. This problem involves modeling and estimating the underlying spatial dependence structure of the observed data. Commonly, this is accomplished through statistical tools such as the variogram or covariogram computed on the whole domain of interest, under the stationarity assumption. However, in practice the stationarity assumption can be doubtful due to many factors, including specific landscape and topographic features of the region of interest or other localized effects. These local influences can be reflected computing local variograms, whose characteristics may vary across the domain of study. In such cases, carry out predictions based on a stationary approach could produce less accurate predictions, including an incorrect assessment of the estimation error (Stein, 1999).

Several approaches have been proposed for modeling and estimating non-stationary dependence structure (see Guttorp and Schmidt (2013), for a brief review). One of the most interesting is the explicit non-stationary covariances class proposed by Paciorek and Schervish (2006). However, the parameter estimation of these latter remains a crucial problem. In this work, we develop a procedure of estimating parameters that govern this class of closed-form non-stationary covariances under a single realization and local stationarity framework, through a step by step approach. First, we compute local variograms by a non-parametric kernel estimator. Then, it is used in a weighted local least squares procedure for estimating the parameters at a reduced set of representative points referred to as anchor points. Finally, a kernel smoothing method is used to interpolate the parameters at any location of interest. Then, the estimated non-stationary parameters are integrated naturally into a kriging procedure or conditional simulations.

As a motivating example in this work we consider the prediction of the breccia pipe elevation named "Braden" of the El Teniente mine in Chile. This latter is one of the largest known porphyry-copper ore bodies. The pipe is poorly mineralized and surrounded by different kinds of mineralized geological units. Knowing the exact location of the pipe surface is important, as it constitutes the internal limit of the deposit. Previous approaches have been applied on this dataset by Séguret and Celhay (2013).
The paper is structured as follows: the model formulation is described in Section 2. In Section 3, the statistical inference is detailed. Spatial predictions and conditional simulations are presented in Section 4. In Section 5 the proposed approach is applied on the breccia pipe datasets. Section 6 is a concluding remarks.

## 2 Model Formulation

Let $Y = \{Y(\mathbf{x}) : \mathbf{x} \in G \subseteq \mathbb{R}^p, p \ge 1\}$ be a random field defined on a fixed continuous domain of interest $G$ of the Euclidean space $\mathbb{R}^p$ and reflecting the underlying studied phenomenon. We consider that $Y$ is governed by the following model:

$$Y(\mathbf{x}) = m(\mathbf{x}) + \sigma(\mathbf{x})Z(\mathbf{x}), \quad \forall \mathbf{x} \in G, \qquad (1)$$

where: $m: \mathbb{R}^p \to \mathbb{R}$ is an unknown fixed function; $\sigma: \mathbb{R}^p \to \mathbb{R}^+$ is an unknown positive fixed function; $Z$ is a zero-expectation, unit variance random field with correlation function defined by:

$$R^{NS}(\mathbf{x}, \mathbf{y}) = \phi_{\mathbf{xy}} R^S \left( \sqrt{Q_{\mathbf{xy}}(\mathbf{x} - \mathbf{y})} \right), \quad \forall (\mathbf{x}, \mathbf{y}) \in G \times G, \qquad (2)$$

with: $\phi_{\mathbf{xy}} = |\Sigma_{\mathbf{x}}|^{\frac{1}{4}} |\Sigma_{\mathbf{y}}|^{\frac{1}{4}} \left| \frac{\Sigma_{\mathbf{x}} + \Sigma_{\mathbf{y}}}{2} \right|^{-\frac{1}{2}}, Q_{\mathbf{xy}}(\mathbf{h}) = \mathbf{h}^T \left( \frac{\Sigma_{\mathbf{x}} + \Sigma_{\mathbf{y}}}{2} \right)^{-1} \mathbf{h}, \forall \mathbf{h} \in \mathbb{R}^p$;

$\Sigma: \mathbb{R}^p \to PD_p(\mathbb{R})$, $\mathbf{x} \mapsto \Sigma_{\mathbf{x}}$ is a mapping from $\mathbb{R}^p$ to $PD_p(\mathbb{R})$ the set of real-valued positive definite $p$-dimensional square matrices; $R^S(.)$ is a continuous isotropic stationary correlation function, positive definite on $\mathbb{R}^p$, for all $p \in \mathbb{N}^*$.

The expression (2) represents the closed form non-stationary covariances class proposed by Paciorek and Schervish (2006). The construction of this class is based on a convolution of an orthogonal random measure with a spatially varying random weighting function (see Fouedjio et al. (2014) for more details). We can notice that this class gives non-stationary versions of some well known stationary correlation functions (Gaussian, exponential, Matérn and Cauchy) for a specific choice of $R^S(.)$. The intuition behind this class is that to each location $\mathbf{x}$ is assigned a local Gaussian kernel matrix $\Sigma_{\mathbf{x}}$ and the correlation between two locations $\mathbf{x}$ and $\mathbf{y}$ is calculated by averaging between the two local kernels at $\mathbf{x}$ and $\mathbf{y}$. In this way, the local characteristics at both locations influence the correlation of the corresponding target values. Thus, it is possible to account for non-stationarity. It is done by specifying the mapping $\Sigma(.)$ which models the anisotropy of the correlation function. The resulting kernel matrix $\Sigma_{\mathbf{x}}$ at each point $\mathbf{x}$ is interpreted as a locally varying geometric anisotropy matrix. It controls the anisotropic behavior of the random field in a small neighborhood around $\mathbf{x}$.

From model defined in (1), the two first moments of the random field $Y$ is given by:

$$E(Y(\mathbf{x})) = m(\mathbf{x}), \qquad (3)$$

$$\mathrm{Cov}(Y(\mathbf{x}), Y(\mathbf{y})) = \sigma(\mathbf{x})\sigma(\mathbf{y})R^{NS}(\mathbf{x}, \mathbf{y}) \equiv C^{NS}(\mathbf{x}, \mathbf{y}). \qquad (4)$$

Then, the non-stationarity of the random field $Y$ is characterized by the non-stationary parameters $m(.)$, $\sigma(.)$ and $\Sigma(.)$ defined at any location of the region of interest.

## 3 Statistical Inference

Let $\mathbf{Y} = (Y(\mathbf{s}_1), \dots, Y(\mathbf{s}_n))^T$ be a $(n \times 1)$ vector of observations from a unique realization of the random field $Y$, associated to known locations $\{\mathbf{s}_1, \dots, \mathbf{s}_n\} \subset G \subseteq \mathbb{R}^p$. The objective is to use the data $\mathbf{Y}$ to estimate the mean function $m(.)$, the standard deviation function $\sigma(.)$ and the correlation function determined by $\Sigma(.)$. The estimation of these parameters relies on the slightly local stationarity assumption which allows certain simplifications.
## 3.1 Local Stationarity

The local stationarity assumption Matheron (1971) implies that at any location $\mathbf{x}_0 \in G$ there exists a neighborhood $\mathcal{V}_{\mathbf{x}_0} = \{\mathbf{x} \in G, \|\mathbf{x} - \mathbf{x}_0\| \le b\}$ where the random field $\mathbf{Y}$ can be approximated by a stationary random field. Thus, $\forall (\mathbf{x}, \mathbf{y}) \in \mathcal{V}_{\mathbf{x}_0} \times \mathcal{V}_{\mathbf{x}_0}, m(\mathbf{x}) \approx m(\mathbf{y}) \approx m(\mathbf{x}_0)$ and $C^{NS}(\mathbf{x}, \mathbf{y}) \approx C^S(\mathbf{x} - \mathbf{y}; \mathbf{x}_0) = C^S(\mathbf{h}; \mathbf{x}_0), \|\mathbf{h}\| \le b$; where $C^S(\cdot)$ is a stationary covariance and the limit $b$ represents the radius of the local stationarity neighborhood $\mathcal{V}_{\mathbf{x}_0}$. In this way, the parameters are assumed to be very smooth functions which vary slowly over the domain. The expectation of the random field $\mathbf{Y}$ being approximately equal to a constant inside the local stationarity neighborhood, the resulting local covariance structure at any location $\mathbf{x}_0$ is written as follows:

$$C^S(\mathbf{h}; \mathbf{x}_0) = \sigma^2(\mathbf{x}_0) R^S \left( \sqrt{\mathbf{h}^T \mathbf{\Sigma}_{\mathbf{x}_0}^{-1} \mathbf{h}} \right), \quad \|\mathbf{h}\| \le b. \qquad (5)$$

Locally, the non-stationary covariance $C^{NS}(\cdot, \cdot)$ (4) is thus reduced to an anisotropic stationary one $C^S(\cdot)$ (5). The anisotropy function $\mathbf{\Sigma}(\cdot)$ is parametrized through the spectral decomposition and then the positive definiteness is guaranteed. Precisely, at any location $\mathbf{x}_0 \in G$, $\mathbf{\Sigma}_{\mathbf{x}_0} = \mathbf{\Psi}_{\mathbf{x}_0} \mathbf{\Lambda}_{\mathbf{x}_0} \mathbf{\Psi}_{\mathbf{x}_0}^T$, where $\mathbf{\Lambda}_{\mathbf{x}_0}$ is the diagonal matrix of eigenvalues and $\mathbf{\Psi}_{\mathbf{x}_0}$ is the eigenvector matrix. We assume working in 2D ($p = 2$) from now and we have:

$$\mathbf{\Lambda}_{\mathbf{x}_0} = \begin{pmatrix} \lambda_1^2(\mathbf{x}_0) & 0 \\ 0 & \lambda_2^2(\mathbf{x}_0) \end{pmatrix}, \quad \mathbf{\Psi}_{\mathbf{x}_0} = \begin{pmatrix} \cos \psi(\mathbf{x}_0) & \sin \psi(\mathbf{x}_0) \\ -\sin \psi(\mathbf{x}_0) & \cos \psi(\mathbf{x}_0) \end{pmatrix}, \quad \lambda_1(\mathbf{x}_0), \lambda_2(\mathbf{x}_0) > 0 \text{ and } \psi(\mathbf{x}_0) \in [0, \pi).$$

At each point, the square roots of the eigenvalues control the local ranges and the eigenvector matrix specify the local orientations. Thus, the anisotropy function $\mathbf{\Sigma}(\cdot)$ is characterized by the functions $\lambda_1(\cdot)$, $\lambda_2(\cdot)$ and $\psi(\cdot)$.

## 3.2 Local Variogram Kernel Estimator

Under the local stationarity assumption, we define a non-parametric kernel moment estimator of the stationary local variogram at a fixed location $\mathbf{x}_0 \in G$ and lag $\mathbf{h} \in \mathbb{R}^p$, $\gamma(\mathbf{h}; \mathbf{x}_0) = \sigma^2(\mathbf{x}_0) - C^S(\mathbf{h}; \mathbf{x}_0), \|\mathbf{h}\| \le b$ as follows:

$$\hat{\gamma}_{\epsilon}(\mathbf{h}; \mathbf{x}_0) = \frac{\sum_{V(\mathbf{h})} K_{\epsilon}^{*}(\mathbf{x}_0, \mathbf{s}_i) K_{\epsilon}^{*}(\mathbf{x}_0, \mathbf{s}_j) [Y(\mathbf{s}_i) - Y(\mathbf{s}_j)]^2}{2 \sum_{V(\mathbf{h})} K_{\epsilon}^{*}(\mathbf{x}_0, \mathbf{s}_i) K_{\epsilon}^{*}(\mathbf{x}_0, \mathbf{s}_j)}, \quad \|\mathbf{h}\| \le b, \qquad (6)$$

where the average (6) is taken over $V(\mathbf{h}) = \{(\mathbf{s}_i, \mathbf{s}_j) : \mathbf{s}_i - \mathbf{s}_j = \mathbf{h}\}$, the set of all pairs of locations separated by vector $\mathbf{h}$; $K_{\epsilon}^{*}(\mathbf{x}_0, \mathbf{s}_i) = K_{\epsilon}(\mathbf{x}_0, \mathbf{s}_i) / \sum_{l=1}^n K_{\epsilon}(\mathbf{x}_0, \mathbf{s}_l)$ are standardized weights; $K_{\epsilon}(\cdot, \cdot)$ is a non-negative, symmetric kernel on $\mathbb{R}^p \times \mathbb{R}^p$ with bandwidth parameter $\epsilon > 0$.

This moment estimator of the local variogram at any location $\mathbf{x}_0 \in G$ is a kernel weighted local average of squared differences of the regionalized variable. The kernel function is used to smoothly down-weight the squared differences (for each lag interval) according to the distance of these paired values from a target location. We assign to each data pair a weight proportional to the product of the individual weights. Observation pairs near to the target location $\mathbf{x}_0$ have more influence on the local variogram estimator than those which are distant.

To calculate the non-parametric kernel estimator (6), we choose an isotropic stationary Gaussian kernel: $K_{\epsilon}(\mathbf{x}, \mathbf{y}) \propto \exp(-\frac{1}{2\epsilon^2} \|\mathbf{x} - \mathbf{y}\|^2)$, $\forall (\mathbf{x}, \mathbf{y}) \in G \times G$. The latter has a non-compact support and therefore considers all observations. Thus, the local variogram estimator is not limited only to the local information, distant points are also considered. This avoids artefacts caused by the only use of observations close to the target location. It also reduces instability of the obtained local variogram at regions with low sampling density. Furthermore, it provides a smooth parameter estimate and then is compatible with the quasi-stationarity assumption. The size of the quasi-stationarity neighborhood $b$, it is set with respect to the bandwidth $\epsilon$. We take $b = \sqrt{3\epsilon}$ such that the standard deviation of the isotropic stationary Gaussian kernel matches the isotropic stationary uniform kernel (with compact support).
## 3.3 Parameter Estimation

The estimation of the parameters vector $\theta(\mathbf{x}_0) = (\sigma(\mathbf{x}_0), \lambda_1(\mathbf{x}_0), \lambda_2(\mathbf{x}_0), \psi(\mathbf{x}_0))$ which characterizes the stationary local variogram $\gamma(\cdot; \mathbf{x}_0) \equiv \gamma(\cdot; \theta(\mathbf{x}_0))$ at a fixed location $\mathbf{x}_0$ are found via the following minimization problem:

$$ \hat{\theta}(\mathbf{x}_0) = \arg\min_{\theta(\mathbf{x}_0) \in \Theta} \|\mathbf{w}_{\epsilon}(\mathbf{x}_0) \odot (\gamma(\theta(\mathbf{x}_0)) - \hat{\gamma}_{\epsilon}(\mathbf{x}_0))\|, \quad (7) $$

where $\odot$ is the product term by term; $\gamma^T(\theta(\mathbf{x}_0)) = [\gamma(\mathbf{h}_j)]_{j=1...J}$; $\hat{\gamma}_{\epsilon}^T(\mathbf{x}_0) = [\hat{\gamma}_{\epsilon}(\mathbf{h}_j; \mathbf{x}_0)]_{j=1...J}$; $\mathbf{w}_{\epsilon}^T(\mathbf{x}_0) = [w_{\epsilon}(\mathbf{h}_j; \mathbf{x}_0)]_{j=1...J}$, $w_{\epsilon}(\mathbf{h}; \mathbf{x}_0) = [(\sum_{V(\mathbf{h})} K_{\epsilon}^*(\mathbf{x}_0, \mathbf{s}_i) K_{\epsilon}^*(\mathbf{x}_0, \mathbf{s}_j))/\|\mathbf{h}\|]^{1/2}$; $\{\mathbf{h}_j \in \mathbb{R}^p, j = 1, ..., J\}$ are given lag vectors; $\theta(\mathbf{x}_0) \in \Theta$ is the vector of unknown parameters and $\Theta$ is an open parameter space.

Using the estimate of the vector of structural parameters $\hat{\theta}(\mathbf{x}_0)$ obtained in (7), the mean parameter $m(\mathbf{x}_0)$ is estimated explicitly by a local stationary kriging of the mean (Matheron, 1971).

For the prediction purpose, one needs to compute the non-stationary parameters $m(\cdot)$, $\sigma(\cdot)$ and $\Sigma(\cdot)$ at prediction and observation locations. In practice, it is unnecessary to solve the minimization problem (7) at each target location. Indeed, doing so is computationally intensive and redundant for close locations, since these estimates are highly correlated. To reduce the computational burden, the proposed idea consists in obtaining the parameter estimates only at some reduced set of $m \ll n$ representative points referred to as anchor points defined over the domain of interest. Then, using the estimates obtained at anchor points, a kernel smoothing method is used to make available estimates at any location of interest. We work with the Nadaraya-Watson kernel smoother which is appropriate and relatively simple. However, other smoothers can be used as well (local polynomials, splines, etc.). The choice of the smoothing bandwidth associated to Nadaraya-Watson kernel smoother is done through the generalized cross-validation criteria (Wand and Jones, 1995).

The estimation of the non-stationary parameters depends on the bandwidth parameter $\epsilon$ used in the computation of the local variogram non-parametric kernel estimator defined in (6). Indeed, the size of the local stationarity neighborhood is expressed in terms of this bandwidth parameter. The data-driven method used to select the bandwidth $\epsilon$ consists of leaving out one data location and using a form of cross-validation. Because the estimation of the spatial dependence structure is rarely a goal per se but an intermediate step before kriging, we want to choose the bandwidth that gives the best cross-validation mean square error.

# 4 Prediction

The main purposes of modelling and estimating the spatial dependence structure is to spatially interpolate data and perform conditional simulations. The expected benefit using the closed-form non-stationary covariances class (4) is to obtain spatial predictions and variance estimation errors more realistic than those based on an inadequate stationary covariances.

## 4.1 Kriging

Let $C^{NS}(\cdot, \cdot)$ the non-stationary covariance of the random field $Y$ and $m(\cdot)$ its mean. Given the vector of observations $\mathbf{Y} = (Y(\mathbf{s}_1), \dots, Y(\mathbf{s}_n))^T$ at $n$ fixed locations $\mathbf{s}_1, \dots, \mathbf{s}_n \in G$, the point predictor for the unknown value of $Y$ at unsampled location $\mathbf{s}_0 \in G$ is given by the optimal linear predictor:

$$ \hat{Y}(\mathbf{s}_0) = m(\mathbf{s}_0) + \sum_{i=1}^{n} \eta_i(\mathbf{s}_0)(Y(\mathbf{s}_i) - m(\mathbf{s}_i)). \quad (8) $$
The kriging weight vector $\eta = [\eta_i(s_0)]$ and the corresponding kriging variance $Q(s_0)$ are given by:

$$ \eta = C^{-1}C_0 \quad \text{et} \quad Q(s_0) = \sigma^2(s_0) - C_0^T C^{-1}C_0. \qquad (9) $$

where $C_0 = [C^{NS}(s_i, s_0)]$; $C = [C^{NS}(s_i, s_j)]$.

## 4.2 Conditional Simulations

Here we assume that the random field $Y$ is Gaussian with mean $m(.)$ and non-stationary covariance structure $C^{NS}(.,.)$. We want to simulate at a large number of locations a Gaussian random field with same mean and covariance, and ensure that the realization honors the observed values $Y(s_1), \dots, Y(s_n)$. This can be achieved from an unconditional simulation of the random field $Y$ as follows (Lantuejoul, 2002):

1. realize a unconditional simulation $\{X(s), s \in G\}$ of the random field $Y$;

2. carried out a simple kriging of $\{X(s) - Y(s), s \in G\}$ from its values taken at the data points $\{s_i, i = 1, \dots, n\}$, using $m(.)$ and $C^{NS}(.,.)$;

3. add the unconditional simulation and the result of kriging.

We have $Y(\mathbf{x}) = m(\mathbf{x}) + \sigma(\mathbf{x})Z(\mathbf{x}), \forall \mathbf{x} \in G$, where $Z$ is a Gaussian random field with zero expectation, unit variance and non-stationary correlation function $R^{NS}(.,.)$. Thus, to simulate the Gaussian random field $Y$ (step 1 of the previous algorithm), we need to known how we can simulate $Z$. Simulation of the Gaussian random field $Z$ can be carried out using a propagative version of the Gibbs sampler proposed by Lantuejoul and Desassis (2012). This algorithm allows to simulate a Gaussian vector at a large number of locations (comparatively to the existing classical algorithms such as Cholesky method or Gibbs sampler) without relying on a Markov assumption (it does not need to have a sparse precision matrix). The algorithm proposed in (Lantuejoul and Desassis, 2012) requires neither the inversion nor the factorization of a covariance matrix. Note that simulation methods such as spectral method or turning bands method are not adapted to the non-stationary case (Lantuejoul, 2002). The representation that underlies these methods relies on the stationarity assumption.

# 5 Application

The methodology presented in Section 4 has been applied to the elevation data of the breccia pipe called "braden" of the El Teniente mine in Chile. We have a training data (616 observations) which serves to calibrate the model and a validation data (200 observations) which serves only to assess the prediction performances. A comparison scheme of kriging under stationary and non-stationary models is carried out through a validation sample.

Raw estimates of non-stationarity parameters $m(.)$, $\sigma^2(.)$ and $\Sigma(.)$ at anchor points are shown respectively on Figures 1b, 1c and 1d. They are based on the non-stationary exponential covariance function. Concerning the estimated anisotropy function $\hat{\Sigma}(.)$ at anchor points, it is represented by ellipses as shown in Figure 1d. Based on these estimates, non-stationarity in the data is quite visible. Especially, from Figure 1d where we can clearly see the spatially varying azimuth. Such directional effects are also quite apparent on data (Figure 1a). Note that the stationary approach has not detected a global geometric anisotropy.
<figure><img src="image_179.png" /><figcaption>Figure 1: (a) Training data; (b) Estimated mean function $\hat{m}(.)$ at anchor points; (c) Estimated variance function $\hat{\sigma}^2(.)$ at anchor points; (d) Estimated anisotropy function $\hat{\Sigma}(.)$ at anchor points where the ellipses were scaled to ease vizualisation.</figcaption></figure>

Figure 2 shows the maps of smoothed parameters over the whole domain of observations: mean, variance, anisotropy ratio and azimuth. A visualization of the covariance at certain points (with all other points) via the level contours for estimated stationary and non-stationary models is presented in Figure 3. We can see how the non-stationary spatial dependence structure changes the shape from one place to another as compared to the stationary one. The stationary model is a nested isotropic model (nugget effect, exponential and spherical) while the non-stationary model corresponds to the non-stationary exponential covariance function.
<figure><img src="image_180.png" /><figcaption>Figure 2: Smoothed parameters over the domain of observations: (a) mean, (b) variance, (c) anisotropy ratio, (d) azimuth.</figcaption></figure>

<figure><img src="image_181.png" /><figcaption>Figure 3: Covariance level contours at few points for the estimated stationary and non-stationary models (a, b). Level contours correspond to the values: 30000 (black), 20000 (red) et 10000 (green).</figcaption></figure>

Table 1 presents the summary statistics for the external validation (200 hold-out sample) results using the classical stationary approach and the non-stationary proposed one. Some well-known discrepancy measures are used (Chilès (2012)), namely the Mean Absolute Error (MAE), the Root Mean Square Error (RMSE), the Normalized Mean Square Error (NMSE), the Logarithmic Score (LogS) and the Continued Rank Probability Score (CRPS). For RMSE, LogS and CRPS, the smaller the better; for MAE, the nearer to zero the better; for NMSE the nearer to one the better. Table 1 shows that the proposed approach outperforms the stationary one with respect to all the measures. The cost of non-using the non-stationary approach in this
case is substantial: in average the prediction at validation locations is about 23% better for the non-stationary approach than for the stationary one, in terms of RMSE.

<table><caption>Table 1: External validation on a set of 200 observations.</caption><thead><tr><th></th><th>Stationary</th><th>Non-stationary</th></tr></thead><tbody><tr><th>MAE</th><td>79.74</td><td>61.39</td></tr><tr><th>RMSE</th><td>154.41</td><td>117.99</td></tr><tr><th>NMSE</th><td>0.98</td><td>0.74</td></tr><tr><th>LogS</th><td>2439</td><td>2315</td></tr><tr><th>CRPS</th><td>123.59</td><td>121.11</td></tr></tbody></table>

The kriged values and the kriging standard deviations for the estimated stationary and non-stationary models are shown in Figure 4. The overall look of the predicted values and prediction standard deviations associated with each model differ notably. In particular, the proposed method takes into account certain local characteristics (such as locally varying anisotropy) of the regionalization that the stationary approach is unable to retrieve. Figure 5 shows some conditional simulations in the Gaussian framework, based on the estimated non-stationary model.

<figure><img src="image_182.png" /><figcaption>Figure 4: (a,b) Predictions and prediction standard deviations for the estimated stationary model. (c,d) Predictions and prediction standard deviations for the estimated non-stationary model.</figcaption></figure>
<figure><img src="image_183.png" /><figcaption>Figure 5: Conditional simulations based on the estimated non-stationary model.</figcaption></figure>

# 6 Conclusion

In this paper we are proposed a statistical methodology based on a non-stationary covariances class to predict the elevation of the breccia pipe elevation named "Braden" of the El Teniente mine in Chile. The estimation method offers an integrated treatment of all aspects of non-stationarity (mean, variance, covariance) in the modeling process and relies on the mild hypothesis of quasi-stationarity. The proposed method has revealed an increased prediction accuracy when compared to the standard stationary method, and demonstrated the ability to extract the underlying non-stationarity from a single realization. It also provides an exploratory analysis tool for the non-stationarity. Beyond the spatial predictions, we also show how conditional simulations can be carried out in this non-stationary framework.

# Acknowledgements

The authors would like to thank the company Codelco, Chile for providing the data used in this paper.
References

J. P. Chilès, P. Delfiner (2012) *Geostatistics: modeling spatial uncertainty*. Springer.

F. Fouedjio, N. Desassis, J. Rivoirard (2014). *A Generalized Convolution Model and Estimation for Non-stationary Random Functions*. arXiv:1412.1373.

P. Guttorp, A. M. Schmidt (2013). *Covariance Structure of Spatial and Spatio-temporal Processes*. Wiley Interdisciplinary Reviews: Computational Statistics 5, 279-287.

C. Lantuejoul, C. and N. Desassis (2012). *Simulation of a Gaussian Random Vector: A Propagative Version of the Gibbs Sampler*. The 9th International Geostatistics Congress.

C. Lantuejoul (2002). *Geostatistical Simulation: Models and Algorithms*. Springer, 256 p.

G. F. Matheron (1971). *The Theory of Regionalized Variables and its Applications*. Les cahiers du Centre de Morhologie Mathématique de Fontainebleau 5, Ecole Nationale Supérieure des Mines de Paris.

C. J. Paciorek, M. J. Schervish (2006). *Spatial Modelling Using a New Class of Nonstationary Covariance Functions*. Environmetrics 17, 483-506.

S. Séguret, F. Celhay (2013). *Geometric Modeling of a Breccia Pipe - Comparing Five Approaches*. Apcom - Application of Computers and Operations research in the Mineral Industry 1, 257-266.

M. L. Stein (1999). *Interpolation of Spatial Data: Some Theory for Kriging*. Springer.

M. Wand, C. Jones (1995). *Kernel Smoothing*. Monographs on Statistics and Applied Probability, Chapman and Hall.