# Estimating Quartz Reserves Using Compositional Kriging

Authors: J. Taboada; Á. Saavedra; C. Iglesias; E. Giráldez
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

The aim of this study was to determine spatial distribution and volume of four commercial quartz grades, namely, silicon metal, ferrosilicon, aggregate, and kaolin (depending on content in impurities) in a quartz seam. The chemical and mineralogical composition of the reserves in the seam were determined from samples collected from outcrops, blasting operations, and exploratory drilling, and compositional kriging was used to calculate the volume and distribution of the reserves. A more accurate knowledge of the deposit ensures better mine planning, leading to higher profitability and an improved relationship with the environment.


---

## Body

## 1. Introduction

A problem in exploiting mineral deposits is the limited knowledge available regarding the continuity, size, and grade distribution of the deposit, all of which are factors that determine the exploitability and profitability of the reserves. This explains why studies of grade and volume have long been the focus of considerable interest, with exploration efforts typically extending over the life of the operation.Different approaches have been taken to reduce the uncertainty as to the characteristics of a deposit, based on using various mathematical techniques. These techniques include kriging—particularly fuzzy kriging [1] and compositional kriging [2–4] as the two most advanced prediction techniques. In these methods, the grade distribution of a deposit is predicted by dividing the deposit into cells of a specific size and calculating the percentages for the different grades present in each cell.Several works in the literature describe the use of kriging techniques to calculate reserves in mineral deposits. Badel et al. [5] applied median indicator kriging to estimate grades of iron in a deposit located in southern Iran. Sohrabian and Ozcelik [6] studied the exploitability of andesite blocks in a quarry in Turkey using independent component kriging. Kriging techniques have been used to evaluate mineral reserves in Spanish deposits in previous works by our research group. Thus, Taboada et al. [7] evaluated reserves of Rosa Porriño ornamental granite in a deposit in southern Pontevedra (NW Spain) using fuzzy kriging techniques and Saavedra et al. [8], subsequently modelled the same site using compositional kriging, obtaining improved results.The deposit used for this study contains not granite, but a seam of quartz. Saavedra et al. [9] used fuzzy kriging techniques to estimate the volume and the spatial distribution of four commercial grades in the seam, marketed according to alumina (Al2O3) content in silicon metal, ferrosilicon, kaolin, or aggregate. In the present study, compositional kriging is applied in order to improve the results obtained with fuzzy kriging. The objective is to determine the volume and spatial distribution of the quartz reserves so as enable better planning of mining operations. Better planning would optimize resource use and increase mining yield, with the resulting positive impact on the environment (reduced extraction and processing waste) and the company (increased profits).The article is structured as follows: first, we describe the methodology, the mathematical concepts used, and the field of study and explain how the data were processed. In the last two sections, we present the results obtained and the conclusions drawn from the research.


## 2. Materials and Methods

### 2.1. Compositional Kriging

