Serge Antoine Séguret
email: serge.seguret@mines-paristech.fr

Francky Fouedjio

For Codelco

Ramon Freire

Sebastian De La Fuente

Cristian Guajardo

Claudio Rojas

N° R151215sseg

Jacques Rivoirard

C Séguret

Rojas

S A Seguret

Francky Fouedjio Kameni
email: francky.fouedjiokameni@csiro.au

MODERN MINING GEOSTATISTICS CONVENTION WITH CODELCO YEAR 10 (2015) ANNIVERSARY REPORT (June 2016 version) For the Ecole des Mines

Global report of year 10 (2015). Contrat ARMINES

Introduction

During the 10 th 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

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, 7 th 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. 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 [START_REF] Matheron | Les Variables Régionalisées et leur estimation[END_REF]:

0 0 3 35 3 10.7 a a = -® = If ()
g denotes a variogram normalized by its sill, the underlying point sill C 0 of the exponential structure is produced by [START_REF] Journel | Mining Geostatistics[END_REF]:

0 (1 (3,3))

drill C C = -g
For the exponential structure, the charts in 1 yield: 0 [START_REF] Matheron | Les Variables Régionalisées et leur estimation[END_REF][START_REF] Matheron | Les Variables Régionalisées et leur estimation[END_REF] 0.087 0.055

C g = ® =
For l=15m, we deduce: 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 [START_REF] Séguret | Spatial sampling effect of laboratory practices in a porphyry copper deposit[END_REF].

Linear structure

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

( , ) 3

l l l b g = (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: ' ( , ) ( ', ') 3 l l l l l l b g -g = obtained by ( 4), the attenuation is 0.00015, a negligible quantity. When l=3m and l'=15m, 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.

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 [START_REF] Journel | Mining Geostatistics[END_REF] is not acceptable and we have to use charts that produce the exact calculation (see 1 , 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 S. A. Séguret "Blast & Drill holes..." Draft author. Final paper in: Proceedings of the 7th International Conference on Sampling and Blending , Kim H. Esbensen Editor, IMP publications, DOI 10.1255/t0sf.51

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 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

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:

2 2 ' ' ' ' ( ( ) ( )) [( ( ) ( )) ]) 2 ( ) 2 ( , ) ( , ) ( ', ') l l l l ll h D Y x Y x h E Y x Y x h h l l l l l l - + = - + = g = g -g -g (8)
S. A. Séguret "Blast & Drill holes..." Draft author. Final paper in: Proceedings of the 7th International Conference on Sampling and Blending , Kim H. Esbensen Editor, IMP publications, DOI 10.1255/t0sf.51 In (8), l' h represents the translation of the support l' by a vector h. ( , ) l l g and ( ', ') l l g 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 l l E l l g = g (9) ' ' [ ( )] [ ( , )] [ ( , )] ll h E h E l l E l l g = g -g (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.

' [ ( , )]

h E l l g 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:

17.5 12.5 1 [ ( , )] ( , ) E l l l l dl l g = g ò (11)
To evaluate the range of variations, the integral ( 11) is approximated by a finite sum:

13,14,15,16,17

1 [ ( , )] ( , ) 5 l E l l l l = g » g å (12)
We use the same charts as previously to calculate the values of ( , ) l l g involved and finally (12) yields: 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:

( ) ( ) ( ) blast m R h h h g = g + g (15)
with ( ) R h g , the nugget effect due to the blast error, having the variance ( ) h g , the point variogram, assumed to be isotropic;

15 [0,15] 2 1 ( ) ( ) ( | | 15)1 (| |) 15 m P h p p h h h = * = - + 15 1 ( ) ( ) 1 ( ) ( ) ( )
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" estimation5 . 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:

$ ! %&' + ! * 1 1 0 -. / 2 3 = . ! %&' + 4 * 5 1 3
In this system, R g 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 6! 7 8 %&' 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:

$ ! %&' + ! * 1 1 0 -. / 2 3 = . ! 7 9 %&' : (! 7 9 %&' )(0) + 4 * 5 1 3
The difference with the previous system is that in the second member, ! 7 8 %&' (capital "8 &' ") is replaced by:

15 15 ( , , ) ( , , ) ( ) m x y z x y m z p h h h h h u p h u du +¥ -¥ g * = g - ò (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):

; ! <' ! <',%&' 1 ! <',%&' ! %&' + ! * 1 1 1 0 > ; / /? 2 > = ;
! 7 9 <' 7 9 @ : (! 7 9 <' 7 9 @ )(0) ! 7 9 %&' 7 9 @ : (! 7 9 <' 7 9 @ )(0) + 4 * 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. 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.

Chapter

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?

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

-Comparing the measurement qualities by comparing their variograms -Establishing a formal link between the two measurements -Deducing linear systems enabling us to use the two types of measurements together

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)

