Capping and kriging grades with long-tailed distributions

 

 

M. Maleki; N. Madani; X. Emery

Department of Mining Engineering & Advanced Mining Technology Center, University of Chile, University of Chile, Santiago, Chile

 

 

---

SYNOPSIS

Variogram analysis and kriging lack robustness in the presence of outliers and data with long-tailed distributions, which often arises when estimating grades in precious metal deposits. The capping technique, consisting of truncating the data to some top-cut grade, is widely used in order to mitigate the influence of the values in the upper tail of the distribution. However, this procedure omits part of the grade variability and is likely to provoke a bias in the estimates. To avoid these issues, a recently proposed approach is to decompose the grade of interest into three components (the truncated grade, a weighted indicator above the top-cut grade, and a zero-mean residual) and jointly estimate the truncated grade and the indicator by cokriging. This approach is attractive as it provides unbiased grade estimates, allows choosing the 'optimal' top-cut value, and essentially works with truncated and indicator data, thus avoiding the use of outlying values for calculating sample variograms and performing spatial interpolation.
This work presents an application of this approach to a disseminated gold deposit that has been identified through exploration drilling. The effect of using an indicator covariate is assessed through leave-one-out cross-validation, by comparing the grade estimates with the true grades and with the grade estimates obtained with the conventional capping approach, which considers only the truncated grade as the variable of interest. As a result, cokriging the truncated grade and the indicator above top-cut grade outperforms the conventional capping approach, yielding significantly more accurate estimates. A few complementary guidelines are provided for validating the model hypotheses and for the implementation of cokriging.

Keywords: Top-cut model, high values, outliers, cokriging, indicator.

---

 

 

Introduction

Geostatistics has received a great deal of attention in the mining industry over the last few decades, as it offers varied tools and models to assess mineral resources and ore reserves (Journel and Huijbregts, 1978; David, 1988; Krige, 1999). Geostatistical analyses are based on exploration data and/or production data such as the grade assays in drill-hole or blast-hole samples. Data integrity is essential to obtain accurate estimates of resources and reserves. In this regard, data with extreme high values may be problematic, because of the effect on the sample variograms and resource and reserve estimates (Krige and Magri, 1982; Armstrong, 1984). The modeller is often tempted to consider such extreme data as outliers or erroneous measurements and to omit them from geostatistical analyses. Such a procedure is questionable if the extreme data corresponds to true values, because it may lead to underestimate the resources or reserves. The bias can be severe when estimating the grades of elements with heavy-tailed distributions, such as gold or silver, for which removing the highest values can affect the economic appraisal of the ore deposit.

To mitigate this impediment, a number of sophisticated robust estimation procedures have been proposed to reduce the influence of extreme high values (Journel and Arik, 1988; Parker, 1991; Arik, 1992; Costa, 2003; Machado et at, 2011, 2012). A simple alternative is to truncate the extreme high values to some threshold or top-cut value, a procedure known as 'capping' or 'cutting' (Sinclair and Blackwell, 2002; Rossi and Deutsch, 2014). As an example, Costa (2003) proposes to truncate the values that deviate by more than their cross-validation standard deviation error. With this method, the same high value can result in different truncated values, depending on the the kriging neighbourhood, and no fixed top-cut is used. In practice, a fixed top-cut value is often considered, chosen (with some arbitrariness) in the last percentiles of the data distribution. This practice is widely used in the evaluation of precious metal deposits (David, 1988; Krige, 1999; Dagbert, 2005) and is accepted in current international codes for reporting mineral resources and ore reserves (SAMREC, 2007; JORC, 2012).

However, although less detrimental than removing extreme high values, considering truncated values omits some part of the data information and is likely to provoke a bias in the estimates, which has to be assessed and, if possible, fixed. On this subject, Rivoirard et al. (2013) recently presented a model in which a fixed top-cut value is used and the estimates are corrected for bias. The model is attractive from several points of view, as it allows choosing the 'optimal' top-cut value and working with truncated and indicator data, thus avoiding the use of extreme high values for calculating sample variograms and for performing spatial interpolation. This work presents an application of this top-cut model to a case study in mineral resources evaluation, together with proposed guidelines for practical implementation. Before dealing with the case study, a short review of the model is presented. Mathematical details and proofs can be found in the paper by Rivoirard et al. (2013).

 