Compositional data are understood to be a set of nonnegative vectors such that the sum of their components is constantk, which takes the value k=100 when the components are percentages and k=1 when the components are proportions. Thus, a compositional data point (or composition) can be interpreted as providing information on the relative importance of the various parts of some whole.Denoting the centroid of a primary block asx∈R3, the compositional random variable Z(x)=(Z1(x),…,Zp(x)) is defined such that the one-dimensional variables Zj(x),j=1,…,p, reflect the jth part of the composition. Thus, given a sample realization {z(x1),…,z(xn)} for the compositional random variable Z(x), it can be verified that each z(xi)=(z1(xi),…,zp(xi)) is a composition, and hence, zj(xi)≥0, ∑j=1pzj(xi)=k. The objective is to predict the value of Z(x0) for a new location of interest x0.Classical spatial statistical methods include predictions with kriging, based on one-dimensional variables, by means of linear combinations of the sample realizations:z^j(x0)=∑i=1nλizj(xi), j=1,…,p, or predictions with cokriging, based on the entire set of variables in accordance with the expression z^j(x0)=∑j=1p∑i=1nλi,jzj(xi), where λi and λi,j are weights determined by solving matrix systems. See [10–12] for a more detailed explanation of these classical prediction methods.It is a proven fact that none of these methods of prediction will preserve the characteristics of the compositional data; that is, the various components may be negative and so may not respect the constant sum, [4]. Consequently, these classical prediction techniques may give rise to erroneous results. The best way to analyse spatial dependence and interpolate compositional data is to apply a geostatistical technique that allows all the elements that form the composition to be handled jointly.Walvoort and de Gruijter [4] proposed a prediction method based on classical approaches. These authors include, in a matrix system, the constraints necessary for the predictions to take admissible values for the composite random function. Authors like [2] propose a transformation of the sample data z*(xi)=f(z(xi)),i=1,…,n, before applying any kriging-based forecasting technique. If the transformation function has been appropriately chosen, an admissible composition can be obtained by backtransforming the predicted values.Tolosana-Delgado [3] proposed a compositional kriging approach based on the properties of the simplex, a Euclidean space formed by the sample space of the compositions and endowed with the operations of addition, external product, and internal product. As in any Euclidean space, a vector can be represented by its coordinates relative to an orthonormal reference system of dimensions p-1. The spatial variance-covariance matrices of the coordinates have (p-1)×(p-1) elements and have all the properties of a variance-covariance matrix. The spatial dependence between coordinates and their subsequent interpolation can be tackled using conventional cokriging programs for regionalized variables. With a new change of coordinates, the predictions can be expressed as elements of the original Euclidean space, the simplex. Theoretical aspects and demonstrations of this methodology can be found in [2–4].The procedure can be briefly summarized as follows.(1)
Changing the coordinates, the sample space of composite data is transformed into a new space of dimensionp-1, z*(xi)=Ψln(z(xi)),i=1,…,n, where z*(xi)=(z1*(xi),…,zp-1*(xi))t, Ψ is the coordinate-change matrix, of dimension (p-1)×p, formed of vectors with an orthonormal basis arranged in columns, ln(z(xi))=(ln(z1(xi)),…,ln(zp(xi)))t, and where the superscript t indicates “transposed”.(2)
Obtained in the Euclidean space of dimensionp-1 is the prediction z^*(x0) of Z*(x0) using cokriging techniques.(3)
The prediction value ofZ(x0) is given by z^(x0)=C(exp(Ψtz^*(x0))), with C((ε1,…,εp)t)=(ε1/∑j=1pεj,…,εp/∑j=1pεj) as a normalization operator.This methodology, which ensures admissible compositional predictions, has been applied to geological compositions with spatial dependence in studies by Zhang et al. [13] and Saavedra et al. [8].In the second step, we need to calculate(p-1) autovariograms and (p-1)×(p-2)/2 cross-variograms. Semivariograms should verify that the variance of any linear combination of these variables is not negative. Put another way, it should be ensured that the prediction variance is always not negative. This could be checked by means of the Cauchy-Schwarz inequality. To solve this problem of model selection, a linear coregionalization model is typically used described briefly below.Let there beq variables, Yr, with zero mean, uncorrelated two by two, in such a way that all the variables can be expressed as a linear combination:
(1)Zj*(x)=mj+∑r=1qαjrYr(x),j=1,…,p-1,
where mj is the average of the variable Zj* and αjr coefficients. From the previous equation it follows that the semivariograms that model the spatial dependence of the variables Zj* and Zj* can be expressed as
(2)γjj′(h)=βjj′1γ1(h)+⋯+βjj′qγq(h),j,j′=1,…,p-1,
with βjj′r=αjrαj′r, r=1,…,q, with γr as the semivariogram function that reflects the spatial structure of the regionalized variable Yr. For further details, see [14].


### 2.2. Description of the Studied Deposit