7

Sampling 2015 -Bordeaux -France S A Séguret 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 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 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 10 Sampling 2015 -Bordeaux -France S A Séguret

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… 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.

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 indicatethe 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.

Finally we are able to establish a formal link between blast and drill measurements.

… 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.

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.

Ø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 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 microstructure 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.

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…. (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.

Good morning everyone

It is my pleasure to present the work entitled 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?

DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente

Mines-Planning 2015 -Antofagasta -Chile

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 Ø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 Ø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 depositand 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. Ø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. Ø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 Ø 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 convolutionis 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). Ø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 Ø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.

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. 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. Ø 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 Ø Now, we have to face this problem: the samples have different drilling directions Ø This is a very general problem that occurs when geostatical 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 Ø If, for a given estimation, based on a given amount of samplestypically the measurement used in the kriging neighborhoodwe 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 Ø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.

Ø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 "l1Sinus 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).

Ø 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.

Ø 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 Ø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 Ø 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.

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 openpit 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 [START_REF] Deere | Design of surface and near-surface construction in rock, Failure and Breakage of Rock[END_REF]Stagg & Zienkiesicz, 1968) and the (Linear) Fracture Frequency [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF]. 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 [START_REF] Priest | Estimation of Discontinuity Spacing and Trace Length Using Scanline Survey[END_REF][START_REF] Sen | Discontinuity Spacing and RQD Estimates from Finite Length Scanlines[END_REF]. 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 [START_REF] Terzaghi | Sources of error in joint surveys[END_REF], 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 [START_REF] Barton | Engineering classification of rock mass for the design of tunnel support[END_REF]. 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 [START_REF] Priest | Discontinuity Spacings in Rock[END_REF]. 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 [START_REF] Wallis | Discontinuity spacings in a crystalline rock[END_REF][START_REF] Goodman | RQD and fracture spacing[END_REF]). In the following, we do not address the thorny problem of making measurements made along lines representative of a volume.

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).

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). The spacing measured along direction n°1 is l 1 while it should be 1 sin l q 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( ) q 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.

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 10 , 1 100

cm cm L L L L RQD L L > > + = + 1 2 1 L L 1 2 1 2 = 10 , 1 10 , 2 1 2 1 2 1 2 100 100 cm cm L L L L L L L L > > + = + 1 2 1 2 1 2 L L RQD RQD L L L L + = +
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 rockstrength 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)

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 [START_REF] Seguret | Geostatistical Evaluation of Fracture Frequency and Crushing[END_REF]. 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.

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. [START_REF] Terzaghi | Sources of error in joint surveys[END_REF] proposed to correct this bias by multiplying the number of fractures by 1 sinq . As opposed to RQD, this correction is possible because the fractures are classed. In this paper, the coefficants 3.86, 1.414 and 1.035 are the multiplicative factors used for respectively [0°30°], [30°,60°], [60°,90°] classes.

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 [START_REF] Priest | Discontinuity Spacings in Rock[END_REF] Under the assumption that the distribution of spacings between discontinuities along a line follows a negative exponential one, the authors deduce that: In the theory developed by these authors, a distribution of spacing is associated with each value of FF, as follows:

( )