Methodology

Spatial interpolation with top-cut model

Let x denote the vector of coordinates associated with a spatial location, Z the variable targeted for estimation (typically, the grade of an element of interest) and z_e a top-cut value. The variable (original grade) can be split into two components, a truncated grade and an excess, as follows:

where 1Z_(x)>ze denotes an indicator function:

In turn, the excess can be written as the sum of its regression onto the indicator and a residual:

where m⁺(z_e) is the mean value of Z above z_e, so that m⁺(z_e) - z_e is the mean excess above z_e, and R_ze (x) is the regression residual. The proposed model assumes that this residual has no spatial correlation with the indicator function or with the truncated grade:

Under such an assumption, the estimation of the original grade can be made by jointly estimating the truncated grade and the indicator via cokriging (CK) and separately estimating the residual via kriging (K):

If, furthermore, the residual has no spatial correlation structure, its estimate is equal to its mean value, i.e. zero. This results in the following estimate for the original grade:

The model therefore amounts to cokriging the truncated grade and the indicator function. Both of them no longer have extreme high values, in contrast to the original grade, which lends more robustness to the calculation of sample variograms and to the final estimates.

Validation of model assumptions

As assumed in Equation [4], the residual should be spatially uncorrelated with the truncated grade and the indicator. In other words, the high grade values (above z_e), which account for the non-zero residual values, should be independent of the low grade values (below z_e) and of the geometry of the set A (z_e) of locations with grade values greater than ze. Rivoirard et al. (2013) state that this occurs when there is no edge effects within A (z_e), i.e. when the value observed at a location of A(z_e) does not depend on whether or not this location is close to the boundary of A(z_e).

The absence of edge effects can be verified by examining the indicator variograms. Let us denote by yze (h) the variogram of the indicator function 1_Z(X)>Ze and by yze,z (h) the cross-variogram between the indicators 1_Z(X)>Ze and 1_Z(X)>Z with z greater than ze (under an assumption of second-order stationarity, these variograms are functions of the lag separation vector h). Then, in the absence of edge effects, one obtains (Rivoirard, 1994)

Choice of optimal top-cut value

The previous ratio of indicator variograms is the main tool for finding the most appropriate top-cut value ze. Indeed, one can select several thresholds (candidate top-cut values) {z₁,..., zn} and calculate the ratios between the indicator cross- and direct variograms associated with successive thresholds. The first threshold for which the ratio is approximately constant (independent of the lag separation vector) corresponds to the minimal acceptable value (zmin) for choosing the top-cut value ze.

Also, it is convenient (although not compulsory) to choose the top-cut value so that the residual variability above ze is pure nugget effect. In such a case, there is no advantage in choosing a higher top-cut value, since this would amount to incorporating poorly structured grade values in the truncated variable. In other words, when the residual variability above ze is pure nugget effect, there is little or no loss of information entailed by working with the truncated grade min{Z(x),ze} and the indicator 1_Z(X)>Ze instead of the original grade Z(x). Following Rivoirard et al. (2013), this condition can be checked by calculating a residual indicator variogram defined as a difference between normalized indicator variograms:

where z and z' (with z < z') are thresholds greater than z_mn, while T(z) and T(z') are the proportions of grade values above z and z, respectively. The threshold z for which the above difference is pure nugget provides the maximum acceptable top-cut value (zmax).

In summary, the examination of indicator variograms allows definition of an interval [z_mn, zmax] in which to choose the 'optimal' top-cut value ze. In addition to these considerations (Rivoirard et al., 2013), the following additional checks are proposed once the top-cut value has been chosen:

1) Calculate the truncated grade, indicator, and residual at each data location

2) Calculate the autocorrelation function (correlogram) of the residual, in order to determine whether or not it is a pure nugget effect. If so, the contribution of the residual can be omitted in the estimation of grade (Equation [6]). Otherwise, one may prefer the estimator that accounts for the residual (Equation [5])

3) Calculate the cross-correlogram between the residual and the truncated grade, as well as the cross-correlogram between the residual and the indicator, and check that these are close to zero for every lag separation vector. This corroborates that the residual has no spatial correlation with the truncated grade and with the indicator, as it is assumed in the model (Equation [4]).