The studied deposit is located in the Sierra del Pico Sacro, southeast of the city of Santiago de Compostela (NW Spain), at an altitude of about 530 metres above sea level.This deposit, with a vein of pure, white quartz extending along the marginal fault between Cira and Pico Sacro, has a fine to very fine-grained saccharoidal matrix with some coarser grains and with the quartz in lattices (Figure1, MAGNA Sheet no. 121). The seam runs N301E, with a dip varying from 53°–60° to the NE, and has an average potential of 70 metres.Figure 1
Geological map of the area where the quartz seam is located.The quartz seam is quarried for four products and classified according to alumina content as follows.(i)
Silicon metal (alumina content below 0.35%) is a high-quality product used as raw material in the production of aluminium and silicon alloys and in the manufacture of photovoltaic solar panels.(ii)
Ferrosilicon (alumina content below 0.6%) is a slightly lower-quality product than silicon metal and is used mainly in the production of steels and cast iron.(iii)
Aggregate (alumina content below 3%) is used in construction and in cement and concrete manufacture.(iv)
Kaolin (alumina content above 3%) is a white clay used in the ceramics and paper industry and in the manufacture of paints and plastics.Bearing in mind that the deposit contains products with different values depending on grade and specifications, it is clearly important to know the volume and spatial distribution of each ahead of mining operations.


### 2.3. Data Processing

A total of 29 profiles, equispaced by a distance of 100 metres, were constructed to cut across the seam. The profiles were considered to be representative of portions of the deposit extending 50 metres each way and so spanning a total of 100 metres. Information concerning the geology and mineralogy of the different sampled points was included in order to characterize, in terms of grade, the reserves intersected by each profile. Data were available from three sources.(i)
Outcrop inspection: field reconnaissance was conducted, and different outcrops in the area were sampled for subsequent analysis.(ii)
Borehole logging: 288 surveys for 6,048 drilled metres were performed before 2000, and 148 surveys for 10,038 drilled metres were performed between 2000 and 2009.(iii)
Sampling of blast material: data were collected from 436 blasting operations conducted in 2000 and in 2003.Thus, the profiles, covering the entire deposit, contained information on the grades in cross-sections of the deposit as determined from the aforementioned data (Figure 2).Figure 2
Profiles 100 metres equidistant constructed for application of the method. Grades are indicated as follows: silicon metal (blue), ferrosilicon (red), aggregate (green) and kaolin (yellow).Geostatistical data processing began by depicting a grid of cubic cells with 10-metre edges. The centroid of each cube contained information about the different percentages of each marketable product. Since the seam’s real dip (53°–60°) and direction (N301E) would have caused some of the cells to contain points outside the seam, without information (which would have implied an additional data processing burden for no benefit), it was necessary to reorient the seam was run in a north-south direction with a vertical dip. Boundary conditions as dictated by the characteristics of the deposit and by the operating method were applied to the grid.Using compositional kriging, the vertical grid of the geostatistical model was projected onto the established profiles. The distribution of grades for the entire deposit was then established, and kriging results were fitted to the cloud of points generated by the grid. The percentages of the four grades were calculated for each cell, and the distribution percentage was allocated to the centroids of the nearest cells, generating a database entry with the UTM coordinates (x,y,z) and the grade distributions for the entire mesh.Geostatistical calculations were performed using the R freeware program, which also programs and stores packets of operations [15]. Once the results were available, the deposit was plotted in 3D using Rhinoceros, with different colours assigned to each commercial grade: silicon metal (blue), ferrosilicon (red), aggregate (green), and kaolin (yellow).


## 2.1. Compositional Kriging