FFx f x FFe - =
with f(x), frequency of the spacing x (2) The quality of this assumption has been discussed in detail by de Marsily & Chilès in [START_REF] Bear | Flow and contaminant transport in fractured rock, chapter 4 "Stochastic models of fracture systems and their use in flow and transport modeling[END_REF] and it is not the aim of the present paper to continue this discussion. We use [START_REF] Journel | Mining Geostatistics[END_REF] to better illustrate the distortions of the correlation between FF and RQD, providing a kind of visual reference.

Application

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). All the drill holes belong to vertical West-East oriented planes and the sample directions can be summarized by six angular classes (Table 1):

Angular class Azimuth Dip 1 [245°, 295°] [-57°, -42°] 2 [245°, 295°] [-69°, -58°] 3 [245°, 295°] [-80°, -70°] 4 [66°, 115°] [-80°, -70°] 5 [66°, 115°] [-67°, -57°] 6 [66°, 115°] [-56°, -42°]

RQD

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

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.

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 7. Horizontal cross-section of RQD estimation. From a) to f), RQD maps by sampling directions. g) RQD map using all the samples together

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.

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. 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.

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. We focus on [0, 20] FF interval and Figures 9a and9b 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.

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 geostatical 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. 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 (FF) D

. 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.

Ø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 Ø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 1 and L 2 , 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 Ø 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 Ø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 "l1Sinus 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) Ø 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 an 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

93 GEOSTATISTICAL EVALUATION OF ROCK QUALITY DESIGNATION & ITS LINK WITH FRACTURE FREQUENCY S A Séguret, C Guajardo IAMG 2015 -Freiberg -Germany
Ø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 Ø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 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 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 ØNot necessarily. Right-hand figure presents RQD versus FF without the Terzaghicorrection. 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

Ø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 Ø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 geostatical 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.

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 [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF]. Several approaches have been proposed for modeling and estimating non-stationary dependence structure (see [START_REF] Guttorp | Covariance Structure of Spatial and Spatio-temporal Processes[END_REF], for a brief review ). One of the most interesting is the explicit non-stationary covariances class proposed by [START_REF] Paciorek | Spatial Modelling Using a New Class of Nonstationary Covariance Functions[END_REF]. 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 [START_REF] Séguret | Geometric Modeling of a Breccia Pipe -Comparing Five Approaches[END_REF].

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.

Model Formulation

Let Y = {Y (x) : x ∈ G ⊆ R p , p ≥ 1} be a random field defined on a fixed continuous domain of interest G of the Euclidean space R p and reflecting the underlying studied phenomenon. We consider that Y is governed by the following model:

Y (x) = m(x) + σ(x)Z(x), ∀x ∈ G, (1)
where: m : R p → R is an unknown fixed function; σ : R p → R + is an unknown positive fixed function; Z is a zero-expectation, unit variance random field with correlation function defined by:

R N S (x, y) = φ xy R S Q xy (x -y) , ∀(x, y) ∈ G × G, (2)
with:

φ xy = |Σ x | 1 4 |Σ y | 1 4 Σ x + Σ y 2 -1 2 , Q xy (h) = h T Σ x + Σ y 2 -1 h, ∀h ∈ R p ; Σ : R p → P D p (R), x → Σ
x is a mapping from R p to P D p (R) the set of real-valued positive definite p-dimensional square matrices; R S (.) is a continuous isotropic stationary correlation function, positive definite on R p , for all p ∈ N ⋆ .