The use of sample correlograms and cross-correlograms is suggested because of their robustness to outliers or to data with extreme high values (Isaaks and Srivastava, 1988), although these tools may be biased estimates of the true underlying spatial correlations because of the centring and normalization of the data for each lag separation vector (the magnitude of the bias depends on the number of data and their spatial distribution). Unbiased estimates of the spatial correlation could be obtained with the variogram or with the non-centred covariance, at the price of a loss of robustness (Rivoirard et al., 2000, Chiles and Delfiner, 2012).

Type of cokriging

Let us return to the estimator in Equation [5]. This estimator relies on cokriging the truncated grade and the indicator and, if not nugget, kriging the residual. Rivoirard et al. (2013) consider the use of simple and ordinary kriging/cokriging, depending on whether the mean values of the variables are deemed known or not.

The latter case (unknown mean values) is interesting in practice, insofar as the means can vary in space, provided that they remain approximately constant at the scale of the kriging or cokriging neighbourhood (assumption of local stationarity). When omitting the residual term and using the estimator in Equation [6] with ordinary cokriging, one further assumes that the mean grade above top-cut value, m⁺(z_e), is constant in space and known (but the probability that Z(x) exceeds ze may vary locally and is unknown), so that the local mean value of the residual is zero everywhere (Rivoirard et al., 2013). Now, if one denotes by m^-(ze) the mean grade below top-cut value, it is possible to establish a relationship between the mean values of the truncated grade and indicator, as shown in Table I.

Accordingly, the relationship between the mean values of the indicator (m₁) and truncated grade (m₂) can be modelled in the following fashion:

Even when considering that m₁ and m₂ are unknown, it is reasonable to assume that the previous relationship remains valid and to incorporate this relationship in the ordinary cokriging system (Emery, 2012). This is more restrictive than traditional ordinary cokriging, which assumes that m₁ and m₂ are unknown and unrelated (a somehow naive assumption, because of the strong dependence between the indicator and the truncated grade, as highlighted in Table I), but more versatile than simple cokriging, which assumes that both mean values are known without any uncertainty. For the coefficients of the relationship between m₁ and m₂ to be constant (Equation [9]), it is assumed that the mean grade below top-cut value, m^-(ze), is constant in space and known, while the probability that Z(x) exceeds or falls short of ze may be locally variable and unknown.

In the next section, these tools and models are applied to a case study in mineral resources estimation.

 

Case study: porphyry copper-gold deposit

Presentation of the data set and deposit

The available data consists of samples taken from 167 exploration drill-holes over an area of about 0.24 km² in a porphyry copper-gold deposit (Figure 1A). The mineralization is disseminated and subvertical in orientation, with a main direction dipping about 60° with respect to the horizontal plane. Mineralization comprises chalcopyrite, bornite, digenite, chalcocite, covellite, molybdenite, and pyrite. Gold mineralization occurs as native gold associated with gangue minerals, with particle sizes up to 160 µη, and as blebs within bornite and chalcocite.

In the following, it is of interest to estimate the gold grade within the deposit. To this end, the available samples have been composited to a length of 5 m. The distribution and statistics of the grade data are summarized in Figure 1B and Table II, indicating a long-tailed distribution with a median of 0.30 g/t, a mean of 1.01 g/t, and a maximum of almost 19 g/t. To preserve the confidentiality of the data, the original values have been multiplied by a constant factor.

Choice of top-cut value

A set of thresholds (from 3 g/t to 6 g/t) are considered as candidate top-cut values and the associated indicator variograms are calculated. From these, it is determined that, in every case, the ratio of cross-to-direct variograms (Equation [7]) does not vary significantly with the lag separation distance (Figure 2A, 2C, and 2E), denoting the absence of edge effects within the set of locations with grade greater than the thresholds under consideration. Following Rivoirard et al. (2013), any top-cut value greater than or equal to 3 g/t is therefore eligible. On the other hand, the variograms of indicator residuals (Equation [8]) display some spatial structure at short scales (up to approximately 20 m) (Figure 2B, 2D, and 2F). Accordingly, the maximum eligible top-cut value is not well-defined, since the residual variability above the candidate top-cut values is not pure nugget effect.