Compositional data are understood to be a set of nonnegative vectors such that the sum of their components is constantk, which takes the value k=100 when the components are percentages and k=1 when the components are proportions. Thus, a compositional data point (or composition) can be interpreted as providing information on the relative importance of the various parts of some whole.Denoting the centroid of a primary block asx∈R3, the compositional random variable Z(x)=(Z1(x),…,Zp(x)) is defined such that the one-dimensional variables Zj(x),j=1,…,p, reflect the jth part of the composition. Thus, given a sample realization {z(x1),…,z(xn)} for the compositional random variable Z(x), it can be verified that each z(xi)=(z1(xi),…,zp(xi)) is a composition, and hence, zj(xi)≥0, ∑j=1pzj(xi)=k. The objective is to predict the value of Z(x0) for a new location of interest x0.Classical spatial statistical methods include predictions with kriging, based on one-dimensional variables, by means of linear combinations of the sample realizations:z^j(x0)=∑i=1nλizj(xi), j=1,…,p, or predictions with cokriging, based on the entire set of variables in accordance with the expression z^j(x0)=∑j=1p∑i=1nλi,jzj(xi), where λi and λi,j are weights determined by solving matrix systems. See [10–12] for a more detailed explanation of these classical prediction methods.It is a proven fact that none of these methods of prediction will preserve the characteristics of the compositional data; that is, the various components may be negative and so may not respect the constant sum, [4]. Consequently, these classical prediction techniques may give rise to erroneous results. The best way to analyse spatial dependence and interpolate compositional data is to apply a geostatistical technique that allows all the elements that form the composition to be handled jointly.Walvoort and de Gruijter [4] proposed a prediction method based on classical approaches. These authors include, in a matrix system, the constraints necessary for the predictions to take admissible values for the composite random function. Authors like [2] propose a transformation of the sample data z*(xi)=f(z(xi)),i=1,…,n, before applying any kriging-based forecasting technique. If the transformation function has been appropriately chosen, an admissible composition can be obtained by backtransforming the predicted values.Tolosana-Delgado [3] proposed a compositional kriging approach based on the properties of the simplex, a Euclidean space formed by the sample space of the compositions and endowed with the operations of addition, external product, and internal product. As in any Euclidean space, a vector can be represented by its coordinates relative to an orthonormal reference system of dimensions p-1. The spatial variance-covariance matrices of the coordinates have (p-1)×(p-1) elements and have all the properties of a variance-covariance matrix. The spatial dependence between coordinates and their subsequent interpolation can be tackled using conventional cokriging programs for regionalized variables. With a new change of coordinates, the predictions can be expressed as elements of the original Euclidean space, the simplex. Theoretical aspects and demonstrations of this methodology can be found in [2–4].The procedure can be briefly summarized as follows.(1)
Changing the coordinates, the sample space of composite data is transformed into a new space of dimensionp-1, z*(xi)=Ψln(z(xi)),i=1,…,n, where z*(xi)=(z1*(xi),…,zp-1*(xi))t, Ψ is the coordinate-change matrix, of dimension (p-1)×p, formed of vectors with an orthonormal basis arranged in columns, ln(z(xi))=(ln(z1(xi)),…,ln(zp(xi)))t, and where the superscript t indicates “transposed”.(2)
Obtained in the Euclidean space of dimensionp-1 is the prediction z^*(x0) of Z*(x0) using cokriging techniques.(3)
The prediction value ofZ(x0) is given by z^(x0)=C(exp(Ψtz^*(x0))), with C((ε1,…,εp)t)=(ε1/∑j=1pεj,…,εp/∑j=1pεj) as a normalization operator.This methodology, which ensures admissible compositional predictions, has been applied to geological compositions with spatial dependence in studies by Zhang et al. [13] and Saavedra et al. [8].In the second step, we need to calculate(p-1) autovariograms and (p-1)×(p-2)/2 cross-variograms. Semivariograms should verify that the variance of any linear combination of these variables is not negative. Put another way, it should be ensured that the prediction variance is always not negative. This could be checked by means of the Cauchy-Schwarz inequality. To solve this problem of model selection, a linear coregionalization model is typically used described briefly below.Let there beq variables, Yr, with zero mean, uncorrelated two by two, in such a way that all the variables can be expressed as a linear combination:
(1)Zj*(x)=mj+∑r=1qαjrYr(x),j=1,…,p-1,
where mj is the average of the variable Zj* and αjr coefficients. From the previous equation it follows that the semivariograms that model the spatial dependence of the variables Zj* and Zj* can be expressed as
(2)γjj′(h)=βjj′1γ1(h)+⋯+βjj′qγq(h),j,j′=1,…,p-1,
with βjj′r=αjrαj′r, r=1,…,q, with γr as the semivariogram function that reflects the spatial structure of the regionalized variable Yr. For further details, see [14].


## 2.2. Description of the Studied Deposit