The expression (2) represents the closed form non-stationary covariances class proposed by [START_REF] Paciorek | Spatial Modelling Using a New Class of Nonstationary Covariance Functions[END_REF]. The construction of this class is based on a convolution of an orthogonal random measure with a spatially varying random weighting function (see [START_REF] Fouedjio | A Generalized Convolution Model and Estimation for Non-stationary Random Functions[END_REF] 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 x is assigned a local Gaussian kernel matrix Σ x and the correlation between two locations x and y is calculated by averaging between the two local kernels at x and 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 Σ(.) which models the anisotropy of the correlation function. The resulting kernel matrix Σ x at each point x is interpreted as a locally varying geometric anisotropy matrix. It controls the anisotropic behavior of the random field in a small neighborhood around x.

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

E(Y (x)) = m(x), (3)
Cov(Y (x), Y (y)) = σ(x)σ(y)R N S (x, y) ≡ C N S (x, y). (4)
Then, the non-stationarity of the random field Y is characterized by the non-stationary parameters m(.), σ(.) and Σ(.) defined at any location of the region of interest.

Statistical Inference

Let Y = (Y (s 1 ), . . . , Y (s n )) T be a (n × 1) vector of observations from a unique realization of the random field Y , associated to known locations {s 1 , . . . , s n } ⊂ G ⊆ R p . The objective is to use the data Y to estimate the mean function m(.), the standard deviation function σ(.) and the correlation function determined by Σ(.). The estimation of these parameters relies on the slightly local stationarity assumption which allows certain simplifications.

Local Stationarity

The local stationarity assumption [START_REF] Matheron | The Theory of Regionalized Variables and its Applications[END_REF] implies that at any location x 0 ∈ G there exists a neighborhood V x 0 = {x ∈ G, xx 0 ≤ b} where the random field Y can be approximated by a stationary randon field. Thus, ∀(x, y) ∈ V x 0 × V x 0 , m(x) ≈ m(y) ≈ m(x 0 ) and C N S (x, y) ≈ C S (xy; x 0 ) = C S (h; x 0 ), h ≤ b; where C S (.) is a stationary covariance and the limit b represents the radius of the local stationarity neighborhood V 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 Y being approximately equal to a constant inside the local stationarity neighborhood, the resulting local covariance structure at any location x 0 is written as follows:

C S (h; x 0 ) = σ 2 (x 0 )R S h T Σ -1 x 0 h , h ≤ b. (5)
Locally, the non-stationary covariance C N S (., .) (4) is thus reduced to an anisotropic stationary one C S (.) [START_REF] Chilès | Geostatistics[END_REF]. The anisotropy function Σ(.) is parametrized through the spectral decomposition and then the positive definiteness is guaranteed. Precisely, at any location x 0 ∈ G,

Σ x 0 = Ψ x 0 Λ x 0 Ψ T
x 0 , where Λ x 0 is the diagonal matrix of eigenvalues and Ψ x 0 is the eigenvector matrix. We assume working in 2D (p = 2) from now and we have:

Λ x 0 = λ 2 1 (x 0 ) 0 0 λ 2 2 (x 0 ) , Ψ x 0 = cos ψ(x 0 ) sin ψ(x 0 ) -sin ψ(x 0 ) cos ψ(x 0 ) , λ 1 (x 0 ), λ 2 (x 0 ) > 0 and ψ(x 0 ) ∈ [0, π).
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 Σ(.) is characterized by the functions λ 1 (.), λ 2 (.) and ψ(.).

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 x 0 ∈ G and lag h ∈ R p , γ(h; x 0 ) = σ 2 (x 0 ) -C S (h; x 0 ), h ≤ b as follows:

γ ǫ (h; x 0 ) = V (h) K ⋆ ǫ (x 0 , s i )K ⋆ ǫ (x 0 , s j )[Y (s i ) -Y (s j )] 2 2 V (h) K ⋆ ǫ (x 0 , s i )K ⋆ ǫ (x 0 , s j ) , h ≤ b, (6)
where the average ( 6) is taken over V (h) = {(s i , s j ) : s is j = h}, the set of all pairs of locations separated by vector h; K ⋆ ǫ (x 0 , s i ) = K ǫ (x 0 , s i )/ n l=1 K ǫ (x 0 , s l ) are standardized weights; K ǫ (., .) is a non-negative, symmetric kernel on R p × R p with bandwidth parameter ǫ > 0.

This moment estimator of the local variogram at any location x 0 ∈ 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 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 ǫ (x, y) ∝ exp(-1 2ǫ 2 x -y 2 ), ∀(x, y) ∈ G × G.
The latter has a noncompact 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 ǫ. We take b = √ 3ǫ such that the standard deviation of the isotropic stationary Gaussian kernel matches the isotropic stationary uniform kernel (with compact support).

Parameter Estimation

The estimation of the parameters vector θ(x 0 ) = (σ(x 0 ), λ 1 (x 0 ), λ 2 (x 0 ), ψ(x 0 )) which characterizes the stationary local variogram γ(.; x 0 ) ≡ γ(.; θ(x 0 )) at a fixed location x 0 are found via the following minimization problem:

θ(x 0 ) = arg min θ(x 0 )∈Θ w ǫ (x 0 ) ⊙ (γ(θ(x 0 )) -γ ǫ (x 0 )) , (7)
where ⊙ is the product term by term ; γ T (θ

(x 0 )) = [γ(h j )] j=1...J ; γ T ǫ (x 0 ) = [ γ ǫ (h j ; x 0 )] j=1...J ; w T ǫ (x 0 ) = [w ǫ (h j ; x 0 )] j=1...J , w ǫ (h; x 0 ) = ( V (h) K ⋆ ǫ (x 0 , s i )K ⋆ ǫ (x 0 , s j ))/ h 1/2
; {h j ∈ R p , j = 1, . . . , J} are given lag vectors; θ(x 0 ) ∈ Θ is the vector of unknown parameters and Θ is an open parameter space.

Using the estimate of the vector of structural parameters θ(x 0 ) obtained in ( 7), the mean parameter m(x 0 ) is estimated explicitly by a local stationary kriging of the mean [START_REF] Matheron | The Theory of Regionalized Variables and its Applications[END_REF].

For the prediction purpose, one needs to compute the non-stationary parameters m(.), σ(.) and Σ(.) at prediction and observation locations. In practice, it is unnecessary to solve the minization 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 ≪ 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 works 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 [START_REF] Wand | Kernel Smoothing[END_REF].

The estimation of the non-stationary parameters depends on the bandwidth parameter ǫ 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 ǫ 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.

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 closedform non-stationary covariances class ( 4) is to obtain spatial predictions and variance estimation errors more realistic than those based on a inadequate stationary covariances.

Kriging

Let C N S (., .) the non-stationary covariance of the random field Y and m(.) its mean. Given the vector of observations Y = (Y (s 1 ), . . . , Y (s n )) T at n fixed locations s 1 , . . . , s n ∈ G, the point predictor for the unknown value of Y at unsampled location s 0 ∈ G is given by the optimal linear predictor:

Y (s 0 ) = m(s 0 ) + n i=1 η i (s 0 )(Y (s i ) -m(s i )).
(8)

The kriging weight vector η = [η i (s 0 )] and the corresponding kriging variance Q(s 0 ) are given by: η

= C -1 C 0 et Q(s 0 ) = σ 2 (s 0 ) -C T 0 C -1 C 0 . (9)
where C 0 = [C N S (s i , s 0 )]; C = [C N S (s i , s j )].

Conditional Simulations

Here we assume that the random field Y is Gaussian with mean m(.) and non-stationary covariance structure C N S (., .). 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 ), . . . , Y (s n ). This can be achieved from an unconditional simulation of the random field Y as follows [START_REF] Lantuejoul | Geostatistical Simulation: Models and Algorithms[END_REF]:

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