Based on these premises, a top-cut value of 3 g/t has been finally chosen, corresponding to the 92% percentile of the gold grade distribution. This choice instead of a higher top-cut value (4, 5, or 6 g/t) is motivated by the fact that the truncated gold grade is likely to exhibit a better spatial continuity, as it is no longer 'contaminated' by outlying data. The gold grade can therefore be split into three components (a truncated grade, an indicator associated with the chosen top-cut, and a residual), as per Equation [3]. To validate the model, it is convenient to check that there is no spatial correlation between the truncated grade and the residual, as well as between the indicator and the residual. This is done by calculating the cross-correlograms between the residual and the truncated grade and between the residual and the indicator, and verifying that these correlograms are close to zero for every separation distance (Figure 3).

Variogram modelling

The direct and cross-variograms of the truncated grade and the indicator are calculated along the identified main anisotropy directions and are fitted with a spherical structure (Figure 4A, 4B, and 4C):

The above equation gives a valid coregionalization model, as the eigenvalues of the sill matrix are non-negative (Wackernagel, 2003). The direction of main continuity (with a correlation range of 40 m) is dipping 60° with respect to the horizontal plane, in agreement with the known direction of mineralization, whereas the variograms are found to be isotropic with a correlation range of 7 m in the plane orthogonal to this direction. Note that the fitting relies mainly on the sample variograms at lag distances greater than 5 m (composite length), since few data pairs are involved in the calculation of the very first experimental point along each direction (only 1 pair for the direction of main continuity and 36 pairs for the orthogonal plane, while all the other experimental points involve several hundreds to thousands of data pairs).

It is also interesting to determine whether or not the residual is spatially correlated, in order to determine which estimator (Equation [5] or Equation [6]) is best suited to the data. To this end, instead of the traditional sample variogram of the residual, we calculated its sample correlogram, which is a more robust spatial continuity measure (Isaaks and Srivastava, 1988). This sample correlogram is then converted into a standardized variogram, showing the existence of a spatial correlation structure, although with a shorter range (28 m) than the truncated grade or the indicator. The residual variogram can be modelled by an anisotropic spherical structure, with the same anisotropy directions as the indicator and truncated grade (Figure 4D):

Resources estimation

According to the previous models, three approaches for estimating gold grade are compared:

1) Ordinary kriging of the truncated gold grade, corresponding to the traditional capping approach

2) Cokriging of the truncated gold grade and indicator, as proposed in Equation [6]. As stated previously, the mean values of both variables are assumed unknown but linearly related (Equation [9])

3) Cokriging of the truncated gold grade and indicator, as above, together with ordinary kriging of the residual, as proposed in Equation [5]. This third estimator is tested because the residual exhibits a spatial structure, as reflected in Figure 4D.

In each case, the same search neighbourhood is used, consisting of an ellipsoid with semi-axes 400 m along the main anisotropy direction and 100 m along the orthogonal directions. This ellipsoid is divided into octants and up to eight data are searched for in each octant. The dimensions of the search ellipsoid as well as the large number of searched data (up to 64 for each target block) have been chosen in order to obtain as accurate and conditionally unbiased estimates as possible, even if they inevitably yield smoothed grade models (Rivoirard, 1987; Krige, 1996, 1997b; Vann et al., 2003).

The traditional approach (kriging of truncated grade) leads to an estimate that is mostly less than the chosen top-cut grade (3 g/t) (Figure 5A), although the true grade is likely to exceed this top-cut, as suggested by the data histogram in Figure 1B. In contrast, with the cokriging approaches (Figure 5B and 5C), the grade estimates are no longer limited by the top-cut grade, which looks more realistic. The locations of high grade estimates are controlled by the indicator and by the residual, which correspond to the two corrective terms introduced in the traditional estimate (Equations [5] and [6]).

The differences between the three approaches can be assessed globally, by calculating the mean grade above different cut-off grades (Table III). It is seen that, for all the cut-offs, the traditional approach (kriging of truncated grade) yields biased estimates in comparison with the other two approaches. The bias can be explained because the traditional approach works only with the truncated grade and neglects the excess above the top-cut grade, whereas the other two approaches account for such an excess via the covariates (indicator and residual). Although the differences between the two top-cut models (with and without residual) are distinguishable on the maps drawn in Figure 5, their effect on the global statistics on the estimates is not significant, at least for low cut-offs; this is explained because, by construction, the residual has a zero mean value.

Cross-validation