The studied deposit is located in the Sierra del Pico Sacro, southeast of the city of Santiago de Compostela (NW Spain), at an altitude of about 530 metres above sea level.This deposit, with a vein of pure, white quartz extending along the marginal fault between Cira and Pico Sacro, has a fine to very fine-grained saccharoidal matrix with some coarser grains and with the quartz in lattices (Figure1, MAGNA Sheet no. 121). The seam runs N301E, with a dip varying from 53°–60° to the NE, and has an average potential of 70 metres.Figure 1
Geological map of the area where the quartz seam is located.The quartz seam is quarried for four products and classified according to alumina content as follows.(i)
Silicon metal (alumina content below 0.35%) is a high-quality product used as raw material in the production of aluminium and silicon alloys and in the manufacture of photovoltaic solar panels.(ii)
Ferrosilicon (alumina content below 0.6%) is a slightly lower-quality product than silicon metal and is used mainly in the production of steels and cast iron.(iii)
Aggregate (alumina content below 3%) is used in construction and in cement and concrete manufacture.(iv)
Kaolin (alumina content above 3%) is a white clay used in the ceramics and paper industry and in the manufacture of paints and plastics.Bearing in mind that the deposit contains products with different values depending on grade and specifications, it is clearly important to know the volume and spatial distribution of each ahead of mining operations.


## 2.3. Data Processing

A total of 29 profiles, equispaced by a distance of 100 metres, were constructed to cut across the seam. The profiles were considered to be representative of portions of the deposit extending 50 metres each way and so spanning a total of 100 metres. Information concerning the geology and mineralogy of the different sampled points was included in order to characterize, in terms of grade, the reserves intersected by each profile. Data were available from three sources.(i)
Outcrop inspection: field reconnaissance was conducted, and different outcrops in the area were sampled for subsequent analysis.(ii)
Borehole logging: 288 surveys for 6,048 drilled metres were performed before 2000, and 148 surveys for 10,038 drilled metres were performed between 2000 and 2009.(iii)
Sampling of blast material: data were collected from 436 blasting operations conducted in 2000 and in 2003.Thus, the profiles, covering the entire deposit, contained information on the grades in cross-sections of the deposit as determined from the aforementioned data (Figure 2).Figure 2
Profiles 100 metres equidistant constructed for application of the method. Grades are indicated as follows: silicon metal (blue), ferrosilicon (red), aggregate (green) and kaolin (yellow).Geostatistical data processing began by depicting a grid of cubic cells with 10-metre edges. The centroid of each cube contained information about the different percentages of each marketable product. Since the seam’s real dip (53°–60°) and direction (N301E) would have caused some of the cells to contain points outside the seam, without information (which would have implied an additional data processing burden for no benefit), it was necessary to reorient the seam was run in a north-south direction with a vertical dip. Boundary conditions as dictated by the characteristics of the deposit and by the operating method were applied to the grid.Using compositional kriging, the vertical grid of the geostatistical model was projected onto the established profiles. The distribution of grades for the entire deposit was then established, and kriging results were fitted to the cloud of points generated by the grid. The percentages of the four grades were calculated for each cell, and the distribution percentage was allocated to the centroids of the nearest cells, generating a database entry with the UTM coordinates (x,y,z) and the grade distributions for the entire mesh.Geostatistical calculations were performed using the R freeware program, which also programs and stores packets of operations [15]. Once the results were available, the deposit was plotted in 3D using Rhinoceros, with different colours assigned to each commercial grade: silicon metal (blue), ferrosilicon (red), aggregate (green), and kaolin (yellow).


## 3. Results