2. carried out a simple kriging of {X(s) -Y (s), s ∈ G} from its values taken at the data points {s i , i = 1, . . . , n}, using m(.) and C N S (., .) ;

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

We have Y (x) = m(x) + σ(x)Z(x), ∀x ∈ G, where Z is a Gaussian random field with zero expectation, unit variance and non-stationary correlation function R N S (., .). 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 [START_REF] Lantuejoul | Simulation of a Gaussian Random Vector: A Propagative Version of the Gibbs Sampler[END_REF]. 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 [START_REF] Lantuejoul | Simulation of a Gaussian Random Vector: A Propagative Version of the Gibbs Sampler[END_REF] 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 [START_REF] Lantuejoul | Geostatistical Simulation: Models and Algorithms[END_REF]. The representation that underlies these methods relies on the stationarity assumption.

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(.), σ 2 (.) and Σ(.) at anchor points are shown respectively on Figures 1b, 1c and1d. They are based on the non-stationary exponential covariance function. Concerning the estimated anisotropy function Σ(.) 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 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 nonstationary model corresponds to the non-stationary exponential covariance function. 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. The kriged values and the kriging standard deviations for the estimated stationary and nonstationary 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.

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.

Figure 3 .

 3 Figure 3. Drill hole variogram fitting. Dotted line, the experimental curve; continuous line, the model

  .Séguret "Blast & Drill holes..." Draft author. Final paper in: Proceedings of the 7th International Conference onSampling and Blending , Kim H. Esbensen Editor, IMP publications, DOI 10.1255/t0sf.51 When l=0m and l'=3m, and with