In order to compare the model performances, we realized leave-one-out cross-validation: each data is successively removed and re-estimated from the surrounding data (Journel and Huijbregts 1978). The statistics on the cross-validation errors (mean error, mean absolute error, and mean squared error) are reported in Table IV, while the scatter diagrams between true and estimated gold grades are shown in Figure 6. These results confirm the bias of the traditional approach (mean error of -0.28 g/t), for which the estimates are practically limited by the top-cut grade. These estimates are also conditionally biased, insofar as the regression of the actual gold grade upon the estimated grade is significantly different from the diagonal, with a slope much greater than 1 (Figure 6A). The existence of a conditional bias is a pitfall for the resource/reserve model, as it leads to an erroneous assessment of the mining project (Krige, 1996, 1997a). In contrast, the top-cut model does not suffer any global or conditional bias in both cases (the mean errors are close to zero and the regressions of the actual upon the estimated grades are close to the diagonal line). A greater accuracy is achieved when accounting for the residual, as indicated by the smaller mean absolute and mean squared errors, and by the smaller dispersion of the scatter diagram around the diagonal line.

 

Conclusions

Top-cut modelling is a helpful tool to handle extreme high values in the estimation of variables with long-tailed distributions, as it improves the accuracy of the estimates and reduces global and conditional biases with respect to the traditional capping approach. Some additional guidelines have been proposed, consisting of checking the independence of the residual with the truncated grade and with the indicator above top-cut grade, through the use of cross-correlograms, and in replacing the estimators based on simple kriging and cokriging by ordinary kriging and cokriging with related mean values.

The case study presented shows that the residual should not be discarded when it has a spatial correlation structure. On the contrary, if the residual is pure nugget or if its structure cannot be interpreted properly (due to large fluctuations in the sample variogram or sample correlogram), the estimation can be limited to truncated grade and indicator cokriging, and the modelling can then be made in a convenient manner as it does not involve extreme high data values. In the case study, the resulting estimates turn out to be globally and conditionally unbiased, whether or not the residual is taken into account in the model. The top-cut model therefore achieves a trade-off between accuracy, simplicity of use, and robustness against extreme high values.

 

References

Arik, A. 1992. Outlier restricted kriging: a new kriging algorithm for handling of outlier high grade data in ore reserve estimation. Proceedings of the 23rd International Symposium on the Application of Computers and Operations Research in the Minerals Industry. Kim, Y.C. (ed.). Society of Mining Engineers, Littleton, CO. pp. 181-187.         [ Links ]

Armstrong M. 1984. Common problems seen in variograms. Mathematical Geology, vol. 16, no. 3. pp. 305-313.         [ Links ]

Chilès J.P. and Delfiner, P. 2012. Geostatistics: Modeling Spatial Uncertainty. Wiley, New York.         [ Links ]

Costa, J.F. 2003. Reducing the impact of outliers in ore reserves estimation. Mathematical Geology, vol. 35, no. 3. pp. 323-345.         [ Links ]

Dagbert, M. 2005. Geostatistics in resource/reserve estimation: a survey of the Canadian mining industry practice. Geostatistics Banff2004. Leuangthong, O. and Deutsch, C.V. (eds). Springer, Dordrecht. pp. 345-350.         [ Links ]

David, M. 1988. Handbook of Applied Advanced Geostatistical Ore Reserve Estimation. Elsevier, Amsterdam.         [ Links ]

Emery, X. 2012. Cokriging random fields with means related by known linear combinations. Computers and Geosciences, vol. 38, no. 1. pp. 136-144.         [ Links ]

Isaaks, E.H. and Srivastava, R.M. 1988. Spatial continuity measures for probabilistic and deterministic geostatistics. Mathematical Geology, vol. 20, no. 4. pp. 313-341.         [ Links ]

JORC. 2012. Australasian Code for Reporting of Exploration Results, Minerals Resources and Ore Reserves. Prepared by the Joint Ore Reserves Committee of the Australasian Institute of Mining and Metallurgy, Australian Institute of Geoscientists and Minerals Council of Australia. 44 pp.         [ Links ]

Journel, A.G., Arik, A. 1988. Dealing with outlier high grade data in precious metals deposits. Proceedings of the First Canadian Conference on Computer Applications in the Mineral Industry. Fytas, K., Collins, J.L., and Singhal, R.K. (eds). Balkema, Rotterdam. pp. 161-171.         [ Links ]