The data obtained in the field constituted a set ofn=18267 compositional data, in such a way that their p=4 components totalled 100. These components represent the percentages of silicon metal, ferrosilicon, aggregate, and kaolin for each 10×10×10 cubic metre cell.Following a study of the main components of the sampled data, an orthogonal basis was obtained which determined the main directions of variability in the observations. This basis, previously normalized, constituted the orthonormal basis that would give rise to the coordinate-change matrixψ. As can be confirmed in [16], the orthonormal basis chosen has no great bearing on the final results. However, by choosing directions close to those of maximum variability, we aimed to reflect as accurately as possible the underlying grade distribution behaviour. Principal components method has been applied in order to know the directions close to those of maximum variability. Then, Gram-Schmidt process has been applied in order to get an orthonormal basis formed of the following vectors:
(3)e1=(0.4243,-0.2828,-0.1414,0.8485)t,e2=(0.4108,0.2783,-0.8311,-0.2512)t,e3=(-0.3221,0.8388,-0.011,0.4388)t.
Since several components in the sampled data were recorded with null values, a positive constant was added prior to the change of coordinates z*(xi)=Ψln(z(xi)), i=1,2,…,n. The sum of this constant does not change the compositional condition of the observations, but it should be taken into account in the final stage of the process to obtain predictions whose components sum 100.After transformation of the data, the new random variablesZj*, j=1,2,3, were structurally analysed. To fit the linear coregionalization model, r=2 uncorrelated variables with the following features were used: γγ1 pure nugget effect semivariogram with sill 1 and γ2 spherical semivariogram with range 300 and sill 1. The coefficients βjj′r, j,j′=1,2,3, r=1,2, that completed the linear co-regionalization model (estimated using the R freeware) are shown in Table 1.Table 1
Co-regionalization model coefficients used in the structural analysis of the transformed variablesZj*, j=1,2,3.
β
11
r
β
22
r
β
33
r
β
12
r
β
13
r
β
23
r
r
=
1
0.0517
0.0382
0.0219
0.0064
0.0060
0.0045
r
=
2
0.0917
0.0534
0.0003
−0.0074
0.0047
0.0016Figure3 shows the experimental semivariograms (circles) together with the fitted theoretical semivariograms γjj′, j,j′=1,2,3 (lines). Spatial correlation between the transformed variables Zj*,j=1,2,3, can be observed. Both Table 1 and Figure 3 show how the variable Z3 has negligible spatial correlation, since the coefficient β332 corresponding to the spherical semivariogram is practically zero. The reason for this behaviour is that the orthonormal basis has been chosen according to the main directions of variability in the observations.Figure 3
Experimental semivariograms (circles) together with the fitted theoretical semivariograms (lines) that describe the spatial correlation between the variablesZj*,j=1,2,3.Applying cokriging techniques and inverting the transformation applied to the data, we obtained grade predictions for the cells that were not included in the database. This completed the geological model of the seam, resulting in an estimate of quartz reserves in terms of the four commercial grades. The results of the cubication measurements as percentages representing the four grades are shown in Table2.Table 2
Commercial grade distribution.
Product
As % of all reserves
Silicon metal
17.15
Ferrosilicon
20.94
Aggregate
59.34
Kaolin
2.57
Total
100.00The value of the prediction variance obtained by compositional kriging was significantly lower than that obtained by the fuzzy kriging in [9]. This is because compositional kriging used real data for the percentage distribution of the four grades in the cells, whereas there was an information loss with fuzzy kriging due to the transformation of the input data to a triangular function.Figure4 shows the modelled deposit with the distribution in four grades indicated by different colours. Observe that, although the predominant grade is aggregate, there are some high quality areas of silicon metal and ferrosilicon. There is also evidence of the presence of a portion of the seam whose end use would be kaolin; this area is located in the central region of the seam, to the southeast and hidden by the surrounding areas in the image shown.Figure 4
Three-dimensional model of the site with the predictions of grade distributions [silicon metal (blue), ferrosilicon (red), aggregate (green), and kaolin (yellow)] and the terrain.


## 4. Conclusions

The compositional kriging method took into account that each block extracted from a deposit could contain different grades of quartz. Unlike other prediction methods, such as classical kriging or cokriging, compositional kriging ensures that the sum of percentages for the different grades adds to 100%. Despite the apparent complexity of the method, in practice is quite simple to implement in high-level languages such as program R. It was also possible to generate a three-dimensional model of the deposit using the program Rhinoceros.The quartz deposit was modelled and zoned in terms of four different quality grades—silicon metal, ferrosilicon, aggregate, and kaolin—defined by the market according to different specifications. The results satisfactorily fitted the geological reality of the seam.The knowledge acquired of a deposit in this way improves mine planning and enables selective exploitation of resources while reducing the amount of waste generated, all of which translates into higher profits and a better relationship with the environment.



---

*Source: 716593-2013-05-20.xml*