Figure 4 .

 4 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.

  y,z), the point grade assumed to be isotropic and devoid of any measurement error; "*" denotes a convolution product; 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:

  15

  Presented at WCSB7, 7 th world conference on sampling and blending, 10-12 June, Bordeaux, France Good morning everybody. I am pleased to present you 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.

  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

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.

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

"

  Diamond Drill-Holes, Blast-Holes and Cokriging" … a work supported by Codelco and the Paris School of Mines 36 DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente Mines-Planning 2015 -Antofagasta -Chile

  Bock modelling using blast and drill measurements together We test the different systems on a simulation 38 DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente Ø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. Ø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.

  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. 40 DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente Mines-Planning 2015 -Antofagasta -Chile

Ø

  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 squarred horizontally and 100 meters vertically Ø 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 43 DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente Mines-Planning 2015 -Antofagasta -Chile Ø As we are in a stationary case, matrices are presented using covariance C, it is simplier

  e. the blast without errors which we know because we work on a simulation where everything is known 44 DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente Mines-Planning 2015 -Antofagasta -Chile Ø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

Ø

  In the three scatter diagrams, the horizontal axis represents the true pointsupport 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 isdone. 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 47 DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente Mines-Planning 2015 -Antofagasta -Chile Ø 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.

•

  Block grade estimate by kriging using only blast holes 48 DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente Mines-Planning 2015 -Antofagasta -Chile Ø 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

  50 DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente Mines-Planning 2015 -Antofagasta -Chile 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.

  51 DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente Mines-Planning 2015 -Antofagasta -ChileThe 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. 52 DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente Mines-Planning 2015 -Antofagasta -Chile 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. 58 ANISOTROPY OF THE ROCK QUALITY DESIGNATION (RQD) & ITS GEOSTATISTICAL EVALUATION Ø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 1 and L 2 , each with its own sum of core pieces longer than 10 cm, respectively equal to L1>10 and L2>10... 59 ANISOTROPY OF THE ROCK QUALITY DESIGNATION (RQD) & ITS GEOSTATISTICAL EVALUATION Geomin 2015 -Antofagasta -Chile S A Séguret, C Guajardo, C Rojas Ø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

  ØThese ideas are illustrated below by a case study 67 ANISOTROPY OF THE ROCK QUALITY DESIGNATION (RQD) & ITS GEOSTATISTICAL EVALUATION Ø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 68 ANISOTROPY OF THE ROCK QUALITY DESIGNATION (RQD) & ITS GEOSTATISTICAL EVALUATION Ø All the drill holes belong to vertical West-East oriented planes and the sample directions can be summarized by six angular classes Ø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 71 ANISOTROPY OF THE ROCK QUALITY DESIGNATION (RQD) & ITS GEOSTATISTICAL EVALUATION Ø 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 Ø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 73 ANISOTROPY OF THE ROCK QUALITY DESIGNATION (RQD) & ITS GEOSTATISTICAL EVALUATION

Figure 1 .

 1 Figure 1. Detail of RQD calculation. Property highlighted in red is a cause of directionality