Journel, A.G. and Huijbregts, C.J. 1978. Mining Geostatistics. Academic Press, London.         [ Links ]

Krige, D.G. 1996. A basic perspective on the roles of classical statistics, data search routines, conditional biases and information and smoothing effects in ore block valuations. International Conference on Mining Geostatistics. Unisa Publishers, Pretoria. pp. 1-10.         [ Links ]

Krige, D.G. 1997a. A practical analysis of the effects of spatial structure and of data available and accessed, on conditional biases in ordinary kriging. Geostatistics Wollongong' 96. Baafi, E.Y., and Schofield, N.A. (eds). Kluwer Academic Publishers, Dordrecht. pp. 799-810.         [ Links ]

Krige D,G, 1997b. Block Kriging and the fallacy of endeavouring to reduce or eliminate smoothing. Proceedings of the 2nd Regional Symposium on the Application of Computers and Operations Research in the Minerals Industry. Puchkov, L.A. (ed.). Moscow State Mining University Publishing Center, Moscow.         [ Links ]

Krige, D.G. 1999. Essential basic concepts in mining geostatistics and their links with geology and classical statistics. South African Journal of Geology, vol. 102, no. 2. pp. 147-151.         [ Links ]

Krige, D.G. and Magri, E. 1982. Studies of the effects of outliers and data transformation on variogram estimates for a base metal and a gold ore body. Mathematical Geology, vol. 14, no. 6. pp. 557-567.         [ Links ]

Machado, R.S., Armony, M., Costa, J.F.C.L., and Koppe, J.C. 2011. Field parametric geostatistics - geostatistical rigorous solutions dealing with extreme grades during resource estimation. Proceedings of the 35th International Symposium on the Application of Computers and Operations Research in the Minerals Industry. Baafi, E.Y., Porter, I., and Kininmonth, R.J. (eds). Australasian Institute of Mining and Metallurgy, Melbourne. pp. 191-201.         [ Links ]

Machado, R.S., Armony, M., and Costa, J.F.C.L. 2012. Field parametric geostatistics - a rigorous theory to solve problems of highly skewed distributions. Geostatistics Oslo 2012. Abrahamsen, P., Hauge, R., and Kolbj0rnsen, O. (eds). Springer, Berlin. pp. 383-395.         [ Links ]

Parker, H.M. 1991. Statistical treatment of outlier data in epithermal gold deposit reserve estimation. Mathematical Geology, vol. 23, no. 2. pp. 175-199.         [ Links ]

Rivoirard, J. 1987. Two key parameters when choosing the kriging neighborhood. Mathematical Geology, vol. 19, no. 8. pp. 851-856.         [ Links ]

Rivoirard, J. 1994. Introduction to Disjunctive Kriging and Non-linear Geostatistics. Oxford University Press, Oxford.         [ Links ]

Rivoirard, J., Demange, C., Freulon, X., Lécureuil, A., and Bellot, N. 2013. A top-cut model for deposits with heavy-tailed grade distribution. Mathematical Geosciences, vol. 45, no. 8. pp. 967-982.         [ Links ]

Rivoirard, J., Simmonds, J., Foote, K.G., Fernandes, P., and Bez, N. 2000. Geostatistics for Estimating Fish Abundance. Wiley, New York.         [ Links ]

Rossi, M.E. and Deutsch, C.V. 2014. Mineral Resource Estimation. Springer, Dordrecht.         [ Links ]

Sinclair, A.J., and Blackwell, G.H. 2002. Applied Mineral Inventory Estimation. Cambridge University Press, Cambridge.         [ Links ]

SAMREC. 2007. The South African Code for the Reporting of Exploration Results, Mineral Resources and Mineral Reserves. Prepared by the South African Mineral Resource Committee, Working Group under the joint auspices of the Southern African Institute of Mining and Metallurgy and the Geological Society of South Africa, 50 pp.         [ Links ]

Vann, J., Jackson, S., and Bertoli, O. 2003. Quantitative kriging neighbourhood analysis for the mining geologist - a description of the method with worked case examples. Proceedings of the 5th International Mining Geology Conference. Dominy, S. (ed.). Australasian Institute of Mining and Metallurgy, Melbourne. pp. 215-223.         [ Links ]

Wackernagel, H. 2003. Multivariate Geostatistics: an Introduction with Applications. Springer, Berlin.         [ Links ] ♦