Figure 2 .

 2 Figure 2. RQD directionality, detail of the calculations

Figure 3 .

 3 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.

Figure 4 .

 4 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 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.

Figure 5 .

 5 Figure 5. Base maps of the samples.

Figure 8 .

 8 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

Figure 9 .

 9 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

Figure 10 .

 10 Figure 10. Proposal for a RQD correction. Red curve: Priest&Hudson model; black curve: conditional expection 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 Prist&Hudson's

ØThese

  ideas are illustrated below by a case study 94 GEOSTATISTICAL EVALUATION OF ROCK QUALITY DESIGNATION & ITS LINK WITH FRACTURE FREQUENCY Ø 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 95 GEOSTATISTICAL EVALUATION OF ROCK QUALITY DESIGNATION & ITS LINK WITH FRACTURE FREQUENCY S A Séguret, C Guajardo IAMG 2015 -Freiberg -Germany Ø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 96 GEOSTATISTICAL EVALUATION OF ROCK QUALITY DESIGNATION & ITS LINK WITH FRACTURE FREQUENCY S A Séguret, C Guajardo IAMG 2015 -Freiberg -Germany Ø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 Ø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 ØSo direction 50° differs strongly from direction 105°

  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 99 GEOSTATISTICAL EVALUATION OF ROCK QUALITY DESIGNATION & ITS LINK WITH FRACTURE FREQUENCY Ø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 ? GEOSTATISTICAL EVALUATION OF ROCK QUALITY DESIGNATION & ITS LINK WITH FRACTURE FREQUENCY S A Séguret, C Guajardo IAMG 2015 -Freiberg -Germany Ø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 101 GEOSTATISTICAL EVALUATION OF ROCK QUALITY DESIGNATION & ITS LINK WITH FRACTURE FREQUENCY S A Séguret, C Guajardo IAMG 2015 -Freiberg -Germany Ø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?

Figure 1 :

 1 Figure 1: (a) Training data; (b) Estimated mean function m(.) at anchor points; (c) Estimated variance function σ 2 (.) at anchor points; (d) Estimated anisotropy function Σ(.) at anchor points where the ellipses were scaled to ease vizualisation.

Figure 2 :

 2 Figure 2: Smoothed parameters over the domain of observations: (a) mean, (b) variance, (c) anisotropy ratio, (d) azimuth.

Figure 3 :

 3 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).

Figure 4 :

 4 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.

Figure 5 :

 5 Figure 5: Conditional simulations based on the estimated non-stationary model.

Table 1 .

 1 Sampling direction classes

Ø

  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 Ø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

	Ø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

Table 1 :

 1 External validation on a set of 200 observations.

		Stationary Non-stationary
	MAE	79.74	61.39
	RMSE	154.41	117.99
	NMSE	0.98	0.74
	LogS	2439	2315
	CRPS	123.59	121.11

1 >

Mines-Planning 2015 -Antofagasta -Chile

DIAMOND DRILL HOLES, BLAST HOLES & COKRIGING S A Séguret, S de La Fuente Mines-Planning 2015 -Antofagasta -Chile

S A Séguret, C Guajardo, C Rojas Geomin 2015 -Antofagasta -ChileHere is the program that we propose to follow.

Geomin 2015 -Antofagasta -Chile S A Séguret, C Guajardo, C Rojas

S A Séguret, C Guajardo IAMG 2015 -Freiberg -Germany

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.

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 governement 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.

Acknowledgements

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

(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. A paper presented at IAMG 2015, 17 th 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 porphyrycopper ore bodies. Previous approaches have been applied on this dataset by [START_REF] Séguret | Geometric Modeling of a Breccia Pipe -Comparing Five Approaches[END_REF]. Here the problem is tackeld by using a geostatistical approach based on nonstationary covariances. The proposed methodology offers an integrated treatment of all aspects of non-stationarity: mean, variance and spatial continuity. The estimation of nonstationary 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 nonstationarity 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 a kriging 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.